MeanField/reports/sources/original_audits/MeanField_audit_2026-09-16.md

MeanField architecture and current-state audit

Date: September 16, 2026. Repository examined: /Users/tboudreaux/Programming/4DSTAR/MeanField.

Scope and confidence

This was a read-only source and artifact audit of the active equilibrium path, its operator, mapping, normalization, preconditioning, solver, seed, test, and experiment infrastructure. No repository files were changed, no solver run was launched, and no fresh build or regression suite was executed. This report is stored in the separate ChatGPT project workspace.

The checked-out commit is 75cc63873992984966ea20011cc04d9984f54be4, but the working tree is substantially newer: 49 tracked files differ, with 3,626 additions and 267 deletions, plus untracked experiments, helper code, documentation, and numerical artifacts. The working tree, not HEAD alone, is the object of this audit. Existing modifications were preserved.

The user's observations—roughly 20 seconds per Newton step at n=1 and zero angular momentum, convergence at low rotation, and plateaus at higher index/rotation—are the current operational baseline. The archived measurements establish some of this independently, but do not establish the current rotation threshold or the cause of the plateau. This audit does not claim exhaustive line-by-line verification of every module or numerical correctness of the entire working tree.

Purpose and numerical problem

MeanField is an experimental numerical test bed intended to inform SERiF/4D-STAR development. The immediate implemented problem is a three-dimensional, self-gravitating, rigidly rotating, barotropic free-boundary equilibrium. The presence of older advection, viscosity, Coriolis, and continuity integrators does not mean those equations participate in the current sandbox solve.

The polytropic EOS implements

[ P=K\rho^{1+1/n},\qquad h=(n+1)K\rho^{1/n},\qquad \rho(h)=\left[h/((n+1)K)\right]^n\quad(h\geq0). ]

The active EOS accepts finite n >= 1 and positive K. Its standalone density-from-enthalpy relation extends to zero density for negative h, but the prepared closure rejects negative enthalpy at its quadrature points. This distinction matters when interpreting rejected Newton trials. The independent density field is checked for finiteness in that preparation path; a successful preparation is not itself a global density-positivity certificate.

With the code's positive rotation potential Psi = |Omega cross (x-center)|^2/2, the hydrostatic equation is the weak form of h + Phi - Psi - C = 0. Gravity uses a mixed gradient–potential formulation, with the stored flux representing grad(Phi), not physical acceleration -grad(Phi). Density closure is weakly enforced in the density test space. Pressure in the mechanical operator is evaluated from enthalpy.

The shape residual is obtained by first assembling the volume mechanical virtual-work residual and then applying the transpose of the surface-to-volume deformation Jacobian. Its Jacobian includes both the derivative of the mechanical residual and the derivative of this pullback. Consequently, solving the shape equations does not independently force every unrestricted volume mechanical test row to zero.

The usual sandbox model additionally imposes:

• Fixed total mass, generating the Bernoulli coordinate C.
• Fixed axial angular momentum, generating a scalar angular velocity along a prescribed axis; both density and changing geometry contribute to the inertia derivative.
• Fixed central density, implemented through a central enthalpy point condition and a hydrostatic border coordinate lambda_c.
• Zero surface pressure, compiled to enthalpy boundary rows replacing the corresponding hydrostatic rows.

Key implementations: EOS, physical core residual/Jacobian, central-density border, rotation coupling.

Architecture

The main flow is: physics specifications -> compiled equilibrium problem and block manifest -> prepared physical operator -> normalized operator and scaled inverse -> FGMRES inside damped Newton -> a report with a borrowed accepted-state view.

Layer	Responsibility and current implementation
Model declarations	model::StellarModel canonicalizes a compile-time specification pack. EOS, surface conditions, integral constraints, and phase conditions declare their own dependencies and contributions.
Compilation	The compiler derives state/residual blocks, Jacobian edges, manifests, normalization metadata, and prepared specification runtimes. PreparedVariadicStellarEquilibriumOperator is the active aggregate.
Discretization	MFEM parallel spaces on paired physical/logical STROID meshes. Field-support maps distinguish reduced solver coordinates from full MFEM true DOFs. Nonconforming support resolution is explicitly rejected.
Geometry	Radial surface parameters generate a vector displacement through a default power-law interior extension (power 2) and fixed-infinity vacuum extension. Kelvin compactification represents the exterior through infinity.
Prepared operators	Cache reference basis information, geometry, thermodynamics, local algebraic blocks, and residuals. Dependency identities/revisions determine invalidation. The global Jacobian is applied without assembling a monolithic matrix.
Normalization	Frozen reference-geometry diagonal Riesz factors and physical scales from M, reference R, and G. Distinct primal and dual coordinate maps.
Preconditioning	Material/surface and gravity block approximations, wrapped by a small dense specification-border correction and coordinate permutation/scaling.
Nonlinear runtime	A context owns the problem, accepted/trial state, normalization, inverse, and linear backend. Trial preparation, acceptance, refresh, rollback, and view invalidation are explicit.
Validation	Catch2 operator/runtime tests plus separate research experiments and a saved-field physical-validation/replay harness.

C++23 modules divide interfaces (libmeanfield/interface) from implementations (libmeanfield/impl). The build uses MFEM, MPI, Hypre, UMFPACK, STROID, XAD, and other installed dependencies. Several old APIs remain compatibility paths; design/handoff documents are historical evidence and must be checked against current source.

The default field families are density and potential L2 order 2, gravity RT index 2, and enthalpy/displacement H1 order 3, plus a scalar surface parameter space. A compile-time option increments registered field orders uniformly. The mesh generator uses order-4 geometry; this is not fourth-order approximation for every field.

The saved coarse n=1 model has 1,216 hex elements (832 stellar, 384 vacuum) and 178,074 root values:

Block	Size
Gravity gradient	98,928
Potential	32,832
Density on stellar support	22,464
Surface shape	866
Enthalpy on stellar support	22,981
Three generated scalar coordinates	3

The root geometry unknown is the surface field, not the much larger full vector displacement. Extension operators and coordinate maps are therefore essential numerical components, not incidental mesh utilities.

Entry points: sandbox, equilibrium problem, field registry, default deformation.

Normalization and preconditioning

For a spatial coefficient block with reference Gram diagonal d_i and characteristic physical scale s, normalized state factors are sqrt(d_i)/s and normalized weak-residual factors are 1/(s sqrt(d_i)). Surface pointwise hydrostatic replacement rows receive special treatment; generated global scalar blocks use declared physical scales.

The wrapper implements x = R x_hat, F_hat = L F, and J_hat = L J R. Its inverse wrapper applies R^-1 P^-1 L^-1. The normalization remains fixed during a solve unless explicitly refreshed. It uses reference computational geometry, including finite reference-domain gravity metrics, rather than an evolving exact physical norm over the infinite exterior. This is a deliberate coordinate policy, not evidence of a missing normalization layer.

The current default inverse has:

• Diagonal density and enthalpy approximations, with enthalpy-to-density coupling.
• An uncalibrated surface mass/Riesz diagonal and a surface-then-material triangular application.
• Independent material/surface and gravity subsystems at the structure level.
• Gravity approximate LDU with a fixed order-5 matrix-free Chebyshev H(div) mass inverse and three fixed BoomerAMG cycles for a potential Schur surrogate.
• A directly solved small dense border Schur matrix. Its construction uses the approximate structure inverse, so it is not the exact Schur complement of the full physical Jacobian.

Alternative triangular structure couplings, approximate LDU policies, calibrated surface scalars, and a surface mass-plus-stiffness surrogate already exist. The September 3 campaign found that apparently better isolated surface approximations did not improve the full system. It selected the defaults primarily using nonrotating n=3 linearizations and explicitly left rotating/refinement generalization unresolved. Its recommendation to add normalization is historical: normalization is now implemented.

References: normalization factors, frozen wrapper, default structure policy, surface surrogate, border Schur assembly.

Newton and initialization

The seed is a spherical Lane–Emden profile projected into independent density and enthalpy spaces, with a separate mixed gravity solve. Initial shape is zero, C is initialized from -GM/R, and angular velocity from target J divided by spherical inertia. A seed at nonzero target J is therefore not already a rotating equilibrium. The seed requires a compatible spherical reference radius.

The sandbox currently requests FGMRES restart 40, relative linear tolerance 0.03, at most 200 linear iterations, and at most 30 Newton iterations. The nonlinear criterion is max(atol, rtol times the initial residual), with sandbox atol=0 and rtol=1e-8. This is materially tighter than the saved validation experiment's atol=1e-8: for its initial residual 1.81e-4, the sandbox-style threshold would be about 1.81e-12.

FGMRES reports a separately recomputed true residual. Newton refuses a correction if that linear solve has not converged, even when it might provide some nonlinear decrease. The forcing tolerance is fixed across Newton steps. Backtracking uses the squared normalized-residual merit, sampled geometry admissibility, and thermodynamic rejection. A geometric preflight estimates safe step length from determinant polynomials at selected sample points. The next linear guess retains (1-alpha) times the previous correction after damping; this is vector reuse, not Krylov recycling.

Accepted-state refresh and rollback are carefully implemented. However, a successful report certifies nonlinear residual convergence, not independently verified physical equilibrium. Capture and serialization methods remain explicit unimplemented stubs. The experiment's saved-field files provide a separate persistence/replay mechanism; its replay currently requires exactly zero rotation.

References: Newton loop, linear solve, Newton options, seed projection.

Evidence for present numerical and performance state

The archived n=1, J=0 benchmark reaches normalized residual 6.60533e-9 in three full Newton steps, with 24, 25, and 28 Krylov iterations. G=M=R=1, K=2/pi, and target central density pi/4. These are actual saved data, not just a README claim.

Diagnostic	Archived final value
Relative mass error	2.36e-11
Energy virial error using P(rho)	3.10e-7
Force-based virial error using P(rho)	3.73e-7
Energy virial error using P(h)	2.62e-7
Gravity energy/force consistency	6.30e-8
Density relative volume L2 error	4.70e-4
Enthalpy relative volume L2 error	3.27e-4
Interior potential relative volume L2 error	1.26e-4
Gravity-gradient relative volume L2 error	5.98e-4
Surface-radius RMS error / R	2.75e-4

Thus the saved model comfortably satisfies the user's 1e-5 virial requirement and also the older 1e-6 screen. Eight existing physical screening checks nevertheless fail. The diagnostic replay finds approximately 1.3% near-surface exterior-potential error; this is outside the interior potential norm above. Weak closure permits pointwise EOS mismatch because rho and h occupy different spaces. The experiment quantitatively demonstrates the n=1 projection contribution to the difference between pressure integrals from h and rho.

The old problematic mesh had a documented corner folding mechanism missed by production samples. The later saved n=1 solve has healthy sampled geometry and substantially better reference conditioning; it should not be described as suffering the same demonstrated fold. Nevertheless, current production inspection and preflight still use finite quadrature sample sets, not an element-wide positivity proof.

A one-level refined mesh and analytic geometry control exist. The refined equilibrium output directory contains input mesh, metadata, and analytic checks, but no accepted equilibrium fields or completed equilibrium convergence results. The corresponding report still describes that solve as in progress. I therefore do not infer demonstrated equilibrium h-convergence from it.

The September 9 optimization report records 1.53x whole-process, 1.70x Newton, and 2.02x linear-solve speedups for that particular batch, with roundoff-level agreement in saved physical metrics. Later working-tree optimizations go substantially further: tensor contractions, compact geometry data, fixed-mesh caches, batched field evaluations, specialized H(div) mass action, cheaper Riesz diagonal assembly, and fixed divergence reuse. The isolation harness documents a later first-Newton-step reference of 24.174 seconds. This is consistent with the user's approximate 10x improvement relative to the old 227–255 second steps, but is not a fresh matched end-to-end benchmark performed here.

Seven saved isolation-check CSVs contain 1,328 check rows, all marked passed. The latest batch includes 182 passing checks and measured incremental speedups of about 2.12x for pressure displacement action, 6.32x for source displacement action, 4.89x for hydrostatic preparation, and 6.61x for Schur preparation. These are operation-specific measurements with different baselines; they must not be multiplied into a claimed global speedup.

Validation provenance matters: the archived full-suite XML contains 592 tests, 588 passing and four failing. The failures are centrifugal virial h-refinement geometry, the n=3 projected-seed closure screen, and two-/four-rank MPI groups. The historical report reproduced those failures in preserved earlier binaries; it did not certify the newer September 14 working tree. The current H(div) constructor still requires both a local stellar and a local vacuum representative on every rank, confirming the partition limitation remains in source.

At audit time, 12 production source files were newer than the release sandbox executable and 30 newer than the release test executable. The isolation harness executable was newer than the inspected production sources. Timestamps are provenance clues, not hashes proving build identity, but they rule out treating the existing sandbox/tests as a demonstrated fresh validation of the latest source.

Evidence: saved physical findings, actual screening results, refinement status, historical full-suite report, latest microbenchmarks.

Main concerns for the next investigation

These distinguish observed implementation limitations from hypotheses about the plateau.

1. Physical acceptance is separate from solver success. The default Newton stopping criterion does not require virial balance. The dedicated validation program is currently specialized to single-rank, nonrotating n=1. General physical diagnostics and a 1e-5 acceptance gate are needed for the requested campaign.

2. Central-density border and branch compatibility deserve early scrutiny. The implemented equations are R_c=h(0)-h(rho_c,target) and R_h <- R_h+lambda_c e_c. This is not a direct constraint on the independent DG density at the center. At the saved equilibrium lambda_c is approximately -4.64e-13, with negligible effect on the unbordered residual. Along a rotating sequence, do not assume independently retained K, M, J, and rho_c targets remain compatible with the unmodified physical equations. Measure lambda_c, its normalized residual contribution, and the unbordered hydrostatic residual. This is a formulation question to investigate, not a proven diagnosis of present stagnation.

3. The default inverse omits structure-level cross coupling and gives shape only a reference mass surrogate. Exact Jacobian couplings are present; the approximation is in the inverse. Its adequacy at higher index/deformation must be measured on the actual troublesome Newton right-hand sides. Existing nonrotating policy rankings need not transfer.

4. Higher n changes surface sensitivity. The implemented d rho/d h tends to zero at h=0 for n>1, whereas n=1 has a nonzero one-sided value. Negative quadrature enthalpy triggers trial rejection. Inspect where rejected trials and weak material/shape responses occur rather than assuming the issue is only a Krylov iteration cap.

5. Globalization is basic and failure modes need separation. Fixed linear forcing, restart length, geometry caps, positivity rejection, and merit decrease can each limit progress. A plateau needs block-resolved residuals, true linear residuals, step lengths, and rejection reasons before selecting a remedy. The much tighter sandbox stopping threshold must be recorded in any comparison.

6. The critical-rotation target is not yet an implemented acceptance coordinate. The active sandbox takes angular momentum; some experiments use a reference Keplerian frequency. A declared Omega/Omega_crit convention, equatorial radius, and equatorial effective-gravity diagnostic are needed before labeling a model 50% or 99% critical. These goals are not established by setting a reference frequency fraction alone.

7. Geometry safety and spatial accuracy remain independent requirements. The present sampled geometry can be healthy while the surface and local gravity field are inaccurate. Near-corner/edge inspection and off-interface gravity sampling should accompany deformed models. A low scalar virial error alone cannot certify the spatial solution.

8. Current-source regression status must be refreshed. Rebuild explicit targets before comparing new solves. Maintain the known failure ledger without weakening assertions. Distributed scalability is presently constrained by local-material ownership assumptions, separate from serial nonlinear convergence.

Validation additions with the highest immediate value

For zero surface pressure, preserve the independent requirement

[ \epsilon_{vir}=\frac{|2T+W+3\Pi|}{|W|}\leq10^{-5},\qquad W=\tfrac12\int\rho\Phi,dV,\quad\Pi=\int P,dV. ]

Evaluate energy- and force-based versions, compare P(h) with P(rho), and vary diagnostic quadrature independently of solve quadrature. If nonzero surface pressure is later enabled, include the boundary-pressure contribution rather than applying the zero-pressure formula unchanged. Do not loosen the existing stricter baseline screens to meet the new minimum requirement.

The useful next tests are:

• A reproducible end-to-end baseline and rotating/index sweep recording state, mesh, source/build identity, solver controls, per-block residuals, unbordered residuals, border contribution, geometry minima, thermodynamic minima, and virial diagnostics.
• Independent radial Lane–Emden comparisons at additional indices, including a noninteger index, plus completed h/p convergence of accepted fields, shape, and global quantities. Separate reference integration from production seed generation.
• Full coupled Jacobian directional checks at accepted nonzero-rotation and higher-index states, including simultaneous angular-velocity, density, enthalpy, and surface perturbations. Use multiple admissible finite-difference step sizes and inspect each block, not just a global norm dominated by gravity.
• Rotating sanity checks: mass and angular momentum closure, T=Omega J/2 for the implemented rigid rotation, Bernoulli variation, surface pressure, symmetry/centering, and equatorial effective gravity. Positive/negative spin partners provide a useful symmetry comparison.
• Independent finite-radius gravitational flux/Poisson balance and exterior monopole checks, with sampling slightly off mesh/block interfaces. This helps distinguish localized trace effects from finite-volume field errors.
• Tensor virial or componentwise force-balance checks to expose anisotropic errors hidden by cancellation in the scalar virial.
• A targeted regression for the historical between-sample fold, plus stronger geometric sampling at large deformation. This is distinct from the existing test that rejects an obvious fold.

Suggested order of subsequent work

First establish a freshly built baseline and preserve one reproducible plateau state with complete block and line-search telemetry. Audit target compatibility and border contamination at that state. Then check its full Jacobian and normalization balance before comparing existing and new inverse policies. Introduce continuation through accepted states and adaptive solve controls as measured experiments, retaining cost-to-solution and physical acceptance as the outcomes.

The intended milestones are understood: robust higher-index equilibria, at least 50% critical rotation, an eventual attempt beyond 99%, and every accepted model virialized to 1e-5 or better. The current infrastructure supports a disciplined investigation, but neither high-rotation convergence nor those models' physical accuracy has yet been demonstrated by the artifacts inspected here.