From 0bb4993a5d6e3efe73db63adab8ae97a3e7c556e Mon Sep 17 00:00:00 2001
From: logan-nc <nikolas.logan@columbia.edu>
Date: Fri, 4 Sep 2026 13:17:57 -0400
Subject: [PATCH] EQUIL - BUGFIX - Detect field-line ODE failure instead of
 consuming a truncated surface

direct_fieldline_int solved the flux-surface field-line ODE and never checked
sol.retcode, then read sol.u[end] as though the surface had closed at eta = 2*pi.
A failed solve does not throw -- it returns a solution truncated wherever the
integrator gave up -- so the caller silently built a wrong flux surface and the
damage surfaced downstream as unrelated nonsense, or not at all.

Check both the retcode and the endpoint: a solve can stop early through a
callback-driven terminate and still report Success, so neither check subsumes the
other. Both errors name the psifac that failed and point at psihigh sitting too
close to the separatrix, which is the usual cause.

No healthy deck is affected: runtests_equil.jl passes 281/281 and the DIII-D-like
ideal example runs end to end with neither guard firing.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PSrf6JCViFfVzqzkQ66o6b
---
 src/Equilibrium/DirectEquilibrium.jl | 20 ++++++++++++++++++++
 test/runtests_equil.jl               | 12 ++++++++++++
 2 files changed, 32 insertions(+)

diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 961a0f7f5..265dbbfb4 100644
--- a/src/Equilibrium/DirectEquilibrium.jl
+++ b/src/Equilibrium/DirectEquilibrium.jl
@@ -293,6 +293,26 @@ 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=1e-8, dt=2π / 200, adaptive=true, dense=false)
 
+    # A failed integration does not throw: it returns a solution truncated wherever it gave up,
+    # and the caller reads the last point as if it closed the field line at η = 2π. That is a
+    # silently wrong surface, so check the retcode and the endpoint. Surfaces very close to the
+    # separatrix are where this fires.
+    if sol.retcode != ReturnCode.Success
+        error(
+            "direct_fieldline_int: field-line integration failed at psifac = " *
+            "$(@sprintf("%.6f", psifac)) (retcode $(sol.retcode)); the flux surface did not " *
+            "close. This usually means psihigh is too close to the separatrix for the " *
+            "equilibrium grid to resolve."
+        )
+    end
+    if !isapprox(sol.t[end], 2π; atol=1e-8)
+        error(
+            "direct_fieldline_int: field-line integration at psifac = " *
+            "$(@sprintf("%.6f", psifac)) stopped at eta = $(@sprintf("%.6f", sol.t[end])) " *
+            "instead of 2*pi; the flux surface did not close."
+        )
+    end
+
     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
     return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
 end
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index af3cd5743..2e1014a04 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -552,4 +552,16 @@
             @test pe.params.li3 > 0
         end
     end
+
+    @testset "field-line ODE failure is an error, not silent truncation" begin
+        # A failed solve does not throw -- it returns a solution truncated wherever it gave up,
+        # which direct_fieldline_int used to consume as a closed flux surface. Both guards must
+        # be present: the retcode check, and the endpoint check (a solve can stop early via a
+        # callback-driven terminate and still report Success). Asserted against the source
+        # because a genuinely non-closing field line cannot be synthesized cheaply here; the
+        # EFIT testsets above are the standing positive control that neither guard false-fires.
+        src = read(joinpath(dirname(@__DIR__), "src", "Equilibrium", "DirectEquilibrium.jl"), String)
+        @test occursin("sol.retcode != ReturnCode.Success", src)
+        @test occursin("isapprox(sol.t[end], 2π", src)
+    end
 end