PerturbedEquilibrium - FEATURE - Solve the perturbed equilibrium for several toroidal modes in one run

[jhalpern30]

Release note

• Audience: users
• Numerical impact: single-n cases unchanged; new diiid_multi_n case (harness @ 5e65f13)
• Migration: none. New outputs PerturbedEquilibrium/mode_m and mode_n; KineticForces now errors on multi-n runs.

PerturbedEquilibrium now runs with nn_low < nn_high: response matrices, energies, torque, field profiles and singular coupling cover every n in one run, and each n block reproduces the matching single-n run. On the DIII-D example, n > 1 is currently limited by a separate ForceFreeStates fix (see Findings).

Review summary (3 min read, figures for both Solovev and DIII-D): https://claude.ai/artifact/5X3EewrkkJmo4QgQdSXnqK

Regression report

PR head 5e65f13 against develop e88b9b0:

diiid_n1         53 unchanged
solovev_n1       22 unchanged
solovev_multi_n  15 unchanged
diiid_multi_n    new case (example absent on develop); baseline on this branch:
                 PE plasma energy 2.937414e+00, vacuum 3.114145e+00, surface 3.539610e+00,
                 torque 4.817613e-03, total energy Re(et[1]) -5.752607e-01,
                 7 resonant (surface, n) rows, mpert 29 × npert 2

Verification against single-n

Each case was run as n = 1, 2 together and as n = 1 and n = 2 alone, with the m range matched per n.

Solovev (Solovev_ideal_example, q₀ = 2.2, mixed-n ASCII forcing; W healthy, anti-Hermitian at 1e-7):

• n = 1: Λ, P, ϱ and the response field agree to 2e-9 to 2e-8.
• n = 2: they agree to 4e-5 to 9e-5, which is set by where the ODE crosses the q = 3 surface shared with n = 1, not by PerturbedEquilibrium.
• Both n: L and the forcing field are exact, resonant field and island width agree to ≤ 1.4e-3 per row, and vacuum, surface and plasma energies equal the sum of the single-n runs to ≤ 3e-6. The torque is zero to roundoff in both, and the cross-n blocks are exactly zero.
• Each n = 2 singular row is compared with the single-n run that crosses its surface at the same distance. A surface shared by several n is crossed at singfac_min/min(n).

DIII-D (DIIID-like_ideal_multi_n_example):

• n = 1 block: Λ, P, ϱ and the response field agree to 1e-9 to 3e-8. Resonant field and island width agree to about 5e-6. L, the forcing field and vacuum energy agree exactly. The cross-n blocks of Λ, L, P and ϱ are exactly zero.
• n = 2 block: not verifiable on current develop. The ForceFreeStates W is 33–65% anti-Hermitian for n = 2 in every run, single-n included, so the comparison is dominated by that upstream error. Once W is healthy, multi-n and single-n agree to 4e-5 in Λ and 1e-4 in the response field.

Findings (upstream, not fixed in this PR)

• The forward Euler-Lagrange integrator has no per-column absolute tolerance. EulerLagrange.jl passes only reltol, so abstol falls back to 1e-6. With the fixed axis start (U₁ = 0), near-axis entries of U₁ stay far below that and are integrated with no error control. The symplectic invariant U₁†U₂ − U₂†U₁ drifts from ψ ≈ 0.01. Fortran ode_step sets atol = max|u(:,isol,ieq)|·tol before every step, and the Julia Riccati propagator already does the same.
• Tested fix: a per-column abstol together with ucrit = 1e3 (the value in Fortran's shipped dcon.in). It takes W from 0.2–1.0 anti-Hermitian to ≤ 2e-6 in every DIII-D case tried: single n = 2 and multi-n n = 1–2 and 1–3, each across nudges of singfac_min and psilow. It also roughly halves the ODE steps. Solovev cases are unchanged.
• Issue Large non-hermitian component of plasma response matrix when calculating crit #93 case (DIIID-like_ideal_example, fixed start): the least-stable energy goes from −5426 − 3126i to 0.8014 + 8e-5i, and the non-Hermitian flags from 2243 steps to 7. The remaining ~1e-4 is introduced at the first singular-surface crossing.
• Column selection at crossings: changing which column is zeroed at a crossing proved unnecessary once the tolerance is fixed (bit-identical results in all 15 test runs). The Fortran index(1) rule stays.

Follow-ups: a ForceFreeStates BUGFIX PR for the tolerance and ucrit changes, retitling #93 to the crossing residual, an issue for the Frobenius m = 0 start, and a re-baseline of diiid_multi_n once the fix lands.

Notes for reviewers

• Mode axes span numpert_total = mpert·npert (m fastest, one block per n), labelled by PerturbedEquilibrium/mode_m and mode_n. Poloidal convolutions and DFTs act within each n block, since the axisymmetric equilibrium does not couple different n.
• The surface inductance L is built block-diagonal in n from one Vacuum call over nlow:nhigh.
• Plasma energy and torque use only the diagonal n blocks of Λ: T = Σₙ −2n·Im(Φₙ†·Λₙₙ⁻¹·Φₙ)/4 (Park 2011 PoP 18 110702, eq. 1).
• Singular-coupling rows are sorted by ψ across all n. An integer-q surface appears once per resonant n, and those rows are not counted as each other's Chirikov neighbours.
• Vacuum: the Legendre quadrature cache is looked up once per kernel call and passed down, fixing a thread race when several n run concurrently.
• set_psilim_via_dmlim is ignored for multi-n; the new example sets the edge with qhigh.
• Equilibrium: warns when grid_type = "auto" or "log_asymptotic" is used with psihigh ≤ 0.98, which gives a non-monotonic ψ grid.

[ebursch]

Should there be a way for the user to change the core/mid/edge breakdown or a physics motivation for these psi values?

[ebursch]

I guess that's unrelated to this PR though

[jhalpern30]

Yeah I really only added that here because it was silently breaking some of my runs, so at least this way it dumps an output. I am not familiar with the auto/log_asymptotic parts of the code, so I'll just flag @logan-nc on this. I think its relatively moot since most humans would not set psihigh < 0.98 anyway, so this acts more as a check on an LLM setting an unphysically low psihigh.

[logan-nc]

Agreed. This is fine for now as there is no known physical reason psihigh should ever be less than 0.98. If one comes up, we can improve this to be more robust.

NO, the user should not be exposed to this. The biggest learning hurdle of the Fortran code was the ballooned name lists. We are working on the premise that we should strive to automate things robustly instead of turning everything into a manual knob.

[ebursch]

This looks ready to go.