Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 20 additions & 0 deletions src/Equilibrium/DirectEquilibrium.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 12 additions & 0 deletions test/runtests_equil.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading