EQUIL + GPEC - Two-pass measured-curvature auto psi grid with rational-surface knots

[logan-nc]

Closes #302, closes #308. Implements the rational-surface half of #309 (kinetic-resonance knots remain a follow-up; the mandatory-node hook is in place).

Problem

The one-pass auto grid (log_asymptotic, mpsi=0; the auto grid_type is renamed "auto" in this PR, with "log_asymptotic" kept as a working alias) sized knots from a-priori models blind to the formed solution: ~5 edge knots at the default psi_accuracy=1e-3 (spurious dW/Δ′, #302), while tightening the tolerance over-packed the mid-region 4× because N_mid ∝ τ^(−1/2) vs N_edge ∝ τ^(−1/4).

Approach

Two-pass measured-curvature refinement (src/Equilibrium/GridRefinement.jl + driver orchestration):

1. Pass 1 forms the equilibrium on a coarse fixed-128 log_asymptotic layout.
2. A knot density is built from measured nodal data — divided-difference fourth derivatives of the 1D profiles, the rzphi geometry channels along sampled θ-lines (the bicubic ψ-axis shares these knots), and the kinetic profiles when loaded (steep pedestal gradients attract knots, Ballooning: guarantee pedestal resolution in auto-gridding, and make a default psi-window the standard for ballooning scans #308 part 1). Nodal values come from independent field-line integrations, so the estimate is immune to inter-knot spline ringing.
3. Mandatory knots are pinned on every rational surface q=m/n in the requested n range (ForceFreeStates.rational_psi_nodes, extracted from sing_find!), each with local geometric packing; a δ_min snap guard prevents near-duplicate nodes (Kinetic-aware equilibrium grid packing: place knots on rational and kinetic-resonance surfaces #309).
4. Pass 2 re-forms via the existing in-memory rerun path (override_psi_nodes keyword through all solver paths) — no file re-read.

Key modeling decisions (each validated on the DIII-D example):

• h³ derivative error model err(f′) ≈ h³|f''''|/24: the stability physics consumes spline derivatives (q′ at rational surfaces, p′/V′ in the EL and ballooning coefficients). Every region scales as τ^(−1/3), fixing the auto mpsi under-resolves the edge, giving spurious dW and Delta-prime #302 imbalance.
• Curvature measured in ρ=√ψ, where the equilibrium is regular at the axis; the ψ-space geometry channels diverge there and made the implied knot count grow without bound.
• Local packing around rational surfaces (spacing 0.06·τ^(1/3), radius 0.05) converges the Δ′ BVP, which samples the ψ-splined coefficient matrices around each surface — a pinned knot alone left Δ′@q=2 2× off.
• The core below ψ=0.03 uses the a-priori geometric density only (smallest-surface nodal data is integration/extrapolation noise).

Also included:

• Ballooning α-boundary scan window ψ_N ∈ [0.1, min(0.99, ψ_edge)] (hard-coded per the Ballooning: guarantee pedestal resolution in auto-gridding, and make a default psi-window the standard for ballooning scans #308 discussion, NaN-sentinel index mask, locstab path untouched) — 93 of 200 surfaces skipped on the DIII-D case.
• q-vs-ι benchmark (benchmarks/benchmark_q_vs_iota_edge.jl): splining ι=1/q gives a ~1.2–2× constant-factor gain but the same convergence order, documenting that wholesale ι replacement is not warranted for grids ending inside the separatrix (PR EQUIL + ForceFreeStates - Extend equilibrium splines to separatrix for diverted plasmas #179 lineage owns separatrix-extension ι work).

Validation (DIII-D-like example, n=1, defaults)

• 200 knots (vs ~1012 the old model needed for a converged edge): top eigenvalue within 0.01% of ldp-512; Δ′ inner surfaces (q=2,3,4) within ~5%; near-separatrix surfaces q=5,6 match the dense references where ldp-256 was 35% off / wrong-signed — the two-pass grids are internally converged there (−217.5/−217.6/−218.1 across 268/293/454 knots) while ldp is still wandering (+205 → −146).
• fortran-physics-reviewer audit of the density mathematics: PASS (divided-difference table, sharp Hall–Meyer h³/24 constant, ρ-Jacobian, edge-floor algebra dlog=(4τ)^(1/3), equidistribution, merge semantics, rerun dispatch); its two MINOR robustness suggestions are applied.
• Full test suite green, including new test/runtests_grid_refinement.jl (38 assertions: merge invariants, equidistribution, quartic-exact d4, τ^(−1/3) scaling, Solovev two-pass round-trip).
• Regression harness vs develop: solovev_n1 and solovev_multi_n bitwise-unchanged (fixed-grid path untouched). diiid_n1 moves as expected — equilibrium scalars at spline-resolution level (≤1e-5 relative); grid-sensitive rows (energies, Δ′, island widths, resonant field, PE torque) shift because the develop baseline was computed on the defective under-resolved grid; the new values match the dense references above. The diiid_n1 baseline should be re-pinned on merge.
• The DIIID parallel-integration testsets now build the two-pass auto grid exactly as the driver does and pin its values — the auto grid is the production default, so that is what regression guards (near-separatrix Δ′ entries asserted finite-only, per the grid-convergence evidence; the ldp-512 equivalence table is recorded in the PR comments).

Notes for reviewers

• PR Add D3D-like kinetic profiles and update example equilibrium #295's ldp+mpsi=256 workaround can flip back to mpsi=0 auto after this merges (and re-pin its regression).
• Standalone setup_equilibrium callers with mpsi=0 now get the coarse pass-1 grid unless they refine via refined_psi_grid + override_psi_nodes (documented in the docstrings and docs/src/equilibrium.md).
• Kinetic-aware equilibrium grid packing: place knots on rational and kinetic-resonance surfaces #309's kinetic-resonance knots plug into merge_mandatory_nodes/refined_psi_grid(...; mandatory=...) without further Equilibrium changes.

🤖 Generated with Claude Code

[logan-nc]

Knot placement: auto two-pass vs ldp (for the record)

Top: cumulative knot fraction (index/mpsi) vs ψ_N — where each grid spends its budget. Bottom: local spacing Δψ_N (log scale). Dotted verticals are the n=1 rational surfaces. In the top panel the two ldp curves overlap exactly (normalized sin² packing is N-independent). Regenerate with benchmarks/plot_grid_knot_placement.jl.

Reading the figure: the auto grid (blue, 199 intervals) coasts across the smooth mid-radius at ~10× coarser spacing than ldp, dives at every rational surface (the spacing notches), and packs the core, pedestal, and edge — which is exactly where ldp either wastes knots (mid) or runs out of them (near-separatrix rationals).

The grid logic, in plain bullets

• Pass 1: form the equilibrium once on a coarse fixed layout (~128 intervals). Cheap; its only job is to measure the real solution.
• Mid-region (0.03 ≤ ψ_N ≤ 0.9): measured spline error. Take 4th divided differences of the pass-1 node values (profiles F, P, V′, q; the four rzphi geometry channels along θ-lines; kinetic n/T/ω_E when loaded) and size the local spacing so the cubic-spline derivative error stays below psi_accuracy (τ). Node values come from independent field-line integrations, so this measurement can't be fooled by inter-knot ringing. Curvature is measured against ρ=√ψ, where the geometry is regular at the axis.
• Core (ψ_N ≤ 0.03): a-priori only — NOT the same as the mid. The smallest flux surfaces' nodal data is integration/axis-extrapolation noise, so measurement is not trusted there; the density is the geometric-in-log(ψ) form outright.
• Edge (ψ_N ≥ 0.9): measured, with an a-priori floor — also not the same as the core. The measured density applies, but is floored by geometric-in-log(1−ψ) packing so a pass-1 undersampled edge can't hide. The floor constant comes from uniform relative q′ error and is independent of the separatrix slope A.
• Rational surfaces: pinned + locally packed. Every q=m/n in the requested n range gets a knot exactly on the surface plus geometric packing around it (this, not the pinned knot alone, is what converges Δ′).
• Equidistribute the resulting density into the final knot vector; pass 2 re-forms the equilibrium on it from the in-memory input. Every region scales as τ^(−1/3), so one knob refines everything proportionally.

Rename: log_asymptotic → auto

Agreed the old name misdescribed the method — there are no q→∞ asymptotics anywhere on the grid (all finite q; the log-form survives only as the edge/core floor shape). grid_type = "auto" is now the default; "log_asymptotic" remains a working alias (two-pass when mpsi=0, the three-region log layout when mpsi>0), so no existing TOML breaks.

Reproduces ldp mpsi=512, and test pinning

Convergence evidence on this example (n=1; Δ′ = BVP diagonal, PEST3 convention as stored in gpec.h5):

grid	knots	et[end]	Δ′(q=2)	Δ′(q=3)	Δ′(q=4)	Δ′(q=5)	Δ′(q=6)
auto, τ=1e-3 (default)	200	0.078112	8.47	−4.03	−10.23	12.35	−144.9
auto, τ=1e-4	454	0.078086	8.98	−3.86	−10.23	12.35	−217.5
ldp, mpsi=256	256	0.077800	8.79	−3.88	−10.21	16.62	+205.2
ldp, mpsi=512 (reference)	512	0.078104	8.90	−3.46	−10.16	12.23	−146.1

• The default auto grid matches ldp-512 on the eigenvalue to 0.01% and on Δ′ to ~5% on the inner surfaces — at 200 knots vs 512. On the near-separatrix surfaces it agrees with ldp-512 where ldp-256 is 35% off (q=5) or wrong-signed (q=6). Tightening τ converges Δ′ further (τ=1e-4 column).
• Test pinning is now on the auto grid (commit 42d9cf4): the DIIID parallel-FM and STRIDE-BVP testsets build the two-pass grid exactly as the driver does (rational_psi_nodes → refined_psi_grid → ingest re-form) and pin its values — the auto grid is the production default, so that is what regression must guard. The table above is the recorded proof of ldp-512 equivalence; the near-separatrix Δ′ entries are asserted finite-only since they are not grid-converged in any grid family at these counts.

🤖 Generated with Claude Code

[logan-nc]

Follow-up: default-convergence scan, core ψ→0 handling, and edge-logic assessment

Should the default converge everything? Scan says: no default change is justified

Full-pipeline scan over (psi_accuracy, sing_pack_coef) on the DIII-D example (identical BVP path for every row; anchor row densest):

τ	pack coef	knots	et[1]	Δ′(q=2)	Δ′(q=3)	Δ′(q=4)	Δ′(q=5)	Δ′(q=6)
5e-5	0.02	574	1.31456	9.80	−3.15	−10.18	11.86	−206.2
1e-3	0.06 (default)	199	1.31163	8.46	−4.04	−10.24	11.70	−138.3
1e-3	0.04	214	1.31147	8.44	−4.01	−10.13	11.63	−140.3
1e-3	0.03	221	1.31156	8.71	−3.40	−10.31	11.63	−140.2
1e-3	0.02	235	1.31147	9.25	−3.45	−10.16	11.76	−137.0
5e-4	0.03	272	1.31291	9.30	−3.09	−10.25	11.34	−201.7
3e-4	0.03	313	1.31523	9.13	−3.28	−10.22	11.48	−194.6

What the table shows:

• Converged at the 199-knot default: eigenvalue (0.2% of the 574-knot anchor), Δ′(q=4) (±0.5%), Δ′(q=5) (±2%).
• Not convergeable by grid density: Δ′(q=2) wanders 8.44 → 9.80 (±8%) and Δ′(q=3) −4.04 → −3.09 (±13%) without settling — adjacent settings (221 vs 235 knots) flip q=2 by 6%. The wandering is correlated with the q=6 edge state: rows with Δ′(q=6) ≈ −140 form one branch, rows with ≈ −200 another. The residual inner-Δ′ spread tracks the edge-truncation branch (where the ψ-limit/dW-peak lands for that grid), the same mechanism as the long-known q=6 bistability — it is not a knot-count effect, so spending 300–600 knots buys nothing.
• Conclusion: the lean τ=1e-3 default stands; the existing test pins (exact default grid) remain the right regression guard. The Δ′-vs-edge-truncation sensitivity deserves its own issue (BVP outer boundary condition / psilim landing), separate from gridding — happy to file it.

sing_pack_coef is now exposed as a keyword on refined_psi_grid for future convergence studies (default unchanged).

How does the auto grid treat the bicubic splines blowing up near ψ→0?

As a coordinate artifact, which is what it is: the rzphi channels behave as ψ^(k/2) near the axis (R−R₀ ~ √ψ), so their ψ-derivatives diverge under refinement while their ρ=√ψ derivatives converge. The grid (a) measures curvature in ρ, (b) converts density with the dψ/dρ = 2√ψ Jacobian (√ψ-tight packing toward the axis by itself), and (c) below ψ=0.03 switches to pure geometric-in-log(ψ) packing — constant-ratio spacing ∝ ψ down to psilow (~37 knots at psilow=1e-4), i.e. yes, asymptotically tight — because nodal data from cm-scale flux surfaces is integration/extrapolation noise that refinement amplifies rather than resolves. Rational-surface pinning+packing still applies at any ψ, including inside the core region.

Is the different edge treatment needed? Honest assessment

Mostly no — one consistent logic already governs, and the docs now say so (commit 053526a). The measured ρ-space density supplies ~3× more edge knots than the edge model demands on this example (84 vs 29 in ψ>0.95), because the pass-1 layout is already log-packed at the edge — so the edge "floor" is inactive insurance, not a distinct treatment. The honest framing, now in the module docstring and docs/src/equilibrium.md: one sizing rule (cubic derivative error h³|f''''|/24 ≤ τ·scale) applied to three curvature sources — measured divided differences (mid + edge), the analytic separatrix log model (edge floor, normally inactive), and the analytic power-law axis model (core, where it replaces measurement because tiny-surface nodal noise is irreducible in a two-pass scheme). The core is the only place a genuinely different rule applies, and that is forced by data, not aesthetics.

🤖 Generated with Claude Code

[logan-nc]

"Curvature of what?" — pedestal-packing attribution, with data

Definition: the density measures |d⁴f/dρ⁴| (ρ=√ψ_N) — the quantity in the cubic-spline derivative error bound err(f′) ≈ h³|f''''|/24 ≤ τ·scale — for each of: the 1D profiles F = R·B_φ, μ₀·P, dV/dψ, q; the four rzphi bicubic channels (r², η-offset, ν, Jacobian) along every 8th θ-line (worst case over θ); and the kinetic profiles (nᵢ, nₑ, Tᵢ, Tₑ, ω_E) when loaded. The final density is the max over all of them. Note the develop example loads no kinetic profiles (no [KineticForces] section until #295 merges), so the flagship's pedestal packing rests on the equilibrium quantities alone today; the kinetic n/T/ω_E channels will add their own pedestal pull after #295.

Is the pressure pedestal gradient appropriately included? Yes — verified by per-source decomposition rather than asserted:

Band integrals (knots contributed per band if that source acted alone), DIII-D example, τ=1e-3:

source	[0.30,0.85]	[0.85,0.90]	[0.90,0.95]	[0.95,0.9995]
F	0.0	0.3	0.6	2.1
P	5.1	4.2	3.6	16.9
dV/dψ	2.4	3.9	3.6	53.5
q	5.1	3.0	3.5	55.3
rzphi r²	6.8	3.6	3.2	65.2
rzphi offset	16.3	5.0	4.2	76.6
rzphi ν	9.1	4.3	3.4	66.4
rzphi jac	2.4	3.9	3.6	53.5
floors (core/edge/min)	2.8	0.3	1.2	32.5
rational packing	25.6	9.0	4.7	12.6

Reading:

• The pedestal band is measured-curvature-driven, not floor-driven. At every node in ψ∈[0.8, 0.98] the winning source is measured (the η-offset geometry channel, 62–588 knots/unit-ψ), well above the edge floor (5–315). The floor only matters as insurance.
• P contributes materially and comparably to q and dV/dψ (~8 knots across [0.85,0.95] on its own; h≈0.012 pedestal spacing from P alone, scaling as τ^(−1/3)). The geometry channels contribute more — which is physically right: the offset/Shafranov-shift steepening is the Grad-Shafranov response to that same pedestal p′, and its splines live on the same knot vector. The pedestal is therefore multiply covered, and since everything is measured from the formed solution, the packing follows the pedestal wherever it sits — no hardcoded 0.9–0.95 band anywhere.
• One flagged (theoretical) weakness stays on record: f_scale = max|f| normalization suppresses large-offset/low-variation quantities — visible as F's tiny contributions above. Harmless here because F's pedestal information is carried redundantly by P/q/V′/geometry; if a future case ever shows a pedestal that only F sees (it shouldn't, physically), range-based normalization is the one-line fix.

Figure + one-paragraph summary committed to the docs (a145a47).

🤖 Generated with Claude Code

[logan-nc]

@d-burg it seems to me that all the auto grid logic checks out. Default settings give ~200 points for DIIID-like examples. I adde dsub issues to yours to make sure the logic captures the pedestal region well for ballooning and kinetic resonance regions well for kinetic terms. The logic also now explicitly knots on rationals, which I think should be a big win for all applications but I hope for you in particular.

Please test on the cases that prompted #302 and confirm all is good. If yes, then merge away.