MGLevelGlobalTransfer/MGTransferMatrixFree: allow user to init vectors

doc/news/changes/minor/20210309Munch

@@ -0,0 +1,7 @@
+New: The class MGTransferMatrixFree::build() now also
+accepts an optional function for initializing the internal level vectors. 
+This is useful if one uses the the transfer operators in the context of 
+smoothers that are built around MatrixFree objects.
+<br>
+(Peter Munch, 2021/03/09)
+

include/deal.II/multigrid/mg_transfer.h

@@ -611,6 +611,13 @@ class MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>
   mutable MGLevelObject<LinearAlgebra::distributed::Vector<Number>>
     solution_ghosted_level_vector;
 
+  /**
+   * Function to initialize internal level vectors.
+   */
+  std::function<void(const unsigned int,
+                     LinearAlgebra::distributed::Vector<Number> &)>
+    initialize_dof_vector;
+
 private:
   /**
    * This function is called to make sure that build() has been invoked.

include/deal.II/multigrid/mg_transfer.templates.h

@@ -438,14 +438,15 @@ MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>::copy_to_mg(
         // In case a ghosted level vector has been initialized, we can simply
         // use that as a template for the vector partitioning. If not, we
         // resort to the locally owned range of the dof handler.
-        if (level <= ghosted_level_vector.max_level() &&
-            ghosted_level_vector[level].size() == dof_handler.n_dofs(level))
+        if (initialize_dof_vector)
+          initialize_dof_vector(level, dst[level]);
+        else if (level <= ghosted_level_vector.max_level() &&
+                 ghosted_level_vector[level].size() ==
+                   dof_handler.n_dofs(level))
           dst[level].reinit(ghosted_level_vector[level], false);
         else
-          {
-            dst[level].reinit(dof_handler.locally_owned_mg_dofs(level),
-                              dof_handler.get_communicator());
-          }
+          dst[level].reinit(dof_handler.locally_owned_mg_dofs(level),
+                            dof_handler.get_communicator());
       }
     else if ((perform_plain_copy == false &&
               perform_renumbered_plain_copy == false) ||
@@ -561,20 +562,29 @@ MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number>>::copy_from_mg(
     {
       // the ghosted vector should already have the correct local size (but
       // different parallel layout)
-      AssertDimension(ghosted_level_vector[level].locally_owned_size(),
-                      src[level].locally_owned_size());
+      if (ghosted_level_vector[level].size() > 0)
+        AssertDimension(ghosted_level_vector[level].locally_owned_size(),
+                        src[level].locally_owned_size());
 
       // the first time around, we copy the source vector to the temporary
       // vector that we hold for the purpose of data exchange
       LinearAlgebra::distributed::Vector<Number> &ghosted_vector =
         ghosted_level_vector[level];
-      ghosted_vector = src[level];
-      ghosted_vector.update_ghost_values();
+
+      if (ghosted_level_vector[level].size() > 0)
+        ghosted_vector = src[level];
+
+      const auto ghosted_vector_ptr = (ghosted_level_vector[level].size() > 0) ?
+                                        &ghosted_vector :
+                                        &src[level];
+
+      ghosted_vector_ptr->update_ghost_values();
+
 
       auto copy_unknowns = [&](const Table<2, unsigned int> &indices) {
         for (unsigned int i = 0; i < indices.n_cols(); ++i)
           dst.local_element(indices(0, i)) =
-            ghosted_vector.local_element(indices(1, i));
+            ghosted_vector_ptr->local_element(indices(1, i));
       };
 
       // first copy local unknowns

include/deal.II/multigrid/mg_transfer_matrix_free.h

@@ -105,6 +105,16 @@ class MGTransferMatrixFree
           &external_partitioners =
             std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>());
 
+  /**
+   * Same as above but taking a lambda for initializing vector instead of
+   * partitioners.
+   */
+  void
+  build(const DoFHandler<dim, dim> &dof_handler,
+        const std::function<void(const unsigned int,
+                                 LinearAlgebra::distributed::Vector<Number> &)>
+          &initialize_dof_vector);
+
   /**
    * Prolongate a vector from level <tt>to_level-1</tt> to level
    * <tt>to_level</tt> using the embedding matrices of the underlying finite

source/multigrid/mg_transfer_matrix_free.cc

@@ -93,6 +93,7 @@ MGTransferMatrixFree<dim, Number>::clear()
 }
 
 
+
 template <int dim, typename Number>
 void
 MGTransferMatrixFree<dim, Number>::build(
@@ -106,6 +107,20 @@ MGTransferMatrixFree<dim, Number>::build(
            "distribute_mg_dofs() function called, but this is a prerequisite "
            "for multigrid transfers. You will need to call this function, "
            "probably close to where you already call distribute_dofs()."));
+  if (external_partitioners.size() > 0)
+    {
+      Assert(
+        this->initialize_dof_vector == nullptr,
+        ExcMessage(
+          "A initialize_dof_vector function has already been registered in the constructor!"));
+
+      this->initialize_dof_vector =
+        [external_partitioners](
+          const unsigned int                          level,
+          LinearAlgebra::distributed::Vector<Number> &vec) {
+          vec.reinit(external_partitioners[level]);
+        };
+    }
 
   this->fill_and_communicate_copy_indices(dof_handler);
 
@@ -186,6 +201,39 @@ MGTransferMatrixFree<dim, Number>::build(
 
 
 
+template <int dim, typename Number>
+void
+MGTransferMatrixFree<dim, Number>::build(
+  const DoFHandler<dim, dim> &dof_handler,
+  const std::function<void(const unsigned int,
+                           LinearAlgebra::distributed::Vector<Number> &)>
+    &initialize_dof_vector)
+{
+  if (initialize_dof_vector)
+    {
+      const unsigned int n_levels =
+        dof_handler.get_triangulation().n_global_levels();
+
+      std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+        external_partitioners(n_levels);
+
+      for (unsigned int level = 0; level < n_levels; ++level)
+        {
+          LinearAlgebra::distributed::Vector<Number> vector;
+          initialize_dof_vector(level, vector);
+          external_partitioners[level] = vector.get_partitioner();
+        }
+
+      build(dof_handler, external_partitioners);
+    }
+  else
+    {
+      build(dof_handler);
+    }
+}
+
+
+
 template <int dim, typename Number>
 void
 MGTransferMatrixFree<dim, Number>::prolongate(

tests/matrix_free/multigrid_dg_sip_01.cc

@@ -596,17 +596,19 @@ do_test(const DoFHandler<dim> &dof, const unsigned n_q_points_1d)
   mg_constrained_dofs.initialize(dof);
   mg_constrained_dofs.make_zero_boundary_constraints(dof, {0});
 
-  MGTransferMF<dim, LevelMatrixType> mg_transfer(mg_matrices,
-                                                 mg_constrained_dofs);
-  mg_transfer.build(dof);
+  MGTransferMatrixFree<dim, number> mg_transfer(mg_constrained_dofs);
+
+  mg_transfer.build(dof, [&](const unsigned int level, auto &vec) {
+    mg_matrices[level].initialize_dof_vector(vec);
+  });
 
   mg::Matrix<LinearAlgebra::distributed::Vector<double>> mg_matrix(mg_matrices);
 
   Multigrid<LinearAlgebra::distributed::Vector<double>> mg(
     mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
   PreconditionMG<dim,
                  LinearAlgebra::distributed::Vector<double>,
-                 MGTransferMF<dim, LevelMatrixType>>
+                 MGTransferMatrixFree<dim, number>>
     preconditioner(dof, mg, mg_transfer);
 
   {