EQUIL - BUGFIX - Detect field-line ODE failure; add psiedge floor and guard the trial loop

Section 0 of the plan first: triple-checked the metres->flux conversion before
building on it. All three checks pass on the DIII-D-like deck.
  - d_i = 0.0509 m vs the NRL c/omega_pi form 0.0510 m (0.2%)
  - |dels_db| = 0.67..1.45, i.e. genuinely O(1); delta_s = 1.8..4.4 mm;
    delta_s/r_s = 0.3..1.1%
  - da/dpsi agrees to <=3.6% across surface_da_dpsi, a different-h difference,
    and the secant between adjacent surfaces
  - the set of overlapping pairs is IDENTICAL computed all-in-flux and
    all-in-metres, so the conversion cannot invent an overlap
Decisive on its own, needing no conversion: q=6 and q=7 sit 1.821 mm apart with
half-widths summing to 2.338 mm. Every other adjacent pair is 8.7 mm or more
apart. The truncation removes 3.59 mm, 0.66% of the minor radius -- far out, as
expected; the psi number reads deep only because psi is compressed at the edge.
The original mixed-units comparison finds NO overlap at all on the in-grid
surfaces, so it would not have truncated.

DirectEquilibrium.jl: check sol.retcode and the eta=2*pi endpoint after the
field-line solve. A failed integration does not throw -- it returns a solution
truncated wherever it gave up, and the caller read the last point as if the
surface closed. That was a silently wrong equilibrium; it is now a named error
that says which psifac failed.

Lower gate: the resolved psihigh may not fall below the top of the rational-q
window psiedge sits in -- the first rational surface OUTBOARD of psiedge --
because the dW edge scan runs over [psiedge, psilim] and needs room for its
maximum. psiedge >= psihigh (the documented scan-disabled sentinel) means no
floor, else decks that switch the scan off could never truncate. Independently,
refuse to truncate if fewer than two rational surfaces would survive.

On the SLAYER deck the floor binds: the criterion asks 0.98919, psiedge=0.99
sits in the q=6 window, so the domain is held at 0.99328 and the warning says
the retained surfaces are set by psiedge rather than by physics.

Trial loop: catch any exception from setup_equilibrium (axis finder, separatrix
finder, field-line ODE, Roots convergence) as a failed trial and back off, and
restore eq_config.psihigh, which was mutated before the attempt. Also thread
DirectRunInput.psi_in through so the round-trip tracing residual runs on plain
efit decks, not only efit_by_inversion -- the capability existed, it was unwired.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PSrf6JCViFfVzqzkQ66o6b

Author: logan-nc

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; callers trialling a large psihigh catch it and back off.
+    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

src/GeneralizedPerturbedEquilibrium.jl

@@ -120,12 +120,25 @@ function _layer_overlap_profiles(inputs, dir_path::AbstractString)
     return nothing
 end
 
+# Lower gate. The resolved psihigh may not fall below the top of the rational-q window that
+# `psiedge` sits in -- i.e. the first rational surface outboard of psiedge -- because the dW edge
+# scan runs over [psiedge, psilim] and needs room for its maximum to be findable. `min` with the
+# deck value keeps the floor from ever raising psihigh above what was asked for. Returns nothing
+# when psiedge >= psihigh (the documented "scan disabled" sentinel: no dW band to protect) or
+# when no rational surface lies outboard of psiedge.
+function _psiedge_floor(psis, psiedge::Real, psihigh_deck::Real)
+    psiedge >= psihigh_deck && return nothing
+    idx = findfirst(>(psiedge), psis)
+    idx === nothing && return nothing
+    return (floor=min(psihigh_deck, psis[idx]), index=idx)
+end
+
 # Run the overlap scan on the pass-1 equilibrium, gate the result against what the
 # equilibrium reconstruction can actually deliver, and set `eq_config.psihigh`.
 # Returns the equilibrium to carry forward: the pass-1 one for an inward move (it spans a
 # superset, so its measured knot density is still usable), or the accepted trial for an
 # outward move. `eq_config.psihigh` is left untouched when the scan cannot run.
-function _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, intr;
+function _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, intr, ctrl;
     max_trials::Int=4)
     if intr.nlow != intr.nhigh
         @warn "psihigh_from_layer_overlap is ignored for multi-n runs (nlow=$(intr.nlow), nhigh=$(intr.nhigh)); " *
@@ -153,21 +166,66 @@ function _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, in
     k = scan.first_overlap
     why = "layers of the q=$(scan.m[k])/$(scan.n[k]) and q=$(scan.m[k-1])/$(scan.n[k-1]) surfaces overlap"
 
+    # Refuse outright if the domain would keep fewer than two rational surfaces: a one-surface
+    # domain cannot support the analysis this feature serves, and widths asking for it are far
+    # more likely wrong than the plasma is that pathological.
+    n_kept = count(<=(target), scan.psi)
+    if n_kept < 2
+        @warn "psihigh_from_layer_overlap: the criterion asks for psihigh = $(@sprintf("%.6f", target)), " *
+              "which would keep only $n_kept rational surface(s) ($why). Refusing to truncate — " *
+              "check the layer widths. Keeping psihigh = $(@sprintf("%.6f", psihigh_safe))."
+        return equil
+    end
+
+    # Lower gate: never cut inside the rational window psiedge sits in.
+    fl = _psiedge_floor(scan.psi, ctrl.psiedge, psihigh_safe)
+    if fl !== nothing && target < fl.floor
+        eq_config.psihigh = fl.floor
+        @warn "psihigh_from_layer_overlap: the layer criterion asks for psihigh = " *
+              "$(@sprintf("%.6f", target)) ($why), but psiedge = $(@sprintf("%.6f", ctrl.psiedge)) " *
+              "sits inside the q=$(scan.m[fl.index])/$(scan.n[fl.index]) window, so the domain is " *
+              "held at $(@sprintf("%.6f", fl.floor)) to keep the dW edge scan band open. The layer " *
+              "criterion is overridden here — the retained surfaces are set by psiedge, not by physics."
+        return equil
+    end
+
     # Inward: the pass-1 equilibrium already covers the domain, so nothing to verify.
     if target <= psihigh_safe
         eq_config.psihigh = target
-        @info "psihigh_from_layer_overlap: $(@sprintf("%.6f", psihigh_safe)) -> $(@sprintf("%.6f", target)) ($why)"
+        a_top = Tearing.InnerLayer.SLAYER.surface_minor_radius(equil, psihigh_safe)
+        a_cut = Tearing.InnerLayer.SLAYER.surface_minor_radius(equil, target)
+        @info "psihigh_from_layer_overlap: $(@sprintf("%.6f", psihigh_safe)) -> $(@sprintf("%.6f", target)) " *
+              "($why); removes $(@sprintf("%.2f", 1e3 * (a_top - a_cut))) mm of minor radius, " *
+              "$(@sprintf("%.2f", 100 * (a_top - a_cut) / a_top))% of a, keeping $n_kept of $(length(scan.psi)) surfaces"
         return equil
     end
 
-    # Outward: the domain has to be formed and checked before it can be trusted.
+    # Outward: the domain has to be formed and checked before it can be trusted. The raw psi(R,Z)
+    # map lives on DirectRunInput, which plain `efit` uses too -- pass it so the round-trip
+    # tracing residual runs on every direct/EFIT deck, not only on efit_by_inversion.
+    psi_map = rerun_input isa Equilibrium.DirectRunInput ? rerun_input.psi_in : nothing
+    psio_ref = rerun_input isa Equilibrium.DirectRunInput ? rerun_input.psio : nothing
     accepted, accepted_equil = psihigh_safe, equil
     lo, hi = psihigh_safe, target
     for trial in 1:max_trials
         eq_config.psihigh = hi
-        trial_equil = Equilibrium.setup_equilibrium(eq_config, rerun_input)
+        # Forming an equilibrium this far out can fail hard: the axis and separatrix finders,
+        # the field-line ODE, and the Roots calls inside it all throw near the separatrix. Any
+        # of those is a failed trial, not a failed run -- and psihigh must be restored, since
+        # it was mutated before the attempt.
+        trial_equil = try
+            Equilibrium.setup_equilibrium(eq_config, rerun_input)
+        catch err
+            err isa InterruptException && rethrow()
+            eq_config.psihigh = accepted
+            @warn "psihigh_from_layer_overlap: equilibrium formation threw at psihigh = " *
+                  "$(@sprintf("%.6f", hi)) (trial $trial/$max_trials): $(sprint(showerror, err)); " *
+                  "backing off toward $(@sprintf("%.6f", lo))."
+            hi = (lo + hi) / 2
+            continue
+        end
         hi = eq_config.psihigh   # setup_equilibrium may clamp back to the closed-flux region
-        q = Equilibrium.equilibrium_edge_quality(trial_equil)
+        q = Equilibrium.equilibrium_edge_quality(trial_equil; psi_map=psi_map, psio=psio_ref)
         if q.ok
             accepted, accepted_equil = hi, trial_equil
             @info "psihigh_from_layer_overlap: trial at psihigh = $(@sprintf("%.6f", hi)) accepted " *
@@ -392,7 +450,7 @@ function main_from_inputs(
     # moves after the nodes are built.
     psihigh_pass1 = eq_config.psihigh
     if eq_config.psihigh_from_layer_overlap
-        equil = _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, intr)
+        equil = _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, intr, ctrl)
     end
     psihigh_moved = eq_config.psihigh != psihigh_pass1
 

test/runtests_layer_overlap.jl

@@ -118,6 +118,43 @@
         @test q.ringing_ok
     end
 
+    @testset "psiedge floor: window top, not bottom" begin
+        # Surfaces bracketing psiedge=0.99: the floor must be the first surface OUTBOARD of it,
+        # so the dW scan band [psiedge, psilim] stays open.
+        psis = [0.52, 0.77, 0.89, 0.968, 0.9933, 0.9985]
+        fl = GeneralizedPerturbedEquilibrium._psiedge_floor(psis, 0.99, 0.9995)
+        @test fl !== nothing
+        @test fl.floor ≈ 0.9933          # q=6, the window top -- NOT 0.968
+        @test fl.index == 5
+        # min() with the deck value: the floor can never raise psihigh above what was asked for.
+        fl2 = GeneralizedPerturbedEquilibrium._psiedge_floor(psis, 0.99, 0.992)
+        @test fl2.floor ≈ 0.992
+    end
+
+    @testset "psiedge floor: disabled sentinel and no-surface cases" begin
+        psis = [0.52, 0.77, 0.89, 0.968, 0.9933, 0.9985]
+        # psiedge >= psihigh is the documented "scan disabled" sentinel -> no floor at all,
+        # otherwise decks that switch the dW scan off could never truncate.
+        @test GeneralizedPerturbedEquilibrium._psiedge_floor(psis, 1.0, 0.9995) === nothing
+        @test GeneralizedPerturbedEquilibrium._psiedge_floor(psis, 0.9995, 0.9995) === nothing
+        # No rational surface outboard of psiedge -> nothing to floor against.
+        @test GeneralizedPerturbedEquilibrium._psiedge_floor(psis, 0.999, 0.9995) === nothing
+        @test GeneralizedPerturbedEquilibrium._psiedge_floor(Float64[], 0.99, 0.9995) === nothing
+    end
+
+    @testset "field-line ODE failure is an error, not silent truncation" begin
+        # A failed solve returns a solution truncated short of eta = 2*pi, which used to be
+        # consumed as a closed flux surface. Both guards must be present in direct_fieldline_int:
+        # the retcode check and the endpoint check (a solve can stop early and still report
+        # Success via a callback-driven terminate).
+        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)
+        # Both must name the surface so a failure is attributable to a psihigh, not just a stack trace.
+        @test count("direct_fieldline_int: field-line integration", src) == 2
+        @test occursin("psihigh is too close to the separatrix", src)
+    end
+
     @testset "psihigh_from_layer_overlap defaults off" begin
         @test Equilibrium.EquilibriumConfig().psihigh_from_layer_overlap == false
         @test eq_cfg.psihigh_from_layer_overlap == false