Record density-order study, convergence limits, and log-density implications

This commit is contained in:
2026-09-19 13:12:41 -04:00
parent c0a496e9d6
commit ef76e0c786

View File

@@ -0,0 +1,446 @@
# Independent density order: controlled study
Branch: `research/density-order-study-20260919`.
Parent: `7ed4d23`, the completed negative-density investigation.
The separate log-coordinate assessment is on `research/log-density-assessment-20260919` at `90580b2`.
**Finding:** Q3 removes the resolved n=1 density undershoot, and Q4 substantially
improves higher-index density accuracy. Neither provides positivity at higher n
or resolves the n=3 coupled-solver stall in the common eight-step budget. All
n=1 and n=1.5 controls converge and remain well virialized. Density remains an
independent unknown. Log density was assessed mathematically and in small 1D
controls, not implemented in the coupled stellar solver.
## Question and controls
Does increasing only the independent density approximation order reduce the
observed unphysical undershoots, improve internal EOS consistency, and preserve
virial accuracy at acceptable coupled-solver cost?
The density field remains an independent unknown in every equilibrium solve.
The frozen EOS projections are observations of approximation spaces; none is
installed as a replacement field or used to alter a solver state.
The baseline is discontinuous tensor Q2 density. The alternatives are Q3 and
Q4 in the same Gauss-Legendre nodal family. Enthalpy remains continuous Q3;
gravity retains its RT2/Q2 pair; displacement remains Q3; the mesh and geometry
are unchanged. The surface uses the existing degree-four spherical harmonics
with all azimuthal orders and simultaneous interior/surface solution.
The active density block has 832 times `(p+1)^3` coefficients: 22,464 at Q2,
53,248 at Q3 and 104,000 at Q4. Other field blocks retain their dimensions.
The central-density constraint is present only in the n=3 cases.
## Implementation and mathematical scope
`MEAN_FIELD_DENSITY_POLYNOMIAL_ORDER` selects the base density order at build
time; its default remains two. The existing uniform order increment still
applies, and is zero throughout this study. The field registry supplies the
selected order to typed field spaces, quadrature form queries, adapters,
normalization and compiled block layouts. The physics-facing API is unchanged.
Under the study's default policy the closure selects a rule covering both its
density mass form and EOS source form. Their nominal degrees differ once density can have higher order
than enthalpy. Selecting the greater resolved rule preserves policy overrides
and the Q2 baseline. Adjacent requested orders may use the same MFEM Gauss rule,
so a missing nominal bound does not by itself prove actual underintegration.
Q3 is a useful n=1 control because the linear EOS image of Q3 enthalpy lies in
the density space. For n=3, a general Q3 enthalpy cubed is Q9. For n=1.5 it is
nonpolynomial. Raising density order therefore does not promise positivity for
the higher-index models.
There is also a discrete consistency benefit at Q3 and above. On each cell,
Q3 enthalpy belongs to the DG density test space, even through the common
curved VALUE mapping. If `a` represents that enthalpy in the density basis and
`R` is the physical closure residual, the implemented polytropic EOS gives
`sum_q w_q h_h rho_h - (n+1) sum_q w_q P(h_h) = a^T R`.
A zero closure residual therefore enforces this enthalpy-weighted EOS moment
cell by cell under the closure quadrature. Q2 does not generally enforce it.
This algebraic property does not establish positivity or virial balance;
other quadrature rules and finite residuals introduce their own defects, and
the central-density border does not change this identity.
## Frozen and coupled observations
Frozen observations retain the saved enthalpy and geometry from the n=1 and
n=1.5 nonrotating equilibria and the n=3 initial seed. Each target space receives
an independently assembled physical-volume-weighted L2 projection, with one
common q26 quadrature, and a nodal interpolant of `rho_EOS(h_h)`. The latter is
explicitly not the original independently tabulated Lane-Emden seed.
The squared EOS error has a higher degree than the projection equations. A
separate shared q30 evaluation checks it against q26 without changing or
re-solving the projected fields. Both rules are also used on coupled final
states. At fractional n this is a numerical sensitivity check rather than a
polynomial exactness argument.
The coupled runs retain the same production quadrature policy, whose resolved
rule can change when density order changes. Actual closure rules are saved for
each case. The frozen common-rule controls isolate approximation order more
strictly; the coupled comparison measures the complete consistent
discretization, including its policy-selected quadrature adjustments.
Density order does not enter the pressure-force rule or its enthalpy
admissibility guard: their field operands are enthalpy and the displacement
gradient. Pressure-preparation rejection counts therefore measure changed
reached states and coupled directions, rather than a density-order change in
pressure sampling. The current preparation order evaluates displacement and
pressure before the barotropic closure. Closure therefore cannot preempt that
pressure stage in these runs; all rejection sources are nevertheless retained
separately rather than interpreting a stage count as total admissibility.
Coupled runs use the original Lane-Emden initialization, the same normalization
and default preconditioner, an FGMRES restart of 160, at most 200 linear
iterations per step, relative linear tolerance 0.03, and at most eight Newton
steps. Nonlinear relative tolerance is 1e-8 and absolute tolerance is zero.
Cases cover n=1,1.5,3 without rotation and the established n=1,J=0.01 rotating
control. A failure within this budget is a recorded trajectory outcome, not a
proof that the discretization has no root.
The initializer tabulates Lane–Emden density and enthalpy before constructing
either FEM field, then interpolates them independently into their spaces
(`libmeanfield/impl/seed/lane_emden.cpp`, lines 231–236, and
`libmeanfield/impl/seed/stellar_equilibrium_projection.cpp`, lines 164–180).
Gravity uses the represented density; it does not replace enthalpy. Thus the
production density seed is `I_p(rho_LE_tabulated)`, whereas the frozen nodal
and L2 controls are `I_p(rho_EOS(h_h))` and `Pi_p(rho_EOS(h_h))`.
Nonlinear EOS evaluation does not commute with either radial or FEM
interpolation. Increasing density order does not remove the discrepancy
between the tabulated density target and the EOS image of represented Q3
enthalpy. Source structure supports unchanged seeded enthalpy across these
orders; no cross-run bitwise identity of the initial states is asserted.
For n=3, the existing central-density phase condition constrains central
enthalpy to its EOS value at the requested central density; it does not directly
pin the independent DG density's point value. Bordered convergence is assessed
alongside the residual with artificial central forcing removed and the
normalized size of that forcing. That unbordered diagnostic retains the scalar
constraint rows. A small bordered residual or virial alone cannot establish
that the artificial forcing has vanished.
All expensive builds, numerical cases and artifact audits run serially.
Reported solver costs must exclude the post-solve observer and distinguish
Newton steps, Krylov work, preparation and preconditioner refresh.
The normalized residual and its relative stopping threshold belong to each
order's discrete system. Changing density order changes its test space,
diagonal Riesz weights and initial interpolant as well as possibly its resolved
quadrature. A smaller final residual across orders is not a physical-error
ranking. Krylov totals also reflect different Newton trajectories, warm starts
and per-action costs; they do not by themselves measure preconditioner quality.
Physical diagnostics and measured Newton time are therefore reported alongside
the iteration counts and density block sizes.
Timings are single observations on this workstation, with serial workloads;
they are not repeated benchmark means or estimates with error bars.
Density sign regions are integrated by isolating polynomial roots of degree
two through four along reference-cell lines. The observer also checks surface
traces, EOS mismatch, signed and negative mass, moments and unchanged source
fields. It separates resolved negative witnesses from values at reconstruction
uncertainty. Tangential sampling and mapped integration remain numerical
observations, not global sign certificates.
Here `M = integral rho_h dV` is signed mass,
`M_- = -integral min(rho_h,0) dV`, and `V_-` is the volume where density is
negative. Reported fractions are `M_-/abs(M)` and `V_-/V`; the denominator of
the mass fraction is not the integral of absolute density. EOS L2 error means
`sqrt(integral (rho_h-rho_EOS(h_h))^2 dV)` in the model's units, without an
additional relative normalization.
If the stored-state observer encounters negative enthalpy, complete EOS error
integrals are marked undefined while density integrals remain reported. No
density is clipped. Both virial definitions from pressure evaluated at enthalpy
and, where admissible, at density are retained; the required virial tolerance
remains 1e-5.
The summary enthalpy minimum and EOS-availability flag refer to integration
nodes. Surface enthalpy samples and the enthalpy at exported negative-density
witnesses are reported separately. Sliced mass and volume are cross-checked
against the physical observer; first moments and inertia are only compared,
because their additional coordinate factors need not be integrated exactly by
the fixed eight-point normal rule. Tangential refinement alone cannot establish
accuracy in that normal direction.
## Results
### Frozen approximation controls
The frozen campaign completed all three cases and its final artifact audit
passed. The following uses the common-q26 L2 projection; EOS L2 errors use the
separate unsliced q30 measurement. Negative measures use the finest root-sliced
tangential rule. These are approximation controls, not new equilibrium solves.
| Frozen model | Density order | M_- / abs(M) | Negative volume fraction | EOS L2 error |
|---|---:|---:|---:|---:|
| n=1 equilibrium | Q2 | 6.09585e-8 | 3.27346e-4 | 1.25266e-4 |
| n=1 equilibrium | Q3 | 1.05e-31 (unresolved roundoff) | 2.49e-16 | 1.56e-15 |
| n=1 equilibrium | Q4 | 1.25e-31 (unresolved roundoff) | 2.73e-16 | 1.60e-15 |
| n=1.5 equilibrium | Q2 | 2.24299e-5 | 1.49181e-2 | 3.67965e-4 |
| n=1.5 equilibrium | Q3 | 9.33602e-6 | 1.04476e-2 | 1.12325e-4 |
| n=1.5 equilibrium | Q4 | 3.18362e-6 | 6.83440e-3 | 4.80663e-5 |
| n=3 seed | Q2 | 2.92063e-5 | 1.21064e-1 | 1.06409e-2 |
| n=3 seed | Q3 | 2.63486e-6 | 3.73200e-2 | 6.86223e-4 |
| n=3 seed | Q4 | 1.02066e-8 | 1.69628e-2 | 3.84980e-5 |
The n=1 control matches the representability argument: Q3 and Q4 contain the
linear EOS image of Q3 enthalpy, and no resolved negative witnesses remain.
Higher-order density does not guarantee positivity for either higher-index
case. The n=3 Q4 negative mass is about 2,862 times smaller than Q2, yet negative
density still occupies about 1.70% of the sampled stellar volume.
The difference between n=1.5 and n=3 is consistent with surface regularity.
Where enthalpy vanishes linearly, the fractional-index density behaves as
distance-to-surface raised to 1.5 and has a singular second derivative. For
n=3, cubing the represented enthalpy gives a polynomial Q9 target. Increasing
density order can therefore improve the integer-index approximation very
rapidly without resolving the conditioning of the coupled nonlinear equations.
At n=1.5, changing the error rule from q26 to q30 changes the Q4 EOS L2 norm by
about 0.18%; the corresponding fixed eight-point root-sliced norm differs by
about 2.2%. This is why the unsliced q30 values are used above. Tangential
refinement gives stable negative masses but cannot establish convergence of
the fractional EOS error in the normal direction. The integral-index norms
agree across q26/q30 to arithmetic accuracy, except for relative comparisons of
n=1 Q3/Q4 errors that are already at roundoff.
### Coupled results
All nine n=1 and n=1.5 controls converged with full Newton steps and no rejected
trials. Their energy and force virial errors from `P(h)` are all below 6.6e-7.
The following times sum Newton iterations and exclude initialization and the
post-solve density observer.
| Model | Density | Newton updates | Krylov iterations | Newton seconds | M_- / abs(M) | EOS L2 error, q30 |
|---|---:|---:|---:|---:|---:|---:|
| n=1, J=0 | Q2 | 5 | 212 | 147.5 | 6.09585e-8 | 1.25266e-4 |
| n=1, J=0 | Q3 | 6 | 286 | 204.3 | unresolved roundoff | 3.13581e-16 |
| n=1, J=0 | Q4 | 6 | 286 | 215.1 | unresolved roundoff | 3.51021e-16 |
| n=1, J=0.01 | Q2 | 6 | 296 | 196.1 | 6.11256e-8 | 1.25327e-4 |
| n=1, J=0.01 | Q3 | 6 | 305 | 217.8 | unresolved roundoff | 7.26203e-16 |
| n=1, J=0.01 | Q4 | 6 | 304 | 229.1 | unresolved roundoff | 3.90198e-16 |
| n=1.5, J=0 | Q2 | 6 | 335 | 222.0 | 2.24535e-5 | 3.67965e-4 |
| n=1.5, J=0 | Q3 | 6 | 329 | 241.7 | 9.35454e-6 | 1.12195e-4 |
| n=1.5, J=0 | Q4 | 7 | 397 | 303.7 | 3.20115e-6 | 4.80640e-5 |
For n=1, Q3/Q4 remove the resolved undershoot and reduce the EOS mismatch to
roundoff. Q3 and Q4 also agree in their virial diagnostics to about 1e-14. The
nonrotating Q3/Q4 total times are 39%/46% above Q2, but both took an extra
Newton step: their fifth residuals were 1.91e-12 and 1.80e-12 against thresholds
of about 1.55e-12 and 1.54e-12. Q2's fifth residual narrowly passed its own
threshold. Those total-time differences should not be interpreted as pure
per-iteration order overhead. The rotating control takes six updates at every
order; Q3/Q4 cost about 11%/17% more there.
For n=1.5, Q3 reduces negative mass by 2.40 times and EOS L2 error by 3.28 times
for 8.9% more Newton time. Q4 reduces them by 7.01 and 7.66 times, respectively,
for 36.8% more time. Its sixth residual, 2.26e-12, narrowly missed its 1.93e-12
threshold; the seventh update cost 58.8 seconds. Negative volume fractions
remain 1.493%, 1.046% and 0.685% at Q2, Q3 and Q4. These are converged discrete
roots with a resolved density defect, not merely incomplete Newton iterates.
Density/EOS consistency is not the only internal accuracy measure. At n=1.5,
`hydrostatic_gradient_rms` is 0.00727475, 0.00731970 and 0.00729901 at Q2/Q3/Q4.
It is the unscaled volume RMS of `grad(h) + g_h - grad(Psi)` in model
acceleration units, using the independent gravity-gradient field `g_h`.
It is neither a relative error nor the normalized algebraic residual. Thus
density enrichment has not appreciably improved this local balance diagnostic.
None of the three n=3 cases converged within eight updates. Their endpoints
must not be treated as equally accurate equilibrium models or their times as
time to solution.
| Density | F_final / F_initial | Smallest accepted alpha | Krylov iterations | Newton seconds | Energy virial, P(h) | Force virial, P(h) |
|---|---:|---:|---:|---:|---:|---:|
| Q2 | 0.430062 | 1/64 | 221 | 301.6 | 1.77163e-4 | 1.34597e-4 |
| Q3 | 0.985301 | 1/256 | 452 | 478.8 | 8.07504e-6 | 2.28169e-5 |
| Q4 | 0.987949 | 1/512 | 458 | 508.2 | 1.37586e-5 | 2.84372e-5 |
Q3 reaches the energy-virial milestone alone; its force virial and nonlinear
residual still fail their respective criteria. Q4's better approximation
space does not produce a better virial endpoint than Q3. Q2 recovers larger
steps during this window, while Q3/Q4 make only about 1–1.5% residual progress.
This is evidence of a stalled short trajectory, not a proof against eventual
convergence with a longer run or different initialization.
All recorded n=3 linear solves meet their true-residual tolerance. Q2 has four
pressure-force rejections and four insufficient-decrease rejections; Q3 and
Q4 each have five and nine. No other rejection stages occur in these runs.
Both enriched cases first accept only 1/256 of their Newton direction. The
unchanged pressure quadrature makes this a change in reached states and
directions, not a change in pressure sampling. Increasing density order alone
has not removed pressure-admissibility and insufficient-decrease failures.
This campaign records those failures but does not independently localize the
nonlinear mechanism responsible for insufficient decrease.
The n=3 endpoint negative-mass fractions are 2.22117e-6, 9.13593e-7 and
5.21562e-7, with negative volume fractions 3.440%, 2.931% and 2.619%.
EOS L2 errors are 9.78088e-3, 8.99701e-4 and 5.08742e-4. These still-unconverged
fields show much less improvement than the frozen projections; their
different trajectories prevent interpreting those ratios as equilibrium
convergence rates.
The unbordered n=3 residuals remain 1.33327e-3, 7.18640e-4 and 7.26009e-4.
Normalized central-border forcing is 3.17039e-5, 1.27669e-5 and 1.20966e-5.
Thus neither a physical root nor vanishing artificial forcing has been
demonstrated. Signed multipliers and both virial definitions are retained in
the full tables.
The initializer distinction also matters. At n=3, the production Q4 seed has
an observed EOS L2 mismatch of about 5.26e-4 under q18, while the frozen Q4
projection of represented enthalpy has 3.85e-5 under q30. These are different
fields and quadratures, not a controlled improvement ratio. Together with the
initializer's source structure, they show why refining the density space need
not make the independently interpolated seed close to the weak closure
manifold. A closure-consistent initial guess is a separate hypothesis; it
could populate an independent density unknown without deriving it from
enthalpy during the subsequent coupled solve.
The actual closure rules also differ only at Q4 in this study:
| Model index | Q2 rule: achieved order / points per cell | Q3 rule | Q4 rule |
|---|---:|---:|---:|
| n=1 | 17 / 729 | 17 / 729 | 19 / 1,000 |
| n=1.5 | 19 / 1,000 | 19 / 1,000 | 21 / 1,331 |
| n=3 | 23 / 1,728 | 23 / 1,728 | 25 / 2,197 |
All 832 stellar elements use the listed rule in each case. Q2/Q3 therefore
share closure sampling, while Q4's coupled accuracy/cost includes the policy-selected
quadrature enrichment as well as the larger space.
### Implications and recommendation
Q3 is the economical consistency control for n=1: it contains the EOS image
and Q4 adds no resolved accuracy benefit there. Q4 is a useful higher-index
accuracy reference, but neither higher order is a positivity mechanism or a
demonstrated n=3 convergence remedy. The production default remains Q2; these
are isolated, reproducible research variants.
The next density-representation experiment should test a pointwise exponential
with an independently chosen vacuum reference, retaining density as an
independent unknown and solving surface and interior together. Judge it by
EOS consistency, accepted Newton steps, both virial diagnostics and cost—not
positivity alone. A separate conservative, closure-consistent seed control
would distinguish initialization effects from representation effects. Neither
of these further stellar experiments was performed here.
The complete tables and checked PNG/PDF figure are in
`../../diagnostics/density_order_2026-09-19/final_analysis/` relative to the
worktree. `trial_stages.csv` retains every source/disposition count, including
accepted trials; `coupled_outcomes.csv` retains all endpoint diagnostics.
[Comparison figure](../../../diagnostics/density_order_2026-09-19/final_analysis/density_order_comparison.png)
and [coupled-result table](../../../diagnostics/density_order_2026-09-19/final_analysis/coupled_outcomes.csv).
## Relation to the independent log-density hypothesis
The separate `LOG_DENSITY_ASSESSMENT_2026-09-19.md` in the log-density worktree
contains the derivation, primary-source MESA comparison and reproducible
one-dimensional controls. It does not contain a coupled stellar log solve.
Exponentiating polynomial nodal coefficients cannot remove the observed
between-node undershoots: the current negative fields can already have positive
nodal values. A pointwise representation `rho = rho_scale * exp(u_h)` changes
the trial family and gives positive interior density in exact arithmetic while
keeping density independent. A finite polynomial log field cannot also have an
exact zero-density trace.
For the present vacuum polytrope, the more suitable hypothesis is
`rho = rho_ref * exp(u_h)`, where a reference chosen independently of live
enthalpy is positive inside and vanishes at the surface with the correct
leading exponent. That exponent is a polytropic assumption to test, not a
universal prescription for the future general EOS. Geometry dependencies of
the reference and all mapped operators must be differentiated consistently.
The log-coordinate closure derivative is a density-weighted mass matrix.
Consequently normalization, local solves and globalization must account for
its small response near vacuum. Exponentiation alone does not fix those
conditioning issues, nor does it make the separately represented enthalpy
admissible in every Newton trial. Floating-point underflow also needs an
explicit rejection/range policy; silently clipping density would change the
equations. A typed physical-density evaluator and its directional action would
allow polynomial and log representations to be exchanged internally without
exposing representation details in the physics-facing API.
## Validation
All three builds passed the same 48 focused contract tests with the same RNG
seed. They cover field registration and spaces, quadrature policy, prepared
closure values and derivatives, gravity coupling, physical mass diagonals, and
the normalized complete operator.
| Density order | Test cases passed | Assertions passed |
|---|---:|---:|
| Q2 | 48 | 9,766 |
| Q3 | 48 | 9,822 |
| Q4 | 48 | 9,894 |
The final frozen artifact audit and all three coupled artifact audits pass.
All twelve coupled cases completed with their declared outcome. The audit helper
has nine passing small tests; the comparison helper has seven, including a check
that separately audited orders retain the same mesh and non-density blocks.
The independent log-coordinate control has four passing small tests.
No ESTER comparison or moderately rotating n=3 validation is part of this
study. The n=1,J=0.01 case is a rotation control; its numerical J value is not a
fraction of critical rotation.
## Reproduction and retained evidence
Raw campaigns, immutable diagnostic executables, input/output hashes, build
logs and test logs are retained under
`../../diagnostics/density_order_2026-09-19/` relative to this worktree.
Every campaign records its source commit, mesh hash, complete CMake cache,
arguments and runtime environment. The mesh SHA256 is
`ce11ae99e21e6c3bbfbee402acd2e191c1da0d8261d2227b4203f73f2f337a74`.
Producer commits are `241621c` for frozen controls, `c0e8551` for coupled Q2,
`cd71f37` for coupled Q3 and `e054c74` for coupled Q4. Later changes to reporting
scripts and this document do not replace those immutable executables or their
recorded provenance. The final comparison generator is from `c0a496e`.
The log control's source and CSV hashes were separately verified in
`log_reference_validation.json`.
The committed presets are `density_order_frozen_cases_2026-09-19.json` and
`density_order_coupled_cases_2026-09-19.json` in `experiments/`. Frozen replay
requires the earlier harmonic-surface campaign retained alongside these data;
the runner verifies the saved state, layout, metadata and original provenance.
For example, the Q3 configuration uses:
```sh
/Applications/CLion.app/Contents/bin/cmake/mac/aarch64/bin/cmake -S . -B cmake-build-release-density-q3 -G Ninja \
-DCMAKE_BUILD_TYPE=Release \
-DCMAKE_CXX_COMPILER=/opt/homebrew/opt/llvm/bin/clang++ \
-DCMAKE_C_COMPILER=/opt/homebrew/opt/llvm/bin/clang \
-DCMAKE_MAKE_PROGRAM=/opt/homebrew/bin/ninja \
-DBoost_DIR=/opt/homebrew/anaconda3/lib/cmake/Boost-1.82.0 \
-DMEAN_FIELD_DENSITY_POLYNOMIAL_ORDER=3 \
-DMEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT=0 -DBUILD_TESTING=ON
```
Build `plateau_diagnostic` and `density_order_contract_tests` with two compile
workers. The focused test filter includes all non-normalization cases in that
target plus the reference-mass-diagonal comparison and complete normalized
operator test:
```sh
env OMPI_MCA_btl=self,sm cmake-build-release-density-q3/density_order_contract_tests \
--catch2 --reporter compact --rng-seed 1218739852 \
'~[normalization],Direct Reference Mass Diagonals Match Full True DOF Assembly,Normalized Stellar Equilibrium Preserves Physical Residuals Jacobians And Frozen Solve Coordinates'
```
The same fixed seed is used for the other density orders. Campaign runtime settings are
`OMPI_MCA_btl=self,sm`, `OPENBLAS_NUM_THREADS=1`, `OMP_NUM_THREADS=1` and one MPI
rank. The campaign runner uses a fresh output directory for each invocation;
its `--timeout 1800` applies per case. Use the order-specific build directory
and the coupled preset with `experiments/run_plateau_campaign.py`; all other
arguments are preserved in each `provenance.json` and `*_run.json`.
`experiments/audit_density_order.py` checks the retained evidence independently
of the producer's summary tables. It validates saved polynomial algebra,
bindings, sums, trial histories and declared gates, but cannot reconstruct the
full FEM integrals because local matrices, RHS vectors and all mapped weights
are not exported. `experiments/summarize_density_order.py` consumes only passing
audit JSON files and emits comparison tables and scientific figures. It keeps
undefined values and unresolved sign witnesses explicit.