EQUIL - IMPROVEMENT - Sample every flux surface at the same straight-fieldline angles

[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.

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, baseline develop @ 9491f89
(branch is current with origin/develop; 1 ahead, 0 behind).

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 ~1e-6-level white-in-ψ geometry error (measured; see the mechanism section), 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.

⚠️ Known regression this PR introduces: DIII-D kinetic-calculated EL stepping

Bisected cleanly (same deck, worktrees per commit): the DIII-D kinetic_source="calculated" forward integration takes 223,271 steps on this branch vs 7,972 on develop (FFS wall 1079 s vs 427 s), with physics agreeing to ~0.1% (et[1] 1.005545 vs 1.006441). This PR alone reproduces the full-stack number bit-for-bit; the stacked follow-ups add nothing. The proximate signature is an integrator grind at ψ ≈ 0.0012–0.0015 (median step 9e-9 vs develop's 2e-6); the deep mechanism is unresolved at the h5 level — the kinetic matrices, their spline third-derivative jumps, condition numbers, crit, integration start point, and solution norms are all identical at every stored digit between the two runs (open follow-up: the unstored F-matrix spline / step-control internals). No harness case covers this configuration — a DIII-D kinetic-calculated FFS-only case should be added. The stacked PR #414 (near-axis validity suppression) removes the sensitive region on physics grounds and cures the regression outright (4,246 steps — below develop's own baseline).

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_equil.jl — 279 pass (includes the arclength path; its return-type annotation was updated to match)
• 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 covered: InverseEquilibrium.jl (chease / lar / tj_analytic / efit_by_inversion), which does its own two-stage SFL resample. Follow-up.
• 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.

[github-actions]

This pull request is missing an assignee and a reviewer.

If you are not ready to name them, mark this pull request as a draft.
docs/development/contributors.md suggests lead developers to ask.
Merging is not blocked here, but no pull request may be merged without human review.