Cleanup debris and formatting fixes

palace/drivers/drivensolver.cpp

@@ -120,6 +120,7 @@ ErrorIndicators DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator
   // Because the Dirichlet BC is always homogenous, no special elimination is required on
   // the RHS. Assemble the linear system for the initial frequency (so we can call
   // KspSolver::SetOperators). Compute everything at the first frequency step.
+  BlockTimer bt0(Timer::CONSTRUCT);
   auto K = spaceop.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
   auto C = spaceop.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
   auto M = spaceop.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
@@ -165,11 +166,12 @@ ErrorIndicators DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator
 
   // Main frequency sweep loop.
   double omega = omega0;
+  auto t0 = Timer::Now();
   for (int step = step0; step < nstep; step++)
   {
-    // const double freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
-    // Mpi::Print("\nIt {:d}/{:d}: ω/2π = {:.3e} GHz (elapsed time = {:.2e} s)\n", step + 1,
-    //            nstep, freq, Timer::Duration(Timer::Now() - t0).count());
+    const double freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
+    Mpi::Print("\nIt {:d}/{:d}: ω/2π = {:.3e} GHz (elapsed time = {:.2e} s)\n", step + 1,
+               nstep, freq, Timer::Duration(Timer::Now() - t0).count());
 
     // Assemble the linear system.
     if (step > step0)
@@ -191,8 +193,7 @@ ErrorIndicators DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator
     ksp.Mult(RHS, E);
 
     // Compute the error indicators, and post process the indicator field.
-    UpdateErrorIndicators(E, step,
-                          iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega));
+    UpdateErrorIndicators(E, step, freq);
 
     // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
     // PostOperator for all postprocessing operations.
@@ -233,6 +234,7 @@ ErrorIndicators DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
                                             double delta_omega) const
 {
   // Configure default parameters if not specified.
+  BlockTimer bt0(Timer::CONSTRUCT);
   double offline_tol = iodata.solver.driven.adaptive_tol;
   int nmax = iodata.solver.driven.adaptive_nmax;
   int ncand = iodata.solver.driven.adaptive_ncand;
@@ -339,6 +341,8 @@ ErrorIndicators DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
              (iter == nmax) ? " reached maximum" : " converged with", iter,
              prom.GetReducedDimension(), max_error, offline_tol);
   utils::PrettyPrint(prom.GetSampleFrequencies(), f0, " Sampled frequencies (GHz):");
+  Mpi::Print(" Total offline phase elapsed time: {:.2e} s\n",
+             Timer::Duration(Timer::Now() - t0).count());  // Timing on root
   SaveMetadata(prom.GetLinearSolver());
 
   // Set the indicator field to the combined field for postprocessing.

palace/drivers/eigensolver.cpp

@@ -251,7 +251,6 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
   ksp->SetOperators(*A, *P);
   eigen->SetLinearSolver(*ksp);
 
-  // TODO: Add a block timer for estimate construction
   auto estimator = [&]()
   {
     BlockTimer bt(Timer::ESTCONSTRUCT);

palace/drivers/electrostaticsolver.cpp

@@ -70,8 +70,6 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &me
     // Next terminal.
     step++;
   }
-
-  // Construct estimator and reducer for error indicators.
   auto estimator = [&]()
   {
     BlockTimer bt(Timer::ESTCONSTRUCT);

palace/drivers/transientsolver.hpp

@@ -51,7 +51,7 @@ class TransientSolver : public BaseSolver
   using BaseSolver::BaseSolver;
 
   ErrorIndicators
-  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const override;
+  Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const final;
 };
 
 }  // namespace palace

palace/fem/coefficient.hpp

@@ -484,7 +484,6 @@ class CurlFluxCoefficient : public mfem::VectorCoefficient
 private:
   const mfem::ParGridFunction &X;
   MaterialPropertyCoefficient<MaterialPropertyType::INV_PERMEABILITY> coef;
-
   mfem::DenseMatrix muinv;
   mfem::Vector curl;
 
@@ -511,7 +510,6 @@ class GradFluxCoefficient : public mfem::VectorCoefficient
 private:
   const mfem::ParGridFunction &phi;
   MaterialPropertyCoefficient<MaterialPropertyType::PERMITTIVITY_REAL> coef;
-
   mfem::Vector grad;
   mfem::DenseMatrix eps;
 

palace/fem/integrator.hpp

@@ -66,7 +66,6 @@ class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator
       fe.CalcVShape(ip, vshape);
 
       Q.Eval(f_loc, Tr, ip);
-
       Tr.InverseJacobian().Mult(f_loc, f_hat);
       f_hat *= ip.weight * Tr.Weight();
       vshape.AddMult(f_hat, elvect);

palace/linalg/errorestimator.cpp

@@ -53,7 +53,6 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const ComplexVecto
   mfem::ParComplexGridFunction field(&fes);
   field.real().SetFromTrueDofs(v.Real());
   field.imag().SetFromTrueDofs(v.Imag());
-
   const int nelem = smooth_flux_fes.GetFinestFESpace().GetNE();
 
   // Coefficients for computing the discontinuous flux., i.e. (W, μ⁻¹∇ × V).
@@ -97,11 +96,8 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const ComplexVecto
     flux.imag().ExchangeFaceNbrData();
     return flux;
   };
-
   auto smooth_flux_func = build_func(smooth_flux, smooth_flux_fes.GetFinestFESpace());
-
   mfem::ParComplexGridFunction coarse_flux_func(&coarse_flux_fes);
-
   coarse_flux_func.real().ProjectCoefficient(real_coef);
   coarse_flux_func.imag().ProjectCoefficient(imag_coef);
 
@@ -166,7 +162,6 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const Vector &v,
 {
   mfem::ParGridFunction field(&fes);
   field.SetFromTrueDofs(v);
-
   const int nelem = smooth_flux_fes.GetFinestFESpace().GetNE();
 
   // Coefficients for computing the discontinuous flux., i.e. (W, μ⁻¹∇ × V).
@@ -183,7 +178,6 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const Vector &v,
 
     return RHS;
   };
-
   const auto smooth_flux_rhs = rhs_from_coef(smooth_flux_fes.GetFinestFESpace(), coef);
 
   // Given the RHS vector of non-smooth flux, construct a flux projector and perform mass
@@ -207,7 +201,6 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const Vector &v,
     return flux;
   };
   auto smooth_flux_func = build_func(smooth_flux, smooth_flux_fes.GetFinestFESpace());
-
   mfem::ParGridFunction coarse_flux_func(&coarse_flux_fes);
   coarse_flux_func.ProjectCoefficient(coef);
 
@@ -289,7 +282,6 @@ IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v,
 {
   mfem::ParGridFunction field(&fes);
   field.SetFromTrueDofs(v);
-
   const int nelem = smooth_flux_fes.GetNE();
 
   // Coefficients for computing the discontinuous flux., i.e. (V, ϵ ∇ ϕ).
@@ -328,7 +320,6 @@ IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v,
       rhs_comp.MakeRef(rhs, i * stride, stride);
       proj.Mult(rhs_comp, flux_comp);
     }
-
     return flux;
   };
   auto smooth_flux = build_flux(smooth_projector, smooth_flux_rhs);
@@ -362,7 +353,6 @@ IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v,
     const int ndof = coarse_vec.Size() / 3;
     coarse_sub_vec.SetSize(ndof);
     smooth_sub_vec.SetSize(ndof);
-
     for (int c = 0; c < 3; c++)
     {
       coarse_sub_vec.MakeRef(coarse_vec, c * ndof);
@@ -373,7 +363,6 @@ IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v,
       estimates[e] += scalar_mass_matrices[e].InnerProduct(coarse_sub_vec,
                                                            coarse_sub_vec);  // Integrate
     }
-
     estimates[e] = std::sqrt(estimates[e]);
   }
   Mpi::GlobalSum(1, &normalization, field.ParFESpace()->GetComm());
@@ -394,7 +383,6 @@ IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v,
     est_field.SetFromTrueDofs(estimates);
 
     paraview.RegisterField("ErrorIndicator", &est_field);
-
     paraview.Save();
   }
   if (normalize)

palace/linalg/errorestimator.hpp

@@ -21,8 +21,8 @@ class PetscParVector;
 }  // namespace petsc
 
 // Class used for computing curl flux error estimate,
-// i.e. || μ⁻¹∇ × V - F ||_K
-// where F denotes a smooth reconstruction of μ⁻¹∇ × V.
+// i.e. || μ⁻¹∇ × Vₕ - F ||_K
+// where F denotes a smooth reconstruction of μ⁻¹∇ × Vₕ.
 class CurlFluxErrorEstimator
 {
   const MaterialOperator &mat_op;
@@ -54,8 +54,8 @@ class CurlFluxErrorEstimator
 };
 
 // Class used for computing grad flux error estimate,
-// i.e. || ϵ ∇ ϕ - F ||_K
-// where F denotes a smooth reconstruction of ϵ ∇ ϕ.
+// i.e. || ϵ ∇ ϕₕ - F ||_K
+// where F denotes a smooth reconstruction of ϵ ∇ ϕₕ.
 class GradFluxErrorEstimator
 {
   const MaterialOperator &mat_op;

palace/linalg/fluxprojector.hpp

@@ -16,20 +16,20 @@ namespace palace
 class MaterialOperator;
 
 //
-// This solver implements a solver to compute a smooth reconstruction of a
-// discontinuous flux. The difference between this resulting smooth flux and the
-// original non-smooth flux provides a localizable error estimate. An instance
-// of FluxProjector can be reused across solutions, thus the construction of the
-// operator is separated from the construction of the flux RHS.
+// This solver computes a smooth reconstruction of a discontinuous flux. The difference
+// between this resulting smooth flux and the original non-smooth flux provides a
+// localizable error estimate. An instance  of FluxProjector can be reused across solutions,
+// thus the construction of the operator is separated from the construction of the flux RHS.
 class FluxProjector
 {
 private:
-  // Operator for the mass matrix inversion
+  // Operator for the mass matrix inversion.
   std::unique_ptr<Operator> M;
 
-  // Linear solver and preconditioner for the projected linear system M σ = σ̂
+  // Linear solver and preconditioner for the projected linear system M σ = σ̂.
   std::unique_ptr<KspSolver> ksp;
 
+  // Intermediate storage vector used in Mult.
   mutable Vector tmp;
 
 public:

palace/main.cpp

@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
 #endif
 
   // Initialize the problem driver.
-  const std::unique_ptr<BaseSolver> solver = [&]() -> std::unique_ptr<BaseSolver>
+  const auto solver = [&]() -> std::unique_ptr<BaseSolver>
   {
     switch (iodata.problem.type)
     {

palace/models/postoperator.cpp

@@ -47,7 +47,7 @@ PostOperator::PostOperator(const IoData &iodata, SpaceOperator &spaceop,
     E(&spaceop.GetNDSpace()), B(&spaceop.GetRTSpace()), V(std::nullopt), A(std::nullopt),
     indicator_fec(0, spaceop.GetNDSpace().GetParMesh()->Dimension()),
     indicator_fes(spaceop.GetNDSpace().GetParMesh(), &indicator_fec),
-    indicator_field(std::nullopt), lumped_port_init(false), wave_port_init(false),
+    indicator_field(&indicator_fes), lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), spaceop.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  spaceop.GetNDSpace().GetParMesh()),
@@ -104,7 +104,7 @@ PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplaceop,
     V(&laplaceop.GetH1Space()), A(std::nullopt),
     indicator_fec(0, laplaceop.GetH1Space().GetParMesh()->Dimension()),
     indicator_fes(laplaceop.GetH1Space().GetParMesh(), &indicator_fec),
-    indicator_field(std::nullopt), lumped_port_init(false), wave_port_init(false),
+    indicator_field(&indicator_fes), lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), laplaceop.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  laplaceop.GetNDSpace().GetParMesh()),
@@ -134,7 +134,7 @@ PostOperator::PostOperator(const IoData &iodata, CurlCurlOperator &curlcurlop,
     A(&curlcurlop.GetNDSpace()),
     indicator_fec(0, curlcurlop.GetNDSpace().GetParMesh()->Dimension()),
     indicator_fes(curlcurlop.GetNDSpace().GetParMesh(), &indicator_fec),
-    indicator_field(std::nullopt), lumped_port_init(false), wave_port_init(false),
+    indicator_field(&indicator_fes), lumped_port_init(false), wave_port_init(false),
     paraview(CreateParaviewPath(iodata, name), curlcurlop.GetNDSpace().GetParMesh()),
     paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
                  curlcurlop.GetNDSpace().GetParMesh()),
@@ -231,6 +231,7 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
     paraview.RegisterField("A", &*A);
     paraview_bdr.RegisterVCoeffField("A", As.get());
   }
+  paraview.RegisterField("ErrorIndicator", &*indicator_field);
 
   // Extract surface charge from normally discontinuous ND E-field. Also extract surface
   // currents from tangentially discontinuous RT B-field The surface charge and surface
@@ -340,10 +341,9 @@ void PostOperator::SetAGridFunction(const Vector &a)
 
 void PostOperator::SetIndicatorGridFunction(const mfem::Vector &i)
 {
-  indicator_field = mfem::ParGridFunction(&indicator_fes);
+  MFEM_VERIFY(indicator_field,
+              "Incorrect usage of PostOperator::SetIndicatorGridFunction!");
   indicator_field->SetFromTrueDofs(i);
-  // Reregistration overwrites the underlying.
-  paraview.RegisterField("ErrorIndicator", std::addressof(indicator_field.value()));
 }
 
 void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega)

palace/models/postoperator.hpp

@@ -168,9 +168,8 @@ class PostOperator
 
   // 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 *
-  // (Z₀ H₀), V <= V * (Z₀ H₀ L₀), etc. Optionally also write error indicator field.
+  // (Z₀ H₀), V <= V * (Z₀ H₀ L₀), etc.
   void WriteFields(int step, double time) const;
-  void WriteFields(int step, double time, ErrorIndicators &indicators) const;
 
   // Probe the E- and B-fields for their vector-values at speceified locations in space.
   // Locations of probes are set up in constructor from configuration file data. If

palace/utils/timer.hpp

@@ -81,6 +81,7 @@ class Timer
 
   // Return the time elapsed since timer creation.
   Duration TimeFromStart() const { return Now() - start_time; }
+
   // Save a timing step by adding a duration, without lapping; optionally, count it.
   Duration SaveTime(Index idx, Duration time, bool count_it = true)
   {