From e7e97f9efbff61cf4ee6873e6f43ee083bfe8e23 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 13 Oct 2023 12:40:40 -0700 Subject: [PATCH 1/9] Set default smoother order for multigrid based on the solver order, and increase smooth order for mass-matrix solve in error estimation --- docs/src/config/solver.md | 6 ++++-- palace/linalg/errorestimator.cpp | 22 +++++++++++++--------- palace/linalg/hcurl.cpp | 4 +++- palace/utils/configfile.hpp | 2 +- palace/utils/iodata.cpp | 4 ++++ 5 files changed, 25 insertions(+), 13 deletions(-) diff --git a/docs/src/config/solver.md b/docs/src/config/solver.md index cc3e5357d..620b185d0 100644 --- a/docs/src/config/solver.md +++ b/docs/src/config/solver.md @@ -413,8 +413,10 @@ multigrid preconditioners (when `"UseMultigrid"` is `true` or `"Type"` is `"AMS" `"MGSmoothIts" [1]` : Number of pre- and post-smooth iterations used for multigrid preconditioners (when `"UseMultigrid"` is `true` or `"Type"` is `"AMS"` or `"BoomerAMG"`). -`"MGSmoothOrder" [4]` : Order of polynomial smoothing for geometric multigrid -preconditioning (when `"UseMultigrid"` is `true`). +`"MGSmoothOrder" [0]` : Order of polynomial smoothing for geometric multigrid +preconditioning (when `"UseMultigrid"` is `true`). A value less than 1 defaults to twice +the solution order given in [`config["Solver"]["Order"]`] +(problem.md#config%5B%22Solver%22%5D) or 4, whichever is larger. `"PCMatReal" [false]` : When set to `true`, constructs the preconditioner for frequency domain problems using a real-valued approximation of the system matrix. This is always diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index b7774fbf8..56008e4d8 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -63,9 +63,15 @@ ConfigureLinearSolver(const mfem::ParFiniteElementSpaceHierarchy &fespaces, doub std::unique_ptr> pc; if (fespaces.GetNumLevels() > 1) { + const int dim = fespaces.GetFinestFESpace().GetParMesh()->Dimension(); + const auto type = fespaces.GetFinestFESpace().FEColl()->GetRangeType(dim); + const int mg_smooth_order = + (type == mfem::FiniteElement::SCALAR) + ? 2 + : std::max(fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( - std::move(amg), fespaces, nullptr, 1, 1, 2, 1.0, 0.0, true, pa_order_threshold, - pa_discrete_interp); + std::move(amg), fespaces, nullptr, 1, 1, mg_smooth_order, 1.0, 0.0, true, + pa_order_threshold, pa_discrete_interp); } else { @@ -90,14 +96,13 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, bool pa_discrete_interp) { BlockTimer bt(Timer::CONSTRUCTESTIMATOR); - constexpr int skip_zeros = 0; + constexpr int skip_zeros = false; { constexpr auto MatType = MaterialPropertyType::INV_PERMEABILITY; MaterialPropertyCoefficient muinv_func(mat_op); BilinearForm flux(nd_fespaces.GetFinestFESpace()); flux.AddDomainIntegrator(muinv_func); - Flux = std::make_unique(flux.Assemble(pa_order_threshold, skip_zeros), - nd_fespaces.GetFinestFESpace()); + Flux = std::make_unique(flux.Assemble(), nd_fespaces.GetFinestFESpace()); } M = GetMassMatrix(nd_fespaces, pa_order_threshold, skip_zeros); @@ -115,15 +120,14 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, bool pa_discrete_interp) { BlockTimer bt(Timer::CONSTRUCTESTIMATOR); - constexpr int skip_zeros = 0; + constexpr int skip_zeros = false; { constexpr auto MatType = MaterialPropertyType::PERMITTIVITY_REAL; MaterialPropertyCoefficient epsilon_func(mat_op); BilinearForm flux(h1_fespaces.GetFinestFESpace(), h1d_fespace); flux.AddDomainIntegrator(epsilon_func); - Flux = - std::make_unique(flux.Assemble(pa_order_threshold, skip_zeros), - h1_fespaces.GetFinestFESpace(), h1d_fespace, false); + Flux = std::make_unique(flux.Assemble(), h1_fespaces.GetFinestFESpace(), + h1d_fespace, false); } M = GetMassMatrix(h1_fespaces, pa_order_threshold, skip_zeros); diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp index 792f0cf37..8f6b89aa4 100644 --- a/palace/linalg/hcurl.cpp +++ b/palace/linalg/hcurl.cpp @@ -71,8 +71,10 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( std::unique_ptr> pc; if (nd_fespaces.GetNumLevels() > 1) { + const int mg_smooth_order = + std::max(nd_fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( - std::move(ams), nd_fespaces, &h1_fespaces, 1, 1, 2, 1.0, 0.0, true, + std::move(ams), nd_fespaces, &h1_fespaces, 1, 1, mg_smooth_order, 1.0, 0.0, true, pa_order_threshold, pa_discrete_interp); } else diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index b77fffa9f..a7d8a4e55 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -746,7 +746,7 @@ struct LinearSolverData int mg_smooth_it = 1; // Order of polynomial smoothing for geometric multigrid. - int mg_smooth_order = 4; + int mg_smooth_order = -1; // Safety factors for eigenvalue estimates associated with Chebyshev smoothing for // geometric multigrid. diff --git a/palace/utils/iodata.cpp b/palace/utils/iodata.cpp index a4805e10d..6839fe54d 100644 --- a/palace/utils/iodata.cpp +++ b/palace/utils/iodata.cpp @@ -420,6 +420,10 @@ void IoData::CheckConfiguration() solver.linear.mg_smooth_aux = 1; } } + if (solver.linear.mg_smooth_order < 0) + { + solver.linear.mg_smooth_order = std::max(2 * solver.order, 4); + } } namespace From d670273ab7f9b23964977f45cfa0464381fe7fac Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 27 Oct 2023 09:04:10 -0700 Subject: [PATCH 2/9] Address PR feedback: Fix type --- palace/linalg/errorestimator.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index 56008e4d8..030a55b67 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -96,7 +96,7 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, bool pa_discrete_interp) { BlockTimer bt(Timer::CONSTRUCTESTIMATOR); - constexpr int skip_zeros = false; + constexpr bool skip_zeros = false; { constexpr auto MatType = MaterialPropertyType::INV_PERMEABILITY; MaterialPropertyCoefficient muinv_func(mat_op); @@ -120,7 +120,7 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, bool pa_discrete_interp) { BlockTimer bt(Timer::CONSTRUCTESTIMATOR); - constexpr int skip_zeros = false; + constexpr bool skip_zeros = false; { constexpr auto MatType = MaterialPropertyType::PERMITTIVITY_REAL; MaterialPropertyCoefficient epsilon_func(mat_op); From b1754fda2f04a321c5945c4b9d14899456c8ea70 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 27 Oct 2023 09:12:04 -0700 Subject: [PATCH 3/9] Address PR feedback: Increase polynomial smoothing order regardless of H1 vs. ND for consistency even if performance is similar --- palace/linalg/divfree.cpp | 6 ++++-- palace/linalg/errorestimator.cpp | 6 +----- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/palace/linalg/divfree.cpp b/palace/linalg/divfree.cpp index 44e940554..41dbda7a6 100644 --- a/palace/linalg/divfree.cpp +++ b/palace/linalg/divfree.cpp @@ -65,9 +65,11 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, std::unique_ptr> pc; if (h1_fespaces.GetNumLevels() > 1) { + const int mg_smooth_order = + std::max(h1_fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( - std::move(amg), h1_fespaces, nullptr, 1, 1, 2, 1.0, 0.0, true, pa_order_threshold, - pa_discrete_interp); + std::move(amg), h1_fespaces, nullptr, 1, 1, mg_smooth_order, 1.0, 0.0, true, + pa_order_threshold, pa_discrete_interp); } else { diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index 030a55b67..60b5e3034 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -63,12 +63,8 @@ ConfigureLinearSolver(const mfem::ParFiniteElementSpaceHierarchy &fespaces, doub std::unique_ptr> pc; if (fespaces.GetNumLevels() > 1) { - const int dim = fespaces.GetFinestFESpace().GetParMesh()->Dimension(); - const auto type = fespaces.GetFinestFESpace().FEColl()->GetRangeType(dim); const int mg_smooth_order = - (type == mfem::FiniteElement::SCALAR) - ? 2 - : std::max(fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); + std::max(fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( std::move(amg), fespaces, nullptr, 1, 1, mg_smooth_order, 1.0, 0.0, true, pa_order_threshold, pa_discrete_interp); From f6cb4e7c6d02c8223ffbdb38a4c18e6c4631d1a6 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 23 Oct 2023 12:25:09 -0700 Subject: [PATCH 4/9] Refactor construction of multigrid prolongation operators and discrete linear operators for discrete gradient, curl, etc. to be lazy and allow for reuse between different use cases --- docs/src/config/solver.md | 1 - palace/drivers/basesolver.cpp | 7 +- palace/drivers/basesolver.hpp | 4 +- palace/drivers/drivensolver.cpp | 20 +-- palace/drivers/eigensolver.cpp | 17 +-- palace/drivers/electrostaticsolver.cpp | 16 +-- palace/drivers/magnetostaticsolver.cpp | 10 +- palace/drivers/transientsolver.cpp | 3 +- palace/fem/fespace.cpp | 89 ++++++++++++++ palace/fem/fespace.hpp | 164 +++++++++++++++++++++++++ palace/fem/integrator.hpp | 3 + palace/fem/multigrid.hpp | 114 +++++++++++------ palace/linalg/ams.cpp | 26 ++-- palace/linalg/ams.hpp | 19 +-- palace/linalg/distrelaxation.cpp | 19 +-- palace/linalg/distrelaxation.hpp | 15 +-- palace/linalg/divfree.cpp | 23 ++-- palace/linalg/divfree.hpp | 13 +- palace/linalg/errorestimator.cpp | 61 +++++---- palace/linalg/errorestimator.hpp | 27 ++-- palace/linalg/gmg.cpp | 40 +++--- palace/linalg/gmg.hpp | 29 +++-- palace/linalg/hcurl.cpp | 36 +++--- palace/linalg/hcurl.hpp | 10 +- palace/linalg/ksp.cpp | 26 ++-- palace/linalg/ksp.hpp | 13 +- palace/linalg/operator.hpp | 10 +- palace/linalg/rap.cpp | 17 ++- palace/linalg/rap.hpp | 8 +- palace/models/curlcurloperator.cpp | 24 +--- palace/models/curlcurloperator.hpp | 26 ++-- palace/models/laplaceoperator.cpp | 19 +-- palace/models/laplaceoperator.hpp | 18 +-- palace/models/romoperator.cpp | 3 +- palace/models/spaceoperator.cpp | 138 +++++---------------- palace/models/spaceoperator.hpp | 39 +++--- palace/models/timeoperator.cpp | 2 +- palace/models/timeoperator.hpp | 4 +- palace/utils/configfile.cpp | 3 - palace/utils/configfile.hpp | 4 - scripts/schema/config/solver.json | 1 - test/unit/test-libceed.cpp | 23 ++-- 42 files changed, 661 insertions(+), 483 deletions(-) diff --git a/docs/src/config/solver.md b/docs/src/config/solver.md index 620b185d0..686e611d6 100644 --- a/docs/src/config/solver.md +++ b/docs/src/config/solver.md @@ -456,7 +456,6 @@ vectors in Krylov subspace methods or other parts of the code. ### Advanced linear solver options - `"InitialGuess" [true]` - - `"MGLegacyTransfer" [false]` - `"MGAuxiliarySmoother" [true]` - `"MGSmoothEigScaleMax" [1.0]` - `"MGSmoothEigScaleMin" [0.0]` diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp index 8c66ee91c..8c3e8770f 100644 --- a/palace/drivers/basesolver.cpp +++ b/palace/drivers/basesolver.cpp @@ -8,6 +8,7 @@ #include #include #include "fem/errorindicator.hpp" +#include "fem/fespace.hpp" #include "linalg/ksp.hpp" #include "models/domainpostoperator.hpp" #include "models/postoperator.hpp" @@ -84,17 +85,17 @@ BaseSolver::BaseSolver(const IoData &iodata, bool root, int size, int num_thread } } -void BaseSolver::SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const +void BaseSolver::SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const { if (post_dir.length() == 0) { return; } - const mfem::ParFiniteElementSpace &fespace = fespaces.GetFinestFESpace(); + const auto &fespace = fespaces.GetFinestFESpace(); HYPRE_BigInt ne = fespace.GetParMesh()->GetNE(); Mpi::GlobalSum(1, &ne, fespace.GetComm()); std::vector ndofs(fespaces.GetNumLevels()); - for (int l = 0; l < fespaces.GetNumLevels(); l++) + for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++) { ndofs[l] = fespaces.GetFESpaceAtLevel(l).GlobalTrueVSize(); } diff --git a/palace/drivers/basesolver.hpp b/palace/drivers/basesolver.hpp index 786c99ed4..f8068dbb3 100644 --- a/palace/drivers/basesolver.hpp +++ b/palace/drivers/basesolver.hpp @@ -12,7 +12,6 @@ namespace mfem { -class ParFiniteElementSpaceHierarchy; class ParMesh; } // namespace mfem @@ -21,6 +20,7 @@ namespace palace { class ErrorIndicator; +class FiniteElementSpaceHierarchy; class IoData; class PostOperator; class Timer; @@ -90,7 +90,7 @@ class BaseSolver Solve(const std::vector> &mesh) const = 0; // These methods write different simulation metadata to a JSON file in post_dir. - void SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const; + void SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const; template void SaveMetadata(const SolverType &ksp) const; void SaveMetadata(const Timer &timer) const; diff --git a/palace/drivers/drivensolver.cpp b/palace/drivers/drivensolver.cpp index a6bbfc6f8..169fe19be 100644 --- a/palace/drivers/drivensolver.cpp +++ b/palace/drivers/drivensolver.cpp @@ -116,7 +116,7 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator & auto C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); auto M = spaceop.GetMassMatrix(Operator::DIAG_ZERO); auto A2 = spaceop.GetExtraSystemMatrix(omega0, Operator::DIAG_ZERO); - auto Curl = spaceop.GetCurlMatrix(); + const auto &Curl = spaceop.GetCurlMatrix(); // Set up the linear solver and set operators for the first frequency step. The // preconditioner for the complex linear system is constructed from a real approximation @@ -132,15 +132,14 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator & // Set up RHS vector for the incident field at port boundaries, and the vector for the // first frequency step. - ComplexVector RHS(Curl->Width()), E(Curl->Width()), B(Curl->Height()); + ComplexVector RHS(Curl.Width()), E(Curl.Width()), B(Curl.Height()); E = 0.0; B = 0.0; // Initialize structures for storing and reducing the results of error estimation. CurlFluxErrorEstimator estimator( spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol, - iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold, - iodata.solver.pa_discrete_interp); + iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold); ErrorIndicator indicator; // Main frequency sweep loop. @@ -176,7 +175,8 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator & // PostOperator for all postprocessing operations. BlockTimer bt2(Timer::POSTPRO); double E_elec = 0.0, E_mag = 0.0; - Curl->Mult(E, B); + Curl.Mult(E.Real(), B.Real()); + Curl.Mult(E.Imag(), B.Imag()); B *= -1.0 / (1i * omega); postop.SetEGridFunction(E); postop.SetBGridFunction(B); @@ -244,16 +244,15 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator // Allocate negative curl matrix for postprocessing the B-field and vectors for the // high-dimensional field solution. - auto Curl = spaceop.GetCurlMatrix(); - ComplexVector E(Curl->Width()), B(Curl->Height()); + const auto &Curl = spaceop.GetCurlMatrix(); + ComplexVector E(Curl.Width()), B(Curl.Height()); E = 0.0; B = 0.0; // Initialize structures for storing and reducing the results of error estimation. CurlFluxErrorEstimator estimator( spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol, - iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold, - iodata.solver.pa_discrete_interp); + iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold); ErrorIndicator indicator; // Configure the PROM operator which performs the parameter space sampling and basis @@ -337,7 +336,8 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator // PostOperator for all postprocessing operations. BlockTimer bt4(Timer::POSTPRO); double E_elec = 0.0, E_mag = 0.0; - Curl->Mult(E, B); + Curl.Mult(E.Real(), B.Real()); + Curl.Mult(E.Imag(), B.Imag()); B *= -1.0 / (1i * omega); postop.SetEGridFunction(E); postop.SetBGridFunction(B); diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index 63ad9a074..6b45a209a 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -35,12 +35,12 @@ EigenSolver::Solve(const std::vector> &mesh) cons auto K = spaceop.GetStiffnessMatrix(Operator::DIAG_ONE); auto C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); auto M = spaceop.GetMassMatrix(Operator::DIAG_ZERO); - auto Curl = spaceop.GetCurlMatrix(); + const auto &Curl = spaceop.GetCurlMatrix(); SaveMetadata(spaceop.GetNDSpaces()); // Configure objects for postprocessing. PostOperator postop(iodata, spaceop, "eigenmode"); - ComplexVector E(Curl->Width()), B(Curl->Height()); + ComplexVector E(Curl.Width()), B(Curl.Height()); // Define and configure the eigensolver to solve the eigenvalue problem: // (K + λ C + λ² M) u = 0 or K u = -λ² M u @@ -166,7 +166,7 @@ EigenSolver::Solve(const std::vector> &mesh) cons spaceop.GetMaterialOp(), spaceop.GetNDSpace(), spaceop.GetH1Spaces(), spaceop.GetAuxBdrTDofLists(), iodata.solver.linear.divfree_tol, iodata.solver.linear.divfree_max_it, divfree_verbose, - iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp); + iodata.solver.pa_order_threshold); eigen->SetDivFreeProjector(*divfree); } @@ -192,9 +192,10 @@ EigenSolver::Solve(const std::vector> &mesh) cons eigen->SetInitialSpace(v0); // Copies the vector // Debug - // auto Grad = spaceop.GetGradMatrix(); + // const auto &Grad = spaceop.GetGradMatrix(); // ComplexVector r0(Grad->Width()); - // Grad->MultTranspose(v0, r0); + // Grad.MultTranspose(v0.Real(), r0.Real()); + // Grad.MultTranspose(v0.Imag(), r0.Imag()); // r0.Print(); } @@ -260,8 +261,7 @@ EigenSolver::Solve(const std::vector> &mesh) cons // Calculate and record the error indicators. CurlFluxErrorEstimator estimator( spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol, - iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold, - iodata.solver.pa_discrete_interp); + iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold); ErrorIndicator indicator; for (int i = 0; i < iodata.solver.eigenmode.n; i++) { @@ -296,7 +296,8 @@ EigenSolver::Solve(const std::vector> &mesh) cons // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in // PostOperator for all postprocessing operations. eigen->GetEigenvector(i, E); - Curl->Mult(E, B); + Curl.Mult(E.Real(), B.Real()); + Curl.Mult(E.Imag(), B.Imag()); B *= -1.0 / (1i * omega); postop.SetEGridFunction(E); postop.SetBGridFunction(B); diff --git a/palace/drivers/electrostaticsolver.cpp b/palace/drivers/electrostaticsolver.cpp index 19de2fcbc..4da1be60d 100644 --- a/palace/drivers/electrostaticsolver.cpp +++ b/palace/drivers/electrostaticsolver.cpp @@ -86,21 +86,21 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop, // charges from the prescribed voltage to get C directly as: // Q_i = ∫ ρ dV = ∫ ∇ ⋅ (ε E) dV = ∫ (ε E) ⋅ n dS // and C_ij = Q_i/V_j. The energy formulation avoids having to locally integrate E = -∇V. - auto Grad = laplaceop.GetGradMatrix(); + const auto &Grad = laplaceop.GetGradMatrix(); const std::map> &terminal_sources = laplaceop.GetSources(); int nstep = static_cast(terminal_sources.size()); mfem::DenseMatrix C(nstep), Cm(nstep); - Vector E(Grad->Height()), Vij(Grad->Width()); + Vector E(Grad.Height()), Vij(Grad.Width()); if (iodata.solver.electrostatic.n_post > 0) { Mpi::Print("\n"); } // Calculate and record the error indicators. - GradFluxErrorEstimator estimator( - laplaceop.GetMaterialOp(), laplaceop.GetH1Spaces(), - iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0, - iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp); + GradFluxErrorEstimator estimator(laplaceop.GetMaterialOp(), laplaceop.GetH1Spaces(), + iodata.solver.linear.estimator_tol, + iodata.solver.linear.estimator_max_it, 0, + iodata.solver.pa_order_threshold); ErrorIndicator indicator; for (int i = 0; i < nstep; i++) { @@ -113,7 +113,7 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop, // Compute E = -∇V on the true dofs, and set the internal GridFunctions in PostOperator // for all postprocessing operations. E = 0.0; - Grad->AddMult(V[i], E, -1.0); + Grad.AddMult(V[i], E, -1.0); postop.SetEGridFunction(E); postop.SetVGridFunction(V[i]); double Ue = postop.GetEFieldEnergy(); @@ -151,7 +151,7 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop, { linalg::AXPBYPCZ(1.0, V[i], 1.0, V[j], 0.0, Vij); E = 0.0; - Grad->AddMult(Vij, E, -1.0); + Grad.AddMult(Vij, E, -1.0); postop.SetEGridFunction(E); double Ue = postop.GetEFieldEnergy(); C(i, j) = Ue - 0.5 * (C(i, i) + C(j, j)); diff --git a/palace/drivers/magnetostaticsolver.cpp b/palace/drivers/magnetostaticsolver.cpp index 1fd9dfa14..00c40afd2 100644 --- a/palace/drivers/magnetostaticsolver.cpp +++ b/palace/drivers/magnetostaticsolver.cpp @@ -88,11 +88,11 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, // Φ_i = ∫ B ⋅ n_j dS // and M_ij = Φ_i/I_j. The energy formulation avoids having to locally integrate B = // ∇ x A. - auto Curl = curlcurlop.GetCurlMatrix(); + const auto &Curl = curlcurlop.GetCurlMatrix(); const SurfaceCurrentOperator &surf_j_op = curlcurlop.GetSurfaceCurrentOp(); int nstep = static_cast(surf_j_op.Size()); mfem::DenseMatrix M(nstep), Mm(nstep); - Vector B(Curl->Height()), Aij(Curl->Width()); + Vector B(Curl.Height()), Aij(Curl.Width()); Vector Iinc(nstep); if (iodata.solver.magnetostatic.n_post > 0) { @@ -103,7 +103,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, CurlFluxErrorEstimator estimator( curlcurlop.GetMaterialOp(), curlcurlop.GetNDSpaces(), iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0, - iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp); + iodata.solver.pa_order_threshold); ErrorIndicator indicator; for (int i = 0; i < nstep; i++) { @@ -120,7 +120,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, // Compute B = ∇ x A on the true dofs, and set the internal GridFunctions in // PostOperator for all postprocessing operations. - Curl->Mult(A[i], B); + Curl.Mult(A[i], B); postop.SetBGridFunction(B); postop.SetAGridFunction(A[i]); double Um = postop.GetHFieldEnergy(); @@ -157,7 +157,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, else if (j > i) { linalg::AXPBYPCZ(1.0, A[i], 1.0, A[j], 0.0, Aij); - Curl->Mult(Aij, B); + Curl.Mult(Aij, B); postop.SetBGridFunction(B); double Um = postop.GetHFieldEnergy(); M(i, j) = Um / (Iinc(i) * Iinc(j)) - diff --git a/palace/drivers/transientsolver.cpp b/palace/drivers/transientsolver.cpp index f1d192776..3c3a50a1a 100644 --- a/palace/drivers/transientsolver.cpp +++ b/palace/drivers/transientsolver.cpp @@ -80,8 +80,7 @@ TransientSolver::Solve(const std::vector> &mesh) // Initialize structures for storing and reducing the results of error estimation. CurlFluxErrorEstimator estimator( spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol, - iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold, - iodata.solver.pa_discrete_interp); + iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold); ErrorIndicator indicator; // Main time integration loop. diff --git a/palace/fem/fespace.cpp b/palace/fem/fespace.cpp index 09bb83c8a..9617f1436 100644 --- a/palace/fem/fespace.cpp +++ b/palace/fem/fespace.cpp @@ -3,6 +3,9 @@ #include "fespace.hpp" +#include "fem/bilinearform.hpp" +#include "fem/integrator.hpp" +#include "linalg/rap.hpp" #include "utils/omp.hpp" namespace palace @@ -24,4 +27,90 @@ std::size_t FiniteElementSpace::GetId() const return id; } +const Operator &AuxiliaryFiniteElementSpace::BuildDiscreteInterpolator() const +{ + // G is always partially assembled. + const int dim = GetParMesh()->Dimension(); + const auto aux_map_type = FEColl()->GetMapType(dim); + const auto primal_map_type = primal_fespace.FEColl()->GetMapType(dim); + if (aux_map_type == mfem::FiniteElement::VALUE && + primal_map_type == mfem::FiniteElement::H_CURL) + { + // Discrete gradient interpolator + DiscreteLinearOperator interp(*this, primal_fespace); + interp.AddDomainInterpolator(); + G = std::make_unique(interp.Assemble(), *this, primal_fespace, true); + } + else if (primal_map_type == mfem::FiniteElement::VALUE && + aux_map_type == mfem::FiniteElement::H_CURL) + { + // Discrete gradient interpolator (spaces reversed) + DiscreteLinearOperator interp(primal_fespace, *this); + interp.AddDomainInterpolator(); + G = std::make_unique(interp.Assemble(), primal_fespace, *this, true); + } + else if (aux_map_type == mfem::FiniteElement::H_CURL && + primal_map_type == mfem::FiniteElement::H_DIV) + { + // Discrete curl interpolator + DiscreteLinearOperator interp(*this, primal_fespace); + interp.AddDomainInterpolator(); + G = std::make_unique(interp.Assemble(), *this, primal_fespace, true); + } + else if (primal_map_type == mfem::FiniteElement::H_CURL && + aux_map_type == mfem::FiniteElement::H_DIV) + { + // Discrete curl interpolator (spaces reversed) + DiscreteLinearOperator interp(primal_fespace, *this); + interp.AddDomainInterpolator(); + G = std::make_unique(interp.Assemble(), primal_fespace, *this, true); + } + else if (aux_map_type == mfem::FiniteElement::H_DIV && + primal_map_type == mfem::FiniteElement::INTEGRAL) + { + // Discrete divergence interpolator + DiscreteLinearOperator interp(*this, primal_fespace); + interp.AddDomainInterpolator(); + G = std::make_unique(interp.Assemble(), *this, primal_fespace, true); + } + else if (primal_map_type == mfem::FiniteElement::H_DIV && + aux_map_type == mfem::FiniteElement::INTEGRAL) + { + // Discrete divergence interpolator (spaces reversed) + DiscreteLinearOperator interp(primal_fespace, *this); + interp.AddDomainInterpolator(); + G = std::make_unique(interp.Assemble(), primal_fespace, *this, true); + } + else + { + MFEM_ABORT("Unsupported trial/test FE spaces for AuxiliaryFiniteElementSpace discrete " + "interpolator!"); + } + + return *G; +} + +const Operator &FiniteElementSpaceHierarchy::BuildProlongationAtLevel(std::size_t l) const +{ + // P is always partially assembled. + MFEM_VERIFY(l >= 0 && l < GetNumLevels() - 1, + "Can only construct a finite element space prolongation with more than one " + "space in the hierarchy!"); + if (fespaces[l]->GetParMesh() != fespaces[l + 1]->GetParMesh()) + { + P[l] = std::make_unique( + std::make_unique(*fespaces[l], *fespaces[l + 1]), + *fespaces[l], *fespaces[l + 1], true); + } + else + { + DiscreteLinearOperator p(*fespaces[l], *fespaces[l + 1]); + p.AddDomainInterpolator(); + P[l] = + std::make_unique(p.Assemble(), *fespaces[l], *fespaces[l + 1], true); + } + + return *P[l]; +} + } // namespace palace diff --git a/palace/fem/fespace.hpp b/palace/fem/fespace.hpp index d6bcdea44..82748938d 100644 --- a/palace/fem/fespace.hpp +++ b/palace/fem/fespace.hpp @@ -4,7 +4,10 @@ #ifndef PALACE_FEM_FESPACE_HPP #define PALACE_FEM_FESPACE_HPP +#include +#include #include +#include "linalg/operator.hpp" namespace palace { @@ -35,6 +38,167 @@ class FiniteElementSpace : public mfem::ParFiniteElementSpace std::size_t GetId() const; }; +// +// An AuxiliaryFiniteElement space is just a FiniteElementSpace which allows for lazy +// construction of the interpolation operator (discrete gradient or curl) from the primal +// space to this one. +// +class AuxiliaryFiniteElementSpace : public FiniteElementSpace +{ +private: + const FiniteElementSpace &primal_fespace; + mutable std::unique_ptr G; + + const Operator &BuildDiscreteInterpolator() const; + +public: + template + AuxiliaryFiniteElementSpace(const FiniteElementSpace &primal_fespace, T &&...args) + : FiniteElementSpace(std::forward(args)...), primal_fespace(primal_fespace) + { + } + + // Return the discrete gradient or discrete curl matrix interpolating from the auxiliary + // to the primal space, constructing it on the fly as necessary. + const Operator &GetDiscreteInterpolator() const + { + return G ? *G : BuildDiscreteInterpolator(); + } +}; + +// +// A collection of FiniteElementSpace objects constructed on the same mesh with the ability +// to construct the prolongation operators between them as needed. +// +class FiniteElementSpaceHierarchy +{ +protected: + std::vector> fespaces; + mutable std::vector> P; + + const Operator &BuildProlongationAtLevel(std::size_t l) const; + +public: + FiniteElementSpaceHierarchy() = default; + FiniteElementSpaceHierarchy(std::unique_ptr &&fespace) + { + AddLevel(std::move(fespace)); + } + + auto GetNumLevels() const { return fespaces.size(); } + auto size() const { return GetNumLevels(); } + bool empty() const { return GetNumLevels() == 0; } + + virtual void AddLevel(std::unique_ptr &&fespace) + { + fespaces.push_back(std::move(fespace)); + P.push_back(nullptr); + } + + FiniteElementSpace &GetFESpaceAtLevel(std::size_t l) + { + MFEM_ASSERT(l >= 0 && l < GetNumLevels(), + "Out of bounds request for finite element space at level " << l << "!"); + return *fespaces[l]; + } + const FiniteElementSpace &GetFESpaceAtLevel(std::size_t l) const + { + MFEM_ASSERT(l >= 0 && l < GetNumLevels(), + "Out of bounds request for finite element space at level " << l << "!"); + return *fespaces[l]; + } + + FiniteElementSpace &GetFinestFESpace() + { + MFEM_ASSERT(!empty(), "Out of bounds request for finite element space at level 0!"); + return *fespaces.back(); + } + const FiniteElementSpace &GetFinestFESpace() const + { + MFEM_ASSERT(!empty(), "Out of bounds request for finite element space at level 0!"); + return *fespaces.back(); + } + + const Operator &GetProlongationAtLevel(std::size_t l) const + { + MFEM_ASSERT(l >= 0 && l < GetNumLevels() - 1, + "Out of bounds request for finite element space prolongation at level " + << l << "!"); + return P[l] ? *P[l] : BuildProlongationAtLevel(l); + } + + std::vector GetProlongationOperators() const + { + std::vector P_(GetNumLevels() - 1); + for (std::size_t l = 0; l < P_.size(); l++) + { + P_[l] = &GetProlongationAtLevel(l); + } + return P_; + } +}; + +// +// A special type of FiniteElementSpaceHierarchy where all members are auxiliary finite +// element spaces. +// +class AuxiliaryFiniteElementSpaceHierarchy : public FiniteElementSpaceHierarchy +{ +public: + AuxiliaryFiniteElementSpaceHierarchy() = default; + AuxiliaryFiniteElementSpaceHierarchy( + std::unique_ptr &&fespace) + { + AddLevel(std::move(fespace)); + } + + void AddLevel(std::unique_ptr &&fespace) override + { + MFEM_ABORT("All finite element spaces in an AuxiliaryFiniteElementSpaceHierarchy must " + "inherit from AuxiliaryFiniteElementSpace!"); + } + + void AddLevel(std::unique_ptr &&fespace) + { + // Guarantees that every object in fespaces is an AuxiliaryFiniteElementSpace. + fespaces.push_back(std::move(fespace)); + P.push_back(nullptr); + } + + AuxiliaryFiniteElementSpace &GetAuxiliaryFESpaceAtLevel(std::size_t l) + { + return *static_cast(&GetFESpaceAtLevel(l)); + } + const AuxiliaryFiniteElementSpace &GetAuxiliaryFESpaceAtLevel(std::size_t l) const + { + return *static_cast(&GetFESpaceAtLevel(l)); + } + + AuxiliaryFiniteElementSpace &GetFinestAuxiliaryFESpace() + { + return *static_cast(&GetFinestFESpace()); + } + const AuxiliaryFiniteElementSpace &GetFinestAuxiliaryFESpace() const + { + return *static_cast(&GetFinestFESpace()); + } + + const Operator &GetDiscreteInterpolatorAtLevel(std::size_t l) const + { + return GetAuxiliaryFESpaceAtLevel(l).GetDiscreteInterpolator(); + } + + std::vector GetDiscreteInterpolators() const + { + std::vector G_(GetNumLevels()); + for (std::size_t l = 0; l < G_.size(); l++) + { + G_[l] = &GetDiscreteInterpolatorAtLevel(l); + } + return G_; + } +}; + } // namespace palace #endif // PALACE_FEM_FESPACE_HPP diff --git a/palace/fem/integrator.hpp b/palace/fem/integrator.hpp index f2a40a173..2f232033a 100644 --- a/palace/fem/integrator.hpp +++ b/palace/fem/integrator.hpp @@ -506,6 +506,9 @@ using GradientInterpolator = DiscreteInterpolator; // Interpolator for the discrete curl map from an H(curl) space to an H(div) space. using CurlInterpolator = DiscreteInterpolator; +// Interpolator for the discrete divergence map from an H(div) space to an L2 space. +using DivergenceInterpolator = DiscreteInterpolator; + // Similar to MFEM's VectorFEBoundaryTangentLFIntegrator for ND spaces, but instead of // computing (n x f, v), this just computes (f, v). Also eliminates the a and b quadrature // parameters and uses fem::GetDefaultIntegrationOrder instead. diff --git a/palace/fem/multigrid.hpp b/palace/fem/multigrid.hpp index 48306c875..e8c6fbaef 100644 --- a/palace/fem/multigrid.hpp +++ b/palace/fem/multigrid.hpp @@ -7,11 +7,7 @@ #include #include #include -#include "fem/bilinearform.hpp" #include "fem/fespace.hpp" -#include "fem/integrator.hpp" -#include "linalg/operator.hpp" -#include "linalg/rap.hpp" #include "utils/iodata.hpp" namespace palace::fem @@ -23,9 +19,10 @@ namespace palace::fem // Construct sequence of FECollection objects. template -std::vector> inline ConstructFECollections( - int p, int dim, int mg_max_levels, - config::LinearSolverData::MultigridCoarsenType mg_coarsen_type, bool mat_lor) +inline std::vector> +ConstructFECollections(int p, int dim, int mg_max_levels, + config::LinearSolverData::MultigridCoarsenType mg_coarsen_type, + bool mat_lor) { // If the solver will use a LOR preconditioner, we need to construct with a specific basis // type. @@ -70,17 +67,15 @@ std::vector> inline ConstructFECollections( } } std::reverse(fecs.begin(), fecs.end()); + return fecs; } -// Construct a hierarchy of finite element spaces given a sequence of meshes and -// finite element collections. Dirichlet boundary conditions are additionally -// marked. +// Construct a hierarchy of finite element spaces given a sequence of meshes and finite +// element collections. Additionally, Dirichlet boundary conditions are marked. template -inline std::unique_ptr -ConstructFiniteElementSpaceHierarchy( - int mg_max_levels, bool mg_legacy_transfer, int pa_order_threshold, - bool pa_discrete_interp, const std::vector> &mesh, +inline FiniteElementSpaceHierarchy ConstructFiniteElementSpaceHierarchy( + int mg_max_levels, const std::vector> &mesh, const std::vector> &fecs, const mfem::Array *dbc_marker = nullptr, std::vector> *dbc_tdof_lists = nullptr) @@ -90,54 +85,95 @@ ConstructFiniteElementSpaceHierarchy( "Empty mesh or FE collection for FE space construction!"); int coarse_mesh_l = std::max(0, static_cast(mesh.size() + fecs.size()) - 1 - mg_max_levels); - auto *fespace = new FiniteElementSpace(mesh[coarse_mesh_l].get(), fecs[0].get()); + FiniteElementSpaceHierarchy fespaces( + std::make_unique(mesh[coarse_mesh_l].get(), fecs[0].get())); if (dbc_marker && dbc_tdof_lists) { - fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back()); + fespaces.GetFinestFESpace().GetEssentialTrueDofs(*dbc_marker, + dbc_tdof_lists->emplace_back()); } - auto fespaces = std::make_unique( - mesh[coarse_mesh_l].get(), fespace, false, true); // h-refinement for (std::size_t l = coarse_mesh_l + 1; l < mesh.size(); l++) { - fespace = new FiniteElementSpace(mesh[l].get(), fecs[0].get()); + fespaces.AddLevel(std::make_unique(mesh[l].get(), fecs[0].get())); if (dbc_marker && dbc_tdof_lists) { - fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back()); + fespaces.GetFinestFESpace().GetEssentialTrueDofs(*dbc_marker, + dbc_tdof_lists->emplace_back()); } - auto *P = new ParOperator( - std::make_unique(fespaces->GetFinestFESpace(), *fespace), - fespaces->GetFinestFESpace(), *fespace, true); - fespaces->AddLevel(mesh[l].get(), fespace, P, false, true, true); } // p-refinement for (std::size_t l = 1; l < fecs.size(); l++) { - fespace = new FiniteElementSpace(mesh.back().get(), fecs[l].get()); + fespaces.AddLevel( + std::make_unique(mesh.back().get(), fecs[l].get())); if (dbc_marker && dbc_tdof_lists) { - fespace->GetEssentialTrueDofs(*dbc_marker, dbc_tdof_lists->emplace_back()); + fespaces.GetFinestFESpace().GetEssentialTrueDofs(*dbc_marker, + dbc_tdof_lists->emplace_back()); + } + } + + return fespaces; +} + +// Similar to ConstructFiniteElementSpaceHierarchy above, but in this case the finite +// element space at each level is an auxiliary space associated with the coresponding level +// of the provided finite element space objects. +template +inline AuxiliaryFiniteElementSpaceHierarchy ConstructAuxiliaryFiniteElementSpaceHierarchy( + const FiniteElementSpaceHierarchy &primal_fespaces, + const std::vector> &fecs, + const mfem::Array *dbc_marker = nullptr, + std::vector> *dbc_tdof_lists = nullptr) +{ + MFEM_VERIFY(!primal_fespaces.empty() && !fecs.empty() && + (!dbc_tdof_lists || dbc_tdof_lists->empty()), + "Empty mesh or FE collection for FE space construction!"); + mfem::ParMesh *mesh = primal_fespaces.GetFESpaceAtLevel(0).GetParMesh(); + AuxiliaryFiniteElementSpaceHierarchy fespaces( + std::make_unique(primal_fespaces.GetFESpaceAtLevel(0), + mesh, fecs[0].get())); + if (dbc_marker && dbc_tdof_lists) + { + fespaces.GetFinestFESpace().GetEssentialTrueDofs(*dbc_marker, + dbc_tdof_lists->emplace_back()); + } + + // h-refinement + std::size_t l; + for (l = 1; l < primal_fespaces.size(); l++) + { + if (primal_fespaces.GetFESpaceAtLevel(l).GetParMesh() == mesh) + { + break; } - ParOperator *P; - if (!mg_legacy_transfer) + fespaces.AddLevel(std::make_unique( + primal_fespaces.GetFESpaceAtLevel(l), + primal_fespaces.GetFESpaceAtLevel(l).GetParMesh(), fecs[0].get())); + if (dbc_marker && dbc_tdof_lists) { - constexpr bool skip_zeros_interp = true; - DiscreteLinearOperator p(fespaces->GetFinestFESpace(), *fespace); - p.AddDomainInterpolator(); - P = new ParOperator( - p.Assemble(pa_discrete_interp ? pa_order_threshold : 99, skip_zeros_interp), - fespaces->GetFinestFESpace(), *fespace, true); + fespaces.GetFinestFESpace().GetEssentialTrueDofs(*dbc_marker, + dbc_tdof_lists->emplace_back()); } - else + mesh = primal_fespaces.GetFESpaceAtLevel(l).GetParMesh(); + } + + // p-refinement + const auto l0 = l - 1; + for (; l < primal_fespaces.size(); l++) + { + fespaces.AddLevel(std::make_unique( + primal_fespaces.GetFESpaceAtLevel(l), mesh, fecs[l - l0].get())); + if (dbc_marker && dbc_tdof_lists) { - P = new ParOperator( - std::make_unique(fespaces->GetFinestFESpace(), *fespace), - fespaces->GetFinestFESpace(), *fespace, true); + fespaces.GetFinestFESpace().GetEssentialTrueDofs(*dbc_marker, + dbc_tdof_lists->emplace_back()); } - fespaces->AddLevel(mesh.back().get(), fespace, P, false, true, true); } + return fespaces; } diff --git a/palace/linalg/ams.cpp b/palace/linalg/ams.cpp index 770ac6637..21925df15 100644 --- a/palace/linalg/ams.cpp +++ b/palace/linalg/ams.cpp @@ -4,14 +4,15 @@ #include "ams.hpp" #include "fem/bilinearform.hpp" +#include "fem/fespace.hpp" #include "fem/integrator.hpp" #include "linalg/rap.hpp" namespace palace { -HypreAmsSolver::HypreAmsSolver(const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace, int cycle_it, +HypreAmsSolver::HypreAmsSolver(const FiniteElementSpace &nd_fespace, + const AuxiliaryFiniteElementSpace &h1_fespace, int cycle_it, int smooth_it, int agg_coarsen, bool vector_interp, bool op_singular, int print) : mfem::HypreSolver(), @@ -47,17 +48,18 @@ HypreAmsSolver::~HypreAmsSolver() } void HypreAmsSolver::ConstructAuxiliaryMatrices( - const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace) + const FiniteElementSpace &nd_fespace, const AuxiliaryFiniteElementSpace &h1_fespace) { // Set up the auxiliary space objects for the preconditioner. Mostly the same as MFEM's // HypreAMS:Init. Start with the discrete gradient matrix. { constexpr bool skip_zeros_interp = true; - DiscreteLinearOperator grad(h1_fespace, nd_fespace); - grad.AddDomainInterpolator(); - ParOperator RAP_G(grad.FullAssemble(skip_zeros_interp), h1_fespace, nd_fespace, true); - G = RAP_G.StealParallelAssemble(); + const auto *PtGP = + dynamic_cast(&h1_fespace.GetDiscreteInterpolator()); + MFEM_VERIFY( + PtGP, + "HypreAmsSolver requires the discrete gradient matrix as a ParOperator operator!"); + G = &PtGP->ParallelAssemble(skip_zeros_interp); } // Vertex coordinates for the lowest order case, or Nedelec interpolation matrix or @@ -65,9 +67,9 @@ void HypreAmsSolver::ConstructAuxiliaryMatrices( mfem::ParMesh &mesh = *h1_fespace.GetParMesh(); if (h1_fespace.GetMaxElementOrder() == 1) { - mfem::ParGridFunction x_coord(const_cast(&h1_fespace)), - y_coord(const_cast(&h1_fespace)), - z_coord(const_cast(&h1_fespace)); + mfem::ParGridFunction x_coord(const_cast(&h1_fespace)), + y_coord(const_cast(&h1_fespace)), + z_coord(const_cast(&h1_fespace)); if (mesh.GetNodes()) { mesh.GetNodes()->GetNodalValues(x_coord, 1); @@ -118,7 +120,7 @@ void HypreAmsSolver::ConstructAuxiliaryMatrices( mfem::ParFiniteElementSpace h1d_fespace(&mesh, h1_fespace.FEColl(), space_dim, mfem::Ordering::byVDIM); mfem::DiscreteLinearOperator pi(&h1d_fespace, - const_cast(&nd_fespace)); + const_cast(&nd_fespace)); pi.AddDomainInterpolator(new mfem::IdentityInterpolator); pi.SetAssemblyLevel(mfem::AssemblyLevel::LEGACY); pi.Assemble(); diff --git a/palace/linalg/ams.hpp b/palace/linalg/ams.hpp index dc0e95e22..50386bd2a 100644 --- a/palace/linalg/ams.hpp +++ b/palace/linalg/ams.hpp @@ -12,6 +12,9 @@ namespace palace { +class AuxiliaryFiniteElementSpace; +class FiniteElementSpace; + // // A wrapper for Hypre's AMS solver. // @@ -28,8 +31,8 @@ class HypreAmsSolver : public mfem::HypreSolver // Control print level for debugging. const int print; - // Discrete gradient matrix. - std::unique_ptr G; + // Discrete gradient matrix (not owned). + const mfem::HypreParMatrix *G; // Nedelec interpolation matrix and its components, or, for p = 1, the mesh vertex // coordinates. @@ -37,8 +40,8 @@ class HypreAmsSolver : public mfem::HypreSolver std::unique_ptr x, y, z; // Helper function to set up the auxiliary objects required by the AMS solver. - void ConstructAuxiliaryMatrices(const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace); + void ConstructAuxiliaryMatrices(const FiniteElementSpace &nd_fespace, + const AuxiliaryFiniteElementSpace &h1_fespace); // Helper function to construct and configure the AMS solver. void InitializeSolver(); @@ -46,12 +49,12 @@ class HypreAmsSolver : public mfem::HypreSolver public: // Constructor requires the ND space, but will construct the H1 and (H1)ᵈ spaces // internally as needed. - HypreAmsSolver(const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace, int cycle_it, int smooth_it, + HypreAmsSolver(const FiniteElementSpace &nd_fespace, + const AuxiliaryFiniteElementSpace &h1_fespace, int cycle_it, int smooth_it, int agg_coarsen, bool vector_interp, bool op_singular, int print); HypreAmsSolver(const IoData &iodata, bool coarse_solver, - const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace, int print) + const FiniteElementSpace &nd_fespace, + const AuxiliaryFiniteElementSpace &h1_fespace, int print) : HypreAmsSolver( nd_fespace, h1_fespace, coarse_solver ? 1 : iodata.solver.linear.mg_cycle_it, iodata.solver.linear.mg_smooth_it, diff --git a/palace/linalg/distrelaxation.cpp b/palace/linalg/distrelaxation.cpp index 91a9c23cf..7ef9d3663 100644 --- a/palace/linalg/distrelaxation.cpp +++ b/palace/linalg/distrelaxation.cpp @@ -4,8 +4,6 @@ #include "distrelaxation.hpp" #include -#include "fem/bilinearform.hpp" -#include "fem/integrator.hpp" #include "linalg/chebyshev.hpp" #include "linalg/rap.hpp" @@ -14,20 +12,11 @@ namespace palace template DistRelaxationSmoother::DistRelaxationSmoother( - const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace, int smooth_it, int cheby_smooth_it, - int cheby_order, double cheby_sf_max, double cheby_sf_min, bool cheby_4th_kind, - int pa_order_threshold, bool pa_discrete_interp) - : Solver(), pc_it(smooth_it), A(nullptr), A_G(nullptr), dbc_tdof_list_G(nullptr) + const Operator &G, int smooth_it, int cheby_smooth_it, int cheby_order, + double cheby_sf_max, double cheby_sf_min, bool cheby_4th_kind) + : Solver(), pc_it(smooth_it), G(&G), A(nullptr), A_G(nullptr), + dbc_tdof_list_G(nullptr) { - // Construct discrete gradient matrix for the auxiliary space. - constexpr bool skip_zeros_interp = true; - DiscreteLinearOperator grad(h1_fespace, nd_fespace); - grad.AddDomainInterpolator(); - G = std::make_unique( - grad.Assemble(pa_discrete_interp ? pa_order_threshold : 99, skip_zeros_interp), - h1_fespace, nd_fespace, true); - // Initialize smoothers. if (cheby_4th_kind) { diff --git a/palace/linalg/distrelaxation.hpp b/palace/linalg/distrelaxation.hpp index dca92ee0b..c447f7186 100644 --- a/palace/linalg/distrelaxation.hpp +++ b/palace/linalg/distrelaxation.hpp @@ -14,7 +14,6 @@ namespace mfem template class Array; -class ParFiniteElementSpace; } // namespace mfem @@ -36,13 +35,13 @@ class DistRelaxationSmoother : public Solver // Number of smoother iterations. const int pc_it; + // Discrete gradient matrix (not owned). + const Operator *G; + // System matrix and its projection GᵀAG (not owned). const OperType *A, *A_G; const mfem::Array *dbc_tdof_list_G; - // Discrete gradient matrix. - std::unique_ptr G; - // Point smoother objects for each matrix. mutable std::unique_ptr> B; std::unique_ptr> B_G; @@ -51,11 +50,9 @@ class DistRelaxationSmoother : public Solver mutable VecType r, x_G, y_G; public: - DistRelaxationSmoother(const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &h1_fespace, int smooth_it, - int cheby_smooth_it, int cheby_order, double cheby_sf_max, - double cheby_sf_min, bool cheby_4th_kind, int pa_order_threshold, - bool pa_discrete_interp); + DistRelaxationSmoother(const Operator &G, int smooth_it, int cheby_smooth_it, + int cheby_order, double cheby_sf_max, double cheby_sf_min, + bool cheby_4th_kind); void SetOperator(const OperType &op) override { diff --git a/palace/linalg/divfree.cpp b/palace/linalg/divfree.cpp index 41dbda7a6..4266bb856 100644 --- a/palace/linalg/divfree.cpp +++ b/palace/linalg/divfree.cpp @@ -7,6 +7,7 @@ #include #include "fem/bilinearform.hpp" #include "fem/coefficient.hpp" +#include "fem/fespace.hpp" #include "fem/integrator.hpp" #include "linalg/amg.hpp" #include "linalg/gmg.hpp" @@ -18,18 +19,17 @@ namespace palace { DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, + const FiniteElementSpace &nd_fespace, + const AuxiliaryFiniteElementSpaceHierarchy &h1_fespaces, const std::vector> &h1_bdr_tdof_lists, - double tol, int max_it, int print, int pa_order_threshold, - bool pa_discrete_interp) + double tol, int max_it, int print, int pa_order_threshold) { constexpr bool skip_zeros = false; constexpr auto MatType = MaterialPropertyType::PERMITTIVITY_REAL; MaterialPropertyCoefficient epsilon_func(mat_op); { auto M_mg = std::make_unique(h1_fespaces.GetNumLevels()); - for (int l = 0; l < h1_fespaces.GetNumLevels(); l++) + for (std::size_t l = 0; l < h1_fespaces.GetNumLevels(); l++) { // Force coarse level operator to be fully assembled always. const auto &h1_fespace_l = h1_fespaces.GetFESpaceAtLevel(l); @@ -49,14 +49,7 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, std::make_unique(weakdiv.Assemble(pa_order_threshold, skip_zeros), nd_fespace, h1_fespaces.GetFinestFESpace(), false); } - { - constexpr bool skip_zeros_interp = true; - DiscreteLinearOperator grad(h1_fespaces.GetFinestFESpace(), nd_fespace); - grad.AddDomainInterpolator(); - Grad = std::make_unique( - grad.Assemble(pa_discrete_interp ? pa_order_threshold : 99, skip_zeros_interp), - h1_fespaces.GetFinestFESpace(), nd_fespace, true); - } + Grad = &h1_fespaces.GetFinestAuxiliaryFESpace().GetDiscreteInterpolator(); bdr_tdof_list_M = &h1_bdr_tdof_lists.back(); // The system matrix for the projection is real and SPD. @@ -68,8 +61,8 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, const int mg_smooth_order = std::max(h1_fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( - std::move(amg), h1_fespaces, nullptr, 1, 1, mg_smooth_order, 1.0, 0.0, true, - pa_order_threshold, pa_discrete_interp); + std::move(amg), h1_fespaces.GetProlongationOperators(), nullptr, 1, 1, + mg_smooth_order, 1.0, 0.0, true); } else { diff --git a/palace/linalg/divfree.hpp b/palace/linalg/divfree.hpp index fe5b360ba..ead820895 100644 --- a/palace/linalg/divfree.hpp +++ b/palace/linalg/divfree.hpp @@ -15,13 +15,14 @@ namespace mfem template class Array; -class ParFiniteElementSpaceHierarchy; } // namespace mfem namespace palace { +class AuxiliaryFiniteElementSpaceHierarchy; +class FiniteElementSpace; class MaterialOperator; // @@ -33,7 +34,8 @@ class DivFreeSolver { private: // Operators for the divergence-free projection. - std::unique_ptr WeakDiv, Grad, M; + std::unique_ptr WeakDiv, M; + const Operator *Grad; const mfem::Array *bdr_tdof_list_M; // Linear solver for the projected linear system (Gᵀ M G) y = x. @@ -43,11 +45,10 @@ class DivFreeSolver mutable Vector psi, rhs; public: - DivFreeSolver(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, + DivFreeSolver(const MaterialOperator &mat_op, const FiniteElementSpace &nd_fespace, + const AuxiliaryFiniteElementSpaceHierarchy &h1_fespaces, const std::vector> &h1_bdr_tdof_lists, double tol, - int max_it, int print, int pa_order_threshold, bool pa_discrete_interp); + int max_it, int print, int pa_order_threshold); // Given a vector of Nedelec dofs for an arbitrary vector field, compute the Nedelec dofs // of the irrotational portion of this vector field. The resulting vector will satisfy diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index 60b5e3034..93242ecf4 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -6,7 +6,6 @@ #include #include "fem/bilinearform.hpp" #include "fem/coefficient.hpp" -#include "fem/fespace.hpp" #include "fem/integrator.hpp" #include "linalg/amg.hpp" #include "linalg/gmg.hpp" @@ -23,14 +22,13 @@ namespace palace namespace { -std::unique_ptr -GetMassMatrix(const mfem::ParFiniteElementSpaceHierarchy &fespaces, int pa_order_threshold, - int skip_zeros) +std::unique_ptr GetMassMatrix(const FiniteElementSpaceHierarchy &fespaces, + int pa_order_threshold, int skip_zeros) { const int dim = fespaces.GetFinestFESpace().GetParMesh()->Dimension(); const auto type = fespaces.GetFinestFESpace().FEColl()->GetRangeType(dim); auto M = std::make_unique(fespaces.GetNumLevels()); - for (int l = 0; l < fespaces.GetNumLevels(); l++) + for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++) { // Force coarse level operator to be fully assembled always. const auto &fespace_l = fespaces.GetFESpaceAtLevel(l); @@ -53,9 +51,8 @@ GetMassMatrix(const mfem::ParFiniteElementSpaceHierarchy &fespaces, int pa_order } std::unique_ptr -ConfigureLinearSolver(const mfem::ParFiniteElementSpaceHierarchy &fespaces, double tol, - int max_it, int print, int pa_order_threshold, - bool pa_discrete_interp) +ConfigureLinearSolver(const FiniteElementSpaceHierarchy &fespaces, double tol, int max_it, + int print) { // The system matrix for the projection is real and SPD. auto amg = @@ -66,8 +63,8 @@ ConfigureLinearSolver(const mfem::ParFiniteElementSpaceHierarchy &fespaces, doub const int mg_smooth_order = std::max(fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( - std::move(amg), fespaces, nullptr, 1, 1, mg_smooth_order, 1.0, 0.0, true, - pa_order_threshold, pa_discrete_interp); + std::move(amg), fespaces.GetProlongationOperators(), nullptr, 1, 1, mg_smooth_order, + 1.0, 0.0, true); } else { @@ -87,37 +84,35 @@ ConfigureLinearSolver(const mfem::ParFiniteElementSpaceHierarchy &fespaces, doub } // namespace FluxProjector::FluxProjector(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, - double tol, int max_it, int print, int pa_order_threshold, - bool pa_discrete_interp) + const FiniteElementSpaceHierarchy &nd_fespaces, double tol, + int max_it, int print, int pa_order_threshold) { BlockTimer bt(Timer::CONSTRUCTESTIMATOR); - constexpr bool skip_zeros = false; { + // Flux operator is always partially assembled. constexpr auto MatType = MaterialPropertyType::INV_PERMEABILITY; MaterialPropertyCoefficient muinv_func(mat_op); BilinearForm flux(nd_fespaces.GetFinestFESpace()); flux.AddDomainIntegrator(muinv_func); Flux = std::make_unique(flux.Assemble(), nd_fespaces.GetFinestFESpace()); } + constexpr int skip_zeros = false; M = GetMassMatrix(nd_fespaces, pa_order_threshold, skip_zeros); - ksp = ConfigureLinearSolver(nd_fespaces, tol, max_it, print, pa_order_threshold, - pa_discrete_interp); + ksp = ConfigureLinearSolver(nd_fespaces, tol, max_it, print); ksp->SetOperators(*M, *M); rhs.SetSize(nd_fespaces.GetFinestFESpace().GetTrueVSize()); } FluxProjector::FluxProjector(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, - const mfem::ParFiniteElementSpace &h1d_fespace, double tol, - int max_it, int print, int pa_order_threshold, - bool pa_discrete_interp) + const FiniteElementSpaceHierarchy &h1_fespaces, + const FiniteElementSpace &h1d_fespace, double tol, int max_it, + int print, int pa_order_threshold) { BlockTimer bt(Timer::CONSTRUCTESTIMATOR); - constexpr bool skip_zeros = false; { + // Flux operator is always partially assembled. constexpr auto MatType = MaterialPropertyType::PERMITTIVITY_REAL; MaterialPropertyCoefficient epsilon_func(mat_op); BilinearForm flux(h1_fespaces.GetFinestFESpace(), h1d_fespace); @@ -125,10 +120,10 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, Flux = std::make_unique(flux.Assemble(), h1_fespaces.GetFinestFESpace(), h1d_fespace, false); } + constexpr int skip_zeros = false; M = GetMassMatrix(h1_fespaces, pa_order_threshold, skip_zeros); - ksp = ConfigureLinearSolver(h1_fespaces, tol, max_it, print, pa_order_threshold, - pa_discrete_interp); + ksp = ConfigureLinearSolver(h1_fespaces, tol, max_it, print); ksp->SetOperators(*M, *M); rhs.SetSize(h1d_fespace.GetTrueVSize()); @@ -177,12 +172,12 @@ void FluxProjector::Mult(const VecType &x, VecType &y) const template CurlFluxErrorEstimator::CurlFluxErrorEstimator( - const MaterialOperator &mat_op, mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, - double tol, int max_it, int print, int pa_order_threshold, bool pa_discrete_interp) + const MaterialOperator &mat_op, const FiniteElementSpaceHierarchy &nd_fespaces, + double tol, int max_it, int print, int pa_order_threshold) : mat_op(mat_op), nd_fespace(nd_fespaces.GetFinestFESpace()), - projector(mat_op, nd_fespaces, tol, max_it, print, pa_order_threshold, - pa_discrete_interp), - F(nd_fespace.GetTrueVSize()), F_gf(&nd_fespace), U_gf(&nd_fespace) + projector(mat_op, nd_fespaces, tol, max_it, print, pa_order_threshold), + F(nd_fespace.GetTrueVSize()), F_gf(const_cast(&nd_fespace)), + U_gf(const_cast(&nd_fespace)) { } @@ -295,15 +290,15 @@ ErrorIndicator CurlFluxErrorEstimator::ComputeIndicators(const VecType } GradFluxErrorEstimator::GradFluxErrorEstimator( - const MaterialOperator &mat_op, mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, - double tol, int max_it, int print, int pa_order_threshold, bool pa_discrete_interp) + const MaterialOperator &mat_op, const FiniteElementSpaceHierarchy &h1_fespaces, + double tol, int max_it, int print, int pa_order_threshold) : mat_op(mat_op), h1_fespace(h1_fespaces.GetFinestFESpace()), h1d_fespace(std::make_unique( h1_fespace.GetParMesh(), h1_fespace.FEColl(), h1_fespace.GetParMesh()->SpaceDimension(), mfem::Ordering::byNODES)), - projector(mat_op, h1_fespaces, *h1d_fespace, tol, max_it, print, pa_order_threshold, - pa_discrete_interp), - F(h1d_fespace->GetTrueVSize()), F_gf(h1d_fespace.get()), U_gf(&h1_fespace) + projector(mat_op, h1_fespaces, *h1d_fespace, tol, max_it, print, pa_order_threshold), + F(h1d_fespace->GetTrueVSize()), F_gf(h1d_fespace.get()), + U_gf(const_cast(&h1_fespace)) { } diff --git a/palace/linalg/errorestimator.hpp b/palace/linalg/errorestimator.hpp index ad7c54835..3ab5333db 100644 --- a/palace/linalg/errorestimator.hpp +++ b/palace/linalg/errorestimator.hpp @@ -7,6 +7,7 @@ #include #include #include "fem/errorindicator.hpp" +#include "fem/fespace.hpp" #include "linalg/ksp.hpp" #include "linalg/operator.hpp" #include "linalg/vector.hpp" @@ -38,12 +39,12 @@ class FluxProjector public: FluxProjector(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, double tol, - int max_it, int print, int pa_order_threshold, bool pa_discrete_interp); + const FiniteElementSpaceHierarchy &nd_fespaces, double tol, int max_it, + int print, int pa_order_threshold); FluxProjector(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, - const mfem::ParFiniteElementSpace &h1d_fespace, double tol, int max_it, - int print, int pa_order_threshold, bool pa_discrete_interp); + const FiniteElementSpaceHierarchy &h1_fespaces, + const FiniteElementSpace &h1d_fespace, double tol, int max_it, int print, + int pa_order_threshold); template void Mult(const VecType &x, VecType &y) const; @@ -62,7 +63,7 @@ class CurlFluxErrorEstimator const MaterialOperator &mat_op; // Finite element space used to represent U and F. - mfem::ParFiniteElementSpace &nd_fespace; + const FiniteElementSpace &nd_fespace; // Global L2 projection solver. FluxProjector projector; @@ -73,9 +74,8 @@ class CurlFluxErrorEstimator public: CurlFluxErrorEstimator(const MaterialOperator &mat_op, - mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, double tol, - int max_it, int print, int pa_order_threshold, - bool pa_discrete_interp); + const FiniteElementSpaceHierarchy &nd_fespaces, double tol, + int max_it, int print, int pa_order_threshold); // Compute elemental error indicators given a vector of true DOF. ErrorIndicator ComputeIndicators(const VecType &U) const; @@ -96,10 +96,10 @@ class GradFluxErrorEstimator const MaterialOperator &mat_op; // Finite element space used to represent U. - mfem::ParFiniteElementSpace &h1_fespace; + const FiniteElementSpace &h1_fespace; // Vector H1 space used to represent the components of F, ordered by component. - std::unique_ptr h1d_fespace; + std::unique_ptr h1d_fespace; // Global L2 projection solver. FluxProjector projector; @@ -110,9 +110,8 @@ class GradFluxErrorEstimator public: GradFluxErrorEstimator(const MaterialOperator &mat_op, - mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, double tol, - int max_it, int print, int pa_order_threshold, - bool pa_discrete_interp); + const FiniteElementSpaceHierarchy &h1_fespaces, double tol, + int max_it, int print, int pa_order_threshold); // Compute elemental error indicators given a vector of true DOF. ErrorIndicator ComputeIndicators(const Vector &U) const; diff --git a/palace/linalg/gmg.cpp b/palace/linalg/gmg.cpp index 209b03636..9abcf6508 100644 --- a/palace/linalg/gmg.cpp +++ b/palace/linalg/gmg.cpp @@ -15,41 +15,35 @@ namespace palace template GeometricMultigridSolver::GeometricMultigridSolver( std::unique_ptr> &&coarse_solver, - const mfem::ParFiniteElementSpaceHierarchy &fespaces, - const mfem::ParFiniteElementSpaceHierarchy *aux_fespaces, int cycle_it, int smooth_it, - int cheby_order, double cheby_sf_max, double cheby_sf_min, bool cheby_4th_kind, - int pa_order_threshold, bool pa_discrete_interp) - : Solver(), pc_it(cycle_it), A(fespaces.GetNumLevels()), - P(fespaces.GetNumLevels() - 1), dbc_tdof_lists(fespaces.GetNumLevels() - 1), - B(fespaces.GetNumLevels()), X(fespaces.GetNumLevels()), Y(fespaces.GetNumLevels()), - R(fespaces.GetNumLevels()) + const std::vector &P, const std::vector *G, + int cycle_it, int smooth_it, int cheby_order, double cheby_sf_max, double cheby_sf_min, + bool cheby_4th_kind) + : Solver(), pc_it(cycle_it), P(P.begin(), P.end()), A(P.size() + 1), + dbc_tdof_lists(P.size()), B(P.size() + 1), X(P.size() + 1), Y(P.size() + 1), + R(P.size() + 1) { // Configure levels of geometric coarsening. Multigrid vectors will be configured at first // call to Mult. The multigrid operator size is set based on the finest space dimension. - const int n_levels = fespaces.GetNumLevels(); + const auto n_levels = P.size() + 1; MFEM_VERIFY(n_levels > 0, "Empty finite element space hierarchy during multigrid solver setup!"); - - // Configure prolongation operators. - for (int l = 0; l < n_levels - 1; l++) - { - P[l] = fespaces.GetProlongationAtLevel(l); - } + MFEM_VERIFY(!G || G->size() == n_levels, + "Invalid input for distributive relaxation smoother auxiliary space transfer " + "operators (mismatch in number of levels)!"); // Use the supplied level 0 (coarse) solver. B[0] = std::move(coarse_solver); // Configure level smoothers. Use distributive relaxation smoothing if an auxiliary // finite element space was provided. - for (int l = 1; l < n_levels; l++) + for (std::size_t l = 1; l < n_levels; l++) { - if (aux_fespaces) + if (G) { const int cheby_smooth_it = 1; B[l] = std::make_unique>( - fespaces.GetFESpaceAtLevel(l), aux_fespaces->GetFESpaceAtLevel(l), smooth_it, - cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min, cheby_4th_kind, - pa_order_threshold, pa_discrete_interp); + *(*G)[l], smooth_it, cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min, + cheby_4th_kind); } else { @@ -79,12 +73,12 @@ void GeometricMultigridSolver::SetOperator(const OperType &op) MFEM_VERIFY(mg_op, "GeometricMultigridSolver requires a MultigridOperator or " "ComplexMultigridOperator argument provided to SetOperator!"); - const int n_levels = static_cast(A.size()); + const auto n_levels = A.size(); MFEM_VERIFY( mg_op->GetNumLevels() == n_levels && (!mg_op->HasAuxiliaryOperators() || mg_op->GetNumAuxiliaryLevels() == n_levels), "Invalid number of levels for operators in multigrid solver setup!"); - for (int l = 0; l < n_levels; l++) + for (std::size_t l = 0; l < n_levels; l++) { A[l] = &mg_op->GetOperatorAtLevel(l); MFEM_VERIFY( @@ -128,7 +122,7 @@ template void GeometricMultigridSolver::Mult(const VecType &x, VecType &y) const { // Initialize. - const int n_levels = static_cast(A.size()); + const auto n_levels = A.size(); MFEM_ASSERT(!this->initial_guess, "Geometric multigrid solver does not use initial guess!"); MFEM_ASSERT(n_levels > 1 || pc_it == 1, diff --git a/palace/linalg/gmg.hpp b/palace/linalg/gmg.hpp index a912dd05b..958cc23c6 100644 --- a/palace/linalg/gmg.hpp +++ b/palace/linalg/gmg.hpp @@ -16,7 +16,6 @@ namespace mfem template class Array; -class ParFiniteElementSpaceHierarchy; } // namespace mfem @@ -37,9 +36,11 @@ class GeometricMultigridSolver : public Solver // Number of V-cycles per preconditioner application. const int pc_it; - // System matrices at each multigrid level and prolongation operators (not owned). - std::vector A; + // Prolongation operators (not owned). std::vector P; + + // System matrices at each multigrid level (not owned). + std::vector A; std::vector *> dbc_tdof_lists; // Smoothers for each level. Coarse level solver is B[0]. @@ -54,21 +55,19 @@ class GeometricMultigridSolver : public Solver public: GeometricMultigridSolver(std::unique_ptr> &&coarse_solver, - const mfem::ParFiniteElementSpaceHierarchy &fespaces, - const mfem::ParFiniteElementSpaceHierarchy *aux_fespaces, - int cycle_it, int smooth_it, int cheby_order, - double cheby_sf_max, double cheby_sf_min, bool cheby_4th_kind, - int pa_order_threshold, bool pa_discrete_interp); + const std::vector &P, + const std::vector *G, int cycle_it, + int smooth_it, int cheby_order, double cheby_sf_max, + double cheby_sf_min, bool cheby_4th_kind); GeometricMultigridSolver(const IoData &iodata, std::unique_ptr> &&coarse_solver, - const mfem::ParFiniteElementSpaceHierarchy &fespaces, - const mfem::ParFiniteElementSpaceHierarchy *aux_fespaces) + const std::vector &P, + const std::vector *G = nullptr) : GeometricMultigridSolver( - std::move(coarse_solver), fespaces, aux_fespaces, - iodata.solver.linear.mg_cycle_it, iodata.solver.linear.mg_smooth_it, - iodata.solver.linear.mg_smooth_order, iodata.solver.linear.mg_smooth_sf_max, - iodata.solver.linear.mg_smooth_sf_min, iodata.solver.linear.mg_smooth_cheby_4th, - iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp) + std::move(coarse_solver), P, G, iodata.solver.linear.mg_cycle_it, + iodata.solver.linear.mg_smooth_it, iodata.solver.linear.mg_smooth_order, + iodata.solver.linear.mg_smooth_sf_max, iodata.solver.linear.mg_smooth_sf_min, + iodata.solver.linear.mg_smooth_cheby_4th) { } diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp index 8f6b89aa4..a35fd9e99 100644 --- a/palace/linalg/hcurl.cpp +++ b/palace/linalg/hcurl.cpp @@ -6,6 +6,7 @@ #include #include "fem/bilinearform.hpp" #include "fem/coefficient.hpp" +#include "fem/fespace.hpp" #include "fem/integrator.hpp" #include "linalg/ams.hpp" #include "linalg/gmg.hpp" @@ -17,11 +18,11 @@ namespace palace { WeightedHCurlNormSolver::WeightedHCurlNormSolver( - const MaterialOperator &mat_op, const mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, - const mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, + const MaterialOperator &mat_op, const FiniteElementSpaceHierarchy &nd_fespaces, + const AuxiliaryFiniteElementSpaceHierarchy &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, bool pa_discrete_interp) + int print, int pa_order_threshold) { constexpr bool skip_zeros = false; constexpr auto MatTypeMuInv = MaterialPropertyType::INV_PERMEABILITY; @@ -30,33 +31,33 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( MaterialPropertyCoefficient epsilon_func(mat_op); { auto A_mg = std::make_unique(nd_fespaces.GetNumLevels()); - for (int s = 0; s < 2; s++) + for (bool aux : {false, true}) { - const auto &fespaces = (s == 0) ? nd_fespaces : h1_fespaces; - const auto &dbc_tdof_lists = (s == 0) ? nd_dbc_tdof_lists : h1_dbc_tdof_lists; - for (int l = 0; l < fespaces.GetNumLevels(); l++) + const auto &fespaces = aux ? h1_fespaces : nd_fespaces; + const auto &dbc_tdof_lists = aux ? h1_dbc_tdof_lists : nd_dbc_tdof_lists; + for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++) { // Force coarse level operator to be fully assembled always. const auto &fespace_l = fespaces.GetFESpaceAtLevel(l); BilinearForm a(fespace_l); - if (s == 0) + if (aux) { - a.AddDomainIntegrator(muinv_func, epsilon_func); + a.AddDomainIntegrator(epsilon_func); } else { - a.AddDomainIntegrator(epsilon_func); + a.AddDomainIntegrator(muinv_func, epsilon_func); } auto A_l = std::make_unique( a.Assemble((l > 0) ? pa_order_threshold : 99, skip_zeros), fespace_l); A_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE); - if (s == 0) + if (aux) { - A_mg->AddOperator(std::move(A_l)); + A_mg->AddAuxiliaryOperator(std::move(A_l)); } else { - A_mg->AddAuxiliaryOperator(std::move(A_l)); + A_mg->AddOperator(std::move(A_l)); } } } @@ -66,16 +67,17 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( // The system matrix K + M is real and SPD. We use Hypre's AMS solver as the coarse-level // multigrid solve. auto ams = std::make_unique>(std::make_unique( - nd_fespaces.GetFESpaceAtLevel(0), h1_fespaces.GetFESpaceAtLevel(0), 1, 1, 1, false, - false, 0)); + nd_fespaces.GetFESpaceAtLevel(0), h1_fespaces.GetAuxiliaryFESpaceAtLevel(0), 1, 1, 1, + false, false, 0)); std::unique_ptr> pc; if (nd_fespaces.GetNumLevels() > 1) { + const auto G = h1_fespaces.GetDiscreteInterpolators(); const int mg_smooth_order = std::max(nd_fespaces.GetFinestFESpace().GetMaxElementOrder(), 2); pc = std::make_unique>( - std::move(ams), nd_fespaces, &h1_fespaces, 1, 1, mg_smooth_order, 1.0, 0.0, true, - pa_order_threshold, pa_discrete_interp); + std::move(ams), nd_fespaces.GetProlongationOperators(), &G, 1, 1, mg_smooth_order, + 1.0, 0.0, true); } else { diff --git a/palace/linalg/hcurl.hpp b/palace/linalg/hcurl.hpp index c841bcab7..579d2b126 100644 --- a/palace/linalg/hcurl.hpp +++ b/palace/linalg/hcurl.hpp @@ -15,13 +15,14 @@ namespace mfem template class Array; -class ParFiniteElementSpaceHierarchy; } // namespace mfem namespace palace { +class AuxiliaryFiniteElementSpaceHierarchy; +class FiniteElementSpaceHierarchy; class MaterialOperator; // @@ -38,12 +39,11 @@ class WeightedHCurlNormSolver public: WeightedHCurlNormSolver(const MaterialOperator &mat_op, - const mfem::ParFiniteElementSpaceHierarchy &nd_fespaces, - const mfem::ParFiniteElementSpaceHierarchy &h1_fespaces, + const FiniteElementSpaceHierarchy &nd_fespaces, + const AuxiliaryFiniteElementSpaceHierarchy &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, - bool pa_discrete_interp); + double tol, int max_it, int print, int pa_order_threshold); const Operator &GetOperator() { return *A; } diff --git a/palace/linalg/ksp.cpp b/palace/linalg/ksp.cpp index 404c316cd..8414adb47 100644 --- a/palace/linalg/ksp.cpp +++ b/palace/linalg/ksp.cpp @@ -4,6 +4,7 @@ #include "ksp.hpp" #include +#include "fem/fespace.hpp" #include "linalg/amg.hpp" #include "linalg/ams.hpp" #include "linalg/gmg.hpp" @@ -106,8 +107,8 @@ std::unique_ptr> ConfigureKrylovSolver(MPI_Comm comm, template std::unique_ptr> ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata, - mfem::ParFiniteElementSpaceHierarchy &fespaces, - mfem::ParFiniteElementSpaceHierarchy *aux_fespaces) + const FiniteElementSpaceHierarchy &fespaces, + const AuxiliaryFiniteElementSpaceHierarchy *aux_fespaces) { // Create the real-valued solver first. std::unique_ptr pc0; @@ -120,9 +121,9 @@ ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata, // space (in which case fespaces.GetNumLevels() == 1). MFEM_VERIFY(aux_fespaces, "AMS solver relies on both primary space " "and auxiliary spaces for construction!"); - pc0 = std::make_unique(iodata, fespaces.GetNumLevels() > 1, - fespaces.GetFESpaceAtLevel(0), - aux_fespaces->GetFESpaceAtLevel(0), print); + pc0 = std::make_unique( + iodata, fespaces.GetNumLevels() > 1, fespaces.GetFESpaceAtLevel(0), + aux_fespaces->GetAuxiliaryFESpaceAtLevel(0), print); break; case config::LinearSolverData::Type::BOOMER_AMG: pc0 = std::make_unique(iodata, fespaces.GetNumLevels() > 1, print); @@ -176,13 +177,14 @@ ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata, { MFEM_VERIFY(aux_fespaces, "Multigrid with auxiliary space smoothers requires both " "primary space and auxiliary spaces for construction!"); - return std::make_unique>(iodata, std::move(pc), - fespaces, aux_fespaces); + const auto G = aux_fespaces->GetDiscreteInterpolators(); + return std::make_unique>( + iodata, std::move(pc), fespaces.GetProlongationOperators(), &G); } else { - return std::make_unique>(iodata, std::move(pc), - fespaces, nullptr); + return std::make_unique>( + iodata, std::move(pc), fespaces.GetProlongationOperators()); } } else @@ -194,9 +196,9 @@ ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata, } // namespace template -BaseKspSolver::BaseKspSolver(const IoData &iodata, - mfem::ParFiniteElementSpaceHierarchy &fespaces, - mfem::ParFiniteElementSpaceHierarchy *aux_fespaces) +BaseKspSolver::BaseKspSolver( + const IoData &iodata, const FiniteElementSpaceHierarchy &fespaces, + const AuxiliaryFiniteElementSpaceHierarchy *aux_fespaces) : BaseKspSolver( ConfigureKrylovSolver(fespaces.GetFinestFESpace().GetComm(), iodata), ConfigurePreconditionerSolver(fespaces.GetFinestFESpace().GetComm(), diff --git a/palace/linalg/ksp.hpp b/palace/linalg/ksp.hpp index b9580e3bc..4af3f4da6 100644 --- a/palace/linalg/ksp.hpp +++ b/palace/linalg/ksp.hpp @@ -10,16 +10,11 @@ #include "linalg/operator.hpp" #include "linalg/solver.hpp" -namespace mfem -{ - -class ParFiniteElementSpaceHierarchy; - -} // namespace mfem - namespace palace { +class AuxiliaryFiniteElementSpaceHierarchy; +class FiniteElementSpaceHierarchy; class IoData; // @@ -45,8 +40,8 @@ class BaseKspSolver mutable int ksp_mult, ksp_mult_it; public: - BaseKspSolver(const IoData &iodata, mfem::ParFiniteElementSpaceHierarchy &fespaces, - mfem::ParFiniteElementSpaceHierarchy *aux_fespaces = nullptr); + BaseKspSolver(const IoData &iodata, const FiniteElementSpaceHierarchy &fespaces, + const AuxiliaryFiniteElementSpaceHierarchy *aux_fespaces = nullptr); BaseKspSolver(std::unique_ptr> &&ksp, std::unique_ptr> &&pc); diff --git a/palace/linalg/operator.hpp b/palace/linalg/operator.hpp index bf8a9ce56..9ba7cec89 100644 --- a/palace/linalg/operator.hpp +++ b/palace/linalg/operator.hpp @@ -264,7 +264,7 @@ class BaseMultigridOperator : public OperType std::vector> ops, aux_ops; public: - BaseMultigridOperator(int l) : OperType(0) + BaseMultigridOperator(std::size_t l) : OperType(0) { ops.reserve(l); aux_ops.reserve(l); @@ -284,19 +284,19 @@ class BaseMultigridOperator : public OperType bool HasAuxiliaryOperators() const { return !aux_ops.empty(); } - int GetNumLevels() const { return static_cast(ops.size()); } - int GetNumAuxiliaryLevels() const { return static_cast(aux_ops.size()); } + auto GetNumLevels() const { return ops.size(); } + auto GetNumAuxiliaryLevels() const { return aux_ops.size(); } const OperType &GetFinestOperator() const { return *ops.back(); } const OperType &GetFinestAuxiliaryOperator() const { return *aux_ops.back(); } - const OperType &GetOperatorAtLevel(int l) const + const OperType &GetOperatorAtLevel(std::size_t l) const { MFEM_ASSERT(l >= 0 && l < GetNumLevels(), "Out of bounds multigrid level operator requested!"); return *ops[l]; } - const OperType &GetAuxiliaryOperatorAtLevel(int l) const + const OperType &GetAuxiliaryOperatorAtLevel(std::size_t l) const { MFEM_ASSERT(l < GetNumAuxiliaryLevels(), "Out of bounds multigrid level auxiliary operator requested!"); diff --git a/palace/linalg/rap.cpp b/palace/linalg/rap.cpp index 7cd50cf17..30ee39246 100644 --- a/palace/linalg/rap.cpp +++ b/palace/linalg/rap.cpp @@ -3,6 +3,8 @@ #include "rap.hpp" +#include "fem/bilinearform.hpp" + namespace palace { @@ -104,7 +106,7 @@ void ParOperator::AssembleDiagonal(Vector &diag) const } } -mfem::HypreParMatrix &ParOperator::ParallelAssemble() const +mfem::HypreParMatrix &ParOperator::ParallelAssemble(bool skip_zeros) const { if (RAP) { @@ -112,8 +114,17 @@ mfem::HypreParMatrix &ParOperator::ParallelAssemble() const } // Build the square or rectangular assembled HypreParMatrix. - mfem::SparseMatrix *sA = dynamic_cast(A); - MFEM_VERIFY(sA, "ParOperator::ParallelAssemble requires A as a SparseMatrix!"); + auto *sA = dynamic_cast(A); + std::unique_ptr data_sA; + if (!sA) + { + auto *cA = dynamic_cast(A); + MFEM_VERIFY(cA, "ParOperator::ParallelAssemble requires A as an mfem::SparseMatrix or " + "ceed::Operator!"); + data_sA = use_R ? DiscreteLinearOperator::FullAssemble(*cA, skip_zeros) + : BilinearForm::FullAssemble(*cA, skip_zeros); + sA = data_sA.get(); + } if (&trial_fespace == &test_fespace) { mfem::HypreParMatrix *hA = diff --git a/palace/linalg/rap.hpp b/palace/linalg/rap.hpp index e45300304..cbf3e73e8 100644 --- a/palace/linalg/rap.hpp +++ b/palace/linalg/rap.hpp @@ -88,13 +88,13 @@ class ParOperator : public Operator void AssembleDiagonal(Vector &diag) const override; // Assemble the operator as a parallel sparse matrix. The memory associated with the - // local operator is not freed. - mfem::HypreParMatrix &ParallelAssemble() const; + // local operator is free'd. + mfem::HypreParMatrix &ParallelAssemble(bool skip_zeros = false) const; // Steal the assembled parallel sparse matrix. - std::unique_ptr StealParallelAssemble() const + std::unique_ptr StealParallelAssemble(bool skip_zeros = false) const { - ParallelAssemble(); + ParallelAssemble(skip_zeros); return std::move(RAP); } diff --git a/palace/models/curlcurloperator.cpp b/palace/models/curlcurloperator.cpp index 76371d84e..de0bc4973 100644 --- a/palace/models/curlcurloperator.cpp +++ b/palace/models/curlcurloperator.cpp @@ -5,7 +5,6 @@ #include "fem/bilinearform.hpp" #include "fem/coefficient.hpp" -#include "fem/fespace.hpp" #include "fem/integrator.hpp" #include "fem/multigrid.hpp" #include "linalg/rap.hpp" @@ -84,13 +83,10 @@ CurlCurlOperator::CurlCurlOperator(const IoData &iodata, rt_fec(std::make_unique(iodata.solver.order - 1, mesh.back()->Dimension())), nd_fespaces(fem::ConstructFiniteElementSpaceHierarchy( - iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer, - pa_order_threshold, pa_discrete_interp, mesh, nd_fecs, &dbc_marker, - &dbc_tdof_lists)), - h1_fespaces(fem::ConstructFiniteElementSpaceHierarchy( - iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer, - pa_order_threshold, pa_discrete_interp, mesh, h1_fecs, nullptr, nullptr)), - rt_fespace(std::make_unique(mesh.back().get(), rt_fec.get())), + iodata.solver.linear.mg_max_levels, mesh, nd_fecs, &dbc_marker, &dbc_tdof_lists)), + h1_fespaces(fem::ConstructAuxiliaryFiniteElementSpaceHierarchy( + nd_fespaces, h1_fecs)), + rt_fespace(nd_fespaces.GetFinestFESpace(), mesh.back().get(), rt_fec.get()), mat_op(iodata, *mesh.back()), surf_j_op(iodata, GetH1Space()) { // Finalize setup. @@ -127,7 +123,7 @@ std::unique_ptr CurlCurlOperator::GetStiffnessMatrix() Mpi::Print("\nAssembling multigrid hierarchy:\n"); } auto K = std::make_unique(GetNDSpaces().GetNumLevels()); - for (int l = 0; l < GetNDSpaces().GetNumLevels(); l++) + for (std::size_t l = 0; l < GetNDSpaces().GetNumLevels(); l++) { // Force coarse level operator to be fully assembled always. const auto &nd_fespace_l = GetNDSpaces().GetFESpaceAtLevel(l); @@ -163,16 +159,6 @@ std::unique_ptr CurlCurlOperator::GetStiffnessMatrix() return K; } -std::unique_ptr CurlCurlOperator::GetCurlMatrix() -{ - constexpr bool skip_zeros_interp = true; - DiscreteLinearOperator curl(GetNDSpace(), GetRTSpace()); - curl.AddDomainInterpolator(); - return std::make_unique( - curl.Assemble(pa_discrete_interp ? pa_order_threshold : 99, skip_zeros_interp), - GetNDSpace(), GetRTSpace(), true); -} - void CurlCurlOperator::GetExcitationVector(int idx, Vector &RHS) { // Assemble the surface current excitation +J. The SurfaceCurrentOperator assembles -J diff --git a/palace/models/curlcurloperator.hpp b/palace/models/curlcurloperator.hpp index 0331317d1..773f1be07 100644 --- a/palace/models/curlcurloperator.hpp +++ b/palace/models/curlcurloperator.hpp @@ -7,6 +7,7 @@ #include #include #include +#include "fem/fespace.hpp" #include "linalg/operator.hpp" #include "linalg/vector.hpp" #include "models/materialoperator.hpp" @@ -41,8 +42,9 @@ class CurlCurlOperator std::vector> nd_fecs; std::vector> h1_fecs; std::unique_ptr rt_fec; - std::unique_ptr nd_fespaces, h1_fespaces; - std::unique_ptr rt_fespace; + FiniteElementSpaceHierarchy nd_fespaces; + AuxiliaryFiniteElementSpaceHierarchy h1_fespaces; + AuxiliaryFiniteElementSpace rt_fespace; // Operator for domain material properties. MaterialOperator mat_op; @@ -61,14 +63,16 @@ class CurlCurlOperator const auto &GetSurfaceCurrentOp() const { return surf_j_op; } // Return the parallel finite element space objects. - auto &GetNDSpaces() { return *nd_fespaces; } - auto &GetNDSpace() { return nd_fespaces->GetFinestFESpace(); } - const auto &GetNDSpace() const { return nd_fespaces->GetFinestFESpace(); } - auto &GetH1Spaces() { return *h1_fespaces; } - auto &GetH1Space() { return h1_fespaces->GetFinestFESpace(); } - const auto &GetH1Space() const { return h1_fespaces->GetFinestFESpace(); } - auto &GetRTSpace() { return *rt_fespace; } - const auto &GetRTSpace() const { return *rt_fespace; } + auto &GetNDSpaces() { return nd_fespaces; } + const auto &GetNDSpaces() const { return nd_fespaces; } + auto &GetNDSpace() { return nd_fespaces.GetFinestFESpace(); } + const auto &GetNDSpace() const { return nd_fespaces.GetFinestFESpace(); } + auto &GetH1Spaces() { return h1_fespaces; } + const auto &GetH1Spaces() const { return h1_fespaces; } + auto &GetH1Space() { return h1_fespaces.GetFinestFESpace(); } + const auto &GetH1Space() const { return h1_fespaces.GetFinestFESpace(); } + auto &GetRTSpace() { return rt_fespace; } + const auto &GetRTSpace() const { return rt_fespace; } // Return the number of true (conforming) dofs on the finest ND space. auto GlobalTrueVSize() { return GetNDSpace().GlobalTrueVSize(); } @@ -78,7 +82,7 @@ class CurlCurlOperator std::unique_ptr GetStiffnessMatrix(); // Construct and return the discrete curl matrix. - std::unique_ptr GetCurlMatrix(); + const Operator &GetCurlMatrix() const { return GetRTSpace().GetDiscreteInterpolator(); } // Assemble the right-hand side source term vector for a current source applied on // specified excited boundaries. diff --git a/palace/models/laplaceoperator.cpp b/palace/models/laplaceoperator.cpp index 3a0d1ce6a..2b648dc78 100644 --- a/palace/models/laplaceoperator.cpp +++ b/palace/models/laplaceoperator.cpp @@ -5,7 +5,6 @@ #include "fem/bilinearform.hpp" #include "fem/coefficient.hpp" -#include "fem/fespace.hpp" #include "fem/integrator.hpp" #include "fem/multigrid.hpp" #include "linalg/rap.hpp" @@ -125,10 +124,8 @@ LaplaceOperator::LaplaceOperator(const IoData &iodata, nd_fec(std::make_unique(iodata.solver.order, mesh.back()->Dimension())), h1_fespaces(fem::ConstructFiniteElementSpaceHierarchy( - iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer, - pa_order_threshold, pa_discrete_interp, mesh, h1_fecs, &dbc_marker, - &dbc_tdof_lists)), - nd_fespace(std::make_unique(mesh.back().get(), nd_fec.get())), + iodata.solver.linear.mg_max_levels, mesh, h1_fecs, &dbc_marker, &dbc_tdof_lists)), + nd_fespace(h1_fespaces.GetFinestFESpace(), mesh.back().get(), nd_fec.get()), mat_op(iodata, *mesh.back()), source_attr_lists(ConstructSources(iodata)) { // Print essential BC information. @@ -150,7 +147,7 @@ std::unique_ptr LaplaceOperator::GetStiffnessMatrix() Mpi::Print("\nAssembling multigrid hierarchy:\n"); } auto K = std::make_unique(GetH1Spaces().GetNumLevels()); - for (int l = 0; l < GetH1Spaces().GetNumLevels(); l++) + for (std::size_t l = 0; l < GetH1Spaces().GetNumLevels(); l++) { // Force coarse level operator to be fully assembled always. const auto &h1_fespace_l = GetH1Spaces().GetFESpaceAtLevel(l); @@ -186,16 +183,6 @@ std::unique_ptr LaplaceOperator::GetStiffnessMatrix() return K; } -std::unique_ptr LaplaceOperator::GetGradMatrix() -{ - constexpr bool skip_zeros_interp = true; - DiscreteLinearOperator grad(GetH1Space(), GetNDSpace()); - grad.AddDomainInterpolator(); - return std::make_unique( - grad.Assemble(pa_discrete_interp ? pa_order_threshold : 99, skip_zeros_interp), - GetH1Space(), GetNDSpace(), true); -} - void LaplaceOperator::GetExcitationVector(int idx, const Operator &K, Vector &X, Vector &RHS) { diff --git a/palace/models/laplaceoperator.hpp b/palace/models/laplaceoperator.hpp index 7a58340a4..603c98414 100644 --- a/palace/models/laplaceoperator.hpp +++ b/palace/models/laplaceoperator.hpp @@ -8,6 +8,7 @@ #include #include #include +#include "fem/fespace.hpp" #include "linalg/operator.hpp" #include "linalg/vector.hpp" #include "models/materialoperator.hpp" @@ -38,8 +39,8 @@ class LaplaceOperator // electric field (Nedelec) on the given mesh. std::vector> h1_fecs; std::unique_ptr nd_fec; - std::unique_ptr h1_fespaces; - std::unique_ptr nd_fespace; + FiniteElementSpaceHierarchy h1_fespaces; + AuxiliaryFiniteElementSpace nd_fespace; // Operator for domain material properties. MaterialOperator mat_op; @@ -58,11 +59,12 @@ class LaplaceOperator const auto &GetSources() const { return source_attr_lists; } // Return the parallel finite element space objects. - auto &GetH1Spaces() { return *h1_fespaces; } - auto &GetH1Space() { return h1_fespaces->GetFinestFESpace(); } - const auto &GetH1Space() const { return h1_fespaces->GetFinestFESpace(); } - auto &GetNDSpace() { return *nd_fespace; } - const auto &GetNDSpace() const { return *nd_fespace; } + auto &GetH1Spaces() { return h1_fespaces; } + const auto &GetH1Spaces() const { return h1_fespaces; } + auto &GetH1Space() { return h1_fespaces.GetFinestFESpace(); } + const auto &GetH1Space() const { return h1_fespaces.GetFinestFESpace(); } + auto &GetNDSpace() { return nd_fespace; } + const auto &GetNDSpace() const { return nd_fespace; } // Return the number of true (conforming) dofs on the finest H1 space. auto GlobalTrueVSize() { return GetH1Space().GlobalTrueVSize(); } @@ -72,7 +74,7 @@ class LaplaceOperator std::unique_ptr GetStiffnessMatrix(); // Construct and return the discrete gradient matrix. - std::unique_ptr GetGradMatrix(); + const Operator &GetGradMatrix() const { return GetNDSpace().GetDiscreteInterpolator(); } // Assemble the solution boundary conditions and right-hand side vector for a nonzero // prescribed voltage on the specified surface index. diff --git a/palace/models/romoperator.cpp b/palace/models/romoperator.cpp index ee465c79e..78c9f5ea3 100644 --- a/palace/models/romoperator.cpp +++ b/palace/models/romoperator.cpp @@ -118,8 +118,7 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &spaceop) : spaceop 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, - iodata.solver.pa_discrete_interp); + iodata.solver.linear.max_it, curlcurl_verbose, iodata.solver.pa_order_threshold); } // The initial PROM basis is empty. Orthogonalization uses MGS by default, else CGS2. diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index a406117b6..cb3c22292 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -5,7 +5,7 @@ #include #include "fem/bilinearform.hpp" -#include "fem/fespace.hpp" +#include "fem/coefficient.hpp" #include "fem/integrator.hpp" #include "fem/multigrid.hpp" #include "linalg/rap.hpp" @@ -88,14 +88,11 @@ SpaceOperator::SpaceOperator(const IoData &iodata, rt_fec(std::make_unique(iodata.solver.order - 1, mesh.back()->Dimension())), nd_fespaces(fem::ConstructFiniteElementSpaceHierarchy( - iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer, - pa_order_threshold, pa_discrete_interp, mesh, nd_fecs, &dbc_marker, + iodata.solver.linear.mg_max_levels, mesh, nd_fecs, &dbc_marker, &nd_dbc_tdof_lists)), - h1_fespaces(fem::ConstructFiniteElementSpaceHierarchy( - iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_legacy_transfer, - pa_order_threshold, pa_discrete_interp, mesh, h1_fecs, &dbc_marker, - &h1_dbc_tdof_lists)), - rt_fespace(std::make_unique(mesh.back().get(), rt_fec.get())), + h1_fespaces(fem::ConstructAuxiliaryFiniteElementSpaceHierarchy( + nd_fespaces, h1_fecs, &dbc_marker, &h1_dbc_tdof_lists)), + rt_fespace(nd_fespaces.GetFinestFESpace(), mesh.back().get(), rt_fec.get()), mat_op(iodata, *mesh.back()), farfield_op(iodata, mat_op, *mesh.back()), surf_sigma_op(iodata, *mesh.back()), surf_z_op(iodata, *mesh.back()), lumped_port_op(iodata, GetH1Space()), @@ -135,7 +132,7 @@ void SpaceOperator::CheckBoundaryProperties() // aux_bdr_marker = 1; // Mark all boundaries (including material interfaces // // added during mesh preprocessing) // // As tested, this does not eliminate all DC modes! - for (int l = 0; l < GetH1Spaces().GetNumLevels(); l++) + for (std::size_t l = 0; l < GetH1Spaces().GetNumLevels(); l++) { GetH1Spaces().GetFESpaceAtLevel(l).GetEssentialTrueDofs( aux_bdr_marker, aux_bdr_tdof_lists.emplace_back()); @@ -178,9 +175,8 @@ void SpaceOperator::CheckBoundaryProperties() namespace { -void PrintHeader(mfem::ParFiniteElementSpace &h1_fespace, - mfem::ParFiniteElementSpace &nd_fespace, - mfem::ParFiniteElementSpace &rt_fespace, int pa_order_threshold, +void PrintHeader(const FiniteElementSpace &h1_fespace, const FiniteElementSpace &nd_fespace, + const FiniteElementSpace &rt_fespace, int pa_order_threshold, bool &print_hdr) { if (print_hdr) @@ -195,8 +191,8 @@ void PrintHeader(mfem::ParFiniteElementSpace &h1_fespace, } template -std::unique_ptr BuildOperator(const mfem::ParFiniteElementSpace &fespace, T1 *df, - T2 *f, T3 *dfb, T4 *fb, int pa_order_threshold, +std::unique_ptr BuildOperator(const FiniteElementSpace &fespace, T1 *df, T2 *f, + T3 *dfb, T4 *fb, int pa_order_threshold, int skip_zeros) { BilinearForm a(fespace); @@ -234,9 +230,8 @@ std::unique_ptr BuildOperator(const mfem::ParFiniteElementSpace &fespa } template -std::unique_ptr BuildAuxOperator(const mfem::ParFiniteElementSpace &fespace, - T1 *f, T2 *fb, int pa_order_threshold, - int skip_zeros) +std::unique_ptr BuildAuxOperator(const FiniteElementSpace &fespace, T1 *f, T2 *fb, + int pa_order_threshold, int skip_zeros) { BilinearForm a(fespace); if (f && !f->empty()) @@ -402,7 +397,7 @@ namespace auto BuildParSumOperator(int h, int w, double a0, double a1, double a2, const ParOperator *K, const ParOperator *C, const ParOperator *M, - const ParOperator *A2, const mfem::ParFiniteElementSpace &fespace) + const ParOperator *A2, const FiniteElementSpace &fespace) { auto sum = std::make_unique(h, w); if (K && a0 != 0.0) @@ -427,8 +422,7 @@ auto BuildParSumOperator(int h, int w, double a0, double a1, double a2, auto BuildParSumOperator(int h, int w, std::complex a0, std::complex a1, std::complex a2, const ComplexParOperator *K, const ComplexParOperator *C, const ComplexParOperator *M, - const ComplexParOperator *A2, - const mfem::ParFiniteElementSpace &fespace) + const ComplexParOperator *A2, const FiniteElementSpace &fespace) { // Block 2 x 2 equivalent-real formulation for each term in the sum: // [ sumr ] += [ ar -ai ] [ Ar ] @@ -619,15 +613,13 @@ namespace { auto BuildLevelOperator(const MultigridOperator &B, std::unique_ptr &&br, - std::unique_ptr &&bi, - const mfem::ParFiniteElementSpace &fespace) + std::unique_ptr &&bi, const FiniteElementSpace &fespace) { return std::make_unique(std::move(br), fespace); } auto BuildLevelOperator(const ComplexMultigridOperator &B, std::unique_ptr &&br, - std::unique_ptr &&bi, - const mfem::ParFiniteElementSpace &fespace) + std::unique_ptr &&bi, const FiniteElementSpace &fespace) { return std::make_unique(std::move(br), std::move(bi), fespace); } @@ -647,19 +639,18 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub MFEM_VERIFY(GetH1Spaces().GetNumLevels() == GetNDSpaces().GetNumLevels(), "Multigrid hierarchy mismatch for auxiliary space preconditioning!"); auto B = std::make_unique>(GetNDSpaces().GetNumLevels()); - for (int s = 0; s < 2; s++) + for (bool aux : {false, true}) { - auto &fespaces = (s == 0) ? GetNDSpaces() : GetH1Spaces(); - auto &dbc_tdof_lists = (s == 0) ? nd_dbc_tdof_lists : h1_dbc_tdof_lists; - for (int l = 0; l < fespaces.GetNumLevels(); l++) + auto &fespaces = aux ? GetH1Spaces() : GetNDSpaces(); + auto &dbc_tdof_lists = aux ? h1_dbc_tdof_lists : nd_dbc_tdof_lists; + for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++) { // Force coarse level operator to be fully assembled always. auto &fespace_l = fespaces.GetFESpaceAtLevel(l); if (print_prec_hdr) { - Mpi::Print(" Level {:d}{} (p = {:d}): {:d} unknowns", l, - (s == 0) ? "" : " (auxiliary)", fespace_l.GetMaxElementOrder(), - fespace_l.GlobalTrueVSize()); + Mpi::Print(" Level {:d}{} (p = {:d}): {:d} unknowns", l, aux ? " (auxiliary)" : "", + fespace_l.GetMaxElementOrder(), fespace_l.GlobalTrueVSize()); } const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); SumMatrixCoefficient dfr(sdim), fr(sdim), fi(sdim), fbr(sdim), fbi(sdim); @@ -691,17 +682,17 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub std::unique_ptr br, bi; if (!dfr.empty() || !fr.empty() || !dfbr.empty() || !fbr.empty()) { - br = (s == 0) ? BuildOperator(fespace_l, &dfr, &fr, &dfbr, &fbr, - (l > 0) ? pa_order_threshold : 99, skip_zeros) - : BuildAuxOperator(fespace_l, &fr, &fbr, - (l > 0) ? pa_order_threshold : 99, skip_zeros); + br = aux ? BuildAuxOperator(fespace_l, &fr, &fbr, (l > 0) ? pa_order_threshold : 99, + skip_zeros) + : BuildOperator(fespace_l, &dfr, &fr, &dfbr, &fbr, + (l > 0) ? pa_order_threshold : 99, skip_zeros); } if (!fi.empty() || !dfbi.empty() || !fbi.empty()) { - bi = (s == 0) ? BuildOperator(fespace_l, (SumCoefficient *)nullptr, &fi, &dfbi, - &fbi, (l > 0) ? pa_order_threshold : 99, skip_zeros) - : BuildAuxOperator(fespace_l, &fi, &fbi, - (l > 0) ? pa_order_threshold : 99, skip_zeros); + bi = aux ? BuildAuxOperator(fespace_l, &fi, &fbi, (l > 0) ? pa_order_threshold : 99, + skip_zeros) + : BuildOperator(fespace_l, (SumCoefficient *)nullptr, &fi, &dfbi, &fbi, + (l > 0) ? pa_order_threshold : 99, skip_zeros); } if (print_prec_hdr) { @@ -718,13 +709,13 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub } auto B_l = BuildLevelOperator(*B, std::move(br), std::move(bi), fespace_l); B_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE); - if (s == 0) + if (aux) { - B->AddOperator(std::move(B_l)); + B->AddAuxiliaryOperator(std::move(B_l)); } else { - B->AddAuxiliaryOperator(std::move(B_l)); + B->AddOperator(std::move(B_l)); } } } @@ -732,69 +723,6 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub return B; } -namespace -{ - -auto BuildCurl(const mfem::ParFiniteElementSpace &nd_fespace, - const mfem::ParFiniteElementSpace &rt_fespace, int pa_order_threshold, - bool skip_zeros) -{ - DiscreteLinearOperator curl(nd_fespace, rt_fespace); - curl.AddDomainInterpolator(); - return curl.Assemble(pa_order_threshold, skip_zeros); -} - -auto BuildGrad(const mfem::ParFiniteElementSpace &h1_fespace, - const mfem::ParFiniteElementSpace &nd_fespace, int pa_order_threshold, - bool skip_zeros) -{ - DiscreteLinearOperator grad(h1_fespace, nd_fespace); - grad.AddDomainInterpolator(); - return grad.Assemble(pa_order_threshold, skip_zeros); -} - -} // namespace - -template <> -std::unique_ptr SpaceOperator::GetCurlMatrix() -{ - constexpr bool skip_zeros_interp = true; - return std::make_unique( - BuildCurl(GetNDSpace(), GetRTSpace(), pa_discrete_interp ? pa_order_threshold : 99, - skip_zeros_interp), - GetNDSpace(), GetRTSpace(), true); -} - -template <> -std::unique_ptr SpaceOperator::GetCurlMatrix() -{ - constexpr bool skip_zeros_interp = true; - return std::make_unique( - BuildCurl(GetNDSpace(), GetRTSpace(), pa_discrete_interp ? pa_order_threshold : 99, - skip_zeros_interp), - nullptr, GetNDSpace(), GetRTSpace(), true); -} - -template <> -std::unique_ptr SpaceOperator::GetGradMatrix() -{ - constexpr bool skip_zeros_interp = true; - return std::make_unique( - BuildGrad(GetH1Space(), GetNDSpace(), pa_discrete_interp ? pa_order_threshold : 99, - skip_zeros_interp), - GetH1Space(), GetNDSpace(), true); -} - -template <> -std::unique_ptr SpaceOperator::GetGradMatrix() -{ - constexpr bool skip_zeros_interp = true; - return std::make_unique( - BuildGrad(GetH1Space(), GetNDSpace(), pa_discrete_interp ? pa_order_threshold : 99, - skip_zeros_interp), - nullptr, GetH1Space(), GetNDSpace(), true); -} - void SpaceOperator::AddStiffnessCoefficients(double coef, SumMatrixCoefficient &df, SumMatrixCoefficient &f) { diff --git a/palace/models/spaceoperator.hpp b/palace/models/spaceoperator.hpp index 3e5c817f9..dd6ced2f5 100644 --- a/palace/models/spaceoperator.hpp +++ b/palace/models/spaceoperator.hpp @@ -8,7 +8,7 @@ #include #include #include -#include "fem/coefficient.hpp" +#include "fem/fespace.hpp" #include "linalg/operator.hpp" #include "linalg/vector.hpp" #include "models/farfieldboundaryoperator.hpp" @@ -23,6 +23,8 @@ namespace palace { class IoData; +class SumCoefficient; +class SumMatrixCoefficient; // // A class handling spatial discretization of the governing equations. @@ -50,8 +52,9 @@ class SpaceOperator std::vector> nd_fecs; std::vector> h1_fecs; std::unique_ptr rt_fec; - std::unique_ptr nd_fespaces, h1_fespaces; - std::unique_ptr rt_fespace; + FiniteElementSpaceHierarchy nd_fespaces; + AuxiliaryFiniteElementSpaceHierarchy h1_fespaces; + AuxiliaryFiniteElementSpace rt_fespace; // Operator for domain material properties. MaterialOperator mat_op; @@ -117,14 +120,16 @@ class SpaceOperator const auto &GetSurfaceCurrentOp() const { return surf_j_op; } // Return the parallel finite element space objects. - auto &GetNDSpaces() { return *nd_fespaces; } - auto &GetNDSpace() { return nd_fespaces->GetFinestFESpace(); } - const auto &GetNDSpace() const { return nd_fespaces->GetFinestFESpace(); } - auto &GetH1Spaces() { return *h1_fespaces; } - auto &GetH1Space() { return h1_fespaces->GetFinestFESpace(); } - const auto &GetH1Space() const { return h1_fespaces->GetFinestFESpace(); } - auto &GetRTSpace() { return *rt_fespace; } - const auto &GetRTSpace() const { return *rt_fespace; } + auto &GetNDSpaces() { return nd_fespaces; } + const auto &GetNDSpaces() const { return nd_fespaces; } + auto &GetNDSpace() { return nd_fespaces.GetFinestFESpace(); } + const auto &GetNDSpace() const { return nd_fespaces.GetFinestFESpace(); } + auto &GetH1Spaces() { return h1_fespaces; } + const auto &GetH1Spaces() const { return h1_fespaces; } + auto &GetH1Space() { return h1_fespaces.GetFinestFESpace(); } + const auto &GetH1Space() const { return h1_fespaces.GetFinestFESpace(); } + auto &GetRTSpace() { return rt_fespace; } + const auto &GetRTSpace() const { return rt_fespace; } // Return the number of true (conforming) dofs on the finest ND space. auto GlobalTrueVSize() { return GetNDSpace().GlobalTrueVSize(); } @@ -171,12 +176,12 @@ class SpaceOperator std::unique_ptr GetPreconditionerMatrix(double a0, double a1, double a2, double a3); - // Construct and return the discrete curl or gradient matrices. The complex variants - // return a matrix suitable for applying to complex-valued vectors. - template - std::unique_ptr GetCurlMatrix(); - template - std::unique_ptr GetGradMatrix(); + // Construct and return the discrete curl or gradient matrices. + const Operator &GetGradMatrix() const + { + return GetH1Spaces().GetFinestAuxiliaryFESpace().GetDiscreteInterpolator(); + } + const Operator &GetCurlMatrix() const { return GetRTSpace().GetDiscreteInterpolator(); } // Assemble the right-hand side source term vector for an incident field or current source // applied on specified excited boundaries. The return value indicates whether or not the diff --git a/palace/models/timeoperator.cpp b/palace/models/timeoperator.cpp index ac7c77b04..741b4ef97 100644 --- a/palace/models/timeoperator.cpp +++ b/palace/models/timeoperator.cpp @@ -143,7 +143,7 @@ TimeOperator::TimeOperator(const IoData &iodata, SpaceOperator &spaceop, std::function &djcoef) { // Construct discrete curl matrix for B-field time integration. - Curl = spaceop.GetCurlMatrix(); + Curl = &spaceop.GetCurlMatrix(); // Allocate space for solution vectors. E.SetSize(Curl->Width()); diff --git a/palace/models/timeoperator.hpp b/palace/models/timeoperator.hpp index 5373b3f2a..f43cb274d 100644 --- a/palace/models/timeoperator.hpp +++ b/palace/models/timeoperator.hpp @@ -32,8 +32,8 @@ class TimeOperator // Time-dependent operator for the E-field. std::unique_ptr op; - // Discrete curl for B-field time integration. - std::unique_ptr Curl; + // Discrete curl for B-field time integration (not owned). + const Operator *Curl; public: TimeOperator(const IoData &iodata, SpaceOperator &spaceop, diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 471d53bca..f245b2da0 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -1642,7 +1642,6 @@ void LinearSolverData::SetUp(json &solver) // Options related to multigrid. mg_max_levels = linear->value("MGMaxLevels", mg_max_levels); mg_coarsen_type = linear->value("MGCoarsenType", mg_coarsen_type); - mg_legacy_transfer = linear->value("MGLegacyTransfer", mg_legacy_transfer); mg_cycle_it = linear->value("MGCycleIts", mg_cycle_it); mg_smooth_aux = linear->value("MGAuxiliarySmoother", mg_smooth_aux); mg_smooth_it = linear->value("MGSmoothIts", mg_smooth_it); @@ -1682,7 +1681,6 @@ void LinearSolverData::SetUp(json &solver) linear->erase("MGMaxLevels"); linear->erase("MGCoarsenType"); - linear->erase("MGLegacyTransfer"); linear->erase("MGCycleIts"); linear->erase("MGAuxiliarySmoother"); linear->erase("MGSmoothIts"); @@ -1721,7 +1719,6 @@ void LinearSolverData::SetUp(json &solver) // std::cout << "MGMaxLevels: " << mg_max_levels << '\n'; // std::cout << "MGCoarsenType: " << mg_coarsen_type << '\n'; - // std::cout << "MGLegacyTransfer: " << mg_legacy_transfer << '\n'; // std::cout << "MGCycleIts: " << mg_cycle_it << '\n'; // std::cout << "MGAuxiliarySmoother: " << mg_smooth_aux << '\n'; // std::cout << "MGSmoothIts: " << mg_smooth_it << '\n'; diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index a7d8a4e55..beff580be 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -730,10 +730,6 @@ struct LinearSolverData }; MultigridCoarsenType mg_coarsen_type = MultigridCoarsenType::LOGARITHMIC; - // Switch to use mfem::TransferOperator or enable partial assembly for the multigrid - // transfer operators. - bool mg_legacy_transfer = false; - // Number of iterations for preconditioners which support it. For multigrid, this is the // number of V-cycles per Krylov solver iteration. int mg_cycle_it = 1; diff --git a/scripts/schema/config/solver.json b/scripts/schema/config/solver.json index 3d3e501b9..4981213e7 100644 --- a/scripts/schema/config/solver.json +++ b/scripts/schema/config/solver.json @@ -107,7 +107,6 @@ "InitialGuess": { "type": "boolean" }, "MGMaxLevels": { "type": "int", "minimum": 1 }, "MGCoarsenType": { "type": "string" }, - "MGLegacyTransfer": { "type": "boolean" }, "MGAuxiliarySmoother": { "type": "boolean" }, "MGCycleIts": { "type": "integer", "exclusiveMinimum": 0 }, "MGSmoothIts": { "type": "integer", "exclusiveMinimum": 0 }, diff --git a/test/unit/test-libceed.cpp b/test/unit/test-libceed.cpp index 875f7a479..5f847eb54 100644 --- a/test/unit/test-libceed.cpp +++ b/test/unit/test-libceed.cpp @@ -366,7 +366,7 @@ void BenchmarkCeedInterpolator(FiniteElementSpace &trial_fespace, FiniteElementSpace &test_fespace, T1 AssembleTest, T2 AssembleRef) { - const bool skip_zeros = true; + const bool skip_zeros_interp = true; Vector x(trial_fespace.GetVSize()), y_ref(test_fespace.GetVSize()), y_test(test_fespace.GetVSize()); x.UseDevice(true); @@ -379,9 +379,9 @@ void BenchmarkCeedInterpolator(FiniteElementSpace &trial_fespace, if (!benchmark_no_fa) { const auto op_test = AssembleTest(trial_fespace, test_fespace); - const auto op_ref = - AssembleRef(trial_fespace, test_fespace, mfem::AssemblyLevel::LEGACY, skip_zeros); - const auto mat_test = DiscreteLinearOperator::FullAssemble(*op_test, skip_zeros); + const auto op_ref = AssembleRef(trial_fespace, test_fespace, + mfem::AssemblyLevel::LEGACY, skip_zeros_interp); + const auto mat_test = DiscreteLinearOperator::FullAssemble(*op_test, skip_zeros_interp); const auto *mat_ref = &op_ref->SpMat(); nnz = mat_test->NumNonZeroElems(); TestCeedOperatorFullAssemble(*mat_test, *mat_ref); @@ -392,13 +392,13 @@ void BenchmarkCeedInterpolator(FiniteElementSpace &trial_fespace, { BENCHMARK("Assemble (MFEM Legacy)") { - const auto op_ref = - AssembleRef(trial_fespace, test_fespace, mfem::AssemblyLevel::LEGACY, skip_zeros); + const auto op_ref = AssembleRef(trial_fespace, test_fespace, + mfem::AssemblyLevel::LEGACY, skip_zeros_interp); return op_ref->Height(); }; { - const auto op_ref = - AssembleRef(trial_fespace, test_fespace, mfem::AssemblyLevel::LEGACY, skip_zeros); + const auto op_ref = AssembleRef(trial_fespace, test_fespace, + mfem::AssemblyLevel::LEGACY, skip_zeros_interp); y_ref = 0.0; BENCHMARK("AddMult (MFEM Legacy)") { @@ -414,12 +414,12 @@ void BenchmarkCeedInterpolator(FiniteElementSpace &trial_fespace, BENCHMARK("Assemble (MFEM Partial)") { const auto op_ref = AssembleRef(trial_fespace, test_fespace, - mfem::AssemblyLevel::PARTIAL, skip_zeros); + mfem::AssemblyLevel::PARTIAL, skip_zeros_interp); return op_ref->Height(); }; { const auto op_ref = AssembleRef(trial_fespace, test_fespace, - mfem::AssemblyLevel::PARTIAL, skip_zeros); + mfem::AssemblyLevel::PARTIAL, skip_zeros_interp); y_ref = 0.0; BENCHMARK("AddMult (MFEM Partial)") { @@ -449,7 +449,8 @@ void BenchmarkCeedInterpolator(FiniteElementSpace &trial_fespace, BENCHMARK("Full Assemble (libCEED)") { const auto op_test = AssembleTest(trial_fespace, test_fespace); - const auto mat_test = DiscreteLinearOperator::FullAssemble(*op_test, skip_zeros); + const auto mat_test = + DiscreteLinearOperator::FullAssemble(*op_test, skip_zeros_interp); return mat_test->NumNonZeroElems(); }; } From 89b567df6475098f26cd71573da5a91bccf748cd Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 23 Oct 2023 12:26:46 -0700 Subject: [PATCH 5/9] Always partially assemble mass matrices for domain postprocessing --- palace/fem/libceed/operator.cpp | 2 +- palace/models/domainpostoperator.cpp | 12 +++++------- palace/models/domainpostoperator.hpp | 2 +- palace/models/postoperator.cpp | 8 +++----- 4 files changed, 10 insertions(+), 14 deletions(-) diff --git a/palace/fem/libceed/operator.cpp b/palace/fem/libceed/operator.cpp index d30e64e30..bfdb7ad62 100644 --- a/palace/fem/libceed/operator.cpp +++ b/palace/fem/libceed/operator.cpp @@ -349,7 +349,7 @@ void CeedOperatorAssembleCOO(const Operator &op, bool skip_zeros, CeedSize *nnz, PalaceCeedCall(ceed, CeedVectorRestoreArray(*vals, &vals_array)); } - // std::cout << " Operator full assembly has " << *nnz << " NNZ"; + // std::cout << " Operator full assembly (COO) has " << *nnz << " NNZ"; if (skip_zeros && *nnz > 0) { CeedOperatorAssembleCOORemoveZeros(ceed, nnz, rows, cols, vals, mem); diff --git a/palace/models/domainpostoperator.cpp b/palace/models/domainpostoperator.cpp index fdac89a26..9cc9b2a36 100644 --- a/palace/models/domainpostoperator.cpp +++ b/palace/models/domainpostoperator.cpp @@ -15,10 +15,9 @@ namespace palace DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOperator &mat_op, const mfem::ParFiniteElementSpace *nd_fespace, - const mfem::ParFiniteElementSpace *rt_fespace, - int pa_order_threshold) + const mfem::ParFiniteElementSpace *rt_fespace) { - constexpr bool skip_zeros = false; + // Mass operators are always partially assembled. if (nd_fespace) { // Construct ND mass matrix to compute the electric field energy integral as: @@ -30,7 +29,7 @@ DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOpera MaterialPropertyCoefficient epsilon_func(mat_op); BilinearForm m_nd(*nd_fespace); m_nd.AddDomainIntegrator(epsilon_func); - M_ND = m_nd.Assemble(pa_order_threshold, skip_zeros); + M_ND = m_nd.Assemble(); D.SetSize(M_ND->Height()); // Use the provided domain postprocessing indices to group for postprocessing bulk @@ -55,8 +54,7 @@ DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOpera BilinearForm mr_nd(*nd_fespace), mi_nd(*nd_fespace); mr_nd.AddDomainIntegrator(epsilon_func_r); mi_nd.AddDomainIntegrator(epsilon_func_i); - M_NDi.emplace(idx, std::make_pair(mr_nd.Assemble(pa_order_threshold, skip_zeros), - mi_nd.Assemble(pa_order_threshold, skip_zeros))); + M_NDi.emplace(idx, std::make_pair(mr_nd.Assemble(), mi_nd.Assemble())); } } @@ -68,7 +66,7 @@ DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOpera MaterialPropertyCoefficient muinv_func(mat_op); BilinearForm m_rt(*rt_fespace); m_rt.AddDomainIntegrator(muinv_func); - M_RT = m_rt.Assemble(pa_order_threshold - 1, skip_zeros); + M_RT = m_rt.Assemble(); H.SetSize(M_RT->Height()); } } diff --git a/palace/models/domainpostoperator.hpp b/palace/models/domainpostoperator.hpp index 9c0e6ab2f..7c31014cb 100644 --- a/palace/models/domainpostoperator.hpp +++ b/palace/models/domainpostoperator.hpp @@ -32,7 +32,7 @@ class DomainPostOperator public: DomainPostOperator(const IoData &iodata, const MaterialOperator &mat_op, const mfem::ParFiniteElementSpace *nd_fespace, - const mfem::ParFiniteElementSpace *rt_fespace, int pa_order_threshold); + const mfem::ParFiniteElementSpace *rt_fespace); // Access underlying bulk loss postprocessing data structures (for keys). const auto &GetEps() const { return M_NDi; } diff --git a/palace/models/postoperator.cpp b/palace/models/postoperator.cpp index bfa1c10c3..5da999346 100644 --- a/palace/models/postoperator.cpp +++ b/palace/models/postoperator.cpp @@ -42,7 +42,7 @@ PostOperator::PostOperator(const IoData &iodata, SpaceOperator &spaceop, : mat_op(spaceop.GetMaterialOp()), surf_post_op(iodata, spaceop.GetMaterialOp(), spaceop.GetH1Space()), dom_post_op(iodata, spaceop.GetMaterialOp(), &spaceop.GetNDSpace(), - &spaceop.GetRTSpace(), iodata.solver.pa_order_threshold), + &spaceop.GetRTSpace()), has_imaginary(iodata.problem.type != config::ProblemData::Type::TRANSIENT), E(&spaceop.GetNDSpace()), B(&spaceop.GetRTSpace()), V(std::nullopt), A(std::nullopt), lumped_port_init(false), wave_port_init(false), @@ -96,8 +96,7 @@ PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplaceop, const std::string &name) : mat_op(laplaceop.GetMaterialOp()), surf_post_op(iodata, laplaceop.GetMaterialOp(), laplaceop.GetH1Space()), - dom_post_op(iodata, laplaceop.GetMaterialOp(), &laplaceop.GetNDSpace(), nullptr, - iodata.solver.pa_order_threshold), + dom_post_op(iodata, laplaceop.GetMaterialOp(), &laplaceop.GetNDSpace(), nullptr), has_imaginary(false), E(&laplaceop.GetNDSpace()), B(std::nullopt), V(&laplaceop.GetH1Space()), A(std::nullopt), lumped_port_init(false), wave_port_init(false), @@ -124,8 +123,7 @@ PostOperator::PostOperator(const IoData &iodata, CurlCurlOperator &curlcurlop, const std::string &name) : mat_op(curlcurlop.GetMaterialOp()), surf_post_op(iodata, curlcurlop.GetMaterialOp(), curlcurlop.GetH1Space()), - dom_post_op(iodata, curlcurlop.GetMaterialOp(), nullptr, &curlcurlop.GetRTSpace(), - iodata.solver.pa_order_threshold), + dom_post_op(iodata, curlcurlop.GetMaterialOp(), nullptr, &curlcurlop.GetRTSpace()), has_imaginary(false), E(std::nullopt), B(&curlcurlop.GetRTSpace()), V(std::nullopt), A(&curlcurlop.GetNDSpace()), lumped_port_init(false), wave_port_init(false), paraview(CreateParaviewPath(iodata, name), curlcurlop.GetNDSpace().GetParMesh()), From 7df6a19cd2dfda85a1821beaabc814b79952de3c Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Fri, 27 Oct 2023 11:47:02 -0400 Subject: [PATCH 6/9] Cherry-pick hughcars/discrete-interpolator-refactor: Separate the space hierarchies --- palace/fem/fespace.cpp | 6 ++- palace/fem/fespace.hpp | 75 ++++++++++++--------------------- palace/linalg/divfree.cpp | 2 +- palace/linalg/hcurl.cpp | 13 +++--- palace/linalg/ksp.cpp | 6 +-- palace/models/spaceoperator.cpp | 9 ++-- palace/models/spaceoperator.hpp | 2 +- 7 files changed, 50 insertions(+), 63 deletions(-) diff --git a/palace/fem/fespace.cpp b/palace/fem/fespace.cpp index 9617f1436..ff841f672 100644 --- a/palace/fem/fespace.cpp +++ b/palace/fem/fespace.cpp @@ -90,7 +90,8 @@ const Operator &AuxiliaryFiniteElementSpace::BuildDiscreteInterpolator() const return *G; } -const Operator &FiniteElementSpaceHierarchy::BuildProlongationAtLevel(std::size_t l) const +template +const Operator &Hierarchy::BuildProlongationAtLevel(std::size_t l) const { // P is always partially assembled. MFEM_VERIFY(l >= 0 && l < GetNumLevels() - 1, @@ -113,4 +114,7 @@ const Operator &FiniteElementSpaceHierarchy::BuildProlongationAtLevel(std::size_ return *P[l]; } +template class Hierarchy; +template class Hierarchy; + } // namespace palace diff --git a/palace/fem/fespace.hpp b/palace/fem/fespace.hpp index 82748938d..ff03306c9 100644 --- a/palace/fem/fespace.hpp +++ b/palace/fem/fespace.hpp @@ -70,50 +70,51 @@ class AuxiliaryFiniteElementSpace : public FiniteElementSpace // A collection of FiniteElementSpace objects constructed on the same mesh with the ability // to construct the prolongation operators between them as needed. // -class FiniteElementSpaceHierarchy +template +class Hierarchy { protected: - std::vector> fespaces; + std::vector> fespaces; mutable std::vector> P; + static_assert(std::is_convertible_v, + "A hierarchy can only be constructed of FiniteElementSpaces"); + const Operator &BuildProlongationAtLevel(std::size_t l) const; public: - FiniteElementSpaceHierarchy() = default; - FiniteElementSpaceHierarchy(std::unique_ptr &&fespace) - { - AddLevel(std::move(fespace)); - } + Hierarchy() = default; + Hierarchy(std::unique_ptr &&fespace) { AddLevel(std::move(fespace)); } auto GetNumLevels() const { return fespaces.size(); } auto size() const { return GetNumLevels(); } bool empty() const { return GetNumLevels() == 0; } - virtual void AddLevel(std::unique_ptr &&fespace) + virtual void AddLevel(std::unique_ptr &&fespace) { fespaces.push_back(std::move(fespace)); P.push_back(nullptr); } - FiniteElementSpace &GetFESpaceAtLevel(std::size_t l) + FESpace &GetFESpaceAtLevel(std::size_t l) { MFEM_ASSERT(l >= 0 && l < GetNumLevels(), "Out of bounds request for finite element space at level " << l << "!"); return *fespaces[l]; } - const FiniteElementSpace &GetFESpaceAtLevel(std::size_t l) const + const FESpace &GetFESpaceAtLevel(std::size_t l) const { MFEM_ASSERT(l >= 0 && l < GetNumLevels(), "Out of bounds request for finite element space at level " << l << "!"); return *fespaces[l]; } - FiniteElementSpace &GetFinestFESpace() + FESpace &GetFinestFESpace() { MFEM_ASSERT(!empty(), "Out of bounds request for finite element space at level 0!"); return *fespaces.back(); } - const FiniteElementSpace &GetFinestFESpace() const + const FESpace &GetFinestFESpace() const { MFEM_ASSERT(!empty(), "Out of bounds request for finite element space at level 0!"); return *fespaces.back(); @@ -138,54 +139,30 @@ class FiniteElementSpaceHierarchy } }; +class FiniteElementSpaceHierarchy : public Hierarchy +{ +public: + using Hierarchy::Hierarchy; +}; + // // A special type of FiniteElementSpaceHierarchy where all members are auxiliary finite // element spaces. // -class AuxiliaryFiniteElementSpaceHierarchy : public FiniteElementSpaceHierarchy +class AuxiliaryFiniteElementSpaceHierarchy : public Hierarchy { -public: - AuxiliaryFiniteElementSpaceHierarchy() = default; - AuxiliaryFiniteElementSpaceHierarchy( - std::unique_ptr &&fespace) - { - AddLevel(std::move(fespace)); - } - void AddLevel(std::unique_ptr &&fespace) override - { - MFEM_ABORT("All finite element spaces in an AuxiliaryFiniteElementSpaceHierarchy must " - "inherit from AuxiliaryFiniteElementSpace!"); - } - - void AddLevel(std::unique_ptr &&fespace) - { - // Guarantees that every object in fespaces is an AuxiliaryFiniteElementSpace. - fespaces.push_back(std::move(fespace)); - P.push_back(nullptr); - } - - AuxiliaryFiniteElementSpace &GetAuxiliaryFESpaceAtLevel(std::size_t l) - { - return *static_cast(&GetFESpaceAtLevel(l)); - } - const AuxiliaryFiniteElementSpace &GetAuxiliaryFESpaceAtLevel(std::size_t l) const - { - return *static_cast(&GetFESpaceAtLevel(l)); - } +public: + using Hierarchy::Hierarchy; - AuxiliaryFiniteElementSpace &GetFinestAuxiliaryFESpace() - { - return *static_cast(&GetFinestFESpace()); - } - const AuxiliaryFiniteElementSpace &GetFinestAuxiliaryFESpace() const + const Operator &GetDiscreteInterpolatorAtLevel(std::size_t l) const { - return *static_cast(&GetFinestFESpace()); + return GetFESpaceAtLevel(l).GetDiscreteInterpolator(); } - const Operator &GetDiscreteInterpolatorAtLevel(std::size_t l) const + const Operator &GetDiscreteInterpolatorAtFinestLevel() const { - return GetAuxiliaryFESpaceAtLevel(l).GetDiscreteInterpolator(); + return GetFinestFESpace().GetDiscreteInterpolator(); } std::vector GetDiscreteInterpolators() const diff --git a/palace/linalg/divfree.cpp b/palace/linalg/divfree.cpp index 4266bb856..988441cf1 100644 --- a/palace/linalg/divfree.cpp +++ b/palace/linalg/divfree.cpp @@ -49,7 +49,7 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, std::make_unique(weakdiv.Assemble(pa_order_threshold, skip_zeros), nd_fespace, h1_fespaces.GetFinestFESpace(), false); } - Grad = &h1_fespaces.GetFinestAuxiliaryFESpace().GetDiscreteInterpolator(); + Grad = &h1_fespaces.GetDiscreteInterpolatorAtFinestLevel(); bdr_tdof_list_M = &h1_bdr_tdof_lists.back(); // The system matrix for the projection is real and SPD. diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp index a35fd9e99..6d3bb5762 100644 --- a/palace/linalg/hcurl.cpp +++ b/palace/linalg/hcurl.cpp @@ -31,14 +31,17 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( MaterialPropertyCoefficient epsilon_func(mat_op); { auto A_mg = std::make_unique(nd_fespaces.GetNumLevels()); + MFEM_VERIFY(h1_fespaces.GetNumLevels() == nd_fespaces.GetNumLevels(), "!"); + auto num_levels = h1_fespaces.GetNumLevels(); for (bool aux : {false, true}) { - const auto &fespaces = aux ? h1_fespaces : nd_fespaces; + // const auto &fespaces = aux ? h1_fespaces : nd_fespaces; const auto &dbc_tdof_lists = aux ? h1_dbc_tdof_lists : nd_dbc_tdof_lists; - for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++) + for (std::size_t l = 0; l < num_levels; l++) { // Force coarse level operator to be fully assembled always. - const auto &fespace_l = fespaces.GetFESpaceAtLevel(l); + const auto &fespace_l = + aux ? h1_fespaces.GetFESpaceAtLevel(l) : nd_fespaces.GetFESpaceAtLevel(l); BilinearForm a(fespace_l); if (aux) { @@ -67,8 +70,8 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( // The system matrix K + M is real and SPD. We use Hypre's AMS solver as the coarse-level // multigrid solve. auto ams = std::make_unique>(std::make_unique( - nd_fespaces.GetFESpaceAtLevel(0), h1_fespaces.GetAuxiliaryFESpaceAtLevel(0), 1, 1, 1, - false, false, 0)); + nd_fespaces.GetFESpaceAtLevel(0), h1_fespaces.GetFESpaceAtLevel(0), 1, 1, 1, false, + false, 0)); std::unique_ptr> pc; if (nd_fespaces.GetNumLevels() > 1) { diff --git a/palace/linalg/ksp.cpp b/palace/linalg/ksp.cpp index 8414adb47..0c59949be 100644 --- a/palace/linalg/ksp.cpp +++ b/palace/linalg/ksp.cpp @@ -121,9 +121,9 @@ ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata, // space (in which case fespaces.GetNumLevels() == 1). MFEM_VERIFY(aux_fespaces, "AMS solver relies on both primary space " "and auxiliary spaces for construction!"); - pc0 = std::make_unique( - iodata, fespaces.GetNumLevels() > 1, fespaces.GetFESpaceAtLevel(0), - aux_fespaces->GetAuxiliaryFESpaceAtLevel(0), print); + pc0 = std::make_unique(iodata, fespaces.GetNumLevels() > 1, + fespaces.GetFESpaceAtLevel(0), + aux_fespaces->GetFESpaceAtLevel(0), print); break; case config::LinearSolverData::Type::BOOMER_AMG: pc0 = std::make_unique(iodata, fespaces.GetNumLevels() > 1, print); diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index cb3c22292..d2b80c760 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -639,14 +639,17 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub MFEM_VERIFY(GetH1Spaces().GetNumLevels() == GetNDSpaces().GetNumLevels(), "Multigrid hierarchy mismatch for auxiliary space preconditioning!"); auto B = std::make_unique>(GetNDSpaces().GetNumLevels()); + MFEM_VERIFY(GetH1Spaces().GetNumLevels() == GetNDSpaces().GetNumLevels(), "!"); + auto num_levels = GetNDSpaces().GetNumLevels(); for (bool aux : {false, true}) { - auto &fespaces = aux ? GetH1Spaces() : GetNDSpaces(); + // auto &fespaces = aux ? GetH1Spaces() : GetNDSpaces(); auto &dbc_tdof_lists = aux ? h1_dbc_tdof_lists : nd_dbc_tdof_lists; - for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++) + for (std::size_t l = 0; l < num_levels; l++) { // Force coarse level operator to be fully assembled always. - auto &fespace_l = fespaces.GetFESpaceAtLevel(l); + auto &fespace_l = + aux ? GetH1Spaces().GetFESpaceAtLevel(l) : GetNDSpaces().GetFESpaceAtLevel(l); if (print_prec_hdr) { Mpi::Print(" Level {:d}{} (p = {:d}): {:d} unknowns", l, aux ? " (auxiliary)" : "", diff --git a/palace/models/spaceoperator.hpp b/palace/models/spaceoperator.hpp index dd6ced2f5..4f5a737a1 100644 --- a/palace/models/spaceoperator.hpp +++ b/palace/models/spaceoperator.hpp @@ -179,7 +179,7 @@ class SpaceOperator // Construct and return the discrete curl or gradient matrices. const Operator &GetGradMatrix() const { - return GetH1Spaces().GetFinestAuxiliaryFESpace().GetDiscreteInterpolator(); + return GetH1Spaces().GetDiscreteInterpolatorAtFinestLevel(); } const Operator &GetCurlMatrix() const { return GetRTSpace().GetDiscreteInterpolator(); } From ff69f5e9a085305d20aec90e30838a7141344ceb Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 27 Oct 2023 10:17:54 -0700 Subject: [PATCH 7/9] Some amendments to PR feedback in prior commit --- palace/fem/fespace.cpp | 7 ++--- palace/fem/fespace.hpp | 45 +++++++++++++++++---------------- palace/fem/multigrid.hpp | 6 ++--- palace/linalg/divfree.cpp | 2 +- palace/linalg/hcurl.cpp | 14 +++++----- palace/models/spaceoperator.cpp | 14 +++++----- palace/models/spaceoperator.hpp | 2 +- 7 files changed, 45 insertions(+), 45 deletions(-) diff --git a/palace/fem/fespace.cpp b/palace/fem/fespace.cpp index ff841f672..98cd0a44d 100644 --- a/palace/fem/fespace.cpp +++ b/palace/fem/fespace.cpp @@ -91,7 +91,8 @@ const Operator &AuxiliaryFiniteElementSpace::BuildDiscreteInterpolator() const } template -const Operator &Hierarchy::BuildProlongationAtLevel(std::size_t l) const +const Operator & +BaseFiniteElementSpaceHierarchy::BuildProlongationAtLevel(std::size_t l) const { // P is always partially assembled. MFEM_VERIFY(l >= 0 && l < GetNumLevels() - 1, @@ -114,7 +115,7 @@ const Operator &Hierarchy::BuildProlongationAtLevel(std::size_t l) cons return *P[l]; } -template class Hierarchy; -template class Hierarchy; +template class BaseFiniteElementSpaceHierarchy; +template class BaseFiniteElementSpaceHierarchy; } // namespace palace diff --git a/palace/fem/fespace.hpp b/palace/fem/fespace.hpp index ff03306c9..05d495e89 100644 --- a/palace/fem/fespace.hpp +++ b/palace/fem/fespace.hpp @@ -39,7 +39,7 @@ class FiniteElementSpace : public mfem::ParFiniteElementSpace }; // -// An AuxiliaryFiniteElement space is just a FiniteElementSpace which allows for lazy +// An AuxiliaryFiniteElement space is a FiniteElementSpace which allows for lazy // construction of the interpolation operator (discrete gradient or curl) from the primal // space to this one. // @@ -71,26 +71,27 @@ class AuxiliaryFiniteElementSpace : public FiniteElementSpace // to construct the prolongation operators between them as needed. // template -class Hierarchy +class BaseFiniteElementSpaceHierarchy { + static_assert(std::is_base_of::value, + "A space hierarchy can only be constructed of FiniteElementSpace objects!"); + protected: std::vector> fespaces; mutable std::vector> P; - static_assert(std::is_convertible_v, - "A hierarchy can only be constructed of FiniteElementSpaces"); - const Operator &BuildProlongationAtLevel(std::size_t l) const; public: - Hierarchy() = default; - Hierarchy(std::unique_ptr &&fespace) { AddLevel(std::move(fespace)); } + BaseFiniteElementSpaceHierarchy() = default; + BaseFiniteElementSpaceHierarchy(std::unique_ptr &&fespace) + { + AddLevel(std::move(fespace)); + } auto GetNumLevels() const { return fespaces.size(); } - auto size() const { return GetNumLevels(); } - bool empty() const { return GetNumLevels() == 0; } - virtual void AddLevel(std::unique_ptr &&fespace) + void AddLevel(std::unique_ptr &&fespace) { fespaces.push_back(std::move(fespace)); P.push_back(nullptr); @@ -111,12 +112,14 @@ class Hierarchy FESpace &GetFinestFESpace() { - MFEM_ASSERT(!empty(), "Out of bounds request for finite element space at level 0!"); + MFEM_ASSERT(GetNumLevels() > 0, + "Out of bounds request for finite element space at level 0!"); return *fespaces.back(); } const FESpace &GetFinestFESpace() const { - MFEM_ASSERT(!empty(), "Out of bounds request for finite element space at level 0!"); + MFEM_ASSERT(GetNumLevels() > 0, + "Out of bounds request for finite element space at level 0!"); return *fespaces.back(); } @@ -139,32 +142,30 @@ class Hierarchy } }; -class FiniteElementSpaceHierarchy : public Hierarchy +class FiniteElementSpaceHierarchy + : public BaseFiniteElementSpaceHierarchy { public: - using Hierarchy::Hierarchy; + using BaseFiniteElementSpaceHierarchy< + FiniteElementSpace>::BaseFiniteElementSpaceHierarchy; }; // // A special type of FiniteElementSpaceHierarchy where all members are auxiliary finite // element spaces. // -class AuxiliaryFiniteElementSpaceHierarchy : public Hierarchy +class AuxiliaryFiniteElementSpaceHierarchy + : public BaseFiniteElementSpaceHierarchy { - public: - using Hierarchy::Hierarchy; + using BaseFiniteElementSpaceHierarchy< + AuxiliaryFiniteElementSpace>::BaseFiniteElementSpaceHierarchy; const Operator &GetDiscreteInterpolatorAtLevel(std::size_t l) const { return GetFESpaceAtLevel(l).GetDiscreteInterpolator(); } - const Operator &GetDiscreteInterpolatorAtFinestLevel() const - { - return GetFinestFESpace().GetDiscreteInterpolator(); - } - std::vector GetDiscreteInterpolators() const { std::vector G_(GetNumLevels()); diff --git a/palace/fem/multigrid.hpp b/palace/fem/multigrid.hpp index e8c6fbaef..106420186 100644 --- a/palace/fem/multigrid.hpp +++ b/palace/fem/multigrid.hpp @@ -129,7 +129,7 @@ inline AuxiliaryFiniteElementSpaceHierarchy ConstructAuxiliaryFiniteElementSpace const mfem::Array *dbc_marker = nullptr, std::vector> *dbc_tdof_lists = nullptr) { - MFEM_VERIFY(!primal_fespaces.empty() && !fecs.empty() && + MFEM_VERIFY((primal_fespaces.GetNumLevels() > 0) && !fecs.empty() && (!dbc_tdof_lists || dbc_tdof_lists->empty()), "Empty mesh or FE collection for FE space construction!"); mfem::ParMesh *mesh = primal_fespaces.GetFESpaceAtLevel(0).GetParMesh(); @@ -144,7 +144,7 @@ inline AuxiliaryFiniteElementSpaceHierarchy ConstructAuxiliaryFiniteElementSpace // h-refinement std::size_t l; - for (l = 1; l < primal_fespaces.size(); l++) + for (l = 1; l < primal_fespaces.GetNumLevels(); l++) { if (primal_fespaces.GetFESpaceAtLevel(l).GetParMesh() == mesh) { @@ -163,7 +163,7 @@ inline AuxiliaryFiniteElementSpaceHierarchy ConstructAuxiliaryFiniteElementSpace // p-refinement const auto l0 = l - 1; - for (; l < primal_fespaces.size(); l++) + for (; l < primal_fespaces.GetNumLevels(); l++) { fespaces.AddLevel(std::make_unique( primal_fespaces.GetFESpaceAtLevel(l), mesh, fecs[l - l0].get())); diff --git a/palace/linalg/divfree.cpp b/palace/linalg/divfree.cpp index 988441cf1..e997e06c9 100644 --- a/palace/linalg/divfree.cpp +++ b/palace/linalg/divfree.cpp @@ -49,7 +49,7 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, std::make_unique(weakdiv.Assemble(pa_order_threshold, skip_zeros), nd_fespace, h1_fespaces.GetFinestFESpace(), false); } - Grad = &h1_fespaces.GetDiscreteInterpolatorAtFinestLevel(); + Grad = &h1_fespaces.GetFinestFESpace().GetDiscreteInterpolator(); bdr_tdof_list_M = &h1_bdr_tdof_lists.back(); // The system matrix for the projection is real and SPD. diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp index 6d3bb5762..3a8b1395b 100644 --- a/palace/linalg/hcurl.cpp +++ b/palace/linalg/hcurl.cpp @@ -30,18 +30,18 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( MaterialPropertyCoefficient muinv_func(mat_op); MaterialPropertyCoefficient epsilon_func(mat_op); { - auto A_mg = std::make_unique(nd_fespaces.GetNumLevels()); - MFEM_VERIFY(h1_fespaces.GetNumLevels() == nd_fespaces.GetNumLevels(), "!"); - auto num_levels = h1_fespaces.GetNumLevels(); + MFEM_VERIFY(h1_fespaces.GetNumLevels() == nd_fespaces.GetNumLevels(), + "Multigrid hierarchy mismatch for auxiliary space preconditioning!"); + const auto n_levels = nd_fespaces.GetNumLevels(); + auto A_mg = std::make_unique(n_levels); for (bool aux : {false, true}) { - // const auto &fespaces = aux ? h1_fespaces : nd_fespaces; - const auto &dbc_tdof_lists = aux ? h1_dbc_tdof_lists : nd_dbc_tdof_lists; - for (std::size_t l = 0; l < num_levels; l++) + for (std::size_t l = 0; l < n_levels; l++) { // Force coarse level operator to be fully assembled always. const auto &fespace_l = aux ? h1_fespaces.GetFESpaceAtLevel(l) : nd_fespaces.GetFESpaceAtLevel(l); + const auto &dbc_tdof_lists_l = aux ? h1_dbc_tdof_lists[l] : nd_dbc_tdof_lists[l]; BilinearForm a(fespace_l); if (aux) { @@ -53,7 +53,7 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( } auto A_l = std::make_unique( a.Assemble((l > 0) ? pa_order_threshold : 99, skip_zeros), fespace_l); - A_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE); + A_l->SetEssentialTrueDofs(dbc_tdof_lists_l, Operator::DiagonalPolicy::DIAG_ONE); if (aux) { A_mg->AddAuxiliaryOperator(std::move(A_l)); diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index d2b80c760..15b389ac2 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -638,18 +638,16 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub } MFEM_VERIFY(GetH1Spaces().GetNumLevels() == GetNDSpaces().GetNumLevels(), "Multigrid hierarchy mismatch for auxiliary space preconditioning!"); - auto B = std::make_unique>(GetNDSpaces().GetNumLevels()); - MFEM_VERIFY(GetH1Spaces().GetNumLevels() == GetNDSpaces().GetNumLevels(), "!"); - auto num_levels = GetNDSpaces().GetNumLevels(); + const auto n_levels = GetNDSpaces().GetNumLevels(); + auto B = std::make_unique>(n_levels); for (bool aux : {false, true}) { - // auto &fespaces = aux ? GetH1Spaces() : GetNDSpaces(); - auto &dbc_tdof_lists = aux ? h1_dbc_tdof_lists : nd_dbc_tdof_lists; - for (std::size_t l = 0; l < num_levels; l++) + for (std::size_t l = 0; l < n_levels; l++) { // Force coarse level operator to be fully assembled always. - auto &fespace_l = + const auto &fespace_l = aux ? GetH1Spaces().GetFESpaceAtLevel(l) : GetNDSpaces().GetFESpaceAtLevel(l); + const auto &dbc_tdof_lists_l = aux ? h1_dbc_tdof_lists[l] : nd_dbc_tdof_lists[l]; if (print_prec_hdr) { Mpi::Print(" Level {:d}{} (p = {:d}): {:d} unknowns", l, aux ? " (auxiliary)" : "", @@ -711,7 +709,7 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub } } auto B_l = BuildLevelOperator(*B, std::move(br), std::move(bi), fespace_l); - B_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE); + B_l->SetEssentialTrueDofs(dbc_tdof_lists_l, Operator::DiagonalPolicy::DIAG_ONE); if (aux) { B->AddAuxiliaryOperator(std::move(B_l)); diff --git a/palace/models/spaceoperator.hpp b/palace/models/spaceoperator.hpp index 4f5a737a1..35dabfac2 100644 --- a/palace/models/spaceoperator.hpp +++ b/palace/models/spaceoperator.hpp @@ -179,7 +179,7 @@ class SpaceOperator // Construct and return the discrete curl or gradient matrices. const Operator &GetGradMatrix() const { - return GetH1Spaces().GetDiscreteInterpolatorAtFinestLevel(); + return GetH1Spaces().GetFinestFESpace().GetDiscreteInterpolator(); } const Operator &GetCurlMatrix() const { return GetRTSpace().GetDiscreteInterpolator(); } From 726ca649c504f3492c2c69fada66c439634cba89 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 27 Oct 2023 10:20:37 -0700 Subject: [PATCH 8/9] Fix compiler warnings due to unusued variable --- docs/src/config/solver.md | 4 ---- palace/models/curlcurloperator.cpp | 3 +-- palace/models/curlcurloperator.hpp | 5 ++--- palace/models/laplaceoperator.cpp | 3 +-- palace/models/laplaceoperator.hpp | 5 ++--- palace/models/spaceoperator.cpp | 3 +-- palace/models/spaceoperator.hpp | 9 ++++----- palace/utils/configfile.cpp | 3 --- palace/utils/configfile.hpp | 3 --- scripts/schema/config/solver.json | 1 - 10 files changed, 11 insertions(+), 28 deletions(-) diff --git a/docs/src/config/solver.md b/docs/src/config/solver.md index 686e611d6..4735c612d 100644 --- a/docs/src/config/solver.md +++ b/docs/src/config/solver.md @@ -89,10 +89,6 @@ Thus, this object is only relevant for [`config["Problem"]["Type"]: "Magnetostat `"Linear"` : Top-level object for configuring the linear solver employed by all simulation types. -### Advanced solver options - - - `"PartialAssemblyInterpolators" [true]` - ## `solver["Eigenmode"]` ```json diff --git a/palace/models/curlcurloperator.cpp b/palace/models/curlcurloperator.cpp index de0bc4973..4aba3849a 100644 --- a/palace/models/curlcurloperator.cpp +++ b/palace/models/curlcurloperator.cpp @@ -71,8 +71,7 @@ mfem::Array SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMe CurlCurlOperator::CurlCurlOperator(const IoData &iodata, const std::vector> &mesh) - : pa_order_threshold(iodata.solver.pa_order_threshold), - pa_discrete_interp(iodata.solver.pa_discrete_interp), skip_zeros(false), + : pa_order_threshold(iodata.solver.pa_order_threshold), skip_zeros(false), print_hdr(true), dbc_marker(SetUpBoundaryProperties(iodata, *mesh.back())), nd_fecs(fem::ConstructFECollections( iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels, diff --git a/palace/models/curlcurloperator.hpp b/palace/models/curlcurloperator.hpp index 773f1be07..7393ab0db 100644 --- a/palace/models/curlcurloperator.hpp +++ b/palace/models/curlcurloperator.hpp @@ -24,9 +24,8 @@ class IoData; class CurlCurlOperator { private: - const int pa_order_threshold; // Order above which to use partial assembly vs. full - const bool pa_discrete_interp; // Use partial assembly for discrete interpolators - const bool skip_zeros; // Skip zeros during full assembly of matrices + const int pa_order_threshold; // Order above which to use partial assembly vs. full + const bool skip_zeros; // Skip zeros during full assembly of matrices // Helper variable for log file printing. bool print_hdr; diff --git a/palace/models/laplaceoperator.cpp b/palace/models/laplaceoperator.cpp index 2b648dc78..6cde0f370 100644 --- a/palace/models/laplaceoperator.cpp +++ b/palace/models/laplaceoperator.cpp @@ -115,8 +115,7 @@ std::map> ConstructSources(const IoData &iodata) LaplaceOperator::LaplaceOperator(const IoData &iodata, const std::vector> &mesh) - : pa_order_threshold(iodata.solver.pa_order_threshold), - pa_discrete_interp(iodata.solver.pa_discrete_interp), skip_zeros(false), + : pa_order_threshold(iodata.solver.pa_order_threshold), skip_zeros(false), print_hdr(true), dbc_marker(SetUpBoundaryProperties(iodata, *mesh.back())), h1_fecs(fem::ConstructFECollections( iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels, diff --git a/palace/models/laplaceoperator.hpp b/palace/models/laplaceoperator.hpp index 603c98414..28e71d818 100644 --- a/palace/models/laplaceoperator.hpp +++ b/palace/models/laplaceoperator.hpp @@ -24,9 +24,8 @@ class IoData; class LaplaceOperator { private: - const int pa_order_threshold; // Order above which to use partial assembly vs. full - const bool pa_discrete_interp; // Use partial assembly for discrete interpolators - const bool skip_zeros; // Skip zeros during full assembly of matrices + const int pa_order_threshold; // Order above which to use partial assembly vs. full + const bool skip_zeros; // Skip zeros during full assembly of matrices // Helper variable for log file printing. bool print_hdr; diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index 15b389ac2..ad8a5493d 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -74,8 +74,7 @@ mfem::Array SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMe SpaceOperator::SpaceOperator(const IoData &iodata, const std::vector> &mesh) - : pa_order_threshold(iodata.solver.pa_order_threshold), - pa_discrete_interp(iodata.solver.pa_discrete_interp), skip_zeros(false), + : pa_order_threshold(iodata.solver.pa_order_threshold), skip_zeros(false), pc_mat_real(iodata.solver.linear.pc_mat_real), pc_mat_shifted(iodata.solver.linear.pc_mat_shifted), print_hdr(true), print_prec_hdr(true), dbc_marker(SetUpBoundaryProperties(iodata, *mesh.back())), diff --git a/palace/models/spaceoperator.hpp b/palace/models/spaceoperator.hpp index 35dabfac2..3cce968d5 100644 --- a/palace/models/spaceoperator.hpp +++ b/palace/models/spaceoperator.hpp @@ -32,11 +32,10 @@ class SumMatrixCoefficient; class SpaceOperator { private: - const int pa_order_threshold; // Order above which to use partial assembly vs. full - const bool pa_discrete_interp; // Use partial assembly for discrete interpolators - const bool skip_zeros; // Skip zeros during full assembly of matrices - const bool pc_mat_real; // Use real-valued matrix for preconditioner - const bool pc_mat_shifted; // Use shifted mass matrix for preconditioner + const int pa_order_threshold; // Order above which to use partial assembly vs. full + const bool skip_zeros; // Skip zeros during full assembly of matrices + const bool pc_mat_real; // Use real-valued matrix for preconditioner + const bool pc_mat_shifted; // Use shifted mass matrix for preconditioner // Helper variables for log file printing. bool print_hdr, print_prec_hdr; diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index f245b2da0..03ac4171a 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -1756,7 +1756,6 @@ void SolverData::SetUp(json &config) } order = solver->value("Order", order); pa_order_threshold = solver->value("PartialAssemblyOrder", pa_order_threshold); - pa_discrete_interp = solver->value("PartialAssemblyInterpolators", pa_discrete_interp); device = solver->value("Device", device); ceed_backend = solver->value("Backend", ceed_backend); @@ -1770,7 +1769,6 @@ void SolverData::SetUp(json &config) // Cleanup solver->erase("Order"); solver->erase("PartialAssemblyOrder"); - solver->erase("PartialAssemblyInterpolators"); solver->erase("Device"); solver->erase("Backend"); @@ -1787,7 +1785,6 @@ void SolverData::SetUp(json &config) // Debug // std::cout << "Order: " << order << '\n'; // std::cout << "PartialAssemblyOrder: " << pa_order_threshold << '\n'; - // std::cout << "PartialAssemblyInterpolators: " << pa_discrete_interp << '\n'; // std::cout << "Device: " << device << '\n'; // std::cout << "Backend: " << ceed_backend << '\n'; } diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index beff580be..2907cff3b 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -839,9 +839,6 @@ struct SolverData // Order above which to use partial assembly instead of full assembly. int pa_order_threshold = 100; - // Enable partial assembly for discrete linear operators. - bool pa_discrete_interp = true; - // Device used to configure MFEM. enum class Device { diff --git a/scripts/schema/config/solver.json b/scripts/schema/config/solver.json index 4981213e7..d69609cdb 100644 --- a/scripts/schema/config/solver.json +++ b/scripts/schema/config/solver.json @@ -8,7 +8,6 @@ { "Order": { "type": "integer", "minimum": 1 }, "PartialAssemblyOrder": { "type": "integer", "minimum": 1 }, - "PartialAssemblyInterpolators": { "type": "boolean" }, "Device": { "type": "string", "enum": ["CPU", "GPU", "Debug"] }, "Backend": { "type": "string" }, "Eigenmode": From 8c05eaf6638ac19cd287a7c347987c3f96a8fb54 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 31 Oct 2023 11:10:22 -0700 Subject: [PATCH 9/9] Fix a few missed int skip_zeros from rebase (caught in review) --- palace/linalg/errorestimator.cpp | 9 ++++----- palace/models/spaceoperator.cpp | 4 ++-- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp index 93242ecf4..76d75c403 100644 --- a/palace/linalg/errorestimator.cpp +++ b/palace/linalg/errorestimator.cpp @@ -23,8 +23,9 @@ namespace { std::unique_ptr GetMassMatrix(const FiniteElementSpaceHierarchy &fespaces, - int pa_order_threshold, int skip_zeros) + int pa_order_threshold) { + constexpr bool skip_zeros = false; const int dim = fespaces.GetFinestFESpace().GetParMesh()->Dimension(); const auto type = fespaces.GetFinestFESpace().FEColl()->GetRangeType(dim); auto M = std::make_unique(fespaces.GetNumLevels()); @@ -96,8 +97,7 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, flux.AddDomainIntegrator(muinv_func); Flux = std::make_unique(flux.Assemble(), nd_fespaces.GetFinestFESpace()); } - constexpr int skip_zeros = false; - M = GetMassMatrix(nd_fespaces, pa_order_threshold, skip_zeros); + M = GetMassMatrix(nd_fespaces, pa_order_threshold); ksp = ConfigureLinearSolver(nd_fespaces, tol, max_it, print); ksp->SetOperators(*M, *M); @@ -120,8 +120,7 @@ FluxProjector::FluxProjector(const MaterialOperator &mat_op, Flux = std::make_unique(flux.Assemble(), h1_fespaces.GetFinestFESpace(), h1d_fespace, false); } - constexpr int skip_zeros = false; - M = GetMassMatrix(h1_fespaces, pa_order_threshold, skip_zeros); + M = GetMassMatrix(h1_fespaces, pa_order_threshold); ksp = ConfigureLinearSolver(h1_fespaces, tol, max_it, print); ksp->SetOperators(*M, *M); diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index ad8a5493d..cc2a174e1 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -192,7 +192,7 @@ void PrintHeader(const FiniteElementSpace &h1_fespace, const FiniteElementSpace template std::unique_ptr BuildOperator(const FiniteElementSpace &fespace, T1 *df, T2 *f, T3 *dfb, T4 *fb, int pa_order_threshold, - int skip_zeros) + bool skip_zeros) { BilinearForm a(fespace); if (df && !df->empty() && f && !f->empty()) @@ -230,7 +230,7 @@ std::unique_ptr BuildOperator(const FiniteElementSpace &fespace, T1 *d template std::unique_ptr BuildAuxOperator(const FiniteElementSpace &fespace, T1 *f, T2 *fb, - int pa_order_threshold, int skip_zeros) + int pa_order_threshold, bool skip_zeros) { BilinearForm a(fespace); if (f && !f->empty())