From 6628e237d3c4d837b4931e8f952b49cd9c5a66bf Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 09:23:19 -0400 Subject: [PATCH] Add prepared Schur gravity reduction and contract tests --- ...ITY_REDUCTION_IMPLEMENTATION_2026-09-17.md | 93 ++++ .../interface/solver/gravity_reduction.cppm | 327 ++++++++++++++ tests/solver/gravity_reduction.cpp | 411 ++++++++++++++++++ 3 files changed, 831 insertions(+) create mode 100644 experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md create mode 100644 libmeanfield/interface/solver/gravity_reduction.cppm create mode 100644 tests/solver/gravity_reduction.cpp diff --git a/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md b/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md new file mode 100644 index 0000000..9138be2 --- /dev/null +++ b/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md @@ -0,0 +1,93 @@ +# Reduced gravity linearization + +This branch supplies the linear algebra for the gravity-reduced nonlinear +formulation. It does not change the physics equations, normalization, default +Newton strategy, or model specification API. Gravity completion and accepted/trial +state ownership belong to the nonlinear context; the reduced operator has no +preparation or state-mutation entry point. + +## Contract + +Write the full normalized root equations as `F(u,z)=0`, `G(u,z)=0`, where `z` +contains the gravity gradient and potential. The typed gravity partition selects +those physical blocks from the compiled form. All other root coordinates and +rows, including optional specification borders, stay in the retained system. +Row and column orderings need not agree. + +At a prepared state, `GravityReducedOperator` applies + +``` +S p = F_u p - F_z (G_z^{-1} G_u p). +``` + +`TryLiftTangent` returns `[p, -G_z^{-1} G_u p]`. On the gravity manifold this +is the derivative of the completed state with respect to retained coordinates. +`TryCondenseRightHandSide` and `TryLiftCorrection` also support a nonzero gravity +right-hand side: + +``` +b_reduced = b_F - F_z G_z^{-1} b_G, +delta_z = G_z^{-1} (b_G - G_u delta_u). +``` + +Consequently a small remaining gravity defect is retained in the Newton +equations rather than silently dropped. These formulas require an invertible +gravity block; they do not regularize a singular full stellar problem or remove +the n=3 homology. + +The production gravity subsystem owns its inner backend and restricted full +inverse. The reduced operator borrows it, allocates workspaces once, and checks +its preparation generation through `IsCurrent()` before and after actions. Each +inner solve starts from zero to avoid importing a previous Krylov direction into +the next action. A successful backend report is checked against the requested +true-residual threshold. The `Try` APIs return numerical failures without +replacing the caller's output; the MFEM `Mult` adapter throws the typed +`GravityReductionFailure` because a void action cannot return failure data. +The outer owner must translate this into a failed numerical attempt. + +`GravityReducedPreconditioner` injects a retained residual into the full root, +applies the supplied full inverse, and extracts retained coordinates. The +retained block of the exact full inverse is exactly `S^{-1}`. For an approximate +inverse this provides an interchangeable preconditioner without assembling a +second approximate physics operator. Its stationary/flexible application +contract is propagated at compile time. It never rebinds the full inverse to a +reduced matrix. + +## Accuracy and cost + +The Schur formula is exact; an iterative inner solve evaluates it approximately. +Tolerance-stopped inner Krylov is not necessarily a linear map of its right-hand +side. FGMRES permits varying preconditioners, not arbitrary errors in its +operator. The nonlinear owner must choose sufficiently tight inner tolerances, +verify the lifted correction against the full Jacobian and full right-hand side, +and refuse an inaccurate correction. This gate is also necessary when the state +has only approximately completed gravity. + +An action costs one gravity solve and two full Jacobian applications outside the +inner solve. Projection preconditioning costs one full inverse application. +This first implementation deliberately reuses existing operators and inverses; +its cost should be measured before replacing them with specialized gravity +kernels. All scratch vectors persist across applications. Objects borrow stable +addresses, are not reentrant, and must not outlive their prepared subsystem or +full inverse. + +## Verification + +The focused test source `tests/solver/gravity_reduction.cpp` covers: + +- Central finite differences of a genuinely curved, analytically completed + gravity manifold; freezing gravity produces a detectably different derivative. +- Tangent lifting and vanishing gravity action. +- Schur condensation/lifting versus a direct monolithic solve, including a + nonzero gravity right-hand side and a coupled optional border. +- The exact inverse projection identity and compile-time propagation of both + stationary and flexible inverse contracts. +- Reduced and monolithic Newton convergence to the same known equilibrium. +- Numerical failure propagation without partially published results, true-defect + gate enforcement, stale preparation rejection, and dimension checks. +- Linearity, repeatability, and input/output alias safety with an exact inner + solve. + +These algebraic tests establish the reduction contract. They do not establish +FEM convergence, a faster implementation, or an additional virial-qualified +stellar equilibrium; those require the integrated runtime experiments. diff --git a/libmeanfield/interface/solver/gravity_reduction.cppm b/libmeanfield/interface/solver/gravity_reduction.cppm new file mode 100644 index 0000000..d478518 --- /dev/null +++ b/libmeanfield/interface/solver/gravity_reduction.cppm @@ -0,0 +1,327 @@ +module; + +#include +#include +#include +#include +#include +#include + +#include +#include + +export module mean_field:solver.gravity_reduction; + +export import :solver.gravity_subsystem; +export import :solver.linear_backend; + +export namespace mean_field::solver { + template + concept PreparedGravityReductionSubsystem = requires( + Candidate &subsystem, + const Candidate &constantSubsystem, + const mfem::Vector &rightHandSide, + mfem::Vector &correction, + const LinearSolveControl &control + ) { + { constantSubsystem.GetPartition() } -> std::same_as; + { constantSubsystem.IsCurrent() } -> std::same_as; + { constantSubsystem.GetCommunicator() } -> std::same_as; + requires std::derived_from< + std::remove_cvref_t, mfem::Operator>; + { subsystem.Solve(rightHandSide, correction, control) } -> std::same_as; + }; + + /* + * MFEM's void Mult boundary cannot return a failed inner solve. Keep that + * failure distinguishable from a contract violation so an outer backend + * can stop the numerical attempt without accepting an inaccurate Schur + * action. The Try methods below retain the ordinary report-based path. + */ + class GravityReductionFailure final : public std::runtime_error { + public: + explicit GravityReductionFailure(LinearSolveReport report) + : std::runtime_error("A reduced-gravity action failed its inner true-residual tolerance."), + m_report(std::move(report)) { + } + + [[nodiscard]] const LinearSolveReport &GetReport() const noexcept { + return m_report; + } + + private: + LinearSolveReport m_report; + }; + + /* + * For J = [F_u F_z; G_u G_z], this applies F_u - F_z G_z^{-1} G_u. + * Neither the full operator nor its preparation may change during an + * application. All coordinates, including specification borders, come + * from the same partition used by the gravity subsystem. + * + * A tolerance-stopped inner iteration approximates this linear operator; + * flexible outer preconditioning does not make an inexact operator exact. + * Its owner must verify a lifted correction with the full Jacobian and + * choose an inner tolerance commensurate with the outer accuracy target. + */ + template + class GravityReducedOperator final : public mfem::Operator { + public: + GravityReducedOperator( + Subsystem &subsystem, + LinearSolveControl innerControl + ) + : mfem::Operator( + subsystem.GetPartition().RetainedResidualSize(), + subsystem.GetPartition().RetainedValueSize() + ), + m_subsystem(std::addressof(subsystem)), + m_innerControl(std::move(innerControl)), + m_fullDirection(subsystem.GetPartition().RootValueSize()), + m_fullGravityDirection(subsystem.GetPartition().RootValueSize()), + m_fullAction(subsystem.GetPartition().RootResidualSize()), + m_fullGravityAction(subsystem.GetPartition().RootResidualSize()), + m_gravityRightHandSide(subsystem.GetPartition().GravityResidualSize()), + m_gravityCorrection(subsystem.GetPartition().GravityValueSize()), + m_retainedResult(subsystem.GetPartition().RetainedResidualSize()), + m_retainedGravityAction(subsystem.GetPartition().RetainedResidualSize()) { + m_innerControl.Validate(); + const auto &partition = GetPartition(); + const auto &operation = subsystem.GetFullOperator(); + if (operation.Width() != partition.RootValueSize() || + operation.Height() != partition.RootResidualSize()) { + throw std::invalid_argument("The reduced-gravity partition does not match its full operator."); + } + } + + GravityReducedOperator(const GravityReducedOperator &) = delete; + GravityReducedOperator &operator=(const GravityReducedOperator &) = delete; + GravityReducedOperator(GravityReducedOperator &&) = delete; + GravityReducedOperator &operator=(GravityReducedOperator &&) = delete; + + [[nodiscard]] const GravityBlockPartition &GetPartition() const noexcept { + return m_subsystem->GetPartition(); + } + + [[nodiscard]] const Subsystem &GetSubsystem() const noexcept { + return *m_subsystem; + } + + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { + return m_subsystem->GetCommunicator(); + } + + [[nodiscard]] bool IsCurrent() const noexcept { + return m_subsystem->IsCurrent(); + } + + [[nodiscard]] const LinearSolveControl &GetInnerSolveControl() const noexcept { + return m_innerControl; + } + + [[nodiscard]] const LinearSolveReport &GetLastInnerSolve() const noexcept { + return m_lastInnerSolve; + } + + [[nodiscard]] LinearSolveReport TryMult( + const mfem::Vector &input, + mfem::Vector &output + ) const { + RequireCurrent(); + const auto &partition = GetPartition(); + partition.InjectRetainedValues(input, m_fullDirection); + FullMult(m_fullDirection, m_fullAction); + partition.GatherGravityResiduals(m_fullAction, m_gravityRightHandSide); + if (!SolveGravity()) { + return m_lastInnerSolve; + } + partition.InjectGravityValues(m_gravityCorrection, m_fullGravityDirection); + FullMult(m_fullGravityDirection, m_fullGravityAction); + partition.GatherRetainedResiduals(m_fullAction, m_retainedResult); + partition.GatherRetainedResiduals(m_fullGravityAction, m_retainedGravityAction); + m_retainedResult -= m_retainedGravityAction; + output = m_retainedResult; + return m_lastInnerSolve; + } + + void Mult( + const mfem::Vector &input, + mfem::Vector &output + ) const override { + const auto report = TryMult(input, output); + if (!report.Converged()) { + throw GravityReductionFailure(report); + } + } + + /* + * Condensing a nonzero gravity residual preserves equivalence to the + * monolithic Newton equations. Dropping it merely because gravity was + * approximately completed would turn inner error into an outer bias. + */ + [[nodiscard]] LinearSolveReport TryCondenseRightHandSide( + const mfem::Vector &fullRightHandSide, + mfem::Vector &reducedRightHandSide + ) const { + RequireCurrent(); + const auto &partition = GetPartition(); + partition.GatherGravityResiduals(fullRightHandSide, m_gravityRightHandSide); + if (!SolveGravity()) { + return m_lastInnerSolve; + } + partition.InjectGravityValues(m_gravityCorrection, m_fullGravityDirection); + FullMult(m_fullGravityDirection, m_fullGravityAction); + partition.GatherRetainedResiduals(fullRightHandSide, m_retainedResult); + partition.GatherRetainedResiduals(m_fullGravityAction, m_retainedGravityAction); + m_retainedResult -= m_retainedGravityAction; + reducedRightHandSide = m_retainedResult; + return m_lastInnerSolve; + } + + [[nodiscard]] LinearSolveReport TryLiftCorrection( + const mfem::Vector &retainedCorrection, + const mfem::Vector &fullRightHandSide, + mfem::Vector &fullCorrection + ) const { + RequireCurrent(); + const auto &partition = GetPartition(); + partition.InjectRetainedValues(retainedCorrection, m_fullDirection); + FullMult(m_fullDirection, m_fullAction); + m_fullAction *= -1.0; + if (fullRightHandSide.Size() != m_fullAction.Size()) { + throw std::invalid_argument("A lifted gravity correction requires a full residual right-hand side."); + } + m_fullAction += fullRightHandSide; + partition.GatherGravityResiduals(m_fullAction, m_gravityRightHandSide); + return FinishLift(fullCorrection); + } + + [[nodiscard]] LinearSolveReport TryLiftTangent( + const mfem::Vector &retainedDirection, + mfem::Vector &fullDirection + ) const { + RequireCurrent(); + const auto &partition = GetPartition(); + partition.InjectRetainedValues(retainedDirection, m_fullDirection); + FullMult(m_fullDirection, m_fullAction); + partition.GatherGravityResiduals(m_fullAction, m_gravityRightHandSide); + m_gravityRightHandSide *= -1.0; + return FinishLift(fullDirection); + } + + private: + void RequireCurrent() const { + if (!IsCurrent()) { + throw std::logic_error("A reduced-gravity action requires the current prepared gravity subsystem."); + } + } + + void FullMult( + const mfem::Vector &input, + mfem::Vector &output + ) const { + m_subsystem->GetFullOperator().Mult(input, output); + RequireCurrent(); + } + + [[nodiscard]] bool SolveGravity() const { + // A previous Schur action is not a suitable initial guess for an + // unrelated Krylov direction and can make the action history dependent. + m_gravityCorrection = 0.0; + m_lastInnerSolve = m_subsystem->Solve(m_gravityRightHandSide, m_gravityCorrection, m_innerControl); + RequireCurrent(); + if (m_lastInnerSolve.Converged()) { + if (!std::isfinite(m_lastInnerSolve.rightHandSideNorm) || + !std::isfinite(m_lastInnerSolve.trueResidualNorm)) { + m_lastInnerSolve.status = LinearSolveStatus::non_finite; + } else if (m_lastInnerSolve.rightHandSideNorm < 0.0 || m_lastInnerSolve.trueResidualNorm < 0.0 || + m_lastInnerSolve.trueResidualNorm > + m_innerControl.ConvergenceThreshold(m_lastInnerSolve.rightHandSideNorm)) { + m_lastInnerSolve.status = LinearSolveStatus::backend_failure; + } + } + return m_lastInnerSolve.Converged(); + } + + [[nodiscard]] LinearSolveReport FinishLift(mfem::Vector &fullCorrection) const { + if (!SolveGravity()) { + return m_lastInnerSolve; + } + GetPartition().AddGravityValues(m_gravityCorrection, m_fullDirection); + fullCorrection = m_fullDirection; + return m_lastInnerSolve; + } + + Subsystem *m_subsystem; + LinearSolveControl m_innerControl; + mutable mfem::Vector m_fullDirection; + mutable mfem::Vector m_fullGravityDirection; + mutable mfem::Vector m_fullAction; + mutable mfem::Vector m_fullGravityAction; + mutable mfem::Vector m_gravityRightHandSide; + mutable mfem::Vector m_gravityCorrection; + mutable mfem::Vector m_retainedResult; + mutable mfem::Vector m_retainedGravityAction; + mutable LinearSolveReport m_lastInnerSolve; + }; + + /* + * The retained block of an exact full inverse equals the inverse Schur + * complement. Restricting a supplied approximate inverse preserves that + * relationship without assembling a second approximate physics operator. + */ + template + requires std::derived_from && + LinearPreconditionerApplicationContractAvailable + class GravityReducedPreconditioner final : public mfem::Solver { + public: + static constexpr auto applicationContract = linearPreconditionerApplicationContract; + + GravityReducedPreconditioner( + const GravityReducedOperator &operation, + FullInverse &fullInverse + ) + : mfem::Solver(operation.Width(), operation.Height(), false), + m_operation(std::addressof(operation)), + m_fullInverse(std::addressof(fullInverse)), + m_fullRightHandSide(operation.GetPartition().RootResidualSize()), + m_fullCorrection(operation.GetPartition().RootValueSize()) { + if (fullInverse.Width() != m_fullRightHandSide.Size() || + fullInverse.Height() != m_fullCorrection.Size()) { + throw std::invalid_argument("A reduced-gravity preconditioner requires a compatible full inverse."); + } + } + + void SetOperator(const mfem::Operator &operation) override { + if (std::addressof(operation) != m_operation) { + throw std::invalid_argument("The reduced-gravity preconditioner cannot bind a different operator."); + } + } + + void Mult( + const mfem::Vector &input, + mfem::Vector &output + ) const override { + if (!m_operation->IsCurrent()) { + throw std::logic_error("The reduced-gravity preconditioner requires the current prepared subsystem."); + } + const auto &partition = m_operation->GetPartition(); + partition.InjectRetainedResiduals(input, m_fullRightHandSide); + m_fullInverse->Mult(m_fullRightHandSide, m_fullCorrection); + if (!m_operation->IsCurrent()) { + throw std::logic_error("The full inverse changed the reduced-gravity preparation."); + } + partition.GatherRetainedValues(m_fullCorrection, output); + } + + [[nodiscard]] const FullInverse &GetFullInverse() const noexcept { + return *m_fullInverse; + } + + private: + const GravityReducedOperator *m_operation; + FullInverse *m_fullInverse; + mutable mfem::Vector m_fullRightHandSide; + mutable mfem::Vector m_fullCorrection; + }; +} // namespace mean_field::solver diff --git a/tests/solver/gravity_reduction.cpp b/tests/solver/gravity_reduction.cpp new file mode 100644 index 0000000..b4fe0f2 --- /dev/null +++ b/tests/solver/gravity_reduction.cpp @@ -0,0 +1,411 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +import mean_field; + +namespace { + using namespace mean_field::solver; + + class DenseOperation final : public mfem::Operator { + public: + explicit DenseOperation(mfem::DenseMatrix matrix) + : mfem::Operator(matrix.Height(), matrix.Width()), + m_matrix(std::move(matrix)) { + } + + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + output.SetSize(Height()); + m_matrix.Mult(input, output); + ++applications; + } + + [[nodiscard]] const mfem::DenseMatrix &Matrix() const noexcept { + return m_matrix; + } + + mutable int applications{0}; + + private: + mfem::DenseMatrix m_matrix; + }; + + class DenseGravitySubsystem final { + public: + DenseGravitySubsystem(mfem::DenseMatrix jacobian, GravityBlockPartition partition) + : m_operation(std::move(jacobian)), + m_partition(std::move(partition)), + m_gravityMatrix(m_partition.GravityResidualSize(), m_partition.GravityValueSize()) { + mfem::Vector gravityDirection(m_partition.GravityValueSize()); + mfem::Vector fullDirection(m_partition.RootValueSize()); + mfem::Vector fullAction(m_partition.RootResidualSize()); + mfem::Vector gravityAction(m_partition.GravityResidualSize()); + for (int column = 0; column < gravityDirection.Size(); ++column) { + gravityDirection = 0.0; + gravityDirection(column) = 1.0; + m_partition.InjectGravityValues(gravityDirection, fullDirection); + m_operation.Mult(fullDirection, fullAction); + m_partition.GatherGravityResiduals(fullAction, gravityAction); + for (int row = 0; row < gravityAction.Size(); ++row) { + m_gravityMatrix(row, column) = gravityAction(row); + } + } + } + + [[nodiscard]] const GravityBlockPartition &GetPartition() const noexcept { return m_partition; } + [[nodiscard]] const DenseOperation &GetFullOperator() const noexcept { return m_operation; } + [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { return MPI_COMM_SELF; } + [[nodiscard]] bool IsCurrent() const noexcept { return current; } + + [[nodiscard]] LinearSolveReport Solve( + const mfem::Vector &rightHandSide, + mfem::Vector &correction, + const LinearSolveControl &control + ) { + CHECK(correction.Norml2() == 0.0); + ++solves; + mfem::DenseMatrixInverse inverse(m_gravityMatrix); + inverse.Mult(rightHandSide, correction); + mfem::Vector defect(rightHandSide.Size()); + m_gravityMatrix.Mult(correction, defect); + defect -= rightHandSide; + const double rhsNorm = rightHandSide.Norml2(); + const double trueNorm = forcedTrueNorm >= 0.0 ? forcedTrueNorm : defect.Norml2(); + return { + .status = status, + .control = control, + .iterations = 1, + .rightHandSideNorm = rhsNorm, + .trueResidualNorm = trueNorm, + .relativeTrueResidualNorm = rhsNorm > 0.0 ? trueNorm / rhsNorm : 0.0 + }; + } + + bool current{true}; + int solves{0}; + LinearSolveStatus status{LinearSolveStatus::converged}; + double forcedTrueNorm{-1.0}; + + private: + DenseOperation m_operation; + GravityBlockPartition m_partition; + mfem::DenseMatrix m_gravityMatrix; + }; + + template + class ExactFullInverse final : public mfem::Solver { + public: + static constexpr auto applicationContract = Contract; + + explicit ExactFullInverse(const mfem::DenseMatrix &matrix) + : mfem::Solver(matrix.Width(), matrix.Height(), false), m_matrix(matrix) { + } + + void SetOperator(const mfem::Operator &operation) override { + if (operation.Height() != Width() || operation.Width() != Height()) { + throw std::invalid_argument("The dense test inverse requires matching dimensions."); + } + } + + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + output.SetSize(Height()); + mfem::DenseMatrixInverse inverse(m_matrix); + inverse.Mult(input, output); + } + + private: + mfem::DenseMatrix m_matrix; + }; + + struct MissingGravityContract { }; + static_assert(PreparedGravityReductionSubsystem); + static_assert(!PreparedGravityReductionSubsystem); + using StationaryFullInverse = ExactFullInverse; + using FlexibleFullInverse = ExactFullInverse; + using StationaryReducedInverse = GravityReducedPreconditioner; + using FlexibleReducedInverse = GravityReducedPreconditioner; + static_assert(linearPreconditionerApplicationContract == + mean_field::preconditioning::ApplicationContract::stationary_linear); + static_assert(linearPreconditionerApplicationContract == + mean_field::preconditioning::ApplicationContract::flexible); + + constexpr LinearSolveControl innerControl{.relativeTolerance = 1.0e-12, .absoluteTolerance = 1.0e-14, + .maximumIterations = 20}; + + [[nodiscard]] GravityBlockPartition testPartition(const int size = 4) { + // Different row and column permutations expose accidental reliance on + // the physical root's present ordering, including its appended borders. + return GravityBlockPartition(size, size, {2, 0}, {1, 3}); + } + + [[nodiscard]] mfem::Vector manifoldState(const double first, const double second) { + mfem::Vector state(4); + state(0) = std::sin(second) / (1.0 + first * first); + state(1) = first; + state(2) = first * second; + state(3) = second; + return state; + } + + [[nodiscard]] mfem::Vector fullResidual(const mfem::Vector &state) { + const double first = state(1); + const double second = state(3); + const double gravityFirst = state(2); + const double gravitySecond = state(0); + const auto equilibrium = manifoldState(0.5, 0.25); + mfem::Vector residual(4); + residual(0) = second + 2.0 * gravityFirst - gravitySecond - + (0.25 + 2.0 * equilibrium(2) - equilibrium(0)); + residual(1) = gravityFirst - first * second; + residual(2) = first * first + gravityFirst + 3.0 * gravitySecond - + (0.25 + equilibrium(2) + 3.0 * equilibrium(0)); + residual(3) = (1.0 + first * first) * gravitySecond - std::sin(second); + return residual; + } + + [[nodiscard]] mfem::DenseMatrix fullJacobian(const mfem::Vector &state, const bool border = false) { + const double first = state(1); + const double second = state(3); + mfem::DenseMatrix matrix(border ? 5 : 4); + matrix = 0.0; + matrix(0, 0) = -1.0; + matrix(0, 2) = 2.0; + matrix(0, 3) = 1.0; + matrix(1, 1) = -second; + matrix(1, 2) = 1.0; + matrix(1, 3) = -first; + matrix(2, 0) = 3.0; + matrix(2, 1) = 2.0 * first; + matrix(2, 2) = 1.0; + matrix(3, 0) = 1.0 + first * first; + matrix(3, 1) = 2.0 * first * state(0); + matrix(3, 3) = -std::cos(second); + if (border) { + matrix(4, 4) = 2.0; + matrix(4, 0) = 0.25; + matrix(1, 4) = -0.3; + matrix(2, 4) = 0.1; + matrix(4, 3) = 0.4; + } + return matrix; + } + + [[nodiscard]] double relativeError(const mfem::Vector &actual, const mfem::Vector &expected) { + mfem::Vector difference(actual); + difference -= expected; + return difference.Norml2() / std::max(1.0, expected.Norml2()); + } + + [[nodiscard]] mfem::DenseMatrix explicitReduced(const auto &operation) { + mfem::DenseMatrix matrix(operation.Height(), operation.Width()); + mfem::Vector basis(operation.Width()); + mfem::Vector action(operation.Height()); + for (int column = 0; column < operation.Width(); ++column) { + basis = 0.0; + basis(column) = 1.0; + operation.Mult(basis, action); + for (int row = 0; row < operation.Height(); ++row) { + matrix(row, column) = action(row); + } + } + return matrix; + } +} // namespace + +TEST_CASE("Gravity reduction follows the curved gravity manifold", "[solver][gravity_reduction]") { + const auto state = manifoldState(0.7, -0.4); + DenseGravitySubsystem subsystem(fullJacobian(state), testPartition()); + GravityReducedOperator operation(subsystem, innerControl); + mfem::Vector direction(2); + direction(0) = 0.3; + direction(1) = -0.2; + mfem::Vector action(2); + operation.Mult(direction, action); + + double previousError = 0.0; + for (const double step : {1.0e-2, 1.0e-3, 1.0e-4}) { + const auto plus = fullResidual(manifoldState(state(1) + step * direction(0), + state(3) + step * direction(1))); + const auto minus = fullResidual(manifoldState(state(1) - step * direction(0), + state(3) - step * direction(1))); + mfem::Vector difference(plus); + difference -= minus; + difference /= 2.0 * step; + mfem::Vector reducedDifference(2); + subsystem.GetPartition().GatherRetainedResiduals(difference, reducedDifference); + const double error = relativeError(reducedDifference, action); + if (previousError > 0.0) { + CHECK(error < previousError / 80.0); + } + previousError = error; + } + CHECK(previousError < 5.0e-10); + + mfem::Vector fullDirection(4); + REQUIRE(operation.TryLiftTangent(direction, fullDirection).Converged()); + mfem::Vector fullAction(4); + subsystem.GetFullOperator().Mult(fullDirection, fullAction); + mfem::Vector gravityAction(2), retainedAction(2); + subsystem.GetPartition().GatherGravityResiduals(fullAction, gravityAction); + subsystem.GetPartition().GatherRetainedResiduals(fullAction, retainedAction); + CHECK(gravityAction.Norml2() < 1.0e-14); + CHECK(relativeError(retainedAction, action) < 1.0e-14); + + mfem::Vector frozenDirection(4), frozenAction(4), frozenRetained(2); + subsystem.GetPartition().InjectRetainedValues(direction, frozenDirection); + subsystem.GetFullOperator().Mult(frozenDirection, frozenAction); + subsystem.GetPartition().GatherRetainedResiduals(frozenAction, frozenRetained); + CHECK(relativeError(frozenRetained, action) > 0.1); +} + +TEST_CASE("Gravity condensation and lifting recover the monolithic Newton correction", "[solver][gravity_reduction]") { + for (const bool border : {false, true}) { + CAPTURE(border); + const auto matrix = fullJacobian(manifoldState(0.7, -0.4), border); + const int size = matrix.Height(); + DenseGravitySubsystem subsystem(matrix, testPartition(size)); + GravityReducedOperator operation(subsystem, innerControl); + const auto reducedMatrix = explicitReduced(operation); + mfem::DenseMatrixInverse reducedInverse(reducedMatrix); + mfem::DenseMatrixInverse fullInverse(matrix); + mfem::Vector fullRightHandSide(size), condensed(size - 2), retainedCorrection(size - 2); + mfem::Vector fullCorrection(size), reference(size); + for (int index = 0; index < size; ++index) { + fullRightHandSide(index) = 0.4 - 0.17 * static_cast(index); + } + REQUIRE(operation.TryCondenseRightHandSide(fullRightHandSide, condensed).Converged()); + reducedInverse.Mult(condensed, retainedCorrection); + REQUIRE(operation.TryLiftCorrection(retainedCorrection, fullRightHandSide, fullCorrection).Converged()); + fullInverse.Mult(fullRightHandSide, reference); + CHECK(relativeError(fullCorrection, reference) < 2.0e-13); + mfem::Vector checked(size); + matrix.Mult(fullCorrection, checked); + CHECK(relativeError(checked, fullRightHandSide) < 2.0e-14); + } +} + +TEST_CASE("The projected full inverse preconditions the reduced Schur operator", "[solver][gravity_reduction]") { + const auto matrix = fullJacobian(manifoldState(0.7, -0.4), true); + DenseGravitySubsystem subsystem(matrix, testPartition(5)); + GravityReducedOperator operation(subsystem, innerControl); + StationaryFullInverse fullInverse(matrix); + GravityReducedPreconditioner inverse(operation, fullInverse); + inverse.SetOperator(operation); + mfem::Vector rightHandSide(3), correction(3), actual(3); + rightHandSide(0) = 0.4; + rightHandSide(1) = -0.3; + rightHandSide(2) = 0.8; + inverse.Mult(rightHandSide, correction); + operation.Mult(correction, actual); + CHECK(relativeError(actual, rightHandSide) < 2.0e-14); + CHECK_THROWS_AS(inverse.SetOperator(subsystem.GetFullOperator()), std::invalid_argument); +} + +TEST_CASE("Reduced and monolithic Newton converge to the same curved-system equilibrium", "[solver][gravity_reduction]") { + auto reducedState = manifoldState(0.6, 0.2); + auto monolithicState = reducedState; + for (int iteration = 0; iteration < 8; ++iteration) { + DenseGravitySubsystem subsystem(fullJacobian(reducedState), testPartition()); + GravityReducedOperator operation(subsystem, innerControl); + auto rhs = fullResidual(reducedState); + rhs *= -1.0; + mfem::Vector condensed(2), reducedCorrection(2); + REQUIRE(operation.TryCondenseRightHandSide(rhs, condensed).Converged()); + const auto reducedMatrix = explicitReduced(operation); + mfem::DenseMatrixInverse inverse(reducedMatrix); + inverse.Mult(condensed, reducedCorrection); + reducedState = manifoldState(reducedState(1) + reducedCorrection(0), + reducedState(3) + reducedCorrection(1)); + + const auto matrix = fullJacobian(monolithicState); + mfem::DenseMatrixInverse fullInverse(matrix); + auto monolithicRhs = fullResidual(monolithicState); + monolithicRhs *= -1.0; + mfem::Vector fullCorrection(4); + fullInverse.Mult(monolithicRhs, fullCorrection); + monolithicState += fullCorrection; + } + const auto expected = manifoldState(0.5, 0.25); + CHECK(fullResidual(reducedState).Norml2() < 2.0e-14); + CHECK(fullResidual(monolithicState).Norml2() < 2.0e-14); + CHECK(relativeError(reducedState, expected) < 2.0e-13); + CHECK(relativeError(monolithicState, expected) < 2.0e-13); +} + +TEST_CASE("Reduced gravity actions fail without exposing partial results", "[solver][gravity_reduction]") { + DenseGravitySubsystem subsystem(fullJacobian(manifoldState(0.7, -0.4)), testPartition()); + GravityReducedOperator operation(subsystem, innerControl); + mfem::Vector direction(2), output(2), fullOutput(4), fullRhs(4); + direction = 0.3; + output = 7.0; + fullOutput = 9.0; + fullRhs = 0.4; + subsystem.status = LinearSolveStatus::maximum_iterations; + CHECK_FALSE(operation.TryMult(direction, output).Converged()); + CHECK(output(0) == 7.0); + CHECK(output(1) == 7.0); + CHECK_THROWS_AS(operation.Mult(direction, output), GravityReductionFailure); + CHECK_FALSE(operation.TryCondenseRightHandSide(fullRhs, output).Converged()); + CHECK_FALSE(operation.TryLiftCorrection(direction, fullRhs, fullOutput).Converged()); + CHECK_FALSE(operation.TryLiftTangent(direction, fullOutput).Converged()); + for (int index = 0; index < fullOutput.Size(); ++index) { + CHECK(fullOutput(index) == 9.0); + } + subsystem.status = LinearSolveStatus::converged; + subsystem.forcedTrueNorm = 1.0; + CHECK(operation.TryMult(direction, output).status == LinearSolveStatus::backend_failure); + subsystem.forcedTrueNorm = std::numeric_limits::infinity(); + CHECK(operation.TryMult(direction, output).status == LinearSolveStatus::non_finite); + CHECK(output(0) == 7.0); +} + +TEST_CASE("Reduced gravity actions enforce preparation and dimension contracts", "[solver][gravity_reduction]") { + DenseGravitySubsystem subsystem(fullJacobian(manifoldState(0.7, -0.4)), testPartition()); + GravityReducedOperator operation(subsystem, innerControl); + StationaryFullInverse fullInverse(subsystem.GetFullOperator().Matrix()); + GravityReducedPreconditioner inverse(operation, fullInverse); + mfem::Vector direction(2), output(2), badDirection(3), badRhs(3), fullOutput(4); + direction = 0.2; + badDirection = 0.2; + badRhs = 0.2; + CHECK_THROWS_AS(operation.TryMult(badDirection, output), std::invalid_argument); + CHECK_THROWS_AS(operation.TryCondenseRightHandSide(badRhs, output), std::invalid_argument); + CHECK_THROWS_AS(operation.TryLiftCorrection(direction, badRhs, fullOutput), std::invalid_argument); + subsystem.current = false; + const int applications = subsystem.GetFullOperator().applications; + CHECK_THROWS_AS(operation.TryMult(direction, output), std::logic_error); + CHECK_THROWS_AS(operation.TryLiftTangent(direction, fullOutput), std::logic_error); + CHECK_THROWS_AS(inverse.Mult(direction, output), std::logic_error); + CHECK(subsystem.GetFullOperator().applications == applications); +} + +TEST_CASE("Reduced gravity action is repeatable, linear, and alias safe with an exact inner solve", "[solver][gravity_reduction]") { + DenseGravitySubsystem subsystem(fullJacobian(manifoldState(0.7, -0.4)), testPartition()); + GravityReducedOperator operation(subsystem, innerControl); + mfem::Vector first(2), second(2), sum(2), firstAction(2), secondAction(2), combinedAction(2); + first(0) = 0.3; + first(1) = -0.2; + second(0) = -0.1; + second(1) = 0.5; + operation.Mult(first, firstAction); + operation.Mult(second, secondAction); + sum = first; + sum.Add(1.7, second); + operation.Mult(sum, combinedAction); + firstAction.Add(1.7, secondAction); + CHECK(relativeError(combinedAction, firstAction) < 1.0e-14); + operation.Mult(sum, sum); + CHECK(relativeError(sum, combinedAction) < 1.0e-14); + operation.Mult(first, firstAction); + operation.Mult(first, secondAction); + CHECK(relativeError(firstAction, secondAction) == 0.0); +}