Equilibrium - BUGFIX! - 🚨 Tie the equilibrium ODE absolute tolerance to etol

[d-burg]

Release note

• Audience: users
• Numerical impact: DIII-D-like results move; Solovev, GGJ and kinetic Solovev cases are unchanged. Riccati Δ′ real parts move 0.1–1.3 % on the auto grid and at most 0.02 % on fixed grids (0.2 % at the 6/1). et[1] moves 0.07 %, the 2/1 SLAYER growth rate 0.6 % (3/1 0.20 %, 4/1 0.28 %), the NTV torque 0.12 % and the Galerkin Δ_coil block 0.99 %. The harness's worst-element figure for the Riccati Δ′ diagonal is 3.5–15.6 %; that is the imaginary part of the 5/1 element, which is not converged on either side. (harness @ 6d6a493)
• Migration: a non-positive or non-finite etol is now a construction-time error. Decks with a valid etol need no change; DIII-D-like results will differ from earlier runs. An equilibrium ODE solve that stops before the end of its range (step limit, step-size underflow, instability) is now an error instead of a silently truncated solution, including inside the TJ-analytic ν root-find.

The equilibrium field-line, LAR and TJ-analytic integrations used reltol = etol with a hardcoded abstol, so tightening etol past 1e-8 did nothing where the solution is small. The absolute tolerance is now min(etol, previous hardcode): tighter decks converge as requested, and looser decks keep exactly their previous value.

Equilibrium of the DIII-D-like example against etol, each relative to the tightest run. On develop nothing improves past 1e-8; with this PR β_N and the q profile keep converging. At the default etol = 1e-10 the equilibrium setup takes 3.5 s instead of 2.5 s.

Left: change in Re Δ′ at each rational surface, develop to this PR, on three grids. Right: other tracked quantities on the same case.

Regression report

regress --refs 6fc89d021,6d6a4932d --force (all cases) on feynman: current develop against this branch with develop merged.

cases	result
Solovev, GGJ, kinetic Solovev (8 cases)	unchanged
diiid_n1_riccati (auto grid)	Re Δ′ moves 0.14–1.33 % per surface
fixed grids, 128 and 256 points	Re Δ′ moves ≤ 0.02 %; 0.22 % at the 6/1 on 256 points
diiid_n1	et[1] 0.07 %, NTV torque 0.12 %, resonant field magnitudes ≤ 0.06 %, equilibrium scalars ≤ 2e-5
diiid_slayer_n1	γ: 2/1 0.61 %, 3/1 0.20 %, 4/1 0.28 %
gal_resistive_diiid	Δ′ matrix norm 0.05 %, Δ_coil block 0.99 %
diiid_error_field	correction overlap 0.01 %, NTV torques 0.02–0.12 %
diiid_multi_n, efit_fixedbdy_separatrix	≤ 0.03 %
diiid_n1_riccati_precompiled	the same lines as diiid_n1_riccati
gal_resistive_pe	extracts nothing on either ref (a fault of that case, not this PR)

The Riccati Δ′ real parts are bit-identical, on both refs, to the run before the latest develop merge and the review fixes. The 2/1 γ figure depends on develop's unseeded triangulation, which moves it by about 0.1 % from run to run.

An earlier run, regress --refs 79ab50427,e875c3cbb --force (all cases), isolates the first return-code commit: every case is unchanged apart from the 0.01 % SLAYER γ line that is the same triangulation noise. So the check moves nothing, and no equilibrium solve in any harness case stops early.

Notes for reviewers

• Not covered by ForceFreeStates - BUGFIX! - Scale the ODE absolute tolerance to each solution column #480: that PR fixes the ForceFreeStates solves only; this is the equilibrium solve.
• Fortran comparison: Fortran's field-line integration already scales its absolute tolerance with etol. On matched fixed ldp grids (mpsi 128 and 256), this branch moves the converged Δ′ elements by 0.003–0.2 %, each toward Fortran STRIDE.
• Overlap: Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots #398 edits the same two files and may conflict.
• Return-code check added, and it found a bug in this branch. The LAR profile and TJ-analytic shaping solves now error unless the integrator finished. The arc-length field-line solve must end on its midplane callback; reaching the end of its range means the field line never closed. The direct field-line solve keeps develop's own guard, which is stricter.
  • With abstol following etol, the LAR profile solve needs about 23,500 Rosenbrock23 steps at the default etol = 1e-10 (2,463 with the old fixed 1e-8). It was stopping at its hard-coded 10,000-step limit at r = 0.42 a, and the equilibrium was built from the truncated table; the existing LAR unit test still passed.
  • The fixed 10,000-step limit is removed, so the integrator's default of 1e6 applies. No tolerance changed. No harness case uses eq_type = "lar", which is why the earlier report did not show it.
• Changes since the review (commit 6d6a4932d). All seven review suggestions are in.
  • The ν root-find falls back to ν = qa/qc only when the root-find itself has no bracket or does not converge. Any other error, including a trial shaping solve that stopped early, propagates. A test forces a failed trial solve and checks that it does.
  • check_equil_solve accepts only Success and Terminated; a test covers a StalledSuccess solution.
  • The TJ shaping solve loses its 10,000-step limit, as the LAR solve did. The redundant abstol arguments and the single-use constant are gone, and the etol docstring is shortened.
  • The tolerance helper and the solve check moved from EquilibriumTypes.jl to a new SolveTolerances.jl; no existing Equilibrium file was a natural home.
• The torque sign flip reported earlier is gone. It came from an order-dependent tie-break in the m = 0 axis start on develop, which ForceFreeStates - BUGFIX! - Default the axis start to the fixed initialization #514 has since removed. Against current develop the perturbed-equilibrium toroidal torque moves 0.03 %, the PE surface energy 0.04 % and Im(et[1]) 0.07 %. The 50 % and 65 % rotation-scan lines on diiid_error_field are gone as well; that line now moves 0.03 %.

---

NO MERGE WITHOUT HUMAN REVIEW. The approval on this PR predates the develop merge and the review fixes, so the reviewer needs to see the new head before it is merged.

🤖 Generated with Claude Code

[Copilot]

Copilot review overview

🟡 Changes recommended

TJ root finding currently catches and suppresses newly detected ODE integration failures.

Review effort: Balanced
Findings: 1

Open (1)

• Broad catch swallows TJ solve failures during fallback · New

What changed in this PR

Ties equilibrium ODE absolute tolerances to etol, validates tolerance inputs, and rejects incomplete integrations.

Changes:

• Adds shared tolerance and solve-status helpers.
• Applies tolerance scaling across direct, arc-length, LAR, and TJ integrations.
• Adds validation and failure-handling tests.

One unresolved issue remains: TJ root finding can swallow ODE failures and silently fall back.

NO MERGE WITHOUT THIRD-PARTY HUMAN REVIEW

File	Description
src/​Equilibrium/​EquilibriumTypes.jl	Adds validation and shared helpers.
src/​Equilibrium/​DirectEquilibrium.jl	Scales direct integration tolerance.
src/​Equilibrium/​DirectEquilibriumArcLength.jl	Scales positional tolerance and validates closure.
src/​Equilibrium/​AnalyticEquilibrium.jl	Updates LAR/TJ tolerances and completion checks.
test/​runtests_equil.jl	Tests validation and solve status.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[jhalpern30]

Claude's addition to this - a bit in the weeds for me, probably just combine the recommendation:
This catch swallows the new solve check. check_equil_solve throws ErrorException on MaxIters, step-size underflow, or an unstable right-hand side, and find_zero lets that exception out. Roots only raises Roots.ConvergenceFailed (no convergence) or ArgumentError whose message contains "not a bracketing interval" (the bracket has no sign change). Every other exception should propagate.

As written, one trial ν that stops early is logged as a root-find miss and replaced with ν = qa/qc. The final tj_analytic_shape_solve at that fallback can then succeed, so the equilibrium is built at the wrong ν.

catch err
    fallback = err isa Roots.ConvergenceFailed ||
               (err isa ArgumentError && occursin("not a bracketing interval", err.msg))
    fallback || rethrow()
    @warn "ν root-find failed for TJ-analytic equilibrium; falling back to lowest-order ν = qa/qc" exception=(err, catch_backtrace())
    nu_guess
end

exception=(err, catch_backtrace()) is what prints the stack; error=err only stores the object.

[d-burg]

Done with the patch as written. The fallback is taken only for Roots.ConvergenceFailed or the no-bracket ArgumentError, everything else rethrows, and the warning now carries the backtrace. I added a test that forces a failed trial solve (NaN abstol) and checks the error comes out of tj_analytic_find_nu.

One caveat: the no-bracket case is matched on Roots' message text, which I checked against the pinned Roots 2.3.0. If that text ever changes, the case will rethrow, so it fails loudly and not silently.

[jhalpern30]

Looks good, only minor suggestions. Interesting how this is now the second abstol that was dropped but then turned out to be important to results

[jhalpern30]

Claude says to drop the maxiters here, since the default of 1e6 should be fine - just like you did for the main LAR integrator above this.

[d-burg]

Done in 6d6a493. It's gone from both solve calls in tj_analytic_shape_solve, so the integrator's default applies, same as the LAR solve.

[jhalpern30]

This extra argument isn't actually doing anything, you should be able to just pass in the reltol and get the default keyword argument behavior. This applies to 461, 461, and 601, and more clearly differentiates the alternate behavior you have for 608

[d-burg]

Done. The three call sites now pass only reltol. The dense solve keeps its explicit abstol because its ceiling is tighter.

[jhalpern30]

This is only used once below. Why not just hardcode 1e-10 there with a comment for why its there? Or something more obvious

[d-burg]

Done. The constant is removed and the 1e-10 is inline, with a comment saying why that solve needs a tighter ceiling than the others.

[jhalpern30]

From Claude, seems to make sense:
successful_retcode is also true for FloatingPointLimit, StalledSuccess, ExactSolutionLeft, and ExactSolutionRight. Those are nonlinear-solver exits; an equilibrium integral that stops on one of them has not reached the end of its range. ReturnCode is already in scope in this module:

retcode = sol.retcode
retcode == ReturnCode.Success || retcode == ReturnCode.Terminated ||
    error("Equilibrium ODE solve ($what) stopped early at t = $(sol.t[end]) with retcode $(retcode).")

Leave the arc-length check (DirectEquilibriumArcLength.jl:119) as a separate == ReturnCode.Terminated test. There, Success means the midplane callback never fired.

[d-burg]

Done as suggested. check_equil_solve accepts only Success and Terminated, and the arc-length check is still its own == Terminated test. I added a test that a StalledSuccess solution is rejected.

[jhalpern30]

Sure, I get having the check below to ensure a reasonable input. But does this really need to be included in the docstring lol, I think its fairly straightforward that a tolerance is positive and finite

[d-burg]

Agreed, removed. The constructor check stays.

[jhalpern30]

These additions seem a bit out of place in EquilibriumTypes... is there a better place for them?

[d-burg]

Moved both to a new src/Equilibrium/SolveTolerances.jl, included right after EquilibriumTypes.jl. I didn't find an existing Equilibrium file they fit in. Happy to move them again if you have one in mind.

[jhalpern30]

Claude's addition to this - a bit in the weeds for me, probably just combine the recommendation:
This catch swallows the new solve check. check_equil_solve throws ErrorException on MaxIters, step-size underflow, or an unstable right-hand side, and find_zero lets that exception out. Roots only raises Roots.ConvergenceFailed (no convergence) or ArgumentError whose message contains "not a bracketing interval" (the bracket has no sign change). Every other exception should propagate.

As written, one trial ν that stops early is logged as a root-find miss and replaced with ν = qa/qc. The final tj_analytic_shape_solve at that fallback can then succeed, so the equilibrium is built at the wrong ν.

catch err
    fallback = err isa Roots.ConvergenceFailed ||
               (err isa ArgumentError && occursin("not a bracketing interval", err.msg))
    fallback || rethrow()
    @warn "ν root-find failed for TJ-analytic equilibrium; falling back to lowest-order ν = qa/qc" exception=(err, catch_backtrace())
    nu_guess
end

exception=(err, catch_backtrace()) is what prints the stack; error=err only stores the object.

[d-burg]

All seven suggestions are in 6d6a493, together with a merge of current develop. I reran the full harness against develop and updated the numbers in the description.

The Riccati Δ′ real parts are bit-identical to the run you reviewed. Several other numbers shrank (NTV torque 0.88 % → 0.12 %, and the torque sign flip is gone) because #514 removed the axis-start tie-break that was inflating them. Could you take another look at the new head?