Equilibrium - BUGFIX - Detect field-line ODE failure instead of consuming a truncated surface

[logan-nc]

Release note

• Audience: users
• Numerical impact: none (harness @ d65375d)
• Migration: none

A flux surface whose field-line integration failed used to be consumed as if it had closed, so GPEC
carried on and produced a silently wrong equilibrium. It now stops with an error naming the psifac
that failed. This most often means psihigh is too close to the separatrix for the equilibrium grid
to resolve.

Regression report

Run at d65375dff, before the 2026-10-01 develop merge. Not re-run after it: the only later changes on this branch are that merge and comment-only edits, so the stamp is deliberately left where the harness actually ran.

Regression Report: diiid_n1
==============================================================================================================
Ref 1: develop  @ f36ecef86 (2026-09-11)
       env: julia 1.11.7, x86_64-linux-gnu, manifest e0a28574 (pinned), 16 threads/16 BLAS
Ref 2: bugfix/fieldline-ode-retcode  @ d65375dff (2026-09-04)
       env: julia 1.11.7, x86_64-linux-gnu, manifest e0a28574 (pinned), 16 threads/16 BLAS
--------------------------------------------------------------------------------------------------------------
Quantity                                      develop          bugfix/fieldline-ode-retcode  Diff       Status
--------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.012318e-01     8.012318e-01                  0.0e+00    OK    
total energy Im(et[1])                        4.142529e-05     4.142529e-05                  0.0e+00    OK    
plasma energy Re(ep[1])                       -1.348486e+00    -1.348486e+00                 0.0e+00    OK    
vacuum energy Re(ev[1])                       2.149718e+00     2.149718e+00                  0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873976e-01     1.873976e-01                  0.0e+00    OK    
plasma energy (all)                           [35 elem]        [35 elem]                     0.0e+00    OK    
vacuum energy (all)                           [35 elem]        [35 elem]                     0.0e+00    OK    
total energy (all)                            [35 elem]        [35 elem]                     0.0e+00    OK    
ODE steps (saved)                             2576             2576                          0.0e+00    OK    
ODE steps (total)                             4572             4572                          0.0e+00    OK    
q0                                            1.204212e+00     1.204212e+00                  0.0e+00    OK    
q95                                           4.781723e+00     4.781723e+00                  0.0e+00    OK    
beta_t                                        1.327024e-02     1.327024e-02                  0.0e+00    OK    
beta_n                                        1.372511e+00     1.372511e+00                  0.0e+00    OK    
internal inductance li1                       8.842392e-01     8.842392e-01                  0.0e+00    OK    
internal inductance li2                       7.080847e-01     7.080847e-01                  0.0e+00    OK    
internal inductance li3                       7.304433e-01     7.304433e-01                  0.0e+00    OK    
poloidal beta betap1                          6.680744e-01     6.680744e-01                  0.0e+00    OK    
poloidal beta betap2                          5.349834e-01     5.349834e-01                  0.0e+00    OK    
poloidal beta betap3                          5.518761e-01     5.518761e-01                  0.0e+00    OK    
# singular surfaces                           5                5                             0.0e+00    OK    
singular psi locations                        [5 elem]         [5 elem]                      0.0e+00    OK    
singular q values                             [5 elem]         [5 elem]                      0.0e+00    OK    
current beta betaj                            4.236478e-01     4.236478e-01                  0.0e+00    OK    
plasma volume                                 1.829472e+01     1.829472e+01                  0.0e+00    OK    
plasma current                                1.152130e+00     1.152130e+00                  0.0e+00    OK    
mpert                                         35               35                            0.0e+00    OK    
npert                                         1                1                             0.0e+00    OK    
toroidal field bt0                            2.006573e+00     2.006573e+00                  0.0e+00    OK    
wall field bwall                              3.880145e-01     3.880145e-01                  0.0e+00    OK    
aspect ratio                                  2.845746e+00     2.845746e+00                  0.0e+00    OK    
elongation kappa                              1.708350e+00     1.708350e+00                  0.0e+00    OK    
q profile (checksum)                          0cd285cea88d...  0cd285cea88d...               identical  OK    
pressure profile (checksum)                   a1c48b266622...  a1c48b266622...               identical  OK    
Mercier D_I profile (checksum)                eeb06744e795...  eeb06744e795...               identical  OK    
resistive interchange D_R profile (checksum)  fa37296851f6...  fa37296851f6...               identical  OK    
ballooning Delta' profile (checksum)          bb713caf14eb...  bb713caf14eb...               identical  OK    
island half-widths                            [5 elem]         [5 elem]                      0.0e+00    OK    
Chirikov parameter                            [5 elem]         [5 elem]                      0.0e+00    OK    
||resonant area-weighted field||              5.189179e-04     5.189179e-04                  0.0e+00    OK    
PE plasma energy                              3.422677e+00     3.422677e+00                  0.0e+00    OK    
PE vacuum energy                              3.174510e+00     3.174510e+00                  0.0e+00    OK    
PE surface energy                             5.826684e+00     5.826684e+00                  0.0e+00    OK    
PE toroidal torque                            5.087465e-02     5.087465e-02                  0.0e+00    OK    
NTV torque FGAR [N·m]                         5.296762e-01     5.296762e-01                  0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                7.924971e-02     7.924971e-02                  0.0e+00    OK    
Runtime (s)                                   302.9s           280.9s                                   --    
resonant area-weighted field b^r              [5 elem]         [5 elem]                      0.0e+00    OK    
==============================================================================================================
Summary: 47 unchanged

The bug

direct_fieldline_int solves the flux-surface field-line ODE and never checks sol.retcode. The
very next line consumes sol.u[end] as though the field line had closed at eta = 2*pi. But a
failed solve does not throw — it returns a solution truncated wherever the integrator gave up. So
the caller silently builds a flux surface from a partial trace, and the damage surfaces downstream as
unrelated nonsense, or not at all.

The fix

Two checks, because neither subsumes the other:

• sol.retcode != ReturnCode.Success — the ordinary failure.
• sol.t[end] did not reach 2*pi — a solve can stop early through a callback-driven terminate and
  still report Success. direct_fieldline_int installs a DiscreteCallback, so this is reachable.

Both errors name the psifac that failed, and the first points at psihigh sitting too close to the
separatrix, which is the usual cause. The failure becomes attributable to a surface rather than a
stack trace somewhere further down the pipeline.

Verification

• Pkg.test() sandbox, runtests_equil.jl — 286/286 pass (the Equilibrium Unit Tests set;
  281 of those predate this PR, the rest come from sibling work on this base).
• DIII-D-like ideal example runs end to end (equilibrium → FFS → PE → KineticForces), neither guard
  firing.
• Regression harness table above: 47 quantities, all unchanged, rebased onto and compared against
  develop @ f36ecef86. Expected — on a healthy equilibrium this adds two predicates and changes
  nothing.

A note on the test

The added testset asserts the guards against the source text rather than triggering them. A genuinely
non-closing field line is not cheap to synthesize here, and forcing psifac toward the separatrix
would likely hit a different error first and be flaky across machines. The standing positive control
that neither guard false-fires is the EFIT testsets already in that file, which trace real surfaces
through this function. Happy to swap in a behavioral trigger if a reviewer sees a clean way to build one.

Scope

Equilibrium only; no SLAYER/Tearing code touched. Split out of the FKR-width work as an
independent fix, same shape as #399.

---

⚠️ NO MERGE WITHOUT THIRD-PARTY HUMAN REVIEW — NON-NEGOTIABLE ⚠️

🤖 Generated with Claude Code

https://claude.ai/code/session_01PSrf6JCViFfVzqzkQ66o6b