Addressing initial PR feedback

- Removing MFEM exceptions
- Use long long int for true dof count
- SolveEstimateMarkRefine returns void
- Move adaptation configuration options to Refinement level
- Make dimensionalization of the mesh in PostOperator consistent

cmake/ExternalMFEM.cmake

@@ -38,7 +38,6 @@ endif()
 # required
 set(PALACE_MFEM_WITH_ZLIB NO)
 set(PALACE_MFEM_WITH_LIBUNWIND NO)
-set(PALACE_MFEM_WITH_EXCEPTIONS NO)
 find_package(ZLIB)
 if(ZLIB_FOUND)
   message(STATUS "Building MFEM with zlib support for binary output compression")
@@ -53,7 +52,6 @@ if(CMAKE_BUILD_TYPE MATCHES "Debug|debug|DEBUG")
     message(STATUS "Building MFEM with libunwind support")
     set(PALACE_MFEM_WITH_LIBUNWIND YES)
   endif()
-  set(PALACE_MFEM_WITH_EXCEPTIONS YES)
 endif()
 
 set(MFEM_OPTIONS ${PALACE_SUPERBUILD_DEFAULT_ARGS})
@@ -68,7 +66,6 @@ list(APPEND MFEM_OPTIONS
   "-DMFEM_USE_MUMPS=${PALACE_WITH_MUMPS}"
   "-DMFEM_USE_ZLIB=${PALACE_MFEM_WITH_ZLIB}"
   "-DMFEM_USE_LIBUNWIND=${PALACE_MFEM_WITH_LIBUNWIND}"
-  "-DMFEM_USE_EXCEPTIONS=${PALACE_MFEM_WITH_EXCEPTIONS}"
   "-DMFEM_USE_METIS_5=YES"
   "-DMFEM_USE_CEED=NO"
 )

docs/src/config/model.md

@@ -36,6 +36,13 @@ scale based on the bounding box of the computational domain.
 "Refinement":
 {
     "UniformLevels": <int>,
+    "AdaptTol": <float>,
+    "AdaptMaxIts": <int>,
+    "UpdateFraction": <float>,
+    "DOFLimit": <int>,
+    "MaximumImbalance": <float>,
+    "Nonconformal": <bool>,
+    "UseCoarsening": <bool>,
     "Boxes":
     [
         {
@@ -54,26 +61,30 @@ scale based on the bounding box of the computational domain.
             "Radius": float
         },
         ...
-    ],
-    "Adaptation":
-    {
-        "Tol": <float>,
-        "MaxIts": <int>,
-        "UpdateFraction": <float>,
-        "DOFLimit": <int>,
-        "SaveStep": <int>,
-        "MaximumImbalance": <float>,
-        "Nonconformal": <bool>,
-        "UseCoarsening": <bool>,
-        "WriteSerialMesh": <bool>
-    }
+    ]
 }
 ```
 
 with
 
 `"UniformLevels" [0]` :  Levels of uniform parallel mesh refinement to be performed on the
-input mesh.
+input mesh. If not performing AMR, these may be used as levels within a geometric multigrid
+scheme.
+
+`"AdaptTol" [1e-2]` : Relative error convergence tolerance for mesh adaptation.
+
+`"AdaptMaxIts" [0]` : Maximum number of iterations to perform of adaptive mesh refinement.
+
+`"UpdateFraction" [0.4]` : Dörfler marking fraction used to specify which elements to
+refine. This marking strategy will mark the smallest number of elements that make up
+`"UpdateFraction"` of the total error in the mesh. A larger value will refine more elements
+per iteration, at the cost of the final mesh being less efficient.
+
+`"DOFLimit" [0]` : The maximum allowable number of DOF within an AMR loop, if an adapted
+mesh exceeds this value no further adaptation will occur.
+
+`"Nonconformal" [false]` : Whether the adaptation should use nonconformal refinement.
+`"Nonconformal"` is necessary to enable `"UseCoarsening"`.
 
 `"Boxes"` :  Array of box region refinement objects. All elements with a node inside the box
 region will be marked for refinement.
@@ -102,39 +113,13 @@ Specified in mesh length units.
 `"Radius" [None]` : The radius of the sphere for this sphere refinement region, specified in
 mesh length units.
 
-`"Tol" [1e-3]` : Relative error convergence tolerance for mesh adaptation.
-
-`"MaxIts" [0]` : Maximum number of iterations to perform of adaptive mesh refinement.
-
-`"UpdateFraction" [0.4]` : Dörfler marking fraction used to specify which elements to
-refine. This marking strategy will mark the smallest number of elements that make up
-`"UpdateFraction"` of the total error in the mesh.
-
-`"DOFLimit" [0]` : The maximum allowable number of DOF within an AMR loop, if an adapted
-mesh exceeds this value no further adaptation will occur.
-
-`"SaveStep" [1]` : Specify which adaptive iterations to save off within the post processing
-output folder. `"SaveStep"` of 1 specifies to save all iterations, while for example
-`"SaveStep"` of 3 would save every third iteration.
-
-`"MaximumImbalance" [1.0]` : The ratio between maximum number of elements on a processor and
-minimum number of elements on a processor before a rebalance is performed. A value of `2.0`
-would result in rebalancing occurring if one processor exceeds twice the number of elements
-on another processor.
-
-`"Nonconformal" [false]` : Whether the adaptation should use nonconformal refinement.
-`"Nonconformal"` is necessary to enable `"UseCoarsening"`.
-
-`"UseCoarsening" [false]` : Whether or not to perform coarsening if the total number of DOF
-exceeds the `"DOFLimit"`. Coarsening may be useful to improve the efficiency of the mesh if
-a large value of `"UpdateFraction"` is used. Requires `"Nonconformal"` mesh refinement to be
-enabled.
-
-`"WriteSerialMesh" [true]` : Whether to write a serialized mesh to file after adaptation.
-
 ### Advanced model options
 
+  - `"MaximumImbalance" [1.0]`
   - `"MaxNCLevels" [0]`
   - `"Partition" [""]`
-  - `"ReorientTetMesh" [false]`
   - `"RemoveCurvature" [false]`
+  - `"ReorientTetMesh" [false]`
+  - `"SaveStep" [1]`
+  - `"UseCoarsening" [false]`
+  - `"WriteSerialMesh" [true]`

palace/drivers/basesolver.cpp

@@ -123,7 +123,7 @@ mfem::Array<int> MarkedElements(double threshold, const Vector &e)
 void RebalanceMesh(std::unique_ptr<mfem::ParMesh> &mesh, const IoData &iodata,
                    std::string output_dir)
 {
-  const auto &param = iodata.model.refinement.adaptation;
+  const auto &param = iodata.model.refinement;
   auto comm = mesh->GetComm();
   const int width = 1 + static_cast<int>(std::log10(Mpi::Size(comm) - 1));
   if (output_dir.back() != '/')
@@ -137,9 +137,10 @@ void RebalanceMesh(std::unique_ptr<mfem::ParMesh> &mesh, const IoData &iodata,
     if (Mpi::Size(comm) == 1)
     {
       std::ofstream serial(serial_mesh_filename);
-      iodata.DimensionalizeMesh(*mesh);
-      mesh->Mesh::Print(serial);
-      iodata.NondimensionalizeMesh(*mesh);
+      serial.precision(mesh::MSH_FLT_PRECISION);
+      mesh::DimensionalizeMesh(*mesh, iodata.GetMeshScaleFactor());
+      mesh->mfem::Mesh::Print(serial);
+      mesh::NondimensionalizeMesh(*mesh, iodata.GetMeshScaleFactor());
     }
   }
 
@@ -153,7 +154,7 @@ void RebalanceMesh(std::unique_ptr<mfem::ParMesh> &mesh, const IoData &iodata,
     Mpi::GlobalMin(1, &min_elem, comm);
     Mpi::GlobalMax(1, &max_elem, comm);
     const double ratio = double(max_elem) / min_elem;
-    Mpi::Print("Min Elem per processor: {}, Max Elem per processor: {}, Ratio: {:.3e}\n",
+    Mpi::Print("Min Elem per processor: {}, Max Elem per processor: {}, Ratio: {:.3f}\n",
                min_elem, max_elem, double(max_elem) / min_elem);
     if (ratio > param.maximum_imbalance)
     {
@@ -167,9 +168,10 @@ void RebalanceMesh(std::unique_ptr<mfem::ParMesh> &mesh, const IoData &iodata,
         if (Mpi::Root(comm))
         {
           std::ofstream serial(serial_mesh_filename);
-          iodata.DimensionalizeMesh(*mesh);
+          serial.precision(mesh::MSH_FLT_PRECISION);
+          mesh::DimensionalizeMesh(*mesh, iodata.GetMeshScaleFactor());
           mesh->Mesh::Print(serial);
-          iodata.NondimensionalizeMesh(*mesh);
+          mesh::NondimensionalizeMesh(*mesh, iodata.GetMeshScaleFactor());
         }
         Mpi::Barrier(comm);
       }
@@ -195,27 +197,29 @@ void RebalanceMesh(std::unique_ptr<mfem::ParMesh> &mesh, const IoData &iodata,
         // carefully. For NC, this requires creating a duplicate mesh.
         mfem::ParMesh smesh(*mesh);
         mfem::Array<int> serial_partition(mesh->GetNE());
-        smesh.FinalizeTopology();
-        smesh.Finalize();
-        smesh.ExchangeFaceNbrData();
         serial_partition = 0;
         smesh.Rebalance(serial_partition);
         BlockTimer bt_io(Timer::IO);
         if (Mpi::Root(comm))
         {
           std::ofstream serial(serial_mesh_filename);
-          iodata.DimensionalizeMesh(smesh);
-          smesh.Mesh::Print(serial);
+          serial.precision(mesh::MSH_FLT_PRECISION);
+          mesh::DimensionalizeMesh(smesh, iodata.GetMeshScaleFactor());
+          smesh.Mesh::Print(serial);  // Do not need to nondimensionalize the temporary
         }
         Mpi::Barrier(comm);
       }
       else
       {
-        std::ofstream serial(serial_mesh_filename);
-        serial.precision(std::numeric_limits<double>::max_digits10);
-        iodata.DimensionalizeMesh(*mesh);
-        mesh->PrintAsSerial(serial);
-        iodata.NondimensionalizeMesh(*mesh);
+        auto smesh = std::make_unique<mfem::Mesh>(mesh->GetSerialMesh(0));
+        BlockTimer bt(Timer::IO);
+        if (Mpi::Rank(comm) == 0)
+        {
+          std::ofstream serial(serial_mesh_filename);
+          serial.precision(mesh::MSH_FLT_PRECISION);
+          mesh::DimensionalizeMesh(*smesh, iodata.GetMeshScaleFactor());
+          smesh->Print(serial);  // Do not need to nondimensionalize the temporary
+        }
       }
     }
     mesh->ExchangeFaceNbrData();
@@ -224,16 +228,17 @@ void RebalanceMesh(std::unique_ptr<mfem::ParMesh> &mesh, const IoData &iodata,
 
 }  // namespace
 
-ErrorIndicator
-BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+void BaseSolver::SolveEstimateMarkRefine(
+    std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
+  const auto &param = iodata.model.refinement;
   auto comm = mesh.back()->GetComm();
   // Helper to save off postprocess data.
   auto SavePostProcess = [&](int iter)
   {
     Mpi::Barrier(comm);
     BlockTimer bt_io(Timer::IO);
-    if (Mpi::Root(comm))
+    if (Mpi::Root(comm) && param.save_adapt_iterations)
     {
       // Create a subfolder for the results of this adaptation.
       namespace fs = std::filesystem;
@@ -256,8 +261,7 @@ BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>>
     Mpi::Barrier(comm);
   };
 
-  const auto &param = iodata.model.refinement.adaptation;
-  const bool use_amr = param.max_its > 0;
+  const bool use_amr = param.adapt_max_its > 0;
   const bool use_coarsening = mesh.back()->Nonconforming() && param.use_coarsening;
   if (use_amr && mesh.size() > 1)
   {
@@ -282,7 +286,8 @@ BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>>
     else
     {
       Mpi::Print("\nAdaptive Mesh Refinement Parameters:\n");
-      Mpi::Print("MaxIter: {}, Tolerance: {:.3e}\n\n", param.max_its, param.tolerance);
+      Mpi::Print("MaxIter: {}, Tolerance: {:.3e}\n\n", param.adapt_max_its,
+                 param.adapt_tolerance);
       SavePostProcess(iter);  // Save the initial solution
     }
   }
@@ -294,11 +299,11 @@ BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>>
     // Run out of DOFs, and coarsening isn't allowed.
     ret |= (!use_coarsening && ntdof > param.dof_limit);
     // Run out of iterations.
-    ret |= iter >= param.max_its;
+    ret |= iter >= param.adapt_max_its;
     return ret;
   };
 
-  while (indicators.Norml2(comm) > param.tolerance && !exhausted_resources())
+  while (indicators.Norml2(comm) > param.adapt_tolerance && !exhausted_resources())
   {
     BlockTimer bt_adapt(Timer::ADAPTATION);
     Mpi::Print("Adaptation iteration {}: Initial Error Indicator: {:.3e}, DOF: {}, DOF "
@@ -357,15 +362,9 @@ BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>>
 
     // Solve + estimate.
     indicators_and_ntdof = Solve(mesh);
-
-    // Optionally save solution off.
-    if (param.save_step > 0 && iter % param.save_step == 0)
-    {
-      SavePostProcess(iter);
-    }
+    SavePostProcess(iter);
   }
   Mpi::Print("\nFinal Error Indicator: {:.3e}, DOF: {}\n", indicators.Norml2(comm), ntdof);
-  return indicators;
 }
 void BaseSolver::SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const
 {

palace/drivers/basesolver.hpp

@@ -86,7 +86,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, int>
+  virtual std::pair<ErrorIndicator, long long int>
   Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const = 0;
 
 public:
@@ -96,8 +96,7 @@ class BaseSolver
 
   // Performs adaptive mesh refinement using the solve-estimate-mark-refine paradigm.
   // Dispatches to the Solve method for the driver specific calculations.
-  ErrorIndicator
-  SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const;
+  void SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const;
 
   // These methods write different simulation metadata to a JSON file in post_dir.
   void SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const;

palace/drivers/drivensolver.cpp

@@ -26,7 +26,7 @@ namespace palace
 
 using namespace std::complex_literals;
 
-std::pair<ErrorIndicator, int>
+std::pair<ErrorIndicator, long long int>
 DrivenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
   // Set up the spatial discretization and frequency sweep.
@@ -98,10 +98,9 @@ DrivenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) con
   Mpi::Print("\n");
 
   // Main frequency sweep loop.
-  return std::make_pair(
-      adaptive ? SweepAdaptive(spaceop, postop, nstep, step0, omega0, delta_omega)
-               : SweepUniform(spaceop, postop, nstep, step0, omega0, delta_omega),
-      spaceop.GlobalTrueVSize());
+  return {adaptive ? SweepAdaptive(spaceop, postop, nstep, step0, omega0, delta_omega)
+                   : SweepUniform(spaceop, postop, nstep, step0, omega0, delta_omega),
+          spaceop.GlobalTrueVSize()};
 }
 
 ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop,

palace/drivers/drivensolver.hpp

@@ -60,7 +60,7 @@ class DrivenSolver : public BaseSolver
                               const WavePortOperator &wave_port_op, int step,
                               double omega) const;
 
-  std::pair<ErrorIndicator, int>
+  std::pair<ErrorIndicator, long long int>
   Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
 
 public:

palace/drivers/eigensolver.cpp

@@ -24,7 +24,7 @@ namespace palace
 
 using namespace std::complex_literals;
 
-std::pair<ErrorIndicator, int>
+std::pair<ErrorIndicator, long long int>
 EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
   // Construct and extract the system matrices defining the eigenvalue problem. The diagonal

palace/drivers/eigensolver.hpp

@@ -41,7 +41,7 @@ class EigenSolver : public BaseSolver
   void PostprocessEPR(const PostOperator &postop, const LumpedPortOperator &lumped_port_op,
                       int i, std::complex<double> omega, double Em) const;
 
-  std::pair<ErrorIndicator, int>
+  std::pair<ErrorIndicator, long long int>
   Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
 
 public:

palace/drivers/electrostaticsolver.cpp

@@ -17,7 +17,7 @@
 namespace palace
 {
 
-std::pair<ErrorIndicator, int>
+std::pair<ErrorIndicator, long long int>
 ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
   // Construct the system matrix defining the linear operator. Dirichlet boundaries are

palace/drivers/electrostaticsolver.hpp

@@ -42,7 +42,7 @@ class ElectrostaticSolver : public BaseSolver
                             const mfem::DenseMatrix &C, const mfem::DenseMatrix &Cinv,
                             const mfem::DenseMatrix &Cm) const;
 
-  std::pair<ErrorIndicator, int>
+  std::pair<ErrorIndicator, long long int>
   Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
 
 public:

palace/drivers/magnetostaticsolver.cpp

@@ -18,7 +18,7 @@
 namespace palace
 {
 
-std::pair<ErrorIndicator, int>
+std::pair<ErrorIndicator, long long int>
 MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
   // Construct the system matrix defining the linear operator. Dirichlet boundaries are

palace/drivers/magnetostaticsolver.hpp

@@ -40,7 +40,7 @@ class MagnetostaticSolver : public BaseSolver
                             const mfem::DenseMatrix &M, const mfem::DenseMatrix &Minv,
                             const mfem::DenseMatrix &Mm) const;
 
-  std::pair<ErrorIndicator, int>
+  std::pair<ErrorIndicator, long long int>
   Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
 
 public:

palace/drivers/transientsolver.cpp

@@ -20,7 +20,7 @@
 namespace palace
 {
 
-std::pair<ErrorIndicator, int>
+std::pair<ErrorIndicator, long long int>
 TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
   // Set up the spatial discretization and time integrators for the E and B fields.