diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp
index c6246e050..08bf90c28 100644
--- a/palace/drivers/basesolver.cpp
+++ b/palace/drivers/basesolver.cpp
@@ -287,41 +287,9 @@ void BaseSolver::SaveMetadata(const Timer &timer) const
   }
 }
 
-namespace
-{
-
-struct EnergyData
-{
-  const int idx;        // Domain index
-  const double E_elec;  // Electric field energy
-  const double E_mag;   // Magnetic field energy
-};
-
-struct FluxData
-{
-  const int idx;                   // Surface index
-  const std::complex<double> Phi;  // Integrated flux
-  const SurfaceFluxType type;      // Flux type
-};
-
-struct EpsData
-{
-  const int idx;   // Interface index
-  const double p;  // Participation ratio
-  const double Q;  // Quality factor
-};
-
-struct ProbeData
-{
-  const int idx;                          // Probe index
-  const std::complex<double> Fx, Fy, Fz;  // Field values at probe location
-};
-
-}  // namespace
-
 BaseSolver::DomainsPostPrinter::DomainsPostPrinter(bool do_measurement, bool root,
                                                    const fs::path &post_dir,
-                                                   const DomainPostOperator &dom_post_op,
+                                                   const PostOperator &post_op,
                                                    const std::string &idx_col_name,
                                                    int n_expected_rows)
   : do_measurement_{do_measurement}, root_{root}
@@ -333,7 +301,7 @@ BaseSolver::DomainsPostPrinter::DomainsPostPrinter(bool do_measurement, bool roo
   using fmt::format;
 
   domain_E = TableWithCSVFile(post_dir / "domain-E.csv");
-  domain_E.table.reserve(n_expected_rows, 4 + dom_post_op.M_i.size());
+  domain_E.table.reserve(n_expected_rows, 4 + post_op.GetDomainPostOp().M_i.size());
   domain_E.table.insert_column(Column("idx", idx_col_name, 0, {}, {}, ""));
 
   domain_E.table.insert_column("Ee", "E_elec (J)");
@@ -341,7 +309,7 @@ BaseSolver::DomainsPostPrinter::DomainsPostPrinter(bool do_measurement, bool roo
   domain_E.table.insert_column("Ec", "E_cap (J)");
   domain_E.table.insert_column("Ei", "E_ind (J)");
 
-  for (const auto &[idx, data] : dom_post_op.M_i)
+  for (const auto &[idx, data] : post_op.GetDomainPostOp().M_i)
   {
     domain_E.table.insert_column(format("Ee{}", idx), format("E_elec[{}] (J)", idx));
     domain_E.table.insert_column(format("pe{}", idx), format("p_elec[{}]", idx));
@@ -353,47 +321,36 @@ BaseSolver::DomainsPostPrinter::DomainsPostPrinter(bool do_measurement, bool roo
 
 void BaseSolver::DomainsPostPrinter::AddMeasurement(double idx_value_dimensionful,
                                                     const PostOperator &post_op,
-                                                    double E_elec, double E_mag,
-                                                    double E_cap, double E_ind,
                                                     const IoData &iodata)
 {
-  if (!do_measurement_)
-  {
-    return;
-  }
-
-  // MPI Gather
-  std::vector<EnergyData> energy_data;
-  energy_data.reserve(post_op.GetDomainPostOp().M_i.size());
-  for (const auto &[idx, data] : post_op.GetDomainPostOp().M_i)
-  {
-    double E_elec_i = (E_elec > 0.0) ? post_op.GetEFieldEnergy(idx) : 0.0;
-    double E_mag_i = (E_mag > 0.0) ? post_op.GetHFieldEnergy(idx) : 0.0;
-    energy_data.emplace_back(EnergyData{idx, E_elec_i, E_mag_i});
-  }
-
-  if (!root_)
+  if (!do_measurement_ || !root_)
   {
     return;
   }
-
   using VT = IoData::ValueType;
   using fmt::format;
 
+  double oneJ = iodata.DimensionalizeValue(VT::ENERGY, 1.0);
+
   domain_E.table["idx"] << idx_value_dimensionful;
 
-  domain_E.table["Ee"] << iodata.DimensionalizeValue(VT::ENERGY, E_elec);
-  domain_E.table["Em"] << iodata.DimensionalizeValue(VT::ENERGY, E_mag);
-  domain_E.table["Ec"] << iodata.DimensionalizeValue(VT::ENERGY, E_cap);
-  domain_E.table["Ei"] << iodata.DimensionalizeValue(VT::ENERGY, E_ind);
+  double E_elec = post_op.GetEFieldEnergy();
+  double E_mag = post_op.GetHFieldEnergy();
+
+  domain_E.table["Ee"] << E_elec * oneJ;
+  domain_E.table["Em"] << E_mag * oneJ;
+  domain_E.table["Ec"] << post_op.GetLumpedCapacitorEnergy() * oneJ;
+  domain_E.table["Ei"] << post_op.GetLumpedInductorEnergy() * oneJ;
 
   // Write the field and lumped element energies.
-  for (const auto &[idx, E_e, E_m] : energy_data)
+  for (const auto &[idx, data] : post_op.GetDomainPostOp().M_i)
   {
-    domain_E.table[format("Ee{}", idx)] << iodata.DimensionalizeValue(VT::ENERGY, E_e);
+    double E_e = post_op.GetEFieldEnergy(idx);
+    double E_m = post_op.GetHFieldEnergy(idx);
+    domain_E.table[format("Ee{}", idx)] << E_e * oneJ;
     domain_E.table[format("pe{}", idx)]
         << ((std::abs(E_elec) > 0.0) ? (E_e / E_elec) : 0.0);
-    domain_E.table[format("Em{}", idx)] << iodata.DimensionalizeValue(VT::ENERGY, E_m);
+    domain_E.table[format("Em{}", idx)] << E_m * oneJ;
     domain_E.table[format("pm{}", idx)] << ((std::abs(E_mag) > 0.0) ? E_m / E_mag : 0.0);
   }
 
@@ -485,24 +442,20 @@ BaseSolver::SurfacesPostPrinter::SurfacesPostPrinter(bool do_measurement, bool r
 
 void BaseSolver::SurfacesPostPrinter::AddMeasurementFlux(double idx_value_dimensionful,
                                                          const PostOperator &post_op,
-                                                         double E_elec, double E_mag,
                                                          const IoData &iodata)
 {
-  if (!do_measurement_flux_)
+  if (!do_measurement_flux_ || !root_)
   {
     return;
   }
   using VT = IoData::ValueType;
   using fmt::format;
 
-  // Gather
   const bool has_imaginary = post_op.HasImag();
-  std::vector<FluxData> flux_data;
-  flux_data.reserve(post_op.GetSurfacePostOp().flux_surfs.size());
-  for (const auto &[idx, data] : post_op.GetSurfacePostOp().flux_surfs)
+  auto flux_data_vec = post_op.GetSurfaceFluxAll();
+  auto dimensionlize_flux = [&iodata](auto Phi, SurfaceFluxType flux_type)
   {
-    auto Phi = post_op.GetSurfaceFlux(idx);
-    switch (data.type)
+    switch (flux_type)
     {
       case SurfaceFluxType::ELECTRIC:
         Phi *= iodata.DimensionalizeValue(VT::CAPACITANCE, 1.0);
@@ -516,23 +469,17 @@ void BaseSolver::SurfacesPostPrinter::AddMeasurementFlux(double idx_value_dimens
         Phi *= iodata.DimensionalizeValue(VT::POWER, 1.0);
         break;
     }
-    flux_data.emplace_back(FluxData{idx, Phi, data.type});
-  }
-
-  if (!root_)
-  {
-    return;
-  }
-
+    return Phi;
+  };
   surface_F.table["idx"] << idx_value_dimensionful;
-
-  for (const auto &[idx, Phi, data_type] : flux_data)
+  for (const auto &flux_data : flux_data_vec)
   {
-    surface_F.table[format("F_{}_re", idx)] << Phi.real();
-    if (has_imaginary &&
-        (data_type == SurfaceFluxType::ELECTRIC || data_type == SurfaceFluxType::MAGNETIC))
+    auto Phi_unitful = dimensionlize_flux(flux_data.Phi, flux_data.type);
+    surface_F.table[format("F_{}_re", flux_data.idx)] << Phi_unitful.real();
+    if (has_imaginary && (flux_data.type == SurfaceFluxType::ELECTRIC ||
+                          flux_data.type == SurfaceFluxType::MAGNETIC))
     {
-      surface_F.table[format("F_{}_im", idx)] << Phi.imag();
+      surface_F.table[format("F_{}_im", flux_data.idx)] << Phi_unitful.imag();
     }
   }
   surface_F.WriteFullTableTrunc();
@@ -540,51 +487,41 @@ void BaseSolver::SurfacesPostPrinter::AddMeasurementFlux(double idx_value_dimens
 
 void BaseSolver::SurfacesPostPrinter::AddMeasurementEps(double idx_value_dimensionful,
                                                         const PostOperator &post_op,
-                                                        double E_elec, double E_mag,
                                                         const IoData &iodata)
 {
-  if (!do_measurement_eps_)
+  if (!do_measurement_eps_ || !root_)
   {
     return;
   }
   using VT = IoData::ValueType;
   using fmt::format;
 
-  // Gather
-  std::vector<EpsData> eps_data;
-  eps_data.reserve(post_op.GetSurfacePostOp().eps_surfs.size());
-  for (const auto &[idx, data] : post_op.GetSurfacePostOp().eps_surfs)
-  {
-    double p = post_op.GetInterfaceParticipation(idx, E_elec);
-    double tandelta = post_op.GetSurfacePostOp().GetInterfaceLossTangent(idx);
-    double Q = (p == 0.0 || tandelta == 0.0) ? mfem::infinity() : 1.0 / (tandelta * p);
-    eps_data.emplace_back(EpsData{idx, p, Q});
-  }
-
-  if (!root_)
-  {
-    return;
-  }
+  // Interface Participation adds energy contriutions E_elec + E_cap
+  // E_cap returns zero if the solver does not supprot lumped ports.
+  double E_elec = post_op.GetEFieldEnergy() + post_op.GetLumpedCapacitorEnergy();
+  auto eps_data_vec = post_op.GetInterfaceEFieldEnergyAll();
 
   surface_Q.table["idx"] << idx_value_dimensionful;
-  for (const auto &[idx, p, Q] : eps_data)
+  for (const auto &eps_data : eps_data_vec)
   {
-    surface_Q.table[format("p_{}", idx)] << p;
-    surface_Q.table[format("Q_{}", idx)] << Q;
+    double p = post_op.GetInterfaceParticipation(eps_data.idx, E_elec);
+    double tandelta = eps_data.tandelta;
+    double Q = (p == 0.0 || tandelta == 0.0) ? mfem::infinity() : 1.0 / (tandelta * p);
+    surface_Q.table[format("p_{}", eps_data.idx)] << p;
+    surface_Q.table[format("Q_{}", eps_data.idx)] << Q;
   }
   surface_Q.WriteFullTableTrunc();
 }
 
 void BaseSolver::SurfacesPostPrinter::AddMeasurement(double idx_value_dimensionful,
                                                      const PostOperator &post_op,
-                                                     double E_elec, double E_mag,
                                                      const IoData &iodata)
 {
   // If surfaces have been specified for postprocessing, compute the corresponding values
   // and write out to disk. The passed in E_elec is the sum of the E-field and lumped
   // capacitor energies, and E_mag is the same for the B-field and lumped inductors.
-  AddMeasurementFlux(idx_value_dimensionful, post_op, E_elec, E_mag, iodata);
-  AddMeasurementEps(idx_value_dimensionful, post_op, E_elec, E_mag, iodata);
+  AddMeasurementFlux(idx_value_dimensionful, post_op, iodata);
+  AddMeasurementEps(idx_value_dimensionful, post_op, iodata);
 }
 
 BaseSolver::ProbePostPrinter::ProbePostPrinter(bool do_measurement, bool root,
@@ -690,7 +627,7 @@ void BaseSolver::ProbePostPrinter::AddMeasurementE(double idx_value_dimensionful
                                                    const PostOperator &post_op,
                                                    const IoData &iodata)
 {
-  if (!do_measurement_E_)
+  if (!do_measurement_E_ || !root_)
   {
     return;
   }
@@ -703,11 +640,6 @@ void BaseSolver::ProbePostPrinter::AddMeasurementE(double idx_value_dimensionful
                      v_dim, post_op.GetProbes().size(), v_dim * post_op.GetProbes().size(),
                      probe_field.size()))
 
-  if (!root_)
-  {
-    return;
-  }
-
   probe_E.table["idx"] << idx_value_dimensionful;
   size_t i = 0;
   for (const auto &idx : post_op.GetProbes())
@@ -730,7 +662,7 @@ void BaseSolver::ProbePostPrinter::AddMeasurementB(double idx_value_dimensionful
                                                    const PostOperator &post_op,
                                                    const IoData &iodata)
 {
-  if (!do_measurement_B_)
+  if (!do_measurement_B_ || !root_)
   {
     return;
   }
@@ -743,11 +675,6 @@ void BaseSolver::ProbePostPrinter::AddMeasurementB(double idx_value_dimensionful
                      v_dim, post_op.GetProbes().size(), v_dim * post_op.GetProbes().size(),
                      probe_field.size()))
 
-  if (!root_)
-  {
-    return;
-  }
-
   probe_B.table["idx"] << idx_value_dimensionful;
   size_t i = 0;
   for (const auto &idx : post_op.GetProbes())
@@ -785,7 +712,6 @@ BaseSolver::ErrorIndicatorPostPrinter::ErrorIndicatorPostPrinter(bool do_measure
   {
     return;
   }
-
   error_indicator = TableWithCSVFile(post_dir / "error-indicators.csv");
   error_indicator.table.reserve(1, 4);
 
@@ -796,27 +722,16 @@ BaseSolver::ErrorIndicatorPostPrinter::ErrorIndicatorPostPrinter(bool do_measure
 }
 
 void BaseSolver::ErrorIndicatorPostPrinter::PrintIndicatorStatistics(
-    const PostOperator &post_op, const ErrorIndicator &indicator)
+    const PostOperator &post_op, const ErrorIndicator::SummaryStatistics &indicator_stats)
 {
-  if (!do_measurement_)
-  {
-    return;
-  }
-
-  // MPI Gather
-  MPI_Comm comm = post_op.GetComm();
-  std::array<double, 4> data = {indicator.Norml2(comm), indicator.Min(comm),
-                                indicator.Max(comm), indicator.Mean(comm)};
-
-  if (!root_)
+  if (!do_measurement_ || !root_)
   {
     return;
   }
-
-  error_indicator.table["norm"] << data[0];
-  error_indicator.table["min"] << data[1];
-  error_indicator.table["max"] << data[2];
-  error_indicator.table["mean"] << data[3];
+  error_indicator.table["norm"] << indicator_stats.norm;
+  error_indicator.table["min"] << indicator_stats.min;
+  error_indicator.table["max"] << indicator_stats.max;
+  error_indicator.table["mean"] << indicator_stats.mean;
 
   error_indicator.WriteFullTableTrunc();
 }
diff --git a/palace/drivers/basesolver.hpp b/palace/drivers/basesolver.hpp
index f425789bb..d951cebbc 100644
--- a/palace/drivers/basesolver.hpp
+++ b/palace/drivers/basesolver.hpp
@@ -8,6 +8,7 @@
 #include <string>
 #include <vector>
 #include <fmt/os.h>
+#include "fem/errorindicator.hpp"
 #include "utils/filesystem.hpp"
 #include "utils/tablecsv.hpp"
 
@@ -15,7 +16,6 @@ namespace palace
 {
 
 class DomainPostOperator;
-class ErrorIndicator;
 class FiniteElementSpaceHierarchy;
 class IoData;
 class Mesh;
@@ -45,10 +45,9 @@ class BaseSolver
   public:
     DomainsPostPrinter() = default;
     DomainsPostPrinter(bool do_measurement, bool root, const fs::path &post_dir,
-                       const DomainPostOperator &dom_post_op,
-                       const std::string &idx_col_name, int n_expected_rows);
+                       const PostOperator &post_op, const std::string &idx_col_name,
+                       int n_expected_rows);
     void AddMeasurement(double idx_value_dimensionful, const PostOperator &post_op,
-                        double E_elec, double E_mag, double E_cap, double E_ind,
                         const IoData &iodata);
   };
 
@@ -67,11 +66,11 @@ class BaseSolver
                         const PostOperator &post_op, const std::string &idx_col_name,
                         int n_expected_rows);
     void AddMeasurement(double idx_value_dimensionful, const PostOperator &post_op,
-                        double E_elec, double E_mag, const IoData &iodata);
+                        const IoData &iodata);
     void AddMeasurementFlux(double idx_value_dimensionful, const PostOperator &post_op,
-                            double E_elec, double E_mag, const IoData &iodata);
+                            const IoData &iodata);
     void AddMeasurementEps(double idx_value_dimensionful, const PostOperator &post_op,
-                           double E_elec, double E_mag, const IoData &iodata);
+                           const IoData &iodata);
   };
 
   // Common probe postprocessing for all simulation types.
@@ -114,7 +113,7 @@ class BaseSolver
     ErrorIndicatorPostPrinter(bool do_measurement, bool root, const fs::path &post_dir);
 
     void PrintIndicatorStatistics(const PostOperator &post_op,
-                                  const ErrorIndicator &indicator);
+                                  const ErrorIndicator::SummaryStatistics &indicator_stats);
   };
 
   // Performs a solve using the mesh sequence, then reports error indicators and the number
diff --git a/palace/drivers/drivensolver.cpp b/palace/drivers/drivensolver.cpp
index 5d80de0f4..2ce50e37c 100644
--- a/palace/drivers/drivensolver.cpp
+++ b/palace/drivers/drivensolver.cpp
@@ -187,22 +187,25 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &space_op, PostOperator
     B *= -1.0 / (1i * omega);
     post_op.SetEGridFunction(E);
     post_op.SetBGridFunction(B);
-    post_op.UpdatePorts(space_op.GetLumpedPortOp(), space_op.GetWavePortOp(), omega);
-    const double E_elec = post_op.GetEFieldEnergy();
-    const double E_mag = post_op.GetHFieldEnergy();
+    post_op.SetFrequency(omega);
+    post_op.MeasureAll();
+
     Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(space_op.GetComm(), E),
                linalg::Norml2(space_op.GetComm(), RHS));
-    {
-      const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
-      Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", E_elec * J,
-                 E_mag * J, (E_elec + E_mag) * J);
-    }
+
+    const double E_elec = post_op.GetEFieldEnergy();
+    const double E_mag = post_op.GetHFieldEnergy();
+
+    const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
+    Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", E_elec * J,
+               E_mag * J, (E_elec + E_mag) * J);
+
     // Calculate and record the error indicators.
     Mpi::Print(" Updating solution error estimates\n");
     estimator.AddErrorIndicator(E, B, E_elec + E_mag, indicator);
 
-    post_results.PostprocessStep(iodata, post_op, space_op, step, omega, E_elec, E_mag);
+    post_results.PostprocessStep(iodata, post_op, space_op, step);
   }
   // Final postprocessing & printing
   BlockTimer bt0(Timer::POSTPRO);
@@ -363,16 +366,18 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &space_op, PostOperator
     B *= -1.0 / (1i * omega);
     post_op.SetEGridFunction(E);
     post_op.SetBGridFunction(B);
-    post_op.UpdatePorts(space_op.GetLumpedPortOp(), space_op.GetWavePortOp(), omega);
+    post_op.SetFrequency(omega);
+    post_op.MeasureAll();
+
+    Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(space_op.GetComm(), E));
+
     const double E_elec = post_op.GetEFieldEnergy();
     const double E_mag = post_op.GetHFieldEnergy();
-    Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(space_op.GetComm(), E));
-    {
-      const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
-      Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", E_elec * J,
-                 E_mag * J, (E_elec + E_mag) * J);
-    }
-    post_results.PostprocessStep(iodata, post_op, space_op, step, omega, E_elec, E_mag);
+    const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
+    Mpi::Print(" Field energy E ({:.3e} J) + H ({:.3e} J) = {:.3e} J\n", E_elec * J,
+               E_mag * J, (E_elec + E_mag) * J);
+
+    post_results.PostprocessStep(iodata, post_op, space_op, step);
   }
   // Final postprocessing & printing
   BlockTimer bt0(Timer::POSTPRO);
@@ -423,7 +428,7 @@ DrivenSolver::CurrentsPostPrinter::CurrentsPostPrinter(
 }
 
 void DrivenSolver::CurrentsPostPrinter::AddMeasurement(
-    double omega, const SurfaceCurrentOperator &surf_j_op, const IoData &iodata)
+    double freq, const SurfaceCurrentOperator &surf_j_op, const IoData &iodata)
 {
   if (!do_measurement_ || !root_)
   {
@@ -432,7 +437,7 @@ void DrivenSolver::CurrentsPostPrinter::AddMeasurement(
   using VT = IoData::ValueType;
   using fmt::format;
 
-  surface_I.table["idx"] << iodata.DimensionalizeValue(VT::FREQUENCY, omega);
+  surface_I.table["idx"] << freq;
   for (const auto &[idx, data] : surf_j_op)
   {
     auto I_inc = data.GetExcitationCurrent();
@@ -482,7 +487,7 @@ DrivenSolver::PortsPostPrinter::PortsPostPrinter(bool do_measurement, bool root,
 }
 
 void DrivenSolver::PortsPostPrinter::AddMeasurement(
-    double omega, const PostOperator &post_op, const LumpedPortOperator &lumped_port_op,
+    double freq, const PostOperator &post_op, const LumpedPortOperator &lumped_port_op,
     const IoData &iodata)
 {
   if (!do_measurement_ || !root_)
@@ -493,7 +498,6 @@ void DrivenSolver::PortsPostPrinter::AddMeasurement(
 
   // Postprocess the frequency domain lumped port voltages and currents (complex magnitude
   // = sqrt(2) * RMS).
-  auto freq = iodata.DimensionalizeValue(VT::FREQUENCY, omega);
   port_V.table["idx"] << freq;
   port_I.table["idx"] << freq;
 
@@ -511,8 +515,8 @@ void DrivenSolver::PortsPostPrinter::AddMeasurement(
       port_I.table[fmt::format("inc{}", idx)] << I_inc * unit_A;
     }
 
-    std::complex<double> V_i = post_op.GetPortVoltage(lumped_port_op, idx);
-    std::complex<double> I_i = post_op.GetPortCurrent(lumped_port_op, idx);
+    std::complex<double> V_i = post_op.GetPortVoltage(idx);
+    std::complex<double> I_i = post_op.GetPortCurrent(idx);
 
     port_V.table[fmt::format("re{}", idx)] << V_i.real() * unit_V;
     port_V.table[fmt::format("im{}", idx)] << V_i.imag() * unit_V;
@@ -541,7 +545,6 @@ DrivenSolver::SParametersPostPrinter::SParametersPostPrinter(
   {
     return;
   }
-
   // Get excitation index as is currently done: if -1 then no excitation
   // Already ensured that one of lumped or wave ports are empty
   for (const auto &[idx, data] : lumped_port_op)
@@ -565,7 +568,6 @@ DrivenSolver::SParametersPostPrinter::SParametersPostPrinter(
   {
     return;
   }
-
   using fmt::format;
   port_S = TableWithCSVFile(post_dir / "port-S.csv");
   port_S.table.reserve(n_expected_rows, lumped_port_op.Size());
@@ -590,7 +592,7 @@ DrivenSolver::SParametersPostPrinter::SParametersPostPrinter(
 }
 
 void DrivenSolver::SParametersPostPrinter::AddMeasurement(
-    double omega, const PostOperator &post_op, const LumpedPortOperator &lumped_port_op,
+    double freq, const PostOperator &post_op, const LumpedPortOperator &lumped_port_op,
     const WavePortOperator &wave_port_op, const IoData &iodata)
 {
   if (!do_measurement_ || !root_)
@@ -601,7 +603,7 @@ void DrivenSolver::SParametersPostPrinter::AddMeasurement(
   using fmt::format;
 
   // Add frequencies
-  port_S.table["idx"] << iodata.DimensionalizeValue(VT::FREQUENCY, omega);
+  port_S.table["idx"] << freq;
 
   std::vector<int> all_port_indices;
   for (const auto &[idx, data] : lumped_port_op)
@@ -615,15 +617,8 @@ void DrivenSolver::SParametersPostPrinter::AddMeasurement(
 
   for (const auto o_idx : all_port_indices)
   {
-    std::complex<double> S_ij;
-    if (src_lumped_port)
-    {
-      S_ij = post_op.GetSParameter(lumped_port_op, o_idx, source_idx);
-    }
-    else
-    {
-      S_ij = post_op.GetSParameter(wave_port_op, o_idx, source_idx);
-    }
+    std::complex<double> S_ij = post_op.GetSParameter(src_lumped_port, o_idx, source_idx);
+
     auto abs_S_ij = 20.0 * std::log10(std::abs(S_ij));
     auto arg_S_ij = std::arg(S_ij) * 180.8 / M_PI;
 
@@ -642,7 +637,7 @@ DrivenSolver::PostprocessPrintResults::PostprocessPrintResults(
     bool root, const fs::path &post_dir, const PostOperator &post_op,
     const SpaceOperator &space_op, int n_expected_rows, int delta_post_)
   : delta_post{delta_post_}, write_paraview_fields{delta_post_ > 0},
-    domains{true, root, post_dir, post_op.GetDomainPostOp(), "f (GHz)", n_expected_rows},
+    domains{true, root, post_dir, post_op, "f (GHz)", n_expected_rows},
     surfaces{true, root, post_dir, post_op, "f (GHz)", n_expected_rows},
     currents{true, root, post_dir, space_op.GetSurfaceCurrentOp(), n_expected_rows},
     probes{true, root, post_dir, post_op, "f (GHz)", n_expected_rows},
@@ -660,19 +655,17 @@ DrivenSolver::PostprocessPrintResults::PostprocessPrintResults(
 void DrivenSolver::PostprocessPrintResults::PostprocessStep(const IoData &iodata,
                                                             const PostOperator &post_op,
                                                             const SpaceOperator &space_op,
-                                                            int step, double omega,
-                                                            double E_elec, double E_mag)
+                                                            int step)
 {
+  double omega = post_op.GetFrequency().real();
   auto freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
-  auto E_cap = post_op.GetLumpedCapacitorEnergy(space_op.GetLumpedPortOp());
-  auto E_ind = post_op.GetLumpedInductorEnergy(space_op.GetLumpedPortOp());
 
-  domains.AddMeasurement(freq, post_op, E_elec, E_mag, E_cap, E_ind, iodata);
-  surfaces.AddMeasurement(freq, post_op, E_elec + E_cap, E_mag + E_ind, iodata);
-  currents.AddMeasurement(omega, space_op.GetSurfaceCurrentOp(), iodata);
+  domains.AddMeasurement(freq, post_op, iodata);
+  surfaces.AddMeasurement(freq, post_op, iodata);
+  currents.AddMeasurement(freq, space_op.GetSurfaceCurrentOp(), iodata);
   probes.AddMeasurement(freq, post_op, iodata);
-  ports.AddMeasurement(omega, post_op, space_op.GetLumpedPortOp(), iodata);
-  s_parameters.AddMeasurement(omega, post_op, space_op.GetLumpedPortOp(),
+  ports.AddMeasurement(freq, post_op, space_op.GetLumpedPortOp(), iodata);
+  s_parameters.AddMeasurement(freq, post_op, space_op.GetLumpedPortOp(),
                               space_op.GetWavePortOp(), iodata);
   // The internal GridFunctions in PostOperator have already been set:
   if (write_paraview_fields && (step % delta_post == 0))
@@ -687,7 +680,8 @@ void DrivenSolver::PostprocessPrintResults::PostprocessFinal(
     const PostOperator &post_op, const ErrorIndicator &indicator)
 {
   BlockTimer bt0(Timer::POSTPRO);
-  error_indicator.PrintIndicatorStatistics(post_op, indicator);
+  auto indicator_stats = indicator.GetSummaryStatistics(post_op.GetComm());
+  error_indicator.PrintIndicatorStatistics(post_op, indicator_stats);
   if (write_paraview_fields)
   {
     post_op.WriteFieldsFinal(&indicator);
diff --git a/palace/drivers/drivensolver.hpp b/palace/drivers/drivensolver.hpp
index cc1eb8208..1c8bd51f8 100644
--- a/palace/drivers/drivensolver.hpp
+++ b/palace/drivers/drivensolver.hpp
@@ -40,7 +40,7 @@ class DrivenSolver : public BaseSolver
     CurrentsPostPrinter() = default;
     CurrentsPostPrinter(bool do_measurement, bool root, const fs::path &post_dir,
                         const SurfaceCurrentOperator &surf_j_op, int n_expected_rows);
-    void AddMeasurement(double omega, const SurfaceCurrentOperator &surf_j_op,
+    void AddMeasurement(double freq, const SurfaceCurrentOperator &surf_j_op,
                         const IoData &iodata);
   };
 
@@ -55,7 +55,7 @@ class DrivenSolver : public BaseSolver
     PortsPostPrinter() = default;
     PortsPostPrinter(bool do_measurement, bool root, const fs::path &post_dir,
                      const LumpedPortOperator &lumped_port_op, int n_expected_rows);
-    void AddMeasurement(double omega, const PostOperator &post_op,
+    void AddMeasurement(double freq, const PostOperator &post_op,
                         const LumpedPortOperator &lumped_port_op, const IoData &iodata);
   };
 
@@ -79,7 +79,7 @@ class DrivenSolver : public BaseSolver
     SParametersPostPrinter(bool do_measurement, bool root, const fs::path &post_dir,
                            const LumpedPortOperator &lumped_port_op,
                            const WavePortOperator &wave_port_op, int n_expected_rows);
-    void AddMeasurement(double omega, const PostOperator &post_op,
+    void AddMeasurement(double freq, const PostOperator &post_op,
                         const LumpedPortOperator &lumped_port_op,
                         const WavePortOperator &wave_port_op, const IoData &iodata);
   };
@@ -102,8 +102,7 @@ class DrivenSolver : public BaseSolver
                             const PostOperator &post_op, const SpaceOperator &space_op,
                             int n_expected_rows, int delta_post);
     void PostprocessStep(const IoData &iodata, const PostOperator &post_op,
-                         const SpaceOperator &space_op, int step, double omega,
-                         double E_elec, double E_mag);
+                         const SpaceOperator &space_op, int step);
     void PostprocessFinal(const PostOperator &post_op, const ErrorIndicator &indicator);
   };
 
diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp
index 2d0e9e3bf..1304a031e 100644
--- a/palace/drivers/eigensolver.cpp
+++ b/palace/drivers/eigensolver.cpp
@@ -313,7 +313,9 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     B *= -1.0 / (1i * omega);
     post_op.SetEGridFunction(E);
     post_op.SetBGridFunction(B);
-    post_op.UpdatePorts(space_op.GetLumpedPortOp(), omega.real());
+    post_op.SetFrequency(omega);
+    post_op.MeasureAll();
+
     const double E_elec = post_op.GetEFieldEnergy();
     const double E_mag = post_op.GetHFieldEnergy();
 
@@ -324,8 +326,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     }
 
     // Postprocess state and write fields to file
-    post_results.PostprocessStep(iodata, post_op, space_op, i, omega, E_elec, E_mag,
-                                 error_abs, error_bkwd);
+    post_results.PostprocessStep(iodata, post_op, space_op, i, error_abs, error_bkwd);
 
     // Final write: Different condition than end of loop (i = num_conv - 1)
     if (i == iodata.solver.eigenmode.n - 1)
@@ -336,24 +337,6 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   return {indicator, space_op.GlobalTrueVSize()};
 }
 
-namespace
-{
-
-struct EprLData
-{
-  int idx;    // Lumped port index
-  double pj;  // Inductor energy-participation ratio
-};
-
-struct EprIOData
-{
-  int idx;    // Lumped port index
-  double Ql;  // Quality factor
-  double Kl;  // κ for loss rate
-};
-
-}  // namespace
-
 void EigenSolver::EigenPostPrinter::PrintStdoutHeader()
 {
   if (!root_)
@@ -417,9 +400,9 @@ void EigenSolver::EigenPostPrinter::PrintStdoutRow(size_t j)
 }
 
 EigenSolver::EigenPostPrinter::EigenPostPrinter(bool do_measurement, bool root,
-                                                const fs::path &post_dir, int n_eig)
+                                                const fs::path &post_dir, int n_post)
   : root_{root}, do_measurement_(do_measurement),
-    stdout_int_print_width(1 + static_cast<int>(std::log10(n_eig)))
+    stdout_int_print_width(1 + static_cast<int>(std::log10(n_post)))
 {
   // Note: we switch to n_eig rather than n_conv for padding since we don't know n_conv
   // until solve
@@ -428,7 +411,7 @@ EigenSolver::EigenPostPrinter::EigenPostPrinter(bool do_measurement, bool root,
     return;
   }
   eig = TableWithCSVFile(post_dir / "eig.csv");
-  eig.table.reserve(n_eig, 6);
+  eig.table.reserve(n_post, 6);
   eig.table.insert_column(Column("idx", "m", 0, {}, {}, ""));
   eig.table.insert_column("f_re", "Re{f} (GHz)");
   eig.table.insert_column("f_im", "Im{f} (GHz)");
@@ -439,7 +422,7 @@ EigenSolver::EigenPostPrinter::EigenPostPrinter(bool do_measurement, bool root,
 }
 
 void EigenSolver::EigenPostPrinter::AddMeasurement(int eigen_print_idx,
-                                                   std::complex<double> omega,
+                                                   const PostOperator &post_op,
                                                    double error_bkwd, double error_abs,
                                                    const IoData &iodata)
 {
@@ -450,7 +433,8 @@ void EigenSolver::EigenPostPrinter::AddMeasurement(int eigen_print_idx,
   using VT = IoData::ValueType;
   using fmt::format;
 
-  std::complex<double> f = iodata.DimensionalizeValue(VT::FREQUENCY, omega);
+  std::complex<double> f =
+      iodata.DimensionalizeValue(VT::FREQUENCY, post_op.GetFrequency());
   double Q = (f.imag() == 0.0) ? mfem::infinity() : 0.5 * std::abs(f) / std::abs(f.imag());
 
   eig.table["idx"] << eigen_print_idx;
@@ -519,8 +503,8 @@ void EigenSolver::PortsPostPrinter::AddMeasurement(int eigen_print_idx,
 
   for (const auto &[idx, data] : lumped_port_op)
   {
-    std::complex<double> V_i = post_op.GetPortVoltage(lumped_port_op, idx);
-    std::complex<double> I_i = post_op.GetPortCurrent(lumped_port_op, idx);
+    std::complex<double> V_i = post_op.GetPortVoltage(idx);
+    std::complex<double> I_i = post_op.GetPortCurrent(idx);
 
     port_V.table[fmt::format("re{}", idx)] << V_i.real() * unit_V;
     port_V.table[fmt::format("im{}", idx)] << V_i.imag() * unit_V;
@@ -590,33 +574,22 @@ EigenSolver::EPRPostPrinter::EPRPostPrinter(bool do_measurement, bool root,
 
 void EigenSolver::EPRPostPrinter::AddMeasurementEPR(
     double eigen_print_idx, const PostOperator &post_op,
-    const LumpedPortOperator &lumped_port_op, double E_m, const IoData &iodata)
+    const LumpedPortOperator &lumped_port_op, const IoData &iodata)
 {
-  if (!do_measurement_EPR_)
-  {
-    return;
-  }
-
-  // TODO: Does this include a reduce?
-  // Write the mode EPR for lumped inductor elements.
-  std::vector<EprLData> epr_L_data;
-  epr_L_data.reserve(ports_with_L.size());
-  for (const auto idx : ports_with_L)
-  {
-    const double pj = post_op.GetInductorParticipation(lumped_port_op, idx, E_m);
-    epr_L_data.push_back({idx, pj});
-  }
-
-  if (!root_)
+  if (!do_measurement_EPR_ || !root_)
   {
     return;
   }
   using fmt::format;
 
+  double E_elec = post_op.GetEFieldEnergy();
+  double E_cap = post_op.GetLumpedCapacitorEnergy();
+  double E_m = E_elec + E_cap;
+
   port_EPR.table["idx"] << eigen_print_idx;
-  for (const auto &[idx, pj] : epr_L_data)
+  for (const auto idx : ports_with_L)
   {
-    port_EPR.table[format("p_{}", idx)] << pj;
+    port_EPR.table[format("p_{}", idx)] << post_op.GetInductorParticipation(idx, E_m);
   }
   port_EPR.AppendRow();
 }
@@ -624,37 +597,28 @@ void EigenSolver::EPRPostPrinter::AddMeasurementEPR(
 void EigenSolver::EPRPostPrinter::AddMeasurementQ(double eigen_print_idx,
                                                   const PostOperator &post_op,
                                                   const LumpedPortOperator &lumped_port_op,
-                                                  std::complex<double> omega, double E_m,
                                                   const IoData &iodata)
 {
-  if (!do_measurement_Q_)
-  {
-    return;
-  }
-
-  // TODO: Does this include a reduce?
-  // Write the mode EPR for lumped resistor elements.
-  std::vector<EprIOData> epr_IO_data;
-  epr_IO_data.reserve(ports_with_R.size());
-  for (const auto idx : ports_with_R)
-  {
-    const double Kl = post_op.GetExternalKappa(lumped_port_op, idx, E_m);
-    const double Ql = (Kl == 0.0) ? mfem::infinity() : omega.real() / std::abs(Kl);
-    epr_IO_data.push_back(
-        {idx, Ql, iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, Kl)});
-  }
-
-  if (!root_)
+  if (!do_measurement_Q_ || !root_)
   {
     return;
   }
   using fmt::format;
+  using VT = IoData::ValueType;
+
+  auto omega = post_op.GetFrequency();
+  double E_elec = post_op.GetEFieldEnergy();
+  double E_cap = post_op.GetLumpedCapacitorEnergy();
+  double E_m = E_elec + E_cap;
 
   port_EPR.table["idx"] << eigen_print_idx;
-  for (const auto &[idx, Ql, Kl] : epr_IO_data)
+  for (const auto idx : ports_with_R)
   {
+    double Kl = post_op.GetExternalKappa(idx, E_m);
+    double Ql = (Kl == 0.0) ? mfem::infinity() : omega.real() / std::abs(Kl);
+
     port_Q.table[format("Ql_{}", idx)] << Ql;
-    port_Q.table[format("Kl_{}", idx)] << Kl;
+    port_Q.table[format("Kl_{}", idx)] << iodata.DimensionalizeValue(VT::FREQUENCY, Kl);
   }
   port_Q.AppendRow();
 }
@@ -662,11 +626,10 @@ void EigenSolver::EPRPostPrinter::AddMeasurementQ(double eigen_print_idx,
 void EigenSolver::EPRPostPrinter::AddMeasurement(double eigen_print_idx,
                                                  const PostOperator &post_op,
                                                  const LumpedPortOperator &lumped_port_op,
-                                                 std::complex<double> omega, double E_m,
                                                  const IoData &iodata)
 {
-  AddMeasurementEPR(eigen_print_idx, post_op, lumped_port_op, E_m, iodata);
-  AddMeasurementQ(eigen_print_idx, post_op, lumped_port_op, omega, E_m, iodata);
+  AddMeasurementEPR(eigen_print_idx, post_op, lumped_port_op, iodata);
+  AddMeasurementQ(eigen_print_idx, post_op, lumped_port_op, iodata);
 }
 
 EigenSolver::PostprocessPrintResults::PostprocessPrintResults(bool root,
@@ -675,7 +638,7 @@ EigenSolver::PostprocessPrintResults::PostprocessPrintResults(bool root,
                                                               const SpaceOperator &space_op,
                                                               int n_post_)
   : n_post(n_post_), write_paraview_fields(n_post_ > 0),
-    domains{true, root, post_dir, post_op.GetDomainPostOp(), "m", n_post},
+    domains{true, root, post_dir, post_op, "m", n_post},
     surfaces{true, root, post_dir, post_op, "m", n_post},
     probes{true, root, post_dir, post_op, "m", n_post}, eigen{true, root, post_dir, n_post},
     epr{true, root, post_dir, space_op.GetLumpedPortOp(), n_post},
@@ -683,22 +646,19 @@ EigenSolver::PostprocessPrintResults::PostprocessPrintResults(bool root,
 {
 }
 
-void EigenSolver::PostprocessPrintResults::PostprocessStep(
-    const IoData &iodata, const PostOperator &post_op, const SpaceOperator &space_op,
-    int step, std::complex<double> omega, double E_elec, double E_mag, double error_abs,
-    double error_bkward)
+void EigenSolver::PostprocessPrintResults::PostprocessStep(const IoData &iodata,
+                                                           const PostOperator &post_op,
+                                                           const SpaceOperator &space_op,
+                                                           int step, double error_abs,
+                                                           double error_bkward)
 {
   int eigen_print_idx = step + 1;
 
-  auto E_cap = post_op.GetLumpedCapacitorEnergy(space_op.GetLumpedPortOp());
-  auto E_ind = post_op.GetLumpedInductorEnergy(space_op.GetLumpedPortOp());
-
-  domains.AddMeasurement(eigen_print_idx, post_op, E_elec, E_mag, E_cap, E_ind, iodata);
-  surfaces.AddMeasurement(eigen_print_idx, post_op, E_elec + E_cap, E_mag + E_ind, iodata);
+  domains.AddMeasurement(eigen_print_idx, post_op, iodata);
+  surfaces.AddMeasurement(eigen_print_idx, post_op, iodata);
   probes.AddMeasurement(eigen_print_idx, post_op, iodata);
-  eigen.AddMeasurement(eigen_print_idx, omega, error_bkward, error_abs, iodata);
-  epr.AddMeasurement(eigen_print_idx, post_op, space_op.GetLumpedPortOp(), omega,
-                     E_elec + E_cap, iodata);
+  eigen.AddMeasurement(eigen_print_idx, post_op, error_bkward, error_abs, iodata);
+  epr.AddMeasurement(eigen_print_idx, post_op, space_op.GetLumpedPortOp(), iodata);
   // The internal GridFunctions in PostOperator have already been set:
   if (write_paraview_fields && step < n_post)
   {
@@ -712,7 +672,8 @@ void EigenSolver::PostprocessPrintResults::PostprocessFinal(const PostOperator &
                                                             const ErrorIndicator &indicator)
 {
   BlockTimer bt0(Timer::POSTPRO);
-  error_indicator.PrintIndicatorStatistics(post_op, indicator);
+  auto indicator_stats = indicator.GetSummaryStatistics(post_op.GetComm());
+  error_indicator.PrintIndicatorStatistics(post_op, indicator_stats);
   if (write_paraview_fields)
   {
     post_op.WriteFieldsFinal(&indicator);
diff --git a/palace/drivers/eigensolver.hpp b/palace/drivers/eigensolver.hpp
index 4fe4623c6..66ab809bf 100644
--- a/palace/drivers/eigensolver.hpp
+++ b/palace/drivers/eigensolver.hpp
@@ -39,7 +39,7 @@ class EigenSolver : public BaseSolver
 
     EigenPostPrinter() = default;
     EigenPostPrinter(bool do_measurement, bool root, const fs::path &post_dir, int n_post);
-    void AddMeasurement(int eigen_print_idx, std::complex<double> omega, double error_bkwd,
+    void AddMeasurement(int eigen_print_idx, const PostOperator &post_op, double error_bkwd,
                         double error_abs, const IoData &iodata);
   };
 
@@ -76,15 +76,12 @@ class EigenSolver : public BaseSolver
                    const LumpedPortOperator &lumped_port_op, int n_expected_rows);
 
     void AddMeasurementEPR(double eigen_print_idx, const PostOperator &post_op,
-                           const LumpedPortOperator &lumped_port_op, double E_m,
-                           const IoData &iodata);
+                           const LumpedPortOperator &lumped_port_op, const IoData &iodata);
     void AddMeasurementQ(double eigen_print_idx, const PostOperator &post_op,
-                         const LumpedPortOperator &lumped_port_op,
-                         std::complex<double> omega, double E_m, const IoData &iodata);
+                         const LumpedPortOperator &lumped_port_op, const IoData &iodata);
 
     void AddMeasurement(double eigen_print_idx, const PostOperator &post_op,
-                        const LumpedPortOperator &lumped_port_op,
-                        std::complex<double> omega, double E_m, const IoData &iodata);
+                        const LumpedPortOperator &lumped_port_op, const IoData &iodata);
   };
 
   struct PostprocessPrintResults
@@ -104,9 +101,8 @@ class EigenSolver : public BaseSolver
                             const PostOperator &post_op, const SpaceOperator &space_op,
                             int n_post_);
     void PostprocessStep(const IoData &iodata, const PostOperator &post_op,
-                         const SpaceOperator &space_op, int step,
-                         std::complex<double> omega, double E_elec, double E_mag,
-                         double error_abs, double error_bkward);
+                         const SpaceOperator &space_op, int step, double error_abs,
+                         double error_bkward);
     void PostprocessFinal(const PostOperator &post_op, const ErrorIndicator &indicator);
   };
 
diff --git a/palace/drivers/electrostaticsolver.cpp b/palace/drivers/electrostaticsolver.cpp
index 412e97287..5efa5af9a 100644
--- a/palace/drivers/electrostaticsolver.cpp
+++ b/palace/drivers/electrostaticsolver.cpp
@@ -77,6 +77,7 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Grad.AddMult(V[step], E, -1.0);
     post_op.SetVGridFunction(V[step]);
     post_op.SetEGridFunction(E);
+    post_op.MeasureAll();
     const double E_elec = post_op.GetEFieldEnergy();
     Mpi::Print(" Sol. ||V|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(laplace_op.GetComm(), V[step]),
@@ -91,7 +92,7 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     estimator.AddErrorIndicator(E, E_elec, indicator);
 
     // Postprocess field solutions and optionally write solution to disk.
-    post_results.PostprocessStep(iodata, post_op, step, idx, E_elec);
+    post_results.PostprocessStep(iodata, post_op, step, idx);
 
     // Next terminal.
     step++;
@@ -203,7 +204,7 @@ void ElectrostaticSolver::PostprocessTerminals(
 ElectrostaticSolver::PostprocessPrintResults::PostprocessPrintResults(
     bool root, const fs::path &post_dir, const PostOperator &post_op, int n_post_)
   : n_post(n_post_), write_paraview_fields(n_post_ > 0),
-    domains{true, root, post_dir, post_op.GetDomainPostOp(), "i", n_post},
+    domains{true, root, post_dir, post_op, "i", n_post},
     surfaces{true, root, post_dir, post_op, "i", n_post},
     probes{true, root, post_dir, post_op, "i", n_post},
     error_indicator{true, root, post_dir}
@@ -211,10 +212,10 @@ ElectrostaticSolver::PostprocessPrintResults::PostprocessPrintResults(
 }
 
 void ElectrostaticSolver::PostprocessPrintResults::PostprocessStep(
-    const IoData &iodata, const PostOperator &post_op, int step, int idx, double E_elec)
+    const IoData &iodata, const PostOperator &post_op, int step, int idx)
 {
-  domains.AddMeasurement(idx, post_op, E_elec, 0.0, 0.0, 0.0, iodata);
-  surfaces.AddMeasurement(idx, post_op, E_elec, 0.0, iodata);
+  domains.AddMeasurement(idx, post_op, iodata);
+  surfaces.AddMeasurement(idx, post_op, iodata);
   probes.AddMeasurement(idx, post_op, iodata);
   // The internal GridFunctions in PostOperator have already been set from V:
   if (write_paraview_fields && step < n_post)
@@ -229,7 +230,8 @@ void ElectrostaticSolver::PostprocessPrintResults::PostprocessFinal(
     const PostOperator &post_op, const ErrorIndicator &indicator)
 {
   BlockTimer bt0(Timer::POSTPRO);
-  error_indicator.PrintIndicatorStatistics(post_op, indicator);
+  auto indicator_stats = indicator.GetSummaryStatistics(post_op.GetComm());
+  error_indicator.PrintIndicatorStatistics(post_op, indicator_stats);
   if (write_paraview_fields)
   {
     post_op.WriteFieldsFinal(&indicator);
diff --git a/palace/drivers/electrostaticsolver.hpp b/palace/drivers/electrostaticsolver.hpp
index 0a71e8fbd..283d6c540 100644
--- a/palace/drivers/electrostaticsolver.hpp
+++ b/palace/drivers/electrostaticsolver.hpp
@@ -45,7 +45,7 @@ class ElectrostaticSolver : public BaseSolver
     PostprocessPrintResults(bool is_mpi_root, const fs::path &post_dir,
                             const PostOperator &post_op, int n_post_);
     void PostprocessStep(const IoData &iodata, const PostOperator &post_op, int step,
-                         int idx, double E_elec);
+                         int idx);
     void PostprocessFinal(const PostOperator &post_op, const ErrorIndicator &indicator);
   };
 
diff --git a/palace/drivers/magnetostaticsolver.cpp b/palace/drivers/magnetostaticsolver.cpp
index d5fac30f9..f4321daf2 100644
--- a/palace/drivers/magnetostaticsolver.cpp
+++ b/palace/drivers/magnetostaticsolver.cpp
@@ -79,6 +79,7 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Curl.Mult(A[step], B);
     post_op.SetAGridFunction(A[step]);
     post_op.SetBGridFunction(B);
+    post_op.MeasureAll();
     const double E_mag = post_op.GetHFieldEnergy();
     Mpi::Print(" Sol. ||A|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(curlcurl_op.GetComm(), A[step]),
@@ -94,7 +95,7 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     estimator.AddErrorIndicator(B, E_mag, indicator);
 
     // Postprocess field solutions and optionally write solution to disk.
-    post_results.PostprocessStep(iodata, post_op, step, idx, E_mag);
+    post_results.PostprocessStep(iodata, post_op, step, idx);
 
     // Next source.
     step++;
@@ -208,7 +209,7 @@ void MagnetostaticSolver::PostprocessTerminals(PostOperator &post_op,
 MagnetostaticSolver::PostprocessPrintResults::PostprocessPrintResults(
     bool root, const fs::path &post_dir, const PostOperator &post_op, int n_post_)
   : n_post(n_post_), write_paraview_fields(n_post_ > 0),
-    domains{true, root, post_dir, post_op.GetDomainPostOp(), "i", n_post},
+    domains{true, root, post_dir, post_op, "i", n_post},
     surfaces{true, root, post_dir, post_op, "i", n_post},
     probes{true, root, post_dir, post_op, "i", n_post},
     error_indicator{true, root, post_dir}
@@ -216,10 +217,10 @@ MagnetostaticSolver::PostprocessPrintResults::PostprocessPrintResults(
 }
 
 void MagnetostaticSolver::PostprocessPrintResults::PostprocessStep(
-    const IoData &iodata, const PostOperator &post_op, int step, int idx, double E_mag)
+    const IoData &iodata, const PostOperator &post_op, int step, int idx)
 {
-  domains.AddMeasurement(idx, post_op, 0.0, E_mag, 0.0, 0.0, iodata);
-  surfaces.AddMeasurement(idx, post_op, 0.0, E_mag, iodata);
+  domains.AddMeasurement(idx, post_op, iodata);
+  surfaces.AddMeasurement(idx, post_op, iodata);
   probes.AddMeasurement(idx, post_op, iodata);
   // The internal GridFunctions in PostOperator have already been set from A:
   if (write_paraview_fields && step < n_post)
@@ -234,7 +235,8 @@ void MagnetostaticSolver::PostprocessPrintResults::PostprocessFinal(
     const PostOperator &post_op, const ErrorIndicator &indicator)
 {
   BlockTimer bt0(Timer::POSTPRO);
-  error_indicator.PrintIndicatorStatistics(post_op, indicator);
+  auto indicator_stats = indicator.GetSummaryStatistics(post_op.GetComm());
+  error_indicator.PrintIndicatorStatistics(post_op, indicator_stats);
   if (write_paraview_fields)
   {
     post_op.WriteFieldsFinal(&indicator);
diff --git a/palace/drivers/magnetostaticsolver.hpp b/palace/drivers/magnetostaticsolver.hpp
index a95cf91ef..608526c24 100644
--- a/palace/drivers/magnetostaticsolver.hpp
+++ b/palace/drivers/magnetostaticsolver.hpp
@@ -37,7 +37,7 @@ class MagnetostaticSolver : public BaseSolver
     PostprocessPrintResults(bool is_mpi_root, const fs::path &post_dir,
                             const PostOperator &post_op, int n_post_);
     void PostprocessStep(const IoData &iodata, const PostOperator &post_op, int step,
-                         int idx, double E_mag);
+                         int idx);
     void PostprocessFinal(const PostOperator &post_op, const ErrorIndicator &indicator);
   };
 
diff --git a/palace/drivers/transientsolver.cpp b/palace/drivers/transientsolver.cpp
index c5d47ca85..88bd7f8dd 100644
--- a/palace/drivers/transientsolver.cpp
+++ b/palace/drivers/transientsolver.cpp
@@ -116,7 +116,8 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     const Vector &B = time_op.GetB();
     post_op.SetEGridFunction(E);
     post_op.SetBGridFunction(B);
-    post_op.UpdatePorts(space_op.GetLumpedPortOp());
+    post_op.MeasureAll();
+
     const double E_elec = post_op.GetEFieldEnergy();
     const double E_mag = post_op.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e}, ||B|| = {:.6e}\n",
@@ -131,8 +132,7 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Mpi::Print(" Updating solution error estimates\n");
     estimator.AddErrorIndicator(E, B, E_elec + E_mag, indicator);
 
-    post_results.PostprocessStep(iodata, post_op, space_op, step, t, J_coef(t), E_elec,
-                                 E_mag);
+    post_results.PostprocessStep(iodata, post_op, space_op, step, t, J_coef(t));
   }
   // Final postprocessing & printing
   BlockTimer bt1(Timer::POSTPRO);
@@ -364,8 +364,8 @@ void TransientSolver::PortsPostPrinter::AddMeasurement(
       port_I.table[format("inc{}", idx)] << I_inc * unit_A;
     }
 
-    std::complex<double> V_i = post_op.GetPortVoltage(lumped_port_op, idx);
-    std::complex<double> I_i = post_op.GetPortCurrent(lumped_port_op, idx);
+    std::complex<double> V_i = post_op.GetPortVoltage(idx);
+    std::complex<double> I_i = post_op.GetPortCurrent(idx);
 
     port_V.table[format("re{}", idx)] << V_i.real() * unit_V;
     port_I.table[format("re{}", idx)] << I_i.real() * unit_A;
@@ -378,7 +378,7 @@ TransientSolver::PostprocessPrintResults::PostprocessPrintResults(
     bool root, const fs::path &post_dir, const PostOperator &post_op,
     const SpaceOperator &space_op, int n_expected_rows, int delta_post_)
   : delta_post{delta_post_}, write_paraview_fields(delta_post_ > 0),
-    domains{true, root, post_dir, post_op.GetDomainPostOp(), "t (ns)", n_expected_rows},
+    domains{true, root, post_dir, post_op, "t (ns)", n_expected_rows},
     surfaces{true, root, post_dir, post_op, "t (ns)", n_expected_rows},
     currents{true, root, post_dir, space_op.GetSurfaceCurrentOp(), n_expected_rows},
     probes{true, root, post_dir, post_op, "t (ns)", n_expected_rows},
@@ -389,14 +389,12 @@ TransientSolver::PostprocessPrintResults::PostprocessPrintResults(
 
 void TransientSolver::PostprocessPrintResults::PostprocessStep(
     const IoData &iodata, const PostOperator &post_op, const SpaceOperator &space_op,
-    int step, double t, double J_coef, double E_elec, double E_mag)
+    int step, double t, double J_coef)
 {
   auto time = iodata.DimensionalizeValue(IoData::ValueType::TIME, t);
-  auto E_cap = post_op.GetLumpedCapacitorEnergy(space_op.GetLumpedPortOp());
-  auto E_ind = post_op.GetLumpedInductorEnergy(space_op.GetLumpedPortOp());
 
-  domains.AddMeasurement(time, post_op, E_elec, E_mag, E_cap, E_ind, iodata);
-  surfaces.AddMeasurement(time, post_op, E_elec + E_cap, E_mag + E_ind, iodata);
+  domains.AddMeasurement(time, post_op, iodata);
+  surfaces.AddMeasurement(time, post_op, iodata);
   currents.AddMeasurement(t, J_coef, space_op.GetSurfaceCurrentOp(), iodata);
   probes.AddMeasurement(time, post_op, iodata);
   ports.AddMeasurement(t, J_coef, post_op, space_op.GetLumpedPortOp(), iodata);
@@ -413,7 +411,8 @@ void TransientSolver::PostprocessPrintResults::PostprocessFinal(
     const PostOperator &post_op, const ErrorIndicator &indicator)
 {
   BlockTimer bt0(Timer::POSTPRO);
-  error_indicator.PrintIndicatorStatistics(post_op, indicator);
+  auto indicator_stats = indicator.GetSummaryStatistics(post_op.GetComm());
+  error_indicator.PrintIndicatorStatistics(post_op, indicator_stats);
   if (write_paraview_fields)
   {
     post_op.WriteFieldsFinal(&indicator);
diff --git a/palace/drivers/transientsolver.hpp b/palace/drivers/transientsolver.hpp
index ab2413c70..82d6edc97 100644
--- a/palace/drivers/transientsolver.hpp
+++ b/palace/drivers/transientsolver.hpp
@@ -75,8 +75,7 @@ class TransientSolver : public BaseSolver
                             const PostOperator &post_op, const SpaceOperator &space_op,
                             int n_expected_rows, int delta_post);
     void PostprocessStep(const IoData &iodata, const PostOperator &post_op,
-                         const SpaceOperator &space_op, int step, double t, double J_coef,
-                         double E_elec, double E_mag);
+                         const SpaceOperator &space_op, int step, double t, double J_coef);
     void PostprocessFinal(const PostOperator &post_op, const ErrorIndicator &indicator);
   };
 
diff --git a/palace/fem/errorindicator.hpp b/palace/fem/errorindicator.hpp
index c0f301848..23e27f046 100644
--- a/palace/fem/errorindicator.hpp
+++ b/palace/fem/errorindicator.hpp
@@ -60,6 +60,19 @@ class ErrorIndicator
 
   // Return the mean local error indicator.
   auto Mean(MPI_Comm comm) const { return linalg::Mean(comm, local); }
+
+  struct SummaryStatistics
+  {
+    double norm;
+    double min;
+    double max;
+    double mean;
+  };
+
+  SummaryStatistics GetSummaryStatistics(MPI_Comm comm) const
+  {
+    return {Norml2(comm), Min(comm), Max(comm), Mean(comm)};
+  }
 };
 
 }  // namespace palace
diff --git a/palace/fem/interpolator.cpp b/palace/fem/interpolator.cpp
index 01e3468f3..8353815ad 100644
--- a/palace/fem/interpolator.cpp
+++ b/palace/fem/interpolator.cpp
@@ -144,7 +144,7 @@ std::vector<std::complex<double>> InterpolationOperator::ProbeField(const GridFu
   }
   else
   {
-    return std::vector<std::complex<double>>(vr.begin(), vr.end());
+    return {vr.begin(), vr.end()};
   }
 }
 
diff --git a/palace/models/postoperator.cpp b/palace/models/postoperator.cpp
index b737ddec4..80b02819a 100644
--- a/palace/models/postoperator.cpp
+++ b/palace/models/postoperator.cpp
@@ -38,17 +38,17 @@ PostOperator::PostOperator(const IoData &iodata, SpaceOperator &space_op,
     surf_post_op(iodata, space_op.GetMaterialOp(), space_op.GetH1Space()),
     dom_post_op(iodata, space_op.GetMaterialOp(), space_op.GetNDSpace(),
                 space_op.GetRTSpace()),
+    interp_op(iodata, space_op.GetNDSpace()), lumped_port_op(&(space_op.GetLumpedPortOp())),
+    wave_port_op(&(space_op.GetWavePortOp())),
     E(std::make_unique<GridFunction>(space_op.GetNDSpace(),
                                      iodata.problem.type !=
                                          config::ProblemData::Type::TRANSIENT)),
     B(std::make_unique<GridFunction>(space_op.GetRTSpace(),
                                      iodata.problem.type !=
                                          config::ProblemData::Type::TRANSIENT)),
-    lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), &space_op.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
-                 &space_op.GetNDSpace().GetParMesh()),
-    interp_op(iodata, space_op.GetNDSpace())
+                 &space_op.GetNDSpace().GetParMesh())
 {
   U_e = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>>(*E, mat_op);
   U_m = std::make_unique<EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>>(*B, mat_op);
@@ -86,8 +86,7 @@ PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplace_op,
     surf_post_op(iodata, laplace_op.GetMaterialOp(), laplace_op.GetH1Space()),
     dom_post_op(iodata, laplace_op.GetMaterialOp(), laplace_op.GetH1Space()),
     E(std::make_unique<GridFunction>(laplace_op.GetNDSpace())),
-    V(std::make_unique<GridFunction>(laplace_op.GetH1Space())), lumped_port_init(false),
-    wave_port_init(false),
+    V(std::make_unique<GridFunction>(laplace_op.GetH1Space())),
     paraview(CreateParaviewPath(iodata, name), &laplace_op.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  &laplace_op.GetNDSpace().GetParMesh()),
@@ -113,8 +112,7 @@ PostOperator::PostOperator(const IoData &iodata, CurlCurlOperator &curlcurl_op,
     surf_post_op(iodata, curlcurl_op.GetMaterialOp(), curlcurl_op.GetH1Space()),
     dom_post_op(iodata, curlcurl_op.GetMaterialOp(), curlcurl_op.GetNDSpace()),
     B(std::make_unique<GridFunction>(curlcurl_op.GetRTSpace())),
-    A(std::make_unique<GridFunction>(curlcurl_op.GetNDSpace())), lumped_port_init(false),
-    wave_port_init(false),
+    A(std::make_unique<GridFunction>(curlcurl_op.GetNDSpace())),
     paraview(CreateParaviewPath(iodata, name), &curlcurl_op.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  &curlcurl_op.GetNDSpace().GetParMesh()),
@@ -274,7 +272,7 @@ void PostOperator::SetEGridFunction(const ComplexVector &e, bool exchange_face_n
     E->Real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces
     E->Imag().ExchangeFaceNbrData();
   }
-  lumped_port_init = wave_port_init = false;
+  ClearAllMeasurementCache();
 }
 
 void PostOperator::SetBGridFunction(const ComplexVector &b, bool exchange_face_nbr_data)
@@ -289,7 +287,7 @@ void PostOperator::SetBGridFunction(const ComplexVector &b, bool exchange_face_n
     B->Real().ExchangeFaceNbrData();  // Ready for parallel comm on shared faces
     B->Imag().ExchangeFaceNbrData();
   }
-  lumped_port_init = wave_port_init = false;
+  ClearAllMeasurementCache();
 }
 
 void PostOperator::SetEGridFunction(const Vector &e, bool exchange_face_nbr_data)
@@ -302,7 +300,7 @@ void PostOperator::SetEGridFunction(const Vector &e, bool exchange_face_nbr_data
   {
     E->Real().ExchangeFaceNbrData();
   }
-  lumped_port_init = wave_port_init = false;
+  ClearAllMeasurementCache();
 }
 
 void PostOperator::SetBGridFunction(const Vector &b, bool exchange_face_nbr_data)
@@ -315,7 +313,7 @@ void PostOperator::SetBGridFunction(const Vector &b, bool exchange_face_nbr_data
   {
     B->Real().ExchangeFaceNbrData();
   }
-  lumped_port_init = wave_port_init = false;
+  ClearAllMeasurementCache();
 }
 
 void PostOperator::SetVGridFunction(const Vector &v, bool exchange_face_nbr_data)
@@ -328,6 +326,7 @@ void PostOperator::SetVGridFunction(const Vector &v, bool exchange_face_nbr_data
   {
     V->Real().ExchangeFaceNbrData();
   }
+  ClearAllMeasurementCache();
 }
 
 void PostOperator::SetAGridFunction(const Vector &a, bool exchange_face_nbr_data)
@@ -340,66 +339,282 @@ void PostOperator::SetAGridFunction(const Vector &a, bool exchange_face_nbr_data
   {
     A->Real().ExchangeFaceNbrData();
   }
+  ClearAllMeasurementCache();
 }
 
-double PostOperator::GetEFieldEnergy() const
+// Measurements
+
+void PostOperator::ClearAllMeasurementCache()
 {
-  if (V)
+  // Clear Cache: Save omega since this set by hand like fields E,...
+  bool has_omega = measurment_cache.omega.has_value();
+  std::complex<double> omega;
+  if (has_omega)
   {
-    return dom_post_op.GetElectricFieldEnergy(*V);
+    omega = *measurment_cache.omega;
   }
-  else
+
+  measurment_cache = {};
+
+  if (has_omega)
   {
-    MFEM_VERIFY(E, "PostOperator is not configured for electric field energy calculation!");
-    return dom_post_op.GetElectricFieldEnergy(*E);
+    measurment_cache.omega = omega;
   }
 }
 
-double PostOperator::GetHFieldEnergy() const
+void PostOperator::MeasureAll()
 {
-  if (A)
+  ClearAllMeasurementCache();
+
+  // Domain Energy: Electric Field Contribution
+  if (V || E)
   {
-    return dom_post_op.GetMagneticFieldEnergy(*A);
+    GetEFieldEnergy();  // if (dom_post_op.M_elec)
+    if (dom_post_op.M_i.size() > 0)
+    {
+      GetEFieldEnergy(dom_post_op.M_i.begin()->first);  // Measures all domains
+    }
   }
-  else
+  // Domain Energy: Magnetic Field Contribution
+  if (A || B)
+  {
+    GetHFieldEnergy();  // if (dom_post_op.M_mag)
+    if (dom_post_op.M_i.size() > 0)
+    {
+      GetHFieldEnergy(dom_post_op.M_i.begin()->first);  // Measures all domains
+    }
+  }
+
+  if (E || B)
+  {
+    GetSurfaceFluxAll();
+  }
+  if (E)
+  {
+    GetInterfaceEFieldEnergyAll();
+    ProbeEField();
+  }
+  if (B)
+  {
+    ProbeBField();
+  }
+  if (E && B)
+  {
+    if (lumped_port_op != nullptr)
+    {
+      MeasureLumpedPorts();
+    }
+    if (wave_port_op != nullptr)
+    {
+      MeasureWavePorts();
+    }
+  }
+}
+
+void PostOperator::SetFrequency(double omega)
+{
+  measurment_cache.omega = std::complex<double>(omega);
+}
+
+void PostOperator::SetFrequency(std::complex<double> omega)
+{
+  measurment_cache.omega = omega;
+}
+
+std::complex<double> PostOperator::GetFrequency() const
+{
+  MFEM_VERIFY(measurment_cache.omega.has_value(),
+              "Frequency value omega has not been correctly set!");
+  return *measurment_cache.omega;
+}
+
+double PostOperator::GetEFieldEnergy() const
+{
+  if (!measurment_cache.domain_E_field_energy_all.has_value())
+  {
+    if (V)
+    {
+      measurment_cache.domain_E_field_energy_all = dom_post_op.GetElectricFieldEnergy(*V);
+    }
+    else if (E)
+    {
+      measurment_cache.domain_E_field_energy_all = dom_post_op.GetElectricFieldEnergy(*E);
+    }
+    else
+    {
+      // No failure: returns zero
+      measurment_cache.domain_E_field_energy_all = 0.0;
+    }
+  }
+  return *measurment_cache.domain_E_field_energy_all;
+}
+
+double PostOperator::GetHFieldEnergy() const
+{
+  if (!measurment_cache.domain_H_field_energy_all.has_value())
   {
-    MFEM_VERIFY(B, "PostOperator is not configured for magnetic field energy calculation!");
-    return dom_post_op.GetMagneticFieldEnergy(*B);
+    if (A)
+    {
+      measurment_cache.domain_H_field_energy_all = dom_post_op.GetMagneticFieldEnergy(*A);
+    }
+    else if (B)
+    {
+      measurment_cache.domain_H_field_energy_all = dom_post_op.GetMagneticFieldEnergy(*B);
+    }
+    else
+    {
+      // No failure: returns zero
+      measurment_cache.domain_H_field_energy_all = 0.0;
+    }
   }
+  return *measurment_cache.domain_H_field_energy_all;
 }
 
 double PostOperator::GetEFieldEnergy(int idx) const
 {
-  if (V)
+  if (!measurment_cache.domain_E_field_energy_i.has_value())
   {
-    return dom_post_op.GetDomainElectricFieldEnergy(idx, *V);
+    // Do all measurements
+    measurment_cache.domain_E_field_energy_i.emplace();
+    for (const auto &[idx, data] : dom_post_op.M_i)
+    {
+      double out;
+      if (V)
+      {
+        out = dom_post_op.GetDomainElectricFieldEnergy(idx, *V);
+      }
+      else if (E)
+      {
+        out = dom_post_op.GetDomainElectricFieldEnergy(idx, *E);
+      }
+      else
+      {
+        // No failure: returns zero
+        out = 0.0;
+      }
+      measurment_cache.domain_E_field_energy_i->emplace(idx, out);
+    }
   }
-  else
+  auto it = measurment_cache.domain_E_field_energy_i->find(idx);
+  if (it == measurment_cache.domain_E_field_energy_i->end())
   {
-    MFEM_VERIFY(E, "PostOperator is not configured for electric field energy calculation!");
-    return dom_post_op.GetDomainElectricFieldEnergy(idx, *E);
+    MFEM_ABORT(fmt::format("Could not find domain index {} for E field energy!", idx));
   }
+  return it->second;
 }
 
 double PostOperator::GetHFieldEnergy(int idx) const
 {
-  if (A)
+  if (!measurment_cache.domain_H_field_energy_i.has_value())
   {
-    return dom_post_op.GetDomainMagneticFieldEnergy(idx, *A);
+    // Do all measurements
+    measurment_cache.domain_H_field_energy_i.emplace();
+    for (const auto &[idx, data] : dom_post_op.M_i)
+    {
+      double out;
+      if (A)
+      {
+        out = dom_post_op.GetDomainMagneticFieldEnergy(idx, *A);
+      }
+      else if (B)
+      {
+        out = dom_post_op.GetDomainMagneticFieldEnergy(idx, *B);
+      }
+      else
+      {
+        // No failure: returns zero
+        out = 0.0;
+      }
+      measurment_cache.domain_H_field_energy_i->emplace(idx, out);
+    }
   }
-  else
+  auto it = measurment_cache.domain_H_field_energy_i->find(idx);
+  if (it == measurment_cache.domain_H_field_energy_i->end())
   {
-    MFEM_VERIFY(B, "PostOperator is not configured for magnetic field energy calculation!");
-    return dom_post_op.GetDomainMagneticFieldEnergy(idx, *B);
+    MFEM_ABORT(fmt::format("Could not find domain index {} for H field energy!", idx));
   }
+  return it->second;
 }
 
-std::complex<double> PostOperator::GetSurfaceFlux(int idx) const
+// Code Note: for returning the full vector we chose name GetSurfaceFluxAll() rather than
+// GetSurfaceFlux(), to keep consistancy with GetEFieldEnergy(). GetEFieldEnergy() returns
+// the total energy (calculated differently), not the vector of the individual
+// GetHFieldEnergy(int idx).
+std::vector<PostOperator::FluxData> PostOperator::GetSurfaceFluxAll() const
 {
   // Compute the flux through a surface as Φ_j = ∫ F ⋅ n_j dS, with F = B, F = ε D, or F =
   // E x H. The special coefficient is used to avoid issues evaluating MFEM GridFunctions
   // which are discontinuous at interior boundary elements.
-  return surf_post_op.GetSurfaceFlux(idx, E.get(), B.get());
+  if (!measurment_cache.surface_flux_i.has_value())
+  {
+    // Do all measurements
+    MFEM_VERIFY(
+        E || B,
+        "PostOperator needs either electric or magnetic field for flux calculation!");
+    measurment_cache.surface_flux_i.emplace();
+    measurment_cache.surface_flux_i->reserve(surf_post_op.flux_surfs.size());
+    for (const auto &[idx, data] : surf_post_op.flux_surfs)
+    {
+      measurment_cache.surface_flux_i->emplace_back(
+          FluxData{idx, surf_post_op.GetSurfaceFlux(idx, E.get(), B.get()), data.type});
+    }
+  }
+  return *measurment_cache.surface_flux_i;
+}
+
+PostOperator::FluxData PostOperator::GetSurfaceFlux(int idx) const
+{
+  if (!measurment_cache.surface_flux_i.has_value())
+  {
+    GetSurfaceFluxAll();
+  }
+  auto it = std::find_if(measurment_cache.surface_flux_i->begin(),
+                         measurment_cache.surface_flux_i->end(),
+                         [idx](const auto &d) { return d.idx == idx; });
+  if (it == measurment_cache.surface_flux_i->end())
+  {
+    MFEM_ABORT(fmt::format("Could not find surface index {} for flux!", idx));
+  }
+  return *it;
+}
+
+std::vector<PostOperator::InterfaceData> PostOperator::GetInterfaceEFieldEnergyAll() const
+{
+  // Compute the surface dielectric participation ratio and associated quality factor for
+  // the material interface given by index idx. We have:
+  //                            1/Q_mj = p_mj tan(δ)_j
+  // with:
+  //          p_mj = 1/2 t_j Re{∫_{Γ_j} (ε_j E_m)ᴴ E_m dS} /(E_elec + E_cap).
+  if (!measurment_cache.interface_eps_i.has_value())
+  {
+    // Do all measurements
+    MFEM_VERIFY(E, "Electric field solution required for E field interface energy!");
+    measurment_cache.interface_eps_i.emplace();
+    measurment_cache.interface_eps_i->reserve(surf_post_op.eps_surfs.size());
+    for (const auto &[idx, data] : surf_post_op.eps_surfs)
+    {
+      measurment_cache.interface_eps_i->emplace_back(
+          InterfaceData{idx, surf_post_op.GetInterfaceElectricFieldEnergy(idx, *E),
+                        surf_post_op.GetInterfaceLossTangent(idx)});
+    }
+  }
+  return *measurment_cache.interface_eps_i;
+}
+
+PostOperator::InterfaceData PostOperator::GetInterfaceEFieldEnergy(int idx) const
+{
+  if (!measurment_cache.interface_eps_i.has_value())
+  {
+    GetInterfaceEFieldEnergyAll();
+  }
+  auto it = std::find_if(measurment_cache.interface_eps_i->begin(),
+                         measurment_cache.interface_eps_i->end(),
+                         [idx](const auto &d) { return d.idx == idx; });
+  if (it == measurment_cache.interface_eps_i->end())
+  {
+    MFEM_ABORT(fmt::format("Could not find surface index {} for interface energy!", idx));
+  }
+  return *it;
 }
 
 double PostOperator::GetInterfaceParticipation(int idx, double E_m) const
@@ -410,25 +625,33 @@ double PostOperator::GetInterfaceParticipation(int idx, double E_m) const
   // with:
   //          p_mj = 1/2 t_j Re{∫_{Γ_j} (ε_j E_m)ᴴ E_m dS} /(E_elec + E_cap).
   MFEM_VERIFY(E, "Surface Q not defined, no electric field solution found!");
-  return surf_post_op.GetInterfaceElectricFieldEnergy(idx, *E) / E_m;
+  auto data = GetInterfaceEFieldEnergy(idx);
+  return data.energy / E_m;
 }
 
-void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega)
+void PostOperator::MeasureLumpedPorts() const
 {
-  MFEM_VERIFY(E && B, "Incorrect usage of PostOperator::UpdatePorts!");
-  if (lumped_port_init)
+  MFEM_VERIFY(E && B && lumped_port_op != nullptr,
+              "Incorrect usage of PostOperator::MeasureLumpedPorts!");
+  if (measurment_cache.lumped_port_vi.has_value())
   {
-    return;
+    measurment_cache.lumped_port_vi->clear();
   }
-  for (const auto &[idx, data] : lumped_port_op)
+  else
+  {
+    measurment_cache.lumped_port_vi.emplace();
+  }
+  for (const auto &[idx, data] : *lumped_port_op)
   {
-    auto &vi = lumped_port_vi[idx];
+    auto &vi = (*measurment_cache.lumped_port_vi)[idx];
     vi.P = data.GetPower(*E, *B);
     vi.V = data.GetVoltage(*E);
     if (HasImag())
     {
       // Compute current from the port impedance, separate contributions for R, L, C
       // branches.
+      // Get value and make real: Matches current behaviour (even for eigensolver!)
+      auto omega = GetFrequency().real();
       MFEM_VERIFY(
           omega > 0.0,
           "Frequency domain lumped port postprocessing requires nonzero frequency!");
@@ -453,182 +676,238 @@ void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double
       vi.I[1] = vi.I[2] = vi.S = 0.0;
     }
   }
-  lumped_port_init = true;
 }
 
-void PostOperator::UpdatePorts(const WavePortOperator &wave_port_op, double omega)
+void PostOperator::MeasureWavePorts() const
 {
-  MFEM_VERIFY(HasImag() && E && B, "Incorrect usage of PostOperator::UpdatePorts!");
-  if (wave_port_init)
+
+  MFEM_VERIFY(E && B && wave_port_op != nullptr,
+              "Incorrect usage of PostOperator::MeasureWavePorts!");
+  if (measurment_cache.wave_port_vi.has_value())
   {
-    return;
+    measurment_cache.wave_port_vi->clear();
   }
-  for (const auto &[idx, data] : wave_port_op)
+  else
   {
+    measurment_cache.wave_port_vi.emplace();
+  }
+  if (!HasImag())
+  {
+    return;  // Wave ports need Imag; leave empty otherwise // TODO: Fix in long run
+  }
+  for (const auto &[idx, data] : *wave_port_op)
+  {
+    MFEM_VERIFY(measurment_cache.omega.has_value(),
+                "Measuring port currents with Imag fields, requires frequency to be set "
+                "with SetFrequency!");
+
+    // Get value and make real: Matches current behaviour
+    auto omega = measurment_cache.omega->real();
+
     MFEM_VERIFY(omega > 0.0,
                 "Frequency domain wave port postprocessing requires nonzero frequency!");
-    auto &vi = wave_port_vi[idx];
+    auto &vi = (*measurment_cache.wave_port_vi)[idx];
     vi.P = data.GetPower(*E, *B);
     vi.S = data.GetSParameter(*E);
     vi.V = vi.I[0] = vi.I[1] = vi.I[2] = 0.0;  // Not yet implemented
                                                // (Z = V² / P, I = V / Z)
   }
-  wave_port_init = true;
 }
 
-double PostOperator::GetLumpedInductorEnergy(const LumpedPortOperator &lumped_port_op) const
+void PostOperator::ValidateDoPortMeasurement() const
+{
+  if (!measurment_cache.lumped_port_vi.has_value())
+  {
+    if (lumped_port_op != nullptr)
+    {
+      MeasureLumpedPorts();
+    }
+    else
+    {
+      MFEM_ABORT("A lumped port measurement called, but the lumped port operator is not "
+                 "defined by the solver.")
+    }
+  }
+  if (!measurment_cache.wave_port_vi.has_value())
+  {
+    if (wave_port_op != nullptr)
+    {
+      MeasureWavePorts();
+    }
+    else
+    {
+      MFEM_ABORT("A wave port measurement called, but the wave port operator is not "
+                 "defined by the solver.")
+    }
+  }
+}
+
+double PostOperator::GetLumpedInductorEnergy() const
 {
   // Add contribution due to all inductive lumped boundaries in the model:
   //                      E_ind = ∑_j 1/2 L_j I_mj².
-  double U = 0.0;
-  for (const auto &[idx, data] : lumped_port_op)
+  if (!measurment_cache.lumped_port_inductor_energy.has_value())
   {
-    if (std::abs(data.L) > 0.0)
+    // No failure if space has no lumped ports: Returns zero
+    double U = 0.0;
+    if (lumped_port_op != nullptr)
     {
-      std::complex<double> I_j =
-          GetPortCurrent(lumped_port_op, idx, LumpedPortData::Branch::L);
-      U += 0.5 * std::abs(data.L) * std::real(I_j * std::conj(I_j));
+      for (const auto &[idx, data] : *lumped_port_op)
+      {
+        if (std::abs(data.L) > 0.0)
+        {
+          std::complex<double> I_j = GetPortCurrent(idx, LumpedPortData::Branch::L);
+          U += 0.5 * std::abs(data.L) * std::real(I_j * std::conj(I_j));
+        }
+      }
     }
+    measurment_cache.lumped_port_inductor_energy = U;
   }
-  return U;
+  return *measurment_cache.lumped_port_inductor_energy;
 }
 
-double
-PostOperator::GetLumpedCapacitorEnergy(const LumpedPortOperator &lumped_port_op) const
+double PostOperator::GetLumpedCapacitorEnergy() const
 {
   // Add contribution due to all capacitive lumped boundaries in the model:
   //                      E_cap = ∑_j 1/2 C_j V_mj².
-  double U = 0.0;
-  for (const auto &[idx, data] : lumped_port_op)
+  if (!measurment_cache.lumped_port_capacitor_energy.has_value())
   {
-    if (std::abs(data.C) > 0.0)
+    // No failure if space has no lumped ports: Returns zero
+    double U = 0.0;
+    if (lumped_port_op != nullptr)
     {
-      std::complex<double> V_j = GetPortVoltage(lumped_port_op, idx);
-      U += 0.5 * std::abs(data.C) * std::real(V_j * std::conj(V_j));
+      for (const auto &[idx, data] : *lumped_port_op)
+      {
+        if (std::abs(data.C) > 0.0)
+        {
+          std::complex<double> V_j = GetPortVoltage(idx);
+          U += 0.5 * std::abs(data.C) * std::real(V_j * std::conj(V_j));
+        }
+      }
     }
+    measurment_cache.lumped_port_capacitor_energy = U;
   }
-  return U;
+  return *measurment_cache.lumped_port_capacitor_energy;
 }
 
-std::complex<double> PostOperator::GetSParameter(const LumpedPortOperator &lumped_port_op,
-                                                 int idx, int source_idx) const
+std::complex<double> PostOperator::GetSParameter(bool is_lumped_port, int idx,
+                                                 int source_idx) const
 {
-  MFEM_VERIFY(lumped_port_init,
-              "Port S-parameters not defined until ports are initialized!");
-  const LumpedPortData &data = lumped_port_op.GetPort(idx);
-  const LumpedPortData &src_data = lumped_port_op.GetPort(source_idx);
-  const auto it = lumped_port_vi.find(idx);
-  MFEM_VERIFY(src_data.excitation,
-              "Lumped port index " << source_idx << " is not marked for excitation!");
-  MFEM_VERIFY(it != lumped_port_vi.end(),
-              "Could not find lumped port when calculating port S-parameters!");
-  std::complex<double> S_ij = it->second.S;
-  if (idx == source_idx)
+  ValidateDoPortMeasurement();
+  // TODO: In multi-excittion PR we will gurantee that lumped & wave ports have unique idx
+  // TODO: Merge lumped and wave port S_ij calcluations to allow both at same time.
+  if (is_lumped_port)
   {
-    S_ij.real(S_ij.real() - 1.0);
+    const LumpedPortData &data = lumped_port_op->GetPort(idx);
+    const LumpedPortData &src_data = lumped_port_op->GetPort(source_idx);
+    const auto it = measurment_cache.lumped_port_vi->find(idx);
+    MFEM_VERIFY(src_data.excitation,
+                "Lumped port index " << source_idx << " is not marked for excitation!");
+    MFEM_VERIFY(it != measurment_cache.lumped_port_vi->end(),
+                "Could not find lumped port when calculating port S-parameters!");
+    std::complex<double> S_ij = it->second.S;
+    if (idx == source_idx)
+    {
+      S_ij.real(S_ij.real() - 1.0);
+    }
+    // Generalized S-parameters if the ports are resistive (avoids divide-by-zero).
+    if (std::abs(data.R) > 0.0)
+    {
+      S_ij *= std::sqrt(src_data.R / data.R);
+    }
+    return S_ij;
   }
-  // Generalized S-parameters if the ports are resistive (avoids divide-by-zero).
-  if (std::abs(data.R) > 0.0)
+  else
   {
-    S_ij *= std::sqrt(src_data.R / data.R);
+    // Wave port modes are not normalized to a characteristic impedance so no generalized
+    // S-parameters are available.
+    const WavePortData &data = wave_port_op->GetPort(idx);
+    const WavePortData &src_data = wave_port_op->GetPort(source_idx);
+    const auto it = measurment_cache.wave_port_vi->find(idx);
+    MFEM_VERIFY(src_data.excitation,
+                "Wave port index " << source_idx << " is not marked for excitation!");
+    MFEM_VERIFY(it != measurment_cache.wave_port_vi->end(),
+                "Could not find wave port when calculating port S-parameters!");
+    std::complex<double> S_ij = it->second.S;
+    if (idx == source_idx)
+    {
+      S_ij.real(S_ij.real() - 1.0);
+    }
+    // Port de-embedding: S_demb = S exp(ikₙᵢ dᵢ) exp(ikₙⱼ dⱼ) (distance offset is default 0
+    // unless specified).
+    S_ij *= std::exp(1i * src_data.kn0 * src_data.d_offset);
+    S_ij *= std::exp(1i * data.kn0 * data.d_offset);
+    return S_ij;
   }
-  return S_ij;
 }
 
-std::complex<double> PostOperator::GetSParameter(const WavePortOperator &wave_port_op,
-                                                 int idx, int source_idx) const
+std::complex<double> PostOperator::GetPortPower(int idx) const
 {
-  // Wave port modes are not normalized to a characteristic impedance so no generalized
-  // S-parameters are available.
-  MFEM_VERIFY(wave_port_init, "Port S-parameters not defined until ports are initialized!");
-  const WavePortData &data = wave_port_op.GetPort(idx);
-  const WavePortData &src_data = wave_port_op.GetPort(source_idx);
-  const auto it = wave_port_vi.find(idx);
-  MFEM_VERIFY(src_data.excitation,
-              "Wave port index " << source_idx << " is not marked for excitation!");
-  MFEM_VERIFY(it != wave_port_vi.end(),
-              "Could not find wave port when calculating port S-parameters!");
-  std::complex<double> S_ij = it->second.S;
-  if (idx == source_idx)
+  ValidateDoPortMeasurement();
+  // TODO: In multi-excittion PR we will gurantee that lumped & wave ports have unique idx
+  auto it_lumped = measurment_cache.lumped_port_vi->find(idx);
+  if (it_lumped != measurment_cache.lumped_port_vi->end())
   {
-    S_ij.real(S_ij.real() - 1.0);
+    return it_lumped->second.P;
   }
-  // Port de-embedding: S_demb = S exp(ikₙᵢ dᵢ) exp(ikₙⱼ dⱼ) (distance offset is default 0
-  // unless specified).
-  S_ij *= std::exp(1i * src_data.kn0 * src_data.d_offset);
-  S_ij *= std::exp(1i * data.kn0 * data.d_offset);
-  return S_ij;
-}
-
-std::complex<double> PostOperator::GetPortPower(const LumpedPortOperator &lumped_port_op,
-                                                int idx) const
-{
-  MFEM_VERIFY(lumped_port_init,
-              "Lumped port quantities not defined until ports are initialized!");
-  const auto it = lumped_port_vi.find(idx);
-  MFEM_VERIFY(it != lumped_port_vi.end(),
-              "Could not find lumped port when calculating lumped port power!");
-  return it->second.P;
-}
-
-std::complex<double> PostOperator::GetPortPower(const WavePortOperator &wave_port_op,
-                                                int idx) const
-{
-  MFEM_VERIFY(wave_port_init,
-              "Wave port quantities not defined until ports are initialized!");
-  const auto it = wave_port_vi.find(idx);
-  MFEM_VERIFY(it != wave_port_vi.end(),
-              "Could not find wave port when calculating wave port power!");
-  return it->second.P;
-}
-
-std::complex<double> PostOperator::GetPortVoltage(const LumpedPortOperator &lumped_port_op,
-                                                  int idx) const
-{
-  MFEM_VERIFY(lumped_port_init,
-              "Lumped port quantities not defined until ports are initialized!");
-  const auto it = lumped_port_vi.find(idx);
-  MFEM_VERIFY(it != lumped_port_vi.end(),
-              "Could not find lumped port when calculating lumped port voltage!");
-  return it->second.V;
+  auto it_wave = measurment_cache.wave_port_vi->find(idx);
+  if (it_wave != measurment_cache.wave_port_vi->end())
+  {
+    return it_wave->second.P;
+  }
+  MFEM_ABORT(
+      fmt::format("Port Power: Could not find a lumped or wave port with index {}!", idx));
 }
 
-std::complex<double> PostOperator::GetPortVoltage(const WavePortOperator &wave_port_op,
-                                                  int idx) const
+std::complex<double> PostOperator::GetPortVoltage(int idx) const
 {
-  MFEM_ABORT("GetPortVoltage is not yet implemented for wave port boundaries!");
-  return 0.0;
+  ValidateDoPortMeasurement();
+  // TODO: In multi-excittion PR we will gurantee that lumped & wave ports have unique idx
+  auto it_lumped = measurment_cache.lumped_port_vi->find(idx);
+  if (it_lumped != measurment_cache.lumped_port_vi->end())
+  {
+    return it_lumped->second.V;
+  }
+  auto it_wave = measurment_cache.wave_port_vi->find(idx);
+  if (it_wave != measurment_cache.wave_port_vi->end())
+  {
+    MFEM_ABORT("GetPortVoltage is not yet implemented for wave port boundaries!");
+  }
+  MFEM_ABORT(fmt::format(
+      "Port Voltage: Could not find a lumped or wave port with index {}!", idx));
 }
 
-std::complex<double> PostOperator::GetPortCurrent(const LumpedPortOperator &lumped_port_op,
-                                                  int idx,
+std::complex<double> PostOperator::GetPortCurrent(int idx,
                                                   LumpedPortData::Branch branch) const
 {
-  MFEM_VERIFY(lumped_port_init,
-              "Lumped port quantities not defined until ports are initialized!");
-  const auto it = lumped_port_vi.find(idx);
-  MFEM_VERIFY(it != lumped_port_vi.end(),
-              "Could not find lumped port when calculating lumped port current!");
-  return ((branch == LumpedPortData::Branch::TOTAL || branch == LumpedPortData::Branch::R)
-              ? it->second.I[0]
-              : 0.0) +
-         ((branch == LumpedPortData::Branch::TOTAL || branch == LumpedPortData::Branch::L)
-              ? it->second.I[1]
-              : 0.0) +
-         ((branch == LumpedPortData::Branch::TOTAL || branch == LumpedPortData::Branch::C)
-              ? it->second.I[2]
-              : 0.0);
-}
-
-std::complex<double> PostOperator::GetPortCurrent(const WavePortOperator &wave_port_op,
-                                                  int idx) const
-{
-  MFEM_ABORT("GetPortCurrent is not yet implemented for wave port boundaries!");
-  return 0.0;
+  ValidateDoPortMeasurement();
+  // TODO: In multi-excittion PR we will gurantee that lumped & wave ports have unique idx
+  auto it_lumped = measurment_cache.lumped_port_vi->find(idx);
+  if (it_lumped != measurment_cache.lumped_port_vi->end())
+  {
+    auto &I_loc = it_lumped->second.I;
+    switch (branch)
+    {
+      case LumpedPortData::Branch::R:
+        return I_loc[0];
+      case LumpedPortData::Branch::L:
+        return I_loc[1];
+      case LumpedPortData::Branch::C:
+        return I_loc[2];
+      default:
+        return std::accumulate(I_loc.begin(), I_loc.end(), std::complex<double>{0.0, 0.0});
+    }
+  }
+  auto it_wave = measurment_cache.wave_port_vi->find(idx);
+  if (it_wave != measurment_cache.wave_port_vi->end())
+  {
+    MFEM_ABORT("GetPortCurrent is not yet implemented for wave port boundaries!");
+  }
+  MFEM_ABORT(fmt::format(
+      "Port Current: Could not find a lumped or wave port with index {}!", idx));
 }
 
-double PostOperator::GetInductorParticipation(const LumpedPortOperator &lumped_port_op,
-                                              int idx, double E_m) const
+double PostOperator::GetInductorParticipation(int idx, double E_m) const
 {
   // Compute energy-participation ratio of junction given by index idx for the field mode.
   // We first get the port line voltage, and use lumped port circuit impedance to get peak
@@ -638,25 +917,30 @@ double PostOperator::GetInductorParticipation(const LumpedPortOperator &lumped_p
   //                            p_mj = 1/2 L_j I_mj² / E_m.
   // An element with no assigned inductance will be treated as having zero admittance and
   // thus zero current.
-  const LumpedPortData &data = lumped_port_op.GetPort(idx);
-  std::complex<double> I_mj =
-      GetPortCurrent(lumped_port_op, idx, LumpedPortData::Branch::L);
+  if (lumped_port_op == nullptr)
+  {
+    return 0.0;
+  }
+  const LumpedPortData &data = lumped_port_op->GetPort(idx);
+  std::complex<double> I_mj = GetPortCurrent(idx, LumpedPortData::Branch::L);
   return std::copysign(0.5 * std::abs(data.L) * std::real(I_mj * std::conj(I_mj)) / E_m,
                        I_mj.real());  // mean(I²) = (I_r² + I_i²) / 2
 }
 
-double PostOperator::GetExternalKappa(const LumpedPortOperator &lumped_port_op, int idx,
-                                      double E_m) const
+double PostOperator::GetExternalKappa(int idx, double E_m) const
 {
-  // Compute participation ratio of external ports (given as any port boundary with nonzero
-  // resistance). Currently no reactance of the ports is supported. The κ of the port
-  // follows from:
+  // Compute participation ratio of external ports (given as any port boundary with
+  // nonzero resistance). Currently no reactance of the ports is supported. The κ of the
+  // port follows from:
   //                          κ_mj = 1/2 R_j I_mj² / E_m
   // from which the mode coupling quality factor is computed as:
   //                              Q_mj = ω_m / κ_mj.
-  const LumpedPortData &data = lumped_port_op.GetPort(idx);
-  std::complex<double> I_mj =
-      GetPortCurrent(lumped_port_op, idx, LumpedPortData::Branch::R);
+  if (lumped_port_op == nullptr)
+  {
+    return 0.0;
+  }
+  const LumpedPortData &data = lumped_port_op->GetPort(idx);
+  std::complex<double> I_mj = GetPortCurrent(idx, LumpedPortData::Branch::R);
   return std::copysign(0.5 * std::abs(data.R) * std::real(I_mj * std::conj(I_mj)) / E_m,
                        I_mj.real());  // mean(I²) = (I_r² + I_i²) / 2
 }
@@ -811,14 +1095,35 @@ void PostOperator::WriteFieldsFinal(const ErrorIndicator *indicator) const
 
 std::vector<std::complex<double>> PostOperator::ProbeEField() const
 {
-  MFEM_VERIFY(E, "PostOperator is not configured for electric field probes!");
-  return interp_op.ProbeField(*E);
+  if (!measurment_cache.probe_E_field.has_value())
+  {
+    MFEM_VERIFY(E, "PostOperator is not configured for electric field probes!");
+    if (interp_op.GetProbes().size() > 0)
+    {
+      measurment_cache.probe_E_field = interp_op.ProbeField(*E);
+    }
+    else
+    {
+      measurment_cache.probe_E_field.emplace();
+    }
+  }
+  return *measurment_cache.probe_E_field;
 }
-
 std::vector<std::complex<double>> PostOperator::ProbeBField() const
 {
-  MFEM_VERIFY(B, "PostOperator is not configured for magnetic flux density probes!");
-  return interp_op.ProbeField(*B);
+  if (!measurment_cache.probe_B_field.has_value())
+  {
+    MFEM_VERIFY(B, "PostOperator is not configured for magnetic flux density probes!");
+    if (interp_op.GetProbes().size() > 0)
+    {
+      measurment_cache.probe_B_field = interp_op.ProbeField(*B);
+    }
+    else
+    {
+      measurment_cache.probe_B_field.emplace();
+    }
+  }
+  return *measurment_cache.probe_B_field;
 }
 
 }  // namespace palace
diff --git a/palace/models/postoperator.hpp b/palace/models/postoperator.hpp
index b21e79fb0..41e0d95c2 100644
--- a/palace/models/postoperator.hpp
+++ b/palace/models/postoperator.hpp
@@ -39,15 +39,26 @@ class PostOperator
   // Reference to material property operator (not owned).
   const MaterialOperator &mat_op;
 
-  // Surface boundary and domain postprocessors.
-  const SurfacePostOperator surf_post_op;
-  const DomainPostOperator dom_post_op;
-
   // Objects for grid function postprocessing from the FE solution.
-  mutable std::unique_ptr<GridFunction> E, B, V, A;
+  std::unique_ptr<GridFunction> E, B, V, A;
   std::unique_ptr<mfem::VectorCoefficient> S, E_sr, E_si, B_sr, B_si, A_s, J_sr, J_si;
   std::unique_ptr<mfem::Coefficient> U_e, U_m, V_s, Q_sr, Q_si;
 
+  // Data collection for writing fields to disk for visualization
+  // Paraview fields are mutable as the writing is triggered by const solver printer
+  mutable mfem::ParaViewDataCollection paraview, paraview_bdr;
+  double mesh_Lc0;
+
+  // ----- Measurements from Fields -----
+
+  DomainPostOperator dom_post_op;           // Energy in bulk
+  SurfacePostOperator surf_post_op;         //  Dielectric Interface Energy and Flux
+  mutable InterpolationOperator interp_op;  // E & B fields: mutates during measure
+
+  // Port Contributions: not owned, view onto space_op only, it must not go out of scope
+  LumpedPortOperator *lumped_port_op = nullptr;
+  WavePortOperator *wave_port_op = nullptr;
+
   // Wave port boundary mode field postprocessing.
   struct WavePortFieldData
   {
@@ -55,19 +66,56 @@ class PostOperator
   };
   std::map<int, WavePortFieldData> port_E0;
 
-  // Lumped and wave port voltage and current (R, L, and C branches) caches updated when
-  // the grid functions are set.
+public:
+  // Mini storage clases for cache: Definitions are public
+  struct FluxData
+  {
+    int idx;                   // Surface index
+    std::complex<double> Phi;  // Integrated flux
+    SurfaceFluxType type;      // Flux type
+  };
+
+  struct InterfaceData
+  {
+    int idx;          // Interface index
+    double energy;    // Surface ELectric Field Energy
+    double tandelta;  // Dissipation tangent tan(δ)
+  };
+
+  // For both lumped and wave port
   struct PortPostData
   {
-    std::complex<double> P, V, I[3], S;
+    std::complex<double> P, V, S;
+    std::array<std::complex<double>, 3> I;  // Separate R, L, and C branches
   };
-  std::map<int, PortPostData> lumped_port_vi, wave_port_vi;
-  bool lumped_port_init, wave_port_init;
 
-  // Data collection for writing fields to disk for visualization and sampling points.
-  mutable mfem::ParaViewDataCollection paraview, paraview_bdr;
-  mutable InterpolationOperator interp_op;
-  double mesh_Lc0;
+private:
+  struct MeasurementCache
+  {
+    std::optional<std::complex<double>> omega = std::nullopt;
+
+    std::optional<double> domain_E_field_energy_all = std::nullopt;
+    std::optional<double> domain_H_field_energy_all = std::nullopt;
+
+    std::optional<std::map<int, double>> domain_E_field_energy_i = std::nullopt;
+    std::optional<std::map<int, double>> domain_H_field_energy_i = std::nullopt;
+
+    std::optional<std::vector<FluxData>> surface_flux_i = std::nullopt;
+    std::optional<std::vector<InterfaceData>> interface_eps_i = std::nullopt;
+
+    std::optional<std::map<int, PortPostData>> lumped_port_vi = std::nullopt;
+    std::optional<std::map<int, PortPostData>> wave_port_vi = std::nullopt;
+
+    std::optional<double> lumped_port_inductor_energy = std::nullopt;
+    std::optional<double> lumped_port_capacitor_energy = std::nullopt;
+
+    std::optional<std::vector<std::complex<double>>> probe_E_field = std::nullopt;
+    std::optional<std::vector<std::complex<double>>> probe_B_field = std::nullopt;
+  };
+  mutable MeasurementCache measurment_cache = {};
+
+  void ValidateDoPortMeasurement() const;
+
   void InitializeDataCollection(const IoData &iodata);
 
 public:
@@ -118,6 +166,22 @@ class PostOperator
     return *A;
   }
 
+  // Function that triggers all available post-processing measurements and populate cache.
+  void MeasureAll();
+
+  // Clear internal measurement caches
+  void ClearAllMeasurementCache();
+
+  // Treat the frequency, for driven and eigemode solvers, as a "measurement", that other
+  // measurements can depend on. This has to be supplied during the solver loop separate
+  // from the fields.
+  void SetFrequency(double omega);
+  void SetFrequency(std::complex<double> omega);
+
+  // Return stored frequency that was given in SetFrequency. Always promotes to complex
+  // frequency.
+  std::complex<double> GetFrequency() const;
+
   // Postprocess the total electric and magnetic field energies in the electric and magnetic
   // fields.
   double GetEFieldEnergy() const;
@@ -130,55 +194,42 @@ class PostOperator
 
   // Postprocess the electric or magnetic field flux for a surface index using the computed
   // electcric field and/or magnetic flux density field solutions.
-  std::complex<double> GetSurfaceFlux(int idx) const;
+  std::vector<FluxData> GetSurfaceFluxAll() const;
+  FluxData GetSurfaceFlux(int idx) const;
 
   // Postprocess the partitipation ratio for interface lossy dielectric losses in the
   // electric field mode.
   double GetInterfaceParticipation(int idx, double E_m) const;
+  std::vector<InterfaceData> GetInterfaceEFieldEnergyAll() const;
+  InterfaceData GetInterfaceEFieldEnergy(int idx) const;
 
-  // Update cached port voltages and currents for lumped and wave port operators.
-  void UpdatePorts(const LumpedPortOperator &lumped_port_op,
-                   const WavePortOperator &wave_port_op, double omega = 0.0)
-  {
-    UpdatePorts(lumped_port_op, omega);
-    UpdatePorts(wave_port_op, omega);
-  }
-  void UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega = 0.0);
-  void UpdatePorts(const WavePortOperator &wave_port_op, double omega = 0.0);
+  // Measure and cache port voltages and currents for lumped and wave port operators.
+  void MeasureLumpedPorts() const;
+  void MeasureWavePorts() const;
 
   // Postprocess the energy in lumped capacitor or inductor port boundaries with index in
   // the provided set.
-  double GetLumpedInductorEnergy(const LumpedPortOperator &lumped_port_op) const;
-  double GetLumpedCapacitorEnergy(const LumpedPortOperator &lumped_port_op) const;
+  double GetLumpedInductorEnergy() const;
+  double GetLumpedCapacitorEnergy() const;
 
   // Postprocess the S-parameter for recieving lumped or wave port index using the electric
   // field solution.
-  std::complex<double> GetSParameter(const LumpedPortOperator &lumped_port_op, int idx,
-                                     int source_idx) const;
-  std::complex<double> GetSParameter(const WavePortOperator &wave_port_op, int idx,
-                                     int source_idx) const;
+  std::complex<double> GetSParameter(bool is_lumped_port, int idx, int source_idx) const;
 
   // Postprocess the circuit voltage and current across lumped port index using the electric
   // field solution. When the internal grid functions are real-valued, the returned voltage
   // has only a nonzero real part.
-  std::complex<double> GetPortPower(const LumpedPortOperator &lumped_port_op,
-                                    int idx) const;
-  std::complex<double> GetPortPower(const WavePortOperator &wave_port_op, int idx) const;
-  std::complex<double> GetPortVoltage(const LumpedPortOperator &lumped_port_op,
-                                      int idx) const;
-  std::complex<double> GetPortVoltage(const WavePortOperator &wave_port_op, int idx) const;
+  std::complex<double> GetPortPower(int idx) const;
+  std::complex<double> GetPortVoltage(int idx) const;
   std::complex<double>
-  GetPortCurrent(const LumpedPortOperator &lumped_port_op, int idx,
+  GetPortCurrent(int idx,
                  LumpedPortData::Branch branch = LumpedPortData::Branch::TOTAL) const;
-  std::complex<double> GetPortCurrent(const WavePortOperator &wave_port_op, int idx) const;
 
   // Postprocess the EPR for the electric field solution and lumped port index.
-  double GetInductorParticipation(const LumpedPortOperator &lumped_port_op, int idx,
-                                  double E_m) const;
+  double GetInductorParticipation(int idx, double E_m) const;
 
   // Postprocess the coupling rate for radiative loss to the given I-O port index.
-  double GetExternalKappa(const LumpedPortOperator &lumped_port_op, int idx,
-                          double E_m) const;
+  double GetExternalKappa(int idx, double E_m) const;
 
   // Write to disk the E- and B-fields extracted from the solution vectors. Note that fields
   // are not redimensionalized, to do so one needs to compute: B <= B * (μ₀ H₀), E <= E *