Add prepared Schur gravity reduction and contract tests

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.

libmeanfield/interface/solver/gravity_reduction.cppm

@@ -0,0 +1,327 @@
+module;
+
+#include <cmath>
+#include <concepts>
+#include <memory>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+#include <mfem.hpp>
+#include <mpi.h>
+
+export module mean_field:solver.gravity_reduction;
+
+export import :solver.gravity_subsystem;
+export import :solver.linear_backend;
+
+export namespace mean_field::solver {
+    template <typename Candidate>
+    concept PreparedGravityReductionSubsystem = requires(
+        Candidate &subsystem,
+        const Candidate &constantSubsystem,
+        const mfem::Vector &rightHandSide,
+        mfem::Vector &correction,
+        const LinearSolveControl &control
+    ) {
+        { constantSubsystem.GetPartition() } -> std::same_as<const GravityBlockPartition &>;
+        { constantSubsystem.IsCurrent() } -> std::same_as<bool>;
+        { constantSubsystem.GetCommunicator() } -> std::same_as<MPI_Comm>;
+        requires std::derived_from<
+            std::remove_cvref_t<decltype(constantSubsystem.GetFullOperator())>, mfem::Operator>;
+        { subsystem.Solve(rightHandSide, correction, control) } -> std::same_as<LinearSolveReport>;
+    };
+
+    /*
+     * 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 <PreparedGravityReductionSubsystem Subsystem>
+    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 <PreparedGravityReductionSubsystem Subsystem, typename FullInverse>
+        requires std::derived_from<FullInverse, mfem::Solver> &&
+                 LinearPreconditionerApplicationContractAvailable<FullInverse>
+    class GravityReducedPreconditioner final : public mfem::Solver {
+    public:
+        static constexpr auto applicationContract = linearPreconditionerApplicationContract<FullInverse>;
+
+        GravityReducedPreconditioner(
+            const GravityReducedOperator<Subsystem> &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<Subsystem> *m_operation;
+        FullInverse *m_fullInverse;
+        mutable mfem::Vector m_fullRightHandSide;
+        mutable mfem::Vector m_fullCorrection;
+    };
+} // namespace mean_field::solver

tests/solver/gravity_reduction.cpp

@@ -0,0 +1,411 @@
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+#include <mpi.h>
+
+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 <mean_field::preconditioning::ApplicationContract Contract>
+    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<DenseGravitySubsystem>);
+    static_assert(!PreparedGravityReductionSubsystem<MissingGravityContract>);
+    using StationaryFullInverse = ExactFullInverse<mean_field::preconditioning::ApplicationContract::stationary_linear>;
+    using FlexibleFullInverse = ExactFullInverse<mean_field::preconditioning::ApplicationContract::flexible>;
+    using StationaryReducedInverse = GravityReducedPreconditioner<DenseGravitySubsystem, StationaryFullInverse>;
+    using FlexibleReducedInverse = GravityReducedPreconditioner<DenseGravitySubsystem, FlexibleFullInverse>;
+    static_assert(linearPreconditionerApplicationContract<StationaryReducedInverse> ==
+                  mean_field::preconditioning::ApplicationContract::stationary_linear);
+    static_assert(linearPreconditionerApplicationContract<FlexibleReducedInverse> ==
+                  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<double>(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<double>::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);
+}