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

src/Equilibrium/AnalyticEquilibrium.jl

@@ -140,7 +140,8 @@ 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=equil_abstol(equil_input.etol), dense=false)
+    check_equil_solve(sol, "large-aspect-ratio radial profiles")
 
     r_arr = sol.t
     y_mat = reduce(hcat, sol.u)'
@@ -369,18 +370,20 @@ 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=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)
-    if saveat === nothing
-        return solve(prob, Vern9(); reltol, abstol, maxiters=10000, dense=false)
+    sol = if saveat === nothing
+        solve(prob, Vern9(); reltol, abstol, dense=false)
     else
-        return 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
 
 """
@@ -391,11 +394,12 @@ 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)
+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)
@@ -406,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
@@ -597,9 +604,9 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo
     # 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=1e-10, 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)'
 

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=equil_abstol(equil_config.etol), dt=2π / 200, adaptive=true, dense=false)
 
     # 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.

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.
@@ -112,9 +112,15 @@ 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 = [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)
+    # 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)
@@ -144,4 +150,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
-

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")

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
 """
@@ -143,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

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

test/runtests_equil.jl

@@ -136,6 +136,33 @@
         @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 "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
+        # 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
         # 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.
@@ -150,7 +177,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 +276,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 +292,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 +365,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 +394,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 +416,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 +445,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 +470,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 +481,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 +492,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 +507,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 +516,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 +532,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 +549,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