From b18f41f88f54cebc408f40b0650ce59456611a24 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 16:19:22 -0400 Subject: [PATCH] Add fixed-work mass-consistent gravity Schur refinement --- .../GRAVITY_CONSISTENT_SCHUR_2026-09-18.md | 119 +++++++ .../preconditioning/gravity_field.cppm | 121 ++++++- tests/preconditioning/gravity_field.cpp | 298 ++++++++++++++++++ 3 files changed, 531 insertions(+), 7 deletions(-) create mode 100644 experiments/GRAVITY_CONSISTENT_SCHUR_2026-09-18.md diff --git a/experiments/GRAVITY_CONSISTENT_SCHUR_2026-09-18.md b/experiments/GRAVITY_CONSISTENT_SCHUR_2026-09-18.md new file mode 100644 index 0000000..4c2150b --- /dev/null +++ b/experiments/GRAVITY_CONSISTENT_SCHUR_2026-09-18.md @@ -0,0 +1,119 @@ +# Fixed-work mass-consistent gravity Schur refinement + +This optional preconditioner policy leaves the exact full Newton Jacobian, +normalization, nonlinear merit, and defaults unchanged. It is selected through the +existing typed prescription: + +```cpp +preconditioning::GravityFieldBlock( + backend::MatrixFreeChebyshev{.order = 5, .powerIterations = 20}, + backend::HypreBoomerAMG{backend::FixedCycles{.cycles = 3}}, + preconditioning::GravityConsistentSchurLDU{ + .correctionSweeps = 1, + .relaxation = 0.5 + } +) +``` + +The numbers above specify a benchmark candidate, not a recommended production +setting. `GravityApproximateLDU` remains the default. The new policy composes with +all existing stellar-structure policies and specification-border handling. + +## Algebra and sign + +Write the mixed gravity block in its existing reduced field coordinates as + +\[ + D=\begin{bmatrix}M&B^T\\ B&0\end{bmatrix},\qquad + D\begin{bmatrix}q\\\phi\end{bmatrix} + =\begin{bmatrix}f\\g\end{bmatrix}. +\] + +Let \(R\) denote the prepared mass inverse action used by flux recovery, and let +\(Q\) denote the existing reduced potential inverse action. The latter uses AMG +on the assembled diagonal-mass Schur surrogate. In particular, increasing the +accuracy of \(R\) does not automatically change the operator targeted by \(Q\). + +The positive mass-consistent Schur action and its signed right-hand side are + +\[ + S_R=B R B^T,\qquad r=B R f-g. +\] + +The new policy starts with the existing potential answer \(\phi_0=Q r\), performs +exactly the requested number of corrections + +\[ + \phi_{j+1}=\phi_j+\omega Q(r-S_R\phi_j), +\] + +and uses the existing flux recovery \(q=R f-R B^T\phi_k\). For stationary linear +\(R\), this gives + +\[ + Bq-g=r-S_R\phi_k. +\] + +Thus the corrected residual targets the Poisson row of the same approximate +mixed factorization used for recovery. The gradient-row residual still depends +on the quality of \(R\); this is not an exact inverse of \(D\). Even a converged +Schur refinement would invert the mass-approximated block only, unless \(R=M^{-1}\). + +For fixed linear child inverses the potential error propagates with +\(I-\omega Q S_R\). A positive finite \(\omega\) does not prove contraction. +Neither the toy test nor the implementation establishes spectral bounds for the +stellar models. Benchmarking must measure full-J true defects and cost. + +## Work, lifetime, and sparse actions + +Each additional sweep costs one application each of \(R\), \(Q\), \(B\), and +\(B^T\). No Krylov iteration, convergence test, nonlinear preparation, or full-J +application occurs inside the refinement. Two potential-sized work vectors are +allocated once for positive sweep counts. The existing statistics count every +child application, and `schurRefinementSweeps` records completed corrections. + +With positive sweeps, the prepared policy reuses the cached assembled divergence +already required by the Schur surrogate. The same sparse matrix and reduced DOF +maps supply both its forward and transpose actions throughout the factorization. +The cache precedes the borrowed operator in owner lifetime and is rebound after +discretization invalidation. Geometry-only refresh retains the computational-mesh +divergence while rebuilding the geometry-dependent mass inverse and Schur +surrogate. + +With zero sweeps, the policy keeps the old matrix-free divergence binding and the +old LDU arithmetic, including operation order. This separate path permits exact +baseline comparison. Positive-sweep comparisons include both the refinement and +the equivalent sparse implementation of the coupling; sparse/matrix-free +equivalence is tested independently. + +The policy preserves the existing compile-time application contract: it is +stationary linear exactly when both child backends declare that contract. A +tolerance-driven child remains flexible. Fixed work must not be mistaken for +linearity of an arbitrary child solver. Keeping the default fixed-cycle children +allows the specification border to retain its cached structural responses. + +Sweep counts must be nonnegative; relaxation must be finite and positive. The +distributed prepared block validates these conditions collectively and requires +identical choices on every rank before constructing collective child solvers. +No upper bound on relaxation is claimed to ensure stability. + +## Validation and comparison scope + +The added unit tests check rational answers for a noncommuting dense mass/Schur +example, preservation of an exact Schur fixed point, zero-sweep bitwise baseline +equivalence, linearity, repeated-application independence, fixed work statistics, +and invalid options. Contraction is asserted only for the selected dense fixture. + +The added FEM integration test compares sparse and matrix-free reduced divergence +and transpose actions on three dense vectors, compares the full refined action +against the matrix-free coupling implementation, and checks zero-sweep exact +equivalence using identical child inverses. It repeats these checks after geometry +and discretization revisions and tests stale rejection, cache rebuilding, and +collective option mismatch rejection when run on multiple ranks. + +Tests are provided for the coordinated build/test campaign; no benchmark result +or production improvement is claimed by this implementation note. Compare the +new policy against the default using the same mass backend, AMG cycles, material +inverse, full-J operator, saved right-hand sides, and scalar-border sources. +Evaluate independent and two-way stellar-structure policies separately so a +coupling change is not mistaken for a gravity-inverse improvement. diff --git a/libmeanfield/interface/preconditioning/gravity_field.cppm b/libmeanfield/interface/preconditioning/gravity_field.cppm index 04d948c..c7e687f 100644 --- a/libmeanfield/interface/preconditioning/gravity_field.cppm +++ b/libmeanfield/interface/preconditioning/gravity_field.cppm @@ -1,6 +1,7 @@ module; #include +#include #include #include #include @@ -9,6 +10,7 @@ module; #include #include +#include export module mean_field:preconditioning.gravity_field; @@ -23,12 +25,21 @@ export namespace mean_field::preconditioning { struct GravityUpperTriangular final { }; struct GravityApproximateLDU final { }; + // A fixed number of corrections preserves a stationary linear application + // when both child inverses are stationary. Relaxation is a numerical choice; + // positivity alone does not establish Richardson contraction. + struct GravityConsistentSchurLDU final { + int correctionSweeps{1}; + double relaxation{1.0}; + }; + template struct IsGravityFactorizationPolicy : std::false_type { }; template <> struct IsGravityFactorizationPolicy : std::true_type { }; template <> struct IsGravityFactorizationPolicy : std::true_type { }; template <> struct IsGravityFactorizationPolicy : std::true_type { }; template <> struct IsGravityFactorizationPolicy : std::true_type { }; + template <> struct IsGravityFactorizationPolicy : std::true_type { }; template concept GravityFactorizationPolicy = IsGravityFactorizationPolicy>::value; @@ -179,6 +190,7 @@ export namespace mean_field::preconditioning { std::uint64_t potentialSchurApplications{0}; std::uint64_t divergenceApplications{0}; std::uint64_t transposeDivergenceApplications{0}; + std::uint64_t schurRefinementSweeps{0}; }; template class GravityFactorizationOperator final : public mfem::Solver { @@ -197,7 +209,17 @@ export namespace mean_field::preconditioning { m_offsets(3), m_potentialWorkspace(potentialSchurInverse.Height()), m_gradientWorkspace(massInverse.Height()), - m_massCorrection(massInverse.Height()) { + m_massCorrection(massInverse.Height()), + m_schurResidual(RefinementWorkspaceSize(m_policy, potentialSchurInverse.Height())), + m_potentialCorrection(RefinementWorkspaceSize(m_policy, potentialSchurInverse.Height())) { + if constexpr (std::same_as) { + if (m_policy.correctionSweeps < 0 || !std::isfinite(m_policy.relaxation) || + m_policy.relaxation <= 0.0) { + throw std::invalid_argument( + "Gravity Schur refinement requires nonnegative sweeps and finite positive relaxation." + ); + } + } if (massInverse.Height() <= 0 || massInverse.Height() != massInverse.Width()) { throw std::invalid_argument("The gravity factorization requires a square gradient-mass inverse."); } @@ -269,11 +291,36 @@ export namespace mean_field::preconditioning { ++m_statistics.transposeDivergenceApplications; ++m_statistics.massInverseApplications; } else { - static_assert(std::same_as); + static_assert( + std::same_as || std::same_as + ); m_massInverse->Mult(gradientRightHandSide, gradientAction); m_divergence->Mult(gradientAction, m_potentialWorkspace); m_potentialWorkspace -= potentialRightHandSide; m_potentialSchurInverse->Mult(m_potentialWorkspace, potentialAction); + + if constexpr (std::same_as) { + // For [M B^T; B 0], flux recovery uses R(f-B^T phi), so its + // consistent positive Schur operator is B R B^T and its RHS + // is B R f-g. The assembled diagonal-mass Schur inverse Q + // remains a cheap correction, rather than the target operator. + // Zero sweeps deliberately retain the old LDU arithmetic. + for (int sweep = 0; sweep < m_policy.correctionSweeps; ++sweep) { + m_divergence->MultTranspose(potentialAction, m_gradientWorkspace); + m_massInverse->Mult(m_gradientWorkspace, m_massCorrection); + m_divergence->Mult(m_massCorrection, m_schurResidual); + m_schurResidual *= -1.0; + m_schurResidual += m_potentialWorkspace; + m_potentialSchurInverse->Mult(m_schurResidual, m_potentialCorrection); + potentialAction.Add(m_policy.relaxation, m_potentialCorrection); + ++m_statistics.massInverseApplications; + ++m_statistics.potentialSchurApplications; + ++m_statistics.divergenceApplications; + ++m_statistics.transposeDivergenceApplications; + ++m_statistics.schurRefinementSweeps; + } + } + m_divergence->MultTranspose(potentialAction, m_gradientWorkspace); m_massInverse->Mult(m_gradientWorkspace, m_massCorrection); gradientAction -= m_massCorrection; @@ -295,6 +342,14 @@ export namespace mean_field::preconditioning { } private: + [[nodiscard]] static int RefinementWorkspaceSize(const Policy &policy, const int potentialSize) noexcept { + if constexpr (std::same_as) { + return policy.correctionSweeps > 0 ? potentialSize : 0; + } else { + return 0; + } + } + Policy m_policy; const mfem::Solver *m_massInverse; const mfem::Solver *m_potentialSchurInverse; @@ -303,6 +358,8 @@ export namespace mean_field::preconditioning { mutable mfem::Vector m_potentialWorkspace; mutable mfem::Vector m_gradientWorkspace; mutable mfem::Vector m_massCorrection; + mutable mfem::Vector m_schurResidual; + mutable mfem::Vector m_potentialCorrection; mutable GravityFactorizationStatistics m_statistics; }; @@ -474,12 +531,12 @@ export namespace mean_field::preconditioning { Block block ) : mfem::Solver(GravitySize(geometryContext)), - m_block(std::move(block)), + m_block(ValidateConfiguration(f, std::move(block))), m_geometryContext(std::addressof(geometryContext)), m_gradientMap(geometryContext.GetMassOperator().GetFluxMap()), m_potentialMap(geometryContext.GetSourceOperator().GetPotentialMap()), m_divergence( - geometryContext.GetDivergenceOperator(), + DivergenceForFactorization(f, geometryContext), m_gradientMap, m_potentialMap ), @@ -560,13 +617,13 @@ export namespace mean_field::preconditioning { return report; } - m_divergence.Rebind(geometryContext.GetDivergenceOperator()); + if (report.discretizationChanged) m_fixedDivergence.Clear(); + m_divergence.Rebind(DivergenceForFactorization(f, geometryContext)); report.rebuiltDivergenceBinding = report.discretizationChanged; RefreshMassInverse(geometryContext); report.rebuiltMassInverse = true; - if (report.discretizationChanged) m_fixedDivergence.Clear(); auto potentialSchur = AssemblePotentialSchur(f, geometryContext); m_potentialSchurInverse.Refresh(*potentialSchur); m_potentialSchurSurrogate = std::move(potentialSchur); @@ -614,6 +671,56 @@ export namespace mean_field::preconditioning { } private: + [[nodiscard]] static Block ValidateConfiguration(const fem::FEM &f, Block block) { + if constexpr (std::same_as) { + if (f.mesh == nullptr) { + throw std::invalid_argument("Gravity Schur refinement requires a distributed mesh."); + } + const auto &policy = block.factorizationPolicy(); + const int localValid = policy.correctionSweeps >= 0 && std::isfinite(policy.relaxation) && + policy.relaxation > 0.0; + int globallyValid = 0; + const MPI_Comm communicator = f.mesh->GetComm(); + if (MPI_Allreduce(&localValid, &globallyValid, 1, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS) { + throw std::runtime_error("Could not synchronize gravity Schur refinement validity."); + } + if (!globallyValid) { + throw std::invalid_argument( + "Gravity Schur refinement requires nonnegative sweeps and finite positive relaxation." + ); + } + + // Different work counts would enter collective child inverses in + // different orders; different damping would change the global + // operator. Check once before constructing either child inverse. + const double localOptions[2]{static_cast(policy.correctionSweeps), policy.relaxation}; + double minimumOptions[2], maximumOptions[2]; + if (MPI_Allreduce(localOptions, minimumOptions, 2, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(localOptions, maximumOptions, 2, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) { + throw std::runtime_error("Could not synchronize gravity Schur refinement options."); + } + if (minimumOptions[0] != maximumOptions[0] || minimumOptions[1] != maximumOptions[1]) { + throw std::invalid_argument("Gravity Schur refinement options must agree on every rank."); + } + } + return block; + } + + [[nodiscard]] const mfem::Operator &DivergenceForFactorization( + const fem::FEM &f, + const operators::context::gravity_field::GravityFieldGeometryContext &geometryContext + ) { + if constexpr (std::same_as) { + if (m_block.factorizationPolicy().correctionSweeps > 0) { + // The Schur surrogate already uses this cached sparse + // computational-mesh divergence. Reuse it for the extra + // applications; zero sweeps retain the original arithmetic. + return m_fixedDivergence.Get(f); + } + } + return geometryContext.GetDivergenceOperator(); + } + [[nodiscard]] static int GravitySize(const operators::context::gravity_field::GravityFieldGeometryContext &geometryContext) { if (!geometryContext.IsPrepared()) { @@ -668,9 +775,9 @@ export namespace mean_field::preconditioning { const operators::context::gravity_field::GravityFieldGeometryContext *m_geometryContext; field::FieldDofMap m_gradientMap; field::FieldDofMap m_potentialMap; + PreparedGravityDivergenceCache m_fixedDivergence; ReducedGravityDivergenceOperator m_divergence; PreparedMassInverse m_massInverse; - PreparedGravityDivergenceCache m_fixedDivergence; std::unique_ptr m_potentialSchurSurrogate; backend::PreparedHypreBoomerAMG m_potentialSchurInverse; ReducedFieldSolverAdapter m_reducedPotentialSchurInverse; diff --git a/tests/preconditioning/gravity_field.cpp b/tests/preconditioning/gravity_field.cpp index abdcdb1..7675164 100644 --- a/tests/preconditioning/gravity_field.cpp +++ b/tests/preconditioning/gravity_field.cpp @@ -1,6 +1,8 @@ #include #include +#include #include +#include #include #include #include @@ -26,6 +28,10 @@ namespace { GravityFieldBlock; using AdaptiveGravityLDU = preconditioning::GravityFieldBlock; + using ConsistentGravityLDU = preconditioning:: + GravityFieldBlock; + using AdaptiveConsistentGravityLDU = preconditioning:: + GravityFieldBlock; using DensityIdentity = preconditioning::IdentityBlock; @@ -141,6 +147,72 @@ namespace { } } + struct DenseSchurRefinementFixture final { + static mfem::DenseMatrix makeMass() { + mfem::DenseMatrix result(3); + result = 0.0; + result(0, 0) = 4.0; + result(1, 1) = 3.0; + result(2, 2) = 2.0; + result(0, 1) = result(1, 0) = 1.0; + result(1, 2) = result(2, 1) = 1.0; + return result; + } + + static mfem::DenseMatrix makeDivergence() { + mfem::DenseMatrix result(2, 3); + result(0, 0) = 1.0; + result(0, 1) = 2.0; + result(0, 2) = -1.0; + result(1, 0) = 2.0; + result(1, 1) = -1.0; + result(1, 2) = 1.0; + return result; + } + + static mfem::DenseMatrix makeSurrogate() { + // B diag(M)^-1 B^T differs from B M^-1 B^T. Its inverse + // does not commute with the latter, so a reversed correction + // order cannot pass the rational-reference checks below. + mfem::DenseMatrix result(2); + result(0, 0) = 25.0 / 12.0; + result(0, 1) = result(1, 0) = -2.0 / 3.0; + result(1, 1) = 11.0 / 6.0; + return result; + } + + mfem::DenseMatrix mass{makeMass()}; + mfem::DenseMatrix divergence{makeDivergence()}; + mfem::DenseMatrix surrogate{makeSurrogate()}; + backend::PreparedDenseDirect massInverse{backend::DenseDirect{}, mass}; + backend::PreparedDenseDirect schurInverse{backend::DenseDirect{}, surrogate}; + + static mfem::Vector rightHandSide() { + mfem::Vector result(5); + result(0) = 3.0; + result(1) = -2.0; + result(2) = 1.0; + result(3) = 4.0; + result(4) = -3.0; + return result; + } + + mfem::Vector residual(const mfem::Vector &action, const mfem::Vector &rhs) const { + mfem::Vector result(5); + mfem::Vector flux(const_cast(action.GetData()), 3); + mfem::Vector potential(const_cast(action.GetData()) + 3, 2); + mfem::Vector gradientResidual(result.GetData(), 3); + mfem::Vector potentialResidual(result.GetData() + 3, 2); + mfem::Vector transposeAction(3); + mass.Mult(flux, gradientResidual); + divergence.MultTranspose(potential, transposeAction); + gradientResidual += transposeAction; + divergence.Mult(flux, potentialResidual); + result -= rhs; + return result; + } + }; + struct PreparedGeometry final { mean_field::fem::FEM finiteElements; gravity_context::GravityFieldGeometryContext context; @@ -175,9 +247,15 @@ TEST_CASE( STATIC_CHECK(preconditioning::PreconditionerComponent); STATIC_CHECK(preconditioning::PreconditionerComponent); STATIC_CHECK(preconditioning::PreconditionerComponent); + STATIC_CHECK(preconditioning::PreconditionerComponent); + STATIC_CHECK(preconditioning::GravityFactorizationPolicy); STATIC_CHECK(preconditioning::backend::ArnoldiAdmissible); STATIC_CHECK(preconditioning::backend::ArnoldiAdmissible); STATIC_CHECK_FALSE(preconditioning::backend::ArnoldiAdmissible); + STATIC_CHECK(preconditioning::backend::ArnoldiAdmissible); + STATIC_CHECK_FALSE( + preconditioning::backend::ArnoldiAdmissible + ); STATIC_CHECK(preconditioning::CompletePreconditionerFor); STATIC_CHECK(preconditioning::CompatiblePreconditionerFor); STATIC_CHECK(preconditioning::StationaryLinearPreconditionerPlan); @@ -216,6 +294,123 @@ TEST_CASE( tags::preconditioning_gravity_unit ) { checkExactDenseRecovery(preconditioning::GravityApproximateLDU{}); + checkExactDenseRecovery(preconditioning::GravityConsistentSchurLDU{.correctionSweeps = 2, .relaxation = 0.5}); +} + +TEST_CASE( + "Mass-Consistent Gravity Schur Refinement Preserves The Baseline And Signed Algebra", + tags::preconditioning_gravity_unit +) { + DenseSchurRefinementFixture fixture; + const mfem::Vector rhs = fixture.rightHandSide(); + preconditioning::GravityFactorizationOperator baseline( + preconditioning::GravityApproximateLDU{}, fixture.massInverse, fixture.schurInverse, fixture.divergence + ); + mfem::Vector baselineAction(5); + baseline.Mult(rhs, baselineAction); + const double baselineResidual = fixture.residual(baselineAction, rhs).Norml2(); + + // These rational answers were obtained from the two dense saddle-point + // blocks, with S_R=[3,-2;-2,59/18] and omega=1/2. They exercise a genuinely + // inexact Schur inverse rather than an already-exact scalar example. + const std::array, 2> expected{{ + {-458734.0 / 4782969.0, 6011089.0 / 4782969.0, -7518059.0 / 4782969.0, + -646046.0 / 531441.0, 888176.0 / 531441.0}, + {-631839955.0 / 6973568802.0, 4376327111.0 / 3486784401.0, -10946970773.0 / 6973568802.0, + -472857412.0 / 387420489.0, 644634595.0 / 387420489.0} + }}; + + double previousResidual = baselineResidual; + for (const int sweeps : {0, 1, 2}) { + CAPTURE(sweeps); + preconditioning::GravityFactorizationOperator refined( + preconditioning::GravityConsistentSchurLDU{.correctionSweeps = sweeps, .relaxation = 0.5}, + fixture.massInverse, fixture.schurInverse, fixture.divergence + ); + mfem::Vector action(5); + action = std::numeric_limits::quiet_NaN(); + const auto *const actionStorage = action.GetData(); + refined.Mult(rhs, action); + CHECK(action.GetData() == actionStorage); + for (int index = 0; index < action.Size(); ++index) { + if (sweeps == 0) { + CHECK(std::bit_cast(action(index)) == + std::bit_cast(baselineAction(index))); + } else { + CHECK(action(index) == Catch::Approx(expected[sweeps - 1][index]).margin(2.0e-13)); + } + } + const auto &statistics = refined.GetStatistics(); + CHECK(statistics.applications == 1); + CHECK(statistics.massInverseApplications == static_cast(2 + sweeps)); + CHECK(statistics.potentialSchurApplications == static_cast(1 + sweeps)); + CHECK(statistics.divergenceApplications == static_cast(1 + sweeps)); + CHECK(statistics.transposeDivergenceApplications == static_cast(1 + sweeps)); + CHECK(statistics.schurRefinementSweeps == static_cast(sweeps)); + + const mfem::Vector residual = fixture.residual(action, rhs); + for (int index = 0; index < 3; ++index) CHECK(std::abs(residual(index)) < 2.0e-13); + if (sweeps > 0) { + CHECK(residual.Norml2() < previousResidual); + previousResidual = residual.Norml2(); + } + } +} + +TEST_CASE( + "Fixed Gravity Schur Refinement Is Linear And Has No Application History", + tags::preconditioning_gravity_unit +) { + DenseSchurRefinementFixture fixture; + preconditioning::GravityFactorizationOperator refined( + preconditioning::GravityConsistentSchurLDU{.correctionSweeps = 2, .relaxation = 0.5}, + fixture.massInverse, fixture.schurInverse, fixture.divergence + ); + const mfem::Vector x = fixture.rightHandSide(); + mfem::Vector y(5), combination(x), px(5), py(5), actual(5), expected(5), repeat(5); + for (int index = 0; index < y.Size(); ++index) y(index) = std::sin(0.7 * (index + 1)); + combination *= -2.3; + combination.Add(0.4, y); + refined.Mult(x, px); + refined.Mult(y, py); + refined.Mult(combination, actual); + expected = px; + expected *= -2.3; + expected.Add(0.4, py); + actual -= expected; + CHECK(actual.Norml2() < 2.0e-13 * expected.Norml2()); + repeat = std::numeric_limits::quiet_NaN(); + refined.Mult(x, repeat); + for (int index = 0; index < repeat.Size(); ++index) CHECK(repeat(index) == px(index)); + + combination = 0.0; + refined.Mult(combination, actual); + CHECK(actual.Norml2() == 0.0); + CHECK(refined.GetStatistics().applications == 5); + CHECK(refined.GetStatistics().schurRefinementSweeps == 10); +} + +TEST_CASE( + "Gravity Schur Refinement Rejects Invalid Fixed-Work Options", + tags::preconditioning_gravity_unit +) { + DenseSchurRefinementFixture fixture; + const std::array invalid{{ + {.correctionSweeps = -1, .relaxation = 1.0}, + {.correctionSweeps = 0, .relaxation = 0.0}, + {.correctionSweeps = 1, .relaxation = -0.5}, + {.correctionSweeps = 1, .relaxation = std::numeric_limits::infinity()}, + {.correctionSweeps = 1, .relaxation = std::numeric_limits::quiet_NaN()} + }}; + for (const auto policy : invalid) { + CAPTURE(policy.correctionSweeps, policy.relaxation); + CHECK_THROWS_AS( + (preconditioning::GravityFactorizationOperator( + policy, fixture.massInverse, fixture.schurInverse, fixture.divergence + )), + std::invalid_argument + ); + } } TEST_CASE( @@ -398,6 +593,109 @@ TEST_CASE( } } +TEST_CASE( + "Prepared Gravity Schur Refinement Reuses Sparse Couplings Across Geometry Revisions", + tags::preconditioning_gravity_integration +) { + const auto arguments = test_utils::setup_args(); + PreparedGeometry geometry(arguments); + const auto makeBlock = [](const int sweeps) { + return preconditioning::GravityFieldBlock( + backend::MatrixFreeChebyshev{.order = 2, .powerIterations = 10}, + FixedAMG{backend::FixedCycles{.cycles = 1}}, + preconditioning::GravityConsistentSchurLDU{.correctionSweeps = sweeps, .relaxation = 0.5} + ); + }; + auto baseline = preconditioning::prepare(geometry.finiteElements, geometry.context, makeBlock(0)); + auto refined = preconditioning::prepare(geometry.finiteElements, geometry.context, makeBlock(2)); + const MPI_Comm communicator = geometry.finiteElements.mesh->GetComm(); + const auto gradientMap = geometry.context.GetMassOperator().GetFluxMap(); + const auto potentialMap = geometry.context.GetSourceOperator().GetPotentialMap(); + const mfem::Vector rhs = gravity_prepared_test_utils::make_deterministic_vector(refined.Width(), 0.29); + mfem::Vector actual(refined.Height()), expected(refined.Height()); + + const auto compareApplications = [&] { + const auto sparse = preconditioning::assembleGravityDivergenceSurrogate(geometry.finiteElements); + preconditioning::ReducedGravityDivergenceOperator sparseReduced(*sparse, gradientMap, potentialMap); + preconditioning::ReducedGravityDivergenceOperator matrixFreeReduced( + geometry.context.GetDivergenceOperator(), gradientMap, potentialMap + ); + for (const double phase : {0.19, 0.53, 0.87}) { + const auto flux = gravity_prepared_test_utils::make_deterministic_vector(sparseReduced.Width(), phase); + const auto potential = + gravity_prepared_test_utils::make_deterministic_vector(sparseReduced.Height(), phase + 0.17); + mfem::Vector sparseAction(sparseReduced.Height()), matrixFreeAction(matrixFreeReduced.Height()); + mfem::Vector sparseTranspose(sparseReduced.Width()), matrixFreeTranspose(matrixFreeReduced.Width()); + sparseReduced.Mult(flux, sparseAction); + matrixFreeReduced.Mult(flux, matrixFreeAction); + sparseReduced.MultTranspose(potential, sparseTranspose); + matrixFreeReduced.MultTranspose(potential, matrixFreeTranspose); + CHECK(gravity_prepared_test_utils::relative_error(sparseAction, matrixFreeAction, communicator) < 2.0e-12); + CHECK(gravity_prepared_test_utils::relative_error(sparseTranspose, matrixFreeTranspose, communicator) < + 2.0e-12); + } + + // Reuse precisely the same child solvers, so zero-sweep equivalence + // checks the path itself rather than independent AMG setup roundoff. + preconditioning::ReducedFieldSolverAdapter baselineSchur(baseline.GetPotentialSchurInverse(), potentialMap); + preconditioning::GravityFactorizationOperator baselineReference( + preconditioning::GravityApproximateLDU{}, baseline.GetMassInverse(), baselineSchur, matrixFreeReduced + ); + baseline.Mult(rhs, actual); + baselineReference.Mult(rhs, expected); + for (int index = 0; index < actual.Size(); ++index) { + CHECK(std::bit_cast(actual(index)) == std::bit_cast(expected(index))); + } + + preconditioning::ReducedFieldSolverAdapter refinedSchur(refined.GetPotentialSchurInverse(), potentialMap); + preconditioning::GravityFactorizationOperator refinedReference( + makeBlock(2).factorizationPolicy(), refined.GetMassInverse(), refinedSchur, matrixFreeReduced + ); + refined.Mult(rhs, actual); + refinedReference.Mult(rhs, expected); + CHECK(gravity_prepared_test_utils::relative_error(actual, expected, communicator) < 2.0e-10); + }; + + compareApplications(); + CHECK(refined.GetDivergenceAssemblyCount() == 1); + CHECK_FALSE(refined.Refresh(geometry.finiteElements, geometry.context).DidAnyWork()); + const mfem::Vector displacementTrue = gravity_prepared_test_utils::make_displacement(geometry.finiteElements, 0.4); + const mfem::Vector displacement = geometry.context.GetDisplacementMap().gather(displacementTrue); + geometry.context.PreparePrimal(displacement, {.value = 1}, {.value = 2}); + CHECK_THROWS_AS(refined.Mult(rhs, actual), std::logic_error); + const auto geometryRefresh = refined.Refresh(geometry.finiteElements, geometry.context); + (void) baseline.Refresh(geometry.finiteElements, geometry.context); + CHECK(geometryRefresh.geometryChanged); + CHECK_FALSE(geometryRefresh.rebuiltDivergenceBinding); + CHECK(refined.GetDivergenceAssemblyCount() == 1); + compareApplications(); + + geometry.context.PreparePrimal(displacement, {.value = 2}, {.value = 2}); + CHECK_THROWS_AS(refined.Mult(rhs, actual), std::logic_error); + const auto discretizationRefresh = refined.Refresh(geometry.finiteElements, geometry.context); + (void) baseline.Refresh(geometry.finiteElements, geometry.context); + CHECK(discretizationRefresh.rebuiltDivergenceBinding); + CHECK(refined.GetDivergenceAssemblyCount() == 2); + compareApplications(); + + int rank = 0, rankCount = 0; + MPI_Comm_rank(communicator, &rank); + MPI_Comm_size(communicator, &rankCount); + if (rankCount > 1) { + const auto checkCollectiveRejection = [&](const preconditioning::GravityConsistentSchurLDU policy) { + const auto block = preconditioning::GravityFieldBlock( + backend::Diagonal{}, FixedAMG{backend::FixedCycles{.cycles = 1}}, policy + ); + CHECK_THROWS_AS( + preconditioning::prepare(geometry.finiteElements, geometry.context, block), std::invalid_argument + ); + }; + checkCollectiveRejection({.correctionSweeps = rank == 0 ? -1 : 1, .relaxation = 0.5}); + checkCollectiveRejection({.correctionSweeps = rank == 0 ? 0 : 1, .relaxation = 0.5}); + checkCollectiveRejection({.correctionSweeps = 1, .relaxation = rank == 0 ? 0.5 : 1.0}); + } +} + TEST_CASE( "Fixed Gravity Divergence Cache Preserves Changing Mass Weights And Invalidates", tags::preconditioning_gravity_integration