From 4d0dc702f1b51fd3bc4b9b0fd4f34f3603916480 Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Mon, 9 Oct 2023 14:15:26 -0400 Subject: [PATCH] New Eval method to avoid duplicate gradient evaluation --- palace/fem/coefficient.hpp | 13 +++++++-- palace/linalg/errorestimator.cpp | 47 +++++++++++++++----------------- palace/linalg/errorestimator.hpp | 3 +- 3 files changed, 34 insertions(+), 29 deletions(-) diff --git a/palace/fem/coefficient.hpp b/palace/fem/coefficient.hpp index 3ff76e661..97bf711c4 100644 --- a/palace/fem/coefficient.hpp +++ b/palace/fem/coefficient.hpp @@ -508,14 +508,14 @@ class GradFluxCoefficient : public mfem::Coefficient private: const mfem::ParGridFunction &gf; const MaterialOperator &mat_op; - mfem::Vector grad, V; + mfem::Vector grad, tmp; int component; public: GradFluxCoefficient(const mfem::ParGridFunction &gf, const MaterialOperator &mat_op) : mfem::Coefficient(), gf(gf), mat_op(mat_op), grad(gf.ParFESpace()->GetParMesh()->SpaceDimension()), - V(gf.ParFESpace()->GetParMesh()->SpaceDimension()), component(-1) + tmp(gf.ParFESpace()->GetParMesh()->SpaceDimension()), component(-1) { } @@ -532,9 +532,16 @@ class GradFluxCoefficient : public mfem::Coefficient MFEM_ASSERT(component >= 0 && component < gf.ParFESpace()->GetParMesh()->SpaceDimension(), "Invalid component index, try calling SetComponent(int)!"); + Eval(tmp, T, ip); + return tmp(component); + } + + // VectorCoefficient style Eval method for use in elemental integral evaluation. + inline void Eval(mfem::Vector &V, mfem::ElementTransformation &T, + const mfem::IntegrationPoint &ip) + { gf.GetGradient(T, grad); mat_op.GetPermittivityReal(T.Attribute).Mult(grad, V); - return V(component); } }; diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index 34d414821..57e07370f 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -278,7 +278,8 @@ GradFluxErrorEstimator::GradFluxErrorEstimator( smooth_flux(h1_fespaces.GetFinestFESpace().GetTrueVSize()), flux_rhs(h1_fespaces.GetFinestFESpace().GetTrueVSize()), field_gf(&h1_fespaces.GetFinestFESpace()), - smooth_flux_gf(&h1_fespaces.GetFinestFESpace()) + smooth_flux_gf{&h1_fespaces.GetFinestFESpace(), &h1_fespaces.GetFinestFESpace(), + &h1_fespaces.GetFinestFESpace()} { } @@ -310,39 +311,35 @@ ErrorIndicator GradFluxErrorEstimator::ComputeIndicators(const Vector &v) const h1_fespace.GetProlongationMatrix()->MultTranspose(rhs, flux_rhs); smooth_projector.Mult(flux_rhs, smooth_flux); } - smooth_flux_gf.SetFromTrueDofs(smooth_flux); - smooth_flux_gf.ExchangeFaceNbrData(); - for (int e = 0; e < h1_fespace.GetNE(); e++) + smooth_flux_gf[c].SetFromTrueDofs(smooth_flux); + smooth_flux_gf[c].ExchangeFaceNbrData(); + } + + Vector coef_eval(3); + for (int e = 0; e < h1_fespace.GetNE(); e++) + { + auto &T = *h1_fespace.GetElementTransformation(e); + // integration order 2p + q + const auto &ir = mfem::IntRules.Get(T.GetGeometryType(), + 2 * h1_fespace.GetFE(e)->GetOrder() + T.Order()); + for (const auto &ip : ir) { - auto &T = *h1_fespace.GetElementTransformation(e); - // integration order 2p + q - const auto &ir = mfem::IntRules.Get(T.GetGeometryType(), - 2 * h1_fespace.GetFE(e)->GetOrder() + T.Order()); - for (const auto &ip : ir) + T.SetIntPoint(&ip); + coef.Eval(coef_eval, T, ip); + for (int c = 0; c < sdim; c++) { - T.SetIntPoint(&ip); - double smooth_val = smooth_flux_gf.GetValue(e, ip); - double coarse_val = coef.Eval(T, ip); + coef.SetComponent(c); + double smooth_val = smooth_flux_gf[c].GetValue(e, ip); const double w_i = ip.weight * T.Weight(); - estimates[e] += w_i * std::pow(smooth_val - coarse_val, 2.0); - normalization += w_i * std::pow(coarse_val, 2.0); + estimates[e] += w_i * std::pow(smooth_val - coef_eval(c), 2.0); + normalization += w_i * std::pow(coef_eval(c), 2.0); } } - if constexpr (false) - { - // Debugging branch generates some intermediate fields for paraview. - mfem::ParaViewDataCollection paraview("debug_coeff" + std::to_string(c), - h1_fespace.GetParMesh()); - paraview.RegisterCoeffField("Flux", &coef); - paraview.RegisterField("SmoothFlux", &smooth_flux_gf); - paraview.Save(); - } + estimates[e] = std::sqrt(estimates[e]); } - linalg::Sqrt(estimates); Mpi::GlobalSum(1, &normalization, h1_fespace.GetComm()); normalization = std::sqrt(normalization); - if (normalization > 0) { estimates /= normalization; diff --git a/palace/linalg/errorestimator.hpp b/palace/linalg/errorestimator.hpp index 872ec9f58..c4f1b7c52 100644 --- a/palace/linalg/errorestimator.hpp +++ b/palace/linalg/errorestimator.hpp @@ -85,7 +85,8 @@ class GradFluxErrorEstimator FluxProjector smooth_projector; mutable Vector smooth_flux, flux_rhs; - mutable mfem::ParGridFunction field_gf, smooth_flux_gf; + mutable mfem::ParGridFunction field_gf; + mutable std::array smooth_flux_gf; public: // Constructor for using geometric and p multigrid.