From d65375dffb4d61653151b793126fb2790f710672 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 4 Sep 2026 13:17:57 -0400 Subject: [PATCH 1/2] Equilibrium - 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 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 From e13de9c35a7a86fd7618aefd4049b156f6b542d0 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 1 Oct 2026 15:33:39 -0400 Subject: [PATCH 2/2] Equilibrium - MINOR - Trim the field-line guard comments to two lines Both annotation blocks narrated the investigation behind the guard. Keep only the non-obvious reason two checks are needed: a callback can end the solve early and still report Success. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PSrf6JCViFfVzqzkQ66o6b --- src/Equilibrium/DirectEquilibrium.jl | 6 ++---- test/runtests_equil.jl | 8 ++------ 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index cf775e600..3ade63787 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -293,10 +293,8 @@ 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. + # 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. if sol.retcode != ReturnCode.Success error( "direct_fieldline_int: field-line integration failed at psifac = " * diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 2e1014a04..62c9719cf 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -554,12 +554,8 @@ 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. + # Checked in the source, since a non-closing field line is costly to synthesize; the EFIT + # testsets above are the 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)