Field coupling review: higher-index harmonic models

Date: 2026-09-19. This is a read-only source and mathematical review supporting the current native spherical-harmonic convergence experiments. No solver changes or numerical runs were performed for this review. Proposed mass-inverse and closure-path changes below are unimplemented and have no demonstrated convergence or performance benefit. The stellar LDU policies already exist; matched frozen tests, including the newly planned n2/n3 controls, must establish their usefulness here.

The relevant current configuration uses independent Q2 density, H1 Q3 enthalpy and a simultaneous native harmonic surface. Conclusions from the earlier unrestricted nodal-shape campaigns do not establish the bottleneck of this configuration. Density and enthalpy remain separate nonlinear unknowns throughout this discussion.

What is measured, and what is not

The completed n3, J=0 baseline trajectories have different final nonlinear residual distributions after twelve updates:

Outer linear tolerance	Final full residual	Leading shares of squared nonlinear residual
0.03	8.2439198751e-4	Poisson 79.9270%; gravity-gradient relation 13.8543%; closure 6.14387%
0.003	5.5073074058e-4	Closure 71.32895%; Poisson 25.87622%; gravity-gradient relation 1.65855%

These are observations from diagnostics/coupled_convergence_2026-09-19/screen_v1, relative to the workspace containing this worktree. Both cases exhausted the nonlinear update budget, and every recorded linear solve passed its requested gate. They follow different trajectories: the tighter run's smaller final residual is not a same-state comparison of preconditioners or evidence of a residual floor.

A large nonlinear closure residual does not demonstrate a large closure linear defect. It can reflect nonlinear EOS curvature, changes in geometry, or the distribution of the accepted linear error. The current frozen benchmark should separate these with preconditioner_blocks.csv and linear_direction_trial_blocks.csv: compare the closure row of Jp+F with the closure component of F(x+alpha p)-(F+alpha Jp). The latter trial vectors are not saved; the exported block remainder norms remain producer observations rather than independently reassembled FEM evidence.

No claim below identifies one of these mechanisms as the demonstrated cause of the n3 trajectory.

Source facts: exact equations and current inverse

The registry declares density an independent L2 field and enthalpy an independent H1 field; the latter has base order three. See libmeanfield/interface/field/field_registry.cppm:29–42 and :293–317.

On fixed geometry, the closure row is

[ F_\rho=\int_{\Omega_*}v_\rho[\rho-\rho_{\rm EOS}(h)],dV, \qquad \delta F_\rho=M_\rho\delta\rho-C\delta h, ]

where C is the cross-mass operator weighted by d rho_EOS/dh. Preparation caches both the closure mismatch and that derivative on the actual mapped production quadrature. The exact action applies the density mass term minus the weighted enthalpy term. See libmeanfield/impl/operators/prepared_barotropic_closure.cpp:647–675 and :918–931.

The hydrostatic enthalpy block is a mapped H1 mass matrix, accompanied by a cross-mass action from the gravity potential. Surface-condition rows are subsequently replaced. See libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp:710–749 and libmeanfield/interface/preconditioning/material_surface.cppm:1432–1437.

The baseline is not diagonal in all physical fields. Its default material policy is SurfaceThenMaterialTriangular, with diagonal density, enthalpy and surface inverses. Writing a for surface parameters, one application computes

[ \delta a=D_a^{-1}b_a,\qquad \delta h=D_h^{-1}(b_h-A_{ha}\delta a),\qquad \delta\rho=D_\rho^{-1}(b_\rho-A_{\rho a}\delta a-A_{\rho h}\delta h). ]

The h -> rho action is the exact prepared EOS cross derivative. Surface-to-material actions also use the prepared closure and hydrostatic operators. The approximation is in the diagonal mass solves and surface surrogate, together with the missing reverse material-to-surface recovery. References:

Default policy: libmeanfield/interface/preconditioning/material_surface.cppm:679–694.
Exact cross actions: the same file, :783–819.
Triangular application: :997–1009.
Density/enthalpy diagonal assembly: :1425–1438 and libmeanfield/impl/operators/prepared_barotropic_closure.cpp:717–743.

At the next level, IndependentStellarSubsystems applies material/surface and gravity inverses separately. It omits their direct cross actions within the structure inverse. The surrounding generated scalar-border completion is still present, so this does not mean the complete preconditioner lacks every global coupling. See libmeanfield/interface/preconditioning/stellar_structure.cppm:347–351, :861–869, and experiments/plateau_preconditioner_benchmark.hpp:133–137.

The live gravity-schur recipe retains that same independent structure policy and material block. It changes only gravity factorization: fixed correction sweeps target the gravity-only consistent Schur operator B R B^T, using the existing surrogate inverse for corrections. It does not insert EOS/material feedback into that target. See experiments/plateau_newton_preconditioner.hpp:60–65, experiments/plateau_preconditioner_benchmark.hpp:110–137, and libmeanfield/interface/preconditioning/gravity_field.cppm:298–326.

The surface/scalar Schur experiment also depends on the quality of its remaining field inverse B = R_F P0 E_F; exact coarse-row recovery does not make that field inverse exact. See experiments/plateau_surface_border_schur.hpp:147–158.

Near-vacuum derivative: derivation and limits

For s=(n+1)K, the implemented positive-enthalpy relation gives

[ \rho=(h/s)^n,\qquad \chi=\frac{d\rho}{dh}=\frac ns(h/s)^{n-1}. ]

The derivative is finite and constant for n=1, and tends to zero at the surface for n>1. Relative to its central value, a profile with h/h_c=t has chi/chi_c=t^(n-1): at t=0.1, that ratio is 0.1 for n=2 and 0.01 for n=3. These are analytical scaling examples, not measurements of the present endpoints. Source: libmeanfield/interface/eos/polytropic.cppm:145–164.

Neither current material diagonal divides by chi. In particular, M_rho remains the mapped density mass matrix. Vanishing chi therefore does not, by itself, prove that this local material block is singular. It changes the spatial distribution of the thermodynamic coupling; its interaction with approximate mass solves and gravity feedback is a hypothesis requiring frozen evidence.

Relative curvature satisfies rho''/rho'=(n-1)/h for positive h. A direction large relative to local h can consequently have an appreciable nonlinear closure remainder even with an accurate linear solve. At fixed geometry, the exact n=3 closure remainder is

[ F_\rho(x+\alpha p)-F_\rho(x)-\alpha(Jp)\rho =-\int v\rho\frac{3h(\alpha\delta h)^2+(\alpha\delta h)^3}{(4K)^3},dV. ]

Moving geometry adds geometric and mixed terms. The formula is a derivation from the implemented EOS, not a diagnosis of the measured dominant mode.

Conditional field-inverse improvements

Use the existing two-way stellar LDU as a matched control. ApproximateStellarBlockLDU first solves material, propagates its density/surface response into gravity, then propagates the gravity response back into material. Those cross actions include the exact density source, geometry response, potential-to-hydrostatic action, and gravity-gradient-to-mechanical action. See libmeanfield/interface/preconditioning/stellar_structure.cppm:371–384, :453–478, and :484–511. This is an inexpensive implementation choice because the policy already exists. It still approximates the field Schur complement; additional inverse and cross-operator applications must be included in cost comparisons.
Strengthen mapped mass inverses if frozen linear defects implicate them. Independent DG density permits an exact local mapped-mass factorization per element; Q2 hexahedral cells have 27 density coefficients. Retain the exact EOS cross action and strengthen the H1 enthalpy mass inverse with a fixed stationary application. H1 connectivity precludes treating its inverse as independent element solves without assembly/overlap handling. The replaced surface rows also require explicit lifting or a free-node formulation before assuming symmetry and using an SPD solver. Cache factors by the relevant geometry generation and retain typed application-contract checks. These changes are not implemented here; no speedup is assumed.
Consider an EOS-aware field Schur only if simpler controls warrant it. With surface parameters and generated scalars held fixed, let m=(rho,h) and let g denote the mixed gravity fields. The exact algebraic elimination is

[ A_{mm}=\begin{bmatrix}M_\rho&-C\0&H\end{bmatrix},\qquad S_g=A_{gg}-A_{gm}A_{mm}^{-1}A_{mg}. ]

The additional term carries potential -> enthalpy -> density -> gravity feedback missing from gravity-only refinement. The signed operators must come from the actual Jacobian, including boundary rows and normalization. This Schur operator is not guaranteed positive definite, so the existing positive Poisson surrogate cannot simply be declared its exact inverse. This would eliminate correction variables only inside a preconditioner; density remains an independent nonlinear unknown.

The current material prepared implementation specializes its stored inverse types to diagonal backends. A mass upgrade needs an explicit typed prepared-backend extension rather than relabeling the existing diagonal recipe. None of these choices requires additional physics-facing parameters or changed equations.

Conditional affine-target closure trial path

This is a separate nonlinear hypothesis, to consider only if the frozen trial observations identify closure remainder as a substantial obstacle. It is not an implemented or demonstrated remedy.

Keep the full Newton unknowns and direction. Freeze the accepted normalized residual F0 and full action J0 p, and define the closure target

[ L_\rho(\alpha)=R_\rho[F_0+\alpha J_0p]. ]

At an admissible prepared raw trial y=x0+alpha p, solve

[ A_{\rho\rho}(y)c_\rho=L_\rho(\alpha)-F_\rho(y), ]

then alter only density. A_rhorho(y) is the current mapped DG mass block with the actual row/value normalization. It is not the density diagonal approximation. Closure is affine in density at fixed h and geometry, so a fresh preparation can verify the resulting closure residual against L_rho separately from the inner linear solve.

The target must be the affine prediction, not zero closure: it makes the correction zero at alpha=0 and preserves the original Newton tangent. With a locally Lipschitz Jacobian, a bounded mass inverse and sufficiently accurate auxiliary solution, the correction is O(alpha^2) up to arithmetic error. At exact vacuum with 1<n<2, the positive-part power law is C1 with a Hölder derivative; its second derivative is singular. A uniform quadratic remainder must not be claimed there: O(alpha^n) can occur instead. The n=3 polynomial case does not have that particular regularity obstruction.

This operation neither makes density a derived state variable nor guarantees its positivity. It also changes total mass, rotational inertia and mechanical balance. It cannot repair a pressure-preparation rejection when enthalpy and geometry remain unchanged; production independently rejects negative sampled h (libmeanfield/impl/operators/prepared_barotropic_closure.cpp:632–644, with the pressure guard also retained).

Correct composition order: correct density, freshly prepare, then correct gravity against the original frozen target G0+alpha Jg0 p using the changed density source. Gravity-only changes leave h, geometry and closure unchanged, so they preserve the closure target. They can still alter hydrostatic and surface residuals. Full residual acceptance, admissibility, independently checked auxiliary defects and transactional restoration remain mandatory.

The existing strategy lifecycle supports such an experiment but does not provide this ordering automatically. PreparedNewtonStrategy::PrepareTrial invokes gravity first, and third-slot TryImproveTrial follows afterward (libmeanfield/interface/solver/stellar_newton_strategy.cppm:508–521). Inserting closure there would need another gravity correction, or a deliberate typed composite preparation policy ordering closure before gravity. Reuse the fresh affine-target verification pattern in the same file at :192–245; keep the full Newton solve and original direction capture authoritative.

The prepared closure currently exposes complete actions, exact diagonal assembly and its mapped rule visitor, but no reusable exact DG mass inverse (libmeanfield/interface/operators/prepared_barotropic_closure_operator.cppm:73–118). A future policy must supply that inverse or a separately verified equivalent. The decision to build it remains pending the actual frozen linear-defect and trial-remainder data.

13 KiB Raw Permalink Blame History Unescape Escape