ForceFreeStates - BUGFIX! - Give every Euler-Lagrange solve an explicit absolute tolerance and default to Vern7

[d-burg]

Release note

• Audience: users
• Numerical impact: ForceFreeStates (FFS) results move toward the converged solution. Every Euler-Lagrange solve had an unset absolute tolerance (the OrdinaryDiffEq default 1e-6), which capped accuracy regardless of eulerlagrange_tolerance. Against a reltol 1e-12 / abstol 1e-14 reference, NTV torques on the Solovev kinetic decks were 1.5–1.7 % off on develop and are 0.30–0.36 % off here. Δ′ on the DIII-D-like Riccati case moves 0.01 %, and δW moves at the 1e-6 to 1e-9 level. The one metric that gets worse is Re(et[1]) on solovev_n1, from 5.6e-7 to 4.1e-6. (harness @ 9fc00dd)
• Migration: existing decks pick up ode_solver = "Vern7" and ode_abstol = 1e-8 automatically. ode_solver accepts Vern6, Vern7, Vern8, Vern9 or DP8; an unknown name is rejected before the equilibrium solve. To approximate old results, set ode_solver = "Vern9" and ode_abstol = 1e-6.

Every Euler-Lagrange solve (forward sweep, Riccati outer plasma, propagator chunks, Δ′ shooting) now runs with an explicit absolute tolerance and a configurable explicit Runge-Kutta method, defaulting to Vern7. Previously it used Vern9 at the OrdinaryDiffEq default absolute tolerance of 1e-6. BLAS is pinned to one thread inside the integration.

Why this is a bug fix, not a speed-up. With abstol = 1e-6, the solve error stopped responding to eulerlagrange_tolerance below about 1e-8. Two symptoms:

• Vern9 at reltol 1e-13 took the same 1667 steps as at 1e-10.
• The Riccati Δ′ depended on the chunk decomposition at 1.6e-5 (surface 1, auto vs auto+37 chunks). That dependence did not shrink when reltol went from 1e-10 to 1e-12, and it drops to 2e-8 once the chunk propagators get a real abstol.

Step counts rise, 20–97 % in the harness, because the tolerance is now actually enforced. Runtime stays roughly the same.

Regression report

regress --refs 20c466994,9fc00ddd5 --force over every case, on feynman (SLURM, 8 threads; develop tip vs this branch rebased onto it, both fresh, manifest pinned). Changed rows only; gal_resistive_pe is N/A on both refs (a pre-existing extraction gap on develop).

ggj_ray_q500i: 4 unchanged
Case: solovev_kinetic_nuzero
vacuum energy Re(ev[1])  1.072237e+01  1.072237e+01  2.060e-07 (0.00%)  ** CHANGED **
total energy Re(et[1])  2.243336e+00  2.243336e+00  9.054e-08 (0.00%)  ** CHANGED **
total energy Im(et[1])  -2.026025e+00  -2.026025e+00  8.209e-08 (0.00%)  ** CHANGED **
plasma energy Re(ep[1])  -8.479029e+00  -8.479030e+00  2.966e-07 (0.00%)  ** CHANGED **
ODE steps (saved)  608  1091  4.830e+02 (79.44%)  ** CHANGED **
total energy (all)  [32 elem]  [32 elem]  7.304e-06 (0.00%)  ** CHANGED **
ODE steps (total)  752  1356  6.040e+02 (80.32%)  ** CHANGED **
Summary: 7 changed, 7 unchanged

Case: diiid_slayer_n1
SLAYER γ_Hz [2/1,3/1,4/1]  [3 elem]  [3 elem]  1.489e-01 (0.01%)  ** CHANGED **
Summary: 1 changed, 16 unchanged

gal_resistive_pe: 8 missing/N/A
Case: solovev_kinetic_calculated
vacuum energy Re(ev[1])  1.043757e+01  1.043757e+01  1.784e-08 (0.00%)  ** CHANGED **
total energy Re(et[1])  1.873745e+00  1.873745e+00  6.129e-09 (0.00%)  ** CHANGED **
total energy Im(et[1])  -1.434288e+00  -1.434288e+00  2.188e-08 (0.00%)  ** CHANGED **
plasma energy Re(ep[1])  -8.563829e+00  -8.563829e+00  2.397e-08 (0.00%)  ** CHANGED **
ODE steps (saved)  599  1082  4.830e+02 (80.63%)  ** CHANGED **
total energy (all)  [32 elem]  [32 elem]  4.095e-06 (0.00%)  ** CHANGED **
ODE steps (total)  734  1335  6.010e+02 (81.88%)  ** CHANGED **
Summary: 7 changed, 7 unchanged

gal_resistive_diiid: 10 unchanged
Case: solovev_n1
total energy Re(et[1])  6.773398e-01  6.773430e-01  3.171e-06 (0.00%)  ** CHANGED **
total energy Im(et[1])  -1.099802e-03  -1.099238e-03  5.644e-07 (0.05%)  ** CHANGED **
plasma energy Re(ep[1])  -9.736133e+00  -9.736128e+00  4.441e-06 (0.00%)  ** CHANGED **
vacuum energy Re(ev[1])  1.041347e+01  1.041347e+01  1.270e-06 (0.00%)  ** CHANGED **
plasma energy (all)  [32 elem]  [32 elem]  2.124e-03 (0.00%)  ** CHANGED **
vacuum energy (all)  [32 elem]  [32 elem]  1.814e-03 (0.00%)  ** CHANGED **
total energy (all)  [32 elem]  [32 elem]  5.456e-04 (0.00%)  ** CHANGED **
ODE steps (saved)  388  655  2.670e+02 (68.81%)  ** CHANGED **
ODE steps (total)  618  997  3.790e+02 (61.33%)  ** CHANGED **
ca_left (checksum)  320262e6cac8...  7f8bce2b5070...  0.000e+00 (0.00%)  ** CHANGED **
Summary: 10 changed, 12 unchanged

ggj_reference: 4 unchanged
Case: solovev_multi_n
total energy Re(et[1])  -1.868943e+00  -1.868936e+00  6.186e-06 (0.00%)  ** CHANGED **
total energy Im(et[1])  -3.604869e-05  -3.625360e-05  2.049e-07 (0.57%)  ** CHANGED **
plasma energy (all)  [70 elem]  [70 elem]  9.428e-02 (0.00%)  ** CHANGED **
vacuum energy (all)  [70 elem]  [70 elem]  5.237e-02 (0.00%)  ** CHANGED **
total energy (all)  [70 elem]  [70 elem]  7.400e-02 (0.00%)  ** CHANGED **
ODE steps (saved)  36  51  1.500e+01 (41.67%)  ** CHANGED **
ODE steps (total)  1141  2101  9.600e+02 (84.14%)  ** CHANGED **
Summary: 7 changed, 8 unchanged

Case: diiid_n1
plasma energy (all)  [35 elem]  [35 elem]  1.286e-09 (0.00%)  ** CHANGED **
vacuum energy (all)  [35 elem]  [35 elem]  2.419e-10 (0.00%)  ** CHANGED **
total energy (all)  [35 elem]  [35 elem]  1.230e-09 (0.00%)  ** CHANGED **
ODE steps (saved)  2576  3308  7.320e+02 (28.42%)  ** CHANGED **
ODE steps (total)  4572  5491  9.190e+02 (20.10%)  ** CHANGED **
island half-widths  [5 elem]  [5 elem]  9.836e-07 (0.00%)  ** CHANGED **
Chirikov parameter  [5 elem]  [5 elem]  7.861e-05 (0.01%)  ** CHANGED **
dominant-coupling singular values  [3 elem]  [3 elem]  4.101e-05 (0.00%)  ** CHANGED **
PE toroidal torque  5.087465e-02  5.087465e-02  2.298e-12 (0.00%)  ** CHANGED **
resonant area-weighted field b^r  [5 elem]  [5 elem]  1.528e-08 (0.00%)  ** CHANGED **
Summary: 10 changed, 40 unchanged

Case: solovev_kinetic_ntv
NTV torque fgar [Re, Im]  [1 elem]  [1 elem]  1.896e-06 (1.22%)  ** CHANGED **
NTV ψ quadrature evaluations  840  810  3.000e+01 (3.57%)  ** CHANGED **
root-area-weighted total energy Re(et[1])  6.773398e-01  6.773430e-01  3.171e-06 (0.00%)  ** CHANGED **
Summary: 3 changed, 3 unchanged

Case: diiid_n1_riccati
delta prime (BVP diagonal)  [5 elem]  [5 elem]  3.358e-01 (0.01%)  ** CHANGED **
delta prime (raw side-major)  [10 elem]  [10 elem]  1.685e-01 (0.00%)  ** CHANGED **
edge coil response delta_coil  [10 elem]  [10 elem]  1.348e-06 (0.00%)  ** CHANGED **
total energy Re(et[1])  8.037196e-01  8.037194e-01  2.270e-07 (0.00%)  ** CHANGED **
plasma energy Re(ep[1])  -1.344751e+00  -1.344751e+00  5.246e-07 (0.00%)  ** CHANGED **
vacuum energy Re(ev[1])  2.148470e+00  2.148471e+00  2.976e-07 (0.00%)  ** CHANGED **
total energy (all)  [35 elem]  [35 elem]  8.199e-04 (0.00%)  ** CHANGED **
ODE steps (saved)  51  52  1.000e+00 (1.96%)  ** CHANGED **
ODE steps (total)  1647  3246  1.599e+03 (97.09%)  ** CHANGED **
Summary: 9 changed, 8 unchanged

Case: solovev_kinetic_multiion
NTV total (D+T+imp+e) [N·m]  1.315589e-04  1.333826e-04  1.824e-06 (1.39%)  ** CHANGED **
NTV Deuterium [N·m]  7.785958e-05  7.880733e-05  9.478e-07 (1.22%)  ** CHANGED **
NTV Tritium [N·m]  8.195889e-05  8.296629e-05  1.007e-06 (1.23%)  ** CHANGED **
NTV electron [N·m]  -2.825956e-05  -2.839106e-05  1.315e-07 (0.47%)  ** CHANGED **
total energy Re(et[1])  6.773398e-01  6.773430e-01  3.171e-06 (0.00%)  ** CHANGED **
Summary: 5 changed, 1 unchanged

efit_fixedbdy_separatrix: 5 unchanged

The diiid_slayer_n1 γ line (0.01 %) is the unseeded-triangulation tie-break noise that #463 seeds away.

Is it closer to converged?

The decks at eulerlagrange_tolerance = 1e-7 were scored against two references at reltol 1e-12 / abstol 1e-14 from different solver families: Vern9 and DP8. The numbers below are the error against the DP8 reference. The two references agree to 2e-5 or better on every quantity scored here.

quantity	develop	this PR
NTV torque (kinetic NTV; multi-ion D, T, total)	1.5–1.7 %	0.30–0.36 %
NTV electron (multi-ion)	5.8e-3	1.2e-3
multi-n energies (all elements)	3.7e-6 – 1.5e-5	4e-8 – 1.7e-7
multi-n Re(et[1])	3.5e-6	1.8e-7
solovev_n1 energies (all elements)	2.5e-8 – 3.2e-7	8.6e-9 – 1.1e-7
solovev_n1 Re(et[1])	5.6e-7	4.1e-6
kinetic-calculated and ν=0 decks, energies	≤ 4e-8	≤ 3e-8 (mixed)

• Im(et[1]) is not scored on solovev_n1 or multi-n. On solovev_n1 the two references disagree by 3.7e-3, so that value is not converged at this resolution.
• No Vern9 reference on the kinetic-calculated decks. At reltol 1e-12 Vern9 falls into its tiny-step collapse there: over 80 GB of memory before it was stopped, and an HDF5 write failure in the harness. DP8 at the same tolerance ran normally, so those rows use DP8 alone. This is the same Vern9 breakdown the study saw in forward mode, and one more reason not to default to it.
• Residual NTV error. The 0.3 % that remains is the 1e-7 tolerance of those decks. Tightening their eulerlagrange_tolerance would be a separate, deck-level decision.

Evidence from the integrator study

Reproducing the checks with what is in the repository.

• The converged-reference table above. Copy a harness case, give it [overrides] with "ForceFreeStates.ode_solver" = "DP8", "ForceFreeStates.eulerlagrange_tolerance" = 1e-12 and "ForceFreeStates.ode_abstol" = 1e-14, and run both cases at one commit, e.g. regress --cases solovev_kinetic_ntv,solovev_kinetic_ntv_ref --refs <sha>. Run the same pair at a develop commit for the "develop" column.
• The chunk-decomposition floor. On examples/DIIID-like_ideal_example, compare solve(equil, Riccati(; nchunks=0)) with nchunks = 90: delta_prime.matrix[1,1] differs by 1.6e-5 relative on develop, the same at eulerlagrange_tolerance = 1e-12, and by 2e-8 here.

The integrator study behind the method choice ran 20 solvers over 3 reltols × 4 abstols. Its drivers are not part of this PR, and the tables below are its output. On the shipped DIII-D-like case, with errors relative to a Vern9 reltol 1e-13 / abstol 1e-15 reference:

setting	integration	FFS stage	lowest-5 δW error	Δ′ error
develop: Vern9, 1e-10 / 1e-6	1.50 s	9.9 s	5e-8	1e-4
this PR: Vern7, 1e-10 / 1e-8	1.84 s	8.6 s	2e-9	1.5e-5
Vern9, 1e-10 / 1e-8	2.48 s	10.9 s	2e-9	1e-4

At equal tolerances Vern7 is about 1.8× faster than Vern9, because Vern9 rejects most of its steps near the rational surfaces. So the fix costs little time. The 8.6 s vs 9.9 s stage time comes largely from the BLAS pin, not from the method.

Ten IDA_run geqdsks (EFIT, CAKE, TokaMaker; 65×65 to 257×257; H- and L-mode), reltol 1e-10 / abstol 1e-10 for every solver, 4 threads:

shot / time	N	msing	STRIDE wall	Vern9 integ.	Vern7 integ.	Vern7 Δ′ err / Vern9 Δ′ err
148798 / 3306	33	2	3.8 s	4.78 s	2.56 s	1.9e-6 / 3.5e-6
153833 / 3450	34	4	3.5 s	5.51 s	3.00 s	4.8e-6 / 7.0e-6
169510 / 3000 CAKE	34	4	3.8 s	6.06 s	3.34 s	4.1e-6 / 7.9e-6
178305 / 1800 (65×65)	36	5	4.5 s	8.83 s	4.85 s	1.7e-6 / 3.0e-6
179633 / 3200 L-mode	34	4	3.6 s	2.57 s	1.50 s	2.3e-6 / 6.9e-6
189392 / 4080 L-mode	34	4	3.6 s	2.42 s	1.39 s	1.2e-6 / 1.9e-6
201586 / 4200 (65×65)	40	10	13.5 s	13.73 s	8.81 s	3.2e-7 / 3.6e-7
204441 / 4409	37	7	7.8 s	5.80 s	3.92 s	3.3e-6 / 6.7e-6
TkMkr D3D H-mode	35	5	5.5 s	4.33 s	2.86 s	1.5e-5 / 9.7e-5

Vern7 vs Vern9 integration speedup: min 1.48×, median 1.81×, max 4.7×; δW within 2e-10 of the reference on all; Δ′ closer on 10/10.

Notes for reviewers

• BLAS pin measured. It covers the whole integration because the RHS works on mpert×mpert blocks, where multithreaded BLAS costs more in synchronization than it saves. The one large-matrix call inside it, the edge-dW vacuum spline, takes 0.35 s pinned and 0.36–0.40 s with 8 BLAS threads, so the wide pin costs nothing there.
• ode_solver is limited to the benchmarked methods. Tsit5 lost Δ′ accuracy 10–100× on three corpus cases, and VCABM was not robust on the 10-surface case. The Δ′ shooting solves keep their per-column absolute tolerance (scale × reltol) rather than ode_abstol, as the field docstring now says.
• Rebuilt after an adversarial review (2026-09-23).
  • Rebased onto current develop as the original commit plus review fixes.
  • The benchmark scripts and PNGs are dropped from the PR, and the body no longer claims a speed-up or an unfinished corpus sweep.
  • ode_solver is validated early, with tests.
  • The equation blocks the formatter flattened in the Riccati docstrings are fenced.
• Not sequenced after Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots #398. The abstol defect is independent of the equilibrium resampling Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots #398 fixes, and Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots #398 re-baselines these numbers either way.
• Reverse-shear double-rational geqdsks (q crossed twice) are pathological for the singular-surface handling and are a separate bugfix, not an integrator issue.

NO MERGE WITHOUT HUMAN REVIEW. This PR is a draft and needs a named human reviewer and an assignee before it can be considered for merge. Non-negotiable.

🤖 Generated with Claude Code

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

[d-burg]

Closing in favour of #480, which fixes the same unset absolute tolerance with a per-column value refreshed every step and keeps Vern9.

Two things found while diagnosing this PR's tolerance study are tracked separately: the wrong-solution removal at ideal crossings (#489) and the reduction trigger missing near-parallel solutions (#490). The selectable ode_solver option here is not in #480; it can come back as its own PR if a Vern7 comparison at the deck tolerances is wanted.