Preliminary refactoring and minor changes

- Some of these are necessary in the AMR branch - Some of these are minor cleanup/rationalization - Some are aesthetics
awslabs · Aug 23, 2023 · ee3a45b · ee3a45b
1 parent 240f6f1
commit ee3a45b
Show file tree

Hide file tree

Showing 10 changed files with 132 additions and 117 deletions.
diff --git a/palace/fem/integrator.hpp b/palace/fem/integrator.hpp
@@ -9,44 +9,42 @@
 namespace palace
 {
 
-//
-// Derived integrator classes extending the linear and bilinear form integrators of MFEM.
-//
-
-class DefaultIntegrationRule
+namespace fem
 {
-protected:
-  static const mfem::IntegrationRule *GetDefaultRule(const mfem::FiniteElement &trial_fe,
-                                                     const mfem::FiniteElement &test_fe,
-                                                     mfem::ElementTransformation &Tr)
-  {
-    const int ir_order = trial_fe.GetOrder() + test_fe.GetOrder() + Tr.OrderW();
-    return &mfem::IntRules.Get(trial_fe.GetGeomType(), ir_order);
-  }
+// Helper functions for creating an integration rule to exactly integrate 2p + q
+// polynomials. order_increment can be used to raise or lower the order, e.g. in
+// the case of derivative fes.
+inline const mfem::IntegrationRule *GetDefaultRule(const mfem::FiniteElement &trial_fe,
+                                                   const mfem::FiniteElement &test_fe,
+                                                   mfem::ElementTransformation &Tr,
+                                                   int order_increment = 0)
+{
+  const int ir_order =
+      trial_fe.GetOrder() + test_fe.GetOrder() + Tr.OrderW() + order_increment;
+  MFEM_ASSERT(ir_order >= 0, "Negative integration order not allowed");
+  return &mfem::IntRules.Get(trial_fe.GetGeomType(), ir_order);
+}
+inline const mfem::IntegrationRule *GetDefaultRule(const mfem::FiniteElement &fe,
+                                                   mfem::ElementTransformation &Tr,
+                                                   int order_increment = 0)
+{
+  return GetDefaultRule(fe, fe, Tr, order_increment);
+}
 
-  static const mfem::IntegrationRule *GetDefaultRule(const mfem::FiniteElement &fe,
-                                                     mfem::ElementTransformation &Tr)
-  {
-    return GetDefaultRule(fe, fe, Tr);
-  }
-};
+}  // namespace fem
 
 // Similar to MFEM's VectorFEBoundaryTangentLFIntegrator for ND spaces, but instead of
 // computing (n x f, v), this just computes (f, v). Also eliminates the a and b quadrature
-// parameters and uses GetDefaultRule instead.
-class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator,
-                                     public DefaultIntegrationRule
+// parameters and uses fem::GetDefaultRule instead.
+class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator
 {
 private:
   mfem::VectorCoefficient &Q;
   mfem::DenseMatrix vshape;
   mfem::Vector f_loc, f_hat;
 
 public:
-  VectorFEBoundaryLFIntegrator(mfem::VectorCoefficient &QG)
-    : Q(QG), f_loc(QG.GetVDim()), f_hat(QG.GetVDim() - 1)
-  {
-  }
+  VectorFEBoundaryLFIntegrator(mfem::VectorCoefficient &QG) : Q(QG), f_loc(QG.GetVDim()) {}
 
   void AssembleRHSElementVect(const mfem::FiniteElement &fe,
                               mfem::ElementTransformation &Tr,
@@ -55,10 +53,11 @@ class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator,
     const int dof = fe.GetDof();
     const int dim = fe.GetDim();
     const mfem::IntegrationRule *ir =
-        (IntRule != nullptr) ? IntRule : GetDefaultRule(fe, Tr);
+        (IntRule != nullptr) ? IntRule : fem::GetDefaultRule(fe, Tr);
     vshape.SetSize(dof, dim);
     elvect.SetSize(dof);
     elvect = 0.0;
+    f_hat.SetSize(dim);
 
     for (int i = 0; i < ir->GetNPoints(); i++)
     {
@@ -67,6 +66,7 @@ class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator,
       fe.CalcVShape(ip, vshape);
 
       Q.Eval(f_loc, Tr, ip);
+
       Tr.InverseJacobian().Mult(f_loc, f_hat);
       f_hat *= ip.weight * Tr.Weight();
       vshape.AddMult(f_hat, elvect);
@@ -75,9 +75,8 @@ class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator,
 };
 
 // Similar to MFEM's BoundaryLFIntegrator for H1 spaces, but eliminates the a and b
-// quadrature parameters and uses GetDefaultRule instead.
-class BoundaryLFIntegrator : public mfem::LinearFormIntegrator,
-                             public DefaultIntegrationRule
+// quadrature parameters and uses fem::GetDefaultRule instead.
+class BoundaryLFIntegrator : public mfem::LinearFormIntegrator
 {
 private:
   mfem::Coefficient &Q;
@@ -92,7 +91,7 @@ class BoundaryLFIntegrator : public mfem::LinearFormIntegrator,
   {
     const int dof = fe.GetDof();
     const mfem::IntegrationRule *ir =
-        (IntRule != nullptr) ? IntRule : GetDefaultRule(fe, Tr);
+        (IntRule != nullptr) ? IntRule : fem::GetDefaultRule(fe, Tr);
     shape.SetSize(dof);
     elvect.SetSize(dof);
     elvect = 0.0;

diff --git a/palace/fem/multigrid.hpp b/palace/fem/multigrid.hpp
@@ -15,7 +15,7 @@ namespace palace::utils
 {
 
 //
-// Methods for constructing hierarchies of finite element spaces for geometric multigrid.
+// Methods for constructing hierarchies of finite element spaces for multigrid.
 //
 
 // Helper function for getting the order of the finite element space underlying a bilinear
@@ -122,37 +122,27 @@ std::vector<std::unique_ptr<FECollection>> inline ConstructFECollections(
 }
 
 // Construct a hierarchy of finite element spaces given a sequence of meshes and
-// finite element collections. Dirichlet boundary conditions are additionally
-// marked.
+// finite element collections. Uses geometric multigrid and p multigrid.
 template <typename FECollection>
 inline mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
     int mg_max_levels, bool mg_legacy_transfer, int pa_order_threshold,
     const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
-    const std::vector<std::unique_ptr<FECollection>> &fecs,
-    const mfem::Array<int> *dbc_marker = nullptr,
-    std::vector<mfem::Array<int>> *dbc_tdof_lists = nullptr)
+    const std::vector<std::unique_ptr<FECollection>> &fecs, int dim = 1)
 {
-  MFEM_VERIFY(!mesh.empty() && !fecs.empty() &&
-                  (!dbc_tdof_lists || dbc_tdof_lists->empty()),
+  MFEM_VERIFY(!mesh.empty() && !fecs.empty(),
               "Empty mesh or FE collection for FE space construction!");
   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());
-  }
+  auto *fespace =
+      new mfem::ParFiniteElementSpace(mesh[coarse_mesh_l].get(), fecs[0].get(), dim);
   mfem::ParFiniteElementSpaceHierarchy fespaces(mesh[coarse_mesh_l].get(), fespace, false,
                                                 true);
 
   // h-refinement
   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)
-    {
-      fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back());
-    }
+    fespace = new mfem::ParFiniteElementSpace(mesh[l].get(), fecs[0].get(), dim);
+
     auto *P = new ParOperator(
         std::make_unique<mfem::TransferOperator>(fespaces.GetFinestFESpace(), *fespace),
         fespaces.GetFinestFESpace(), *fespace, true);
@@ -162,11 +152,8 @@ inline mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy
   // p-refinement
   for (std::size_t l = 1; l < fecs.size(); l++)
   {
-    fespace = new mfem::ParFiniteElementSpace(mesh.back().get(), fecs[l].get());
-    if (dbc_marker && dbc_tdof_lists)
-    {
-      fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back());
-    }
+    fespace = new mfem::ParFiniteElementSpace(mesh.back().get(), fecs[l].get(), dim);
+
     ParOperator *P;
     if (!mg_legacy_transfer && mfem::DeviceCanUseCeed())
     {
@@ -188,6 +175,27 @@ inline mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy
   return fespaces;
 }
 
+// Overload for treatment of Dirichlet boundary conditions, extracts the true
+// dof vectors for each level within a finite element space hierarchy.
+template <typename FECollection, template <class... C> typename Container,
+          typename... MeshT>
+mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
+    int mg_max_levels, bool mg_legacy_transfer, int pa_order_threshold,
+    const Container<MeshT...> &mesh, const std::vector<std::unique_ptr<FECollection>> &fecs,
+    const mfem::Array<int> &dbc_marker, std::vector<mfem::Array<int>> &dbc_tdof_lists,
+    int dim = 1)
+{
+  dbc_tdof_lists.clear();
+  auto fespaces = ConstructFiniteElementSpaceHierarchy(mg_max_levels, mg_legacy_transfer,
+                                                       pa_order_threshold, mesh, fecs, dim);
+  for (int l = 0; l < fespaces.GetNumLevels(); l++)
+  {
+    fespaces.GetFESpaceAtLevel(l).GetEssentialTrueDofs(dbc_marker,
+                                                       dbc_tdof_lists.emplace_back());
+  }
+  return fespaces;
+}
+
 }  // namespace palace::utils
 
 #endif  // PALACE_FEM_MULTIGRID_HPP
diff --git a/palace/main.cpp b/palace/main.cpp
@@ -7,6 +7,7 @@
 #include <string_view>
 #include <vector>
 #include <mfem.hpp>
+
 #include "drivers/drivensolver.hpp"
 #include "drivers/eigensolver.hpp"
 #include "drivers/electrostaticsolver.hpp"
@@ -153,30 +154,28 @@ int main(int argc, char *argv[])
 #endif
 
   // Initialize the problem driver.
-  std::unique_ptr<BaseSolver> solver;
-  switch (iodata.problem.type)
+  const std::unique_ptr<BaseSolver> solver = [&]() -> std::unique_ptr<BaseSolver>
   {
-    case config::ProblemData::Type::DRIVEN:
-      solver = std::make_unique<DrivenSolver>(iodata, world_root, world_size, num_thread,
+    switch (iodata.problem.type)
+    {
+      case config::ProblemData::Type::DRIVEN:
+        return std::make_unique<DrivenSolver>(iodata, world_root, world_size, num_thread,
                                               git_tag);
-      break;
-    case config::ProblemData::Type::EIGENMODE:
-      solver = std::make_unique<EigenSolver>(iodata, world_root, world_size, num_thread,
+      case config::ProblemData::Type::EIGENMODE:
+        return std::make_unique<EigenSolver>(iodata, world_root, world_size, num_thread,
                                              git_tag);
-      break;
-    case config::ProblemData::Type::ELECTROSTATIC:
-      solver = std::make_unique<ElectrostaticSolver>(iodata, world_root, world_size,
+      case config::ProblemData::Type::ELECTROSTATIC:
+        return std::make_unique<ElectrostaticSolver>(iodata, world_root, world_size,
                                                      num_thread, git_tag);
-      break;
-    case config::ProblemData::Type::MAGNETOSTATIC:
-      solver = std::make_unique<MagnetostaticSolver>(iodata, world_root, world_size,
+      case config::ProblemData::Type::MAGNETOSTATIC:
+        return std::make_unique<MagnetostaticSolver>(iodata, world_root, world_size,
                                                      num_thread, git_tag);
-      break;
-    case config::ProblemData::Type::TRANSIENT:
-      solver = std::make_unique<TransientSolver>(iodata, world_root, world_size, num_thread,
+      case config::ProblemData::Type::TRANSIENT:
+        return std::make_unique<TransientSolver>(iodata, world_root, world_size, num_thread,
                                                  git_tag);
-      break;
-  }
+    }
+    return nullptr;
+  }();
 
   // Read the mesh from file, refine, partition, and distribute it. Then nondimensionalize
   // it and the input parameters.
@@ -188,6 +187,7 @@ int main(int argc, char *argv[])
 
   // Run the problem driver.
   solver->Solve(mesh, timer);
+
   timer.Reduce(world_comm);
   timer.Print(world_comm);
   solver->SaveMetadata(timer);

diff --git a/palace/models/curlcurloperator.cpp b/palace/models/curlcurloperator.cpp
@@ -81,10 +81,10 @@ CurlCurlOperator::CurlCurlOperator(const IoData &iodata,
     rt_fec(iodata.solver.order - 1, mesh.back()->Dimension()),
     nd_fespaces(utils::ConstructFiniteElementSpaceHierarchy<mfem::ND_FECollection>(
         iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer,
-        pa_order_threshold, mesh, nd_fecs, &dbc_marker, &dbc_tdof_lists)),
+        pa_order_threshold, mesh, nd_fecs, dbc_marker, dbc_tdof_lists)),
     h1_fespaces(utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
         iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer,
-        pa_order_threshold, mesh, h1_fecs, nullptr, nullptr)),
+        pa_order_threshold, mesh, h1_fecs)),
     rt_fespace(mesh.back().get(), &rt_fec), mat_op(iodata, *mesh.back()),
     surf_j_op(iodata, GetH1Space())
 {

diff --git a/palace/models/curlcurloperator.hpp b/palace/models/curlcurloperator.hpp
@@ -69,6 +69,9 @@ class CurlCurlOperator
   auto &GetRTSpace() { return rt_fespace; }
   const auto &GetRTSpace() const { return rt_fespace; }
 
+  // Return the number of true (conforming) dofs on the finest ND space.
+  auto GetNDof() { return GetNDSpace().GlobalTrueVSize(); }
+
   // Construct and return system matrix representing discretized curl-curl operator for
   // Ampere's law.
   std::unique_ptr<Operator> GetStiffnessMatrix();

diff --git a/palace/models/laplaceoperator.cpp b/palace/models/laplaceoperator.cpp
@@ -121,7 +121,7 @@ LaplaceOperator::LaplaceOperator(const IoData &iodata,
     nd_fec(iodata.solver.order, mesh.back()->Dimension()),
     h1_fespaces(utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
         iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer,
-        pa_order_threshold, mesh, h1_fecs, &dbc_marker, &dbc_tdof_lists)),
+        pa_order_threshold, mesh, h1_fecs, dbc_marker, dbc_tdof_lists)),
     nd_fespace(mesh.back().get(), &nd_fec), mat_op(iodata, *mesh.back()),
     source_attr_lists(ConstructSources(iodata))
 {

diff --git a/palace/models/laplaceoperator.hpp b/palace/models/laplaceoperator.hpp
@@ -63,6 +63,9 @@ class LaplaceOperator
   auto &GetNDSpace() { return nd_fespace; }
   const auto &GetNDSpace() const { return nd_fespace; }
 
+  // Return the number of true (conforming) dofs on the finest H1 space.
+  auto GetNDof() { return GetH1Space().GlobalTrueVSize(); }
+
   // Construct and return system matrix representing discretized Laplace operator for
   // Gauss's law.
   std::unique_ptr<Operator> GetStiffnessMatrix();

diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp
@@ -86,10 +86,10 @@ SpaceOperator::SpaceOperator(const IoData &iodata,
     rt_fec(iodata.solver.order - 1, mesh.back()->Dimension()),
     nd_fespaces(utils::ConstructFiniteElementSpaceHierarchy<mfem::ND_FECollection>(
         iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer,
-        pa_order_threshold, mesh, nd_fecs, &dbc_marker, &nd_dbc_tdof_lists)),
+        pa_order_threshold, mesh, nd_fecs, dbc_marker, nd_dbc_tdof_lists)),
     h1_fespaces(utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
         iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer,
-        pa_order_threshold, mesh, h1_fecs, &dbc_marker, &h1_dbc_tdof_lists)),
+        pa_order_threshold, mesh, h1_fecs, dbc_marker, h1_dbc_tdof_lists)),
     rt_fespace(mesh.back().get(), &rt_fec), mat_op(iodata, *mesh.back()),
     farfield_op(iodata, mat_op, *mesh.back()), surf_sigma_op(iodata, *mesh.back()),
     surf_z_op(iodata, *mesh.back()), lumped_port_op(iodata, GetH1Space()),

diff --git a/palace/models/spaceoperator.hpp b/palace/models/spaceoperator.hpp
@@ -126,6 +126,9 @@ class SpaceOperator
   auto &GetRTSpace() { return rt_fespace; }
   const auto &GetRTSpace() const { return rt_fespace; }
 
+  // Return the number of true (conforming) dofs on the finest ND space.
+  auto GetNDof() { return GetNDSpace().GlobalTrueVSize(); }
+
   // Construct any part of the frequency-dependent complex linear system matrix:
   //                     A = K + iω C - ω² (Mr + i Mi) + A2(ω) .
   // For time domain problems, any one of K, C, or M = Mr can be constructed. The argument