From 73e70d1e0c939c48194ca833f62ebc9a142850b2 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 31 May 2013 08:14:35 +0000 Subject: [PATCH] Make distribute() work properly with parallel::distributed::Vector. Also cleanup in extracting locally owned dofs of Trilinos vector. Distributing constraints with parallel::distributed::Vector has become simpler now, so use that in step-48. git-svn-id: https://svn.dealii.org/trunk@29678 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-48/step-48.cc | 68 ++++++++++--------- .../deal.II/lac/constraint_matrix.templates.h | 60 ++++++++-------- .../deal.II/lac/trilinos_vector_base.h | 12 +--- tests/mpi/step-48.cc | 16 ++--- 4 files changed, 71 insertions(+), 85 deletions(-) diff --git a/deal.II/examples/step-48/step-48.cc b/deal.II/examples/step-48/step-48.cc index 2b3f7a651a..9e2f4a2466 100644 --- a/deal.II/examples/step-48/step-48.cc +++ b/deal.II/examples/step-48/step-48.cc @@ -280,7 +280,7 @@ namespace Step48 void make_grid_and_dofs (); void oldstyle_operation (); void assemble_system (); - void output_results (const unsigned int timestep_number) const; + void output_results (const unsigned int timestep_number); #ifdef DEAL_II_WITH_P4EST parallel::distributed::Triangulation triangulation; @@ -422,32 +422,34 @@ namespace Step48 // This function prints the norm of the solution and writes the solution // vector to a file. The norm is standard (except for the fact that we need // to be sure to only count norms on locally owned cells), and the second is - // similar to what we did in step-40. However, we first need to generate an - // appropriate vector for output: The ones we used during time stepping - // contained information about ghosts dofs that one needs write access to - // during the loops over cell. However, that is not the same as needed when - // outputting. So we first initialize a vector with locally relevant degrees - // of freedom by copying the solution (note how we use the function @p - // copy_from to transfer data between vectors with the same local range, but - // different layouts of ghosts). Then, we import the values on the ghost - // DoFs and then distribute the constraints (as constraints are zero in the - // vectors during loop over all cells). + // similar to what we did in step-40. Note that we can use the same vector + // for output as we used for computation: The vectors in the matrix-free + // framework always provide full information on all locally owned cells + // (this is what is needed in the local evaluations, too), including ghost + // vector entries on these cells. This is the only data that is needed in + // the integrate_difference function as well as in DataOut. We only need to + // make sure that we tell the vector to update its ghost values before we + // read them. This is a feature present only in the + // parallel::distributed::Vector class. Distributed vectors with PETSc and + // Trilinos, on the other hand, need to be copied to special vectors + // including ghost values (see the relevant section in step-40). If we + // wanted to access all degrees of freedom on ghost cells, too (e.g. when + // computing error estimators that use the jump of solution over cell + // boundaries), we would need more information and create a vector + // initialized with locally relevant dofs just as in step-40. Observe also + // that we need to distribute constraints for output - they are not filled + // during computations (rather, they are distributed on the fly in the + // matrix-free method read_dof_values). template void - SineGordonProblem::output_results (const unsigned int timestep_number) const + SineGordonProblem::output_results (const unsigned int timestep_number) { - parallel::distributed::Vector locally_relevant_solution; - locally_relevant_solution.reinit (dof_handler.locally_owned_dofs(), - locally_relevant_dofs, - MPI_COMM_WORLD); - locally_relevant_solution.copy_from (solution); - locally_relevant_solution.update_ghost_values(); - constraints.distribute (locally_relevant_solution); + constraints.distribute (solution); Vector norm_per_cell (triangulation.n_active_cells()); - locally_relevant_solution.update_ghost_values(); + solution.update_ghost_values(); VectorTools::integrate_difference (dof_handler, - locally_relevant_solution, + solution, ZeroFunction(), norm_per_cell, QGauss(fe_degree+1), @@ -464,7 +466,7 @@ namespace Step48 DataOut data_out; data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (locally_relevant_solution, "solution"); + data_out.add_data_vector (solution, "solution"); data_out.build_patches (); const std::string filename = @@ -502,10 +504,10 @@ namespace Step48 // the finest mesh size. The finest mesh size is computed as the diameter of // the last cell in the triangulation, which is the last cell on the finest // level of the mesh. This is only possible for Cartesian meshes, otherwise, - // one needs to loop over all cells). Note that we need to query all the + // one needs to loop over all cells. Note that we need to query all the // processors for their finest cell since the not all processors might hold // a region where the mesh is at the finest level. Then, we readjust the - // time step a little to hit the final time exactly if necessary. + // time step a little to hit the final time exactly. template void SineGordonProblem::run () @@ -556,15 +558,15 @@ namespace Step48 // // Note how this shift is implemented: We simply call the swap method on // the two vectors which swaps only some pointers without the need to copy - // data around. Obviously, this is a more efficient way to move data - // around. Let us see what happens in more detail: First, we exchange - // old_solution with old_old_solution, which - // means that old_old_solution gets - // old_solution, which is what we expect. Similarly, - // old_solution gets the content from solution - // in the next step. Afterward, solution holds - // old_old_solution, but that will be overwritten during this - // step. + // data around. Obviously, this is a more efficient way to update the + // vectors during time stepping. Let us see what happens in more detail: + // First, we exchange old_solution with + // old_old_solution, which means that + // old_old_solution gets old_solution, which is + // what we expect. Similarly, old_solution gets the content + // from solution in the next step. Afterward, + // solution holds old_old_solution, but that + // will be overwritten during this step. unsigned int timestep_number = 1; Timer timer; diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h index 5758495937..ed19a1c0d6 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h @@ -733,41 +733,36 @@ namespace internal { namespace { - // TODO: in general we should iterate over the constraints and not over all DoFs - // for performance reasons template - void set_zero_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, unsigned int shift = 0) + void set_zero_parallel(const std::vector &cm, VEC &vec, unsigned int shift = 0) { Assert(!vec.has_ghost_elements(), ExcInternalError());//ExcGhostsPresent()); - - const unsigned int - start = vec.local_range().first, - end = vec.local_range().second; - for (unsigned int i=start; i::const_iterator it = cm.begin(); + it != cm.end(); ++it) + if (locally_owned.is_element(*it)) + vec(*it) = 0.; } - // TODO: in general we should iterate over the constraints and not over all DoFs - // for performance reasons template - void set_zero_parallel(const dealii::ConstraintMatrix &cm, parallel::distributed::Vector &vec, unsigned int shift = 0) + void set_zero_parallel(const std::vector &cm, parallel::distributed::Vector &vec, unsigned int shift = 0) { - for (unsigned int i=0; i::const_iterator it = cm.begin(); + it != cm.end(); ++it) + if (vec.in_local_range(*it)) + vec(*it) = 0.; vec.zero_out_ghosts(); } template - void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type) + void set_zero_in_parallel(const std::vector &cm, VEC &vec, internal::bool2type) { set_zero_parallel(cm, vec, 0); } // in parallel for BlockVectors template - void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type) + void set_zero_in_parallel(const std::vector &cm, VEC &vec, internal::bool2type) { unsigned int start_shift = 0; for (unsigned int j=0; j - void set_zero_serial(const dealii::ConstraintMatrix &cm, VEC &vec) + void set_zero_serial(const std::vector &cm, VEC &vec) { - // TODO would be faster: - /* std::vector::const_iterator constraint_line = cm.lines.begin(); - for (; constraint_line!=cm.lines.end(); ++constraint_line) - vec(constraint_line->line) = 0.;*/ - for (unsigned int i=0; i::const_iterator it = cm.begin(); + it != cm.end(); ++it) + vec(*it) = 0.; } template - void set_zero_all(const dealii::ConstraintMatrix &cm, VEC &vec) + void set_zero_all(const std::vector &cm, VEC &vec) { set_zero_in_parallel(cm, vec, internal::bool2type::value>()); vec.compress(VectorOperation::insert); @@ -798,13 +789,13 @@ namespace internal template - void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::Vector &vec) + void set_zero_all(const std::vector &cm, dealii::Vector &vec) { set_zero_serial(cm, vec); } template - void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::BlockVector &vec) + void set_zero_all(const std::vector &cm, dealii::BlockVector &vec) { set_zero_serial(cm, vec); } @@ -817,7 +808,12 @@ template void ConstraintMatrix::set_zero (VectorType &vec) const { - internal::ConstraintMatrix::set_zero_all(*this, vec); + // since we lines is a private member, we cannot pass it to the functions + // above. therefore, copy the content which is cheap + std::vector constrained_lines(lines.size()); + for (unsigned int i=0; i &output, const internal::bool2type /*is_block_vector*/) { + const_cast&>(vec).zero_out_ghosts(); output.reinit (locally_owned_elements, needed_elements, vec.get_mpi_communicator()); output = vec; + output.update_ghost_values(); } diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h index f8d2b168b3..1e31ba8471 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector_base.h +++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h @@ -1257,16 +1257,8 @@ namespace TrilinosWrappers else if (vector->Map().NumMyElements() > 0) { const unsigned int n_indices = vector->Map().NumMyElements(); - int * vector_indices = vector->Map().MyGlobalElements(); - unsigned int range_start = vector_indices[0]; - unsigned int range_end = range_start+1; - for (unsigned int i=1; i(vector_indices[i]) != range_end) - { - is.add_range(range_start, range_end); - range_end = range_start = vector_indices[i]; - } - is.add_range(range_start, range_end); + unsigned int * vector_indices = (unsigned int*)vector->Map().MyGlobalElements(); + is.add_indices(vector_indices, vector_indices+n_indices); is.compress(); } diff --git a/tests/mpi/step-48.cc b/tests/mpi/step-48.cc index 5d4e4e70b5..361289a12c 100644 --- a/tests/mpi/step-48.cc +++ b/tests/mpi/step-48.cc @@ -205,7 +205,7 @@ namespace Step48 void make_grid_and_dofs (); void oldstyle_operation (); void assemble_system (); - void output_results (const unsigned int timestep_number) const; + void output_results (const unsigned int timestep_number); parallel::distributed::Triangulation triangulation; FE_Q fe; @@ -301,20 +301,14 @@ namespace Step48 template void - SineGordonProblem::output_results (const unsigned int timestep_number) const + SineGordonProblem::output_results (const unsigned int timestep_number) { - parallel::distributed::Vector locally_relevant_solution; - locally_relevant_solution.reinit (dof_handler.locally_owned_dofs(), - locally_relevant_dofs, - MPI_COMM_WORLD); - locally_relevant_solution.copy_from (solution); - locally_relevant_solution.update_ghost_values(); - constraints.distribute (locally_relevant_solution); + constraints.distribute (solution); + solution.update_ghost_values(); Vector norm_per_cell (triangulation.n_active_cells()); - locally_relevant_solution.update_ghost_values(); VectorTools::integrate_difference (dof_handler, - locally_relevant_solution, + solution, ZeroFunction(), norm_per_cell, QGauss(fe_degree+1), -- 2.39.5