MeanField/experiments/SPHERICAL_HARMONIC_SURFACE_2026-09-19.md

Native spherical harmonic surface candidate

Branch: research/spherical-harmonic-surface-20260919. Parent: 2171b8b76ad78856323fdbcf3d2c537d9039201f.

The preferred research candidate is degree 4 with all real azimuthal orders: 25 independent surface coefficients. This is a discretization choice, not a change to the physical model. The coefficients and interior fields belong to the same Newton system, including their full cross derivatives. The existing nodal radial prescription remains selectable for comparisons. The immediate physical target is rotation through 50% critical; this integration does not establish that target or choose a converged harmonic cutoff for it.

Coordinates and configuration

For reference boundary nodes with radial directions rhat_i, the scalar radial displacement is

q_i = sum_{l=0..L, |m|<=min(l,M)} a_lm Y_lm(rhat_i).

Omitting M retains every order. The real basis has unit integral of its square over solid angle, includes the Condon–Shortley phase, and uses negative orders for sine and positive orders for cosine. Coefficients have units of length; zero coefficients preserve the reference mesh. The catalog is ordered by increasing degree, then increasing signed order. Degree zero is a genuine one-parameter monopole.

Choice	Independent surface parameters
L=4, all m	25
L=4, maximum absolute m=2	19
L=4, m=0 only	5
L=8, all m	81

No odd-degree or azimuthal modes are silently removed. Restricting to m=0 constrains the shape about the fixed global z axis; it does not impose axisymmetry on the three-dimensional interior fields. The reference center is configurable, but changing the frame is not part of a Newton update.

auto discretization = equilibrium::makeStellarDiscretization(
    std::move(finiteElements),
    normalization::PhysicalRieszDiagonal{dimensions::LengthValue{radius}, G},
    deformation::SphericalHarmonicSurface{
        deformation::SphericalHarmonicSurfaceOptions{.maximumDegree = 4}
    }
);

Set maximumAzimuthalOrder = 0 in the options for the axisymmetric shape comparison. Mode counts and actual sampled independence are checked at compilation onto the mesh. Unresolved modes are rejected, not silently discarded. Prescription/extension compatibility and the coupled operator types remain compile-time constrained; degree, mesh rank and sample dimensions are runtime data.

The sandbox selects this candidate explicitly. Existing two-argument discretization calls retain their nodal behavior. Surface choice lives on the numerical discretization, so the physics-facing EOS, surface-pressure and mass/angular-momentum specification is unchanged.

Mathematical and architectural changes

• SphericalHarmonicSurface satisfies the existing surface-prescription contract and composes with the existing interior and fixed-infinity vacuum extensions. Its cached reference basis has exact forward and transpose actions and a zero derivative of the pullback.
• StellarDiscretizationFor owns typed surface/interior/vacuum prescriptions. The modern equilibrium compiler no longer constructs a hardcoded nodal surface.
• Normalization and the default material/surface inverse assemble a diagonal of the pulled-back reference vector trace mass matrix for modal coordinates. The state/residual Gram factors remain a primal/dual pair. The nodal assembly path is preserved.
• Reference metric assembly checks the FEM identity and an explicit canonical trace-layout capability. The erased bridge supports older volume-only deformation operators without claiming that arbitrary trace layouts can use this metric.
• The calibrated diagonal surface inverse supports the actual ordered parameters. The old nodal H1 surrogate explicitly rejects harmonic coordinates; implementing a projected H1 or a dense small surface Schur block is a later preconditioner experiment.

The mesh geometry and field orders are unchanged: Q4 reference geometry, Q3 displacement. Harmonic nodal samples are interpolated through the existing FE trace. Consequently the resulting surface is not an exact analytic harmonic graph between nodes, and the original reference-mesh interpolation error remains.

This differs deliberately from yesterday's frozen reduced-space diagnostic. That experiment retained a general nodal base q0, area-whitened its radial basis, and used separate normalized QR maps. The native prescription uses fixed raw reference coordinates from the seed onward. Equality of corresponding physical residuals and Jacobians is the validation criterion; identical Krylov trajectories under different diagonal metrics are not expected.

The implementation currently requires one MPI rank, consistent with the existing radial volume extension. Replicated global coefficients are not passed to distributed Newton inner products.

Validation and reproducibility

New tests cover the harmonic catalog and analytic low-order functions, derivative/transpose/pullback identities, genuine sampled alias rejection, equivalence to the same nodal geometry, coupled physical residual/Jacobian equivalence and centered differences with generated rotation, normalization before physical preparation, independent boundary integration of the metric, and default/calibrated inverse recovery. Unsupported metric layouts, a different current FEM and the unsupported H1 combination are rejected explicitly.

The live preset is spherical_harmonic_surface_cases_2026-09-19.json. Native use of the diagnostic driver is explicit:

--surface-prescription harmonics --surface-degree 4

Optionally add --surface-maximum-order 0. Legacy frozen nodal probes are excluded from native mode. Known harmonic snapshot metadata is rejected by legacy loaders; undocumented external snapshots still require the caller to supply compatible coordinates.

Every live case retains q14/q18 virial, EOS, positivity and hydrostatic observations. Additional artifacts contain the exact sampled basis, ordered coefficients, full nodal shape force and its component outside the harmonic span. The excluded force uses the Euclidean nodal-dual norm and is not the solver's normalized residual. Both geometry and force pullback equivalence are checked in the observer.

Evidence directory: ../../diagnostics/spherical_harmonic_surface_2026-09-19 relative to this checkout. Builds, tests, live cases and independent artifact audits run serially. Source reviews can run concurrently.

Implementation validation

• All 10 new harmonic cases passed: 4,273,672 assertions in 174.514 seconds. The large assertion count includes array-element checks; the substantive evidence is the independent analytic, adjoint, coupled finite-difference and nodal-equivalence tests.
• All 64 selected existing regression cases passed: 2,795,921 assertions in 643.882 seconds. These cover nodal/radial/domain deformation, normalization, material/surface preconditioning, modern equilibrium problems and the user API.
• The two-rank invocation exited successfully. It exercised the explicit harmonic rejection paths; tests requiring serial physical preparation were skipped. This does not establish distributed harmonic-solver support.
• The independent artifact auditor passed its valid fixture and rejected six deliberate corruptions (catalog, radial displacement, excluded force, norm, nonfinite data and coefficient count).

Clang initially exhausted its source-location budget when the new interfaces included all of MFEM. Narrowing their includes and moving the non-template metric assembly into an implementation file resolved the test-build failure without changing compiler flags or numerical behavior. The test executable, sandbox and live diagnostic driver all built successfully.

A small modal residual or good virial balance alone is not proof of a physically accurate rotating model.

Live coupled results

All five cases used the same mesh, cold Lane–Emden initialization, no central-density constraint, the baseline line search and a six-update budget. The stopping criterion was ||F|| / ||F_initial|| <= 1e-8 with zero absolute tolerance. Each linear solve used restart 160, at most 200 iterations and a true relative residual gate of 0.03. All 29 accepted updates took a full step with one trial, and all passed that linear gate. Surface coefficients, thermodynamic fields, gravity and global multipliers were solved in the same coupled system throughout.

The runs below establish native integration and a baseline for the next preconditioner investigation. They do not test 50% critical rotation. J=0.01 is the driver's angular-momentum input in its unit-mass/unit-reference-radius/unit-G convention, not a calibrated fraction of critical rotation.

Case	Shape parameters	Newton updates	Final normalized residual	Final / initial residual	Outcome	Newton time (s)	Process wall time (s)
n=1, J=0	25	5	1.52715e-12	9.64952e-9	Converged	146.00	210.55
n=1.5, J=0	25	6	1.29965e-12	5.01751e-9	Converged	221.40	290.70
n=1, J=0.01	25	6	3.86299e-13	4.29274e-10	Converged	197.23	265.89
n=1, J=0.01, m=0	5	6	3.08442e-13	3.42757e-10	Converged	164.90	232.54
n=1.5, J=0, gravity-Schur	25	6	2.19823e-11	8.48662e-8	Update budget exhausted	189.54	257.49

Process time includes setup and diagnostic observations; Newton time is the sum of the recorded iteration timers. These are single-run observations, not repeated timing benchmarks. The gravity-Schur variant used one correction sweep with relaxation 0.5. Its failure to meet the requested tolerance within six updates is not evidence of a plateau: its last update still reduced the residual by about 40 times.

Virial balance and internal correctness

The following are the final q18 measurements. All enthalpy-based energy and force virial errors also pass 1e-5 at q14. The energy form uses the potential binding energy; the force form uses the force binding integral. Both are normalized by the magnitude of the potential binding energy in the existing observer.

Case	Energy virial, P(h)	Force virial, P(h)	Energy virial, P(rho)	Force virial, P(rho)
n=1, J=0	3.87036e-7	2.65938e-7	4.26994e-7	3.05896e-7
n=1.5, J=0	6.55153e-7	6.00565e-7	Undefined	Undefined
n=1, J=0.01	3.86528e-7	2.65956e-7	4.26535e-7	3.05963e-7
n=1, J=0.01, m=0	3.86668e-7	2.65918e-7	4.26626e-7	3.05876e-7
n=1.5, gravity-Schur	6.55158e-7	6.00559e-7	Undefined	Undefined

The n=1 cases have no negative density at the sampled interior quadrature points. For baseline n=1.5, however, the minimum sampled density is -4.07687e-4, the negative-density mass fraction is 2.89318e-5, and the negative-density volume fraction is 0.0186652. The fractional polytropic density EOS is undefined at those samples. Its density-derived pressure integral and both density-based virial measures are therefore undefined, not successful tests. The gravity-Schur endpoint has essentially the same defect.

Case	Scaled EOS density mismatch RMS	Hydrostatic gradient RMS, initial	Hydrostatic gradient RMS, final	Minimum deformation determinant
n=1, J=0	7.79532e-5	1.15196e-3	3.96736e-3	0.998938
n=1.5, J=0	1.25551e-4	1.52810e-3	7.27475e-3	0.998671
n=1, J=0.01	7.80099e-5	1.47640e-3	3.96661e-3	0.997774
n=1, J=0.01, m=0	7.79615e-5	1.47640e-3	3.95930e-3	0.997774

These pointwise diagnostics are not the assembled weak residual. The n=1.5 hydrostatic gradient RMS worsens by 4.76 times while the native residual converges. This reproduces the main physical caveat from yesterday's projected-space experiment; the new implementation has not resolved it. Mass errors are below 5.3e-14 for the four converged cases, and angular-momentum errors are below 3.6e-16.

Quadrature sensitivity is small for the enthalpy energy virial: q18 versus q14 differs by at most 7.7e-14 for these n=1 cases and about 1.83e-9 for n=1.5. The n=1.5 negative-density mass fraction, in contrast, changes by about 30% between these quadratures. Good integral virial stability must not be used as evidence that the positivity defect is resolved or fully measured.

What the restricted surface balance leaves out

The final full nodal mechanical force is almost entirely orthogonal to the selected harmonic span. Its Euclidean norm is approximately 9.22e-6 for n=1, J=0 and 8.99e-6 for n=1.5, J=0. In the latter case, the retained modal force norm is only 3.61e-14. The independent audit reproduces the basis pullback and excluded-force projection from saved vectors.

This is convergence of the selected Galerkin shape problem, not convergence of the former unrestricted nodal shape problem. Orthogonality of the remaining force is expected at a modal root; the absolute excluded force still needs a degree/mesh convergence study and an appropriate physical or dual-norm interpretation. Its dimensional Euclidean norm cannot be compared directly with either the normalized nonlinear residual or the dimensionless virial tolerance.

All signed m modes remain available and coupled in the all-m representation. For example, the final real (l,m)=(4,4) coefficient is 8.97263e-5 for n=1, J=0 and 9.23867e-4 for n=1.5, J=0, in radial-length units. The latter is comparable to that case's (4,0) coefficient, 1.09313e-3. Nonaxisymmetric corrections in these nominally spherical models are a reason to investigate mesh/discretization imprint under refinement, not evidence of a physical nonaxisymmetric branch. L=4 is a preferred starting resolution, not a validated universal cutoff.

Implications for the next preconditioner study

Case	Krylov iterations by Newton update	Total Krylov iterations	Linear-solve time (s)
n=1, J=0	19, 23, 28, 74, 68	212	89.79
n=1.5, J=0	25, 26, 26, 84, 87, 87	335	150.06
n=1, J=0.01	18, 20, 22, 48, 92, 96	296	129.96
n=1, J=0.01, m=0	18, 20, 22, 27, 65, 67	219	97.29
n=1.5, gravity-Schur	18, 20, 18, 22, 70, 58	206	116.25

The m=0 comparison reduces total Krylov iterations by 26.0% and Newton time by 16.4% for this low-rotation case, with nearly identical virial measures. It is a useful selectable control. Restricting m also changes the admissible shape, so this is not a pure preconditioner speedup. All m remain the preferred 3D candidate.

Gravity-Schur reduces Krylov work in the six-update comparison but changes the inexact Newton trajectory and does not reach the common stopping tolerance within that budget. Its lower time is not a demonstrated time-to-convergence improvement. Late iteration counts increase substantially even in the 25-parameter system; the recorded initial-linear-residual/RHS ratios remain essentially one, so this observation is not explained by an increasingly poor recycled initial guess.

The next focused experiment should measure the coupled surface Schur operator at identical saved states. With surface coefficients a and all remaining unknowns u, the target is

S = J_aa - J_au J_uu^{-1} J_ua.

The field response must include gravity, thermodynamics and the relevant mass/angular-momentum borders. The existing diagonal boundary metric supplies coordinate scaling but does not model this response. Twenty-five coefficients make a small dense reference Schur matrix feasible as a diagnostic, although its field solves may be expensive. First verify the inner field-block solve defects and determine whether that block has a usable inverse. Measure the Schur singular values, conditioning, spectrum, scaling and setup cost before choosing a reusable approximation; eigenvalues alone do not characterize weak inverse response for a nonsymmetric operator. Compare default diagonal, calibrated diagonal and any new Schur approximation at the same state, then confirm complete nonlinear time to a common tolerance. A block elimination inside the preconditioner still serves the simultaneous coupled Newton solve; it does not introduce a separate surface-relaxation algorithm.

No new gauge constraint, field order, positivity modification or rotation-criticality model is introduced by this integration. Possible weak directions must be measured before deciding that any mode is a gauge or should be removed. Increasing n or J remains subject to the recorded local validation defects as well as the 1e-5 virial milestone.

Evidence binding

Implementation and live producer source: 7b0eb8d92d0f550f1fe3fb19cd8eecc52dc80c56. The tracked source was clean during the campaign. The report-only follow-up commit does not alter the producer.

• Mesh SHA256: ce11ae99e21e6c3bbfbee402acd2e191c1da0d8261d2227b4203f73f2f337a74.
• Copied live binary SHA256: 855d9b6bc8af326d689840b7eb8bf9c4f85d5f2584ef30c225247803a86ded36.
• implementation_validation.json records build/test evidence and executable hashes.
• campaign/provenance.json records the branch, source commit, toolchain, mesh, binary and exact case arguments. Each case's run record binds output hashes and timing.
• campaign_audit.json records a successful independent audit of all five cases, saved algebra and scalar identities. This is an artifact-consistency result, not a model-validity certificate or a fresh FEM integration.
• results_summary.json contains the compact quantities used in this report; the original CSVs, logs and accepted states remain available in campaign/.

To repeat the artifact audit from this checkout:

/opt/homebrew/anaconda3/envs/general/bin/python experiments/audit_spherical_harmonic_surface.py ../../diagnostics/spherical_harmonic_surface_2026-09-19/campaign --output ../../diagnostics/spherical_harmonic_surface_2026-09-19/campaign_audit.json

The preferred surface candidate is now integrated and tested. Acceptance of a physical higher-n or 50%-critical model remains separate from this implementation milestone.