Introduce Flux Error Estimation

CHANGELOG.md

@@ -42,6 +42,10 @@ The format of this changelog is based on
     based for now).
   - Added improved OpenMP support in `palace` wrapper script and CI tests.
   - Added Apptainer/Singularity container build definition for Palace.
+  - Added flux-based error estimation, reported in `error-estimate.csv`. This computes the
+    difference between the numerical gradient (electrostatics) or curl (otherwise) of the
+    solution, and a smoother approximation obtained through a global mass matrix inversion.
+    The results are reported in `error-estimates.csv` within the `"Output"` folder.
 
 ## [0.11.2] - 2023-07-14
 

docs/src/config/solver.md

@@ -417,9 +417,15 @@ iterative solver choices, and the default choice depends on the iterative solver
 `"DivFreeTol" [1.0e-12]` :  Relative tolerance for divergence-free cleaning used in the
 eigenmode simulation type.
 
-`"DivFreeMaxIts" [100]` :  Maximum number of iterations for divergence-free cleaning use in
+`"DivFreeMaxIts" [1000]` :  Maximum number of iterations for divergence-free cleaning use in
 the eigenmode simulation type.
 
+`"EstimatorTol" [1e-6]` :  Relative tolerance for flux projection used in the
+error estimate calculation.
+
+`"EstimatorMaxIts" [100]` :  Maximum number of iterations for flux projection use in
+the error estimate calculation.
+
 `"GSOrthogonalization" ["MGS"]` :  Gram-Schmidt variant used to explicitly orthogonalize
 vectors in Krylov subspace methods or other parts of the code.
 

palace/drivers/basesolver.cpp

@@ -3,9 +3,11 @@
 
 #include "basesolver.hpp"
 
+#include <array>
 #include <complex>
 #include <mfem.hpp>
 #include <nlohmann/json.hpp>
+#include "fem/errorindicator.hpp"
 #include "linalg/ksp.hpp"
 #include "models/domainpostoperator.hpp"
 #include "models/postoperator.hpp"
@@ -551,18 +553,50 @@ void BaseSolver::PostprocessProbes(const PostOperator &postop, const std::string
 #endif
 }
 
-void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double time) const
+void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double time,
+                                   const ErrorIndicator *indicator) const
 {
   // Save the computed fields in parallel in format for viewing with ParaView.
   BlockTimer bt(Timer::IO);
   if (post_dir.length() == 0)
   {
-    Mpi::Warning("No file specified under [\"Problem\"][\"Output\"]!\nSkipping saving of "
+    Mpi::Warning(postop.GetComm(),
+                 "No file specified under [\"Problem\"][\"Output\"]!\nSkipping saving of "
                  "fields to disk!\n");
     return;
   }
-  postop.WriteFields(step, time);
-  Mpi::Barrier();
+  postop.WriteFields(step, time, indicator);
+  Mpi::Barrier(postop.GetComm());
+}
+
+void BaseSolver::PostprocessErrorIndicator(const PostOperator &postop,
+                                           const ErrorIndicator &indicator) const
+{
+  // Write the indicator statistics.
+  if (post_dir.length() == 0)
+  {
+    return;
+  }
+  MPI_Comm comm = postop.GetComm();
+  std::array<double, 4> data = {indicator.Norml2(comm), indicator.Min(comm),
+                                indicator.Max(comm), indicator.Mean(comm)};
+  if (root)
+  {
+    std::string path = post_dir + "error-indicators.csv";
+    auto output = OutputFile(path, false);
+    // clang-format off
+    output.print("{:>{}s},{:>{}s},{:>{}s},{:>{}s}\n",
+                 "Norm", table.w,
+                 "Minimum", table.w,
+                 "Maximum", table.w,
+                 "Mean", table.w);
+    output.print("{:+{}.{}e},{:+{}.{}e},{:+{}.{}e},{:+{}.{}e}\n",
+                 data[0], table.w, table.p,
+                 data[1], table.w, table.p,
+                 data[2], table.w, table.p,
+                 data[3], table.w, table.p);
+    // clang-format on
+  }
 }
 
 template void BaseSolver::SaveMetadata<KspSolver>(const KspSolver &) const;

palace/drivers/basesolver.hpp

@@ -20,6 +20,7 @@ class ParMesh;
 namespace palace
 {
 
+class ErrorIndicator;
 class IoData;
 class PostOperator;
 class Timer;
@@ -57,23 +58,36 @@ class BaseSolver
                                                fmt::file::TRUNC);
   }
 
-  // Common postprocessing functions for all simulation types.
+  // Common domain postprocessing for all simulation types.
   void PostprocessDomains(const PostOperator &postop, const std::string &name, int step,
                           double time, double E_elec, double E_mag, double E_cap,
                           double E_ind) const;
+
+  // Common surface postprocessing for all simulation types.
   void PostprocessSurfaces(const PostOperator &postop, const std::string &name, int step,
                            double time, double E_elec, double E_mag, double Vinc,
                            double Iinc) const;
+
+  // Common probe postprocessing for all simulation types.
   void PostprocessProbes(const PostOperator &postop, const std::string &name, int step,
                          double time) const;
-  void PostprocessFields(const PostOperator &postop, int step, double time) const;
+
+  // Common field visualization postprocessing for all simulation types.
+  void PostprocessFields(const PostOperator &postop, int step, double time,
+                         const ErrorIndicator *indicator = nullptr) const;
+
+  // Common error indicator postprocessing for all simulation types.
+  void PostprocessErrorIndicator(const PostOperator &postop,
+                                 const ErrorIndicator &indicator) const;
 
 public:
   BaseSolver(const IoData &iodata, bool root, int size = 0, int num_thread = 0,
              const char *git_tag = nullptr);
   virtual ~BaseSolver() = default;
 
-  virtual void Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const = 0;
+  // Performs a solve using the mesh sequence, then reports error indicators.
+  virtual ErrorIndicator
+  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const = 0;
 
   // These methods write different simulation metadata to a JSON file in post_dir.
   void SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const;

palace/drivers/drivensolver.cpp

@@ -5,6 +5,8 @@
 
 #include <complex>
 #include <mfem.hpp>
+#include "fem/errorindicator.hpp"
+#include "linalg/errorestimator.hpp"
 #include "linalg/ksp.hpp"
 #include "linalg/operator.hpp"
 #include "linalg/vector.hpp"
@@ -24,7 +26,8 @@ namespace palace
 
 using namespace std::complex_literals;
 
-void DrivenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
+ErrorIndicator
+DrivenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const
 {
   // Set up the spatial discretization and frequency sweep.
   BlockTimer bt0(Timer::CONSTRUCT);
@@ -95,18 +98,13 @@ void DrivenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
   Mpi::Print("\n");
 
   // Main frequency sweep loop.
-  if (adaptive)
-  {
-    SweepAdaptive(spaceop, postop, nstep, step0, omega0, delta_omega);
-  }
-  else
-  {
-    SweepUniform(spaceop, postop, nstep, step0, omega0, delta_omega);
-  }
+  return adaptive ? SweepAdaptive(spaceop, postop, nstep, step0, omega0, delta_omega)
+                  : SweepUniform(spaceop, postop, nstep, step0, omega0, delta_omega);
 }
 
-void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, int nstep,
-                                int step0, double omega0, double delta_omega) const
+ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop,
+                                          int nstep, int step0, double omega0,
+                                          double delta_omega) const
 {
   // Construct the system matrices defining the linear operator. PEC boundaries are handled
   // simply by setting diagonal entries of the system matrix for the corresponding dofs.
@@ -138,6 +136,12 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
   E = 0.0;
   B = 0.0;
 
+  // Initialize structures for storing and reducing the results of error estimation.
+  CurlFluxErrorEstimator<ComplexVector> estimator(
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold);
+  ErrorIndicator indicator;
+
   // Main frequency sweep loop.
   int step = step0;
   double omega = omega0;
@@ -187,20 +191,26 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
                  E_elec + E_mag);
     }
 
+    // Calculate and record the error indicators.
+    estimator.AddErrorIndicator(E, indicator);
+
     // Postprocess S-parameters and optionally write solution to disk.
     Postprocess(postop, spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(),
                 spaceop.GetSurfaceCurrentOp(), step, omega, E_elec, E_mag,
-                !iodata.solver.driven.only_port_post);
+                !iodata.solver.driven.only_port_post,
+                (step == nstep - 1) ? &indicator : nullptr);
 
     // Increment frequency.
     step++;
     omega += delta_omega;
   }
   SaveMetadata(ksp);
+  return indicator;
 }
 
-void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, int nstep,
-                                 int step0, double omega0, double delta_omega) const
+ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop,
+                                           int nstep, int step0, double omega0,
+                                           double delta_omega) const
 {
   // Configure default parameters if not specified.
   BlockTimer bt0(Timer::CONSTRUCT);
@@ -238,6 +248,12 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i
   E = 0.0;
   B = 0.0;
 
+  // Initialize structures for storing and reducing the results of error estimation.
+  CurlFluxErrorEstimator<ComplexVector> estimator(
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold);
+  ErrorIndicator indicator;
+
   // Configure the PROM operator which performs the parameter space sampling and basis
   // construction during the offline phase as well as the PROM solution during the online
   // phase.
@@ -257,8 +273,10 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i
   BlockTimer bt1(Timer::CONSTRUCTPROM);
   prom.SolveHDM(omega0, E);  // Print matrix stats at first HDM solve
   prom.AddHDMSample(omega0, E);
+  estimator.AddErrorIndicator(E, indicator);
   prom.SolveHDM(omega0 + (nstep - step0 - 1) * delta_omega, E);
   prom.AddHDMSample(omega0 + (nstep - step0 - 1) * delta_omega, E);
+  estimator.AddErrorIndicator(E, indicator);
 
   // Greedy procedure for basis construction (offline phase). Basis is initialized with
   // solutions at frequency sweep endpoints.
@@ -281,6 +299,7 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i
         max_error);
     prom.SolveHDM(omega_star, E);
     prom.AddHDMSample(omega_star, E);
+    estimator.AddErrorIndicator(E, indicator);
     iter++;
   }
   Mpi::Print("\nAdaptive sampling{} {:d} frequency samples:\n"
@@ -333,12 +352,14 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i
     // Postprocess S-parameters and optionally write solution to disk.
     Postprocess(postop, spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(),
                 spaceop.GetSurfaceCurrentOp(), step, omega, E_elec, E_mag,
-                !iodata.solver.driven.only_port_post);
+                !iodata.solver.driven.only_port_post,
+                (step == nstep - 1) ? &indicator : nullptr);
 
     // Increment frequency.
     step++;
     omega += delta_omega;
   }
+  return indicator;
 }
 
 int DrivenSolver::GetNumSteps(double start, double end, double delta) const
@@ -356,7 +377,8 @@ void DrivenSolver::Postprocess(const PostOperator &postop,
                                const LumpedPortOperator &lumped_port_op,
                                const WavePortOperator &wave_port_op,
                                const SurfaceCurrentOperator &surf_j_op, int step,
-                               double omega, double E_elec, double E_mag, bool full) const
+                               double omega, double E_elec, double E_mag, bool full,
+                               const ErrorIndicator *indicator) const
 {
   // The internal GridFunctions for PostOperator have already been set from the E and B
   // solutions in the main frequency sweep loop.
@@ -379,9 +401,13 @@ void DrivenSolver::Postprocess(const PostOperator &postop,
   if (iodata.solver.driven.delta_post > 0 && step % iodata.solver.driven.delta_post == 0)
   {
     Mpi::Print("\n");
-    PostprocessFields(postop, step / iodata.solver.driven.delta_post, freq);
+    PostprocessFields(postop, step / iodata.solver.driven.delta_post, freq, indicator);
     Mpi::Print(" Wrote fields to disk at step {:d}\n", step + 1);
   }
+  if (indicator)
+  {
+    PostprocessErrorIndicator(postop, *indicator);
+  }
 }
 
 namespace

palace/drivers/drivensolver.hpp

@@ -18,6 +18,7 @@ class ParMesh;
 namespace palace
 {
 
+class ErrorIndicator;
 class IoData;
 class LumpedPortOperator;
 class PostOperator;
@@ -34,35 +35,35 @@ class DrivenSolver : public BaseSolver
 private:
   int GetNumSteps(double start, double end, double delta) const;
 
-  void SweepUniform(SpaceOperator &spaceop, PostOperator &postop, int nstep, int step0,
-                    double omega0, double delta_omega) const;
-  void SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, int nstep, int step0,
-                     double omega0, double delta_omega) const;
+  ErrorIndicator SweepUniform(SpaceOperator &spaceop, PostOperator &postop, int nstep,
+                              int step0, double omega0, double delta_omega) const;
+
+  ErrorIndicator SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, int nstep,
+                               int step0, double omega0, double delta_omega) const;
 
   void Postprocess(const PostOperator &postop, const LumpedPortOperator &lumped_port_op,
                    const WavePortOperator &wave_port_op,
                    const SurfaceCurrentOperator &surf_j_op, int step, double omega,
-                   double E_elec, double E_mag, bool full) const;
+                   double E_elec, double E_mag, bool full,
+                   const ErrorIndicator *indicator) const;
 
   void PostprocessCurrents(const PostOperator &postop,
                            const SurfaceCurrentOperator &surf_j_op, int step,
                            double omega) const;
+
   void PostprocessPorts(const PostOperator &postop,
                         const LumpedPortOperator &lumped_port_op, int step,
                         double omega) const;
+
   void PostprocessSParameters(const PostOperator &postop,
                               const LumpedPortOperator &lumped_port_op,
                               const WavePortOperator &wave_port_op, int step,
                               double omega) const;
 
 public:
-  DrivenSolver(const IoData &iodata, bool root, int size = 0, int num_thread = 0,
-               const char *git_tag = nullptr)
-    : BaseSolver(iodata, root, size, num_thread, git_tag)
-  {
-  }
+  using BaseSolver::BaseSolver;
 
-  void Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
+  ErrorIndicator Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const final;
 };
 
 }  // namespace palace