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