Address PR feedback

palace/drivers/drivensolver.cpp

@@ -182,8 +182,8 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(spaceop.GetComm(), E),
                linalg::Norml2(spaceop.GetComm(), RHS));
@@ -273,8 +273,8 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
     B *= -1.0 / (1i * omega);
     postop.SetEGridFunction(E, false);
     postop.SetBGridFunction(B, false);
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     estimator.AddErrorIndicator(E, B, E_elec + E_mag, indicator);
   };
   promop.SolveHDM(omega0, E);
@@ -362,8 +362,8 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(spaceop.GetComm(), E));
     {
       const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
@@ -410,8 +410,8 @@ void DrivenSolver::Postprocess(const PostOperator &postop,
   // The internal GridFunctions for PostOperator have already been set from the E and B
   // solutions in the main frequency sweep loop.
   const double freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
-  double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
-  double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
+  const double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
+  const double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessCurrents(postop, surf_j_op, step, omega);
   PostprocessPorts(postop, lumped_port_op, step, omega);
   if (surf_j_op.Size() == 0)

palace/drivers/eigensolver.cpp

@@ -315,8 +315,8 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), omega.real());
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
 
     // Calculate and record the error indicators.
     if (i < iodata.solver.eigenmode.n)
@@ -340,8 +340,8 @@ void EigenSolver::Postprocess(const PostOperator &postop,
 {
   // The internal GridFunctions for PostOperator have already been set from the E and B
   // solutions in the main loop over converged eigenvalues.
-  double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
-  double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
+  const double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
+  const double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessEigen(i, omega, error_bkwd, error_abs, num_conv);
   PostprocessPorts(postop, lumped_port_op, i);
   PostprocessEPR(postop, lumped_port_op, i, omega, E_elec + E_cap);

palace/drivers/electrostaticsolver.cpp

@@ -75,7 +75,7 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Grad.AddMult(V[step], E, -1.0);
     postop.SetVGridFunction(V[step]);
     postop.SetEGridFunction(E);
-    double E_elec = postop.GetEFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
     Mpi::Print(" Sol. ||V|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(laplaceop.GetComm(), V[step]),
                linalg::Norml2(laplaceop.GetComm(), RHS));

palace/drivers/magnetostaticsolver.cpp

@@ -77,7 +77,7 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Curl.Mult(A[step], B);
     postop.SetAGridFunction(A[step]);
     postop.SetBGridFunction(B);
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||A|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(curlcurlop.GetComm(), A[step]),
                linalg::Norml2(curlcurlop.GetComm(), RHS));

palace/drivers/transientsolver.cpp

@@ -115,8 +115,8 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp());
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e}, ||B|| = {:.6e}\n",
                linalg::Norml2(spaceop.GetComm(), E), linalg::Norml2(spaceop.GetComm(), B));
     {
@@ -259,8 +259,8 @@ void TransientSolver::Postprocess(const PostOperator &postop,
   // The internal GridFunctions for PostOperator have already been set from the E and B
   // solutions in the main time integration loop.
   const double ts = iodata.DimensionalizeValue(IoData::ValueType::TIME, t);
-  double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
-  double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
+  const double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
+  const double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessCurrents(postop, surf_j_op, step, t, J_coef);
   PostprocessPorts(postop, lumped_port_op, step, t, J_coef);
   PostprocessDomains(postop, "t (ns)", step, ts, E_elec, E_mag, E_cap, E_ind);

palace/linalg/densematrix.cpp

@@ -267,7 +267,7 @@ mfem::DenseTensor Mult(const mfem::DenseTensor &A, const mfem::DenseTensor &B)
 {
   MFEM_VERIFY(A.SizeK() == B.SizeK(),
               "Size mismatch for product of two DenseTensor objects!");
-  mfem::DenseTensor C(A);
+  mfem::DenseTensor C(A.SizeI(), B.SizeJ(), A.SizeK());
   for (int k = 0; k < C.SizeK(); k++)
   {
     Mult(A(k), B(k), C(k));

palace/linalg/errorestimator.hpp

@@ -55,7 +55,7 @@ class FluxProjector
   void Mult(const VecType &x, VecType &y) const;
 };
 
-// Class used for computing gradient flux error estimate, i.e. || ε Eₕ - D ||_K, where D
+// Class used for computing gradient flux error estimate, η_K = || ε Eₕ - D ||_K, where D
 // denotes a smooth reconstruction of ε Eₕ = ε ∇Vₕ with continuous normal component.
 template <typename VecType = Vector>
 class GradFluxErrorEstimator
@@ -86,7 +86,7 @@ class GradFluxErrorEstimator
   void AddErrorIndicator(const VecType &E, double Et, ErrorIndicator &indicator) const;
 };
 
-// Class used for computing curl flux error estimate, i.e. || μ⁻¹ Bₕ - H ||_K where H
+// Class used for computing curl flux error estimate, η_K = || μ⁻¹ Bₕ - H ||_K where H
 // denotes a smooth reconstruction of μ⁻¹ Bₕ = μ⁻¹ ∇ × Eₕ with continuous tangential
 // component.
 template <typename VecType = Vector>
@@ -118,8 +118,10 @@ class CurlFluxErrorEstimator
   void AddErrorIndicator(const VecType &B, double Et, ErrorIndicator &indicator) const;
 };
 
-// Class used for computing gradient flux error estimate, i.e. || ε Eₕ - D ||_K, where D
-// denotes a smooth reconstruction of ε Eₕ = ε ∇Vₕ with continuous normal component.
+// Class used for computing sum of the gradient flux and curl flux error estimates,
+// η²_K = || ε Eₕ - D ||²_K + || μ⁻¹ Bₕ - H ||²_K, where D and H denote a smooth
+// reconstructions of ε Eₕ = ε ∇Vₕ with continuous normal component and μ⁻¹ Bₕ = μ⁻¹ ∇ × Eₕ
+// with continuous tangential component.
 template <typename VecType>
 class TimeDependentFluxErrorEstimator
 {