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 1, 2023 · 52147a2 · 52147a2
1 parent c6fcae6
commit 52147a2
Show file tree

Hide file tree

Showing 11 changed files with 169 additions and 116 deletions.
diff --git a/.gitignore b/.gitignore
@@ -2,6 +2,7 @@ build*/
 examples/**/postpro*/
 examples/**/log*/
 spack/local/packages/palace/__pycache__/
+examples/**/*.log
 test/ref/**/*.json
 .gitlab-ci-local
 **/.DS_Store

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
@@ -14,7 +14,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.
 //
 
 // Construct sequence of FECollection objects.
@@ -68,35 +68,25 @@ std::vector<std::unique_ptr<FECollection>> ConstructFECollections(bool pc_pmg, b
 }
 
 // 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>
 mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
     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!");
-  auto *fespace = new mfem::ParFiniteElementSpace(mesh[0].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[0].get(), fecs[0].get(), dim);
+
   mfem::ParFiniteElementSpaceHierarchy fespaces(mesh[0].get(), fespace, false, true);
 
   // XX TODO: LibCEED transfer operators!
 
   // h-refinement
   for (std::size_t 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);
@@ -106,11 +96,8 @@ 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);
+
     auto *P = new ParOperator(
         std::make_unique<mfem::TransferOperator>(fespaces.GetFinestFESpace(), *fespace),
         fespaces.GetFinestFESpace(), *fespace, true);
@@ -119,21 +106,73 @@ mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy(
   return fespaces;
 }
 
-// Construct a single-level finite element space hierarchy from a single mesh and
-// finite element collection. Unnecessary to pass the dirichlet boundary
-// conditions as they need not be incorporated in any inter-space projectors.
+// Construct a hierarchy of finite element spaces given a single mesh and a
+// sequence of finite element collections. Uses p multigrid.
 template <typename FECollection>
 mfem::ParFiniteElementSpaceHierarchy
-ConstructFiniteElementSpaceHierarchy(mfem::ParMesh &mesh, const FECollection &fec,
-                                     const mfem::Array<int> *dbc_marker = nullptr,
-                                     mfem::Array<int> *dbc_tdof_list = nullptr)
+ConstructFiniteElementSpaceHierarchy(const std::unique_ptr<mfem::ParMesh> &mesh,
+                                     const std::vector<std::unique_ptr<FECollection>> &fecs,
+                                     int dim = 1)
+{
+  MFEM_VERIFY(!fecs.empty(), "Empty FE collection for FE space construction!");
+  auto *fespace = new mfem::ParFiniteElementSpace(mesh.get(), fecs[0].get(), dim);
+  mfem::ParFiniteElementSpaceHierarchy fespaces(mesh.get(), fespace, false, true);
+
+  // p-refinement
+  for (std::size_t l = 1; l < fecs.size(); l++)
+  {
+    fespace = new mfem::ParFiniteElementSpace(mesh.get(), fecs[l].get(), dim);
+
+    auto *P = new ParOperator(
+        std::make_unique<mfem::TransferOperator>(fespaces.GetFinestFESpace(), *fespace),
+        fespaces.GetFinestFESpace(), *fespace, true);
+    fespaces.AddLevel(mesh.get(), fespace, P, false, true, true);
+  }
+  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(
+    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)
 {
-  auto *fespace = new mfem::ParFiniteElementSpace(&mesh, &fec);
-  if (dbc_marker && dbc_tdof_list)
+  dbc_tdof_lists.clear();
+  auto fespaces = ConstructFiniteElementSpaceHierarchy(mesh, fecs, dim);
+  for (int l = 0; l < fespaces.GetNumLevels(); ++l)
   {
-    fespace->GetEssentialTrueDofs(*dbc_marker, *dbc_tdof_list);
+    fespaces.GetFESpaceAtLevel(l).GetEssentialTrueDofs(dbc_marker,
+                                                       dbc_tdof_lists.emplace_back());
   }
-  return mfem::ParFiniteElementSpaceHierarchy(&mesh, fespace, false, true);
+  return fespaces;
+}
+
+// Construct a hierarchy of finite element spaces given a single mesh and a
+// single finite element collection.
+template <typename FECollection>
+mfem::ParFiniteElementSpaceHierarchy
+ConstructFiniteElementSpaceHierarchy(mfem::ParMesh &mesh, FECollection &fec, int dim = 1)
+{
+  auto *fespace = new mfem::ParFiniteElementSpace(&mesh, &fec, dim);
+  mfem::ParFiniteElementSpaceHierarchy fespaces(&mesh, fespace, false, true);
+  return fespaces;
+}
+
+// Construct a hierarchy of finite element spaces given a single mesh and a
+// single finite element collection. Overload for treatment of Dirichlet
+// boundary conditions, extracts the true dof vector
+template <typename FECollection>
+mfem::ParFiniteElementSpaceHierarchy
+ConstructFiniteElementSpaceHierarchy(mfem::ParMesh &mesh, FECollection &fec,
+                                     const mfem::Array<int> &dbc_marker,
+                                     mfem::Array<int> &dbc_tdof_lists, int dim = 1)
+{
+  auto fespaces = ConstructFiniteElementSpaceHierarchy(mesh, fec, dim);
+  fespaces.GetFinestFESpace().GetEssentialTrueDofs(dbc_marker, dbc_tdof_lists);
+  return fespaces;
 }
 
 }  // namespace palace::utils

diff --git a/palace/linalg/ksp.cpp b/palace/linalg/ksp.cpp
@@ -61,6 +61,7 @@ std::unique_ptr<IterativeSolver<OperType>> ConfigureKrylovSolver(MPI_Comm comm,
         ksp = std::move(fgmres);
       }
       break;
+    case config::LinearSolverData::KspType::MINRES:
     case config::LinearSolverData::KspType::BICGSTAB:
     case config::LinearSolverData::KspType::DEFAULT:
     case config::LinearSolverData::KspType::INVALID:
@@ -84,12 +85,15 @@ std::unique_ptr<IterativeSolver<OperType>> ConfigureKrylovSolver(MPI_Comm comm,
       auto *gmres = static_cast<GmresSolver<OperType> *>(ksp.get());
       switch (iodata.solver.linear.pc_side_type)
       {
+        case config::LinearSolverData::SideType::DEFAULT:
         case config::LinearSolverData::SideType::LEFT:
           gmres->SetPrecSide(GmresSolver<OperType>::PrecSide::LEFT);
           break;
         case config::LinearSolverData::SideType::RIGHT:
           gmres->SetPrecSide(GmresSolver<OperType>::PrecSide::RIGHT);
           break;
+        case config::LinearSolverData::SideType::INVALID:
+          MFEM_ABORT("Only Left, Right or Default values are allowed for SideType");
       }
     }
   }
@@ -116,6 +120,8 @@ std::unique_ptr<IterativeSolver<OperType>> ConfigureKrylovSolver(MPI_Comm comm,
       case config::LinearSolverData::OrthogType::CGS2:
         gmres->SetOrthogonalization(GmresSolver<OperType>::OrthogType::CGS2);
         break;
+      case config::LinearSolverData::OrthogType::INVALID:
+        MFEM_ABORT("Only MGS, CGS or CGS2 are allowed for OrthogType");
     }
   }
 

diff --git a/palace/models/curlcurloperator.cpp b/palace/models/curlcurloperator.cpp
@@ -80,10 +80,10 @@ CurlCurlOperator::CurlCurlOperator(const IoData &iodata,
         pc_mg, false, iodata.solver.order, mesh.back()->Dimension())),
     rt_fec(iodata.solver.order - 1, mesh.back()->Dimension()),
     nd_fespaces(pc_mg ? utils::ConstructFiniteElementSpaceHierarchy(
-                            mesh, nd_fecs, &dbc_marker, &dbc_tdof_lists)
+                            mesh, nd_fecs, dbc_marker, dbc_tdof_lists)
                       : utils::ConstructFiniteElementSpaceHierarchy(
-                            *mesh.back(), *nd_fecs.back(), &dbc_marker,
-                            &dbc_tdof_lists.emplace_back())),
+                            *mesh.back(), *nd_fecs.back(), dbc_marker,
+                            dbc_tdof_lists.emplace_back())),
     h1_fespaces(pc_mg ? utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
                             mesh, h1_fecs)
                       : utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(

diff --git a/palace/models/curlcurloperator.hpp b/palace/models/curlcurloperator.hpp
@@ -70,6 +70,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,10 +121,10 @@ LaplaceOperator::LaplaceOperator(const IoData &iodata,
         pc_mg, false, iodata.solver.order, mesh.back()->Dimension())),
     nd_fec(iodata.solver.order, mesh.back()->Dimension()),
     h1_fespaces(pc_mg ? utils::ConstructFiniteElementSpaceHierarchy(
-                            mesh, h1_fecs, &dbc_marker, &dbc_tdof_lists)
+                            mesh, h1_fecs, dbc_marker, dbc_tdof_lists)
                       : utils::ConstructFiniteElementSpaceHierarchy(
-                            *mesh.back(), *h1_fecs.back(), &dbc_marker,
-                            &dbc_tdof_lists.emplace_back())),
+                            *mesh.back(), *h1_fecs.back(), dbc_marker,
+                            dbc_tdof_lists.emplace_back())),
     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
@@ -64,6 +64,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
@@ -84,15 +84,15 @@ SpaceOperator::SpaceOperator(const IoData &iodata,
         pc_mg, false, iodata.solver.order, mesh.back()->Dimension())),
     rt_fec(iodata.solver.order - 1, mesh.back()->Dimension()),
     nd_fespaces(pc_mg ? utils::ConstructFiniteElementSpaceHierarchy<mfem::ND_FECollection>(
-                            mesh, nd_fecs, &dbc_marker, &nd_dbc_tdof_lists)
+                            mesh, nd_fecs, dbc_marker, nd_dbc_tdof_lists)
                       : utils::ConstructFiniteElementSpaceHierarchy<mfem::ND_FECollection>(
-                            *mesh.back(), *nd_fecs.back(), &dbc_marker,
-                            &nd_dbc_tdof_lists.emplace_back())),
+                            *mesh.back(), *nd_fecs.back(), dbc_marker,
+                            nd_dbc_tdof_lists.emplace_back())),
     h1_fespaces(pc_mg ? utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
-                            mesh, h1_fecs, &dbc_marker, &h1_dbc_tdof_lists)
+                            mesh, h1_fecs, dbc_marker, h1_dbc_tdof_lists)
                       : utils::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
-                            *mesh.back(), *h1_fecs.back(), &dbc_marker,
-                            &h1_dbc_tdof_lists.emplace_back())),
+                            *mesh.back(), *h1_fecs.back(), dbc_marker,
+                            h1_dbc_tdof_lists.emplace_back())),
     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