Automatically select memory-efficient PtAP algorithms for HMG

framework/doc/content/source/preconditioners/HMGPreconditioner.md

@@ -0,0 +1,65 @@
+# HMG
+
+## Overview
+
+HMG stands for High-performance (Hybrid) MultiGrid method.
+[HMG](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCHMG.html)'s
+essential idea is to separate a multigrid method into two steps;
+matrix coarsening and level solvers. The main motivation of HMG in
+the first place is to use
+[HYPRE](https://computing.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods)
+to do matrix coarsening and generate interpolations
+and use [PETSc](https://www.mcs.anl.gov/petsc/) preconditioners/solvers as level solvers.
+However, the code is more general, and it can use other
+codes such as [GAMG](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCGAMG.html)
+or a user code to generate interpolations.
+
+
+## Example 1
+
+```
+[Executioner]
+  type = Steady
+
+  solve_type = 'PJFNK'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'hmg'
+[]
+```
+
+This configuration uses HYPRE to generate interpolations
+and uses SOR preconditioners from PETSc as level solvers.
+
+## Example 2
+
+```
+[Executioner]
+  type = Steady
+
+  solve_type = 'PJFNK'
+
+  petsc_options_iname = '-pc_type -pc_hmg_use_subspace_coarsening'
+  petsc_options_value = 'hmg true'
+  petsc_options = '-snes_view'
+[]
+```
+
+If there are multiple nonlinear variables, this configuration will reuse interpolations
+generated for the first nonlinear variable for all other variables. This will significantly
+speedup preconditioner setup. A complete toy example is [test/tests/preconditioners/hmg/diffusion_hmg.i].
+We also demonstrate this capability for
+realistic neutron transport calculations in the following paper:
+
+```tex
+@article{kong2020highly,
+  title={{A Highly Parallel Multilevel Newton--Krylov--Schwarz Method with Subspace-Based Coarsening and Partition-Based Balancing for the Multigroup Neutron Transport Equation on Three-Dimensional Unstructured Meshes}},
+  author={Kong, Fande and Wang, Yaqi and Gaston, Derek R and Permann, Cody J and Slaughter, Andrew E and Lindsay, Alexander D and DeHart, Mark D and Martineau, Richard C},
+  journal={SIAM Journal on Scientific Computing},
+  volume={42},
+  number={5},
+  pages={C193--C220},
+  year={2020},
+  publisher={SIAM}
+}
+```

framework/src/utils/PetscSupport.C

@@ -620,8 +620,9 @@ storePetscOptions(FEProblemBase & fe_problem, const InputParameters & params)
      * which happens before the parser is even created.  We'll throw an error if somebody attempts
      * to add this option later.
      */
-    if (option == "-log_summary")
-      mooseError("The PETSc option \"-log_summary\" can only be used on the command line.  Please "
+    if (option == "-log_summary" || option == "-log_view")
+      mooseError("The PETSc option \"-log_summary\" or \"-log_view\" can only be used on the "
+                 "command line.  Please "
                  "remove it from the input file");
 
     // Warn about superseded PETSc options (Note: -snes is not a REAL option, but people used it in
@@ -665,29 +666,80 @@ storePetscOptions(FEProblemBase & fe_problem, const InputParameters & params)
   bool tiny_pivot_found = false;
 #endif
   std::string pc_description = "";
+#if !PETSC_VERSION_LESS_THAN(3, 12, 0)
+  // If users use HMG, we would like to set
+  bool hmg_found = false;
+  bool matptap_found = false;
+  bool hmg_strong_threshold_found = false;
+#endif
+  bool option_changed = false;
+  std::vector<std::string> new_option_names;
+  std::vector<std::string> new_option_values;
   for (unsigned int i = 0; i < petsc_options_inames.size(); i++)
   {
+    new_option_names.clear();
+    new_option_values.clear();
+    option_changed = false;
     // Do not add duplicate settings
     if (find(po.inames.begin(), po.inames.end(), petsc_options_inames[i]) == po.inames.end())
     {
 #if !PETSC_VERSION_LESS_THAN(3, 9, 0)
       if (petsc_options_inames[i] == "-pc_factor_mat_solver_package")
-        po.inames.push_back("-pc_factor_mat_solver_type");
-      else
-        po.inames.push_back(petsc_options_inames[i]);
+      {
+        new_option_names.push_back("-pc_factor_mat_solver_type");
+        new_option_values.push_back(petsc_options_values[i]);
+        option_changed = true;
+      }
 #else
       if (petsc_options_inames[i] == "-pc_factor_mat_solver_type")
-        po.inames.push_back("-pc_factor_mat_solver_package");
-      else
-        po.inames.push_back(petsc_options_inames[i]);
+      {
+        new_option_names.push_back("-pc_factor_mat_solver_package");
+        new_option_values.push_back(petsc_options_values[i]);
+        option_changed = true;
+      }
 #endif
-      po.values.push_back(petsc_options_values[i]);
 
       // Look for a pc description
       if (petsc_options_inames[i] == "-pc_type" || petsc_options_inames[i] == "-pc_sub_type" ||
           petsc_options_inames[i] == "-pc_hypre_type")
         pc_description += petsc_options_values[i] + ' ';
 
+#if !PETSC_VERSION_LESS_THAN(3, 12, 0)
+      if (petsc_options_inames[i] == "-pc_type" && petsc_options_values[i] == "hmg")
+        hmg_found = true;
+
+        // MPIAIJ for PETSc 3.12.0: -matptap_via
+        // MAIJ for PETSc 3.12.0: -matmaijptap_via
+        // MPIAIJ for PETSc 3.13 or higher: -matptap_via, -matproduct_ptap_via
+        // MAIJ for PETSc 3.13 or higher: -matproduct_ptap_via
+#if !PETSC_VERSION_LESS_THAN(3, 13, 0)
+      if (hmg_found && (petsc_options_inames[i] == "-matptap_via" ||
+                        petsc_options_inames[i] == "-matmaijptap_via"))
+      {
+        new_option_names.push_back("-matproduct_ptap_via");
+        new_option_values.push_back(petsc_options_values[i]);
+        option_changed = true;
+      }
+#else
+      if (hmg_found && (petsc_options_inames[i] == "-matproduct_ptap_via"))
+      {
+        new_option_names.push_back("-matptap_via");
+        new_option_values.push_back(petsc_options_values[i]);
+        new_option_names.push_back("-matmaijptap_via");
+        new_option_values.push_back(petsc_options_values[i]);
+        option_changed = true;
+      }
+#endif
+
+      if (petsc_options_inames[i] == "-matptap_via" ||
+          petsc_options_inames[i] == "-matmaijptap_via" ||
+          petsc_options_inames[i] == "-matproduct_ptap_via")
+        matptap_found = true;
+
+      // For 3D problems, we need to set this 0.7
+      if (petsc_options_inames[i] == "-hmg_inner_pc_hypre_boomeramg_strong_threshold")
+        hmg_strong_threshold_found = true;
+#endif
       // This special case is common enough that we'd like to handle it for the user.
       if (petsc_options_inames[i] == "-pc_hypre_type" && petsc_options_values[i] == "boomeramg")
         boomeramg_found = true;
@@ -703,6 +755,20 @@ storePetscOptions(FEProblemBase & fe_problem, const InputParameters & params)
       if (petsc_options_inames[i] == "-mat_superlu_dist_replacetinypivot")
         tiny_pivot_found = true;
 #endif
+
+      if (option_changed)
+      {
+        for (MooseIndex(new_option_names.size()) op = 0; op < new_option_names.size(); op++)
+        {
+          po.inames.push_back(new_option_names[op]);
+          po.values.push_back(new_option_values[op]);
+        }
+      }
+      else
+      {
+        po.inames.push_back(petsc_options_inames[i]);
+        po.values.push_back(petsc_options_values[i]);
+      }
     }
     else
     {
@@ -722,6 +788,30 @@ storePetscOptions(FEProblemBase & fe_problem, const InputParameters & params)
     pc_description += "strong_threshold: 0.7 (auto)";
   }
 
+#if !PETSC_VERSION_LESS_THAN(3, 12, 0)
+  if (hmg_found && !hmg_strong_threshold_found && fe_problem.mesh().dimension() == 3)
+  {
+    po.inames.push_back("-hmg_inner_pc_hypre_boomeramg_strong_threshold");
+    po.values.push_back("0.7");
+    pc_description += "strong_threshold: 0.7 (auto)";
+  }
+
+  // Default PETSc PtAP takes too much memory, and it is not quite useful
+  // Let us switch to use new algorithm
+  if (hmg_found && !matptap_found)
+  {
+#if !PETSC_VERSION_LESS_THAN(3, 13, 0)
+    po.inames.push_back("-matproduct_ptap_via");
+    po.values.push_back("allatonce");
+#else
+    po.inames.push_back("-matptap_via");
+    po.values.push_back("allatonce");
+    po.inames.push_back("-matmaijptap_via");
+    po.values.push_back("allatonce");
+#endif
+  }
+#endif
+
 #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
   // In PETSc-3.7.{0--4}, there is a bug when using superlu_dist, and we have to use
   // SamePattern_SameRowPerm, otherwise we use whatever we have in PETSc

test/tests/preconditioners/hmg/diffusion_hmg.i

@@ -0,0 +1,87 @@
+[Mesh]
+  [./dmg]
+    type = DistributedRectilinearMeshGenerator
+    nx = 10
+    ny = 10
+    dim = 2
+  [../]
+[]
+
+[Variables]
+  [u1][]
+  [u2][]
+  [u3][]
+[]
+
+[Kernels]
+  [./diff_1]
+    type = Diffusion
+    variable = u1
+  [../]
+  [./diff_2]
+    type = Diffusion
+    variable = u2
+  [../]
+  [./diff_3]
+    type = Diffusion
+    variable = u3
+  [../]
+[]
+
+[BCs]
+  [./left_1]
+    type = DirichletBC
+    variable = u1
+    boundary = 'left'
+    value = 0
+  [../]
+
+  [./right_1]
+    type = DirichletBC
+    variable = u1
+    boundary = 'right'
+    value = 1
+  [../]
+
+  [./left_2]
+    type = DirichletBC
+    variable = u2
+    boundary = 'left'
+    value = 0
+  [../]
+
+  [./right_2]
+    type = DirichletBC
+    variable = u2
+    boundary = 'right'
+    value = 2
+  [../]
+
+  [./left_3]
+    type = DirichletBC
+    variable = u3
+    boundary = 'left'
+    value = 0
+  [../]
+
+  [./right_3]
+    type = DirichletBC
+    variable = u3
+    boundary = 'right'
+    value = 3
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+
+  solve_type = 'PJFNK'
+
+  petsc_options_iname = '-pc_type -pc_hmg_use_subspace_coarsening'
+  petsc_options_value = 'hmg true'
+  petsc_options = '-snes_view'
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/preconditioners/hmg/tests

@@ -0,0 +1,29 @@
+[Tests]
+  design = 'HMGPreconditioner.md'
+  issues = '#16210'
+
+  [hmg]
+    type = 'Exodiff'
+    input = 'diffusion_hmg.i'
+    exodiff = 'diffusion_hmg_out.e'
+    min_parallel = 2
+    # HMG was introduced in 3.12.0
+    petsc_version = '>=3.12.0'
+    requirement = "The system shall support the use of HMG (high performance MG)"
+    # Check if PtAP algorithms are setup correctly
+    expect_out = "using\s+allatonce\s+MatPtAP\(\)\s+implementation"
+  []
+
+  [hmg_3D]
+    type = 'Exodiff'
+    input = 'diffusion_hmg.i'
+    exodiff = 'diffusion_hmg_3d_out.e'
+    min_parallel = 2
+    # HMG was introduced in 3.12.0
+    petsc_version = '>=3.12.0'
+    cli_args = 'Mesh/dmg/dim=3 Mesh/dmg/nz=10 Outputs/file_base=diffusion_hmg_3d_out -log_view'
+    requirement = "The system shall support the use of HMG (high performance MG) for 3D problems"
+    # Check if strong_threshold is setup correctly
+    expect_out = "PETSc\s+Preconditioner:\s+hmg\s+strong_threshold:\s+0.7"
+  []
+[]