address comments

Author: tjhei

doc/sphinx/references.bib

@@ -12202,6 +12202,20 @@ @article{vankeken:etal:2008
 doi = {10.1016/j.pepi.2008.04.015}
 }
 
+@Article{munch:globalcoarsening:2023,
+  author    = {Peter Munch and Timo Heister and Laura Prieto Saavedra and Martin Kronbichler},
+  title     = {Efficient Distributed Matrix-free Multigrid Methods on Locally Refined Meshes for {FEM} Computations},
+  journal   = {{ACM} Transactions on Parallel Computing},
+  year      = {2023},
+  volume    = {10},
+  number    = {1},
+  pages     = {1--38},
+  month     = mar,
+  doi       = {10.1145/3580314},
+  publisher = {Association for Computing Machinery ({ACM})},
+  url       = {https://arxiv.org/abs/2203.12292},
+}
+
 @article{mulyukova:bercovici:2018,
   title={Collapse of passive margins by lithospheric damage and plunging grain size},
   author={Mulyukova, Elvira and Bercovici, David},

include/aspect/simulator/solver/matrix_free_operators.h

@@ -358,14 +358,16 @@ namespace aspect
         void clear () override;
 
         /**
-        *
+         * Initialize the MatrixFree object given in @p mf_storage and use that to
+         * initialize this operator.
          */
         void reinit(const Mapping<dim>              &mapping,
                     const DoFHandler<dim>           &dof_handler_v,
                     const DoFHandler<dim>           &dof_handler_p,
                     const AffineConstraints<number> &constraints_v,
                     const AffineConstraints<number> &constraints_p,
-                    std::shared_ptr<MatrixFree<dim,double>> mf_storage);
+                    std::shared_ptr<MatrixFree<dim,double>> mf_storage,
+                    const unsigned int level = numbers::invalid_unsigned_int);
 
         /**
          * Pass in a reference to the problem data.
@@ -432,13 +434,16 @@ namespace aspect
         void clear () override;
 
         /**
-         *
+         * Initialize the MatrixFree object given in @p mf_storage and use that to
+         * initialize this operator.
          */
         void reinit(const Mapping<dim>              &mapping,
-                    const DoFHandler<dim>           &dof_handler,
-                    const DoFHandler<dim>           &dof_handler_other,
-                    const AffineConstraints<number> &constraints,
-                    std::shared_ptr<MatrixFree<dim,double>> mf_storage);
+                    const DoFHandler<dim>           &dof_handler_v,
+                    const DoFHandler<dim>           &dof_handler_p,
+                    const AffineConstraints<number> &constraints_v,
+                    const AffineConstraints<number> &constraints_p,
+                    std::shared_ptr<MatrixFree<dim,double>> mf_storage,
+                    const unsigned int level = numbers::invalid_unsigned_int);
         /**
          * Pass in a reference to the problem data.
          */

source/simulator/core.cc

@@ -125,7 +125,7 @@ namespace aspect
 
     template <int dim>
     typename Triangulation<dim>::MeshSmoothing
-    smoothing_flags(bool global_coarsening)
+    smoothing_flags(const bool global_coarsening)
     {
       if (global_coarsening)
         return Triangulation<dim>::limit_level_difference_at_vertices;

source/simulator/parameters.cc

@@ -435,7 +435,10 @@ namespace aspect
 
         prm.declare_entry ("Stokes GMG type", "local smoothing",
                            Patterns::Selection(StokesGMGType::pattern()),
-                           "local smoothing or global coarsening");
+                           "The choice of geometric multigrid, either 'local smoothing' (the default) "
+                           " or 'global coarsening'. Local smoothing (\\cite{clevenger:heister:2021}) "
+                           "has been extensively tested and works in many more situations, while "
+                           "global coarsening is shown to be up to 3x faster (\\cite{munch:globalcoarsening:2023}).");
 
         prm.declare_entry ("Use direct solver for Stokes system", "false",
                            Patterns::Bool(),
@@ -1117,7 +1120,7 @@ namespace aspect
         prm.declare_entry ("Stabilization method", "entropy viscosity",
                            Patterns::Selection("entropy viscosity|SUPG"),
                            "Select the method for stabilizing the advection equation. The original "
-                           "method implemented is 'entropy viscosity' as described in \\cite {kronbichler:etal:2012}. "
+                           "method implemented is 'entropy viscosity' as described in \\cite{kronbichler:etal:2012}. "
                            "SUPG is currently experimental.");
 
         prm.declare_entry ("List of compositional fields with disabled boundary entropy viscosity", "",

source/simulator/solver/matrix_free_operators.cc

@@ -483,12 +483,13 @@ namespace aspect
                                                                              const DoFHandler<dim>           &dof_handler_p,
                                                                              const AffineConstraints<number> &constraints_v,
                                                                              const AffineConstraints<number> &constraints_p,
-                                                                             std::shared_ptr<MatrixFree<dim,double>> mf_storage)
+                                                                             std::shared_ptr<MatrixFree<dim,double>> mf_storage,
+                                                                             const unsigned int level)
   {
     typename MatrixFree<dim, number>::AdditionalData data;
     data.mapping_update_flags =
       update_quadrature_points /*| update_gradients*/ | update_values;
-    data.mg_level = numbers::invalid_unsigned_int;
+    data.mg_level = level;
 
     data.tasks_parallel_scheme =
       MatrixFree<dim,double>::AdditionalData::none;
@@ -663,15 +664,16 @@ namespace aspect
   template <int dim, int degree_v, typename number>
   void
   MatrixFreeStokesOperators::ABlockOperator<dim,degree_v,number>::reinit(const Mapping<dim>              &mapping,
-                                                                         const DoFHandler<dim>           &dof_handler,
-                                                                         const DoFHandler<dim>           &dof_handler_other,
-                                                                         const AffineConstraints<number> &constraints,
-                                                                         std::shared_ptr<MatrixFree<dim,double>> mf_storage)
+                                                                         const DoFHandler<dim>           &dof_handler_v,
+                                                                         const DoFHandler<dim>           &dof_handler_p,
+                                                                         const AffineConstraints<number> &constraints_v,
+                                                                         const AffineConstraints<number> &constraints_p,
+                                                                         std::shared_ptr<MatrixFree<dim,double>> mf_storage,
+                                                                         const unsigned int level)
   {
     typename MatrixFree<dim, number>::AdditionalData data;
-    data.mapping_update_flags =
-      update_quadrature_points | update_gradients | update_values;
-    data.mg_level = numbers::invalid_unsigned_int;
+    data.mapping_update_flags = update_quadrature_points | update_gradients | update_values;
+    data.mg_level = level;
 
     typename MatrixFree<dim,double>::AdditionalData additional_data;
     additional_data.tasks_parallel_scheme =
@@ -681,8 +683,8 @@ namespace aspect
 
     AffineConstraints<number> dummy;
     mf_storage->reinit(mapping,
-    std::vector< const DoFHandler< dim > *> {&dof_handler, &dof_handler_other},
-    std::vector< const AffineConstraints< number > *> {&constraints, &dummy} ,
+    std::vector< const DoFHandler< dim > *> {&dof_handler_v, &dof_handler_p},
+    std::vector< const AffineConstraints< number > *> {&constraints_v, &constraints_p} ,
     QGauss<1>(degree_v+1), additional_data);
 
     this->initialize(mf_storage, std::vector< unsigned int > {0}, std::vector< unsigned int > {0});

source/simulator/solver/stokes_matrix_free_global_coarsening.cc

@@ -145,7 +145,7 @@ namespace aspect
   std::string
   StokesMatrixFreeHandlerGlobalCoarseningImplementation<dim, velocity_degree>::name () const
   {
-    return "GC-GMG";
+    return "GMG-GC";
   }
 
 
@@ -173,11 +173,10 @@ namespace aspect
   template <int dim, int velocity_degree>
   void StokesMatrixFreeHandlerGlobalCoarseningImplementation<dim, velocity_degree>::evaluate_material_model ()
   {
-    sim.pcout << "evaluate_material_model()" << std::endl;
     const auto &dof_handler_projection = dofhandlers_projection.back();
 
     dealii::LinearAlgebra::distributed::Vector<double> active_viscosity_vector(dof_handler_projection.locally_owned_dofs(),
-                                                                               this->get_triangulation().get_communicator());
+                                                                               this->get_mpi_communicator());
 
     const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.velocities;
 
@@ -248,8 +247,8 @@ namespace aspect
       active_viscosity_vector.compress(VectorOperation::insert);
     }
 
-    minimum_viscosity = dealii::Utilities::MPI::min(minimum_viscosity_local, this->get_triangulation().get_communicator());
-    maximum_viscosity = dealii::Utilities::MPI::max(maximum_viscosity_local, this->get_triangulation().get_communicator());
+    minimum_viscosity = dealii::Utilities::MPI::min(minimum_viscosity_local, this->get_mpi_communicator());
+    maximum_viscosity = dealii::Utilities::MPI::max(maximum_viscosity_local, this->get_mpi_communicator());
 
     FEValues<dim> fe_values_projection (this->get_mapping(),
                                         fe_projection,
@@ -369,7 +368,7 @@ namespace aspect
         AffineConstraints<double> cs;
         std::shared_ptr<MatrixFree<dim,double>>
         mf(new MatrixFree<dim,double>());
-        mf->reinit(*sim.mapping, dofhandlers_projection[l], cs, QGauss<1>(degree+1));
+        mf->reinit(this->get_mapping(), dofhandlers_projection[l], cs, QGauss<1>(degree+1));
         temp_ops[l].initialize(mf);
       }
 
@@ -380,9 +379,9 @@ namespace aspect
       temp_ops[l].initialize_dof_vector(vec);
     });
 
-    transfer.template interpolate_to_mg(dof_handler_projection,
-                                        level_viscosity_vector,
-                                        active_viscosity_vector);
+    transfer.interpolate_to_mg(dof_handler_projection,
+                               level_viscosity_vector,
+                               active_viscosity_vector);
 
     for (unsigned int level=min_level; level<=max_level; ++level)
       {
@@ -1454,14 +1453,6 @@ namespace aspect
     if (solve_newton_system == false)
       outputs.pressure_normalization_adjustment = this->normalize_pressure(solution_vector);
 
-    // convert melt pressures
-    // TODO: We assert in the StokesMatrixFreeHandler constructor that we
-    //       are not including melt transport.
-    // TODO: We have to const_cast, because the compute_melt_variables
-    //       function assembles and solves an equation. To avoid this we either
-    //       have to hand over non-const references to the current function, or
-    //       call the compute_melt_variables function after finishing the current function.
-
     return outputs;
   }
 
@@ -1489,16 +1480,16 @@ namespace aspect
                 }
             }
 
-      have_periodic_hanging_nodes = (dealii::Utilities::MPI::max(have_periodic_hanging_nodes ? 1 : 0, this->get_triangulation().get_communicator())) == 1;
+      have_periodic_hanging_nodes = (dealii::Utilities::MPI::max(have_periodic_hanging_nodes ? 1 : 0, this->get_mpi_communicator())) == 1;
       AssertThrow(have_periodic_hanging_nodes==false, ExcNotImplemented());
     }
 
-    const Mapping<dim> &mapping = *sim.mapping;
+    const Mapping<dim> &mapping = this->get_mapping();
 
     // This vector will be refilled with the new MatrixFree objects below:
     matrix_free_objects.clear();
 
-    trias = dealii::MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence (sim.triangulation);
+    trias = dealii::MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence (this->get_triangulation());
     min_level = 0;
     max_level = trias.size() - 1;
     constraints_v.resize(min_level, max_level);
@@ -1540,38 +1531,57 @@ namespace aspect
             constraint.reinit(locally_relevant_dofs);
 #endif
 
-            sim.compute_initial_velocity_boundary_constraints(constraint); // TODO: this can't work
-            sim.compute_current_velocity_boundary_constraints(constraint); // TODO: this can't work
 
-            for (const auto &p : sim.boundary_velocity_manager.get_zero_boundary_velocity_indicators())
+            std::set<types::boundary_id> dirichlet_boundary = this->get_boundary_velocity_manager().get_zero_boundary_velocity_indicators();
+            for (const auto boundary_id: this->get_boundary_velocity_manager().get_prescribed_boundary_velocity_indicators())
+              {
+                const ComponentMask component_mask = this->get_boundary_velocity_manager().get_component_mask(boundary_id);
+
+                if (component_mask != ComponentMask(this->introspection().n_components, false))
+                  {
+                    ComponentMask velocity_mask(fe_v.n_components(), false);
+
+                    for (unsigned int i=0; i<dim; ++i)
+                      velocity_mask.set(i, component_mask[this->introspection().component_indices.velocities[i]]);
+
+                    VectorTools::interpolate_boundary_values (mapping,
+                                                              dof_handler,
+                                                              boundary_id,
+                                                              Functions::ZeroFunction<dim>(dim),
+                                                              constraint,
+                                                              velocity_mask);
+                  }
+                else
+                  {
+                    // no mask given: add at the end
+                    dirichlet_boundary.insert(boundary_id);
+                  }
+              }
+
+
+            for (const auto id : dirichlet_boundary)
               {
                 VectorTools::interpolate_boundary_values (mapping,
                                                           dof_handler,
-                                                          p,
+                                                          id,
                                                           Functions::ZeroFunction<dim>(dim),
                                                           constraint);
-
               }
 
 
-//            VectorTools::interpolate_boundary_values(
-//                mapping,
-//                dof_handler,
-//                0,
-//                Functions::ZeroFunction<dim, typename LevelOperatorType::value_type>(
-//                    dof_handler_in.get_fe().n_components()),
-//                constraint);
-            Assert(sim.boundary_velocity_manager.get_active_plugins().size() == 0,
+            Assert(this->get_boundary_velocity_manager().get_active_plugins().size() == 0,
                    ExcNotImplemented());
 
-            Assert(sim.boundary_velocity_manager.get_tangential_boundary_velocity_indicators().size() == 0,
-                   ExcNotImplemented());
-//            VectorTools::compute_no_normal_flux_constraints (dof_handler_v,
-//                                                            /* first_vector_component= */
-//                                                            0,
-//                                                            sim.boundary_velocity_manager.get_tangential_boundary_velocity_indicators(),
-//                                                            constraints_v,
-//                                                            *sim.mapping);
+            if (this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators().size() > 0)
+              {
+                Assert(false, dealii::StandardExceptions::ExcNotImplemented("This is not tested."));
+                VectorTools::compute_no_normal_flux_constraints (dof_handler,
+                                                                 0 /* first_vector_component */,
+                                                                 this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators(),
+                                                                 constraint,
+                                                                 mapping);
+
+              }
 
             DoFTools::make_hanging_node_constraints(dof_handler, constraint);
             constraint.close();
@@ -1607,7 +1617,7 @@ namespace aspect
 
           matrix_free = std::make_shared<MatrixFree<dim,double>>();
           matrix_free_objects.push_back(matrix_free);
-          mg_matrices_A_block[l].reinit(mapping, dofhandlers_v[l], dofhandlers_p[l], constraints_v[l], matrix_free);
+          mg_matrices_A_block[l].reinit(mapping, dofhandlers_v[l], dofhandlers_p[l], constraints_v[l], constraints_p[l], matrix_free);
 
           matrix_free = std::make_shared<MatrixFree<dim,double>>();
           matrix_free_objects.push_back(matrix_free);
@@ -1625,52 +1635,6 @@ namespace aspect
         }
     }
 
-
-
-    // Multigrid DoF setup
-#if 0
-    {
-      //Ablock GMG
-      dof_handler_v.distribute_mg_dofs();
-
-      mg_constrained_dofs_A_block.clear();
-      mg_constrained_dofs_A_block.initialize(dof_handler_v);
-
-      std::set<types::boundary_id> dirichlet_boundary = this->get_boundary_velocity_manager().get_zero_boundary_velocity_indicators();
-      for (const auto boundary_id: this->get_boundary_velocity_manager().get_prescribed_boundary_velocity_indicators())
-        {
-          const ComponentMask component_mask = this->get_boundary_velocity_manager().get_component_mask(boundary_id);
-
-          if (component_mask != ComponentMask(this->introspection().n_components, false))
-            {
-              ComponentMask velocity_mask(fe_v.n_components(), false);
-
-              for (unsigned int i=0; i<dim; ++i)
-                velocity_mask.set(i, component_mask[this->introspection().component_indices.velocities[i]]);
-
-              mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, {boundary_id}, velocity_mask);
-            }
-          else
-            {
-              // no mask given: add at the end
-              dirichlet_boundary.insert(boundary_id);
-            }
-        }
-
-      // Unconditionally call this function, even if the set is empty. Otherwise, the data structure
-      // for boundary indices will not be created (if mesh has no Dirichlet conditions).
-      mg_constrained_dofs_A_block.make_zero_boundary_constraints(dof_handler_v, dirichlet_boundary);
-
-      //Schur complement matrix GMG
-      dof_handler_p.distribute_mg_dofs();
-
-      mg_constrained_dofs_Schur_complement.clear();
-      mg_constrained_dofs_Schur_complement.initialize(dof_handler_p);
-
-      dof_handler_projection.distribute_mg_dofs();
-    }
-#endif
-
     // Setup the matrix-free operators
     std::shared_ptr<MatrixFree<dim,double>> matrix_free = std::make_shared<MatrixFree<dim,double>>();
     matrix_free_objects.push_back(matrix_free);
@@ -1758,13 +1722,14 @@ namespace aspect
   template <int dim, int velocity_degree>
   void StokesMatrixFreeHandlerGlobalCoarseningImplementation<dim, velocity_degree>::build_preconditioner()
   {
-    TimerOutput::Scope timer (this->get_computing_timer(), "Build Stokes preconditioner");
+    this->get_computing_timer().enter_subsection("Build Stokes preconditioner");
 
     for (auto l = min_level; l <= max_level; ++l)
       {
         mg_matrices_Schur_complement[l].compute_diagonal();
         mg_matrices_A_block[l].compute_diagonal();
       }
+    this->get_computing_timer().leave_subsection("Build Stokes preconditioner");
   }