From: kronbichler Date: Fri, 5 Sep 2008 17:09:19 +0000 (+0000) Subject: Changed a few things to facilitate parallel run. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a73e01c311f694c0b7c0cb166cda306307de5ae1;p=dealii-svn.git Changed a few things to facilitate parallel run. git-svn-id: https://svn.dealii.org/trunk@16755 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 452a79e58a..01fbca07d0 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -572,8 +572,8 @@ void BoussinesqFlowProblem::setup_dofs () { stokes_dof_handler.distribute_dofs (stokes_fe); DoFRenumbering::Cuthill_McKee (stokes_dof_handler); - DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component); DoFRenumbering::subdomain_wise (stokes_dof_handler); + DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component); stokes_constraints.clear (); DoFTools::make_hanging_node_constraints (stokes_dof_handler, @@ -605,16 +605,16 @@ void BoussinesqFlowProblem::setup_dofs () n_T = temperature_dof_handler.n_dofs(); pcout << "Number of active cells: " - << triangulation.n_active_cells() - << " (on " - << triangulation.n_levels() - << " levels)" - << std::endl - << "Number of degrees of freedom: " - << n_u + n_p + n_T - << " (" << n_u << '+' << n_p << '+'<< n_T <<')' - << std::endl - << std::endl; + << triangulation.n_active_cells() + << " (on " + << triangulation.n_levels() + << " levels)" + << std::endl + << "Number of degrees of freedom: " + << n_u + n_p + n_T + << " (" << n_u << '+' << n_p << '+'<< n_T <<')' + << std::endl + << std::endl; stokes_partitioner.clear(); { @@ -628,7 +628,7 @@ void BoussinesqFlowProblem::setup_dofs () BlockSparsityPattern stokes_sparsity_pattern (2,2); - if (trilinos_communicator.MyPID() == 0) + //if (trilinos_communicator.MyPID() == 0) { BlockCompressedSetSparsityPattern csp (2,2); @@ -664,7 +664,7 @@ void BoussinesqFlowProblem::setup_dofs () BlockSparsityPattern stokes_preconditioner_sparsity_pattern (2,2); - if (trilinos_communicator.MyPID() == 0) + //if (trilinos_communicator.MyPID() == 0) { BlockCompressedSetSparsityPattern csp (2,2); @@ -693,7 +693,6 @@ void BoussinesqFlowProblem::setup_dofs () stokes_preconditioner_sparsity_pattern); } - std::cout << " bla " << std::flush; temperature_partitioner = Epetra_Map (n_T, 0, trilinos_communicator); { @@ -703,7 +702,7 @@ void BoussinesqFlowProblem::setup_dofs () SparsityPattern temperature_sparsity_pattern; - if (trilinos_communicator.MyPID() == 0) + //if (trilinos_communicator.MyPID() == 0) { CompressedSetSparsityPattern csp (n_T, n_T); DoFTools::make_sparsity_pattern (temperature_dof_handler, csp); @@ -836,6 +835,8 @@ void BoussinesqFlowProblem::assemble_stokes_system () stokes_rhs=0; + Vector old_temperatures (old_temperature_solution); + QGauss quadrature_formula(stokes_degree+2); QGauss face_quadrature_formula(stokes_degree+2); @@ -879,7 +880,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () std::vector > grads_phi_u (dofs_per_cell); std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - + const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); @@ -898,7 +899,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () local_matrix = 0; local_rhs = 0; - temperature_fe_values.get_function_values (old_temperature_solution, + temperature_fe_values.get_function_values (old_temperatures, old_temperature_values); for (unsigned int q=0; q::assemble_temperature_system () stokes_fe_values.reinit (stokes_cell); temperature_fe_values.get_function_values (old_temperature_solution, - old_temperature_values); + old_temperature_values); temperature_fe_values.get_function_values (old_old_temperature_solution, - old_old_temperature_values); + old_old_temperature_values); temperature_fe_values.get_function_gradients (old_temperature_solution, old_temperature_grads); @@ -1145,9 +1146,9 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_grads); temperature_fe_values.get_function_hessians (old_temperature_solution, - old_temperature_hessians); + old_temperature_hessians); temperature_fe_values.get_function_hessians (old_old_temperature_solution, - old_old_temperature_hessians); + old_old_temperature_hessians); temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), gamma_values); @@ -1157,17 +1158,17 @@ void BoussinesqFlowProblem::assemble_temperature_system () const double nu = compute_viscosity (old_temperature_values, - old_old_temperature_values, - old_temperature_grads, - old_old_temperature_grads, - old_temperature_hessians, - old_old_temperature_hessians, - present_stokes_values, - gamma_values, - global_u_infty, - global_T_range.second - global_T_range.first, - global_Omega_diameter, cell->diameter(), - old_time_step); + old_old_temperature_values, + old_temperature_grads, + old_old_temperature_grads, + old_temperature_hessians, + old_old_temperature_hessians, + present_stokes_values, + gamma_values, + global_u_infty, + global_T_range.second - global_T_range.first, + global_Omega_diameter, cell->diameter(), + old_time_step); for (unsigned int q=0; q::assemble_temperature_system () { for (unsigned int i=0; i::solve () gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner); pcout << " " - << solver_control.last_step() - << " GMRES iterations for Stokes subsystem." - << std::endl; + << solver_control.last_step() + << " GMRES iterations for Stokes subsystem." + << std::endl; stokes_constraints.distribute (stokes_solution); } @@ -1304,9 +1305,9 @@ void BoussinesqFlowProblem::solve () temperature_constraints.distribute (temperature_solution); pcout << " " - << solver_control.last_step() - << " CG iterations for temperature." - << std::endl; + << solver_control.last_step() + << " CG iterations for temperature." + << std::endl; double min_temperature = temperature_solution(0), max_temperature = temperature_solution(0); @@ -1319,8 +1320,8 @@ void BoussinesqFlowProblem::solve () } pcout << " Temperature range: " - << min_temperature << ' ' << max_temperature - << std::endl; + << min_temperature << ' ' << max_temperature + << std::endl; } } diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 14075ee65e..9ed483cdfb 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -26,10 +26,10 @@ namespace TrilinosWrappers SparseMatrix::const_iterator::Accessor:: visit_present_row () { - // if we are asked to visit the - // past-the-end line, then simply - // release all our caches and go on - // with life + // if we are asked to visit the + // past-the-end line, then simply + // release all our caches and go on + // with life if (this->a_row == matrix->m()) { colnum_cache.reset (); @@ -38,11 +38,11 @@ namespace TrilinosWrappers return; } - // otherwise first flush Trilinos caches + // otherwise first flush Trilinos caches matrix->compress (); - // get a representation of the present - // row + // get a representation of the present + // row int ncols; int colnums = matrix->n(); TrilinosScalar *values = new TrilinosScalar(colnums); @@ -52,12 +52,12 @@ namespace TrilinosWrappers ncols, &(values[0])); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - // copy it into our caches if the line - // isn't empty. if it is, then we've - // done something wrong, since we - // shouldn't have initialized an - // iterator for an empty line (what - // would it point to?) + // copy it into our caches if the + // line isn't empty. if it is, then + // we've done something wrong, since + // we shouldn't have initialized an + // iterator for an empty line (what + // would it point to?) Assert (ncols != 0, ExcInternalError()); colnum_cache.reset (new std::vector (colnums, colnums+ncols)); @@ -66,17 +66,17 @@ namespace TrilinosWrappers } - // The constructor is actually the - // only point where we have to check - // whether we build a serial or - // a parallel Trilinos matrix. - // In the end, it does not even - // matter how many threads there - // are, but only if we use an - // MPI compiler or a standard - // compiler. So, one thread on - // an MPI compiler will still get - // a parallel interface. + // The constructor is actually the + // only point where we have to check + // whether we build a serial or a + // parallel Trilinos matrix. + // Actually, it does not even matter + // how many threads there are, but + // only if we use an MPI compiler or + // a standard compiler. So, even one + // thread on a configuration with + // MPI will still get a parallel + // interface. SparseMatrix::SparseMatrix () : #ifdef DEAL_II_COMPILER_SUPPORTS_MPI @@ -160,21 +160,28 @@ namespace TrilinosWrappers { unsigned int n_rows = sparsity_pattern.n_rows(); + // TODO: As of now, assume that the + // sparsity pattern only sits at the + // zeroth processor (completely), let + // this determine the Trilinos + // sparsity pattern on that + // processor, and only broadcast the + // pattern afterwards. if (row_map.Comm().MyPID() == 0) { Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), ExcDimensionMismatch (matrix->NumGlobalRows(), sparsity_pattern.n_rows())); - // Trilinos seems to have a bug for - // rectangular matrices at this point, - // so do not check for consistent - // column numbers here. - // - // - // this bug is filed in the Sandia - // bugzilla under #4123 and should be - // fixed for version 9.0 + // Trilinos seems to have a bug for + // rectangular matrices at this point, + // so do not check for consistent + // column numbers here. + // + // + // this bug is filed in the Sandia + // bugzilla under #4123 and should be + // fixed for version 9.0 // Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(), // ExcDimensionMismatch (matrix->NumGlobalCols(), // sparsity_pattern.n_cols())); @@ -197,21 +204,17 @@ namespace TrilinosWrappers row_indices[col] = sparsity_pattern.column_number (row, col); matrix->InsertGlobalValues(row, row_length, - &values[0], &row_indices[0]); + &values[0], &row_indices[0]); } } + + last_action = Insert; - compress(); // In the end, the matrix needs to // be compressed in order to be - // really ready. However, that is - // a collective operation, so it - // has to be called on all processes - // by the user, whereas this function - // should only be used on one processor - // since our sparsity pattern data - // types are all serial as of now. - } + // really ready. + compress(); + } @@ -219,7 +222,39 @@ namespace TrilinosWrappers SparseMatrix::reinit (const Epetra_Map &input_map, const SparsityPattern &sparsity_pattern) { - reinit (input_map, input_map, sparsity_pattern); + // TODO: There seems to be problem + // in Epetra when a quadratic matrix + // is initialized with both row and + // column map. Maybe find something + // more out about this... + //reinit (input_map, input_map, sparsity_pattern); + matrix.reset(); + + unsigned int n_rows = sparsity_pattern.n_rows(); + + if (row_map.Comm().MyPID() == 0) + { + Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (input_map.NumGlobalElements(), + sparsity_pattern.n_rows())); + Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), + ExcDimensionMismatch (input_map.NumGlobalElements(), + sparsity_pattern.n_cols())); + } + + row_map = input_map; + col_map = input_map; + + std::vector n_entries_per_row(n_rows); + + for (unsigned int row=0; row + (new Epetra_FECrsMatrix(Copy, row_map, &n_entries_per_row[0], + false)); + + reinit (sparsity_pattern); } @@ -234,14 +269,14 @@ namespace TrilinosWrappers unsigned int n_rows = sparsity_pattern.n_rows(); if (row_map.Comm().MyPID() == 0) - { - Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), - ExcDimensionMismatch (input_row_map.NumGlobalElements(), - sparsity_pattern.n_rows())); - Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), - ExcDimensionMismatch (input_col_map.NumGlobalElements(), - sparsity_pattern.n_cols())); - } + { + Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (input_row_map.NumGlobalElements(), + sparsity_pattern.n_rows())); + Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), + ExcDimensionMismatch (input_col_map.NumGlobalElements(), + sparsity_pattern.n_cols())); + } row_map = input_row_map; col_map = input_col_map; @@ -253,7 +288,7 @@ namespace TrilinosWrappers matrix = std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, col_map, - &n_entries_per_row[0], true)); + &n_entries_per_row[0], false)); reinit (sparsity_pattern); } @@ -267,7 +302,7 @@ namespace TrilinosWrappers row_map = sparse_matrix.row_map; col_map = sparse_matrix.col_map; matrix = std::auto_ptr(new Epetra_FECrsMatrix( - *sparse_matrix.matrix)); + *sparse_matrix.matrix)); } @@ -343,9 +378,9 @@ namespace TrilinosWrappers void SparseMatrix::clear () { - // When we clear the matrix, - // reset the pointer and - // generate an empty matrix. + // When we clear the matrix, reset + // the pointer and generate an + // empty matrix. matrix.reset(); #ifdef DEAL_II_COMPILER_SUPPORTS_MPI row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), @@ -362,7 +397,7 @@ namespace TrilinosWrappers void SparseMatrix::compress () { - // flush buffers + // flush buffers int ierr; if (row_map.SameAs(col_map)) ierr = matrix->GlobalAssemble (); @@ -450,13 +485,13 @@ namespace TrilinosWrappers last_action = Add; } - // we have to do above actions in any - // case to be consistent with the MPI - // communication model (see the - // comments in the documentation of - // TrilinosWrappers::Vector), but we - // can save some work if the addend is - // zero + // we have to do above actions in any + // case to be consistent with the MPI + // communication model (see the + // comments in the documentation of + // TrilinosWrappers::Vector), but we + // can save some work if the addend + // is zero if (value == 0) return; @@ -480,8 +515,8 @@ namespace TrilinosWrappers Assert (matrix->Filled()==true, ExcMessage("Matrix must be compressed before invoking clear_row.")); - // Only do this on the rows - // owned locally on this processor. + // Only do this on the rows owned + // locally on this processor. int local_row = matrix->LRID(row); if (local_row >= 0) { diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 5405f9c9f6..d57c9afd14 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -24,8 +24,7 @@ namespace TrilinosWrappers { namespace internal { - VectorReference::operator TrilinosScalar () const // i believe useless with - // trilinos + VectorReference::operator TrilinosScalar () const { Assert (index < vector.size(), ExcIndexRange (index, 0, vector.size())); @@ -131,6 +130,7 @@ namespace TrilinosWrappers #endif vector = std::auto_ptr (new Epetra_FEVector(map)); + last_action = Insert; }