]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make distribute() work properly with parallel::distributed::Vector. Also cleanup...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 31 May 2013 08:14:35 +0000 (08:14 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 31 May 2013 08:14:35 +0000 (08:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@29678 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-48/step-48.cc
deal.II/include/deal.II/lac/constraint_matrix.templates.h
deal.II/include/deal.II/lac/trilinos_vector_base.h
tests/mpi/step-48.cc

index 2b3f7a651a622cda9c6f22b6969a72308d368d57..9e2f4a24664b4572b858eafbf4b56b27ea21534d 100644 (file)
@@ -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<dim>   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 <int dim>
   void
-  SineGordonProblem<dim>::output_results (const unsigned int timestep_number) const
+  SineGordonProblem<dim>::output_results (const unsigned int timestep_number)
   {
-    parallel::distributed::Vector<double> 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<float> 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<dim>(),
                                        norm_per_cell,
                                        QGauss<dim>(fe_degree+1),
@@ -464,7 +466,7 @@ namespace Step48
     DataOut<dim> 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 <int dim>
   void
   SineGordonProblem<dim>::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
-    // <code>old_solution</code> with <code>old_old_solution</code>, which
-    // means that <code>old_old_solution</code> gets
-    // <code>old_solution</code>, which is what we expect. Similarly,
-    // <code>old_solution</code> gets the content from <code>solution</code>
-    // in the next step. Afterward, <code>solution</code> holds
-    // <code>old_old_solution</code>, 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 <code>old_solution</code> with
+    // <code>old_old_solution</code>, which means that
+    // <code>old_old_solution</code> gets <code>old_solution</code>, which is
+    // what we expect. Similarly, <code>old_solution</code> gets the content
+    // from <code>solution</code> in the next step. Afterward,
+    // <code>solution</code> holds <code>old_old_solution</code>, but that
+    // will be overwritten during this step.
     unsigned int timestep_number = 1;
 
     Timer timer;
index 5758495937ecac4fab9dbb6a043ed161ad6e74c3..ed19a1c0d63167a2923e062feeb1910ca54f82b0 100644 (file)
@@ -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<class VEC>
-      void set_zero_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, unsigned int shift = 0)
+      void set_zero_parallel(const std::vector<unsigned int> &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<end; ++i)
-          if (cm.is_constrained (shift + i))
-            vec(i) = 0;
+        IndexSet locally_owned = vec.locally_owned_elements();
+        for (typename std::vector<unsigned int>::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<typename Number>
-      void set_zero_parallel(const dealii::ConstraintMatrix &cm, parallel::distributed::Vector<Number> &vec, unsigned int shift = 0)
+      void set_zero_parallel(const std::vector<unsigned int> &cm, parallel::distributed::Vector<Number> &vec, unsigned int shift = 0)
       {
-        for (unsigned int i=0; i<vec.local_size(); ++i)
-          if (cm.is_constrained (shift + vec.local_range().first+i))
-            vec.local_element(i) = 0;
+        for (typename std::vector<unsigned int>::const_iterator it = cm.begin();
+             it != cm.end(); ++it)
+          if (vec.in_local_range(*it))
+            vec(*it) = 0.;
         vec.zero_out_ghosts();
       }
 
       template<class VEC>
-      void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type<false>)
+      void set_zero_in_parallel(const std::vector<unsigned int> &cm, VEC &vec, internal::bool2type<false>)
       {
         set_zero_parallel(cm, vec, 0);
       }
 
       // in parallel for BlockVectors
       template<class VEC>
-      void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type<true>)
+      void set_zero_in_parallel(const std::vector<unsigned int> &cm, VEC &vec, internal::bool2type<true>)
       {
         unsigned int start_shift = 0;
         for (unsigned int j=0; j<vec.n_blocks(); ++j)
@@ -778,19 +773,15 @@ namespace internal
       }
 
       template<class VEC>
-      void set_zero_serial(const dealii::ConstraintMatrix &cm, VEC &vec)
+      void set_zero_serial(const std::vector<unsigned int> &cm, VEC &vec)
       {
-        // TODO would be faster:
-        /* std::vector<dealii::ConstraintMatrix::ConstraintLine>::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<vec.size(); ++i)
-          if (cm.is_constrained (i))
-            vec(i) = 0;
+        for (typename std::vector<unsigned int>::const_iterator it = cm.begin();
+             it != cm.end(); ++it)
+           vec(*it) = 0.;
       }
 
       template<class VEC>
-      void set_zero_all(const dealii::ConstraintMatrix &cm, VEC &vec)
+      void set_zero_all(const std::vector<unsigned int> &cm, VEC &vec)
       {
         set_zero_in_parallel<VEC>(cm, vec, internal::bool2type<IsBlockVector<VEC>::value>());
         vec.compress(VectorOperation::insert);
@@ -798,13 +789,13 @@ namespace internal
 
 
       template<class T>
-      void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::Vector<T> &vec)
+      void set_zero_all(const std::vector<unsigned int> &cm, dealii::Vector<T> &vec)
       {
         set_zero_serial(cm, vec);
       }
 
       template<class T>
-      void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::BlockVector<T> &vec)
+      void set_zero_all(const std::vector<unsigned int> &cm, dealii::BlockVector<T> &vec)
       {
         set_zero_serial(cm, vec);
       }
@@ -817,7 +808,12 @@ template <class VectorType>
 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<unsigned int> constrained_lines(lines.size());
+  for (unsigned int i=0; i<lines.size(); ++i)
+    constrained_lines[i] = lines[i].line;
+  internal::ConstraintMatrix::set_zero_all(constrained_lines, vec);
 }
 
 
@@ -943,7 +939,7 @@ ConstraintMatrix::distribute (const VectorType &condensed,
               for (unsigned int i=row+1; i<n_rows; ++i)
                 old_line.push_back (i-shift);
               break;
-            };
+            }
         }
       else
         old_line.push_back (row-shift);
@@ -1020,8 +1016,10 @@ namespace internal
                                        parallel::distributed::Vector<number>       &output,
                                        const internal::bool2type<false>             /*is_block_vector*/)
     {
+      const_cast<parallel::distributed::Vector<number>&>(vec).zero_out_ghosts();
       output.reinit (locally_owned_elements, needed_elements, vec.get_mpi_communicator());
       output = vec;
+      output.update_ghost_values();
     }
 
 
index f8d2b168b387e682011263d9e2f530c26bedaf0a..1e31ba84719b66c0505718fb1f5924e162f1cb0b 100644 (file)
@@ -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<n_indices; ++i, ++range_end)
-          if (static_cast<unsigned int>(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();
       }
 
index 5d4e4e70b5c28a60f495d4990c8efe92daba6cef..361289a12c67909004fa1f03116252be2607c79b 100644 (file)
@@ -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<dim>   triangulation;
     FE_Q<dim>            fe;
@@ -301,20 +301,14 @@ namespace Step48
 
   template <int dim>
   void
-  SineGordonProblem<dim>::output_results (const unsigned int timestep_number) const
+  SineGordonProblem<dim>::output_results (const unsigned int timestep_number)
   {
-    parallel::distributed::Vector<double> 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<float> norm_per_cell (triangulation.n_active_cells());
-    locally_relevant_solution.update_ghost_values();
     VectorTools::integrate_difference (dof_handler,
-                                       locally_relevant_solution,
+                                       solution,
                                        ZeroFunction<dim>(),
                                        norm_per_cell,
                                        QGauss<dim>(fe_degree+1),

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.