From 9334c22853ff81824c19deb01630bd07e45d7232 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Wed, 23 Sep 2026 12:32:57 -0400
Subject: [PATCH 1/8] Equilibrium - BUGFIX! - Tie the equilibrium ODE absolute
 tolerance to etol

The field-line, LAR and TJ-analytic integrations passed reltol=etol but a
hardcoded abstol, which dominates wherever the solution is small, so tightening
etol past 1e-8 left q0 wandering at ~1e-9. Use abstol = min(etol, previous
hardcode): tighter decks now converge as requested, and looser decks (the
Solovev examples at etol=1e-7) keep exactly their previous abstol.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/AnalyticEquilibrium.jl        | 6 +++---
 src/Equilibrium/DirectEquilibrium.jl          | 2 +-
 src/Equilibrium/DirectEquilibriumArcLength.jl | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl
index 767f2b718..f9c19b1ac 100644
--- a/src/Equilibrium/AnalyticEquilibrium.jl
+++ b/src/Equilibrium/AnalyticEquilibrium.jl
@@ -140,7 +140,7 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
 
     prob = ODEProblem(dydr, y0, tspan, p)
 
-    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=1e-8, maxiters=10000, dense=false)
+    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=min(equil_input.etol, 1e-8), maxiters=10000, dense=false)
 
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
@@ -372,7 +372,7 @@ downstream Hₙ / ψ splines sit on uniform nodes); leave it `nothing` for
 the default adaptive save pattern used by `tj_analytic_run`.
 """
 function tj_analytic_shape_solve(p::TJAnalyticShapeParams, nu::Float64;
-    reltol::Float64=1e-7, abstol::Float64=1e-8,
+    reltol::Float64=1e-7, abstol::Float64=min(reltol, 1e-8),
     saveat=nothing)
     rhs_params = (; p.a, p.B0, p.qc, p.mu, p.pc, p.epsa2, nu=nu)
     prob = ODEProblem(tj_analytic_shape_rhs!, tj_analytic_shape_initial(p, nu), (p.r0, p.a), rhs_params)
@@ -599,7 +599,7 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo
     # the (R, Z) → (r, w) Newton iteration hits spline interpolation artifacts.
     dense_r = collect(range(p.r0, p.a; length=1024))
     sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol,
-        abstol=1e-10, saveat=dense_r)
+        abstol=min(equil_input.etol, 1e-10), saveat=dense_r)
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
 
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 0bc9d26fa..008268c1d 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -291,7 +291,7 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
     callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false))
 
     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
-    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=false)
+    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=min(equil_config.etol, 1e-8), dt=2π / 200, adaptive=true, dense=false)
 
     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
     return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
index 32b93a0ac..984906d67 100644
--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
+++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
@@ -112,7 +112,7 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     prob = ODEProblem{true}(arclength_fieldline_der!, u0, (0.0, 1.0e4), params)
     # Tight tolerances on position (y[1:2]); integrals (y[3:5]) effectively unconstrained near x-points
     reltol_vec = [equil_config.etol, equil_config.etol, 1e20, 1e20, 1e20]
-    abstol_vec = [1e-8, 1e-8, 1e20, 1e20, 1e20]
+    abstol_vec = [min(equil_config.etol, 1e-8), min(equil_config.etol, 1e-8), 1e20, 1e20, 1e20]
     sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec,
         dt=2π / 200, adaptive=true, dense=false)
 

From d0859be98654bb0178ac974d9eb276a981da46b2 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 11:02:59 -0400
Subject: [PATCH 2/8] Equilibrium - REFACTOR - Name the equilibrium ODE abstol
 ceiling and pass it explicitly

Replace the repeated min(etol, 1e-8) literal with equil_abstol(etol) capped at
EQUIL_ABSTOL_MAX, name the dense TJ-analytic ceiling TJ_DENSE_ABSTOL_MAX, pass
abstol explicitly from tj_analytic_run and tj_analytic_run_direct into the nu
root-find, and document that etol now governs abstol. Same values at every site.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/AnalyticEquilibrium.jl        | 22 +++++++++++--------
 src/Equilibrium/DirectEquilibrium.jl          |  2 +-
 src/Equilibrium/DirectEquilibriumArcLength.jl |  2 +-
 src/Equilibrium/EquilibriumTypes.jl           | 15 ++++++++++++-
 4 files changed, 29 insertions(+), 12 deletions(-)

diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl
index f9c19b1ac..f0c924a1c 100644
--- a/src/Equilibrium/AnalyticEquilibrium.jl
+++ b/src/Equilibrium/AnalyticEquilibrium.jl
@@ -140,7 +140,7 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
 
     prob = ODEProblem(dydr, y0, tspan, p)
 
-    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=min(equil_input.etol, 1e-8), maxiters=10000, dense=false)
+    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol), maxiters=10000, dense=false)
 
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
@@ -369,10 +369,11 @@ end
 Integrate the TJ-analytic shape ODE for the given ν.  Pass `saveat` to collect
 output on a prescribed dense grid (used by `tj_analytic_run_direct` so the
 downstream Hₙ / ψ splines sit on uniform nodes); leave it `nothing` for
-the default adaptive save pattern used by `tj_analytic_run`.
+the default adaptive save pattern used by `tj_analytic_run`. `reltol` and `abstol` go to the
+Vern9 solve; `abstol` defaults to `equil_abstol(reltol)`.
 """
 function tj_analytic_shape_solve(p::TJAnalyticShapeParams, nu::Float64;
-    reltol::Float64=1e-7, abstol::Float64=min(reltol, 1e-8),
+    reltol::Float64=1e-7, abstol::Float64=equil_abstol(reltol),
     saveat=nothing)
     rhs_params = (; p.a, p.B0, p.qc, p.mu, p.pc, p.epsa2, nu=nu)
     prob = ODEProblem(tj_analytic_shape_rhs!, tj_analytic_shape_initial(p, nu), (p.r0, p.a), rhs_params)
@@ -393,9 +394,9 @@ O(εa²) correction relative to the lowest-order guess ν = qa/qc, which
 matters for the TJ-analytic benchmark at large ε.  Falls back to the
 lowest-order ν if the bracket search diverges.
 """
-function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; reltol::Float64=1e-7)
+function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; reltol::Float64=1e-7, abstol::Float64=equil_abstol(reltol))
     function q2_edge(nu::Float64)
-        sol = tj_analytic_shape_solve(p, nu; reltol)
+        sol = tj_analytic_shape_solve(p, nu; reltol, abstol)
         g2end = sol.u[end][2]
         f3end = sol.u[end][5]
         f1end = tj_analytic_f1(1.0, nu, p.qc)
@@ -455,8 +456,8 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig)
     epsa2 = p.epsa2
     p00_phys = B0^2 * epsa2 * pc          # μ₀P = B₀²·εa²·p₂ at axis
 
-    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol)
-    sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol)
+    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol))
+    sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol))
 
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
@@ -549,6 +550,9 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig)
     return InverseRunInput(equil_input, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, R0, 0.0, psio, nothing)
 end
 
+# Tighter abstol ceiling for the dense-saveat TJ-analytic solve that feeds the (R, Z) → (r, w) Newton inversion.
+const TJ_DENSE_ABSTOL_MAX = 1e-10
+
 """
     tj_analytic_run_direct(equil_input, tj_input; nrbox=257, nzbox=257, rc=1.2)
 
@@ -592,14 +596,14 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo
     p00_phys = B0^2 * epsa2 * pc
 
     # ν root-find (cf. Fitzpatrick TJ's Setnu): q₂(1) = qa_target.
-    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol)
+    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol))
 
     # Dense saveat so the downstream splines (H₁, g₂, f₃, ψ) are evaluated on
     # a fine uniform r grid rather than the ~30 adaptive Vern9 steps — otherwise
     # the (R, Z) → (r, w) Newton iteration hits spline interpolation artifacts.
     dense_r = collect(range(p.r0, p.a; length=1024))
     sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol,
-        abstol=min(equil_input.etol, 1e-10), saveat=dense_r)
+        abstol=equil_abstol(equil_input.etol, TJ_DENSE_ABSTOL_MAX), saveat=dense_r)
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
 
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 008268c1d..266fa3cd5 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -291,7 +291,7 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
     callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false))
 
     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
-    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=min(equil_config.etol, 1e-8), dt=2π / 200, adaptive=true, dense=false)
+    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=equil_abstol(equil_config.etol), dt=2π / 200, adaptive=true, dense=false)
 
     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
     return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
index 984906d67..e587f68bd 100644
--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
+++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
@@ -112,7 +112,7 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     prob = ODEProblem{true}(arclength_fieldline_der!, u0, (0.0, 1.0e4), params)
     # Tight tolerances on position (y[1:2]); integrals (y[3:5]) effectively unconstrained near x-points
     reltol_vec = [equil_config.etol, equil_config.etol, 1e20, 1e20, 1e20]
-    abstol_vec = [min(equil_config.etol, 1e-8), min(equil_config.etol, 1e-8), 1e20, 1e20, 1e20]
+    abstol_vec = [equil_abstol(equil_config.etol), equil_abstol(equil_config.etol), 1e20, 1e20, 1e20]
     sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec,
         dt=2π / 200, adaptive=true, dense=false)
 
diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index ac5a212fc..a110a89fc 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -38,7 +38,8 @@ specified in the input.
   - `mtheta::Int` - Number of poloidal grid points
   - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to
     meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped)
-  - `etol::Float64` - Error tolerance for equilibrium solver
+  - `etol::Float64` - Relative tolerance of the equilibrium ODE solves; their absolute tolerance is
+    `equil_abstol(etol)`, i.e. `etol` capped at `EQUIL_ABSTOL_MAX`
   - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations)
   - `use_galgrid::Bool` - Use the same grid as galerkin method
 """
@@ -154,6 +155,18 @@ specified in the input.
     end
 end
 
+"""
+Ceiling on the absolute tolerance of the equilibrium ODE solves; the historical fixed value, kept so a loose `etol` never loosens abstol.
+"""
+const EQUIL_ABSTOL_MAX = 1e-8
+
+"""
+    equil_abstol(etol, cap=EQUIL_ABSTOL_MAX)
+
+Absolute tolerance of an equilibrium ODE solve: follows `etol` but is never looser than `cap`.
+"""
+equil_abstol(etol::Real, cap::Real=EQUIL_ABSTOL_MAX) = min(etol, cap)
+
 """
 Outer constructor for EquilibriumConfig from a parsed TOML dictionary
 """

From bd4a5fb4bf48d5697413486a406b6a952fe7a579 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 11:09:38 -0400
Subject: [PATCH 3/8] Repo - MINOR - Apply the formatter to the audit fixes

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/DirectEquilibriumArcLength.jl | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
index e587f68bd..1ed2d72bb 100644
--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
+++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
@@ -64,11 +64,11 @@ end
 Arc-length-parameterized flux surface integration. Drop-in replacement for
 `direct_fieldline_int` with identical return format:
 
-- `y_out[:, 1]`: geometric angle η ∈ 0 to 2π (CCW from outboard midplane)
-- `y_out[:, 2]`: accumulated ∫dl/Bp
-- `y_out[:, 3]`: rfac = √((R−ro)² + (Z−zo)²)
-- `y_out[:, 4]`: accumulated ∫dl/(R²Bp)
-- `y_out[:, 5]`: accumulated ∫jac·dl/Bp
+  - `y_out[:, 1]`: geometric angle η ∈ 0 to 2π (CCW from outboard midplane)
+  - `y_out[:, 2]`: accumulated ∫dl/Bp
+  - `y_out[:, 3]`: rfac = √((R−ro)² + (Z−zo)²)
+  - `y_out[:, 4]`: accumulated ∫dl/(R²Bp)
+  - `y_out[:, 5]`: accumulated ∫jac·dl/Bp
 
 The ODE is terminated by a `ContinuousCallback` that detects the return to the
 outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
@@ -144,4 +144,3 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     # bfield at the starting point carries F and P for the surface-averaged quantities
     return y_out, bfield
 end
-

From 665357aa4a0ceb9d70aca820d71946ead3b16edf Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 15:38:40 -0400
Subject: [PATCH 4/8] Equilibrium - BUGFIX - Reject a non-positive or
 non-finite etol when building the equilibrium config

etol sets both the relative and (capped) absolute tolerance of the equilibrium ODE solves, so zero, negative, NaN or Inf is now a construction-time error rather than a failed or meaningless solve. Adds a test in runtests_equil.jl.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/EquilibriumTypes.jl | 5 +++--
 test/runtests_equil.jl              | 7 +++++++
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index a110a89fc..d5e7e5673 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -38,8 +38,8 @@ specified in the input.
   - `mtheta::Int` - Number of poloidal grid points
   - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to
     meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped)
-  - `etol::Float64` - Relative tolerance of the equilibrium ODE solves; their absolute tolerance is
-    `equil_abstol(etol)`, i.e. `etol` capped at `EQUIL_ABSTOL_MAX`
+  - `etol::Float64` - Relative tolerance of the equilibrium ODE solves (positive and finite); their
+    absolute tolerance is `equil_abstol(etol)`, i.e. `etol` capped at `EQUIL_ABSTOL_MAX`
   - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations)
   - `use_galgrid::Bool` - Use the same grid as galerkin method
 """
@@ -144,6 +144,7 @@ specified in the input.
         else
             error("Cannot recognize jac_type = $(jac_type)")
         end
+        (isfinite(etol) && etol > 0) || error("etol = $etol must be a positive, finite relative tolerance")
         if psihigh > 1.0
             @warn "psihigh = $psihigh exceeds 1.0 (separatrix); clamping to 1.0"
         end
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index af3cd5743..836344073 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -129,6 +129,13 @@
         @test isapprox(q_axis, 1.05; rtol=0.02)
     end
 
+    @testset "etol must be positive and finite" begin
+        for etol in (0.0, -1e-8, NaN, Inf)
+            @test_throws ErrorException GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(; etol=etol)
+        end
+        @test GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(; etol=1e-8).etol == 1e-8
+    end
+
     @testset "Deprecated TOML keys are dropped, not fatal" begin
         # Removed control knobs must keep old gpec.toml decks (and older gpec.h5 replays,
         # whose stored TOML blob goes through the same path) parsing with a warning.

From 79ab504277cbd8ae7b9157856adee4fcd3c573d9 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 15:43:23 -0400
Subject: [PATCH 5/8] Repo - MINOR - Apply the formatter to the audit fixes

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 test/runtests_equil.jl | 52 +++++++++++++++++++++---------------------
 1 file changed, 26 insertions(+), 26 deletions(-)

diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 836344073..c1fdebe20 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -150,7 +150,7 @@
     @testset "EFIT Method Consistency" begin
         # All three methods solve the same equilibrium — q-profiles should broadly agree.
         # Tolerance is 10% to allow for method-specific discretisation differences.
-        q_efit      = plasma_eq_efit.profiles.q_spline.y
+        q_efit = plasma_eq_efit.profiles.q_spline.y
         q_arclength = plasma_eq_arclength.profiles.q_spline.y
         q_inversion = plasma_eq_inversion.profiles.q_spline.y
 
@@ -249,7 +249,7 @@
         b0exp = 7.4  # CHEASE normalization field [T]
 
         B_nodes_binary = plasma_eq_binary.eqfun_B.nodal_derivs.partials[1, :, :]
-        B_nodes_ascii  = plasma_eq_ascii.eqfun_B.nodal_derivs.partials[1, :, :]
+        B_nodes_ascii = plasma_eq_ascii.eqfun_B.nodal_derivs.partials[1, :, :]
 
         # B field must be finite and positive everywhere
         @test all(isfinite, B_nodes_binary)
@@ -265,7 +265,7 @@
 
         # q must be finite, positive, and in a physically reasonable range
         q_binary = plasma_eq_binary.profiles.q_spline.y
-        q_ascii  = plasma_eq_ascii.profiles.q_spline.y
+        q_ascii = plasma_eq_ascii.profiles.q_spline.y
         @test all(isfinite, q_binary)
         @test all(isfinite, q_ascii)
         @test all(>(0), q_binary)
@@ -338,7 +338,7 @@
         end
 
         @testset "sol_run clamps p0fac to ≥ 1" begin
-            equil_inputs, sol_inputs = make_inputs(p0fac=0.5)
+            equil_inputs, sol_inputs = make_inputs(; p0fac=0.5)
             dri = GeneralizedPerturbedEquilibrium.Equilibrium.sol_run(equil_inputs, sol_inputs)
             @test all(dri.sq_in.y[:, 2] .>= 0)  # no negative pressures
         end
@@ -367,7 +367,7 @@
         end
 
         @testset "sol_run spline integrity" begin
-            equil_inputs, sol_inputs = make_inputs(mr=6, mz=5, ma=3)
+            equil_inputs, sol_inputs = make_inputs(; mr=6, mz=5, ma=3)
             dri = GeneralizedPerturbedEquilibrium.Equilibrium.sol_run(equil_inputs, sol_inputs)
             sq = dri.sq_in
             psi = dri.psi_in
@@ -389,7 +389,7 @@
         end
 
         @testset "sol_run 2D psi field properties" begin
-            equil_inputs, sol_inputs = make_inputs(mr=3, mz=3)
+            equil_inputs, sol_inputs = make_inputs(; mr=3, mz=3)
             dri = GeneralizedPerturbedEquilibrium.Equilibrium.sol_run(equil_inputs, sol_inputs)
             psi = dri.psi_in
 
@@ -418,13 +418,13 @@
         @testset "sol_run extreme inputs" begin
             # minimal grid (CubicInterpolant requires at least 4 points for extrap BC)
             # mr=3, mz=3 creates 4-point grids (mr+1 points)
-            equil_inputs, sol_inputs = make_inputs(mr=3, mz=3, ma=3)
+            equil_inputs, sol_inputs = make_inputs(; mr=3, mz=3, ma=3)
             dri = GeneralizedPerturbedEquilibrium.Equilibrium.sol_run(equil_inputs, sol_inputs)
             @test length(dri.psi_in_xs) == 4
             @test length(dri.psi_in_ys) == 4
 
             # very high aspect ratio
-            equil_inputs, sol_inputs = make_inputs(e=0.8, a=0.1, r0=10.0)
+            equil_inputs, sol_inputs = make_inputs(; e=0.8, a=0.1, r0=10.0)
             dri = GeneralizedPerturbedEquilibrium.Equilibrium.sol_run(equil_inputs, sol_inputs)
             @test isfinite(dri.psio)
         end
@@ -443,7 +443,7 @@
         Eq = GeneralizedPerturbedEquilibrium.Equilibrium
 
         function build_solovev_equilibrium(; e=1.6, a=0.33, r0=1.0, q0=1.9,
-                mpsi=64, mtheta=128)
+            mpsi=64, mtheta=128)
             eq_config = Eq.EquilibriumConfig(;
                 eq_type="sol", eq_filename="unused",
                 jac_type="pest", grid_type="ldp",
@@ -454,7 +454,7 @@
         end
 
         @testset "Elongated Solovev (e=1.6)" begin
-            pe = build_solovev_equilibrium(e=1.6)
+            pe = build_solovev_equilibrium(; e=1.6)
             rsep, zsep, rext, zext = Eq.equilibrium_separatrix_find!(pe)
 
             # rsep[1] = outboard (R > R₀), rsep[2] = inboard (R < R₀)
@@ -465,8 +465,8 @@
             # rsep values consistent with r0=1.0, a=0.33 (Shafranov shift makes it approximate)
             amean = (rsep[1] - rsep[2]) / 2
             rmean = (rsep[1] + rsep[2]) / 2
-            @test amean ≈ 0.33 rtol=0.15
-            @test rmean ≈ 1.0 rtol=0.15
+            @test amean ≈ 0.33 rtol = 0.15
+            @test rmean ≈ 1.0 rtol = 0.15
 
             # rsep should be on the midplane (Z ≈ 0)
             # (verified indirectly: R at η=0 and η=0.5 are midplane by definition)
@@ -480,7 +480,7 @@
             @test zext ≈ zsep
 
             # Up-down symmetry of Solovev: |zsep_top| ≈ |zsep_bottom|
-            @test abs(zsep[1]) ≈ abs(zsep[2]) rtol=0.01
+            @test abs(zsep[1]) ≈ abs(zsep[2]) rtol = 0.01
 
             # Extremum R should be near the magnetic axis
             @test abs(rext[1] - pe.ro) < 0.2 * (rsep[1] - rsep[2])
@@ -489,11 +489,11 @@
             # kappa ≈ elongation
             kappa = (zsep[1] - zsep[2]) / (rsep[1] - rsep[2])
             @test kappa > 0
-            @test kappa ≈ 1.6 rtol=0.02
+            @test kappa ≈ 1.6 rtol = 0.02
         end
 
         @testset "Circular Solovev (e=1.0)" begin
-            pe = build_solovev_equilibrium(e=1.0)
+            pe = build_solovev_equilibrium(; e=1.0)
             rsep, zsep, rext, zext = Eq.equilibrium_separatrix_find!(pe)
 
             @test rsep[1] > pe.ro
@@ -505,15 +505,15 @@
             @test zsep[1] > zsep[2]
 
             # For circular cross-section, rext[1] ≈ rext[2] (top/bottom at same R)
-            @test rext[1] ≈ rext[2] rtol=0.01
+            @test rext[1] ≈ rext[2] rtol = 0.01
 
             kappa = (zsep[1] - zsep[2]) / (rsep[1] - rsep[2])
             @test kappa > 0
-            @test kappa ≈ 1.0 rtol=0.02
+            @test kappa ≈ 1.0 rtol = 0.02
         end
 
         @testset "Global scalars via equilibrium_global_parameters!" begin
-            pe = build_solovev_equilibrium(e=1.6)
+            pe = build_solovev_equilibrium(; e=1.6)
             Eq.equilibrium_global_parameters!(pe)
 
             # Separatrix convention: rsep[1]=outboard, rsep[2]=inboard,
@@ -522,24 +522,24 @@
             @test pe.params.zsep[1] > pe.params.zsep[2]
 
             # Shape parameters — all physically positive quantities.
-            @test pe.params.amean  > 0
-            @test pe.params.rmean  > 0
+            @test pe.params.amean > 0
+            @test pe.params.rmean > 0
             @test pe.params.aratio > 0
-            @test pe.params.kappa  > 0
-            @test pe.params.kappa  ≈ 1.6 rtol=0.02
+            @test pe.params.kappa > 0
+            @test pe.params.kappa ≈ 1.6 rtol = 0.02
 
             # For Solovev (e=1.6, a=0.33, r0=1.0) the shape is approximately
             # recovered (Shafranov shift loosens the match).
-            @test pe.params.amean ≈ 0.33 rtol=0.15
-            @test pe.params.rmean ≈ 1.0  rtol=0.15
+            @test pe.params.amean ≈ 0.33 rtol = 0.15
+            @test pe.params.rmean ≈ 1.0 rtol = 0.15
 
             # Consistency with separatrix formulae.
             @test pe.params.rmean ≈ (pe.params.rsep[1] + pe.params.rsep[2]) / 2
             @test pe.params.amean ≈ (pe.params.rsep[1] - pe.params.rsep[2]) / 2
 
             # Beta and field quantities — all physically positive.
-            @test pe.params.bt0   > 0
-            @test pe.params.crnt  > 0
+            @test pe.params.bt0 > 0
+            @test pe.params.crnt > 0
             @test pe.params.bwall > 0
             @test pe.params.betat > 0
             @test pe.params.betan > 0

From e875c3cbbd01adf32da23d37cc69f39e0facb030 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Sun, 4 Oct 2026 02:21:22 -0400
Subject: [PATCH 6/8] Equilibrium - BUGFIX - Error when an equilibrium ODE
 solve stops early, and raise the LAR step limit

With abstol following etol, the large-aspect-ratio profile solve needs ~2.3e4 Rosenbrock23 steps at the
default etol = 1e-10 and stopped at its 1e4-step limit at r = 0.42 a, so the equilibrium was built from a
truncated profile table. Every equilibrium ODE solve now checks its return code, and the LAR limit is 1e6.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/AnalyticEquilibrium.jl        | 11 +++++++----
 src/Equilibrium/DirectEquilibrium.jl          |  1 +
 src/Equilibrium/DirectEquilibriumArcLength.jl |  1 +
 src/Equilibrium/EquilibriumTypes.jl           | 12 ++++++++++++
 test/runtests_equil.jl                        | 12 ++++++++++++
 5 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl
index f0c924a1c..83842ca3b 100644
--- a/src/Equilibrium/AnalyticEquilibrium.jl
+++ b/src/Equilibrium/AnalyticEquilibrium.jl
@@ -140,7 +140,9 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
 
     prob = ODEProblem(dydr, y0, tspan, p)
 
-    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol), maxiters=10000, dense=false)
+    # The second-order method takes ~2e4 steps at the default etol = 1e-10, and about 3× more per decade below it.
+    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol), maxiters=10^6, dense=false)
+    check_equil_solve(sol, "large-aspect-ratio radial profiles")
 
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
@@ -377,11 +379,12 @@ function tj_analytic_shape_solve(p::TJAnalyticShapeParams, nu::Float64;
     saveat=nothing)
     rhs_params = (; p.a, p.B0, p.qc, p.mu, p.pc, p.epsa2, nu=nu)
     prob = ODEProblem(tj_analytic_shape_rhs!, tj_analytic_shape_initial(p, nu), (p.r0, p.a), rhs_params)
-    if saveat === nothing
-        return solve(prob, Vern9(); reltol, abstol, maxiters=10000, dense=false)
+    sol = if saveat === nothing
+        solve(prob, Vern9(); reltol, abstol, maxiters=10000, dense=false)
     else
-        return solve(prob, Vern9(); reltol, abstol, maxiters=10000, saveat=saveat)
+        solve(prob, Vern9(); reltol, abstol, maxiters=10000, saveat=saveat)
     end
+    return check_equil_solve(sol, "TJ-analytic shaping functions at ν = $nu")
 end
 
 """
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 266fa3cd5..7fad607b8 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -292,6 +292,7 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
 
     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
     sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=equil_abstol(equil_config.etol), dt=2π / 200, adaptive=true, dense=false)
+    check_equil_solve(sol, "field line at ψ_N = $psifac")
 
     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
     return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
index 1ed2d72bb..e7dd40772 100644
--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
+++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
@@ -115,6 +115,7 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     abstol_vec = [equil_abstol(equil_config.etol), equil_abstol(equil_config.etol), 1e20, 1e20, 1e20]
     sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec,
         dt=2π / 200, adaptive=true, dense=false)
+    check_equil_solve(sol, "arc-length field line at ψ_N = $psifac")
 
     n = length(sol.u)
     y_out = Matrix{Float64}(undef, n, 5)
diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index d5e7e5673..bc217c677 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -168,6 +168,18 @@ Absolute tolerance of an equilibrium ODE solve: follows `etol` but is never loos
 """
 equil_abstol(etol::Real, cap::Real=EQUIL_ABSTOL_MAX) = min(etol, cap)
 
+"""
+    check_equil_solve(sol, what)
+
+Return the equilibrium ODE solution `sol`, or raise an error naming the solve `what` if the integrator
+stopped early (step limit, step-size underflow, instability) instead of finishing or terminating on its callback.
+"""
+function check_equil_solve(sol, what::AbstractString)
+    OrdinaryDiffEq.SciMLBase.successful_retcode(sol) ||
+        error("Equilibrium ODE solve ($what) stopped early at t = $(sol.t[end]) with retcode $(sol.retcode).")
+    return sol
+end
+
 """
 Outer constructor for EquilibriumConfig from a parsed TOML dictionary
 """
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index c1fdebe20..cc1440eaf 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -136,6 +136,18 @@
         @test GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(; etol=1e-8).etol == 1e-8
     end
 
+    @testset "an equilibrium ODE solve that stops early is an error" begin
+        using OrdinaryDiffEq
+        check = GeneralizedPerturbedEquilibrium.Equilibrium.check_equil_solve
+        prob = ODEProblem((u, p, t) -> u, 1.0, (0.0, 10.0))
+        sol = solve(prob, Vern9())
+        @test check(sol, "test") === sol
+        stopped = @test_logs (:warn, r"maxiters") match_mode = :any solve(prob, Vern9(); maxiters=1)
+        @test_throws r"stopped early at t = .* with retcode MaxIters" check(stopped, "test")
+        terminated = solve(prob, Vern9(); callback=ContinuousCallback((u, t, integrator) -> t - 0.5, terminate!))
+        @test check(terminated, "test") === terminated
+    end
+
     @testset "Deprecated TOML keys are dropped, not fatal" begin
         # Removed control knobs must keep old gpec.toml decks (and older gpec.h5 replays,
         # whose stored TOML blob goes through the same path) parsing with a warning.

From 41762de6da1d2c665a3e3a2391aa3d123d881360 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Mon, 5 Oct 2026 20:03:05 -0400
Subject: [PATCH 7/8] Equilibrium - BUGFIX - Require the arc-length field line
 to close and drop two redundant solve guards

The arc-length solve passed the return-code check with Success, which there means the field line never returned
to the midplane; it now requires Terminated. The direct field-line check duplicated develop's stricter guard,
and the LAR step limit equalled the integrator's default.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/AnalyticEquilibrium.jl        | 3 +--
 src/Equilibrium/DirectEquilibrium.jl          | 1 -
 src/Equilibrium/DirectEquilibriumArcLength.jl | 7 ++++++-
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl
index 83842ca3b..633fe70a2 100644
--- a/src/Equilibrium/AnalyticEquilibrium.jl
+++ b/src/Equilibrium/AnalyticEquilibrium.jl
@@ -140,8 +140,7 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
 
     prob = ODEProblem(dydr, y0, tspan, p)
 
-    # The second-order method takes ~2e4 steps at the default etol = 1e-10, and about 3× more per decade below it.
-    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol), maxiters=10^6, dense=false)
+    sol = solve(prob, Rosenbrock23(; autodiff=false); reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol), dense=false)
     check_equil_solve(sol, "large-aspect-ratio radial profiles")
 
     r_arr = sol.t
diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 5b98a31de..f68f949ca 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -292,7 +292,6 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
 
     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
     sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=equil_abstol(equil_config.etol), dt=2π / 200, adaptive=true, dense=false)
-    check_equil_solve(sol, "field line at ψ_N = $psifac")
 
     # A failed solve returns a truncated solution instead of throwing; check both the retcode and
     # that the field line reached η = 2π, since a callback can end it early and still report Success.
diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
index e7dd40772..1ddbdd469 100644
--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
+++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
@@ -115,7 +115,12 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     abstol_vec = [equil_abstol(equil_config.etol), equil_abstol(equil_config.etol), 1e20, 1e20, 1e20]
     sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec,
         dt=2π / 200, adaptive=true, dense=false)
-    check_equil_solve(sol, "arc-length field line at ψ_N = $psifac")
+    # The only complete ending is the callback's terminate! on return to the midplane; Success means it never fired.
+    sol.retcode == ReturnCode.Terminated || error(
+        "arclength_fieldline_int: field line at psifac = $(@sprintf("%.6f", psifac)) did not return to the midplane " *
+        "(retcode $(sol.retcode) at arc length $(sol.t[end])); the flux surface did not close. " *
+        "This usually means psihigh is too close to the separatrix for the equilibrium grid to resolve."
+    )
 
     n = length(sol.u)
     y_out = Matrix{Float64}(undef, n, 5)

From 6d6a4932dd34129d45e0d5c61d276fbc8756ba75 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 9 Oct 2026 12:42:59 -0400
Subject: [PATCH 8/8] =?UTF-8?q?Equilibrium=20-=20BUGFIX=20-=20Let=20a=20fa?=
 =?UTF-8?q?iled=20shaping=20solve=20out=20of=20the=20=CE=BD=20root-find=20?=
 =?UTF-8?q?and=20accept=20only=20finished=20return=20codes?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Review fixes. The ν root-find falls back to ν = qa/qc only when the
root-find itself has no bracket or does not converge. check_equil_solve
accepts Success and Terminated only. The TJ shaping solve drops its step
cap, redundant abstol arguments and a single-use constant go, and the two
solve helpers move to SolveTolerances.jl.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/AnalyticEquilibrium.jl | 25 ++++++++++++-----------
 src/Equilibrium/Equilibrium.jl         |  1 +
 src/Equilibrium/EquilibriumTypes.jl    | 28 ++------------------------
 src/Equilibrium/SolveTolerances.jl     | 28 ++++++++++++++++++++++++++
 test/runtests_equil.jl                 |  8 ++++++++
 5 files changed, 52 insertions(+), 38 deletions(-)
 create mode 100644 src/Equilibrium/SolveTolerances.jl

diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl
index 633fe70a2..eb499dba1 100644
--- a/src/Equilibrium/AnalyticEquilibrium.jl
+++ b/src/Equilibrium/AnalyticEquilibrium.jl
@@ -379,9 +379,9 @@ function tj_analytic_shape_solve(p::TJAnalyticShapeParams, nu::Float64;
     rhs_params = (; p.a, p.B0, p.qc, p.mu, p.pc, p.epsa2, nu=nu)
     prob = ODEProblem(tj_analytic_shape_rhs!, tj_analytic_shape_initial(p, nu), (p.r0, p.a), rhs_params)
     sol = if saveat === nothing
-        solve(prob, Vern9(); reltol, abstol, maxiters=10000, dense=false)
+        solve(prob, Vern9(); reltol, abstol, dense=false)
     else
-        solve(prob, Vern9(); reltol, abstol, maxiters=10000, saveat=saveat)
+        solve(prob, Vern9(); reltol, abstol, saveat=saveat)
     end
     return check_equil_solve(sol, "TJ-analytic shaping functions at ν = $nu")
 end
@@ -394,7 +394,8 @@ https://github.com/rfitzp/TJ): solve for ν so that q₂(x=1) matches
 `q₂ = x²·(1+εa²·g₂)·exp(−εa²·f3/f1)/f1`; at x=1 and low β this picks up an
 O(εa²) correction relative to the lowest-order guess ν = qa/qc, which
 matters for the TJ-analytic benchmark at large ε.  Falls back to the
-lowest-order ν if the bracket search diverges.
+lowest-order ν if the root-find has no bracket or does not converge; an
+error from a trial shaping solve propagates.
 """
 function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; reltol::Float64=1e-7, abstol::Float64=equil_abstol(reltol))
     function q2_edge(nu::Float64)
@@ -409,7 +410,10 @@ function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; relto
         find_zero(nu -> q2_edge(nu) - qa_target, (0.5 * nu_guess, 2 * nu_guess);
             atol=1e-8, rtol=1e-10)
     catch err
-        @warn "ν root-find failed for TJ-analytic equilibrium; falling back to lowest-order ν = qa/qc" error = err
+        # Only a failed root-find falls back; a shaping solve that stopped early must propagate.
+        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
 end
@@ -458,8 +462,8 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig)
     epsa2 = p.epsa2
     p00_phys = B0^2 * epsa2 * pc          # μ₀P = B₀²·εa²·p₂ at axis
 
-    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol))
-    sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol))
+    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol)
+    sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol)
 
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
@@ -552,9 +556,6 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig)
     return InverseRunInput(equil_input, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, R0, 0.0, psio, nothing)
 end
 
-# Tighter abstol ceiling for the dense-saveat TJ-analytic solve that feeds the (R, Z) → (r, w) Newton inversion.
-const TJ_DENSE_ABSTOL_MAX = 1e-10
-
 """
     tj_analytic_run_direct(equil_input, tj_input; nrbox=257, nzbox=257, rc=1.2)
 
@@ -598,14 +599,14 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo
     p00_phys = B0^2 * epsa2 * pc
 
     # ν root-find (cf. Fitzpatrick TJ's Setnu): q₂(1) = qa_target.
-    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol))
+    nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol)
 
     # Dense saveat so the downstream splines (H₁, g₂, f₃, ψ) are evaluated on
     # a fine uniform r grid rather than the ~30 adaptive Vern9 steps — otherwise
     # the (R, Z) → (r, w) Newton iteration hits spline interpolation artifacts.
+    # That inversion also needs a tighter abstol ceiling (1e-10) than the other equilibrium solves.
     dense_r = collect(range(p.r0, p.a; length=1024))
-    sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol,
-        abstol=equil_abstol(equil_input.etol, TJ_DENSE_ABSTOL_MAX), saveat=dense_r)
+    sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol, abstol=equil_abstol(equil_input.etol, 1e-10), saveat=dense_r)
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
 
diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl
index 9022cbd37..9bb0bd3ae 100644
--- a/src/Equilibrium/Equilibrium.jl
+++ b/src/Equilibrium/Equilibrium.jl
@@ -12,6 +12,7 @@ import ..Utilities
 
 # --- Internal Module Structure ---
 include("EquilibriumTypes.jl")
+include("SolveTolerances.jl")
 include("GridRefinement.jl")
 include("FluxSurfaceMetrics.jl")
 include("CoordinateInvariant.jl")
diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index bc217c677..53f2ab8f2 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -38,8 +38,8 @@ specified in the input.
   - `mtheta::Int` - Number of poloidal grid points
   - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to
     meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped)
-  - `etol::Float64` - Relative tolerance of the equilibrium ODE solves (positive and finite); their
-    absolute tolerance is `equil_abstol(etol)`, i.e. `etol` capped at `EQUIL_ABSTOL_MAX`
+  - `etol::Float64` - Relative tolerance of the equilibrium ODE solves; their absolute tolerance
+    is `equil_abstol(etol)`, i.e. `etol` capped at `EQUIL_ABSTOL_MAX`
   - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations)
   - `use_galgrid::Bool` - Use the same grid as galerkin method
 """
@@ -156,30 +156,6 @@ specified in the input.
     end
 end
 
-"""
-Ceiling on the absolute tolerance of the equilibrium ODE solves; the historical fixed value, kept so a loose `etol` never loosens abstol.
-"""
-const EQUIL_ABSTOL_MAX = 1e-8
-
-"""
-    equil_abstol(etol, cap=EQUIL_ABSTOL_MAX)
-
-Absolute tolerance of an equilibrium ODE solve: follows `etol` but is never looser than `cap`.
-"""
-equil_abstol(etol::Real, cap::Real=EQUIL_ABSTOL_MAX) = min(etol, cap)
-
-"""
-    check_equil_solve(sol, what)
-
-Return the equilibrium ODE solution `sol`, or raise an error naming the solve `what` if the integrator
-stopped early (step limit, step-size underflow, instability) instead of finishing or terminating on its callback.
-"""
-function check_equil_solve(sol, what::AbstractString)
-    OrdinaryDiffEq.SciMLBase.successful_retcode(sol) ||
-        error("Equilibrium ODE solve ($what) stopped early at t = $(sol.t[end]) with retcode $(sol.retcode).")
-    return sol
-end
-
 """
 Outer constructor for EquilibriumConfig from a parsed TOML dictionary
 """
diff --git a/src/Equilibrium/SolveTolerances.jl b/src/Equilibrium/SolveTolerances.jl
new file mode 100644
index 000000000..097940e65
--- /dev/null
+++ b/src/Equilibrium/SolveTolerances.jl
@@ -0,0 +1,28 @@
+# SolveTolerances.jl
+#
+# Absolute tolerance and completion check shared by the equilibrium ODE solves
+# (direct field-line, large-aspect-ratio and TJ-analytic integrations).
+
+"""
+Ceiling on the absolute tolerance of the equilibrium ODE solves; the historical fixed value, kept so a loose `etol` never loosens abstol.
+"""
+const EQUIL_ABSTOL_MAX = 1e-8
+
+"""
+    equil_abstol(etol, cap=EQUIL_ABSTOL_MAX)
+
+Absolute tolerance of an equilibrium ODE solve: follows `etol` but is never looser than `cap`.
+"""
+equil_abstol(etol::Real, cap::Real=EQUIL_ABSTOL_MAX) = min(etol, cap)
+
+"""
+    check_equil_solve(sol, what)
+
+Return the equilibrium ODE solution `sol`, or raise an error naming the solve `what` if the integrator
+stopped early (step limit, step-size underflow, instability) instead of finishing or terminating on its callback.
+"""
+function check_equil_solve(sol, what::AbstractString)
+    sol.retcode == ReturnCode.Success || sol.retcode == ReturnCode.Terminated ||
+        error("Equilibrium ODE solve ($what) stopped early at t = $(sol.t[end]) with retcode $(sol.retcode).")
+    return sol
+end
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 53b745e6a..2d33380c6 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -153,6 +153,14 @@
         @test_throws r"stopped early at t = .* with retcode MaxIters" check(stopped, "test")
         terminated = solve(prob, Vern9(); callback=ContinuousCallback((u, t, integrator) -> t - 0.5, terminate!))
         @test check(terminated, "test") === terminated
+        # Only Success and Terminated count as finished; other "successful" codes are nonlinear-solver exits.
+        stalled = OrdinaryDiffEq.SciMLBase.solution_new_retcode(sol, ReturnCode.StalledSuccess)
+        @test_throws r"retcode StalledSuccess" check(stalled, "test")
+
+        # A trial shaping solve that stops early must come out of the ν root-find, not become its ν = qa/qc fallback.
+        Equil = GeneralizedPerturbedEquilibrium.Equilibrium
+        tj = Equil.TJAnalyticConfig(; lar_r0=4.0, lar_a=1.0, qc=1.5, qa=3.6, pc=0.001, mu=2.0, B0=12.0, ma=64, mtau=64)
+        @test_throws r"stopped early at t = .* with retcode DtNaN" Equil.tj_analytic_find_nu(Equil.TJAnalyticShapeParams(tj), tj.qa; abstol=NaN)
     end
 
     @testset "Deprecated TOML keys are dropped, not fatal" begin