Make GMG default preconditioner

Author: gassmoeller

benchmarks/newton_solver_benchmark_set/tosi_et_al_2015/input.prm

@@ -28,6 +28,7 @@ subsection Solver parameters
   subsection Stokes solver parameters
     set Number of cheap Stokes solver steps = 200
     set Linear solver tolerance = 1e-6
+    set Stokes solver type = block AMG
   end
 
   subsection Newton solver parameters
@@ -100,6 +101,7 @@ end
 
 subsection Material model
   set Model name = TosiMaterial
+  set Material averaging =  none
 
   subsection Tosi benchmark
     set Reference density             = 1

doc/modules/changes/20240531_gassmoeller

@@ -0,0 +1,10 @@
+Changed: ASPECT's default Stokes preconditioner has been changed from
+a block AMG preconditioner to the geometric multigrid (GMG) preconditioner
+described in Clevenger and Heister, 2021 (https://doi.org/10.1002/nla.2375).
+The GMG preconditioner generally performs better, but requires that
+the viscosity is averaged using one of the material model averaging functions.
+If models use features that are not supported by the GMG preconditioner
+and the GMG preconditioner has not been explicitly set ASPECT will
+fall back to the AMG preconditioner.
+<br>
+(Rene Gassmoeller, 2024/06/09)

include/aspect/material_model/interface.h

@@ -712,7 +712,8 @@ namespace aspect
         log_average,
         harmonic_average_only_viscosity,
         geometric_average_only_viscosity,
-        project_to_Q1_only_viscosity
+        project_to_Q1_only_viscosity,
+        solver_default
       };
 
 

include/aspect/parameters.h

@@ -323,12 +323,13 @@ namespace aspect
       {
         block_amg,
         direct_solver,
-        block_gmg
+        block_gmg,
+        default_solver
       };
 
       static const std::string pattern()
       {
-        return "block AMG|direct solver|block GMG";
+        return "default solver|block AMG|direct solver|block GMG";
       }
 
       static Kind
@@ -340,6 +341,8 @@ namespace aspect
           return direct_solver;
         else if (input == "block GMG")
           return block_gmg;
+        else if (input == "default solver")
+          return default_solver;
         else
           AssertThrow(false, ExcNotImplemented());
 

include/aspect/simulator.h

@@ -1740,6 +1740,14 @@ namespace aspect
       void
       check_consistency_of_formulation ();
 
+      /**
+      * This function checks if the default solver and/or material
+      * averaging were selected and if so, determines the appropriate
+      * solver and/or averaging option.
+      */
+      void
+      select_default_solver_and_averaging ();
+
       /**
        * This function checks that the user-selected boundary conditions do not
        * contain contradictions. If an incorrect selection is detected it

source/material_model/interface.cc

@@ -430,7 +430,7 @@ namespace aspect
     {
       std::string get_averaging_operation_names ()
       {
-        return "none|arithmetic average|harmonic average|geometric average|pick largest|project to Q1|log average|harmonic average only viscosity|geometric average only viscosity|project to Q1 only viscosity";
+        return "none|solver default|arithmetic average|harmonic average|geometric average|pick largest|project to Q1|log average|harmonic average only viscosity|geometric average only viscosity|project to Q1 only viscosity";
       }
 
 
@@ -456,6 +456,8 @@ namespace aspect
           return geometric_average_only_viscosity;
         else if (s == "project to Q1 only viscosity")
           return project_to_Q1_only_viscosity;
+        else if (s == "solver default")
+          return solver_default;
         else
           AssertThrow (false,
                        ExcMessage ("The value <" + s + "> for a material "

source/simulator/core.cc

@@ -212,7 +212,8 @@ namespace aspect
                       Triangulation<dim>::do_not_produce_unrefined_islands)
                    )
                    ,
-                   (parameters.stokes_solver_type == Parameters<dim>::StokesSolverType::block_gmg
+                   (parameters.stokes_solver_type == Parameters<dim>::StokesSolverType::block_gmg ||
+                    parameters.stokes_solver_type == Parameters<dim>::StokesSolverType::default_solver
                     ?
                     static_cast<typename parallel::distributed::Triangulation<dim>::Settings>
                     (parallel::distributed::Triangulation<dim>::mesh_reconstruction_after_repartitioning |
@@ -423,6 +424,9 @@ namespace aspect
         newton_handler->parameters.parse_parameters(prm);
       }
 
+    // choose the default solver and averaging scheme
+    select_default_solver_and_averaging();
+
     if (parameters.stokes_solver_type == Parameters<dim>::StokesSolverType::block_gmg)
       {
         switch (parameters.stokes_velocity_degree)

source/simulator/helper_functions.cc

@@ -26,6 +26,7 @@
 #include <aspect/global.h>
 
 #include <aspect/geometry_model/interface.h>
+#include <aspect/geometry_model/ellipsoidal_chunk.h>
 #include <aspect/heating_model/interface.h>
 #include <aspect/heating_model/adiabatic_heating.h>
 #include <aspect/material_model/interface.h>
@@ -2131,6 +2132,7 @@ namespace aspect
   }
 
 
+
   template <int dim>
   void
   Simulator<dim>::replace_outflow_boundary_ids(const unsigned int offset)
@@ -2533,7 +2535,53 @@ namespace aspect
       }
     return new_linear_stokes_solver_tolerance;
   }
+
+
+
+  template <int dim>
+  void
+  Simulator<dim>::select_default_solver_and_averaging()
+  {
+    if (parameters.stokes_solver_type == Parameters<dim>::StokesSolverType::default_solver)
+      {
+        // Catch all situations that are not supported by the GMG solver:
+        //   - Melt transport
+        //   - Ellipsoidal geometry
+        //   - Locally conservative discretization
+        //   - Implicit reference density profile
+        //   - Periodic boundaries
+        //   - Stokes velocity degree not 2 or 3
+        //   - Material averaging explicitly disabled
+        if (parameters.include_melt_transport == true ||
+            dynamic_cast<const GeometryModel::EllipsoidalChunk<dim>*>(geometry_model.get()) != nullptr ||
+            parameters.use_locally_conservative_discretization == true ||
+            (material_model->is_compressible() == true && parameters.formulation_mass_conservation ==
+             Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile) ||
+            (geometry_model->get_periodic_boundary_pairs().size()) > 0 ||
+            (parameters.stokes_velocity_degree < 2 || parameters.stokes_velocity_degree > 3) ||
+            parameters.material_averaging == MaterialModel::MaterialAveraging::none)
+          {
+            // GMG is not supported (yet), by default fall back to AMG.
+            parameters.stokes_solver_type = Parameters<dim>::StokesSolverType::block_amg;
+          }
+        else
+          {
+            // GMG is supported for all other cases
+            parameters.stokes_solver_type = Parameters<dim>::StokesSolverType::block_gmg;
+          }
+      }
+
+    // Now pick an appropriate material averaging for the chosen solver
+    if (parameters.material_averaging == MaterialModel::MaterialAveraging::solver_default)
+      {
+        if (parameters.stokes_solver_type == Parameters<dim>::StokesSolverType::block_gmg)
+          parameters.material_averaging = MaterialModel::MaterialAveraging::project_to_Q1_only_viscosity;
+        else
+          parameters.material_averaging = MaterialModel::MaterialAveraging::none;
+      }
+  }
 }
+
 // explicit instantiation of the functions we implement in this file
 namespace aspect
 {
@@ -2573,7 +2621,9 @@ namespace aspect
                                                                             const double linear_stokes_solver_tolerance, \
                                                                             const double stokes_residual, \
                                                                             const double newton_residual, \
-                                                                            const double newton_residual_old);
+                                                                            const double newton_residual_old); \
+  template void Simulator<dim>::select_default_solver_and_averaging();
+
 
   ASPECT_INSTANTIATE(INSTANTIATE)
 

source/simulator/parameters.cc

@@ -376,7 +376,7 @@ namespace aspect
 
       prm.enter_subsection ("Stokes solver parameters");
       {
-        prm.declare_entry ("Stokes solver type", "block AMG",
+        prm.declare_entry ("Stokes solver type", "default solver",
                            Patterns::Selection(StokesSolverType::pattern()),
                            "This is the type of solver used on the Stokes system. The block geometric "
                            "multigrid solver currently has a limited implementation and therefore "
@@ -1342,7 +1342,7 @@ namespace aspect
     // preconditioner
     prm.enter_subsection ("Material model");
     {
-      prm.declare_entry ("Material averaging", "none",
+      prm.declare_entry ("Material averaging", "solver default",
                          Patterns::Selection(MaterialModel::MaterialAveraging::
                                              get_averaging_operation_names()),
                          "Whether or not (and in the first case, how) to do any averaging of "