Add a folding method to error estimator class. Change to only dumping…

… the reduced form to the csv

palace/drivers/basesolver.cpp

@@ -566,8 +566,7 @@ void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double
   Mpi::Barrier();
 }
 
-void BaseSolver::PostprocessErrorIndicator(const std::string &name, int step, double time,
-                                           const std::array<double, 5> &data) const
+void BaseSolver::PostprocessErrorIndicator(const std::array<double, 5> &data) const
 {
   if (post_dir.length() == 0)
   {
@@ -578,23 +577,19 @@ void BaseSolver::PostprocessErrorIndicator(const std::string &name, int step, do
   if (root)
   {
     std::string path = post_dir + "error-indicators.csv";
-    auto output = OutputFile(path, (step > 0));
-    if (step == 0)
-    {
-      // clang-format off
-      output.print("{:>{}s},{:>{}s},{:>{}s},{:>{}s},{:>{}s},{:>{}s}\n",
-                   name, table.w1,
-                   "Sum", table.w,
-                   "Min", table.w,
-                   "Max", table.w,
-                   "Mean", table.w,
-                   "Normalization", table.w);
-      // clang-format on
-    }
+    auto output = OutputFile(path, false);
+
+    // clang-format off
+    output.print("{:>{}s},{:>{}s},{:>{}s},{:>{}s},{:>{}s}\n",
+                  "Sum", table.w,
+                  "Min", table.w,
+                  "Max", table.w,
+                  "Mean", table.w,
+                  "Normalization", table.w);
+    // clang-format on
 
     // clang-format off
-    output.print("{:{}.{}e},{:+{}.{}e},{:+{}.{}e},{:+{}.{}e},{:+{}.{}e},{:+{}.{}e}\n",
-                time, table.w1, table.p1,
+    output.print("{:+{}.{}e},{:+{}.{}e},{:+{}.{}e},{:+{}.{}e},{:+{}.{}e}\n",
                 data[0], table.w, table.p,
                 data[1], table.w, table.p,
                 data[2], table.w, table.p,

palace/drivers/basesolver.hpp

@@ -68,11 +68,7 @@ class BaseSolver
   void PostprocessProbes(const PostOperator &postop, const std::string &name, int step,
                          double time) const;
   void PostprocessFields(const PostOperator &postop, int step, double time) const;
-
-  // Postprocess granular error indicator to file. The argument normalized indicates if
-  // supplied indicators have already been normalized.
-  void PostprocessErrorIndicator(const std::string &name, int step, double time,
-                                 const std::array<double, 5> &data) const;
+  void PostprocessErrorIndicator(const std::array<double, 5> &data) const;
 
 public:
   BaseSolver(const IoData &iodata, bool root, int size = 0, int num_thread = 0,

palace/drivers/drivensolver.cpp

@@ -144,19 +144,8 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
   E = 0.0;
   B = 0.0;
 
-  // Initialize structures for storing and reducing the results of error estimation.
-  ErrorIndicator indicators;
-  auto UpdateErrorIndicator = [this, &estimator, &indicators, &postop,
-                               &spaceop](const auto &E, int step, double frequency)
-  {
-    BlockTimer bt0(Timer::ESTIMATION);
-    ErrorIndicator sample = estimator.ComputeIndicators(E);
-    postop.SetIndicatorGridFunction(sample.Local());
-    BlockTimer bt_post(Timer::POSTPRO);
-    PostprocessErrorIndicator("f (GHz)", step, frequency,
-                              sample.GetPostprocessData(spaceop.GetComm()));
-    indicators.AddIndicator(sample);
-  };
+  // Initialize structures for storing the error indicator.
+  ErrorIndicator indicator;
 
   // Main frequency sweep loop.
   double omega = omega0;
@@ -187,7 +176,7 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     ksp.Mult(RHS, E);
 
     // Compute the error indicators, and post process the indicator field.
-    UpdateErrorIndicator(E, step, freq);
+    estimator.AddErrorIndicator(indicator, postop, E);
 
     // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
     // PostOperator for all postprocessing operations.
@@ -218,7 +207,9 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     omega += delta_omega;
   }
   SaveMetadata(ksp);
-  return indicators;
+  BlockTimer bt_post(Timer::POSTPRO);
+  PostprocessErrorIndicator(indicator.GetPostprocessData(spaceop.GetComm()));
+  return indicator;
 }
 
 ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop,
@@ -276,28 +267,18 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
 
   // The error indicators will be calculated for each HDM sample rather than for
   // the online stage.
-  ErrorIndicator indicators;
-  auto UpdateErrorIndicator =
-      [this, &estimator, &indicators, &spaceop](const auto &E, int step, double frequency)
-  {
-    BlockTimer bt0(Timer::ESTIMATION);
-    auto sample_indicators = estimator.ComputeIndicators(E);
-    BlockTimer bt_post(Timer::POSTPRO);
-    PostprocessErrorIndicator("f (GHz)", step, frequency,
-                              sample_indicators.GetPostprocessData(spaceop.GetComm()));
-    indicators.AddIndicator(sample_indicators);
-  };
+  ErrorIndicator indicator;
 
   // Initialize the basis with samples from the top and bottom of the frequency
   // range of interest. Each call for an HDM solution adds the frequency sample to P_S and
   // removes it from P \ P_S. Timing for the HDM construction and solve is handled inside
   // of the RomOperator.
   BlockTimer bt1(Timer::CONSTRUCTPROM);
   prom.SolveHDM(omega0, E);  // Print matrix stats at first HDM solve
-  UpdateErrorIndicator(E, 0, omega0 * f0);
+  estimator.AddErrorIndicator(indicator, E);
   prom.AddHDMSample(omega0, E);
   prom.SolveHDM(omega0 + (nstep - step0 - 1) * delta_omega, E);
-  UpdateErrorIndicator(E, 1, (omega0 + (nstep - step0 - 1) * delta_omega) * f0);
+  estimator.AddErrorIndicator(indicator, E);
   prom.AddHDMSample(omega0 + (nstep - step0 - 1) * delta_omega, E);
 
   // Greedy procedure for basis construction (offline phase). Basis is initialized with
@@ -321,7 +302,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
         max_error);
     prom.SolveHDM(omega_star, E);
     Mpi::Print("Computing error estimates\n");
-    UpdateErrorIndicator(E, iter - iter0 + 1, omega_star * f0);
+    estimator.AddErrorIndicator(indicator, E);
     prom.AddHDMSample(omega_star, E);
     iter++;
   }
@@ -334,8 +315,10 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
              Timer::Duration(Timer::Now() - t0).count());  // Timing on root
   SaveMetadata(prom.GetLinearSolver());
 
-  // Set the indicator field to the combined field for postprocessing.
-  postop.SetIndicatorGridFunction(indicators.Local());
+  // Postprocess the indicator.
+  BlockTimer bt_post(Timer::POSTPRO);
+  postop.SetIndicatorGridFunction(indicator.Local());
+  PostprocessErrorIndicator(indicator.GetPostprocessData(spaceop.GetComm()));
 
   // Main fast frequency sweep loop (online phase).
   BlockTimer bt2(Timer::CONSTRUCT);
@@ -384,7 +367,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
     step++;
     omega += delta_omega;
   }
-  return indicators;
+  return indicator;
 }
 
 int DrivenSolver::GetNumSteps(double start, double end, double delta) const

palace/drivers/eigensolver.cpp

@@ -266,12 +266,12 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
   // Save the eigenvalue estimates.
   BlockTimer bt_est(Timer::ESTIMATION);
   std::vector<ErrorIndicator> indicators;
-  indicators.reserve(num_conv);
+  indicators.reserve(iodata.solver.eigenmode.n);
   for (int i = 0; i < iodata.solver.eigenmode.n; i++)
   {
-    eigen->GetEigenvector(i, E);
     // Only update the error indicator for targeted modes.
     Mpi::Print("\nComputing error estimates for mode {:d}\n", i);
+    eigen->GetEigenvector(i, E);
     indicators.emplace_back(estimator.ComputeIndicators(E));
   }
 
@@ -309,17 +309,16 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), omega.real());
-
-    // Postprocess the mode.
-    Postprocess(postop, spaceop.GetLumpedPortOp(), i, omega, error1, error2, num_conv);
-
-    if (i < iodata.solver.eigenmode.n)
+    if (i < iodata.solver.eigenmode.n_post)
     {
       postop.SetIndicatorGridFunction(indicators[i].Local());
-      PostprocessErrorIndicator("m", i, i + 1,
-                                indicators[i].GetPostprocessData(spaceop.GetComm()));
     }
+
+    // Postprocess the mode.
+    Postprocess(postop, spaceop.GetLumpedPortOp(), i, omega, error1, error2, num_conv);
   }
+  PostprocessErrorIndicator(
+      ErrorIndicator(indicators).GetPostprocessData(spaceop.GetComm()));
   return indicators;
 }
 

palace/drivers/electrostaticsolver.cpp

@@ -70,11 +70,9 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &me
 
     BlockTimer bt1(Timer::SOLVE);
     ksp.Mult(RHS, V[step]);
-
-    BlockTimer bt2(Timer::ESTIMATION);
     indicators[step] = estimator.ComputeIndicators(V[step]);
 
-    BlockTimer bt3(Timer::POSTPRO);
+    BlockTimer bt2(Timer::POSTPRO);
     Mpi::Print(" Sol. ||V|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(laplaceop.GetComm(), V[step]),
                linalg::Norml2(laplaceop.GetComm(), RHS));
@@ -87,6 +85,8 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &me
   BlockTimer bt1(Timer::POSTPRO);
   SaveMetadata(ksp);
   Postprocess(laplaceop, postop, V, indicators);
+  PostprocessErrorIndicator(
+      ErrorIndicator(indicators).GetPostprocessData(laplaceop.GetComm()));
   return indicators;
 }
 
@@ -123,10 +123,9 @@ void ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop, PostOperator &
     PostprocessDomains(postop, "i", i, idx, Ue, 0.0, 0.0, 0.0);
     PostprocessSurfaces(postop, "i", i, idx, Ue, 0.0, 1.0, 0.0);
     PostprocessProbes(postop, "i", i, idx);
-    PostprocessErrorIndicator("i", i, idx,
-                              indicators[i].GetPostprocessData(laplaceop.GetComm()));
     if (i < iodata.solver.electrostatic.n_post)
     {
+      postop.SetIndicatorGridFunction(indicators[i].Local());
       PostprocessFields(postop, i, idx);
       Mpi::Print(" Wrote fields to disk for terminal {:d}\n", idx);
     }

palace/drivers/magnetostaticsolver.cpp

@@ -47,7 +47,7 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &me
   MFEM_VERIFY(nstep > 0,
               "No surface current boundaries specified for magnetostatic simulation!");
 
-  // Source term and solution vector storage.
+  // Source term, solution vector storage and error indicators storage.
   Vector RHS(K->Height());
   std::vector<Vector> A(nstep);
   std::vector<ErrorIndicator> indicators(nstep);
@@ -70,11 +70,9 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &me
 
     BlockTimer bt1(Timer::SOLVE);
     ksp.Mult(RHS, A[step]);
-
-    BlockTimer bt2(Timer::ESTIMATION);
     indicators[step] = estimator.ComputeIndicators(A[step]);
 
-    BlockTimer bt3(Timer::POSTPRO);
+    BlockTimer bt2(Timer::POSTPRO);
     Mpi::Print(" Sol. ||A|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(curlcurlop.GetComm(), A[step]),
                linalg::Norml2(curlcurlop.GetComm(), RHS));
@@ -87,6 +85,8 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &me
   BlockTimer bt1(Timer::POSTPRO);
   SaveMetadata(ksp);
   Postprocess(curlcurlop, postop, A, indicators);
+  PostprocessErrorIndicator(
+      ErrorIndicator(indicators).GetPostprocessData(curlcurlop.GetComm()));
   return indicators;
 }
 
@@ -128,10 +128,9 @@ void MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, PostOperator
     PostprocessDomains(postop, "i", i, idx, 0.0, Um, 0.0, 0.0);
     PostprocessSurfaces(postop, "i", i, idx, 0.0, Um, 0.0, Iinc(i));
     PostprocessProbes(postop, "i", i, idx);
-    PostprocessErrorIndicator("i", i, idx,
-                              indicators[i].GetPostprocessData(curlcurlop.GetComm()));
     if (i < iodata.solver.magnetostatic.n_post)
     {
+      postop.SetIndicatorGridFunction(indicators[i].Local());
       PostprocessFields(postop, i, idx);
       Mpi::Print(" Wrote fields to disk for terminal {:d}\n", idx);
     }

palace/drivers/transientsolver.cpp

@@ -83,18 +83,7 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh)
     BlockTimer bt(Timer::ESTCONSTRUCT);
     return CurlFluxErrorEstimator(iodata, spaceop.GetMaterialOp(), spaceop.GetNDSpaces());
   }();
-  ErrorIndicator indicators;
-  auto UpdateErrorIndicator = [this, &estimator, &indicators, &postop,
-                               &spaceop](const auto &E, int step, double time)
-  {
-    BlockTimer bt0(Timer::ESTIMATION);
-    auto sample_indicators = estimator.ComputeIndicators(E);
-    postop.SetIndicatorGridFunction(sample_indicators.Local());
-    BlockTimer bt_post(Timer::POSTPRO);
-    PostprocessErrorIndicator("t (ns)", step, time,
-                              sample_indicators.GetPostprocessData(spaceop.GetComm()));
-    indicators.AddIndicator(sample_indicators);
-  };
+  ErrorIndicator indicator;
 
   // Main time integration loop.
   int step = 0;
@@ -138,7 +127,7 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh)
     }
 
     // Calculate and record the error indicators.
-    UpdateErrorIndicator(E, step, t);
+    estimator.AddErrorIndicator(indicator, postop, E);
 
     // Postprocess port voltages/currents and optionally write solution to disk.
     Postprocess(postop, spaceop.GetLumpedPortOp(), spaceop.GetSurfaceCurrentOp(), step, t,
@@ -148,7 +137,8 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh)
     step++;
   }
   SaveMetadata(timeop.GetLinearSolver());
-  return indicators;
+  PostprocessErrorIndicator(indicator.GetPostprocessData(spaceop.GetComm()));
+  return indicator;
 }
 std::function<double(double)> TransientSolver::GetTimeExcitation(bool dot) const
 {

palace/linalg/errorestimator.cpp

@@ -12,6 +12,7 @@
 #include "linalg/iterative.hpp"
 #include "linalg/rap.hpp"
 #include "models/materialoperator.hpp"
+#include "models/postoperator.hpp"
 #include "utils/communication.hpp"
 #include "utils/iodata.hpp"
 #include "utils/timer.hpp"
@@ -111,6 +112,7 @@ CurlFluxErrorEstimator::CurlFluxErrorEstimator(
 template <>
 ErrorIndicator CurlFluxErrorEstimator::ComputeIndicators(const ComplexVector &v) const
 {
+  BlockTimer bt_est(Timer::ESTIMATION);
   auto &nd_fespace = nd_fespaces.GetFinestFESpace();
   const int nelem = nd_fespace.GetNE();
 
@@ -176,6 +178,7 @@ ErrorIndicator CurlFluxErrorEstimator::ComputeIndicators(const ComplexVector &v)
 template <>
 ErrorIndicator CurlFluxErrorEstimator::ComputeIndicators(const Vector &v) const
 {
+  BlockTimer bt_est(Timer::ESTIMATION);
   auto &nd_fespace = nd_fespaces.GetFinestFESpace();
   field_gf.SetFromTrueDofs(v);
   const int nelem = nd_fespace.GetNE();
@@ -234,6 +237,37 @@ ErrorIndicator CurlFluxErrorEstimator::ComputeIndicators(const Vector &v) const
   return ErrorIndicator(std::move(estimates), normalization);
 }
 
+template <typename VectorType>
+void CurlFluxErrorEstimator::AddErrorIndicator(ErrorIndicator &indicator,
+                                               PostOperator &postop,
+                                               const VectorType &v) const
+{
+  auto i = ComputeIndicators(v);
+  BlockTimer bt(Timer::POSTPRO);
+  postop.SetIndicatorGridFunction(i.Local());
+  indicator.AddIndicator(i);
+}
+
+template <typename VectorType>
+void CurlFluxErrorEstimator::AddErrorIndicator(ErrorIndicator &indicator,
+                                               const VectorType &v) const
+{
+  BlockTimer bt(Timer::POSTPRO);
+  indicator.AddIndicator(ComputeIndicators(v));
+}
+
+template void
+CurlFluxErrorEstimator::AddErrorIndicator<ComplexVector>(ErrorIndicator &, PostOperator &,
+                                                         const ComplexVector &) const;
+template void CurlFluxErrorEstimator::AddErrorIndicator<Vector>(ErrorIndicator &,
+                                                                PostOperator &,
+                                                                const Vector &) const;
+template void
+CurlFluxErrorEstimator::AddErrorIndicator<ComplexVector>(ErrorIndicator &,
+                                                         const ComplexVector &) const;
+template void CurlFluxErrorEstimator::AddErrorIndicator<Vector>(ErrorIndicator &,
+                                                                const Vector &) const;
+
 GradFluxErrorEstimator::GradFluxErrorEstimator(
     const IoData &iodata, const MaterialOperator &mat_op,
     mfem::ParFiniteElementSpaceHierarchy &h1_fespaces)
@@ -250,6 +284,7 @@ GradFluxErrorEstimator::GradFluxErrorEstimator(
 
 ErrorIndicator GradFluxErrorEstimator::ComputeIndicators(const Vector &v) const
 {
+  BlockTimer bt_est(Timer::ESTIMATION);
   auto &h1_fespace = h1_fespaces.GetFinestFESpace();
   const int sdim = h1_fespace.GetMesh()->SpaceDimension();
   field_gf.SetFromTrueDofs(v);