MeanField/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md

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 selected gravity preconditioner, using a supported compact action or full-root restriction. 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. The compact gravity applications reduce work inside the inner solve; they do not remove these surrounding full-system actions. The measured improvement to one inner solve therefore does not establish a cheaper reduced outer step. 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. Separately, the integrated serial n=1, J=0 reduced control converged in six full steps and passed all eight virial checks, but took longer than the curvature control. Supported low-level two-rank tests passed; full stellar MPI operation remains unvalidated because of the existing single-rank logical radial surface interpolation restriction.

Resolved FEM difference sampling

The first frozen n=1.5 probe completed all gravity solves and preparations. Its central-difference errors at steps 0.01, 0.001, and 0.0001 were respectively 1.7761e-9, 1.9661e-11, and 7.1208e-11. Only the first exceeded twenty times the recorded arithmetic scale. The observed first-interval order was about 1.956, but the conservative checker correctly left the trend unresolved because the finer endpoint did not meet its resolution criterion.

The diagnostic now uses the explicit, validated step array {0.1, 0.03, 0.01, 0.001} and records it in the summary. The larger steps test whether truncation dominates sufficiently to reveal a central-difference trend; retaining 0.001 also preserves an independent small-step accuracy check. This adds one central pair, hence at most four extra gravity correction solves under the existing two-correction completion budget. The inner tolerance, prepared completion gate, and checker thresholds are unchanged. There is no assumption that larger steps will prepare successfully or automatically pass the trend check.

The updated frozen n=1.5 probe passed all 124 algebraic and invariant checks. Absolute central-difference errors on that grid were approximately 1.7784e-7, 1.5987e-8, 1.7761e-9, and 1.9661e-11, resolving the expected second-order trend without relaxing thresholds. These remain finite-range empirical derivative checks. The recorded arithmetic scale is heuristic, and a gravity residual bound by itself is not a bound on its amplified effect in the retained equations. The result establishes no such unmeasured amplification bound and does not reinterpret floor-limited differences as derivative failures.