From 6fe4ca5e5519e30fb932f8e752818a5e7b552e08 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 14 Aug 2023 18:04:35 -0700 Subject: [PATCH] Address PR comments: Clarify comments, fix some typos, rename HCurlNormSolver -> EnergyNormSolver --- palace/fem/multigrid.hpp | 32 ++++++++++++++++---------------- palace/linalg/hcurl.cpp | 13 +++++++------ palace/linalg/hcurl.hpp | 14 +++++++------- palace/models/romoperator.cpp | 2 +- palace/models/romoperator.hpp | 2 +- palace/utils/configfile.hpp | 2 +- palace/utils/iodata.cpp | 1 - 7 files changed, 33 insertions(+), 33 deletions(-) diff --git a/palace/fem/multigrid.hpp b/palace/fem/multigrid.hpp index 3cb3b942b..188f9164d 100644 --- a/palace/fem/multigrid.hpp +++ b/palace/fem/multigrid.hpp @@ -33,7 +33,7 @@ inline auto GetMaxElementOrder(mfem::MixedBilinearForm &a) a.TrialFESpace()->GetMaxElementOrder()); } -// Assembly a bilinear or mixed bilinear form. If the order is lower than the specified +// Assemble a bilinear or mixed bilinear form. If the order is lower than the specified // threshold, the operator is assembled as a sparse matrix. template inline std::unique_ptr @@ -77,16 +77,16 @@ std::vector> inline ConstructFECollections( { // If the solver will use a LOR preconditioner, we need to construct with a specific basis // type. - MFEM_VERIFY(p >= 1, "FE space order must not be less than 1!"); + constexpr int pmin = (std::is_base_of::value || + std::is_base_of::value) + ? 1 + : 0; + MFEM_VERIFY(p >= pmin, "FE space order must not be less than " << pmin << "!"); int b1 = mfem::BasisType::GaussLobatto, b2 = mfem::BasisType::GaussLegendre; if (mat_lor) { b2 = mfem::BasisType::IntegratedGLL; } - constexpr int pm1 = (std::is_base_of::value || - std::is_base_of::value) - ? 0 - : 1; // Construct the p-multigrid hierarchy, first finest to coarsest and then reverse the // order. @@ -96,14 +96,14 @@ std::vector> inline ConstructFECollections( if constexpr (std::is_base_of::value || std::is_base_of::value) { - fecs.push_back(std::make_unique(p - pm1, dim, b1, b2)); + fecs.push_back(std::make_unique(p, dim, b1, b2)); } else { - fecs.push_back(std::make_unique(p - pm1, dim, b1)); + fecs.push_back(std::make_unique(p, dim, b1)); MFEM_CONTRACT_VAR(b2); } - if (p == 1) + if (p == pmin) { break; } @@ -113,7 +113,7 @@ std::vector> inline ConstructFECollections( p--; break; case config::LinearSolverData::MultigridCoarsenType::LOGARITHMIC: - p = (p + 1) / 2; + p = (p + pmin) / 2; break; case config::LinearSolverData::MultigridCoarsenType::INVALID: MFEM_ABORT("Invalid coarsening type for p-multigrid levels!"); @@ -138,18 +138,18 @@ inline mfem::ParFiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy MFEM_VERIFY(!mesh.empty() && !fecs.empty() && (!dbc_tdof_lists || dbc_tdof_lists->empty()), "Empty mesh or FE collection for FE space construction!"); - auto mesh_levels = std::min(mesh.size() - 1, mg_max_levels - fecs.size()); - auto *fespace = new mfem::ParFiniteElementSpace(mesh[mesh.size() - mesh_levels - 1].get(), - fecs[0].get()); + int coarse_mesh_l = + std::max(0, static_cast(mesh.size() + fecs.size()) - 1 - mg_max_levels); + auto *fespace = new mfem::ParFiniteElementSpace(mesh[coarse_mesh_l].get(), fecs[0].get()); if (dbc_marker && dbc_tdof_lists) { fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back()); } - mfem::ParFiniteElementSpaceHierarchy fespaces(mesh[mesh.size() - mesh_levels - 1].get(), - fespace, false, true); + mfem::ParFiniteElementSpaceHierarchy fespaces(mesh[coarse_mesh_l].get(), fespace, false, + true); // h-refinement - for (std::size_t l = mesh.size() - mesh_levels; l < mesh.size(); l++) + for (std::size_t l = coarse_mesh_l + 1; l < mesh.size(); l++) { fespace = new mfem::ParFiniteElementSpace(mesh[l].get(), fecs[0].get()); if (dbc_marker && dbc_tdof_lists) diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp index ae403f3e6..4ae5c7c25 100644 --- a/palace/linalg/hcurl.cpp +++ b/palace/linalg/hcurl.cpp @@ -15,12 +15,13 @@ namespace palace { -HCurlNormSolver::HCurlNormSolver(const MaterialOperator &mat_op, - mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, - mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, - const std::vector> &nd_dbc_tdof_lists, - const std::vector> &h1_dbc_tdof_lists, - double tol, int max_it, int print, int pa_order_threshold) +EnergyNormSolver::EnergyNormSolver(const MaterialOperator &mat_op, + mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, + mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, + const std::vector> &nd_dbc_tdof_lists, + const std::vector> &h1_dbc_tdof_lists, + double tol, int max_it, int print, + int pa_order_threshold) { constexpr int skip_zeros = 0; constexpr auto MatTypeMuInv = MaterialPropertyType::INV_PERMEABILITY; diff --git a/palace/linalg/hcurl.hpp b/palace/linalg/hcurl.hpp index 8f05a2c3a..43760089c 100644 --- a/palace/linalg/hcurl.hpp +++ b/palace/linalg/hcurl.hpp @@ -27,7 +27,7 @@ class MaterialOperator; // // This solver implements a solver for the operator K + M in a Nedelec space. // -class HCurlNormSolver +class EnergyNormSolver { private: // H(curl) norm operator A = K + M and its projection Gáµ€ A G. @@ -37,12 +37,12 @@ class HCurlNormSolver std::unique_ptr ksp; public: - HCurlNormSolver(const MaterialOperator &mat_op, - mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, - mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, - const std::vector> &nd_dbc_tdof_lists, - const std::vector> &h1_dbc_tdof_lists, double tol, - int max_it, int print, int pa_order_threshold); + EnergyNormSolver(const MaterialOperator &mat_op, + mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, + mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, + const std::vector> &nd_dbc_tdof_lists, + const std::vector> &h1_dbc_tdof_lists, double tol, + int max_it, int print, int pa_order_threshold); const Operator &GetOperator() { return *A; } diff --git a/palace/models/romoperator.cpp b/palace/models/romoperator.cpp index dbff4b392..8f59b08b1 100644 --- a/palace/models/romoperator.cpp +++ b/palace/models/romoperator.cpp @@ -114,7 +114,7 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &spaceop) : spaceop if (iodata.solver.driven.adaptive_metric_aposteriori) { constexpr int curlcurl_verbose = 0; - kspKM = std::make_unique( + kspKM = std::make_unique( spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), spaceop.GetH1Spaces(), spaceop.GetNDDbcTDofLists(), spaceop.GetH1DbcTDofLists(), iodata.solver.linear.tol, iodata.solver.linear.max_it, curlcurl_verbose, iodata.solver.pa_order_threshold); diff --git a/palace/models/romoperator.hpp b/palace/models/romoperator.hpp index ffd3a1013..38a4eb09c 100644 --- a/palace/models/romoperator.hpp +++ b/palace/models/romoperator.hpp @@ -42,7 +42,7 @@ class RomOperator std::unique_ptr ksp; // Linear solver for inner product solves for error metric. - std::unique_ptr kspKM; + std::unique_ptr kspKM; // PROM matrices and vectors. Eigen::MatrixXcd Kr, Mr, Cr, Ar; diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 7bd139277..822b6083c 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -740,7 +740,7 @@ struct LinearSolverData // Maximum number of levels for geometric multigrid (set to 1 to disable multigrid). int mg_max_levels = 100; - // Type of coarsening for geometric multigrid. + // Type of coarsening for p-multigrid. enum class MultigridCoarsenType { LINEAR, diff --git a/palace/utils/iodata.cpp b/palace/utils/iodata.cpp index 2c155c20e..e1dd168da 100644 --- a/palace/utils/iodata.cpp +++ b/palace/utils/iodata.cpp @@ -330,7 +330,6 @@ void IoData::CheckConfiguration() // Resolve default values in configuration file. // XX TODO: Default value for pa_order_threshold if we want PA enabled by default - // XX TODO: Also default value for mg_legacy_transfer, maybe always true for PA? if (solver.linear.max_size < 0) {