Merge pull request #166 from awslabs/sjg/mat-coeff-revisions

`palace::Mesh` and libCEED quadrature data refactor

.clang-format

@@ -52,6 +52,7 @@ SpacesBeforeTrailingComments: 2
 StatementMacros: ['PalacePragmaOmp',
                   'PalacePragmaDiagnosticPush',
                   'PalacePragmaDiagnosticPop',
-                  'PalacePragmaDiagnosticDisableDeprecated']
+                  'PalacePragmaDiagnosticDisableDeprecated',
+                  'PalacePragmaDiagnosticDisableUnused']
 TypenameMacros: ['CEED_QFUNCTION']
 UseTab: Never

docs/src/config/solver.md

@@ -466,7 +466,7 @@ vectors in Krylov subspace methods or other parts of the code.
   - `"MGSmoothEigScaleMin" [0.0]`
   - `"MGSmoothChebyshev4th" [true]`
   - `"ColumnOrdering" ["Default"]` :  `"METIS"`, `"ParMETIS"`,`"Scotch"`, `"PTScotch"`,
-    `"Default"`
+    `"PORD"`, `"AMD"`, `"RCM"`, `"Default"`
   - `"STRUMPACKCompressionType" ["None"]` :  `"None"`, `"BLR"`, `"HSS"`, `"HODLR"`, `"ZFP"`,
     `"BLR-HODLR"`, `"ZFP-BLR-HODLR"`
   - `"STRUMPACKCompressionTol" [1.0e-3]`

palace/drivers/basesolver.cpp

@@ -11,6 +11,7 @@
 #include "drivers/transientsolver.hpp"
 #include "fem/errorindicator.hpp"
 #include "fem/fespace.hpp"
+#include "fem/mesh.hpp"
 #include "linalg/ksp.hpp"
 #include "models/domainpostoperator.hpp"
 #include "models/postoperator.hpp"
@@ -136,8 +137,7 @@ BaseSolver::BaseSolver(const IoData &iodata, bool root, int size, int num_thread
   }
 }
 
-void BaseSolver::SolveEstimateMarkRefine(
-    std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+void BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<Mesh>> &mesh) const
 {
   const auto &refinement = iodata.model.refinement;
   const bool use_amr = [&]()
@@ -155,7 +155,7 @@ void BaseSolver::SolveEstimateMarkRefine(
                "the sequence of a priori refinements\n");
     mesh.erase(mesh.begin(), mesh.end() - 1);
     constexpr bool refine = true, fix_orientation = true;
-    mesh.back()->Finalize(refine, fix_orientation);
+    mesh.back()->Get().Finalize(refine, fix_orientation);
   }
   MPI_Comm comm = mesh.back()->GetComm();
 
@@ -206,7 +206,7 @@ void BaseSolver::SolveEstimateMarkRefine(
         refinement.update_fraction);
 
     // Refine.
-    auto &fine_mesh = *mesh.back();
+    mfem::ParMesh &fine_mesh = *mesh.back();
     const auto initial_elem_count = fine_mesh.GetGlobalNE();
     fine_mesh.GeneralRefinement(marked_elements, -1, refinement.max_nc_levels);
     const auto final_elem_count = fine_mesh.GetGlobalNE();
@@ -227,6 +227,7 @@ void BaseSolver::SolveEstimateMarkRefine(
                  "(new ratio = {:.3f})\n",
                  ratio_pre, refinement.maximum_imbalance, ratio_post);
     }
+    mesh.back()->Update();
 
     // Solve + estimate.
     Mpi::Print("\nProceeding with solve/estimate iteration {}...\n", it + 1);
@@ -249,7 +250,7 @@ void BaseSolver::SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const
     return;
   }
   const auto &fespace = fespaces.GetFinestFESpace();
-  HYPRE_BigInt ne = fespace.GetParMesh()->GetNE();
+  HYPRE_BigInt ne = fespace.GetParMesh().GetNE();
   Mpi::GlobalSum(1, &ne, fespace.GetComm());
   std::vector<HYPRE_BigInt> ndofs(fespaces.GetNumLevels());
   for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++)

palace/drivers/basesolver.hpp

@@ -9,18 +9,12 @@
 #include <vector>
 #include <fmt/os.h>
 
-namespace mfem
-{
-
-class ParMesh;
-
-}  // namespace mfem
-
 namespace palace
 {
 
 class ErrorIndicator;
 class FiniteElementSpaceHierarchy;
+class Mesh;
 class IoData;
 class PostOperator;
 class Timer;
@@ -83,7 +77,7 @@ class BaseSolver
   // Performs a solve using the mesh sequence, then reports error indicators and the number
   // of global true dofs.
   virtual std::pair<ErrorIndicator, long long int>
-  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const = 0;
+  Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const = 0;
 
 public:
   BaseSolver(const IoData &iodata, bool root, int size = 0, int num_thread = 0,
@@ -92,7 +86,7 @@ class BaseSolver
 
   // Performs adaptive mesh refinement using the solve-estimate-mark-refine paradigm.
   // Dispatches to the Solve method for the driver specific calculations.
-  void SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const;
+  void SolveEstimateMarkRefine(std::vector<std::unique_ptr<Mesh>> &mesh) const;
 
   // These methods write different simulation metadata to a JSON file in post_dir.
   void SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const;

palace/drivers/drivensolver.cpp

@@ -6,6 +6,7 @@
 #include <complex>
 #include <mfem.hpp>
 #include "fem/errorindicator.hpp"
+#include "fem/mesh.hpp"
 #include "linalg/errorestimator.hpp"
 #include "linalg/ksp.hpp"
 #include "linalg/operator.hpp"
@@ -27,7 +28,7 @@ namespace palace
 using namespace std::complex_literals;
 
 std::pair<ErrorIndicator, long long int>
-DrivenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+DrivenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 {
   // Set up the spatial discretization and frequency sweep.
   BlockTimer bt0(Timer::CONSTRUCT);
@@ -54,7 +55,7 @@ DrivenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) con
     bool first = true;
     for (const auto &[idx, data] : spaceop.GetLumpedPortOp())
     {
-      if (data.IsExcited())
+      if (data.excitation)
       {
         if (first)
         {
@@ -69,7 +70,7 @@ DrivenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) con
     first = true;
     for (const auto &[idx, data] : spaceop.GetWavePortOp())
     {
-      if (data.IsExcited())
+      if (data.excitation)
       {
         if (first)
         {
@@ -508,7 +509,7 @@ void DrivenSolver::PostprocessPorts(const PostOperator &postop,
     const double Iinc = (std::abs(Vinc) > 0.0) ? data.GetExcitationPower() / Vinc : 0.0;
     const std::complex<double> Vi = postop.GetPortVoltage(lumped_port_op, idx);
     const std::complex<double> Ii = postop.GetPortCurrent(lumped_port_op, idx);
-    port_data.push_back({idx, data.IsExcited(),
+    port_data.push_back({idx, data.excitation,
                          iodata.DimensionalizeValue(IoData::ValueType::VOLTAGE, Vinc),
                          iodata.DimensionalizeValue(IoData::ValueType::CURRENT, Iinc),
                          iodata.DimensionalizeValue(IoData::ValueType::VOLTAGE, Vi),
@@ -642,7 +643,7 @@ void DrivenSolver::PostprocessSParameters(const PostOperator &postop,
   int source_idx = -1;
   for (const auto &[idx, data] : lumped_port_op)
   {
-    if (data.IsExcited())
+    if (data.excitation)
     {
       if (src_lumped_port || src_wave_port)
       {
@@ -654,7 +655,7 @@ void DrivenSolver::PostprocessSParameters(const PostOperator &postop,
   }
   for (const auto &[idx, data] : wave_port_op)
   {
-    if (data.IsExcited())
+    if (data.excitation)
     {
       if (src_lumped_port || src_wave_port)
       {

palace/drivers/drivensolver.hpp

@@ -8,19 +8,13 @@
 #include <vector>
 #include "drivers/basesolver.hpp"
 
-namespace mfem
-{
-
-class ParMesh;
-
-}  // namespace mfem
-
 namespace palace
 {
 
 class ErrorIndicator;
 class IoData;
 class LumpedPortOperator;
+class Mesh;
 class PostOperator;
 class SpaceOperator;
 class SurfaceCurrentOperator;
@@ -61,7 +55,7 @@ class DrivenSolver : public BaseSolver
                               double omega) const;
 
   std::pair<ErrorIndicator, long long int>
-  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
+  Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const override;
 
 public:
   using BaseSolver::BaseSolver;

palace/drivers/eigensolver.cpp

@@ -5,6 +5,7 @@
 
 #include <mfem.hpp>
 #include "fem/errorindicator.hpp"
+#include "fem/mesh.hpp"
 #include "linalg/arpack.hpp"
 #include "linalg/divfree.hpp"
 #include "linalg/errorestimator.hpp"
@@ -25,7 +26,7 @@ namespace palace
 using namespace std::complex_literals;
 
 std::pair<ErrorIndicator, long long int>
-EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 {
   // Construct and extract the system matrices defining the eigenvalue problem. The diagonal
   // values for the mass matrix PEC dof shift the Dirichlet eigenvalues out of the
@@ -542,7 +543,7 @@ void EigenSolver::PostprocessEPR(const PostOperator &postop,
   epr_L_data.reserve(lumped_port_op.Size());
   for (const auto &[idx, data] : lumped_port_op)
   {
-    if (std::abs(data.GetL()) > 0.0)
+    if (std::abs(data.L) > 0.0)
     {
       const double pj = postop.GetInductorParticipation(lumped_port_op, idx, Em);
       epr_L_data.push_back({idx, pj});
@@ -582,7 +583,7 @@ void EigenSolver::PostprocessEPR(const PostOperator &postop,
   epr_IO_data.reserve(lumped_port_op.Size());
   for (const auto &[idx, data] : lumped_port_op)
   {
-    if (std::abs(data.GetR()) > 0.0)
+    if (std::abs(data.R) > 0.0)
     {
       const double Kl = postop.GetExternalKappa(lumped_port_op, idx, Em);
       const double Ql = (Kl == 0.0) ? mfem::infinity() : omega.real() / std::abs(Kl);

palace/drivers/eigensolver.hpp

@@ -9,19 +9,13 @@
 #include <vector>
 #include "drivers/basesolver.hpp"
 
-namespace mfem
-{
-
-class ParMesh;
-
-}  // namespace mfem
-
 namespace palace
 {
 
 class ErrorIndicator;
 class IoData;
 class LumpedPortOperator;
+class Mesh;
 class PostOperator;
 class Timer;
 
@@ -45,7 +39,7 @@ class EigenSolver : public BaseSolver
                       int i, std::complex<double> omega, double Em) const;
 
   std::pair<ErrorIndicator, long long int>
-  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
+  Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const override;
 
 public:
   using BaseSolver::BaseSolver;

palace/drivers/electrostaticsolver.cpp

@@ -5,6 +5,7 @@
 
 #include <mfem.hpp>
 #include "fem/errorindicator.hpp"
+#include "fem/mesh.hpp"
 #include "linalg/errorestimator.hpp"
 #include "linalg/ksp.hpp"
 #include "linalg/operator.hpp"
@@ -18,7 +19,7 @@ namespace palace
 {
 
 std::pair<ErrorIndicator, long long int>
-ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 {
   // Construct the system matrix defining the linear operator. Dirichlet boundaries are
   // handled eliminating the rows and columns of the system matrix for the corresponding

palace/drivers/electrostaticsolver.hpp

@@ -16,7 +16,6 @@ namespace mfem
 template <typename T>
 class Array;
 class DenseMatrix;
-class ParMesh;
 
 }  // namespace mfem
 
@@ -26,6 +25,7 @@ namespace palace
 class ErrorIndicator;
 class IoData;
 class LaplaceOperator;
+class Mesh;
 class PostOperator;
 class Timer;
 
@@ -43,7 +43,7 @@ class ElectrostaticSolver : public BaseSolver
                             const mfem::DenseMatrix &Cm) const;
 
   std::pair<ErrorIndicator, long long int>
-  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
+  Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const override;
 
 public:
   using BaseSolver::BaseSolver;

palace/drivers/magnetostaticsolver.cpp

@@ -5,6 +5,7 @@
 
 #include <mfem.hpp>
 #include "fem/errorindicator.hpp"
+#include "fem/mesh.hpp"
 #include "linalg/errorestimator.hpp"
 #include "linalg/ksp.hpp"
 #include "linalg/operator.hpp"
@@ -19,7 +20,7 @@ namespace palace
 {
 
 std::pair<ErrorIndicator, long long int>
-MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 {
   // Construct the system matrix defining the linear operator. Dirichlet boundaries are
   // handled eliminating the rows and columns of the system matrix for the corresponding

palace/drivers/magnetostaticsolver.hpp

@@ -13,7 +13,6 @@ namespace mfem
 {
 
 class DenseMatrix;
-class ParMesh;
 
 }  // namespace mfem
 
@@ -23,6 +22,7 @@ namespace palace
 class CurlCurlOperator;
 class ErrorIndicator;
 class IoData;
+class Mesh;
 class PostOperator;
 class SurfaceCurrentOperator;
 class Timer;
@@ -41,7 +41,7 @@ class MagnetostaticSolver : public BaseSolver
                             const mfem::DenseMatrix &Mm) const;
 
   std::pair<ErrorIndicator, long long int>
-  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
+  Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const override;
 
 public:
   using BaseSolver::BaseSolver;

palace/drivers/transientsolver.cpp

@@ -5,6 +5,7 @@
 
 #include <mfem.hpp>
 #include "fem/errorindicator.hpp"
+#include "fem/mesh.hpp"
 #include "linalg/errorestimator.hpp"
 #include "linalg/vector.hpp"
 #include "models/lumpedportoperator.hpp"
@@ -21,7 +22,7 @@ namespace palace
 {
 
 std::pair<ErrorIndicator, long long int>
-TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 {
   // Set up the spatial discretization and time integrators for the E and B fields.
   BlockTimer bt0(Timer::CONSTRUCT);
@@ -49,7 +50,7 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh)
     bool first = true;
     for (const auto &[idx, data] : spaceop.GetLumpedPortOp())
     {
-      if (data.IsExcited())
+      if (data.excitation)
       {
         if (first)
         {
@@ -367,7 +368,7 @@ void TransientSolver::PostprocessPorts(const PostOperator &postop,
         (std::abs(Vinc) > 0.0) ? data.GetExcitationPower() * J_coef * J_coef / Vinc : 0.0;
     const double Vi = postop.GetPortVoltage(lumped_port_op, idx).real();
     const double Ii = postop.GetPortCurrent(lumped_port_op, idx).real();
-    port_data.push_back({idx, data.IsExcited(),
+    port_data.push_back({idx, data.excitation,
                          iodata.DimensionalizeValue(IoData::ValueType::VOLTAGE, Vinc),
                          iodata.DimensionalizeValue(IoData::ValueType::CURRENT, Iinc),
                          iodata.DimensionalizeValue(IoData::ValueType::VOLTAGE, Vi),