ForceFreeStates - BUGFIX! - Use a per-column absolute tolerance in the Riccati shooting integration

[d-burg]

Release note

• Audience: users
• Numerical impact: Riccati Δ' matrices at the default singfac_min move at the 0.00–0.01 % level (float-level tightening); at smaller singfac_min the previous results were biased low (3/1 Δ' on the DIII-D-like case: 0.909 at 1e-4 rising to a 1.055 plateau by 1e-6 after the fix, 4-way consistent with Galerkin, Fortran and TJ). SLAYER growth rates at defaults change by 0.02 %, which is at the per-surface reproducibility floor of the threaded AMR root search (~0.1 Hz discrete steps, see Regression - BUGFIX - Fix two silent failure modes: PATH julia resolution and single-sample result caching #418); attribution of that residual to the Δ' shift is plausible but not resolved by the harness measurement. (harness @ 362662b)
• Migration: no input changes. Riccati Δ' results computed with a non-default singfac_min (below 1e-4) should be regenerated; values at the default move within 0.01 %.

The Riccati shooting batches were solved with a relative tolerance only, so DifferentialEquations' uniform default absolute tolerance (1e-6) let the large column of each batch set the error floor of the small column, biasing Δ' low by an amount that grew as the singular-surface exclusion radius shrank. Each column now carries its own absolute tolerance scaled to its initial magnitude, making Δ' independent of singfac_min to 0.07 % per decade. This is a tightening of the integration tolerance, not a relaxation.

Regression report

regress --cases diiid_n1_riccati,diiid_slayer_n1 --refs c42d558e,362662b9 (branch head 362662b is the develop-merge carrying #339/#399; c42d558 is its develop side, so the diff isolates the Riccati.jl fix. Same pinned manifest, julia 1.11.6):

Regression Report: diiid_n1_riccati
Ref 1: c42d558e  @ c42d558ef (2026-08-20)
Ref 2: 362662b9  @ 362662b9e (2026-08-20)
Quantity                       c42d558e       362662b9       Diff               Status       
delta prime (BVP diagonal)     [5 elem]       [5 elem]       1.644e-02 (0.00%)  ** CHANGED **
delta prime (raw side-major)   [10 elem]      [10 elem]      8.087e-01 (0.01%)  ** CHANGED **
edge coil response delta_coil  [10 elem]      [10 elem]      1.210e-03 (0.01%)  ** CHANGED **
total energy Re(et[1])         8.038421e-01   8.038421e-01   0.0e+00            OK           
plasma energy Re(ep[1])        -1.344588e+00  -1.344588e+00  0.0e+00            OK           
vacuum energy Re(ev[1])        2.148430e+00   2.148430e+00   0.0e+00            OK           
total energy (all)             [35 elem]      [35 elem]      0.0e+00            OK           
# singular surfaces            5              5              0.0e+00            OK           
singular psi locations         [5 elem]       [5 elem]       0.0e+00            OK           
singular q values              [5 elem]       [5 elem]       0.0e+00            OK           
ODE steps (saved)              51             51             0.0e+00            OK           
ODE steps (total)              1671           1671           0.0e+00            OK           
mpert                          35             35             0.0e+00            OK           
npert                          1              1              0.0e+00            OK           
q0                             1.204212e+00   1.204212e+00   0.0e+00            OK           
q95                            4.781723e+00   4.781723e+00   0.0e+00            OK           
beta_n                         1.372511e+00   1.372511e+00   0.0e+00            OK           
Runtime (s)                    109.2s         111.4s                            --           
Summary: 3 changed, 14 unchanged

Regression Report: diiid_slayer_n1
Ref 1: c42d558e  @ c42d558ef (2026-08-20)
Ref 2: 362662b9  @ 362662b9e (2026-08-20)
Quantity                            c42d558e  362662b9  Diff               Status       
SLAYER surface indices              [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER poloidal m                   [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER toroidal n                   [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER minor radius rs              [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER r-based shear                [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER Lundquist S                  [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER D_norm                       [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER P_perp                       [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER tauk                         [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER iota_e                       [6 elem]  [6 elem]  0.0e+00            OK           
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]  [3 elem]  5.4e-06            OK           
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]  [3 elem]  4.3e-03            OK           
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]  [3 elem]  1.623e-01 (0.02%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]  [3 elem]  0.0e+00            OK           
SLAYER enabled flag                 1         1         0.0e+00            OK           
Runtime (s)                         130.0s    131.4s                       --           
Summary: 1 changed, 14 unchanged

Energies, singular-surface locations and ODE step counts are bit-identical: the fix changes only the per-column absolute tolerance of the shooting batches (a tightening; no tolerance was loosened). The only flagged SLAYER quantity is γ_Hz at 0.02 %, which sits at the root search's own per-surface jitter floor (~0.1 Hz) and should not be read as a resolved mechanism; every layer input is bit-identical and the ψ_N-referenced Δ' moves 0.00–0.01 % on deterministic BVP quantities.

Fixes #412.

Root cause

integrate_fm_with_ua_ic solved each shooting batch with reltol only, leaving the ODE library's uniform default abstol = 1e-6 in force across columns spanning ~10⁻³ (resonant small solution) to ~10³ (resonant big solution). The batch's largest column therefore set the absolute error floor of its smallest, and the extracted Δ' inherited an absolute error proportional to the big solution's magnitude at the matching radius — decaying only as singfac_min shrinks (measured ≈ dpsi^1.2, close to the Frobenius exponent gap 2α = 1.0), and invisible to grid-resolution convergence scans. The reference implementation avoids this by construction: it rebuilds a per-column atol from each column's own magnitude every step, so its Δ' shows only −1.4 % sensitivity over the same singfac decade where this code moved −13 %.

Fix

11 lines: build a per-column abstol array scaled to each column's own magnitude (max|column| × reltol) for both shooting batches, mirroring the reference implementation's per-column policy. This tightens error control (replacing an implicit uniform 1e-6 floor); no criterion is relaxed. A structural alternative (confining the ua IC to the surface-adjacent chunk, composing the rest from identity-IC propagators) was implemented and tested first: it only half-fixed the benchmark and regressed the small-radius plateau, so it was reverted — error control, not integration structure, was the load-bearing difference.

Evidence (TJ circular benchmark, ε-scan point, ldp mpsi=256, n=1)

	Δ'_31 @ singfac 1e-4 (default)	Δ'_31 @ 1e-6	Δ'_21
develop	0.9087 (−14 % vs consensus)	1.0552	10.868
this PR	1.0544	1.0537	10.875
consensus	julia Galerkin 1.0567 · Fortran STRIDE 1.05 ± 0.02 · TJ (ψ_N) 1.067

Δ' is now matching-radius independent (0.07 %/decade, matching the reference implementation's behavior). The full diagnostic chain (Galerkin cross-check, Fortran/TJ corroboration, per-surface singfac scans, sing_order irrelevance proof, line-by-line series/BVP audits) is in #412.

Validation

• Shipped shaped DIII-D example at default settings: all six Δ' diagonals within 0.05 % of develop (no behavior change where the old error was negligible).
• Test files pass: runtests_sing, runtests_riccati, runtests_parallel_integration (94/94), runtests_slayer_riccati (via earlier suite run).
• Regression harness vs merge base 9491f89 (pre-InnerLayer.SLAYER - BUGFIX! - Use the physical toroidal field, not the b0exp normalization #399, matching the branch base at test time): diiid_slayer_n1 15/15 unchanged; diiid_n1_riccati Δ' diagonal/raw/coil move 0.00–0.01 % (float-level tightening shifts); gal_resistive_diiid 0.01 %. Branch since rebased onto develop (InnerLayer.SLAYER - BUGFIX! - Use the physical toroidal field, not the b0exp normalization #399); the fix touches only Riccati.jl, which InnerLayer.SLAYER - BUGFIX! - Use the physical toroidal field, not the b0exp normalization #399 does not.

What moves downstream

Small Δ' matrix elements (|Δ'| ≲ 1–2) computed at the default singfac_min were biased low by up to ~0.15 absolute; they now agree with the Galerkin path and external codes. Large elements move ≤ 0.1 %. Cases whose tearing analysis hinges on small/near-marginal Δ' should be re-examined after merge.

Important

This changes the effective error-control semantics of the Δ' shooting integration (tightening an implicit library default). Do not merge without human review.

🤖 Generated with Claude Code

[d-burg]

Scan-scale validation: full TJ β + ε scans on this branch, both grids — 196/196 runs, zero failures

A/B against the frozen develop scans (same equilibria, same settings; julia_*_413_{ldp256,auto}.csv vs the develop CSVs):

quantity	β ldp@256 (med / max)	β auto (med / max)	ε ldp@256 (med / max)	ε auto (med / max)
Δ'_21 change	0.000 % / 0.01 %	0.000 % / 0.01 %	0.012 % / 0.11 %	0.012 % / 0.11 %
Δ'_31 change (pole-free window)	0.006 % / 2.4 %	0.006 % / 3.5 %	3.1 % / 16.0 %	3.2 % / 15.4 %

Exactly the bug's fingerprint: the 2/1 (|Δ'| ~ 7–22) is untouched at the 0.1 % level across all 196 points, and the movement is confined to the small ε-scan 3/1 elements (|Δ'| ~ 0.9–1.5) where the absolute error floor lived.

Closure against the independent reference: re-cutting the K^(2μ)-corrected comparison with this branch's data, the ε-scan 3/1 median |difference from TJ| drops 5.7 % → 2.3 % (spread [−14.6 %, +0.4 %] → [−6.3 %, +5.4 %]) — the previously "unexplained low-ε 3/1 drift" is now confirmed to have been this bug. The 2/1 panels are unchanged (median 1.4 %).

Scan runner + data: julia_deltaprime_harness/outputs/develop_benchmark/run_fix413_scans.jl (maintainer's benchmark repo). One incidental find while setting this up: a TOML with psi_accuracy > 0 silently engages the two-pass auto grid builder even when grid_type="ldp" is requested — worth a guard or warning in a follow-up.

🤖 Generated with Claude Code

[d-burg]

Forced-rerun verification of the regression stamp — CORRECTED. An earlier version of this comment claimed the gamma_Hz residual of 1.623e-01 Hz (0.02%) "reproduces exactly, four orders of magnitude above the measured sample-noise ceiling — it is mechanism." That claim was wrong and is withdrawn.

The ceiling it relied on (1.5e-05 Hz) came from two forced samples that happened to land in the same solver state. The actual same-source reproducibility floor, measured by comparing a single commit against itself across forced re-runs (layer-overlap investigation, tracked in #418), is a discrete per-surface step of ~0.08–0.15 Hz in gamma_Hz — the threaded AMR root search lands in one of two nearby states per surface. The 0.16 Hz residual on this PR is the same magnitude as that step, so the measurement cannot distinguish the Delta-prime mechanism from the floor; "reproduces to four digits across two pairs" is consistent with landing in the same state four times.

What still stands: the structural single-channel argument (the entire source diff is src/ForceFreeStates/Riccati.jl, every layer input is bit-identical, and the uncoupled case's only path from that file to the dispersion root is the Delta-prime matrix, which moves 0.00–0.01%) makes a Delta-prime-driven shift of this size plausible, and the 6x6 matrix path is confirmed taken (no stub fallback). But attribution of the specific 0.02% number is not resolved by this measurement — it is at the jitter floor. The headline of this PR is unaffected: the Riccati Delta-prime correction itself (0.909 → 1.055 plateau at small singfac_min) is deterministic ODE output, not a root-search quantity, and the default-deck Delta-prime movement (0.00–0.01%) is measured on bit-reproducible BVP quantities.

The cached-vs-forced ω residual (4.3e-03 → 1.8e-03 Hz) is likewise at or below floor. All layer inputs remain bit-identical. The stamped table stands as a record of what the harness reported; its gamma line should be read as "at the per-surface jitter floor."

Follow-ups: the root-search nondeterminism itself is tracked in #420 (seeded/serial candidate ordering + an in-process run-twice reproducibility test); until it lands, the diiid_slayer_n1 gamma_Hz noise threshold has been raised to 2.5e-1 as a documented stopgap derived from the measured per-surface steps (commit b3d65ec on #298).

[github-actions]

This pull request is missing 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.