From a97349cadff23b236c8ad63ceb1d4a23f15e8971 Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Tue, 3 Oct 2023 10:16:28 -0400 Subject: [PATCH] Fixes to better handle normalization --- palace/drivers/basesolver.cpp | 11 +++++----- palace/drivers/basesolver.hpp | 7 +++--- palace/drivers/drivensolver.cpp | 12 +++++++---- palace/drivers/eigensolver.cpp | 6 ++++-- palace/drivers/electrostaticsolver.cpp | 6 ++++-- palace/drivers/magnetostaticsolver.cpp | 6 ++++-- palace/drivers/transientsolver.cpp | 16 +++++++------- palace/linalg/errorestimator.cpp | 30 ++++++++++++++++++-------- palace/linalg/errorestimator.hpp | 7 +++--- palace/main.cpp | 2 +- 10 files changed, 64 insertions(+), 39 deletions(-) diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp index e9b0c37d0..be62bced6 100644 --- a/palace/drivers/basesolver.cpp +++ b/palace/drivers/basesolver.cpp @@ -567,7 +567,8 @@ void BaseSolver::PostprocessFields(const PostOperator &postop, int step, double } void BaseSolver::PostprocessErrorIndicators(const std::string &name, int step, double time, - const ErrorIndicators &indicators) const + const ErrorIndicators &indicators, + bool normalized) const { if (post_dir.length() == 0) { @@ -584,10 +585,10 @@ void BaseSolver::PostprocessErrorIndicators(const std::string &name, int step, d // clang-format off output.print("{:>{}s},{:>{}s},{:>{}s},{:>{}s},{:>{}s},{:>{}s}\n", name, table.w1, - "Global Rel.", table.w, - "Element Min.", table.w, - "Element Max.", table.w, - "Element Mean", table.w, + normalized ? "Sum Rel." : "Sum Abs.", table.w, + normalized ? "Min Elem Rel." : "Min Elem Abs.", table.w, + normalized ? "Max Elem Rel." : "Max Elem Abs.", table.w, + normalized ? "Mean Elem Rel." : "Mean Elem Abs.", table.w, "Normalization", table.w); // clang-format on } diff --git a/palace/drivers/basesolver.hpp b/palace/drivers/basesolver.hpp index ed415ca89..e411c0f24 100644 --- a/palace/drivers/basesolver.hpp +++ b/palace/drivers/basesolver.hpp @@ -67,10 +67,11 @@ class BaseSolver double time) const; void PostprocessFields(const PostOperator &postop, int step, double time) const; - // Postprocess granular error indicator to file. + // Postprocess granular error indicator to file. The argument normalized indicates if + // supplied indicators have already been normalized. void PostprocessErrorIndicators(const std::string &name, int step, double time, - const ErrorIndicators &indicators) const; - // Append combined error indicator to file. + const ErrorIndicators &indicators, bool normalized) const; + // Write a string labeled error indicator. Used for writing statistics. void PostprocessErrorIndicators(const std::string &name, const ErrorIndicators &indicators) const; diff --git a/palace/drivers/drivensolver.cpp b/palace/drivers/drivensolver.cpp index b5bed2ae9..0c944413c 100644 --- a/palace/drivers/drivensolver.cpp +++ b/palace/drivers/drivensolver.cpp @@ -151,13 +151,15 @@ ErrorIndicators DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop](const auto &E, int step, double f) { BlockTimer bt0(Timer::ESTSOLVE); - auto estimate = estimator(E); + constexpr bool normalized = true; + auto estimate = estimator(E, normalized); BlockTimer bt1(Timer::POSTPRO); // Write the indicator for this mode. postop.SetIndicatorGridFunction(estimate.indicators); PostprocessErrorIndicators( "f (GHz)", step, f, - ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}); + ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}, + normalized); ErrorReducer(indicators, std::move(estimate)); }; @@ -285,12 +287,14 @@ ErrorIndicators DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &ErrorReducer](const auto &E, int step, double frequency) { BlockTimer bt0(Timer::ESTSOLVE); - auto estimate = estimator(E); + constexpr bool normalized = true; + auto estimate = estimator(E, normalized); BlockTimer bt1(Timer::POSTPRO); // Write the indicator for this mode. PostprocessErrorIndicators( "f (GHz)", step, frequency, - ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}); + ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}, + normalized); ErrorReducer(indicators, std::move(estimate)); }; diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index d1923d4bc..6e34bcfac 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -272,13 +272,15 @@ EigenSolver::Solve(const std::vector> &mesh) cons [this, &estimator, &indicators, &ErrorReducer, &postop](const auto &E, int i) { BlockTimer bt0(Timer::ESTSOLVE); - auto ind = estimator(E); + constexpr bool normalized = true; + auto ind = estimator(E, normalized); BlockTimer bt1(Timer::POSTPRO); postop.SetIndicatorGridFunction(ind.indicators); // Write the indicator for this mode. PostprocessErrorIndicators( "m", i, i + 1, - ErrorIndicators{ind, indicators.GlobalTrueVSize(), indicators.GetComm()}); + ErrorIndicators{ind, indicators.GlobalTrueVSize(), indicators.GetComm()}, + normalized); ErrorReducer(indicators, std::move(ind)); }; diff --git a/palace/drivers/electrostaticsolver.cpp b/palace/drivers/electrostaticsolver.cpp index 8a834b490..3f9fcb4c6 100644 --- a/palace/drivers/electrostaticsolver.cpp +++ b/palace/drivers/electrostaticsolver.cpp @@ -108,14 +108,16 @@ ErrorIndicators ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop, &postop](const auto &V, int i, double idx) { BlockTimer bt0(Timer::ESTSOLVE); - auto estimate = estimator(V); + constexpr bool normalized = true; + auto estimate = estimator(V, normalized); ErrorReducer(indicators, estimate); BlockTimer bt1(Timer::POSTPRO); // Write the indicator for this mode. postop.SetIndicatorGridFunction(estimate.indicators); PostprocessErrorIndicators( "i", i, idx, - ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}); + ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}, + normalized); ErrorReducer(indicators, std::move(estimate)); }; if (iodata.solver.electrostatic.n_post > 0) diff --git a/palace/drivers/magnetostaticsolver.cpp b/palace/drivers/magnetostaticsolver.cpp index 1d65f920a..b19389fe1 100644 --- a/palace/drivers/magnetostaticsolver.cpp +++ b/palace/drivers/magnetostaticsolver.cpp @@ -110,13 +110,15 @@ ErrorIndicators MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, &postop](const auto &A, int i, double idx) { BlockTimer bt0(Timer::ESTSOLVE); - auto estimate = estimator(A); + constexpr bool normalized = true; + auto estimate = estimator(A, normalized); BlockTimer bt1(Timer::POSTPRO); // Write the indicator for this mode. postop.SetIndicatorGridFunction(estimate.indicators); PostprocessErrorIndicators( "i", i, idx, - ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}); + ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}, + normalized); ErrorReducer(indicators, std::move(estimate)); }; if (iodata.solver.magnetostatic.n_post > 0) diff --git a/palace/drivers/transientsolver.cpp b/palace/drivers/transientsolver.cpp index ccbb92005..568720de5 100644 --- a/palace/drivers/transientsolver.cpp +++ b/palace/drivers/transientsolver.cpp @@ -90,12 +90,15 @@ TransientSolver::Solve(const std::vector> &mesh) &postop](const auto &E, int step, double time) { BlockTimer bt0(Timer::ESTSOLVE); - auto estimate = estimator(E); + // Initial flux of zero would return nan. + bool constexpr normalized = false; + auto estimate = estimator(E, normalized); BlockTimer bt1(Timer::POSTPRO); postop.SetIndicatorGridFunction(estimate.indicators); PostprocessErrorIndicators( "t (ns)", step, time, - ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}); + ErrorIndicators{estimate, indicators.GlobalTrueVSize(), indicators.GetComm()}, + normalized); ErrorReducer(indicators, std::move(estimate)); }; @@ -140,11 +143,8 @@ TransientSolver::Solve(const std::vector> &mesh) E_elec + E_mag); } - if (step > 0) - { - // Calculate and record the error indicators. - UpdateErrorIndicators(E, step, t); - } + // Calculate and record the error indicators. + UpdateErrorIndicators(E, step, t); // Postprocess port voltages/currents and optionally write solution to disk. Postprocess(postop, spaceop.GetLumpedPortOp(), spaceop.GetSurfaceCurrentOp(), step, t, @@ -154,7 +154,7 @@ TransientSolver::Solve(const std::vector> &mesh) step++; } SaveMetadata(timeop.GetLinearSolver()); - PostprocessErrorIndicators("Mean", step, t, indicators); + PostprocessErrorIndicators("Mean", indicators); return indicators; } std::function TransientSolver::GetTimeExcitation(bool dot) const diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index c98f80a44..f65808136 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -47,7 +47,8 @@ CurlFluxErrorEstimator::CurlFluxErrorEstimator( } } -IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const ComplexVector &v) const +IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const ComplexVector &v, + bool normalize) const { mfem::ParComplexGridFunction field(&fes); field.real().SetFromTrueDofs(v.Real()); @@ -152,12 +153,16 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const ComplexVecto Mpi::GlobalSum(1, &normalization, field.ParFESpace()->GetComm()); normalization = std::sqrt(normalization); - std::for_each(estimates.begin(), estimates.end(), - [&normalization](auto &x) { x /= normalization; }); + if (normalize) + { + std::for_each(estimates.begin(), estimates.end(), + [&normalization](auto &x) { x /= normalization; }); + } return {estimates, normalization}; } -IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const Vector &v) const +IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const Vector &v, + bool normalize) const { mfem::ParGridFunction field(&fes); field.SetFromTrueDofs(v); @@ -236,8 +241,11 @@ IndicatorsAndNormalization CurlFluxErrorEstimator::operator()(const Vector &v) c Mpi::GlobalSum(1, &normalization, field.ParFESpace()->GetComm()); normalization = std::sqrt(normalization); - std::for_each(estimates.begin(), estimates.end(), - [&normalization](auto &x) { x /= normalization; }); + if (normalize) + { + std::for_each(estimates.begin(), estimates.end(), + [&normalization](auto &x) { x /= normalization; }); + } return {estimates, normalization}; } @@ -276,7 +284,8 @@ GradFluxErrorEstimator::GradFluxErrorEstimator( } } -IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v) const +IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v, + bool normalize) const { mfem::ParGridFunction field(&fes); field.SetFromTrueDofs(v); @@ -388,8 +397,11 @@ IndicatorsAndNormalization GradFluxErrorEstimator::operator()(const Vector &v) c paraview.Save(); } - std::for_each(estimates.begin(), estimates.end(), - [&normalization](auto &x) { x /= normalization; }); + if (normalize) + { + std::for_each(estimates.begin(), estimates.end(), + [&normalization](auto &x) { x /= normalization; }); + } return {estimates, normalization}; } } // namespace palace diff --git a/palace/linalg/errorestimator.hpp b/palace/linalg/errorestimator.hpp index 0f7d6286b..8d9b9a672 100644 --- a/palace/linalg/errorestimator.hpp +++ b/palace/linalg/errorestimator.hpp @@ -46,10 +46,11 @@ class CurlFluxErrorEstimator mfem::ParFiniteElementSpace &fes); // Compute elemental error indicators given a complex vector of true DOF. - IndicatorsAndNormalization operator()(const ComplexVector &v) const; + IndicatorsAndNormalization operator()(const ComplexVector &v, + bool normalize = true) const; // Compute elemental error indicators given a vector of true DOF. - IndicatorsAndNormalization operator()(const Vector &v) const; + IndicatorsAndNormalization operator()(const Vector &v, bool normalize = true) const; }; // Class used for computing grad flux error estimate, @@ -82,7 +83,7 @@ class GradFluxErrorEstimator mfem::ParFiniteElementSpace &fes); // Compute elemental error indicators given a vector of true DOF. - IndicatorsAndNormalization operator()(const Vector &v) const; + IndicatorsAndNormalization operator()(const Vector &v, bool normalize = true) const; }; } // namespace palace diff --git a/palace/main.cpp b/palace/main.cpp index f62689334..a19902b58 100644 --- a/palace/main.cpp +++ b/palace/main.cpp @@ -188,7 +188,7 @@ int main(int argc, char *argv[]) // Run the problem driver. auto solver_output = solver->Solve(mesh); - Mpi::Print(world_comm, "Normalized Error Indicator: {:.3e}", + Mpi::Print(world_comm, "Error Estimate: {:.3e}\n", solver_output.GetGlobalErrorIndicator()); // Print timing summary.