Add field-feedback controls and document the actual material factorization
This commit is contained in:
@@ -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.
|
||||
|
||||
124
experiments/FIELD_COUPLING_REVIEW_2026-09-19.md
Normal file
124
experiments/FIELD_COUPLING_REVIEW_2026-09-19.md
Normal file
@@ -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<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.
|
||||
@@ -226,5 +226,69 @@
|
||||
"--no-central",
|
||||
"--benchmark-trials"
|
||||
]
|
||||
},
|
||||
{
|
||||
"id": "n2_late_ldu_gravity",
|
||||
"arguments": [
|
||||
"--index",
|
||||
"2",
|
||||
"--angular-momentum",
|
||||
"0",
|
||||
"--surface-prescription",
|
||||
"harmonics",
|
||||
"--surface-degree",
|
||||
"4",
|
||||
"--restart",
|
||||
"160",
|
||||
"--max-linear",
|
||||
"300",
|
||||
"--rtol",
|
||||
"0",
|
||||
"--atol",
|
||||
"1e-12",
|
||||
"--max-newton",
|
||||
"0",
|
||||
"--linear-tol",
|
||||
".003",
|
||||
"--preconditioner-probe",
|
||||
"--preconditioner-variant",
|
||||
"stellar-ldu-gravity-schur",
|
||||
"--restart-snapshot",
|
||||
"../../diagnostics/coupled_convergence_2026-09-19/screen_v1/n2_J0/capture_3/physical_base.txt",
|
||||
"--allow-unconverged-restart",
|
||||
"--no-central"
|
||||
]
|
||||
},
|
||||
{
|
||||
"id": "n3_late_ldu_gravity",
|
||||
"arguments": [
|
||||
"--index",
|
||||
"3",
|
||||
"--angular-momentum",
|
||||
"0",
|
||||
"--surface-prescription",
|
||||
"harmonics",
|
||||
"--surface-degree",
|
||||
"4",
|
||||
"--restart",
|
||||
"160",
|
||||
"--max-linear",
|
||||
"300",
|
||||
"--rtol",
|
||||
"0",
|
||||
"--atol",
|
||||
"1e-12",
|
||||
"--max-newton",
|
||||
"0",
|
||||
"--linear-tol",
|
||||
".003",
|
||||
"--preconditioner-probe",
|
||||
"--preconditioner-variant",
|
||||
"stellar-ldu-gravity-schur",
|
||||
"--restart-snapshot",
|
||||
"../../diagnostics/coupled_convergence_2026-09-19/screen_v1/n3_J0_control/capture_7/physical_base.txt",
|
||||
"--allow-unconverged-restart",
|
||||
"--benchmark-trials"
|
||||
]
|
||||
}
|
||||
]
|
||||
|
||||
Reference in New Issue
Block a user