diff --git a/examples/LAR_reversed_shear_example/gpec.toml b/examples/LAR_reversed_shear_example/gpec.toml new file mode 100644 index 000000000..6387c228b --- /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.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. + +[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..45188af12 --- /dev/null +++ b/regression-harness/cases/lar_reversed_shear_n1.toml @@ -0,0 +1,201 @@ +# Regression case: large-aspect-ratio reversed-shear n=1 ideal stability. +# 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. +[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 + +[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" +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/src/Equilibrium/Grid/GridRefinement.jl b/src/Equilibrium/Grid/GridRefinement.jl index 9999dcd5f..94cee6742 100644 --- a/src/Equilibrium/Grid/GridRefinement.jl +++ b/src/Equilibrium/Grid/GridRefinement.jl @@ -469,21 +469,34 @@ 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. 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 = -Inf - 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 # duplicate rational (same q via several (m,n)) + (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) + 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) - 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] && 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 @@ -492,8 +505,10 @@ 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 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 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] diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 85ba6881f..b792a6bf2 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 ff8ed1ec5..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 @@ -661,7 +647,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 @@ -704,8 +693,15 @@ 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. " * + "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" push!(chunks, IntegrationChunk(; 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 1719efe76..c21e9ddd9 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 @@ -163,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) @@ -178,13 +183,13 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, end """ - remove_singular_surfs!(intr; qmin=-Inf, psimax=intr.psilim) -> keep + remove_singular_surfs!(intr; psimin=-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)" @@ -194,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 873f112e8..c791fd610 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -647,16 +647,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_equil.jl b/test/runtests_equil.jl index fc7b01980..0f13facdd 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -492,6 +492,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 diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index ef8b3146f..5ab027750 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -556,6 +556,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 ed0fd77bd..2df34b25c 100644 --- a/test/runtests_grid_refinement.jl +++ b/test/runtests_grid_refinement.jl @@ -108,6 +108,23 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium for c in (0.3, 0.7) @test !any(g -> abs(g - c) < 1e-12, b) end + + # 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 @testset "_validate_psi_nodes" begin diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index c5f17a049..bd566bf95 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,56 @@ 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 + + # 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.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() + 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 + # --------------------------------- # sing_lim # ---------------------------------