From 97372d7f72303c0f0be2e19352f84d521de2f812 Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Tue, 3 Oct 2023 11:07:37 -0400 Subject: [PATCH] Cleanup debris and formatting fixes --- palace/drivers/drivensolver.cpp | 14 +++++++++----- palace/drivers/eigensolver.cpp | 1 - palace/drivers/electrostaticsolver.cpp | 2 -- palace/drivers/transientsolver.hpp | 2 +- palace/fem/coefficient.hpp | 2 -- palace/fem/integrator.hpp | 1 - palace/linalg/errorestimator.cpp | 12 ------------ palace/linalg/errorestimator.hpp | 8 ++++---- palace/linalg/fluxprojector.hpp | 14 +++++++------- palace/main.cpp | 2 +- palace/models/postoperator.cpp | 12 ++++++------ palace/models/postoperator.hpp | 3 +-- palace/utils/timer.hpp | 1 + 13 files changed, 30 insertions(+), 44 deletions(-) diff --git a/palace/drivers/drivensolver.cpp b/palace/drivers/drivensolver.cpp index 0c944413c..1155218e8 100644 --- a/palace/drivers/drivensolver.cpp +++ b/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(Operator::DIAG_ONE); auto C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); auto M = spaceop.GetMassMatrix(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. diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index 6e34bcfac..0b085c188 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -251,7 +251,6 @@ EigenSolver::Solve(const std::vector> &mesh) cons ksp->SetOperators(*A, *P); eigen->SetLinearSolver(*ksp); - // TODO: Add a block timer for estimate construction auto estimator = [&]() { BlockTimer bt(Timer::ESTCONSTRUCT); diff --git a/palace/drivers/electrostaticsolver.cpp b/palace/drivers/electrostaticsolver.cpp index 3f9fcb4c6..cdb025b17 100644 --- a/palace/drivers/electrostaticsolver.cpp +++ b/palace/drivers/electrostaticsolver.cpp @@ -70,8 +70,6 @@ ElectrostaticSolver::Solve(const std::vector> &me // Next terminal. step++; } - - // Construct estimator and reducer for error indicators. auto estimator = [&]() { BlockTimer bt(Timer::ESTCONSTRUCT); diff --git a/palace/drivers/transientsolver.hpp b/palace/drivers/transientsolver.hpp index 85dba778d..20a17ac5b 100644 --- a/palace/drivers/transientsolver.hpp +++ b/palace/drivers/transientsolver.hpp @@ -51,7 +51,7 @@ class TransientSolver : public BaseSolver using BaseSolver::BaseSolver; ErrorIndicators - Solve(const std::vector> &mesh) const override; + Solve(const std::vector> &mesh) const final; }; } // namespace palace diff --git a/palace/fem/coefficient.hpp b/palace/fem/coefficient.hpp index 464c9e536..3ad9e5db2 100644 --- a/palace/fem/coefficient.hpp +++ b/palace/fem/coefficient.hpp @@ -484,7 +484,6 @@ class CurlFluxCoefficient : public mfem::VectorCoefficient private: const mfem::ParGridFunction &X; MaterialPropertyCoefficient coef; - mfem::DenseMatrix muinv; mfem::Vector curl; @@ -511,7 +510,6 @@ class GradFluxCoefficient : public mfem::VectorCoefficient private: const mfem::ParGridFunction φ MaterialPropertyCoefficient coef; - mfem::Vector grad; mfem::DenseMatrix eps; diff --git a/palace/fem/integrator.hpp b/palace/fem/integrator.hpp index 367dfa2fb..6e068f32c 100644 --- a/palace/fem/integrator.hpp +++ b/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); diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index f65808136..813b7b363 100644 --- a/palace/linalg/errorestimator.cpp +++ b/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) diff --git a/palace/linalg/errorestimator.hpp b/palace/linalg/errorestimator.hpp index 8d9b9a672..17d948ff9 100644 --- a/palace/linalg/errorestimator.hpp +++ b/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; diff --git a/palace/linalg/fluxprojector.hpp b/palace/linalg/fluxprojector.hpp index eccda9fc0..24194b615 100644 --- a/palace/linalg/fluxprojector.hpp +++ b/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 M; - // Linear solver and preconditioner for the projected linear system M σ = σ̂ + // Linear solver and preconditioner for the projected linear system M σ = σ̂. std::unique_ptr ksp; + // Intermediate storage vector used in Mult. mutable Vector tmp; public: diff --git a/palace/main.cpp b/palace/main.cpp index a19902b58..7be814a39 100644 --- a/palace/main.cpp +++ b/palace/main.cpp @@ -155,7 +155,7 @@ int main(int argc, char *argv[]) #endif // Initialize the problem driver. - const std::unique_ptr solver = [&]() -> std::unique_ptr + const auto solver = [&]() -> std::unique_ptr { switch (iodata.problem.type) { diff --git a/palace/models/postoperator.cpp b/palace/models/postoperator.cpp index a49fc6978..49e7c2f5b 100644 --- a/palace/models/postoperator.cpp +++ b/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) diff --git a/palace/models/postoperator.hpp b/palace/models/postoperator.hpp index 943a1dfa2..ff0e6b898 100644 --- a/palace/models/postoperator.hpp +++ b/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 diff --git a/palace/utils/timer.hpp b/palace/utils/timer.hpp index 15a305164..e55d341d1 100644 --- a/palace/utils/timer.hpp +++ b/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) {