Expose AMS and AMG options in configuration file

docs/src/config/model.md

@@ -64,7 +64,7 @@ scale based on the bounding box of the computational domain.
 
 with
 
-`"Tol" [1e-2]` : Relative error convergence tolerance for adaptive mesh refinement (AMR).
+`"Tol" [1.0e-2]` : Relative error convergence tolerance for adaptive mesh refinement (AMR).
 
 `"MaxIts" [0]` : Maximum number of iterations of AMR to perform.
 

docs/src/config/solver.md

@@ -434,7 +434,7 @@ eigenmode simulation type.
 `"DivFreeMaxIts" [1000]` :  Maximum number of iterations for divergence-free cleaning use in
 the eigenmode simulation type.
 
-`"EstimatorTol" [1e-6]` :  Relative tolerance for flux projection used in the
+`"EstimatorTol" [1.0e-6]` :  Relative tolerance for flux projection used in the
 error estimate calculation.
 
 `"EstimatorMaxIts" [10000]` :  Maximum number of iterations for flux projection use in the
@@ -465,5 +465,7 @@ vectors in Krylov subspace methods or other parts of the code.
   - `"STRUMPACKCompressionTol" [1.0e-3]`
   - `"STRUMPACKLossyPrecision" [16]`
   - `"STRUMPACKButterflyLevels" [1]`
-  - `"SuperLU3D" [false]`
-  - `"AMSVector" [false]`
+  - `"SuperLU3DCommunicator" [false]`
+  - `"AMSVectorInterpolation" [false]`
+  - `"AMSSingularOperator" [false]`
+  - `"AMGAggressiveCoarsening" [false]`

palace/linalg/amg.cpp

@@ -6,15 +6,15 @@
 namespace palace
 {
 
-BoomerAmgSolver::BoomerAmgSolver(int cycle_it, int smooth_it, int print)
+BoomerAmgSolver::BoomerAmgSolver(int cycle_it, int smooth_it, bool agg_coarsen, int print)
   : mfem::HypreBoomerAMG()
 {
   HYPRE_BoomerAMGSetPrintLevel(*this, (print > 1) ? print - 1 : 0);
   HYPRE_BoomerAMGSetMaxIter(*this, cycle_it);
   HYPRE_BoomerAMGSetTol(*this, 0.0);
 
   // Set additional BoomerAMG options.
-  int agg_levels = 1;  // Number of aggressive coarsening levels
+  int agg_levels = agg_coarsen ? 1 : 0;  // Number of aggressive coarsening levels
   double theta = 0.5;  // AMG strength parameter = 0.25 is 2D optimal (0.5-0.8 for 3D)
   int relax_type = 8;  // 8 = l1-symm. GS, 13 = l1-GS, 18 = l1-Jacobi, 16 = Chebyshev
   if (mfem::Device::Allows(mfem::Backend::DEVICE_MASK))

palace/linalg/amg.hpp

@@ -16,10 +16,12 @@ namespace palace
 class BoomerAmgSolver : public mfem::HypreBoomerAMG
 {
 public:
-  BoomerAmgSolver(int cycle_it = 1, int smooth_it = 1, int print = 0);
+  BoomerAmgSolver(int cycle_it = 1, int smooth_it = 1, bool agg_coarsen = true,
+                  int print = 0);
   BoomerAmgSolver(const IoData &iodata, bool coarse_solver, int print)
     : BoomerAmgSolver(coarse_solver ? 1 : iodata.solver.linear.mg_cycle_it,
-                      iodata.solver.linear.mg_smooth_it, print)
+                      iodata.solver.linear.mg_smooth_it,
+                      iodata.solver.linear.amg_agg_coarsen, print)
   {
   }
 };

palace/linalg/ams.cpp

@@ -15,7 +15,8 @@ namespace palace
 
 HypreAmsSolver::HypreAmsSolver(FiniteElementSpace &nd_fespace,
                                FiniteElementSpace &h1_fespace, int cycle_it, int smooth_it,
-                               bool vector_interp, bool op_pos, bool op_singular, int print)
+                               bool vector_interp, bool singular_op, bool agg_coarsen,
+                               int print)
   : mfem::HypreSolver(),
     // From the Hypre docs for AMS: cycles 1, 5, 8, 11, 13 are fastest, 7 yields fewest its
     // (MFEM default is 13). 14 is similar to 11/13 but is cheaper in that is uses additive
@@ -24,12 +25,12 @@ HypreAmsSolver::HypreAmsSolver(FiniteElementSpace &nd_fespace,
     // When used as the coarse solver of geometric multigrid, always do only a single
     // V-cycle.
     ams_it(cycle_it), ams_smooth_it(smooth_it),
-    // For positive (SPD) operators, we will use aggressive coarsening but not for frequency
-    // domain problems when the preconditioner matrix is not SPD.
-    ams_pos(op_pos),
     // If we know the operator is singular (no mass matrix, for magnetostatic problems),
     // internally the AMS solver will avoid G-space corrections.
-    ams_singular(op_singular), print((print > 1) ? print - 1 : 0)
+    ams_singular(singular_op),
+    // For positive (SPD) operators, we will use aggressive coarsening but not for frequency
+    // domain problems when the preconditioner matrix is not SPD.
+    agg_coarsen(agg_coarsen), print((print > 1) ? print - 1 : 0)
 {
   // From MFEM: The AMS preconditioner may sometimes require inverting singular matrices
   // with BoomerAMG, which are handled correctly in Hypre's Solve method, but can produce
@@ -151,8 +152,8 @@ void HypreAmsSolver::InitializeSolver()
   }
 
   // Set additional AMS options.
-  int coarsen_type = 10;                 // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
-  int amg_agg_levels = ams_pos ? 1 : 0;  // Number of aggressive coarsening levels
+  int coarsen_type = 10;                     // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
+  int amg_agg_levels = agg_coarsen ? 1 : 0;  // Number of aggressive coarsening levels
   double theta = 0.5;      // AMG strength parameter = 0.25 is 2D optimal (0.5-0.8 for 3D)
   int amg_relax_type = 8;  // 3 = GS, 6 = symm. GS, 8 = l1-symm. GS, 13 = l1-GS,
                            // 18 = l1-Jacobi, 16 = Chebyshev

palace/linalg/ams.hpp

@@ -25,7 +25,7 @@ class HypreAmsSolver : public mfem::HypreSolver
 
   // Parameters used for preconditioner construction.
   const int cycle_type, space_dim, ams_it, ams_smooth_it;
-  const bool ams_pos, ams_singular;
+  const bool ams_singular, agg_coarsen;
 
   // Control print level for debugging.
   const int print;
@@ -49,16 +49,14 @@ class HypreAmsSolver : public mfem::HypreSolver
   // Constructor requires the ND space, but will construct the H1 and (H1)ᵈ spaces
   // internally as needed.
   HypreAmsSolver(FiniteElementSpace &nd_fespace, FiniteElementSpace &h1_fespace,
-                 int cycle_it, int smooth_it, bool vector_interp, bool op_pos,
-                 bool op_singular, int print);
+                 int cycle_it, int smooth_it, bool vector_interp, bool singular_op,
+                 bool agg_coarsen, int print);
   HypreAmsSolver(const IoData &iodata, bool coarse_solver, FiniteElementSpace &nd_fespace,
                  FiniteElementSpace &h1_fespace, int print)
     : HypreAmsSolver(
           nd_fespace, h1_fespace, coarse_solver ? 1 : iodata.solver.linear.mg_cycle_it,
-          iodata.solver.linear.mg_smooth_it, iodata.solver.linear.ams_vector,
-          (iodata.problem.type == config::ProblemData::Type::TRANSIENT ||
-           iodata.problem.type == config::ProblemData::Type::MAGNETOSTATIC),
-          (iodata.problem.type == config::ProblemData::Type::MAGNETOSTATIC), print)
+          iodata.solver.linear.mg_smooth_it, iodata.solver.linear.ams_vector_interp,
+          iodata.solver.linear.ams_singular_op, iodata.solver.linear.amg_agg_coarsen, print)
   {
   }
   ~HypreAmsSolver() override;

palace/linalg/divfree.cpp

@@ -83,7 +83,7 @@ DivFreeSolver<VecType>::DivFreeSolver(
 
   // The system matrix for the projection is real and SPD.
   auto amg = std::make_unique<MfemWrapperSolver<OperType>>(
-      std::make_unique<BoomerAmgSolver>(1, 1, 0));
+      std::make_unique<BoomerAmgSolver>(1, 1, true, 0));
   std::unique_ptr<Solver<OperType>> pc;
   if (h1_fespaces.GetNumLevels() > 1)
   {

palace/linalg/errorestimator.cpp

@@ -76,7 +76,7 @@ auto ConfigureLinearSolver(const FiniteElementSpaceHierarchy &fespaces, double t
   }
   else
   {
-    auto amg = std::make_unique<BoomerAmgSolver>(1, 1, 0);
+    auto amg = std::make_unique<BoomerAmgSolver>(1, 1, true, 0);
     amg->SetStrengthThresh(0.8);  // More coarsening to save memory
     if (fespaces.GetNumLevels() > 1)
     {

palace/utils/configfile.cpp

@@ -1705,8 +1705,10 @@ void LinearSolverData::SetUp(json &solver)
   strumpack_lossy_precision =
       linear->value("STRUMPACKLossyPrecision", strumpack_lossy_precision);
   strumpack_butterfly_l = linear->value("STRUMPACKButterflyLevels", strumpack_butterfly_l);
-  superlu_3d = linear->value("SuperLU3D", superlu_3d);
-  ams_vector = linear->value("AMSVector", ams_vector);
+  superlu_3d = linear->value("SuperLU3DCommunicator", superlu_3d);
+  ams_vector_interp = linear->value("AMSVectorInterpolation", ams_vector_interp);
+  ams_singular_op = linear->value("AMSSingularOperator", ams_singular_op);
+  amg_agg_coarsen = linear->value("AMGAggressiveCoarsening", amg_agg_coarsen);
 
   // Other linear solver options.
   divfree_tol = linear->value("DivFreeTol", divfree_tol);
@@ -1743,8 +1745,10 @@ void LinearSolverData::SetUp(json &solver)
   linear->erase("STRUMPACKCompressionTol");
   linear->erase("STRUMPACKLossyPrecision");
   linear->erase("STRUMPACKButterflyLevels");
-  linear->erase("SuperLU3D");
-  linear->erase("AMSVector");
+  linear->erase("SuperLU3DCommunicator");
+  linear->erase("AMSVectorInterpolation");
+  linear->erase("AMSSingularOperator");
+  linear->erase("AMGAggressiveCoarsening");
 
   linear->erase("DivFreeTol");
   linear->erase("DivFreeMaxIts");
@@ -1785,8 +1789,11 @@ void LinearSolverData::SetUp(json &solver)
     std::cout << "STRUMPACKCompressionTol: " << strumpack_lr_tol << '\n';
     std::cout << "STRUMPACKLossyPrecision: " << strumpack_lossy_precision << '\n';
     std::cout << "STRUMPACKButterflyLevels: " << strumpack_butterfly_l << '\n';
-    std::cout << "SuperLU3D: " << superlu_3d << '\n';
-    std::cout << "AMSVector: " << ams_vector << '\n';
+    std::cout << "SuperLU3DCommunicator: " << superlu_3d << '\n';
+    std::cout << "AMSVectorInterpolation: " << ams_vector_interp << '\n';
+    std::cout << "AMSSingularOperator: " << ams_singular_op << '\n';
+    std::cout << "AMGAggressiveCoarsening: " << amg_agg_coarsen << '\n';
+
     std::cout << "DivFreeTol: " << divfree_tol << '\n';
     std::cout << "DivFreeMaxIts: " << divfree_max_it << '\n';
     std::cout << "EstimatorTol: " << estimator_tol << '\n';

palace/utils/configfile.hpp

@@ -137,7 +137,7 @@ struct RefinementData
 {
 public:
   // Non-dimensional tolerance used to specify convergence of adaptive mesh refinement.
-  double tol = 1e-2;
+  double tol = 1.0e-2;
 
   // Maximum number of iterations to perform during adaptive mesh refinement.
   int max_it = 0;
@@ -864,7 +864,15 @@ struct LinearSolverData
   bool superlu_3d = false;
 
   // Option to use vector or scalar Pi-space corrections for the AMS preconditioner.
-  bool ams_vector = false;
+  bool ams_vector_interp = false;
+
+  // Option to tell the AMS solver that the operator is singular, like for magnetostatic
+  // problems.
+  int ams_singular_op = -1;
+
+  // Option to use aggressive coarsening for Hypre AMG solves (with BoomerAMG or AMS).
+  // Typically use this when the operator is positive definite.
+  int amg_agg_coarsen = -1;
 
   // Relative tolerance for solving linear systems in divergence-free projector.
   double divfree_tol = 1.0e-12;

palace/utils/iodata.cpp

@@ -415,6 +415,18 @@ void IoData::CheckConfiguration()
   {
     solver.linear.mg_smooth_order = std::max(2 * solver.order, 4);
   }
+  if (solver.linear.ams_singular_op < 0)
+  {
+    solver.linear.ams_singular_op =
+        (problem.type == config::ProblemData::Type::MAGNETOSTATIC);
+  }
+  if (solver.linear.amg_agg_coarsen < 0)
+  {
+    solver.linear.amg_agg_coarsen =
+        (problem.type == config::ProblemData::Type::ELECTROSTATIC ||
+         problem.type == config::ProblemData::Type::MAGNETOSTATIC ||
+         problem.type == config::ProblemData::Type::TRANSIENT);
+  }
 
   // Configure settings for quadrature rules and partial assembly.
   BilinearForm::pa_order_threshold = solver.pa_order_threshold;

scripts/schema/config/solver.json

@@ -121,8 +121,10 @@
         "STRUMPACKCompressionTol": { "type": "number", "minimum": 0.0 },
         "STRUMPACKLossyPrecision": { "type": "integer", "minimum": 0 },
         "STRUMPACKButterflyLevels": { "type": "integer", "minimum": 0 },
-        "SuperLU3D": { "type": "boolean" },
-        "AMSVector": { "type": "boolean" },
+        "SuperLU3DCommunicator": { "type": "boolean" },
+        "AMSVectorInterpolation": { "type": "boolean" },
+        "AMSSingularOperator": { "type": "boolean" },
+        "AMGAggressiveCoarsening": { "type": "boolean" },
         "DivFreeTol": { "type": "number", "minimum": 0.0 },
         "DivFreeMaxIts": { "type": "integer", "minimum": 0 },
         "EstimatorTol": { "type": "number", "minimum": 0.0 },