diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..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,14 +384,22 @@ 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 - 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 || @@ -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/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 345495fc6..815a80b27 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/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(; diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 18c45def7..bed1cf30c 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -42,27 +42,32 @@ 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 @@ -70,10 +75,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 - 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 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 27e29ed54..9e355ed43 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/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) diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 02fa2f07e..e6a95144e 100644 --- a/test/runtests_eulerlagrange.jl +++ b/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 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 diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index d58d11c2b..cf2ad7fe6 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()) @@ -172,6 +172,44 @@ 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. + # 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) + 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 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) + @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 + + # 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 + # --------------------------------- # sing_lim # ---------------------------------