diff --git a/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md b/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md index 9685af7..2a454bc 100644 --- a/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md +++ b/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md @@ -128,6 +128,15 @@ inverse at its final saved state, once its original twelve-update budget is complete. This is the combined higher-index/rotation obstruction; it is not assumed to be a failed linear solve, since its update-eight linear gate passes. +Source review also identifies the existing `stellar-ldu-gravity-schur` policy +as a distinct control: it introduces material/gravity feedback that the +independent subsystem policy omits. Add frozen n=2 and n=3 comparisons at +the same captures and eta=.003, with n=3 raw trial observations. The existing +material inverse already includes the exact prepared h-to-density closure +derivative, so a large nonlinear closure residual must not be described as +proof of a missing EOS cross derivative. Frozen block defects and finite-step +closure remainders distinguish those possibilities. + The small Schur block contains all 25 surface parameters and every generated scalar. Its field inverse is the unbordered structure inverse; restricting an already scalar-bordered inverse could spuriously constrain the field solve. diff --git a/experiments/FIELD_COUPLING_REVIEW_2026-09-19.md b/experiments/FIELD_COUPLING_REVIEW_2026-09-19.md new file mode 100644 index 0000000..c5cf233 --- /dev/null +++ b/experiments/FIELD_COUPLING_REVIEW_2026-09-19.md @@ -0,0 +1,124 @@ +# 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 + +1. **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. + +2. **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. + +3. **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