From ae6b37a52a453802b4215370d181b723607b25e0 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 9 Sep 2026 14:18:01 -0400 Subject: [PATCH 01/16] 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 --- 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 Date: Wed, 23 Sep 2026 12:34:30 -0400 Subject: [PATCH 02/16] 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 --- 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 Date: Fri, 25 Sep 2026 11:04:55 -0400 Subject: [PATCH 03/16] 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 --- 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 Date: Fri, 25 Sep 2026 11:04:55 -0400 Subject: [PATCH 04/16] 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 --- 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 Date: Fri, 25 Sep 2026 11:09:42 -0400 Subject: [PATCH 05/16] Repo - MINOR - Apply the formatter to the audit fixes Co-Authored-By: Claude Opus 5.5 --- 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 Date: Fri, 25 Sep 2026 15:38:52 -0400 Subject: [PATCH 06/16] 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 --- 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 Date: Fri, 25 Sep 2026 15:38:52 -0400 Subject: [PATCH 07/16] 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 --- 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 Date: Fri, 25 Sep 2026 15:38:52 -0400 Subject: [PATCH 08/16] 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 --- 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 Date: Fri, 25 Sep 2026 15:43:21 -0400 Subject: [PATCH 09/16] Repo - MINOR - Apply the formatter to the audit fixes Co-Authored-By: Claude Opus 5.5 --- 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) From 538dcf56957c0f6f308440fe9aae2f6f8e933ced Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 17:33:53 -0400 Subject: [PATCH 10/16] ForceFreeStates - BUGFIX - Seed psilim from the edge inward and check the q < qlow region on the grid itself Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/Surfaces/Finding.jl | 35 ++++++++++++++----------- test/runtests_sing.jl | 7 +++++ 2 files changed, 26 insertions(+), 16 deletions(-) diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 05f467fc9..f05423ed1 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -166,8 +166,10 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, # If truncate_at_rational_offset decreased qlim or qhigh < qmax, we need to find the precise psilim via newton iteration if intr.qlim < equil.params.qmax - # Find nearest ψ index where q ≈ qlim - _, jpsi = findmin(abs.(profiles.q_spline.y .- intr.qlim)) + # Seed from the edge inward at the outermost node with q ≤ qlim, so a reversed-shear profile, which crosses + # qlim on both legs of its minimum, gets the outer crossing. + qy = profiles.q_spline.y + jpsi = something(findlast(<=(intr.qlim), qy), argmin(abs.(qy .- intr.qlim))) jpsi = min(jpsi, length(profiles.xs) - 1) hint = Ref(jpsi) @@ -183,22 +185,23 @@ 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. +Error unless the region `q < qlow` on the ψ grid is empty or one interval starting at the axis. +On a reversed-shear profile with `q_min < qlow ≤ q0` it is instead a ring around q_min, and with +several minima it can be several intervals; no single lower bound excludes either: `sing_min!` +places `psilow` on the outermost leg, while `initialize_el_at_axis!` starts at the axis, 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 + qy = equil.profiles.q_spline.y + below = qy .< ctrl.qlow + # A core: below on the first node and on every node up to the first one at or above qlow, then never again. + (!any(below) || (below[1] && issorted(below; rev=true))) && return nothing + error( + "qlow=$(ctrl.qlow): the region q < qlow is not a single core around the axis (q on the axis is $(qy[1]), " * + "the grid minimum is $(minimum(qy))); on a reversed-shear profile it is a ring around q_min, so no single " * + "lower integration bound excludes it. Set qlow below the q minimum to integrate from the axis, or above " * + "the axis value to start outside it." + ) end """ diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index 056a69924..5b972d9b1 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -208,6 +208,13 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex @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 + # qlow equal to the axis value still leaves a ring below it. + @test_throws ErrorException FFS.check_qlow_domain(FFS.ForceFreeStatesControl(; qlow=q0), equil) + + # A qhigh below the axis value is crossed on both legs; the domain ends on the outer one. + intr_lim = FFS.ForceFreeStatesInternal() + FFS.sing_lim!(intr_lim, FFS.ForceFreeStatesControl(; qhigh=2.0, truncate_at_rational_offset=false), equil) + @test isapprox(intr_lim.psilim, psi_q2[2]; atol=1e-6) end # --------------------------------- From 36cbc5e31b383fa2b4048d75e4c8548dc2d2b1a7 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 17:41:19 -0400 Subject: [PATCH 11/16] Equilibrium - BUGFIX - Cap rational-surface brackets at half the gap to their neighbors and the domain edges instead of erroring Co-Authored-By: Claude Opus 5.5 --- src/Equilibrium/Grid/GridRefinement.jl | 39 ++++++++++++++------------ test/runtests_grid_refinement.jl | 18 +++++++++--- 2 files changed, 35 insertions(+), 22 deletions(-) diff --git a/src/Equilibrium/Grid/GridRefinement.jl b/src/Equilibrium/Grid/GridRefinement.jl index 0b0a78024..ac9ee1778 100644 --- a/src/Equilibrium/Grid/GridRefinement.jl +++ b/src/Equilibrium/Grid/GridRefinement.jl @@ -469,31 +469,35 @@ 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. Errors if adjacent brackets overlap (`c₂ − w₂ ≤ c₁ + w₁`) or a -half-width is non-finite (a tangent surface, q′ = 0). +the grid strictly increasing. Each half-width is capped at half the distance to the neighboring +centers, so adjacent brackets never interleave and every center keeps its own interval, and the +stencil part of it also at half the distance to the nearer domain edge, so a near-tangent surface +(q′ → 0, a huge stencil) cannot clear the grid. A cap that cuts a bracket below its stencil floor +warns, since Δ′ there is sampled on a narrower stencil than intended. """ 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, last_w = -Inf, 0.0 - for idx in order - c, hw = centers[idx], min_half_widths[idx] + # In-domain centers in increasing order; centers closer than collapse_atol share one bracket. + cs, hws = Float64[], Float64[] + for idx in sortperm(centers) + c = centers[idx] (lo < c < hi) || continue - c - last_c < collapse_atol && continue # centers closer than collapse_atol share one bracket + (isempty(cs) || c - cs[end] >= collapse_atol) || continue + push!(cs, c) + push!(hws, min_half_widths[idx]) + end + kept = Tuple{Float64,Bool}[(g, false) for g in grid] # (value, is_bracket_knot) + for (i, c) in enumerate(cs) 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 + cap_nb = 0.5 * min(i > 1 ? c - cs[i-1] : Inf, i < length(cs) ? cs[i+1] - c : Inf) + cap_edge = 0.5 * min(c - lo, hi - c) + # Stencil floor capped by neighbors and edges; blend to half the local spacing capped by neighbors only. + w = max(min(hws[i], cap_nb, cap_edge), min(0.5 * (grid[k+1] - grid[k]), cap_nb)) + w < hws[i] && @warn "Rational surface at ψ=$c: grid bracket half-width capped at $w, below its Δ′ stencil $(hws[i]) " * + "(neighboring surface or domain edge too close, or a near-tangent surface with q′ ≈ 0); Δ′ there is less accurate." 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 || @@ -501,7 +505,6 @@ 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, 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 37f2cc3d7..2df34b25c 100644 --- a/test/runtests_grid_refinement.jl +++ b/test/runtests_grid_refinement.jl @@ -109,10 +109,20 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium @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) + # Close centers (e.g. two reverse-shear crossings near q_min) get brackets capped at half their gap, so each + # still sits alone in its own interval; the cut below the stencil floor warns. + b = @test_logs (:warn, r"capped") match_mode = :any GridRef.bracket_mandatory_nodes(grid, [0.5, 0.51], [0.01, 0.01], 1e-4) + @test all(diff(b) .> 0) + for c in (0.5, 0.51) + k = searchsortedlast(b, c) + @test b[k] < c < b[k+1] && isapprox(c - b[k], b[k+1] - c; rtol=1e-6) + end + @test searchsortedlast(b, 0.5) != searchsortedlast(b, 0.51) + # A near-tangent surface (huge or infinite stencil) is capped by the domain edges instead of clearing the grid. + for hw in (1e6, Inf) + b = @test_logs (:warn, r"capped") match_mode = :any GridRef.bracket_mandatory_nodes(grid, [0.3], [hw], 1e-4) + @test count(<(0.3), b) >= 3 && count(>(0.3), b) >= 3 + end # 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 From a97fe4aa2a921c91dfdd09c3f05a1fa959b3e2f0 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 17:43:31 -0400 Subject: [PATCH 12/16] Equilibrium - FEATURE - Add a hollow current profile to the large-aspect-ratio equilibrium for reversed shear Co-Authored-By: Claude Opus 5.5 --- .../Sources/Analytic/LargeAspectRatio.jl | 25 ++++++++++--------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/Equilibrium/Sources/Analytic/LargeAspectRatio.jl b/src/Equilibrium/Sources/Analytic/LargeAspectRatio.jl index 46fdc827f..20b43bcc4 100644 --- a/src/Equilibrium/Sources/Analytic/LargeAspectRatio.jl +++ b/src/Equilibrium/Sources/Analytic/LargeAspectRatio.jl @@ -11,7 +11,8 @@ A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equi - `q0`: The safety factor on axis. - `p_pres`: The exponent for the pressure profile, defined as `p00 * (1 - (r / a)^2)^p_pres`. - `p_sig`: The exponent that determines the shape of the current-related function profile. - - `sigma_type`: The type of sigma profile, can be "default" or "wesson". If "wesson", the sigma profile is defined as `sigma0 * (1 - (r / a)^2)^p_sig`. + - `sigma_type`: The type of sigma profile, can be "default", "wesson" or "hollow". If "wesson", the sigma profile is defined as `sigma0 * (1 - (r / a)^2)^p_sig`. If "hollow", it is `sigma0 * (1 + hollow * (r / a)^2) * (1 - (r / a)^2)^p_sig`, an off-axis current peak that makes q fall from `q0` to a minimum before rising (reversed shear) when `hollow` is large enough. + - `hollow`: The off-axis current enhancement of the "hollow" sigma profile (unused by the other profiles). - `mtheta`: The number of grid points in the poloidal direction. - `ma`: The number of grid points in the radial direction. - `zeroth`: If set to true, it neglects the Shafranov shift @@ -26,6 +27,7 @@ A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equi p_pres::Float64 = 2.0 p_sig::Float64 = 1.0 sigma_type::String = "default" + hollow::Float64 = 0.0 mtheta::Int = 128 ma::Int = 128 zeroth::Bool = false @@ -86,6 +88,14 @@ function lar_init_conditions(rmin::Float64, lar_input::LargeAspectRatioConfig) return r, y end +# Current profile σ(x), x = r/a, with σ(0) = sigma0 for every `sigma_type`. +function _lar_sigma(x::Float64, sigma0::Float64, lar_input::LargeAspectRatioConfig) + p_sig = lar_input.p_sig + lar_input.sigma_type == "wesson" && return sigma0 * (1 - x^2)^p_sig + lar_input.sigma_type == "hollow" && return sigma0 * (1 + lar_input.hollow * x^2) * (1 - x^2)^p_sig + return sigma0 / (1 + x^(2 * p_sig))^(1 + 1 / p_sig) +end + """ lar_der(dy, r, y, lar_input) @@ -111,19 +121,13 @@ function lar_der(dy::Vector{Float64}, r::Float64, y::Vector{Float64}, lar_input: beta0 = lar_input.beta0 p00 = beta0 / 2.0 p_pres = lar_input.p_pres - p_sig = lar_input.p_sig - sigma_type = lar_input.sigma_type x = r / lar_a xfac = 1 - x^2 pp = -2 * p_pres * p00 * x * xfac^(p_pres - 1) # p' sigma0 = 2.0 / (q0 * lar_r0) - sigma = if sigma_type == "wesson" - sigma0 * xfac^p_sig - else - sigma0 / (1 + x^(2 * p_sig))^(1 + 1 / p_sig) - end + sigma = _lar_sigma(x, sigma0, lar_input) # LAR equations have singularities at r=0; ensure we never evaluate there @assert r > 0.0 "LAR equations require r > 0 (called with r=$r)" @@ -159,7 +163,6 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf lar_r0 = lar_input.r0 q0 = lar_input.q0 beta0 = lar_input.beta0 - sigma_type = lar_input.sigma_type p00 = beta0 / 2.0 lar_input.p_pres = max(lar_input.p_pres, 1.001) @@ -194,9 +197,7 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf x = r / lar_a xfac = 1 - x^2 pval = p00 * xfac^lar_input.p_pres - sigma = (sigma_type == "wesson") ? - sigma0 * xfac^lar_input.p_sig : - sigma0 / (1 + x^(2 * lar_input.p_sig))^(1 + 1 / lar_input.p_sig) + sigma = _lar_sigma(x, sigma0, lar_input) q = r^2 * y[2] / (lar_r0 * y[1]) temp[i, 1] = r temp[i, 2] = y[1] From ca75806a9c452e32a2f4a6699385ca3970d3c1e9 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 18:28:33 -0400 Subject: [PATCH 13/16] Regression - TEST - Add a large-aspect-ratio reversed-shear n=1 case with a double q = 2 crossing Co-Authored-By: Claude Opus 5.5 --- examples/LAR_reversed_shear_example/gpec.toml | 44 ++++ .../cases/lar_reversed_shear_n1.toml | 193 ++++++++++++++++++ test/runtests_equil.jl | 11 + 3 files changed, 248 insertions(+) create mode 100644 examples/LAR_reversed_shear_example/gpec.toml create mode 100644 regression-harness/cases/lar_reversed_shear_n1.toml diff --git a/examples/LAR_reversed_shear_example/gpec.toml b/examples/LAR_reversed_shear_example/gpec.toml new file mode 100644 index 000000000..3bfca08bd --- /dev/null +++ b/examples/LAR_reversed_shear_example/gpec.toml @@ -0,0 +1,44 @@ +# Large-aspect-ratio equilibrium with reversed shear: n=1 ideal stability. +# A hollow current profile makes q fall from 2.35 on the axis to a minimum of 1.82 at ψ_N ≈ 0.37 +# and rise to 2.84 at the edge, so q = 2 is crossed twice, with q′ < 0 on the inner crossing and +# q′ > 0 on the outer one. Regression fixture for keeping both crossings as separate surfaces; the +# equilibrium is generated analytically from the embedded [LAR_INPUT] section, with no wall. + +[Equilibrium] +eq_type = "lar" # Type of the input 2D equilibrium file +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux +mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points + +[Wall] +shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +integrator = "forward" # Integration formalism: "forward" (dense ξ), "riccati" (chunked propagator BVP, singular-surface Δ′), or "galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Large-aspect-ratio analytic equilibrium parameters (eq_type = "lar"); see LargeAspectRatioConfig in src/Equilibrium. +[LAR_INPUT] +r0 = 10.0 # Major radius of the plasma [m] +a = 1.0 # Minor radius of the plasma [m] +beta0 = 1e-3 # Beta on axis (normalized pressure) +q0 = 2.4 # Safety factor on axis +p_pres = 2.0 # Exponent of the pressure profile p00·(1 − (r/a)²)^p_pres +p_sig = 2.0 # Exponent of the current-related profile +sigma_type = "hollow" # Current profile: "default", "wesson" or "hollow" (off-axis peak, reversed shear) +hollow = 6.0 # Off-axis current enhancement of the "hollow" profile +mtheta = 128 # Number of grid points in the poloidal direction +ma = 128 # Number of grid points in the radial direction diff --git a/regression-harness/cases/lar_reversed_shear_n1.toml b/regression-harness/cases/lar_reversed_shear_n1.toml new file mode 100644 index 000000000..7445392ca --- /dev/null +++ b/regression-harness/cases/lar_reversed_shear_n1.toml @@ -0,0 +1,193 @@ +# Regression case: large-aspect-ratio reversed-shear n=1 ideal stability. +# q falls from 2.35 on the axis to 1.82 and rises to 2.84 at the edge, so q = 2 is crossed twice; +# the two crossings must stay separate singular surfaces with opposite-sign q′. Tracks the +# free-boundary energies, the singular surfaces and the integration for that fixture. Each +# [quantities.*] block names an HDF5 path, how to extract it, and its noise floor. +[case] +name = "lar_reversed_shear_n1" +description = "Large-aspect-ratio reversed-shear equilibrium (hollow current), n=1, ideal stability" +example_dir = "examples/LAR_reversed_shear_example" + +# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian. +[quantities.et_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 10 + +[quantities.et_imag] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "imag_first" +label = "total energy Im(et[1])" +noise_threshold = 1e-10 +order = 11 + +[quantities.ep_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "real_first" +label = "plasma energy Re(ep[1])" +noise_threshold = 1e-10 +order = 12 + +[quantities.ev_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" +type = "complex_vector" +extract = "real_first" +label = "vacuum energy Re(ev[1])" +noise_threshold = 1e-10 +order = 13 + +[quantities.vacuum_eigenvalue] +h5path = "ForceFreeStates/FreeBoundaryStability/vacuum_eigenvalue" +type = "real_scalar" +extract = "value" +label = "vacuum matrix min eigenvalue" +noise_threshold = 1e-10 +order = 14 + +# Energies — full eigenvalue arrays +[quantities.ep_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "all_complex" +label = "plasma energy (all)" +noise_threshold = 1e-10 +order = 20 + +[quantities.ev_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" +type = "complex_vector" +extract = "all_complex" +label = "vacuum energy (all)" +noise_threshold = 1e-10 +order = 21 + +[quantities.et_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "all_complex" +label = "total energy (all)" +noise_threshold = 1e-10 +order = 22 + +# Integration +[quantities.nstep] +h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep" +type = "int_scalar" +extract = "value" +label = "ODE steps (saved)" +noise_threshold = 0 +order = 30 + +[quantities.nstep_total] +h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" +type = "int_scalar" +extract = "value" +label = "ODE steps (total)" +noise_threshold = 0 +order = 31 + +# Equilibrium +[quantities.q0] +h5path = "Equilibrium/q_axis" +type = "real_scalar" +extract = "value" +label = "q0" +noise_threshold = 1e-12 +order = 40 + +[quantities.q95] +h5path = "Equilibrium/q_95" +type = "real_scalar" +extract = "value" +label = "q95" +noise_threshold = 1e-12 +order = 41 + +[quantities.betat] +h5path = "Equilibrium/beta_t" +type = "real_scalar" +extract = "value" +label = "beta_t" +noise_threshold = 1e-12 +order = 42 + +[quantities.betan] +h5path = "Equilibrium/beta_N" +type = "real_scalar" +extract = "value" +label = "beta_n" +noise_threshold = 1e-12 +order = 43 + +# Singular surfaces +[quantities.msing] +h5path = "SingularSurfaces/rational_count" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 50 + +[quantities.sing_psi] +h5path = "SingularSurfaces/rational_psi" +type = "real_vector" +extract = "all_real" +label = "singular psi locations" +noise_threshold = 1e-8 +order = 51 + +[quantities.sing_q] +h5path = "SingularSurfaces/rational_q" +type = "real_vector" +extract = "all_real" +label = "singular q values" +noise_threshold = 1e-8 +order = 52 + +# Mode numbers +[quantities.mpert] +h5path = "Info/mpert" +type = "int_scalar" +extract = "value" +label = "mpert" +noise_threshold = 0 +order = 60 + +[quantities.npert] +h5path = "Info/npert" +type = "int_scalar" +extract = "value" +label = "npert" +noise_threshold = 0 +order = 61 + +# Profiles (checksums for fast comparison of large arrays) +[quantities.q_profile] +h5path = "Equilibrium/Profiles/q" +type = "real_vector" +extract = "checksum" +label = "q profile (checksum)" +noise_threshold = 0 +order = 70 + +[quantities.pressure_profile] +h5path = "Equilibrium/Profiles/mu0p" +type = "real_vector" +extract = "checksum" +label = "pressure profile (checksum)" +noise_threshold = 0 +order = 71 + +# Runtime (special: not from H5) +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index f8afac63c..efc0fc770 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -478,6 +478,17 @@ using TOML q_lar = plasma_eq_lar.profiles.q_spline.y @test all(>(0), q_lar) @test issorted(q_lar) + + # The "hollow" current profile keeps q(0) = q0 but drives an off-axis current peak, giving reversed shear: + # q falls to an interior minimum and rises again, crossing q = 2 on both legs. + hollow_input = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(; + r0=10.0, a=1.0, q0=2.4, beta0=1e-3, p_sig=2.0, sigma_type="hollow", hollow=6.0, mtheta=64, ma=64 + ) + q_rs = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(lar_config, hollow_input).profiles.q_spline.y + jmin = argmin(q_rs) + @test 1 < jmin < length(q_rs) + @test q_rs[jmin] < 2 < min(q_rs[1], q_rs[end]) + @test isapprox(q_rs[1], 2.4; rtol=0.05) end @testset "Solovev Equilibrium" begin From 3a899ff9d5b5244fa418f6a9b092c1b0341407d0 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 20:55:59 -0400 Subject: [PATCH 14/16] Regression - TEST - Give the axis q of the reversed-shear fixture as the run reports it Co-Authored-By: Claude Opus 5.5 --- examples/LAR_reversed_shear_example/gpec.toml | 2 +- regression-harness/cases/lar_reversed_shear_n1.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/LAR_reversed_shear_example/gpec.toml b/examples/LAR_reversed_shear_example/gpec.toml index 3bfca08bd..6387c228b 100644 --- a/examples/LAR_reversed_shear_example/gpec.toml +++ b/examples/LAR_reversed_shear_example/gpec.toml @@ -1,5 +1,5 @@ # Large-aspect-ratio equilibrium with reversed shear: n=1 ideal stability. -# A hollow current profile makes q fall from 2.35 on the axis to a minimum of 1.82 at ψ_N ≈ 0.37 +# A hollow current profile makes q fall from 2.40 on the axis to a minimum of 1.82 at ψ_N ≈ 0.37 # and rise to 2.84 at the edge, so q = 2 is crossed twice, with q′ < 0 on the inner crossing and # q′ > 0 on the outer one. Regression fixture for keeping both crossings as separate surfaces; the # equilibrium is generated analytically from the embedded [LAR_INPUT] section, with no wall. diff --git a/regression-harness/cases/lar_reversed_shear_n1.toml b/regression-harness/cases/lar_reversed_shear_n1.toml index 7445392ca..7111cf782 100644 --- a/regression-harness/cases/lar_reversed_shear_n1.toml +++ b/regression-harness/cases/lar_reversed_shear_n1.toml @@ -1,5 +1,5 @@ # Regression case: large-aspect-ratio reversed-shear n=1 ideal stability. -# q falls from 2.35 on the axis to 1.82 and rises to 2.84 at the edge, so q = 2 is crossed twice; +# q falls from 2.40 on the axis to 1.82 and rises to 2.84 at the edge, so q = 2 is crossed twice; # the two crossings must stay separate singular surfaces with opposite-sign q′. Tracks the # free-boundary energies, the singular surfaces and the integration for that fixture. Each # [quantities.*] block names an HDF5 path, how to extract it, and its noise floor. From e200765c2ff701f30e30f9ded7c6cd5524de0182 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 29 Oct 2026 14:34:41 -0400 Subject: [PATCH 15/16] ForceFreeStates - BUGFIX - Start every integrator at one qlow bound from the spline q extrema and filter surfaces by psi Replaces check_qlow_domain: qlow_start_psi (outermost crossing of q = qlow, bracketed on the spline q extrema sing_find! uses) is now the lower bound for forward, Riccati and Galerkin alike, and remove_singular_surfs! drops surfaces by psi against it. A qlow between q_min and q0 now excludes the reversed-shear ring and the core instead of starting at the axis or erroring. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/ForceFreeStates/CoreTypes.jl | 4 +- src/ForceFreeStates/EulerLagrange.jl | 21 +---- src/ForceFreeStates/Galerkin/GalerkinGrid.jl | 4 +- src/ForceFreeStates/Surfaces/Finding.jl | 82 +++++++------------- src/GeneralizedPerturbedEquilibrium.jl | 17 ++-- test/runtests_sing.jl | 19 +++-- 6 files changed, 54 insertions(+), 93 deletions(-) diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5d8565cce..72c23a87c 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -70,7 +70,7 @@ A mutable struct holding internal state variables for stability calculations. kinsing_scan_cond::Vector{Float64} = Float64[] kinsing_scan_threshold::Float64 = 0.0 psilim::Float64 = 0.0 - psilow::Float64 = 0.0 # lower integration bound; raised above the axis by sing_min! when qlow > qmin (RDCON gal) + psilow::Float64 = 0.0 # lower integration bound; the outermost crossing of q = qlow (sing_min!), else the first grid node qlim::Float64 = 0.0 q1lim::Float64 = 0.0 debug_settings::DebugSettings = DebugSettings() @@ -139,7 +139,7 @@ gpec.toml. - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (the drift-kinetic matrices of the KineticForces model, evaluated on the `[KineticForces]` profiles; needs `integrator = "forward"` and is deck-driven through `main`) - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components. - - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil + - `qlow::Float64` - Integration starts at the outermost crossing of q = qlow, excluding the q < qlow region inside it (a core, or a reversed-shear ring); no effect when q ≥ qlow everywhere - `dW_edge_scan_start::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [dW_edge_scan_start, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below). - `truncate_at_dW_peak::Bool` - When `true` and `dW_edge_scan_start < psilim`, the edge-dW scan's peak is adopted as the physical plasma edge, reproducing the ode_record_edge heuristic of Fortran STRIDE. `intr.psilim`/`intr.qlim`/`odet.u` are pulled back to the peak, fixups recorded beyond it are discarded, `q1lim` and the rational-surface list are re-derived at the new edge, and on the Riccati path the Δ' chunks are cut back to it. The result is the same physics as setting `psihigh` to the peak by hand. Δ' still depends on where in the edge band the peak falls, so set this deliberately when you mean to measure against the peak-defined edge (e.g. edge-mode regimes); leave at `false` (default) for the full domain at `qhigh` / `psihigh` / `rational_offset_fraction`. - `write_outputs_to_HDF5::Bool` - Write results to HDF5 format diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index b403ec97a..d1a337cc8 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -536,22 +536,8 @@ below it. function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, mats::MatrixSplines, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) - # Default psifac to minimum equilibrium psi value - odet.psifac = profiles.xs[1] - - # Use Newton iteration to find starting psi if qlow is above q0 - if ctrl.qlow > profiles.q_spline.y[1] - # Find last index where q < qlow - idx = findlast(jpsi -> profiles.q_spline.y[jpsi-1] < ctrl.qlow, 2:profiles.npts) - if idx !== nothing - odet.psifac = profiles.xs[idx] - end - odet.psifac = find_zero( - (psi -> profiles.q_spline(psi) - ctrl.qlow, - psi -> profiles.q_deriv(psi)), - odet.psifac, Roots.Newton() - ) - end + # Start at the qlow bound sing_min! set (`qlow_start_psi`), or the first grid node when it was not run. + odet.psifac = max(intr.psilow, profiles.xs[1]) # Find starting singular surface (where sing.psifac > psi(qlow/q0)) # Note: This logic is kept in initialize_el_at_axis! rather than chunk_el_integration_bounds @@ -712,7 +698,8 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro 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." + "singfac_min/|n q'| = $(sg.psifac - psi_end) reaches back past the current chunk start ψ=$psi_current. " * + "Move the domain start (psilow, or a qlow that lands on this surface) off it, or lower singfac_min if a neighboring crossing is too close." ) end @assert isempty(chunks) || psi_current >= chunks[end].psi_end "Overlapping chunks detected" diff --git a/src/ForceFreeStates/Galerkin/GalerkinGrid.jl b/src/ForceFreeStates/Galerkin/GalerkinGrid.jl index 37e1c6da8..41902df6b 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinGrid.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinGrid.jl @@ -50,8 +50,8 @@ end gal_resonant_surfaces(intr, equil) -> (sings, psilow, psihigh) Resonant surfaces inside the integration domain, in increasing ψ (Fortran `sing(1:msing)`), with the -domain bounds. `psilow` is raised above the axis by `sing_min!` when `qlow > qmin` (excludes the q 0 ? intr.psilow : equil.profiles.xs[1] diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index f05423ed1..c21e9ddd9 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -183,35 +183,13 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, end """ - check_qlow_domain(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) + remove_singular_surfs!(intr; psimin=-Inf, psimax=intr.psilim) -> keep -Error unless the region `q < qlow` on the ψ grid is empty or one interval starting at the axis. -On a reversed-shear profile with `q_min < qlow ≤ q0` it is instead a ring around q_min, and with -several minima it can be several intervals; no single lower bound excludes either: `sing_min!` -places `psilow` on the outermost leg, while `initialize_el_at_axis!` starts at the axis, and the -two disagree about the integration domain. -""" -function check_qlow_domain(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) - qy = equil.profiles.q_spline.y - below = qy .< ctrl.qlow - # A core: below on the first node and on every node up to the first one at or above qlow, then never again. - (!any(below) || (below[1] && issorted(below; rev=true))) && return nothing - error( - "qlow=$(ctrl.qlow): the region q < qlow is not a single core around the axis (q on the axis is $(qy[1]), " * - "the grid minimum is $(minimum(qy))); on a reversed-shear profile it is a ring around q_min, so no single " * - "lower integration bound excludes it. Set qlow below the q minimum to integrate from the axis, or above " * - "the axis value to start outside it." - ) -end - -""" - remove_singular_surfs!(intr; qmin=-Inf, psimax=intr.psilim) -> keep - -Drop rational surfaces with `q < qmin` or `psifac > psimax` from `intr.sing`, logging any dropped, +Drop rational surfaces outside `psimin ≤ psifac ≤ psimax` from `intr.sing`, logging any dropped, and return the indices of those kept. """ -function remove_singular_surfs!(intr::ForceFreeStatesInternal; qmin::Real=-Inf, psimax::Real=intr.psilim) - keep = [j for j in 1:intr.msing if intr.sing[j].q >= qmin && intr.sing[j].psifac <= psimax] +function remove_singular_surfs!(intr::ForceFreeStatesInternal; psimin::Real=-Inf, psimax::Real=intr.psilim) + keep = [j for j in 1:intr.msing if psimin <= intr.sing[j].psifac <= psimax] if length(keep) < intr.msing excluded = [(intr.sing[j].m, intr.sing[j].q) for j in setdiff(1:intr.msing, keep)] @info "Filtered $(length(excluded)) singular surface(s) outside integration domain: $(excluded)" @@ -221,39 +199,37 @@ function remove_singular_surfs!(intr::ForceFreeStatesInternal; qmin::Real=-Inf, return keep end +""" + qlow_start_psi(qlow, equil) -> ψ + +Lower integration bound set by `qlow`: the outermost crossing of q = qlow, so the whole q < qlow +region inside it is excluded, whether a monotonic core (q ≤ 1 sawtooth surfaces) or a reversed-shear +ring around q_min. Read from the same spline q extrema `sing_find!` brackets its roots with, so a dip +between grid nodes counts; the root is bracketed on the monotone stretch past the outermost extremum +below qlow. Returns the first grid node when q ≥ qlow everywhere. +""" +function qlow_start_psi(qlow::Real, equil::Equilibrium.PlasmaEquilibrium) + psiex, qex = equil.params.qextrema_psi, equil.params.qextrema_q + k = findlast(<(qlow), qex) + k === nothing && return equil.profiles.xs[1] + k < length(qex) || error("qlow=$qlow: q < qlow up to the edge of the ψ grid (q there is $(qex[end])), so no integration domain remains.") + return find_zero(psi -> equil.profiles.q_spline(psi) - qlow, (psiex[k], psiex[k+1]), Roots.Brent()) +end + """ sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) -Set the lower integration bound `intr.psilow`. Port of Fortran RDCON `sing_min` (sing.f): -when `qlow > qmin`, the q < qlow core (including any q ≤ 1 sawtooth/internal-kink surfaces) must be -excluded from the outer-region Galerkin domain — otherwise the Hermite FEM integrates through those -ideal singularities without imposing the ideal constraint, contaminating Δ′ at the innermost kept -surface. A Newton iteration locates ψ where q = qlow; scanning starts from the edge inward for -robustness in reverse-shear cores. When `qlow ≤ qmin` the axis value is kept. +Set the lower integration bound `intr.psilow` to `qlow_start_psi`, the one bound every +integrator starts from. Excluding the q < qlow region matters most for the outer-region Galerkin +solve, whose Hermite FEM would otherwise integrate through unhandled q ≤ 1 ideal singularities and +contaminate Δ′ at the innermost kept surface. """ function sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) profiles = equil.profiles - intr.psilow = profiles.xs[1] # default: equilibrium axis-side bound - ctrl.qlow > equil.params.qmin || return intr.psilow - - # Scan from the edge inward for the first node with q < qlow (robust for reverse-shear q). - qy = profiles.q_spline.y - jpsi = 1 - for j in (length(profiles.xs)-1):-1:1 - if qy[j] < ctrl.qlow - jpsi = j - break - end - end - - hint = Ref(jpsi) - intr.psilow = find_zero( - (psi -> profiles.q_spline(psi; hint=hint) - ctrl.qlow, - psi -> profiles.q_deriv(psi; hint=hint)), - profiles.xs[jpsi], Roots.Newton() - ) - @info "sing_min: qlow=$(@sprintf("%.3f", ctrl.qlow)) > qmin=$(@sprintf("%.3f", equil.params.qmin)); " * - "raising psilow from $(@sprintf("%.5f", profiles.xs[1])) to $(@sprintf("%.5f", intr.psilow)) (excludes q profiles.xs[1] && + @info "sing_min: qlow=$(@sprintf("%.3f", ctrl.qlow)); raising psilow from $(@sprintf("%.5f", profiles.xs[1])) " * + "to $(@sprintf("%.5f", intr.psilow)) (excludes the q < qlow region)" return intr.psilow end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 911b0808d..385e673fa 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -78,7 +78,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!, check_qlow_domain, sing_min!, sing_find!, remove_singular_surfs!, resist_eval_all!, resist_geometry, ResistGeometry +using .ForceFreeStates: sing_lim!, 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! @@ -657,7 +657,6 @@ function prepare_force_free_states!( # Determine psilim and qlim (where we will integrate to) sing_lim!(intr, ctrl, equil; psilim_cap=psilim_cap) - check_qlow_domain(ctrl, equil) # Fires whether or not the cap was applied: past the overlap point no surface retains a # well-separated inner region, so a domain reaching beyond it is worth flagging even when @@ -672,16 +671,10 @@ function prepare_force_free_states!( # Find all singular surfaces in the equilibrium sing_find!(intr, equil) - # Keep only surfaces inside the integration domain [qlow, qlim], so the Δ' BVP sees only crossable ones. - remove_singular_surfs!(intr; qmin=max(ctrl.qlow, equil.params.qmin)) - - # For the outer-region Galerkin solve, exclude the q < qlow core (incl. any q≤1 sawtooth - # surfaces) by raising psilow to where q = qlow (RDCON sing_min). Without this, the gal FEM - # integrates through the unhandled q≤1 ideal singularity and contaminates Δ′ at the innermost - # kept surface when q0 < qlow. No-op (keeps the axis bound) when qlow ≤ qmin. - if ctrl.integrator == "galerkin" - sing_min!(intr, ctrl, equil) - end + # Lower bound from qlow for every integrator, then keep only surfaces inside [psilow, psilim], + # so the Δ' BVP sees only crossable ones. + sing_min!(intr, ctrl, equil) + remove_singular_surfs!(intr; psimin=intr.psilow) # Populate Glasser-Greene-Johnson geometric coefficients (E, F, G, H, K, M) for each # surviving singular surface; `ggj_parameters` layers τ_A / τ_R on top. diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index 5b972d9b1..bd566bf95 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -202,14 +202,19 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex @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. + # qlow sets one lower bound, the outermost crossing of q = qlow, for every q_min < qlow: a qlow between q_min and q0 + # excludes the reversed-shear ring and the core inside it, and the surfaces there are dropped by ψ. 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 - # qlow equal to the axis value still leaves a ring below it. - @test_throws ErrorException FFS.check_qlow_domain(FFS.ForceFreeStatesControl(; qlow=q0), equil) + @test equil.params.qmin < 2.1 < q0 + @test FFS.qlow_start_psi(1.5, equil) == equil.profiles.xs[1] + @test isapprox(FFS.qlow_start_psi(2.1, equil), 0.4 + sqrt(0.3 / 3); atol=1e-6) + @test isapprox(FFS.qlow_start_psi(q0 + 0.1, equil), 0.4 + sqrt((q0 + 0.1 - 1.8) / 3); atol=1e-6) + intr_min = FFS.ForceFreeStatesInternal() + FFS.sing_min!(intr_min, FFS.ForceFreeStatesControl(; qlow=2.1), equil) + intr.psilim = 1.0 + FFS.remove_singular_surfs!(intr; psimin=intr_min.psilow) + @test intr.msing > 0 && all(s -> s.psifac >= intr_min.psilow, intr.sing) + @test !any(s -> isapprox(s.q, 2.0; atol=1e-8), intr.sing) # A qhigh below the axis value is crossed on both legs; the domain ends on the outer one. intr_lim = FFS.ForceFreeStatesInternal() From 29fdcaf2b896e7949ed5c0ad90ecdb2ad074b13a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 29 Oct 2026 14:34:31 -0400 Subject: [PATCH 16/16] Equilibrium - MINOR - Summarize capped grid brackets in one warning and pin dq/dpsi in the reversed-shear case Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- regression-harness/cases/lar_reversed_shear_n1.toml | 8 ++++++++ src/Equilibrium/Grid/GridRefinement.jl | 7 +++++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/regression-harness/cases/lar_reversed_shear_n1.toml b/regression-harness/cases/lar_reversed_shear_n1.toml index 7111cf782..45188af12 100644 --- a/regression-harness/cases/lar_reversed_shear_n1.toml +++ b/regression-harness/cases/lar_reversed_shear_n1.toml @@ -149,6 +149,14 @@ label = "singular q values" noise_threshold = 1e-8 order = 52 +[quantities.sing_dqdpsi] +h5path = "SingularSurfaces/dqdpsi" +type = "real_vector" +extract = "all_real" +label = "singular dq/dpsi (opposite signs on the two q = 2 crossings)" +noise_threshold = 1e-8 +order = 53 + # Mode numbers [quantities.mpert] h5path = "Info/mpert" diff --git a/src/Equilibrium/Grid/GridRefinement.jl b/src/Equilibrium/Grid/GridRefinement.jl index ac9ee1778..94cee6742 100644 --- a/src/Equilibrium/Grid/GridRefinement.jl +++ b/src/Equilibrium/Grid/GridRefinement.jl @@ -489,14 +489,14 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64} push!(hws, min_half_widths[idx]) end kept = Tuple{Float64,Bool}[(g, false) for g in grid] # (value, is_bracket_knot) + capped = Tuple{Float64,Float64,Float64}[] # (center, capped half-width, stencil floor) for (i, c) in enumerate(cs) k = clamp(searchsortedlast(grid, c), 1, length(grid) - 1) cap_nb = 0.5 * min(i > 1 ? c - cs[i-1] : Inf, i < length(cs) ? cs[i+1] - c : Inf) cap_edge = 0.5 * min(c - lo, hi - c) # Stencil floor capped by neighbors and edges; blend to half the local spacing capped by neighbors only. w = max(min(hws[i], cap_nb, cap_edge), min(0.5 * (grid[k+1] - grid[k]), cap_nb)) - w < hws[i] && @warn "Rational surface at ψ=$c: grid bracket half-width capped at $w, below its Δ′ stencil $(hws[i]) " * - "(neighboring surface or domain edge too close, or a near-tangent surface with q′ ≈ 0); Δ′ there is less accurate." + w < hws[i] && push!(capped, (c, w, hws[i])) left, right = c - w, c + w # Drop non-bracket, non-endpoint knots inside the zone or crowding a new bracket knot. filter!(kept) do t @@ -506,6 +506,9 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64} left > lo && push!(kept, (left, true)) right < hi && push!(kept, (right, true)) end + isempty(capped) || @warn "$(length(capped)) rational surface(s) have grid bracket half-widths capped below their Δ′ stencil " * + "(neighboring surface or domain edge too close, or a near-tangent surface with q′ ≈ 0); Δ′ there is less accurate. " * + "(ψ, half-width, stencil): $(join(capped, ", "))" sort!(kept; by=first) merged = Float64[] for (v, _) in kept