Suggested improvements for #85

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,15 +3,16 @@
 
 #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"
 #include "models/surfacepostoperator.hpp"
 #include "utils/communication.hpp"
-#include "utils/errorindicators.hpp"
 #include "utils/filesystem.hpp"
 #include "utils/iodata.hpp"
 #include "utils/timer.hpp"
@@ -552,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

@@ -11,14 +11,16 @@
 
 namespace mfem
 {
+
 class ParFiniteElementSpaceHierarchy;
 class ParMesh;
+
 }  // namespace mfem
 
 namespace palace
 {
 
-struct ErrorIndicators;
+class ErrorIndicator;
 class IoData;
 class PostOperator;
 class Timer;
@@ -56,24 +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;
 
-  // Performs a solve using the mesh sequence, and report error indicators.
-  virtual ErrorIndicators Solve(const 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;