Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots

[logan-nc]

Summary

Each flux surface was traced independently and then splined on that surface's own solver-chosen abscissae before being resampled onto the common θ grid (DirectEquilibrium.jl). The resample error was therefore uncorrelated between neighbouring surfaces — white noise in ψ that grid refinement amplifies rather than reduces. Neither reltol nor abstol touched it, because it is remap interpolation error, not integration error.

The trace now returns its dense solution, and equilibrium_solver root-solves (Brent, bracketed by the monotone jac-weighted flux integral) for the angle at which the normalised straight-fieldline angle reaches each target node, evaluating there. Every surface is sampled at identical abscissae and the resample error at the output nodes is zero. The arclength tracer returns nothing and keeps the previous path.

Investigated under issue #376.

Review round (2026-10-03)

Sampling is now part of the tracer contract. Every field-line tracer takes the target
straight-fieldline angles and returns its surface already sampled there, through one shared helper
(sample_trace_at_sfl_angles). equilibrium_solver has a single path. Before this round the fix
covered only direct_fieldline_int: efit_arclength still resampled each surface from its own
adaptive steps and kept the noise, and the sol === nothing branch would have given any future
tracer the same silent fallback.

Pre-existing bug fixed: F and P were read at the solver's last evaluation, not on the surface.
direct_fieldline_int returned its ODE scratch DirectBField, which the right-hand side, the
refinement callback and the dense interpolant overwrite on every call. The arclength tracer already
kept them separate. The correction is ≤1.5e-8 in P (Solovev, at the edge) and ~1e-10 in q and ν.

The restructure is otherwise a pure restructure. With the F/P fix applied to both the old
sampling and the refactored one, the harness is bit-identical on all six cases. Every number
below that moved relative to the pre-review head is therefore attributable to the arclength fix or
the F/P fix.

New test, runtests_sfl_sampling.jl, checks the property rather than the mechanism, and was
mutation-checked against the code it guards:

code	Solovev rzphi, mpsi 64 vs 128	efit ↔ efit_arclength, same grid
develop (old resampling, both tracers)	fail, 2.1e-4 to 1.1e-3	fail, 3.9e-4
this branch before the arclength fix	pass	fail, 2.4e-4
this branch	pass, ≤ 2.2e-7	pass, 2.3e-6

The remaining ~2e-6 tracer disagreement is the arclength trace's ordinary integration error. It
follows the level set by construction, and at its configured position tolerance it sits within ~1e-6 m
of the flux surface. That is the floor of the agreement test; tightening it would mean tightening
that tolerance.

Merged in: #414 (near-axis kinetic validity suppression), which also cures the DIII-D kinetic
step regression listed below. Relation to #480: measured, and it does not make this PR
redundant: comment.

Why Solovev's energy moves 98%: an independent check (2026-10-04)

Solovev is analytic, but GPEC does not trace the formula. It samples ψ on an (R, Z) grid, splines it,
and traces the spline. On develop each traced surface was then resampled from its own solver
steps. On Solovev_ideal_example, which runs at etol = 1e-7, the 9th-order solver crosses a whole
surface in 23–29 steps. A cubic spline through ~25 points, resampled at 257 angles, is off by
2–4e-3 in r² compared with a densely traced surface (measured). The error differs from surface to
surface, so the ψ-derivatives in the Euler-Lagrange coefficients amplify it. et[1] is a
near-cancellation (ep ≈ −10.4, ev ≈ +10.4), so a 6.5 % plasma-energy error became a 98 % et[1] error.

This check does not depend on this PR's code. Tightening etol makes develop's steps denser and walks
it toward this PR's answer, and this PR's answer does not depend on etol:

etol	develop et[1]	develop ep	this PR et[1]	this PR ep
1e-7 (the deck)	0.6773	−9.7361	0.010735	−10.41668
1e-10	0.011719	−10.41523	0.010735	−10.41668
1e-13	0.011520	−10.41590	0.010735	−10.41668

Develop is still 7 % high at 1e-13 because it still takes only 40–70 steps per surface. Correction to
earlier text in this description: the removed error was called "~1e-6". On this deck it is ~1e-3. On
DIII-D the shift is only 0.1 % because that deck runs at etol = 1e-10, where the steps are denser.

@d-burg, this is likely relevant to your work:

• Tearing and Δ′ studies on traced equilibria. On develop, any directly traced equilibrium (sol,
  tj_analytic_direct, efit) carried this etol-dependent geometry error, at roughly the scale
  above. That includes the Cerfon–Freidberg equilibrium in Equilibrium - FEATURE - 🌱 Add the Cerfon-Freidberg diverted analytic equilibrium #495, which returns a DirectRunInput.
  Δ′ depends on exactly the ψ-derivatives that amplify it. The inverse paths (lar, tj_analytic,
  CHEASE) do not have this mechanism.
• Equilibrium - BUGFIX! - 🚨 Tie the equilibrium ODE absolute tolerance to etol #491 changes abstol on the two tracer solves this PR restructures. On develop, part of the
  numerical movement from a tolerance change comes through this resample sparsity rather than through
  integration accuracy. After this PR, step density no longer enters the sampled geometry. The two
  PRs also overlap textually on those solve lines, though the resolution is mechanical.

The headline: the answer stops depending on the grid

Solovev, same commit, only this change differing — the free-boundary energy was not converging before and is now converged to 6 significant figures:

mpsi	et[1] develop	et[1] this PR	EL steps develop	EL steps this PR
256	1.959e-02	1.462068e-02	1074	775
512	5.399e-02	1.462087e-02	1558	806
1024	8.890e-02	1.462084e-02	2689	891

develop's et[1] drifts 4.5× across the ladder and is still moving; the plasma energy drifts with it (−10.4049 → −10.3144). With this change both are grid-invariant, and the step count is nearly flat (1.04×, 1.11× per doubling versus 1.45×, 1.73×).

That is the property we want: once the equilibrium splines resolve the equilibrium, adding knots should change neither the answer nor the work.

Regression harness

regress --cases diiid_n1,solovev_n1 --refs develop,local — original pre-port measurement, baseline develop @ 9491f89. The current numbers, re-run against develop @ 349a0c262, are in the ## Regression report section below; the mechanism discussion here is unaffected.

diiid_n1 — physics moves only in the 3rd–4th digit, cost halves:

quantity	develop	local	diff
total energy Re(et[1])	8.012318e-01	8.004564e-01	0.10%
plasma energy Re(ep[1])	-1.348486e+00	-1.349677e+00	0.09%
vacuum energy Re(ev[1])	2.149718e+00	2.150134e+00	0.02%
plasma volume	1.829472e+01	1.829471e+01	0.00%
ODE steps (total)	4572	1974	−56.8%

q0, q95, beta_t, beta_n and the singular-surface locations/count are unchanged to 0.00%.

solovev_n1 — et[1] moves 98%. This is the fix working, not a regression: as the table above shows, the develop value is grid-dependent and non-convergent. Solovev's ν node data sat at ~19% relative white noise under the old resample, so it is the case this change corrects most. et[1] is also a near-cancellation (ep ≈ −10.4, ev ≈ +10.4), so a 7% plasma-energy correction dominates it.

On solovev_n1, q0, q95, the singular-surface count/locations/q-values and mpert are all
exactly unchanged; beta_t/beta_n move 0.02%/0.08%.

Runtime 210.6s → 184.5s.

Kinetic cases — large moves, all of which enter at this PR

regress --cases solovev_kinetic_calculated,solovev_kinetic_ntv,solovev_kinetic_nuzero --refs develop,local. This PR moves the kinetic harness values, some of them dramatically — that is owned here, not hidden:

case	quantity	develop	local	diff
kinetic_calculated	Re(et[1])	1.845905	1.799588	2.5%
kinetic_calculated	Im(et[1])	−1.421059	−1.309210	7.9%
kinetic_ntv	NTV torque (fgar), magnitude	~1.56e-4 N·m	~3.85e-1 N·m	~2500×
kinetic_ntv	ψ-quadrature integrand evaluations	840	60	−93%
kinetic_ntv	root-area-weighted Re(et[1])	6.773e-1	1.074e-2	−98%
kinetic_nuzero	Im(et[1])	−2.034869	−1.559010	23%
kinetic_nuzero	Re(et[1])	2.245440	2.284431	1.7%

Attribution: the stacked follow-up branches (EL grid cap + certified kinetic grid, knob off) reproduce these local values bit-identically (34/34 tracked quantities unchanged vs this branch's head), so every kinetic delta vs develop enters at this commit — none at the follow-ups.

Interpretation: the old per-surface resample injected white-in-ψ geometry error (~1e-3 on the Solovev deck at etol = 1e-7; see "Why Solovev's energy moves 98%" above), and resonance-dominated kinetic quantities amplify exactly that kind of noise — the Solovev cases are the most sensitive in the suite. The quadrature-cost drop (840 → 60 evaluations for the same tolerance) is direct evidence the develop-side torque integrand carried noise structure the quadrature was chasing. The eigenvalue moves (2–8%, 23% for Im at ν→0) are consistent with the ideal-case finding that develop's Solovev values were grid-dependent while this branch's are grid-convergent.

---

Second change, folded in from #408: cap the core knot density of the EL coefficient splines

(#408 was folded here on 2026-09-04 — both are grid-side fixes aimed at the same symptom, the cap is a no-op on production auto grids so its measurements only mean anything alongside the step numbers above, and the two were benchmarked as a pair throughout. One file, src/ForceFreeStates/Fourfit.jl.)

Summary

The Euler-Lagrange coefficient splines (fmats/kmats/gmats and the primitives) inherited every knot of the equilibrium grid. A cubic spline's third-derivative jumps at knots scale as (node error)/Δψ³, so the equilibrium's near-axis packing (Δψ ~ 1e-6 at high mpsi) amplifies even tolerance-level node error (~1e-9, measured) into huge C² kinks — and the adaptive integrator's step size becomes slaved to the knot spacing. Measured directly: core jump magnitudes grow ~34× per mpsi doubling, and an (tol/J)^¼ step model reproduces the observed step-count ladder.

The coefficients are near-cylindrical in the core and do not need that packing. This PR builds their splines on a subset of the equilibrium grid with core density capped at Δψ ≥ 0.05·ψ below ψ = 0.1. Node values are unchanged — only interpolation density — which is why the physics moves at the 1e-7–1e-8 level.

Stacked on #398 (route (a), same investigation); diff shows only the cap once #398 lands.

Measured (DIII-D stripped decks, route-(a) base, mpsi 512/1024)

	mψ=512	mψ=1024
accepted EL steps	2768 → 2188	4403 → 3034
step growth per doubling	—	1.59× → 1.39×
warm run	12.8 → 10.4 s (−19%)	—
et[1] relative change	3e-8	3e-8
Riccati Δ′ diagonal (all 5 surfaces)	4e-7	—

Not a tuned hack: the rule's basis and its generalization

Form: near the axis every component is a Frobenius power law in ψ; power laws are scale-free, so log-uniform sampling (Δψ ≥ c·ψ) resolves them at constant relative accuracy. Constant: cubic interpolation of ψ^p on a log-uniform grid errs by ~(p·c)⁴/384, so c = 0.05 resolves even the steepest spectrum component (p = m_max/2 = 11) to ~2e-4 — and the physics responds far below that because the steep components carry vanishing solution amplitude. Region: ends at min(0.1, innermost rational − RATIONAL_RES_RADIUS), and no knot inside a rational's resolution window is ever removed — guards for decks (e.g. higher n) whose rationals reach the core; verified no-ops on every current case.

Cross-equilibrium check (four equilibria, two construction paths, three grid families):

case	knots	EL steps	et[1] change
DIII-D efit, traced, log-family m1024	1025→575	−31%	3e-8
tj_analytic_direct, traced, log-family m1024	1025→575	−35%	1.3e-9
LAR, inversion, packed auto m1024	1025→575	−8%	exact to 8 digits
Solovev, traced, ldp m512/m1024	505→468	~0	~4.5e-7 abs (3e-5 rel via the ±10.4 → 0.0146 cancellation)

The LAR row is the anti-over-fit witness: clean geometry with no noise to exploit, and the cap is still harmless.

Two properties reviewers should know

1. No-op on production decks. The two-pass auto grid's core spacing already satisfies the cap — verified: 287 → 287 knots on the DIII-D auto deck, and the full harness (diiid_n1, solovev_n1, diiid_n1_riccati, gal_resistive_diiid) reproduces Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots #398's numbers exactly. The cap engages only on explicitly packed fine grids (large mpsi), which is precisely the issue Performance: Why do large equilibrium splines slow down the code? #376 scenario.
2. A fixed cap beats general curvature selection — measured, not assumed. Removal-error knot selection (drop any knot whose interpolated value matches within δ for every element of all 12 matrices) was implemented and swept: it reaches the physics plateau only at δ=1e-7 where it keeps 504/513 knots and buys nothing, and breaks et[1]/Δ′ at any δ loose enough to matter. The core knots that are safe to remove fail the interpolation test (real coefficient curvature, but the solution components multiplying them vanish as ψ^|m|/2). Interpolation error is the wrong objective; the fixed mechanism-based cap is the right scope for the ideal path. Record: handoff/issue376/RESULTS.md §22–§23 (experiment branch).

Verification

• Tests: runtests_eulerlagrange, runtests_riccati, runtests_sing — all pass.
• Harness vs develop @ 9491f89 across all four cases: identical to Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots #398's report (see property 1); Δ′-tracking cases included.
• Mechanism measurement and step-model validation: handoff/issue376/RESULTS.md §22.

Release note

• Audience: users
• Numerical impact: ideal and kinetic results move. The resample noise that grid refinement amplified is gone, so Solovev's free-boundary energy now converges in mpsi. efit_arclength gets the same fix, and F and P are now evaluated on the traced surface. (harness @ 81006c3)
• Migration: re-pin local regression baselines — ideal and kinetic values move (Solovev's free-boundary energy most of all); no interface, deck-key or output-name changes

Every flux surface is now sampled at the same straight-fieldline abscissae instead of being
splined on its own solver-chosen points and remapped, so the surface-to-surface noise that
refinement used to amplify is gone: energies converge in mpsi where they previously drifted, and
the Euler-Lagrange integration stops chasing the noise (DIII-D ODE steps roughly halve). Folded in
from #408, the Euler-Lagrange coefficient splines no longer inherit the equilibrium grid's
near-axis packing, which removes the third-derivative kinks that slaved the integrator's step size
to the knot spacing on finely packed grids.

Regression report

regress --cases diiid_n1,solovev_n1,diiid_n1_riccati,efit_fixedbdy_separatrix,solovev_kinetic_calculated,solovev_kinetic_ntv,solovev_kinetic_nuzero,solovev_kinetic_multiion --refs 57ec7e16a,81006c3a6.
Baseline develop @ 57ec7e1; branch at 81006c3 (includes #414 and today's develop).

diiid_n1                        develop         this PR         diff
total energy Re(et[1])          8.012318e-01    8.004563e-01    0.10%
NTV torque FGAR [N.m]           5.296762e-01    5.325920e-01    OK (within tolerance)
beta_n / plasma volume          unchanged       unchanged       ≤ 3e-7
ODE steps (total)               4572            1953            -57.3%

solovev_n1                      develop         this PR         diff
total energy Re(et[1])          6.773398e-01    1.073537e-02    98.42%   <- the fix; see the mpsi ladder
ODE steps (total)               618             488             -21.0%

diiid_n1_riccati                develop         this PR         diff
Δ′ diagonal, real parts         see below       see below       ≤ 2.3%
Δ′[4,4] imaginary part          -379.5i         -5.3i           spurious Im removed
total energy Re(et[1])          8.037196e-01    8.029446e-01    0.10%
ODE steps (total)               1645            1350            -17.9%

solovev_kinetic_calculated      develop         this PR         diff
total energy Re(et[1])          1.873745e+00    1.827410e+00    2.47%
total energy Im(et[1])         -1.434288e+00   -1.328439e+00    7.38%

solovev_kinetic_nuzero          develop         this PR         diff
total energy Re(et[1])          2.243336e+00    2.290424e+00    2.10%
total energy Im(et[1])         -2.026025e+00   -1.568157e+00    22.6%

solovev_kinetic_ntv             develop         this PR         diff
NTV torque fgar                 ~1.6e-4         ~3.85e-1        dominated by this PR
NTV psi quadrature evaluations  840             45              -94.6%

The Riccati row's 16.7% headline is the imaginary part of Δ′[4,4], which should be zero for an ideal
Δ′. See "Reading the 16.73%" below. Real parts move ≤2.3%, inside the ~2% grid-to-grid spread
recorded for Δ′ in #379. Two quantities that are physically zero for an ideal run sit at the
numerical floor and change sign: DIII-D Im(et[1]) and the ideal PE toroidal torque (±5e-2). The
kinetic moves are this PR's sampling plus #414's suppression, as attributed in #414.

⚠️ Requires third-party human review before merging — do not merge without an approving review.

Rebased onto develop (349a0c2)

Merged develop on 2026-09-04. The only conflict was Fourfit.jl, where develop's MatrixSplines
refactor now builds an immutable IdealMatrices: the cap logic moved into a named helper,
core_capped_knots(xs, rationals), and the constructor builds on the capped subset. Cap
behaviour is unchanged — only its plumbing follows develop.

Verification after the port: runtests_equil 279/279, runtests_grid_refinement 59/59.
Harness vs develop reproduces the headline effect — diiid_n1 ODE steps 4572 → 1974 (−56.8%),
with energies moving 0.02–0.28% (Re(et[1]) 0.10%) and q0/q95/β/singular-surface locations at 0.00%.

Resolved: DIII-D kinetic-calculated EL stepping

This PR alone takes the DIII-D kinetic_source="calculated" forward integration to 223,271 steps,
against 7,972 on develop. #414, now merged into this branch, removes the near-axis region where the
kinetic model is invalid and cures it outright: 4,246 steps, below develop's baseline. No harness case
covers this configuration yet; a DIII-D kinetic-calculated FFS-only case remains a follow-up.

Caveat, stated plainly: a three-orders-of-magnitude move in the ntv-case torque means the develop baseline for that tracked value was noise-dominated, and the new value has no independent reference yet. This needs a physics reviewer's judgment, not just the attribution argument. If this PR is accepted, the kinetic harness baselines must be re-pinned on the merge commit.

Tests

• test/runtests_sfl_sampling.jl (new) — 4 pass; mutation-checked as tabulated above
• Full suite on the merged branch — 4793 / 4793; CI green on Julia 1.11 and 1.x
• test/runtests_equil.jl — 279 pass (includes the arclength path)
• test/runtests_grid_refinement.jl — pass
• test/runtests_tj_analytic.jl — pass (exercises tj_analytic_direct, the same code path)

Scope / follow-ups

• Covers eq_type ∈ {efit, efit_arclength, imas, sol, tj_analytic_direct} — everything reaching equilibrium_solver in DirectEquilibrium.jl.
• Not affected: InverseEquilibrium.jl (chease / lar / tj_analytic / efit_by_inversion). It samples every surface on one common input θ grid through smooth 2-D interpolants, so its abscissae vary smoothly in ψ and the white-noise mechanism does not arise.
• The arclength trace stays within ~1e-6 m of the flux surface, set by its position tolerance. That is the floor of the tracer-agreement test.
• DIII-D retains a deep-core residual (steps still grow 1.59× per mpsi doubling, concentrated in ψ<0.05); under investigation separately, likely tied to the EFIT input's own resolution rather than the trace.

⚠️ Requires third-party human review before merging — do not merge without an approving review.

🤖 Generated with Claude Code

---

Follow-up verification: Δ′ and wall time (requested in review)

The original evidence used diiid_n1/solovev_n1, which run integrator = "forward" and therefore
do not emit the BVP Δ′ matrix at all. Re-ran the Δ′-tracking cases.

diiid_n1_riccati

quantity	develop	local	diff
delta prime (BVP diagonal)	[5 elem]	[5 elem]	16.73%
delta prime (raw side-major)	[10 elem]	[10 elem]	2.35%
edge coil response delta_coil	[10 elem]	[10 elem]	1.26%
total energy Re(et[1])	8.037196e-01	8.029430e-01	0.10%
# singular surfaces / psi / q	5	5	OK / 2e-09 / 0.0
ODE steps (total)	1649	1358	−17.7%
Runtime	202.5s	192.8s	−4.8%

gal_resistive_diiid — the quantities the tearing/matching path consumes

quantity	diff
gal PEST3 Δ diagonal	0.01%
‖gal Δ′ matrix‖	0.03%
gal D_I per surface	0.00%
gal α per surface	0.00%
‖gal inner-layer Δ‖	0.00%
‖gal Δ_coil block‖	0.80%
‖gal match cout‖	1.38%
gal match residual	OK
Runtime	216.8s → 213.6s

Reading the 16.73%

It is one element. Δ′ BVP diagonal across an mpsi ladder on the riccati deck:

mpsi	develop	this PR
256	[9.7, −1.9, −13.1, 103024, 336.0]	[9.7, −1.9, −13.1, 102768, 336.1]
512	[7.2, −5.2, −16.6, −1319.5, 61.4]	[7.2, −5.2, −16.6, −1293.2, 61.5]
1024	[8.6, −5.5, −16.1, −2368.8, 55.1]	[8.6, −5.5, −16.1, −2368.0, 55.6]
drift vs previous mpsi	174% → 79.5%	174.7% → 83.1%

The 4th entry (q = 5 surface) runs 1e5 → −1319 → −2369 and changes sign: it is not converged in
mpsi on either version, so it cannot discriminate between them. The other four diagonal entries
agree between versions at every grid, and version-to-version agreement on the full diagonal is
0.25% / 2.0% / 0.03% at mpsi 256 / 512 / 1024.

Caveat worth stating plainly: this PR does not improve Δ′ convergence either (drift
174.7%/83.1% versus develop's 174.0%/79.5%). The unconverged q = 5 Δ′ element is a pre-existing
issue that this change neither causes nor fixes, and it deserves separate attention.

Wall time improves on every case measured: forward diiid_n1 210.6 → 184.5 s, diiid_n1_riccati
202.5 → 192.8 s, gal_resistive_diiid 216.8 → 213.6 s.

[logan-nc]

Report with visuals for workflow and results can be found here: https://claude.ai/code/artifact/873c8376-3f75-4fe4-934b-4d39f23e2376

[logan-nc]

Review package (visual companion to this PR — mechanism diagram, the Solovev et[1] grid-convergence before/after, step ladders, geometry-residual evidence, and the Δ′ referee table):

📦 https://claude.ai/code/artifact/873c8376-3f75-4fe4-934b-4d39f23e2376

Self-contained page; complements rather than duplicates the description and diff. If the link does not resolve, ask @logan-nc to enable sharing on it.

[logan-nc]

Review package (self-contained, private artifact): https://claude.ai/code/artifact/183e8712-e9fd-48c7-9ca0-1500f51ac263

Covers both folded changes — the aligned SFL surface sampling and the coefficient-spline core knot-density cap (#408, folded in here) — with the mechanism diagrams, the grid-convergence ladder, the C³-jump/knot-density figures and the cross-equilibrium generalization check. Numbers on the page were measured against develop @ 9491f893d; the post-port harness run reproduces them and the body below carries the current table.

[logan-nc]

Does #480 make this stack redundant? Measured: no

#480 (per-column ODE absolute tolerance) also cuts ODE steps, by fixing the cause of the step growth reported in #379. I checked whether that makes the surface sampling here, or the #408 knot cap folded into this PR, unnecessary.

Surface sampling: still needed. It fixes error in the equilibrium data the ODE integrates. #480 changes how tightly the ODE integrates it, and cannot remove noise already in the coefficients. Harness results on the standard decks:

	develop	#480	this stack	stack + #480
Solovev et[1]	0.6773	0.6773	0.0107	0.0107
DIII-D et[1]	0.80123	0.80123	0.80046	0.80046
DIII-D steps	4572	1424	1974	799
Solovev steps	618	410	495	283
DIII-D Riccati steps	1645	4166	1358	2622

#480 leaves every energy unchanged to about 1e-6. The grid-dependent Solovev answer this PR corrects is therefore still present with #480 alone. The step savings compound.

Knot cap (#408): still needed on packed grids. On the auto grid it is inert: every harness case is bit-identical with it on or off, as stated above. On explicitly packed DIII-D grids it still saves steps on top of #480:

packed grid, with #480	cap off	cap on
mpsi 512	879	794 (−10%)
mpsi 1024	1293	979 (−24%)
Riccati, mpsi 1024	3859	2976 (−23%)

et[1] moves by about 3e-8 and Δ′ by about 4e-7. The two fixes are complementary: #480 addresses loose tolerance, and the cap addresses spline kinks at densely packed knots.

Not yet checked: the Solovev grid-convergence ladder in this PR predates #480. I have not run develop + #480 alone at mpsi 256/512/1024. Since #480 leaves the standard-grid energy unchanged, it very likely still drifts.

_{The "stack" columns are #414's head (2abea51d2), which carries this branch. #414's own changes affect kinetic runs only, so the ideal rows reflect this PR. Baselines: develop 56eafc2a5 and #480 d73206dc7. The cap was disabled by a one-line ablation on local throwaway branches, none of them pushed.}

🤖 Generated with Claude Code

[logan-nc]

@d-burg, I've added a section to the description, "Why Solovev's energy moves 98%: an independent check", that may be relevant to your cylindrical tearing work. On develop, directly traced equilibria (sol, tj_analytic_direct, efit, and the Cerfon–Freidberg equilibrium in #495) carried an etol-dependent geometry error. On the Solovev deck it is about 1e-3, and it feeds straight into the ψ-derivatives that Δ′ depends on. The section also explains how this interacts with your #491, which changes the tolerance on the same tracer solves.

🤖 Generated with Claude Code

[logan-nc]

@d-burg I'm gonna merge this weds latest. A read through before then would be 👍 🙏

[d-burg]

Some final claude notes and proposed to-dos:

Please fix before merge

1. Release note. It says no deck-key or output changes, but #414 came in with this branch: axis_validity_suppression (default on), KineticForces/Validity, and a forced reg_spot = 0 in kinetic runs. The title should mention it too.
2. Solovev mpsi ladder. It shows et[1] = 0.01462 "converged to 6 figures", but the harness deck (mpsi 128) gives 0.010735. Please re-run the ladder at this head, including 128.
3. Two claims in the description.
   • The cap is not a no-op on auto grids: it removed 2–43 knots on five of my equilibria (harmlessly).
   • "Spurious Im removed" is not general: on fixed grids Im Δ′ dropped on one surface and rose on others.
4. Unit test for core_capped_knots — see the inline comment.
5. Harness at head for diiid_slayer_n1, gal_resistive_diiid and gal_resistive_pe; the current report covers 8 of 16 cases.

Follow-up issue is fine

• After this PR the Solovev fixture is a ~1000:1 cancellation (et[1] = 0.0107 against ep = −10.42), so the Solovev cases will amplify every later change. The wall distance needs re-tuning, and the deck comment still says et[1] = +0.24.
• I suspect the 2500× NTV torque move is resonant amplification from sitting that close to marginal: (0.677 / 0.0107)² ≈ 4000. Comparing the response amplitude between the two runs would settle it.
• #414's physics still wants a look from someone with PENTRC experience, as @ebursch asked.

[d-burg]

From claude: nothing is resampled any more; the fit is only needed for the θ-derivative.

Suggested change

	# Evaluate the θ-fit and its derivative at the nodes (the Jacobian needs the derivative)

[d-burg]

From claude: the refinement callback saves pre-refinement states (save_positions=(true, false)), so the dense interpolant is offset by one step's refinement and these samples sit slightly off ψ = ψ₀ (tolerance-level). Not blocking. Calling direct_refine(u[2], eta, psi0, params) on each sample would put them on the surface exactly. Worth a follow-up?

[d-burg]

The 2e-4 figure is for the n = 1 DIII-D spectrum. The cap constant does not scale with mpert.

Suggested change

	steepest component (p = m_max/2) to ~2e-4, well below where the physics responds. The capped region
	steepest component (p = m_max/2) to ~2e-4 at m_max = 22 (~1e-2 at m_max = 60). The capped region

[d-burg]

From claude: This has no test, and no harness case removes a knot, so CI never exercises it. Please add this as test/runtests_core_capped_knots.jl (passes 8/8 at this head) and include it in runtests.jl:

using Test
using GeneralizedPerturbedEquilibrium

@testset "core_capped_knots" begin
    FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
    R = GeneralizedPerturbedEquilibrium.Equilibrium.RATIONAL_RES_RADIUS

    # Densely packed grid, one rational far from the core.
    xs = collect(range(1e-3, 0.99; length=2000))
    keep = FFS.core_capped_knots(xs, [0.5])
    kept = xs[keep]
    @test keep[1] == 1 && keep[end] == length(xs)
    @test issorted(keep) && allunique(keep)
    @test length(keep) < length(xs)
    # Nothing at or above ψ = 0.1 is dropped.
    @test count(>=(0.1), kept) == count(>=(0.1), xs)
    # Below it, kept knots are at least 5% of ψ apart.
    core = kept[kept .< 0.1]
    @test all(diff(core) .>= 0.05 .* core[2:end])

    # A core rational ends the capped region before its resolution window.
    keep_rat = FFS.core_capped_knots(xs, [0.05])
    @test all(i in keep_rat for i in eachindex(xs) if abs(xs[i] - 0.05) <= R)

    # A grid already sparser than the cap is returned whole.
    sparse = collect(range(0.01, 0.99; length=20))
    @test FFS.core_capped_knots(sparse, [0.5]) == collect(1:20)
    @test FFS.core_capped_knots([0.1, 0.2], Float64[]) == [1, 2]
end

[d-burg]

Do we want fortran source citations like this one?

Suggested change

	its validity domain. Returns 0.0 when no criterion is met anywhere. The Fortran precedent
	(`ktanh_flag`, dcon/fourfit.F) suppressed the same region with four hand-tuned knobs; here the
	boundary is derived from the profiles with no user parameters.
	its validity domain. Returns 0.0 when no criterion is met anywhere. The boundary is derived from
	the profiles with no user parameters.
	``

[d-burg]

This is a new key that is on by default and changes kinetic results. I think it needs to be in the release note, along with the forced reg_spot = 0 and the new KineticForces/Validity group

[d-burg]

From claude: neither test covers ν: it is zero here, and the EFIT test checks only the Jacobian. Could the EFIT test also compare rzphi_nu? Solovev is analytic, so a direct comparison of r² against the exact surfaces would also be a stronger check than mpsi 64 vs 128.