Merge pull request #234 from awslabs/sjg/error-estimator-dev-prep

Fix Paraview output for error indicator

docs/src/config/solver.md

@@ -170,7 +170,6 @@ solver. This option is relevant only for `"Type": "FEAST"`.
     "MaxFreq": <float>,
     "FreqStep": <float>,
     "SaveStep": <int>,
-    "SaveOnlyPorts": <bool>,
     "Restart": <int>,
     "AdaptiveTol": <float>,
     "AdaptiveMaxSamples": <int>,
@@ -191,11 +190,6 @@ fields to disk for visualization with [ParaView](https://www.paraview.org/). Fil
 saved in the `paraview/` directory under the directory specified by
 [`config["Problem"]["Output"]`](problem.md#config%5B%22Problem%22%5D).
 
-`"SaveOnlyPorts" [false]` :  If set to `true`, postprocessing is only performed for port
-boundaries and skipped for quantities depending on, for example, field integrals over all
-or part of the interior of the computational domain. This can be useful in speeding up
-simulations if only port boundary quantities are required.
-
 `"Restart" [1]` :  Iteration (1-based) from which to restart for a partial frequency sweep
 simulation. That is, the initial frequency will be computed as
 `"MinFreq" + ("Restart" - 1) * "FreqStep"`.
@@ -226,8 +220,7 @@ error tolerance.
     "ExcitationWidth": <float>,
     "MaxTime": <float>,
     "TimeStep": <float>,
-    "SaveStep": <int>,
-    "SaveOnlyPorts": <bool>
+    "SaveStep": <int>
 }
 ```
 
@@ -279,11 +272,6 @@ disk for visualization with [ParaView](https://www.paraview.org/). Files are sav
 `paraview/` directory under the directory specified by
 [`config["Problem"]["Output"]`](problem.md#config%5B%22Problem%22%5D).
 
-`"SaveOnlyPorts" [false]` :  If set to `true`, postprocessing is only performed for port
-boundaries and skipped for quantities depending on, for example, field integrals over all
-or part of the interior of the computational domain. This can be useful in speeding up
-simulations if only port boundary quantities are required.
-
 ## `solver["Electrostatic"]`
 
 ```json

palace/drivers/basesolver.cpp

@@ -736,8 +736,7 @@ void BaseSolver::PostprocessProbes(const PostOperator &postop, const std::string
 #endif
 }
 
-void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double time,
-                                   const ErrorIndicator *indicator) const
+void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double time) const
 {
   // Save the computed fields in parallel in format for viewing with ParaView.
   BlockTimer bt(Timer::IO);
@@ -748,12 +747,13 @@ void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double
                  "fields to disk!\n");
     return;
   }
-  postop.WriteFields(step, time, indicator);
+  postop.WriteFields(step, time);
   Mpi::Barrier(postop.GetComm());
 }
 
 void BaseSolver::PostprocessErrorIndicator(const PostOperator &postop,
-                                           const ErrorIndicator &indicator) const
+                                           const ErrorIndicator &indicator,
+                                           bool fields) const
 {
   // Write the indicator statistics.
   if (post_dir.length() == 0)
@@ -780,6 +780,12 @@ void BaseSolver::PostprocessErrorIndicator(const PostOperator &postop,
                  data[3], table.w, table.p);
     // clang-format on
   }
+  if (fields)
+  {
+    BlockTimer bt(Timer::IO);
+    postop.WriteFieldsFinal(&indicator);
+    Mpi::Barrier(postop.GetComm());
+  }
 }
 
 template void BaseSolver::SaveMetadata<KspSolver>(const KspSolver &) const;

palace/drivers/basesolver.hpp

@@ -67,12 +67,11 @@ class BaseSolver
                          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;
+  void PostprocessFields(const PostOperator &postop, int step, double time) const;
 
   // Common error indicator postprocessing for all simulation types.
   void PostprocessErrorIndicator(const PostOperator &postop,
-                                 const ErrorIndicator &indicator) const;
+                                 const ErrorIndicator &indicator, bool fields) const;
 
   // Performs a solve using the mesh sequence, then reports error indicators and the number
   // of global true dofs.

palace/drivers/drivensolver.cpp

@@ -175,21 +175,19 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
     // PostOperator for all postprocessing operations.
     BlockTimer bt0(Timer::POSTPRO);
-    double E_elec = 0.0, E_mag = 0.0;
     Curl.Mult(E.Real(), B.Real());
     Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
+    double E_elec = postop.GetEFieldEnergy();
+    double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(spaceop.GetComm(), E),
                linalg::Norml2(spaceop.GetComm(), RHS));
-    if (!iodata.solver.driven.only_port_post)
     {
       const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
-      E_elec = postop.GetEFieldEnergy();
-      E_mag = postop.GetHFieldEnergy();
       Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", E_elec * J,
                  E_mag * J, (E_elec + E_mag) * J);
     }
@@ -201,7 +199,6 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     // 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,
                 (step == nstep - 1) ? &indicator : nullptr);
 
     // Increment frequency.
@@ -261,12 +258,16 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
   // 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.
+  auto UpdatePROM = [&](double omega)
+  {
+    // Add the HDM solution to the PROM reduced basis.
+    promop.UpdatePROM(omega, E);
+    estimator.AddErrorIndicator(E, indicator);
+  };
   promop.SolveHDM(omega0, E);
-  promop.UpdatePROM(omega0, E);
-  estimator.AddErrorIndicator(E, indicator);
+  UpdatePROM(omega0);
   promop.SolveHDM(omega0 + (nstep - step0 - 1) * delta_omega, E);
-  promop.UpdatePROM(omega0 + (nstep - step0 - 1) * delta_omega, E);
-  estimator.AddErrorIndicator(E, indicator);
+  UpdatePROM(omega0 + (nstep - step0 - 1) * delta_omega);
 
   // Greedy procedure for basis construction (offline phase). Basis is initialized with
   // solutions at frequency sweep endpoints.
@@ -308,8 +309,7 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
                (memory == 0)
                    ? ""
                    : fmt::format(", memory = {:d}/{:d}", memory, convergence_memory));
-    promop.UpdatePROM(omega_star, E);
-    estimator.AddErrorIndicator(E, indicator);
+    UpdatePROM(omega_star);
     it++;
   }
   Mpi::Print("\nAdaptive sampling{} {:d} frequency samples:\n"
@@ -343,27 +343,24 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
     // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
     // PostOperator for all postprocessing operations.
     BlockTimer bt0(Timer::POSTPRO);
-    double E_elec = 0.0, E_mag = 0.0;
     Curl.Mult(E.Real(), B.Real());
     Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
+    double E_elec = postop.GetEFieldEnergy();
+    double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(spaceop.GetComm(), E));
-    if (!iodata.solver.driven.only_port_post)
     {
       const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
-      E_elec = postop.GetEFieldEnergy();
-      E_mag = postop.GetHFieldEnergy();
       Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", E_elec * J,
                  E_mag * J, (E_elec + E_mag) * J);
     }
 
     // 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,
                 (step == nstep - 1) ? &indicator : nullptr);
 
     // Increment frequency.
@@ -394,36 +391,33 @@ 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,
+                               double omega, double E_elec, double E_mag,
                                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.
-  double freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
+  const double freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
+  double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
+  double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessCurrents(postop, surf_j_op, step, omega);
   PostprocessPorts(postop, lumped_port_op, step, omega);
   if (surf_j_op.Size() == 0)
   {
     PostprocessSParameters(postop, lumped_port_op, wave_port_op, step, omega);
   }
-  if (full)
-  {
-    double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
-    double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
-    PostprocessDomains(postop, "f (GHz)", step, freq, E_elec, E_mag, E_cap, E_ind);
-    PostprocessSurfaces(postop, "f (GHz)", step, freq, E_elec + E_cap, E_mag + E_ind, 1.0,
-                        1.0);
-    PostprocessProbes(postop, "f (GHz)", step, freq);
-  }
+  PostprocessDomains(postop, "f (GHz)", step, freq, E_elec, E_mag, E_cap, E_ind);
+  PostprocessSurfaces(postop, "f (GHz)", step, freq, E_elec + E_cap, E_mag + E_ind, 1.0,
+                      1.0);
+  PostprocessProbes(postop, "f (GHz)", step, freq);
   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, indicator);
+    PostprocessFields(postop, step / iodata.solver.driven.delta_post, freq);
     Mpi::Print(" Wrote fields to disk at step {:d}\n", step + 1);
   }
   if (indicator)
   {
-    PostprocessErrorIndicator(postop, *indicator);
+    PostprocessErrorIndicator(postop, *indicator, iodata.solver.driven.delta_post > 0);
   }
 }
 

palace/drivers/drivensolver.hpp

@@ -12,13 +12,11 @@ namespace palace
 {
 
 class ErrorIndicator;
-class IoData;
 class LumpedPortOperator;
 class Mesh;
 class PostOperator;
 class SpaceOperator;
 class SurfaceCurrentOperator;
-class Timer;
 class WavePortOperator;
 
 //
@@ -38,8 +36,7 @@ class DrivenSolver : public BaseSolver
   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 ErrorIndicator *indicator) const;
+                   double E_elec, double E_mag, const ErrorIndicator *indicator) const;
 
   void PostprocessCurrents(const PostOperator &postop,
                            const SurfaceCurrentOperator &surf_j_op, int step,

palace/drivers/eigensolver.cpp

@@ -206,8 +206,11 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   // Configure the shift-and-invert strategy is employed to solve for the eigenvalues
   // closest to the specified target, σ.
   const double target = iodata.solver.eigenmode.target;
-  const double f_target = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, target);
-  Mpi::Print(" Shift-and-invert σ = {:.3e} GHz ({:.3e})\n", f_target, target);
+  {
+    const double f_target =
+        iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, target);
+    Mpi::Print(" Shift-and-invert σ = {:.3e} GHz ({:.3e})\n", f_target, target);
+  }
   if (C)
   {
     // Search for eigenvalues closest to λ = iσ.
@@ -270,8 +273,8 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   BlockTimer bt2(Timer::POSTPRO);
   SaveMetadata(*ksp);
 
-  // Calculate and record the error indicators.
-  Mpi::Print("\nComputing solution error estimates\n");
+  // Calculate and record the error indicators, and postprocess the results.
+  Mpi::Print("\nComputing solution error estimates and performing postprocessing\n");
   CurlFluxErrorEstimator<ComplexVector> estimator(
       spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
       iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
@@ -284,13 +287,6 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     eigen->SetBMat(*KM);
     eigen->RescaleEigenvectors(num_conv);
   }
-  for (int i = 0; i < iodata.solver.eigenmode.n; i++)
-  {
-    eigen->GetEigenvector(i, E);
-    estimator.AddErrorIndicator(E, indicator);
-  }
-
-  // Postprocess the results.
   Mpi::Print("\n");
   for (int i = 0; i < num_conv; i++)
   {
@@ -318,27 +314,31 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), omega.real());
+    double E_elec = postop.GetEFieldEnergy();
+    double E_mag = postop.GetHFieldEnergy();
+
+    // Calculate and record the error indicators.
+    if (i < iodata.solver.eigenmode.n)
+    {
+      estimator.AddErrorIndicator(E, indicator);
+    }
 
     // Postprocess the mode.
     Postprocess(postop, spaceop.GetLumpedPortOp(), i, omega, error_bkwd, error_abs,
-                num_conv, (i == 0) ? &indicator : nullptr);
+                num_conv, E_elec, E_mag,
+                (i == iodata.solver.eigenmode.n - 1) ? &indicator : nullptr);
   }
   return {indicator, spaceop.GlobalTrueVSize()};
 }
 
 void EigenSolver::Postprocess(const PostOperator &postop,
                               const LumpedPortOperator &lumped_port_op, int i,
                               std::complex<double> omega, double error_bkwd,
-                              double error_abs, int num_conv,
+                              double error_abs, int num_conv, double E_elec, double E_mag,
                               const ErrorIndicator *indicator) const
 {
   // The internal GridFunctions for PostOperator have already been set from the E and B
-  // solutions in the main loop over converged eigenvalues. Note: The energies output are
-  // nondimensional (they can be dimensionalized using the scaling μ₀ * H₀² * L₀³, which
-  // are the free space permeability, characteristic magnetic field strength, and
-  // characteristic length scale, respectively).
-  double E_elec = postop.GetEFieldEnergy();
-  double E_mag = postop.GetHFieldEnergy();
+  // solutions in the main loop over converged eigenvalues.
   double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
   double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessEigen(i, omega, error_bkwd, error_abs, num_conv);
@@ -349,12 +349,12 @@ void EigenSolver::Postprocess(const PostOperator &postop,
   PostprocessProbes(postop, "m", i, i + 1);
   if (i < iodata.solver.eigenmode.n_post)
   {
-    PostprocessFields(postop, i, i + 1, indicator);
+    PostprocessFields(postop, i, i + 1);
     Mpi::Print(" Wrote mode {:d} to disk\n", i + 1);
   }
   if (indicator)
   {
-    PostprocessErrorIndicator(postop, *indicator);
+    PostprocessErrorIndicator(postop, *indicator, iodata.solver.eigenmode.n_post > 0);
   }
 }
 

palace/drivers/eigensolver.hpp

@@ -13,11 +13,9 @@ namespace palace
 {
 
 class ErrorIndicator;
-class IoData;
 class LumpedPortOperator;
 class Mesh;
 class PostOperator;
-class Timer;
 
 //
 // Driver class for eigenmode simulations.
@@ -27,7 +25,8 @@ class EigenSolver : public BaseSolver
 private:
   void Postprocess(const PostOperator &postop, const LumpedPortOperator &lumped_port_op,
                    int i, std::complex<double> omega, double error_bkwd, double error_abs,
-                   int num_conv, const ErrorIndicator *indicator) const;
+                   int num_conv, double E_elec, double E_mag,
+                   const ErrorIndicator *indicator) const;
 
   void PostprocessEigen(int i, std::complex<double> omega, double error_bkwd,
                         double error_abs, int num_conv) const;