Address PR comments: Clarify comments, fix some typos, rename HCurlNo…

…rmSolver -> EnergyNormSolver
awslabs · Aug 15, 2023 · 6fe4ca5 · 6fe4ca5
1 parent 7f2721b
commit 6fe4ca5
Show file tree

Hide file tree

Showing 7 changed files with 33 additions and 33 deletions.
diff --git a/palace/fem/multigrid.hpp b/palace/fem/multigrid.hpp
@@ -33,7 +33,7 @@ inline auto GetMaxElementOrder(mfem::MixedBilinearForm &a)
                   a.TrialFESpace()->GetMaxElementOrder());
 }
 
-// Assembly a bilinear or mixed bilinear form. If the order is lower than the specified
+// Assemble a bilinear or mixed bilinear form. If the order is lower than the specified
 // threshold, the operator is assembled as a sparse matrix.
 template <typename BilinearForm>
 inline std::unique_ptr<Operator>
@@ -77,16 +77,16 @@ std::vector<std::unique_ptr<FECollection>> inline ConstructFECollections(
 {
   // If the solver will use a LOR preconditioner, we need to construct with a specific basis
   // type.
-  MFEM_VERIFY(p >= 1, "FE space order must not be less than 1!");
+  constexpr int pmin = (std::is_base_of<mfem::H1_FECollection, FECollection>::value ||
+                        std::is_base_of<mfem::ND_FECollection, FECollection>::value)
+                           ? 1
+                           : 0;
+  MFEM_VERIFY(p >= pmin, "FE space order must not be less than " << pmin << "!");
   int b1 = mfem::BasisType::GaussLobatto, b2 = mfem::BasisType::GaussLegendre;
   if (mat_lor)
   {
     b2 = mfem::BasisType::IntegratedGLL;
   }
-  constexpr int pm1 = (std::is_base_of<mfem::H1_FECollection, FECollection>::value ||
-                       std::is_base_of<mfem::ND_FECollection, FECollection>::value)
-                          ? 0
-                          : 1;
 
   // Construct the p-multigrid hierarchy, first finest to coarsest and then reverse the
   // order.
@@ -96,14 +96,14 @@ std::vector<std::unique_ptr<FECollection>> inline ConstructFECollections(
     if constexpr (std::is_base_of<mfem::ND_FECollection, FECollection>::value ||
                   std::is_base_of<mfem::RT_FECollection, FECollection>::value)
     {
-      fecs.push_back(std::make_unique<FECollection>(p - pm1, dim, b1, b2));
+      fecs.push_back(std::make_unique<FECollection>(p, dim, b1, b2));
     }
     else
     {
-      fecs.push_back(std::make_unique<FECollection>(p - pm1, dim, b1));
+      fecs.push_back(std::make_unique<FECollection>(p, dim, b1));
       MFEM_CONTRACT_VAR(b2);
     }
-    if (p == 1)
+    if (p == pmin)
     {
       break;
     }
@@ -113,7 +113,7 @@ std::vector<std::unique_ptr<FECollection>> inline ConstructFECollections(
         p--;
         break;
       case config::LinearSolverData::MultigridCoarsenType::LOGARITHMIC:
-        p = (p + 1) / 2;
+        p = (p + pmin) / 2;
         break;
       case config::LinearSolverData::MultigridCoarsenType::INVALID:
         MFEM_ABORT("Invalid coarsening type for p-multigrid levels!");
@@ -138,18 +138,18 @@ inline mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy
   MFEM_VERIFY(!mesh.empty() && !fecs.empty() &&
                   (!dbc_tdof_lists || dbc_tdof_lists->empty()),
               "Empty mesh or FE collection for FE space construction!");
-  auto mesh_levels = std::min(mesh.size() - 1, mg_max_levels - fecs.size());
-  auto *fespace = new mfem::ParFiniteElementSpace(mesh[mesh.size() - mesh_levels - 1].get(),
-                                                  fecs[0].get());
+  int coarse_mesh_l =
+      std::max(0, static_cast<int>(mesh.size() + fecs.size()) - 1 - mg_max_levels);
+  auto *fespace = new mfem::ParFiniteElementSpace(mesh[coarse_mesh_l].get(), fecs[0].get());
   if (dbc_marker && dbc_tdof_lists)
   {
     fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back());
   }
-  mfem::ParFiniteElementSpaceHierarchy fespaces(mesh[mesh.size() - mesh_levels - 1].get(),
-                                                fespace, false, true);
+  mfem::ParFiniteElementSpaceHierarchy fespaces(mesh[coarse_mesh_l].get(), fespace, false,
+                                                true);
 
   // h-refinement
-  for (std::size_t l = mesh.size() - mesh_levels; l < mesh.size(); l++)
+  for (std::size_t l = coarse_mesh_l + 1; l < mesh.size(); l++)
   {
     fespace = new mfem::ParFiniteElementSpace(mesh[l].get(), fecs[0].get());
     if (dbc_marker && dbc_tdof_lists)

diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp
@@ -15,12 +15,13 @@
 namespace palace
 {
 
-HCurlNormSolver::HCurlNormSolver(const MaterialOperator &mat_op,
-                                 mfem::ParFiniteElementSpaceHierarchy &nd_fespaces,
-                                 mfem::ParFiniteElementSpaceHierarchy &h1_fespaces,
-                                 const std::vector<mfem::Array<int>> &nd_dbc_tdof_lists,
-                                 const std::vector<mfem::Array<int>> &h1_dbc_tdof_lists,
-                                 double tol, int max_it, int print, int pa_order_threshold)
+EnergyNormSolver::EnergyNormSolver(const MaterialOperator &mat_op,
+                                   mfem::ParFiniteElementSpaceHierarchy &nd_fespaces,
+                                   mfem::ParFiniteElementSpaceHierarchy &h1_fespaces,
+                                   const std::vector<mfem::Array<int>> &nd_dbc_tdof_lists,
+                                   const std::vector<mfem::Array<int>> &h1_dbc_tdof_lists,
+                                   double tol, int max_it, int print,
+                                   int pa_order_threshold)
 {
   constexpr int skip_zeros = 0;
   constexpr auto MatTypeMuInv = MaterialPropertyType::INV_PERMEABILITY;

diff --git a/palace/linalg/hcurl.hpp b/palace/linalg/hcurl.hpp
@@ -27,7 +27,7 @@ class MaterialOperator;
 //
 // This solver implements a solver for the operator K + M in a Nedelec space.
 //
-class HCurlNormSolver
+class EnergyNormSolver
 {
 private:
   // H(curl) norm operator A = K + M and its projection Gᵀ A G.
@@ -37,12 +37,12 @@ class HCurlNormSolver
   std::unique_ptr<KspSolver> ksp;
 
 public:
-  HCurlNormSolver(const MaterialOperator &mat_op,
-                  mfem::ParFiniteElementSpaceHierarchy &nd_fespaces,
-                  mfem::ParFiniteElementSpaceHierarchy &h1_fespaces,
-                  const std::vector<mfem::Array<int>> &nd_dbc_tdof_lists,
-                  const std::vector<mfem::Array<int>> &h1_dbc_tdof_lists, double tol,
-                  int max_it, int print, int pa_order_threshold);
+  EnergyNormSolver(const MaterialOperator &mat_op,
+                   mfem::ParFiniteElementSpaceHierarchy &nd_fespaces,
+                   mfem::ParFiniteElementSpaceHierarchy &h1_fespaces,
+                   const std::vector<mfem::Array<int>> &nd_dbc_tdof_lists,
+                   const std::vector<mfem::Array<int>> &h1_dbc_tdof_lists, double tol,
+                   int max_it, int print, int pa_order_threshold);
 
   const Operator &GetOperator() { return *A; }
 

diff --git a/palace/models/romoperator.cpp b/palace/models/romoperator.cpp
@@ -114,7 +114,7 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &spaceop) : spaceop
   if (iodata.solver.driven.adaptive_metric_aposteriori)
   {
     constexpr int curlcurl_verbose = 0;
-    kspKM = std::make_unique<HCurlNormSolver>(
+    kspKM = std::make_unique<EnergyNormSolver>(
         spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), spaceop.GetH1Spaces(),
         spaceop.GetNDDbcTDofLists(), spaceop.GetH1DbcTDofLists(), iodata.solver.linear.tol,
         iodata.solver.linear.max_it, curlcurl_verbose, iodata.solver.pa_order_threshold);

diff --git a/palace/models/romoperator.hpp b/palace/models/romoperator.hpp
@@ -42,7 +42,7 @@ class RomOperator
   std::unique_ptr<ComplexKspSolver> ksp;
 
   // Linear solver for inner product solves for error metric.
-  std::unique_ptr<HCurlNormSolver> kspKM;
+  std::unique_ptr<EnergyNormSolver> kspKM;
 
   // PROM matrices and vectors.
   Eigen::MatrixXcd Kr, Mr, Cr, Ar;

diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp
@@ -740,7 +740,7 @@ struct LinearSolverData
   // Maximum number of levels for geometric multigrid (set to 1 to disable multigrid).
   int mg_max_levels = 100;
 
-  // Type of coarsening for geometric multigrid.
+  // Type of coarsening for p-multigrid.
   enum class MultigridCoarsenType
   {
     LINEAR,

diff --git a/palace/utils/iodata.cpp b/palace/utils/iodata.cpp
@@ -330,7 +330,6 @@ void IoData::CheckConfiguration()
   // Resolve default values in configuration file.
 
   // XX TODO: Default value for pa_order_threshold if we want PA enabled by default
-  // XX TODO: Also default value for mg_legacy_transfer, maybe always true for PA?
 
   if (solver.linear.max_size < 0)
   {