MeanField/experiments/POLYTROPE_VALIDATION.md

Nonrotating n=1 polytrope verification

This experiment separates nonlinear stopping criteria from physical accuracy. It uses a closed-form reference independent of the production Lane–Emden seed, then measures the actual mapped mesh and, optionally, the last accepted production state. It does not modify the physical equations, nonlinear policy, mesh, or existing sandbox.

The target is polytrope_validation_experiment, declared EXCLUDE_FROM_ALL. Build only that target:

cmake --build cmake-build-release-homebrew --target polytrope_validation_experiment -j 6

The executable requires one MPI rank and uses the CPU device. The following sequence is deliberately bounded: analytic self-checks, an analytic-field mesh check, then one production solve. Choose fresh output directories; existing directories are refused.

./cmake-build-release-homebrew/polytrope_validation_experiment --self-check --output polytrope_self_checks
./cmake-build-release-homebrew/polytrope_validation_experiment --analytic-mesh --mesh sandbox.smesh --output polytrope_analytic_mesh
./cmake-build-release-homebrew/polytrope_validation_experiment --solve --mesh sandbox.smesh --output polytrope_solution

Inspect each stage before proceeding. The commands are usage examples, not a claim that the corresponding runs passed. Results belong in a separate findings report. No full test suite or parameter/resolution sweep is part of this workflow.

For an existing completed output with profiles, the standard-library-only summarizer can generate a static SVG and Markdown report without rerunning the solver or changing its CSVs:

python3 experiments/summarize_polytrope_validation.py polytrope_solution --control-directory polytrope_analytic_mesh

One additional uniform h-refinement

Use the saved coarse input, not a newly generated mesh, and retain a separate refined snapshot. The helper calls UniformRefinement(loaded_mesh, 1), verifies eight children per hexahedron in each material region, and checks saved-mesh reload and input immutability. It records mesh and STROID extension hashes.

/opt/homebrew/anaconda3/envs/stroidDev/bin/python3.14 experiments/refine_polytrope_mesh.py polytrope_solution_2026-09-09/input.smesh polytrope_h1_mesh
env OMPI_MCA_btl=self ./cmake-build-release-homebrew/polytrope_validation_experiment --analytic-mesh --mesh polytrope_h1_mesh/refined.smesh --output polytrope_h1_analytic
env OMPI_MCA_btl=self ./cmake-build-release-homebrew/polytrope_validation_experiment --solve --mesh polytrope_h1_mesh/refined.smesh --output polytrope_h1_solution
python3 experiments/compare_polytrope_refinement.py polytrope_solution_2026-09-09 polytrope_h1_solution --coarse-control polytrope_analytic_checked_2026-09-09 --fine-control polytrope_h1_analytic --output polytrope_h_comparison

Inspect the analytic control before launching the solve. Keep field orders, model, quadrature and solver tolerances fixed. The full refined solve may take well over an hour; reuse the coarse solve rather than repeating it.

The Python environment above contains multiblock-capable STROID 0.5.0. The current MeanField build links an older /usr/local STROID that can load the saved geometry but cannot correctly regenerate multiblock geometry from its configuration. Therefore refinement is performed by the current Python binding, not by that old linked library. No installation or relinking is needed. The helper fails if core_mapping support is absent.

STROID refines the logical mesh and reprojects the high-order geometry; this is not subdivision of a fixed physical polynomial geometry. Consequently the comparison includes geometry approximation as well as FE field refinement. Run the analytic control on both levels to quantify the geometry contribution. The refined file is complete: the experiment uses extraRefine=0, and its ordinary input.smesh/GF snapshots can be replayed without additional refinement or the newer STROID library.

Two-level ratios E_coarse/E_fine and log2(E_coarse/E_fine) are observed reductions, not proof of asymptotic order. Prioritize stellar volume field errors, shape and virial balance. Exterior profile accuracy is diagnostic only, relevant insofar as it affects the interior solution. Near-zero integral errors and errors near the nonlinear/quadrature floor do not give reliable h-rates.

Modes and controls

Option	Meaning/default
--self-check	Check the independent reference without loading a mesh or constructing a solver.
--analytic-mesh	Evaluate exact fields at physical points on the actual mesh; exercise volume/surface integration and radial location without a Newton context.
--solve	Default mode. Construct the production n=1 problem, solve within the stated budget, and measure the last accepted state even after nonlinear failure.
--replay DIRECTORY	Re-measure saved nonrotating solve-mode fields using that directory's input.smesh; no Newton context or solve. Requires a new --output directory.
--mesh FILE	Default sandbox.smesh; input is not overwritten.
--output DIRECTORY	Default polytrope_validation_results; must not already exist. Its parent directory must exist.
--absolute-tolerance VALUE	Nonlinear absolute tolerance, default 1e-8.
--relative-tolerance VALUE	Nonlinear relative tolerance, default 1e-8.
--linear-tolerance VALUE	Relative linear tolerance, default 0.03.
--max-newton N	Default 8 nonlinear iterations.
--max-linear-iterations N	Default 80 linear iterations per solve; FGMRES restart length is 40.
--quadrature-order N	Base physical-volume quadrature order, default 14.
--check-quadrature-order N	Independent higher order, default 18; must exceed the base order. Also used for surface integration.
--skip-profiles	Omit physical-shell/ray location and its output; other measurements remain enabled.
--mu-points N, --phi-points N	Angular quadrature sizes, default 6 and 12. Increase during cheap replay to check spherical-mean sampling sensitivity.
--exterior-shells N	Number of finite-exterior radii between 1.001R and 2R, default 8. Increase during replay to resolve radial structure.

Select one mode explicitly. The parser accepts the last mode flag if several are provided. Every mode first runs the independent analytic self-checks.

Exit codes are 0 for all requested checks passing, 1 for nonlinear failure, 2 for an execution/input error, and 3 for verification failure. In solve mode, nonlinear failure takes precedence over a subsequent physical-screen result; consult the saved metrics and metadata as well as the exit code. Replay retains that precedence using the source's historical convergence status.

--analytic-mesh is not an FE projection test or a numerical equilibrium solution. Exact fields are evaluated directly, so volume field-versus-reference errors are zero by construction at the same successfully evaluated points. Profile comparisons instead use the requested physical radius and therefore also test the locator's accuracy. The mode checks integration, physical location, reference-domain geometry, and integral identities on that domain.

Diagnostic replay without another solve

After a solve-mode output has saved its grid functions and completed physical_metrics.csv, physical postprocessing can be repeated independently:

./cmake-build-release-homebrew/polytrope_validation_experiment --replay polytrope_solution --output polytrope_solution_replay

Replay requires an original mode=solve output for this nonrotating n=1 benchmark, matching compiled (G,M,R), and saved angular velocity exactly zero. It loads the saved mesh and the five GF files into compatible FE spaces. The saved accepted coefficients need not have passed the physical screen, but all required source artifacts must be present. Keep the original solve output: a replay output is not itself an accepted replay source.

Volume/surface metrics and requested profiles are recomputed. Solver convergence, bordered/unbordered residual norms, central-border action/value, Bernoulli constant, and angular-velocity norm are historical source diagnostics, not recomputed. Metadata records solver_diagnostics=copied_from_source_not_recomputed. In particular, a replay pass is not a new residual evaluation or convergence claim, and solver-tolerance options do not trigger a fresh solve. The source's nonlinear failure still produces exit code 1 after successful physical postprocessing.

Fixed analytic reference and normalization

The benchmark uses the compiled utils::G, utils::MASS, and utils::RADIUS as fixed positive (G,M,R); their values are saved in metadata.txt. It does not fit mass, radius, central density, or a potential offset to the numerical solution. For (\xi=\pi r/R), the stellar solution is

[ \theta=\frac{\sin\xi}{\xi},\qquad K=\frac{2GR^2}{\pi},\qquad \rho_c=\frac{\pi M}{4R^3},\qquad h_c=\frac{GM}{R}, ] [ \rho=\rho_c\theta,\qquad h=h_c\theta,\qquad P=K\rho^2,\qquad \Phi=-h_c(1+\theta),\qquad m(r)=\frac{M}{\pi}(\sin\xi-\xi\cos\xi). ]

The radial potential gradient is outward-positive (g_r=d\Phi/dr=Gm(r)/r^2); the acceleration is its negative. Outside the star the analytic material fields are zero, (\Phi=-GM/r), and (g_r=GM/r^2). The implementation uses origin series and a small-distance-to-surface expression, including exact values at the center and surface. Reference radii must be finite and nonnegative.

The profile normalizations are fixed:

[ \theta_\rho=\rho/\rho_c,\qquad \theta_h=h/h_c,\qquad \theta_\Phi=-R\Phi/(GM)-1. ]

All three agree with (\theta) inside the analytic star. The normalized vacuum potential is negative outside (R), approaching (-1), and is not clipped. Numerical density and enthalpy are likewise never clipped. Negative samples and their minima are reported. Since (P_\rho=K\rho^2) and (P_h=h^2/(4K)) are positive even for negative arguments, pressure checks alone do not establish physical positivity.

The independent integral references are

[ \Pi=\int P,dV=\frac{GM^2}{4R},\qquad W=-\frac{3GM^2}{4R},\qquad I_z=\frac23\left(1-\frac6{\pi^2}\right)MR^2. ]

The self-checks use independent radial Simpson integration, including nonunit scales, origin regularity, surface/vacuum joins, EOS and hydrostatic identities, and the gravitational-field energy with its exterior contribution.

Physical balances and interpretation

The volume measurements integrate over the current stellar material elements with the full physical Jacobian. They compare fields at their actual physical positions, not at logical radii. Relative volume (L^2) errors use the analytic field's (L^2) norm over that same domain. Surface and volume-equivalent radius errors separately measure the domain discrepancy.

Let (g) denote the reconstructed physical mixed gravity field and (\Psi=|\Omega\times x|^2/2). The reported energies are

[ T=\int\rho\Psi,dV,\quad W_\Phi=\tfrac12\int\rho\Phi,dV,\quad W_g=-\int\rho,x\cdot g,dV,\quad \Pi=\int P_\rho,dV. ]

The scalar virial error is (|2T+W_\Phi+3\Pi|/|W_\Phi|); virial_signed retains its sign and virial_ratio is ((2T+3\Pi)/|W_\Phi|), which should approach one. The force virial replaces (W_\Phi) by (W_g), keeping the same denominator. gravity_energy_consistency measures (|W_\Phi-W_g|/|W_\Phi|). These are complementary checks: an inaccurate potential and an inaccurate mixed gravity field need not fail identically.

The central-density constraint is imposed through central enthalpy and an additional hydrostatic border. Its reported achieved density is inferred from that enthalpy; it is not an independently sampled DG density value. The unbordered residual removes only the artificial central-border action from hydrostatic rows, then uses the same production normalization. It retains the physical Bernoulli constant and all scalar constraint rows. A small bordered solver residual does not by itself establish a small unbordered physical residual.

The potential is discontinuous across elements. Its elementwise, or broken, gradient omits interface jumps; it is not the same discrete object as the mixed H(div) gravity field. Their reported gradient mismatch is a diagnostic, not a requirement of pointwise equality at finite resolution. The strong enthalpy-gradient balance is also distinct from the weak assembled residual.

Bernoulli statistics concern (h+\Phi-\Psi). Both its mean error against the fixed analytic constant (-h_c) and its spatial variation are saved. Computing a centered variance does not fit or subtract a potential gauge from the fields.

Weak EOS closure and the pressure projection floor

The default density space is DG/L2 order 2, while enthalpy is continuous H1 order 3; see the field registry. The prepared closure assembles (F_i=\int q_i(\rho-h/(2K)),dV) at (n=1), using the full physical volume weight and density-space test functions. Thus, at fixed geometry and zero closure residual, density is the quadrature-weighted (L^2) projection of (h/(2K)), not necessarily its pointwise value. The seed's independent coefficient projections do not themselves enforce this orthogonality.

The pressure-force kernel uses (P(h)=h^2/(4K)), rather than the diagnostic's primary (P(\rho)=K\rho^2); these formulas follow from the polytropic EOS for admissible nonnegative enthalpy. With (\delta=h-2K\rho), define

[ E=\int\frac{\delta^2}{4K},dV,\qquad D=\Pi_h-\Pi_\rho. ]

The algebraic identity, using the same domain and quadrature, is

[ D=E+\int\rho\delta,dV,\qquad D-E=-2K,\boldsymbol{\rho}^{,T}\mathbf F_{\rm closure}. ]

Consequently, exact weak closure gives (D=E\ge0), even when the strong EOS mismatch (\delta) is nonzero. A pointwise mismatch can therefore persist at a well-converged discrete solution without indicating an EOS/Jacobian algebra bug. Changing diagnostic quadrature introduces an additional discrepancy in the residual-pairing identity and should be checked separately.

The observer records closure_projection_pressure_gap ((E)), closure_density_inner_product ((\int\rho\delta,dV)), and closure_projection_pressure_gap_relative_defect (((D-E)/\Pi_{\rm reference})). Older metric files also allow reconstruction of (E=Vh_c^2,\text{eos_enthalpy_scaled_rms}^2/(4K)). Inspect these alongside the complete barotropic_closure residual block: a small single global pairing can conceal cancellation and does not prove every closure equation is satisfied.

enthalpy_virial_error and enthalpy_force_virial_error substitute (\Pi_h) into the two virial diagnostics. They complement, not replace, the original density-pressure checks. The predetermined pointwise EOS screening budget is not relaxed: it may identify a finite-resolution projection floor requiring further discretization study. Diagnostic replay can produce the additional metrics from saved fields without another Newton solve.

Shell and ray sampling

Default profiles use the origin, 32 interior radii through (0.99R), one additional radius at (0.999R), and eight exterior radii from (1.001R) to (2R). Each noncentral sphere uses six Gauss points in (\mu=\cos\vartheta) and twelve uniformly spaced azimuths. The weights sum to one; angular moments and exact constant-field weighted mean/variance are self-checked before measurement. Separate 26-ray profiles contain six axes, twelve face diagonals, and eight body diagonals. Rays are not angular quadrature.

The locator inverts the full physical mapping. Sampled element bounds prioritize searches but do not exclude elements; missing points may therefore be expensive. Location error, attempts, and coverage are reported.

Shell means are conditional on successfully located finite values. Density and enthalpy are additionally conditional on stellar material. In particular, their means near a displaced surface are not whole-sphere density/enthalpy averages: inspect material and valid-weight coverage. Missing/exterior material values are not silently replaced by zero. Angular spread and total RMS error against the fixed radial reference are both written. Axis/diagonal DG values can be one-sided traces at element interfaces. The origin is explicitly a single trace, not an angular average; a radial gravity component is undefined there.

Output files

File	Contents
metadata.txt	Mode, input mesh/snapshot, compiled scales/settings, FE orders, and available solve/measurement status and timing.
input.smesh	Exact input-mesh copy used to construct the FE spaces in each mesh-loading mode, including replay.
analytic_self_checks.csv	Independent check, observed/expected values, fixed scale, errors, tolerance, pass flag.
volume_metrics_base.csv	Physical-volume metrics at the base quadrature order.
quadrature_comparison.csv	Base/check values and their absolute/relative changes. Relative changes of nearly zero metrics require caution.
physical_metrics.csv	Higher-order volume metrics, surface/corner metrics, optional profile summaries, and available solver/border diagnostics.
seed_physical_metrics.csv	Production seed's lower-order volume metrics and residuals, measured before Newton to expose any physical degradation during correction.
verification_checks.csv	Declared screening budgets, observations, and individual pass flags.
radial_profiles.csv	Physical-shell means, angular spreads, analytic values, scaled errors, coverage, and normalized profiles.
directional_profiles.csv	Axis/diagonal values, analytic values, element/material identity, and location errors.
newton_history.csv	Solve-mode iteration residuals, accepted steps, trial counts, and linear/nonlinear diagnostics.
residual_blocks.csv	Solve-mode physical and normalized block norms, with and without the central border.
field_reconstruction.csv	Reduced/full field sizes and coefficient round-trip errors.
accepted_state.txt, state_layout.csv	Solve-mode last accepted coefficient vector and block layout.
density.gf, enthalpy.gf, potential.gf, gravity_gradient_reference.gf, displacement.gf	Solve-mode reconstructed grid functions saved before expensive postprocessing.
polytrope_profiles.svg, polytrope_summary.md	Optional summarizer products, generated from existing CSVs; may be regenerated independently.

The text/GF data are experiment artifacts, not a production checkpoint or a supported solver-restart format. Their supported reuse is the constrained diagnostic replay described above, not continuation of Newton iterations. Preserve the original solve output, saved input mesh, and metadata. In particular, gravity_gradient_reference.gf contains the reference-mesh Piola representation, not the already transformed physical gravity field; its interpretation requires the corresponding mapping and displacement. Unsupported material coefficients in saved grid functions are not numerical vacuum data.

Screening budgets and limits

Budgets are declared in the driver before solving. The default screen requires:

• Relative mass, radius, density/enthalpy/potential/gravity (L^2), binding energy, pressure integral, and axial inertia errors at most 1e-4.
• Bernoulli scaled RMS variation at most 1e-4.
• Scalar/force virial errors, gravity-energy disagreement, and EOS enthalpy scaled RMS at most 1e-6.
• Virial quadrature change and binding-energy relative quadrature change at most 1e-8; normalized unbordered residual at most 1e-8 when available (fresh in solve mode, historical in replay).
• Maximum negative density/enthalpy excursions divided by their fixed central scales at most 1e-8; the underlying negative values are not altered.
• Zero invalid stellar corner samples and, when enabled, zero missing profile points; kinetic energy at most 1e-14 in the compiled benchmark units.
• In --analytic-mesh mode with profiles enabled, maximum scaled pointwise profile error at most 1e-8 and maximum scaled angular RMS at most 1e-10 as independent location and spherical-scatter control gates.

All raw metrics remain available, including quantities without pass/fail budgets. The corner checks include vertices and nearby interior points but cannot certify positivity everywhere in a high-order element. Two quadrature orders test integration sensitivity, not spatial-discretization convergence.

A passing single-resolution screen is not a physical convergence certificate; a small nonlinear residual is not one either. Future work should separate quadrature, mesh/order, and nonlinear-tolerance errors using planned, bounded resolution studies, rather than starting broad sweeps automatically.

The user provided an independent ESTER executable at /Users/tboudreaux/Programming/ESTER_polytrope.pub/ester and the command ester < dati_ester_polytrope. This is a future cross-validation path only: ESTER is not integrated or run by this experiment. Resolve the relative input in its intended project working directory, and reconcile units, boundary conditions, rotation, potential gauge, and diagnostic definitions before a future comparison.