From ae6b37a52a453802b4215370d181b723607b25e0 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Wed, 9 Sep 2026 14:18:01 -0400
Subject: [PATCH 1/9] ForceFreeStates - BUGFIX! - Keep both crossings of a
 rational surface on reverse-shear profiles
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The singular-surface finder identified a surface by its q value, so the two
crossings of q = m/n on either side of an off-axis q minimum were merged into
one entry carrying m = [m, m]: the multi-n merge path intended for (2,1)/(4,2)
collisions. The outer crossing then had no chunk boundary, no singfac_min
standoff and no ideal jump, the sweep crawled through an unregularized
1/(m - nq) (hundreds of thousands of steps), the Δ′ BVP was skipped as a
multi-resonance surface, and the free-boundary energy was wrong by a factor
of two against STRIDE on DIII-D 153072 (q₀ = 2.44, q_min = 1.86).

Identify surfaces by ψ instead, in both sing_find! and the auto-grid knot
list, so a genuine multi-(m, n) collision still merges while distinct
crossings stay distinct surfaces with their own ψ and q′ sign.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/ForceFreeStates/Surfaces/Finding.jl | 31 +++++++++++++++----------
 1 file changed, 19 insertions(+), 12 deletions(-)

diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl
index 2c7bdac14..6f0f097dd 100644
--- a/src/ForceFreeStates/Surfaces/Finding.jl
+++ b/src/ForceFreeStates/Surfaces/Finding.jl
@@ -42,21 +42,25 @@ end
     rational_psi_nodes(equil::Equilibrium.PlasmaEquilibrium; nlow::Int, nhigh::Int=nlow)
 
 Unique ψ_N locations of all rational surfaces q = m/n for n in `nlow:nhigh`, sorted
-increasing. Used as mandatory knots for the two-pass equilibrium grid refinement (the
-same physical surface reached through several (m, n) pairs is deduplicated by q value).
+increasing. Used as mandatory knots for the two-pass equilibrium grid refinement. A surface
+is identified by its ψ, not by its q: the same physical surface reached through several
+(m, n) pairs collapses to one node, while a reverse-shear profile that crosses the same
+rational on both sides of its q minimum keeps both crossings.
 """
 function rational_psi_nodes(equil::Equilibrium.PlasmaEquilibrium; nlow::Int, nhigh::Int=nlow)
     surfaces = _find_rational_surfaces(equil, nlow, nhigh)
     nodes = Float64[]
-    qs = Float64[]
     for s in surfaces
-        any(q -> isapprox(q, s.m / s.n; atol=1e-8), qs) && continue
-        push!(qs, s.m / s.n)
+        any(psi -> _same_surface(psi, s.psifac), nodes) && continue
         push!(nodes, s.psifac)
     end
     return sort!(nodes)
 end
 
+# Two rational-surface roots are the same physical surface when they sit at the same ψ.
+const SURFACE_PSI_TOLERANCE = 1e-8
+_same_surface(psi_a::Float64, psi_b::Float64) = isapprox(psi_a, psi_b; atol=SURFACE_PSI_TOLERANCE)
+
 """
     sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
 
@@ -70,10 +74,11 @@ function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEqui
 
     for s in _find_rational_surfaces(equil, intr.nlow, intr.nhigh)
         m, n, psifac = s.m, s.n, s.psifac
-        if any(sg -> isapprox(sg.q, m / n; atol=1e-8), intr.sing)
-            # Rational surface with multiplicity > 1, add this m,n to the resonant mode numbers
-            # Technically only need m or n, but simplifies some later code and cheap to store both
-            idx = findfirst(sg -> isapprox(sg.q, m / n; atol=1e-8), intr.sing)
+        # One physical surface can be resonant with several (m, n) of equal ratio in a multi-n run;
+        # it is identified by ψ, so the two crossings of one rational on a reverse-shear profile
+        # (q1 < 0 inside the q minimum, q1 > 0 outside) stay separate surfaces.
+        idx = findfirst(sg -> _same_surface(sg.psifac, psifac), intr.sing)
+        if idx !== nothing
             push!(intr.sing[idx].m, m)
             push!(intr.sing[idx].n, n)
         else
@@ -125,9 +130,11 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     # strategy. Multi-n runs are not supported — the "outermost rational + dmlim/n" cutoff depends
     # on which n is used — and fall back to qhigh / psihigh truncation with a warning.
     if ctrl.set_psilim_via_dmlim && intr.nlow <= 0
-        error("sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " *
-              "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " *
-              "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead.")
+        error(
+            "sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " *
+            "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " *
+            "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead."
+        )
     elseif ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh
         @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation."
     elseif ctrl.set_psilim_via_dmlim

From b3018d56d705d745a0ce9b70176a9fe89037d383 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Wed, 23 Sep 2026 12:34:30 -0400
Subject: [PATCH 2/9] ForceFreeStates - BUGFIX - Match rational surfaces on (q,
 psi) and name unresolvable crossings

Deduplicate rational-surface roots on both q and psi, so two different
rationals can never merge even where steep q puts them within tolerance in
psi. Replace the chunk-bounds assert with an error naming the surface when its
singfac_min standoff overruns, which happens for crossings near a q extremum.
Add a unit test that swaps a synthetic reverse-shear q into the Solovev fixture
and checks both q=2 crossings, their opposite-sign shear, and that (2,1) and
(4,2) still merge on each.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/ForceFreeStates/EulerLagrange.jl    |  8 +++++--
 src/ForceFreeStates/Surfaces/Finding.jl | 28 +++++++++++------------
 test/runtests_sing.jl                   | 30 +++++++++++++++++++++++++
 3 files changed, 49 insertions(+), 17 deletions(-)

diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index 6245d1565..3f2abf983 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -707,8 +707,12 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro
             psi_end = intr.sing[ising_current].psifac - ctrl.singfac_min /
                                                         abs(minimum(intr.sing[ising_current].n) * intr.sing[ising_current].q1)
 
-            # Validate chunk bounds
-            @assert psi_current < psi_end "Invalid chunk bounds: psi_start=$psi_current >= psi_end=$psi_end"
+            # A standoff reaching past the previous boundary means |q'| is too small here (e.g. near a q extremum).
+            if psi_end <= psi_current
+                sg = intr.sing[ising_current]
+                error("Singular surface q=$(sg.q) at ψ=$(sg.psifac) (q'=$(sg.q1)): its singfac_min standoff reaches past ψ=$psi_current. " *
+                      "Crossings this close to a q extremum cannot be resolved; reduce singfac_min or move q_min away from this rational.")
+            end
             @assert isempty(chunks) || psi_current >= chunks[end].psi_end "Overlapping chunks detected"
 
             push!(chunks, IntegrationChunk(;
diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl
index 6f0f097dd..1c943adb9 100644
--- a/src/ForceFreeStates/Surfaces/Finding.jl
+++ b/src/ForceFreeStates/Surfaces/Finding.jl
@@ -42,24 +42,24 @@ end
     rational_psi_nodes(equil::Equilibrium.PlasmaEquilibrium; nlow::Int, nhigh::Int=nlow)
 
 Unique ψ_N locations of all rational surfaces q = m/n for n in `nlow:nhigh`, sorted
-increasing. Used as mandatory knots for the two-pass equilibrium grid refinement. A surface
-is identified by its ψ, not by its q: the same physical surface reached through several
-(m, n) pairs collapses to one node, while a reverse-shear profile that crosses the same
-rational on both sides of its q minimum keeps both crossings.
+increasing. Used as mandatory knots for the two-pass equilibrium grid refinement. Roots
+are deduplicated by (q, ψ), so several (m, n) of equal ratio collapse to one node while the
+two crossings of one rational on a reverse-shear profile stay distinct.
 """
 function rational_psi_nodes(equil::Equilibrium.PlasmaEquilibrium; nlow::Int, nhigh::Int=nlow)
     surfaces = _find_rational_surfaces(equil, nlow, nhigh)
-    nodes = Float64[]
+    kept = @NamedTuple{q::Float64, psifac::Float64}[]
     for s in surfaces
-        any(psi -> _same_surface(psi, s.psifac), nodes) && continue
-        push!(nodes, s.psifac)
+        q = s.m / s.n
+        any(k -> _same_surface(k.q, k.psifac, q, s.psifac), kept) && continue
+        push!(kept, (q=q, psifac=s.psifac))
     end
-    return sort!(nodes)
+    return sort!([k.psifac for k in kept])
 end
 
-# Two rational-surface roots are the same physical surface when they sit at the same ψ.
-const SURFACE_PSI_TOLERANCE = 1e-8
-_same_surface(psi_a::Float64, psi_b::Float64) = isapprox(psi_a, psi_b; atol=SURFACE_PSI_TOLERANCE)
+# Two rational-surface roots are one physical surface when both q and ψ agree.
+const SURFACE_TOLERANCE = 1e-8
+_same_surface(q_a, psi_a, q_b, psi_b) = isapprox(q_a, q_b; atol=SURFACE_TOLERANCE) && isapprox(psi_a, psi_b; atol=SURFACE_TOLERANCE)
 
 """
     sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
@@ -74,10 +74,8 @@ function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEqui
 
     for s in _find_rational_surfaces(equil, intr.nlow, intr.nhigh)
         m, n, psifac = s.m, s.n, s.psifac
-        # One physical surface can be resonant with several (m, n) of equal ratio in a multi-n run;
-        # it is identified by ψ, so the two crossings of one rational on a reverse-shear profile
-        # (q1 < 0 inside the q minimum, q1 > 0 outside) stay separate surfaces.
-        idx = findfirst(sg -> _same_surface(sg.psifac, psifac), intr.sing)
+        # Merge equal-ratio (m, n) at one surface; reverse-shear crossings of one rational stay separate.
+        idx = findfirst(sg -> _same_surface(sg.q, sg.psifac, m / n, psifac), intr.sing)
         if idx !== nothing
             push!(intr.sing[idx].m, m)
             push!(intr.sing[idx].n, n)
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index d58d11c2b..b6b79fae6 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -172,6 +172,36 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         end
     end
 
+    @testset "sing_find reverse shear" begin
+        FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))
+        # Swap in a synthetic reverse-shear q (q_min = 1.8 at ψ = 0.4) so q = 2 is crossed twice.
+        pr = equil.profiles
+        q_rs = @. 1.8 + 3 * (pr.xs - 0.4)^2
+        equil.profiles = GeneralizedPerturbedEquilibrium.Equilibrium.ProfileSplines(pr.xs, pr.F_spline.y, pr.P_spline.y, pr.dVdpsi_spline.y, q_rs)
+        GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_qfind!(equil)
+
+        intr = FFS.ForceFreeStatesInternal()
+        intr.nlow = 1
+        intr.nhigh = 2
+        intr.msing = 0
+        intr.sing = FFS.SingType[]
+        FFS.sing_find!(intr, equil)
+
+        q2 = filter(s -> isapprox(s.q, 2.0; atol=1e-8), intr.sing)
+        psi_q2 = 0.4 .+ [-1, 1] .* sqrt(0.2 / 3)
+        @test length(q2) == 2
+        @test [s.psifac for s in q2] ≈ psi_q2 atol = 1e-6
+        @test q2[1].q1 < 0 < q2[2].q1
+        # (2,1) and (4,2) still merge onto each crossing.
+        @test all(s -> sort(s.m) == [2, 4] && sort(s.n) == [1, 2], q2)
+        @test length(intr.sing) == intr.msing
+
+        nodes = FFS.rational_psi_nodes(equil; nlow=1, nhigh=2)
+        @test count(p -> any(isapprox.(p, psi_q2; atol=1e-6)), nodes) == 2
+        @test length(nodes) == intr.msing
+    end
+
     # ---------------------------------
     # sing_lim
     # ---------------------------------

From 851ff66c215cfeeae2fbbf11c2056cc023e24f62 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 11:04:55 -0400
Subject: [PATCH 3/9] ForceFreeStates - MINOR - State the fact in the
 crossing-standoff error and document the surface merge rule

The chunk-bounds error no longer asserts a q extremum or recommends reducing singfac_min; it reports the surface, q', the standoff and the chunk start. sing_find! documents the equal-ratio merge and reverse-shear split in its docstring and drops the Fortran citation. Refresh a stale dedupe comment in bracket_mandatory_nodes and tidy the reverse-shear test's isapprox call.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/GridRefinement.jl       | 2 +-
 src/ForceFreeStates/EulerLagrange.jl    | 6 +++---
 src/ForceFreeStates/Surfaces/Finding.jl | 8 ++++----
 test/runtests_sing.jl                   | 3 ++-
 4 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl
index e58eb81dc..651b897e3 100644
--- a/src/Equilibrium/GridRefinement.jl
+++ b/src/Equilibrium/GridRefinement.jl
@@ -387,7 +387,7 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}
     for idx in order
         c, hw = centers[idx], min_half_widths[idx]
         (lo < c < hi) || continue
-        c - last_c < collapse_atol && continue  # duplicate rational (same q via several (m,n))
+        c - last_c < collapse_atol && continue  # centers closer than collapse_atol share one bracket
         k = clamp(searchsortedlast(grid, c), 1, length(grid) - 1)
         w = max(hw, 0.5 * (grid[k+1] - grid[k]))  # blend to half the local spacing
         left, right = c - w, c + w
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index 3f2abf983..b43200939 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -707,11 +707,11 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro
             psi_end = intr.sing[ising_current].psifac - ctrl.singfac_min /
                                                         abs(minimum(intr.sing[ising_current].n) * intr.sing[ising_current].q1)
 
-            # A standoff reaching past the previous boundary means |q'| is too small here (e.g. near a q extremum).
+            # The crossing standoff singfac_min/|n q'| must not reach back past the current chunk start.
             if psi_end <= psi_current
                 sg = intr.sing[ising_current]
-                error("Singular surface q=$(sg.q) at ψ=$(sg.psifac) (q'=$(sg.q1)): its singfac_min standoff reaches past ψ=$psi_current. " *
-                      "Crossings this close to a q extremum cannot be resolved; reduce singfac_min or move q_min away from this rational.")
+                error("Singular surface q=$(sg.q) at ψ=$(sg.psifac) (q'=$(sg.q1)): its crossing standoff " *
+                      "singfac_min/|n q'| = $(sg.psifac - psi_end) reaches back past the current chunk start ψ=$psi_current.")
             end
             @assert isempty(chunks) || psi_current >= chunks[end].psi_end "Overlapping chunks detected"
 
diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl
index 1c943adb9..02a5b252f 100644
--- a/src/ForceFreeStates/Surfaces/Finding.jl
+++ b/src/ForceFreeStates/Surfaces/Finding.jl
@@ -64,9 +64,10 @@ _same_surface(q_a, psi_a, q_b, psi_b) = isapprox(q_a, q_b; atol=SURFACE_TOLERANC
 """
     sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
 
-Locate singular rational q-surfaces (q = m/nn) using a bisection method
-between extrema of the q-profile, and store their properties in `intr.sing`.
-Performs the same function as `sing_find` in the Fortran code.
+Locate singular rational q-surfaces (q = m/n) by bisection between extrema of the q-profile,
+and store their properties in `intr.sing`, sorted by increasing ψ. Roots of equal ratio m/n at
+the same ψ (e.g. (2,1) and (4,2)) merge into one surface carrying all its (m, n); the two
+crossings of one rational on a reverse-shear profile stay separate surfaces with opposite-sign q′.
 """
 function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
     profiles = equil.profiles
@@ -74,7 +75,6 @@ function sing_find!(intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEqui
 
     for s in _find_rational_surfaces(equil, intr.nlow, intr.nhigh)
         m, n, psifac = s.m, s.n, s.psifac
-        # Merge equal-ratio (m, n) at one surface; reverse-shear crossings of one rational stay separate.
         idx = findfirst(sg -> _same_surface(sg.q, sg.psifac, m / n, psifac), intr.sing)
         if idx !== nothing
             push!(intr.sing[idx].m, m)
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index b6b79fae6..e269b3bbb 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -176,6 +176,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
         equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))
         # Swap in a synthetic reverse-shear q (q_min = 1.8 at ψ = 0.4) so q = 2 is crossed twice.
+        # sing_find! reads only q, so F, P and geometry are left inconsistent with it.
         pr = equil.profiles
         q_rs = @. 1.8 + 3 * (pr.xs - 0.4)^2
         equil.profiles = GeneralizedPerturbedEquilibrium.Equilibrium.ProfileSplines(pr.xs, pr.F_spline.y, pr.P_spline.y, pr.dVdpsi_spline.y, q_rs)
@@ -191,7 +192,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         q2 = filter(s -> isapprox(s.q, 2.0; atol=1e-8), intr.sing)
         psi_q2 = 0.4 .+ [-1, 1] .* sqrt(0.2 / 3)
         @test length(q2) == 2
-        @test [s.psifac for s in q2] ≈ psi_q2 atol = 1e-6
+        @test isapprox([s.psifac for s in q2], psi_q2; atol=1e-6)
         @test q2[1].q1 < 0 < q2[2].q1
         # (2,1) and (4,2) still merge onto each crossing.
         @test all(s -> sort(s.m) == [2, 4] && sort(s.n) == [1, 2], q2)

From 7b1a1de188cc7911d60ebac8577ee34f8dcf8880 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 11:04:55 -0400
Subject: [PATCH 4/9] ForceFreeStates - TEST - Cover overlapping and tangent
 reverse-shear crossing standoffs
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Two crossings of q = 2 on an analytic parabola resolve when 2ε > singfac_min and raise the named error when their standoffs overlap; a tangent surface (q' = 0) also raises it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 test/runtests_eulerlagrange.jl | 34 ++++++++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl
index 497e1cff0..4558fa4d1 100644
--- a/test/runtests_eulerlagrange.jl
+++ b/test/runtests_eulerlagrange.jl
@@ -440,6 +440,40 @@ end
         chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr)
         @test length(chunks) == 1
         @test chunks[1].needs_crossing == false
+
+        # Case 5: both crossings of q = 2 on q = (2 - ε) + a(ψ - ψ0)², at ψ0 ∓ √(ε/a) with q' = ∓2√(aε).
+        # With n = 1 their standoffs singfac_min/|q'| overlap exactly when 2ε ≤ singfac_min.
+        FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        ctrl = FFS.ForceFreeStatesControl(; numsteps_init=10, numunorms_init=5, singfac_min=1e-4)
+        function reverse_shear_pair(ε; a=3.0, psi0=0.4)
+            intr = FFS.ForceFreeStatesInternal(; mpert=1, numpert_total=1)
+            dpsi, q1 = sqrt(ε / a), 2 * sqrt(a * ε)
+            intr.sing = [FFS.SingType(; psifac=psi0 - dpsi, n=[1], m=[2], q=2.0, q1=-q1), FFS.SingType(; psifac=psi0 + dpsi, n=[1], m=[2], q=2.0, q1=q1)]
+            intr.msing = 2
+            intr.psilim = 1.0
+            intr.mlow = 1
+            intr.mhigh = 3
+            odet = FFS.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
+            odet.psifac = 0.0
+            return odet, intr
+        end
+        odet, intr = reverse_shear_pair(1e-3)
+        chunks = FFS.chunk_el_integration_bounds(odet, ctrl, intr)
+        @test length(chunks) == 3
+        @test all(c.needs_crossing for c in chunks[1:2])
+        odet, intr = reverse_shear_pair(1e-5)
+        @test_throws ErrorException FFS.chunk_el_integration_bounds(odet, ctrl, intr)
+
+        # Case 6: a surface tangent to q_min (q' = 0) has an unbounded standoff.
+        intr = FFS.ForceFreeStatesInternal(; mpert=1, numpert_total=1)
+        intr.sing = [FFS.SingType(; psifac=0.4, n=[1], m=[2], q=2.0, q1=0.0)]
+        intr.msing = 1
+        intr.psilim = 1.0
+        intr.mlow = 1
+        intr.mhigh = 3
+        odet = FFS.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
+        odet.psifac = 0.0
+        @test_throws ErrorException FFS.chunk_el_integration_bounds(odet, ctrl, intr)
     end
 
     @testset "EdgeScanState" begin

From c920ee7dc9d4b051a04e52b0375893e41125e8cb Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 11:09:42 -0400
Subject: [PATCH 5/9] Repo - MINOR - Apply the formatter to the audit fixes

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/ForceFreeStates/EulerLagrange.jl |  6 ++++--
 test/runtests_sing.jl                | 20 ++++++++++----------
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index b43200939..c3b2b15c6 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -710,8 +710,10 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro
             # The crossing standoff singfac_min/|n q'| must not reach back past the current chunk start.
             if psi_end <= psi_current
                 sg = intr.sing[ising_current]
-                error("Singular surface q=$(sg.q) at ψ=$(sg.psifac) (q'=$(sg.q1)): its crossing standoff " *
-                      "singfac_min/|n q'| = $(sg.psifac - psi_end) reaches back past the current chunk start ψ=$psi_current.")
+                error(
+                    "Singular surface q=$(sg.q) at ψ=$(sg.psifac) (q'=$(sg.q1)): its crossing standoff " *
+                    "singfac_min/|n q'| = $(sg.psifac - psi_end) reaches back past the current chunk start ψ=$psi_current."
+                )
             end
             @assert isempty(chunks) || psi_current >= chunks[end].psi_end "Overlapping chunks detected"
 
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index e269b3bbb..999de8072 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -41,14 +41,14 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
 
     function read_solutions_3d(fname::String)
         lines = readlines(fname)
-        blocks = Vector{Vector{Vector{Float64}}}();
+        blocks = Vector{Vector{Vector{Float64}}}()
         current = Vector{Vector{Float64}}()
         for s in lines
-            t = strip(s);
+            t = strip(s)
             isempty(t) && continue
             if occursin("Solution index", t)
                 if !isempty(current)
-                    push!(blocks, current);
+                    push!(blocks, current)
                     current = Vector{Vector{Float64}}()
                 end
                 continue
@@ -56,24 +56,24 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
             if occursin("=", t)
                 continue
             end
-            vals = [parse(Float64, replace(x, 'D'=>'E')) for x in split(t)]
+            vals = [parse(Float64, replace(x, 'D' => 'E')) for x in split(t)]
             push!(current, vals)
         end
         if !isempty(current)
             push!(blocks, current)
         end
-        mpert = length(blocks[1]);
+        mpert = length(blocks[1])
         nsol = length(blocks)
         result = Array{ComplexF64}(undef, mpert, nsol, 2)
         for (j, block) in enumerate(blocks), (i, row) in enumerate(block)
-            result[i, j, 1] = complex(row[2], row[3]);
+            result[i, j, 1] = complex(row[2], row[3])
             result[i, j, 2] = complex(row[4], row[5])
         end
         return result
     end
 
     function copyForSplines(prev_mat, num_range)
-        mmat = zeros(ComplexF64, length(num_range), size(prev_mat, 1)*size(prev_mat, 2))
+        mmat = zeros(ComplexF64, length(num_range), size(prev_mat, 1) * size(prev_mat, 2))
         for i in 1:length(num_range)
             mmat[i, :] .= vec(prev_mat)
         end
@@ -105,7 +105,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(; numpert_total=intr.numpert_total,
             numsteps_init=ctrl.numsteps_init, numunorms_init=ctrl.numunorms_init, msing=intr.msing)
 
-        psifac_dummy = collect(range(0, 1, 10));
+        psifac_dummy = collect(range(0, 1, 10))
         points = length(psifac_dummy)
         amat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/amat.dat"))
         amats = copyForSplines(amat, psifac_dummy)
@@ -115,7 +115,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         cmats = copyForSplines(cmat, psifac_dummy)
 
         fmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/fmat.dat"))
-        fmat .= cholesky(Hermitian(fmat)).L;
+        fmat .= cholesky(Hermitian(fmat)).L
         fmats = copyForSplines(fmat, psifac_dummy)
         kmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/kmat.dat"))
         kmats = copyForSplines(kmat, psifac_dummy)
@@ -125,7 +125,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         umat_p1 = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/umat_p1.dat"))
         umat_p2 = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/umat_p2.dat"))
         odet.psifac = extract_value(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/sing_der_output_normal.dat"), "psifac")
-        odet.u[:, :, 1] .= umat_p1;
+        odet.u[:, :, 1] .= umat_p1
         odet.u[:, :, 2] .= umat_p2
 
         itp_opts = (; extrap=ExtendExtrap())

From 10d4cc7d469a8584445cb485e704f45ce2273673 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 15:38:52 -0400
Subject: [PATCH 6/9] ForceFreeStates - BUGFIX - Reject a qlow between qmin and
 q0 on reversed-shear profiles

When qmin < qlow < q0 the q < qlow region is a ring around q_min rather than a core, and sing_min! (outer leg) and initialize_el_at_axis! (axis) disagree about the lower integration bound. check_qlow_domain raises a named error at domain setup saying what was detected and to set qlow <= qmin or qlow >= q0. Covered on the analytic reverse-shear fixture in runtests_sing.jl.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/ForceFreeStates/Surfaces/Finding.jl | 21 +++++++++++++++++++++
 src/GeneralizedPerturbedEquilibrium.jl  |  3 ++-
 test/runtests_sing.jl                   |  7 +++++++
 3 files changed, 30 insertions(+), 1 deletion(-)

diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl
index 02a5b252f..5c1e0a078 100644
--- a/src/ForceFreeStates/Surfaces/Finding.jl
+++ b/src/ForceFreeStates/Surfaces/Finding.jl
@@ -163,6 +163,27 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     end
 end
 
+"""
+    check_qlow_domain(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
+
+Error when `qmin < qlow < q0`. On such a reversed-shear profile the q < qlow region is a ring
+around q_min rather than a core, so there is no single lower bound that excludes it: `sing_min!`
+places `psilow` on the outer leg, while `initialize_el_at_axis!` starts at the axis because
+`qlow < q0`, and the two disagree about the integration domain.
+"""
+function check_qlow_domain(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
+    q0 = equil.profiles.q_spline.y[1]
+    qmin = equil.params.qmin
+    if ctrl.qlow > qmin && q0 > ctrl.qlow
+        error(
+            "qlow=$(ctrl.qlow) lies between qmin=$qmin and the axis value q0=$q0 of a reversed-shear q profile: " *
+            "the q < qlow region is a ring around q_min, not a core, so no single lower integration bound excludes it. " *
+            "Set qlow ≤ qmin to integrate from the axis, or qlow ≥ q0 to start outside q_min."
+        )
+    end
+    return nothing
+end
+
 """
     sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
 
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index f2ed4ef65..e2e22244b 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -72,7 +72,7 @@ include("Rerun.jl")
 # Import ForceFreeStates types and functions needed for main
 using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings
 using .ForceFreeStates: ForceFreeStatesResult, build_result
-using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry
+using .ForceFreeStates: sing_lim!, check_qlow_domain, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry
 using .ForceFreeStates: make_metric, build_matrix_splines, build_kinetic_matrix_splines
 using .ForceFreeStates: find_kinetic_singular_surfaces!
 using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions!
@@ -496,6 +496,7 @@ function prepare_force_free_states!(
 )
     # Determine psilim and qlim (where we will integrate to)
     sing_lim!(intr, ctrl, equil)
+    check_qlow_domain(ctrl, equil)
 
     # Find all singular surfaces in the equilibrium
     sing_find!(intr, equil)
diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl
index 999de8072..cf2ad7fe6 100644
--- a/test/runtests_sing.jl
+++ b/test/runtests_sing.jl
@@ -201,6 +201,13 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         nodes = FFS.rational_psi_nodes(equil; nlow=1, nhigh=2)
         @test count(p -> any(isapprox.(p, psi_q2; atol=1e-6)), nodes) == 2
         @test length(nodes) == intr.msing
+
+        # qmin < qlow < q0 makes the q < qlow region a ring around q_min, which no single lower bound excludes.
+        q0 = equil.profiles.q_spline.y[1]
+        @test equil.params.qmin < 2.0 < q0
+        @test_throws ErrorException FFS.check_qlow_domain(FFS.ForceFreeStatesControl(; qlow=2.0), equil)
+        @test FFS.check_qlow_domain(FFS.ForceFreeStatesControl(; qlow=1.5), equil) === nothing
+        @test FFS.check_qlow_domain(FFS.ForceFreeStatesControl(; qlow=q0 + 0.1), equil) === nothing
     end
 
     # ---------------------------------

From 0df203cf335aca05b01084e0ea1c2897565b927b Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 15:38:52 -0400
Subject: [PATCH 7/9] Equilibrium - BUGFIX - Error when adjacent
 rational-surface grid brackets overlap

bracket_mandatory_nodes now errors when c2 - w2 <= c1 + w1 for adjacent brackets, or when a half-width is non-finite (tangent surface), instead of interleaving bracket knots and leaving a surface off-center or sharing a spline interval. The collapse_atol dedupe is kept as a safety net. Geometric condition only; no new tolerance.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/Equilibrium/GridRefinement.jl | 15 ++++++++++++---
 test/runtests_grid_refinement.jl  |  7 +++++++
 2 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl
index 651b897e3..ee771afd8 100644
--- a/src/Equilibrium/GridRefinement.jl
+++ b/src/Equilibrium/GridRefinement.jl
@@ -375,7 +375,8 @@ surface (no knot-on-surface jump) and *stable across* grid refinement — a narr
 wide neighbors would be consistent but noisy. Non-bracket knots within `w` of the center, or
 within `min_spacing` of either new bracket knot, are dropped; endpoints always win; bracket knots
 outside the open domain are dropped, and knots closer than `collapse_atol` are collapsed to keep
-the grid strictly increasing.
+the grid strictly increasing. Errors if adjacent brackets overlap (`c₂ − w₂ ≤ c₁ + w₁`) or a
+half-width is non-finite (a tangent surface, q′ = 0).
 """
 function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}, min_half_widths::Vector{Float64}, min_spacing::Float64; collapse_atol::Float64=1e-7)
     isempty(centers) && return copy(grid)
@@ -383,7 +384,7 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}
     lo, hi = grid[1], grid[end]
     order = sortperm(centers)
     kept = Tuple{Float64,Bool}[(g, false) for g in grid]  # (value, is_bracket_knot)
-    last_c = -Inf
+    last_c, last_w = -Inf, 0.0
     for idx in order
         c, hw = centers[idx], min_half_widths[idx]
         (lo < c < hi) || continue
@@ -391,6 +392,14 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}
         k = clamp(searchsortedlast(grid, c), 1, length(grid) - 1)
         w = max(hw, 0.5 * (grid[k+1] - grid[k]))  # blend to half the local spacing
         left, right = c - w, c + w
+        # Overlapping brackets would interleave their knots and leave a surface off-center or sharing an interval.
+        isfinite(w) || error("Rational surface at ψ=$c: its grid bracket half-width $w is not finite (tangent surface, q′ = 0).")
+        if left <= last_c + last_w
+            error(
+                "Rational surface at ψ=$c: its grid bracket half-width $w reaches back past the bracket " *
+                "of the previous rational surface at ψ=$last_c (half-width $last_w)."
+            )
+        end
         # Drop non-bracket, non-endpoint knots inside the zone or crowding a new bracket knot.
         filter!(kept) do t
             t[2] || t[1] == lo || t[1] == hi ||
@@ -398,7 +407,7 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}
         end
         left > lo && push!(kept, (left, true))
         right < hi && push!(kept, (right, true))
-        last_c = c
+        last_c, last_w = c, w
     end
     sort!(kept; by=first)
     merged = Float64[]
diff --git a/test/runtests_grid_refinement.jl b/test/runtests_grid_refinement.jl
index 35c350c90..55bb2e9e2 100644
--- a/test/runtests_grid_refinement.jl
+++ b/test/runtests_grid_refinement.jl
@@ -108,6 +108,13 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium
         for c in (0.3, 0.7)
             @test !any(g -> abs(g - c) < 1e-12, b)
         end
+
+        # Adjacent brackets that overlap (c₂ − w₂ ≤ c₁ + w₁), e.g. two reverse-shear crossings near q_min, error.
+        @test_throws ErrorException GridRef.bracket_mandatory_nodes(grid, [0.5, 0.51], [0.01, 0.01], 1e-4)
+        # A tangent surface (q′ = 0) has an infinite bracket half-width.
+        @test_throws ErrorException GridRef.bracket_mandatory_nodes(grid, [0.5], [Inf], 1e-4)
+        # Centers within collapse_atol still share one bracket instead of erroring.
+        @test all(diff(GridRef.bracket_mandatory_nodes(grid, [0.5, 0.5 + 1e-9], [0.01, 0.01], 1e-4)) .> 0)
     end
 
     @testset "_validate_psi_nodes" begin

From f44496288e4af177caaa2f5bb8ecbfdd04a6fcd2 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 15:38:52 -0400
Subject: [PATCH 8/9] ForceFreeStates - MINOR - Warn when the kinetic path
 skips a crossing with an overlapping standoff

The kinetic chunking branch still skips a kinetic singular surface whose crossing standoff reaches back past the chunk start, but now warns with its m/n, psi, q' and standoff. No result changes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/ForceFreeStates/EulerLagrange.jl | 5 ++++-
 test/runtests_eulerlagrange.jl       | 8 ++++++++
 2 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index c3b2b15c6..7a58e640b 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -664,7 +664,10 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro
                       ctrl.singfac_min / abs(minimum(intr.kinsing[ising_current].n) * intr.kinsing[ising_current].q1)
 
             if psi_current >= psi_end
-                # Surface too close to current position — skip it
+                # The crossing standoff reaches back past the current chunk start: skip this surface, loudly.
+                ks = intr.kinsing[ising_current]
+                @warn "Skipping kinetic singular surface m/n=$(join(string.(ks.m, "/", ks.n), ", ")) at ψ=$(ks.psifac) (q'=$(ks.q1)): " *
+                      "its crossing standoff singfac_min/|n q'| = $(ks.psifac - psi_end) reaches back past the current chunk start ψ=$psi_current."
                 ising_current = find_next_kinsing!(ising_current, intr)
                 continue
             end
diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl
index 4558fa4d1..f63249ca0 100644
--- a/test/runtests_eulerlagrange.jl
+++ b/test/runtests_eulerlagrange.jl
@@ -474,6 +474,14 @@ end
         odet = FFS.OdeState(1, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
         odet.psifac = 0.0
         @test_throws ErrorException FFS.chunk_el_integration_bounds(odet, ctrl, intr)
+
+        # Case 7: the kinetic path skips the second of two overlapping crossings, with a warning naming it.
+        ctrl_kin = FFS.ForceFreeStatesControl(; numsteps_init=10, numunorms_init=5, singfac_min=1e-4, kinetic_factor=1.0)
+        odet, intr = reverse_shear_pair(1e-5)
+        intr.kinsing, intr.kmsing = intr.sing, intr.msing
+        chunks = @test_logs (:warn, r"Skipping kinetic singular surface m/n=2/1") match_mode = :any FFS.chunk_el_integration_bounds(odet, ctrl_kin, intr)
+        @test length(chunks) == 2
+        @test chunks[1].needs_crossing && !chunks[2].needs_crossing
     end
 
     @testset "EdgeScanState" begin

From aac74ca2a2c1eb298a9bb36339ca4ce11595e8f9 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Fri, 25 Sep 2026 15:43:21 -0400
Subject: [PATCH 9/9] Repo - MINOR - Apply the formatter to the audit fixes

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/GeneralizedPerturbedEquilibrium.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index e2e22244b..6af7821d9 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -85,8 +85,8 @@ using .Equilibrium: PlasmaEquilibrium
 using .ForcingTerms: RMPField
 
 const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim",
-                              "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi",
-                              "gal_flag")
+    "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi",
+    "gal_flag")
 const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc")
 
 # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing
@@ -741,10 +741,10 @@ runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces]
 and errors here.
 
 ```julia
-eq   = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
+eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
 prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true)
-ffs  = solve(prob, Riccati())
-ffs  = solve(eq, Riccati(); nn=1, vac_flag=true)   # equivalent one-line form
+ffs = solve(prob, Riccati())
+ffs = solve(eq, Riccati(); nn=1, vac_flag=true)   # equivalent one-line form
 ```
 """
 function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrator)