MeanField/experiments/GRAVITY_CONSISTENT_SCHUR_2026-09-18.md

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:

preconditioning::GravityFieldBlock(
    backend::MatrixFreeChebyshev{.order = 5, .powerIterations = 20},
    backend::HypreBoomerAMG<backend::FixedCycles>{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 and cache rebuilding. It also contains collective option-mismatch checks for multiple ranks, but those checks have not executed successfully in the actual FEM fixture; see the MPI limit below.

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.

Executed coverage and MPI limit

The retained tests_v1 executable built from commit dc8d4b0501e2775b0d1938005e142411319d78b1 passed the focused serial selection: 20 cases and 441,978 assertions, including the new real-FEM Schur refinement and cached stellar cross-action checks. Its SHA-256 is e11a49ebadc85206fad5a998b5cc6ee3ac074b498d09b3578d9ef5e0c5d0c04b.

The same executable passed 13 cases and 213 assertions per rank on two ranks with selection [gravity_block][unit],[preconditioning][stellar_structure][unit],[coupled_shape_coarse]. This verifies the selected algebra and distributed coarse tests. It does not establish distributed FEM Schur assembly, refresh, or option-mismatch handling.

The broader two-rank selection aborted before reaching the new Schur FEM test. The existing Prepared Gravity Block Diagonal Is Legacy Equivalent And Allocation Stable fixture constructs PreparedMappedHDivMassOperator, whose constructor requires both a stellar element and a compactified vacuum element on each rank. Rank 1 has no vacuum element in the production partition of sandbox.smesh, so construction calls MPI_ABORT at that pre-existing check.

fem::setup_fem currently calls GeneratePartitioning(rankCount, 1) internally; there is no partition override in the existing setup interface. A lightweight, read-only reconstruction of the mesh dual graph and installed METIS 4 partition call, using MFEM's descending adjacency order, reproduced the following split:

Mesh	Rank	Core	Envelope	Vacuum
sandbox.smesh	0	0	220	384
sandbox.smesh	1	448	164	0
Existing one-level refined mesh	0	0	1,673	3,072
Existing one-level refined mesh	1	3,584	1,399	0

The refined mesh therefore has the same obstruction; it was inspected through its partition graph, not subjected to another FEM run. Selecting only the new test cannot remove its shared geometry-construction prerequisite. A deliberate partition distributing every material across both ranks could exercise that fixture, but adding a partition-injection seam is outside this preconditioner experiment. Production mesh setup and partitioning remain unchanged.

The precise current validation claim is serial FEM plus two-rank algebra/coarse coverage. Actual FEM Schur refinement, its distributed sparse-cache lifecycle, and its collective option-rejection branches remain unvalidated on multiple ranks. Logs, executable hashes, and commands are retained in the workspace under diagnostics/preconditioner_improvement_2026-09-18/tests_v1/, specifically provenance.json, serial.log, mpi2.log, mpi2_units.json, and mpi2_units.log.