ForceFreeStates - BUGFIX! - ⚠️ Keep both crossings of a rational surface on reverse-shear profiles

src/Equilibrium/GridRefinement.jl

@@ -375,30 +375,39 @@ 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)
     length(centers) == length(min_half_widths) || error("bracket_mandatory_nodes: centers and min_half_widths length mismatch")
     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
-        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
+        # 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 ||
                 !(left < t[1] < right || abs(t[1] - left) < min_spacing || abs(t[1] - right) < min_spacing)
         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[]

src/ForceFreeStates/EulerLagrange.jl

@@ -634,7 +634,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
@@ -677,8 +680,14 @@ 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"
+            # 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."
+                )
+            end
             @assert isempty(chunks) || psi_current >= chunks[end].psi_end "Overlapping chunks detected"
 
             push!(chunks, IntegrationChunk(;

src/ForceFreeStates/Surfaces/Finding.jl

@@ -42,38 +42,41 @@ 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. 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[]
-    qs = Float64[]
+    kept = @NamedTuple{q::Float64, psifac::Float64}[]
     for s in surfaces
-        any(q -> isapprox(q, s.m / s.n; atol=1e-8), qs) && continue
-        push!(qs, s.m / s.n)
-        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 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)
 
-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
     hint = Ref(1)
 
     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)
+        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)
         else
@@ -125,9 +128,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
@@ -158,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
+
 """
     remove_singular_surfs!(intr; qmin=-Inf, psimax=intr.psilim) -> keep
 

src/GeneralizedPerturbedEquilibrium.jl

@@ -77,7 +77,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!, remove_singular_surfs!, resist_eval_all!, resist_geometry, ResistGeometry
+using .ForceFreeStates: sing_lim!, check_qlow_domain, sing_min!, sing_find!, remove_singular_surfs!, 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!
@@ -558,6 +558,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)

test/runtests_eulerlagrange.jl

@@ -494,6 +494,48 @@ 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)
+
+        # 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

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