[matt-pharr]

Release note

• Audience: users
• Numerical impact: none in the harness cases; q profiles with exactly quadratic spline cells now keep their extrema and rational surfaces (harness @ 131c83d)
• Migration: scripts that pass a mode list that is not contiguous and ascending to the 2D vacuum response now get an ArgumentError; pass mlow:mhigh. Deck runs are unaffected.

Fixes six numerical edge cases found in a review of v2.0.0. An exactly quadratic q profile lost its interior extremum, and with it the rational surfaces on either side. The rational-surface search could loop forever on a flat rational segment. The NTV-limited correction current could come out as zero for a weak residual torque. Rank-deficient joint correction solves threw an error or returned meaningless currents. allowable_tolerance returned NaN for scans that start at scale 0. The 2D vacuum response silently computed the wrong poloidal modes for a mode list that is not contiguous.

Changes

• Equilibrium (equilibrium_qfind!): each cell's q′ roots now come from the cancellation-safe quadratic formula, which also covers the linear cell (q‴ = 0) that was skipped before. Cubic splines fit a quadratic exactly on a uniform grid, so most of its cells have q‴ = 0 and the minimum was lost.
• ForceFreeStates (_find_rational_surfaces): an interval whose two extrema have the same q is skipped. If that q was rational, dm = 0 and the loop never ended.
• ErrorFields (correction_current, correction_current_upper): the closed-form roots are rewritten in an algebraically identical form that does not cancel when 4a·excess ≪ C².
• ErrorFields (correction_requirement): the joint least-squares solve uses the pseudoinverse. Identical or collinear arrays now give the minimum-norm currents; before, a square system threw SingularException and a near-singular one returned currents around 1e16 kAt.
• ErrorFields (allowable_tolerance): a target equal to a sampled risk returns that sample's scale exactly, and an interval that starts at scale 0, which has no logarithm, is interpolated linearly.
• Vacuum (2D response): a mode list that is not contiguous and ascending is rejected. Before, [1, 3] computed modes 1 and 2. Deck runs always pass contiguous ranges; only the scripting API could reach this.

Each fix has a test, and each new test fails on the unfixed code. The exception is the flat-segment test, which would hang there.

Regression report

Run against the branch's base 5503af4fa rather than develop, so that #626, which merged after this branch was cut, does not appear in the diff. The branch merges into develop without conflicts.

Run:  regress-fanout numfix_131c83d1b (SDCC sirius), julia 1.12.6, pinned manifest
Ref 1: 5503af4fa
Ref 2: bugfix/numerics-review-fixes @ 131c83d1b

diiid_error_field               10 unchanged
diiid_kinetic_multi_n            8 unchanged
diiid_multi_n                   19 unchanged
diiid_n1                        53 unchanged
diiid_n1_riccati                17 unchanged
diiid_n1_riccati_precompiled    17 unchanged
diiid_slayer_n1                 22 unchanged
diiid_slayer_n1_chi_perp_e      22 unchanged
diiid_slayer_n1_coupled          7 unchanged
diiid_slayer_n1_toroidal         5 unchanged
diiid_slayer_n1_winding         22 unchanged
diiid_truncation_n1             12 unchanged
efit_fixedbdy_separatrix         5 unchanged
gal_resistive_diiid             11 unchanged
gal_resistive_diiid_galinner    11 unchanged
gal_resistive_pe                 8 missing/N/A (N/A on both refs, as on develop)
ggj_ray_q500i                    4 unchanged
ggj_reference                    4 unchanged
lar_ggj_tearing_n1               3 unchanged
lar_reversed_shear_n1           22 unchanged
solovev_kinetic_nuzero          14 unchanged
solovev_multi_n                 15 unchanged
solovev_n1                      22 unchanged
tj_kinetic_calculated           20 unchanged
tj_kinetic_multiion              6 unchanged
tj_kinetic_ntv                   7 unchanged
tj_n1                           25 unchanged

Notes for reviewers

• The equilibria in the harness have no exactly quadratic spline cells, so nothing moves. The q-extrema fix only matters for analytic or synthetic q profiles on uniform grids.
• The same review found further issues that are left for separate PRs:
  • RLAR torque is a linear-amplitude placeholder, and CLAR returns zero.
  • The regression harness counts a dataset that disappears from one ref as N/A instead of a change, so --fail-on-change misses it.
  • A rerun still reads the kinetic and wall files from disk instead of from the snapshot.
  • The rerun overwrite guard compares paths as strings, so an alias of the source file gets through.
  • The 3D vacuum quadrature cache has no lock.

🤖 Generated with Claude Code

https://claude.ai/code/session_014wKUyAaxKdTJ5ch1boMFE1

[logan-nc]

Review at 131c83d1b. I checked items marked ran by running them, and the rest by reading the code (traced).

What held up

• q extrema (traced and ran). In each cell q′ = b + c·x + (d/2)·x². The roots 2s/d and b/s are the cancellation-safe pair for that quadratic.
  • When q‴ = 0, 2s/d is ±∞ and gets rejected, while b/s = −b/c is the linear root.
  • s never reaches 0 on the two-root branch.
  • The change also fixes q‴ that is tiny but not zero, where −c/d ± √… used to cancel.
  • Ran: a minimum exactly on a knot is found exactly once on 15 combinations of grid size and knot position. It is neither lost nor duplicated.
• Flat segment (traced). The loop can only fail to end when dq = 0 exactly, so the guard is sufficient.
• NTV correction roots (traced). Both rewrites give the same roots in a form that avoids cancellation, given delta_per_kat ≥ 0 (it is |δ| per kAt). With zero residual torque, the upper edge used to be NaN; it is now the linear limit.
• pinv (traced). pinv and rank(R) share the default tolerance, so the reported rank matches the solution. Only collinearity at round-off level is cut: nearly collinear arrays still give large, opposing least-squares currents. That is correct least squares and fine for v2.0.
• 2D vacuum (traced). The 3D path builds its basis from m_modes directly, so only the 2D path needs the guard.
• Tests. Each test is targeted, checks against an independent answer, and is fast.
• Repo rules. The stamp is at the head commit, the release note is present and the milestone is set. No PR numbers appear in the source.

Pushed: 517613152

• allowable_tolerance's docstring still said the interpolation is "in the logarithms of both the scale and the risk". It now notes the linear case for a bracket that starts at scale 0.
• The scale-0 test checked 0 < x < 1. The exact answer is 0.5, so it now uses ≈ 0.5.

These are a docstring and a test only, so the stamp stands. Rerunning runtests_risk, _error_fields, _ntv_limits, _sing and _vacuum on that commit: pending, I'll update here.

This is not an approval. A third-party human reviewer's approval is still required before merge.