From 833bd49a2580a92e7b91bc7344500a968370ea63 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 20 Aug 2026 20:38:38 -0400 Subject: [PATCH 1/6] InnerLayer.SLAYER - FEATURE! - Add the Fitzpatrick diffusive-resistive layer width Adds three algebraic scales to `LayerWidths` alongside the Riccati `delta_s`: `delta_norm` (r_s S^-1/3), `delta_visco` (the P^1/6 viscous-resistive broadening) and `delta_dr`, the diffusive-resistive width of Fitzpatrick, Nucl. Fusion 2025 (doi 10.1088/1741-4326/ae4fdd) Eq. (100), which sets the edge layer-overlap criterion. The mapping of Eq. (100) into SLAYER variables is not literal. The paper's tau_A (its Eq. 74) carries no shear, whereas SLAYER's tau_h divides by n*s, so lu = tau_R/tau_h = (n|s|) * (tau_R/tau_A). Substituting that cancels Eq. (100)'s explicit (n|s|)^(-1/2) exactly, leaving delta_dr = r_s^(3/2) * P_perp^(1/4) / ( lu^(1/2) * d_beta^(1/2) ) with no shear factor. Carrying the (n|s|)^(-1/2) on top of SLAYER's lu double-counts the shear and understates the width by sqrt(n|s|) -- a factor of 6.6 at q=7 on the DIII-D-like deck, enough to report no overlap anywhere. Validated against the paper's own Sect. 5.8 JET model, which this reproduces at Psi = 0.9985 for n = 1 (paper: 0.9985) and 0.9936 for n = 4 (paper: 0.9952). Marked `!` because `LayerWidths` gains fields. --- src/InnerLayer/SLAYER/LayerThickness.jl | 41 ++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl index 0b6db9761..9a70a2d38 100644 --- a/src/InnerLayer/SLAYER/LayerThickness.jl +++ b/src/InnerLayer/SLAYER/LayerThickness.jl @@ -137,6 +137,26 @@ is built from, retained as a drift-scale reference. - `delta_s` -- complex layer thickness `δ_s = dels_db · d_β` [m] - `delta_s_m` -- `|δ_s|`, the resistive layer thickness in meters (primary) - `d_beta` -- β-weighted ion scale `c_β·d_i` in meters (drift reference) + - `delta_norm` -- layer normalization length `r_s · S^(-1/3)` in meters, the + length by which Δ' is made dimensionless (`Δ̂' = Δ'·delta_norm`). Shared by all + four regimes of Burgess et al. (2026) Eqs. (10)-(13), not specific to any one + of them, and **not** the classic non-rotating Furth-Killeen-Rosenbluth (1963) + constant-ψ tearing width, which scales as `S^(-2/5)` and is not computed here. + Exactly the reciprocal of the `delta_n` Δ-normalization already carried on + `SLAYERParameters`, restated as a length for comparison with `delta_s_m` + - `delta_dr` -- diffusive-resistive layer thickness in meters, Fitzpatrick (2025) + Eq. (100). Strong magnetic shear near the separatrix forces every resonant layer in + that region into the diffusive-resistive (`nu = 1/4`) regime, so this is the width the + edge overlap criterion compares against surface spacing. In the paper's normalized + radius, `delta_hat = tau_A^(1/2) / (tau_R^(1/4) tau_E^(1/4) d_beta_hat^(1/2) (n|s|)^(1/2))`; + with `lu = tau_R/tau_A` and `P_perp = tau_R/tau_E` that is + `lu^(-1/2) P_perp^(1/4) / (d_beta_hat^(1/2) (n|s|)^(1/2))`, scaled by `rs` for meters. + Distinct from `delta_visco`, which is the 2/3 power of the same timescale grouping and + carries neither the `d_beta` nor the shear factor + - `delta_visco` -- viscous-resistive scale `delta_norm · P_perp^(1/6)` in meters, + from the `P^(1/6)` broadening of the VR-regime growth rate, Burgess et al. + (2026) Eq. (11). The paper gives the growth rate rather than an explicit width; + this is the corresponding length `delta_s_m` should sit within a few orders of magnitude of `d_beta` for a well-posed surface (`dels_db` is O(1)); a large gap flags a normalisation @@ -150,6 +170,9 @@ struct LayerWidths delta_s::ComplexF64 delta_s_m::Float64 d_beta::Float64 + delta_norm::Float64 + delta_visco::Float64 + delta_dr::Float64 end """ @@ -161,11 +184,27 @@ surface. Runs [`riccati_del_s`](@ref) for the dimensionless `δ_s / d_β` and scales by `p.d_beta` to obtain `δ_s` in meters. Keyword arguments are forwarded to `riccati_del_s`. + +The two algebraic comparison scales `delta_norm` and `delta_visco` come from +`p` directly and do not depend on the Riccati solve. """ function slayer_layer_thickness(p::SLAYERParameters; kwargs...) dels_db = riccati_del_s(p; kwargs...) delta_s = dels_db * p.d_beta + # Layer normalization length: exactly 1/delta_n (delta_n = S^(1/3)/r_s), computed + # from rs and lu so it reads as a length. Burgess et al. (2026), the Δ̂' = Δ'/S^(1/3) + # normalization preceding Eq. (10). + delta_norm = p.rs * p.lu^(-1.0 / 3.0) + # Viscous-resistive broadening, P^(1/6) coefficient of Burgess et al. (2026) Eq. (11). + delta_visco = delta_norm * p.P_perp^(1.0 / 6.0) + # Diffusive-resistive width, Fitzpatrick (2025) Eq. (100), in meters. d_beta_hat is the + # normalized ion scale d_beta/rs; the shear is the r-based |s| SLAYER already carries. + # The paper's tau_A (its Eq. 74) carries no shear, whereas SLAYER's tau_h divides by n*s, + # so lu = tau_R/tau_h = (n|s|) * (tau_R/tau_A). Substituting that into Eq. (100) cancels its + # explicit (n|s|)^(-1/2) exactly, leaving no shear dependence in terms of lu. + d_beta_hat = p.d_beta / p.rs + delta_dr = d_beta_hat > 0 ? p.rs * p.lu^(-0.5) * p.P_perp^(0.25) / sqrt(d_beta_hat) : NaN return LayerWidths(p.ising, p.m, p.n, dels_db, delta_s, abs(delta_s), - p.d_beta) + p.d_beta, delta_norm, delta_visco, delta_dr) end From 960dc38fc9165edc90ef7a90bf0858d415acc3f5 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 20 Aug 2026 20:38:55 -0400 Subject: [PATCH 2/6] Equilibrium - FEATURE! - Share the separatrix edge q-law with the layer-overlap scan Factors the diverging edge model q ~ -A*ln(1-psi) out of the grid-refinement density floor into `edge_q_law`, so one fit serves both that floor and the resistive-layer overlap scan's search for rational surfaces beyond the equilibrium grid. A cubic in psi (on q or on iota) saturates instead of diverging and is not usable for the extrapolation; measured on the DIII-D-like deck, extrapolating from a grid ending at psi=0.97 to the true q=8 surface at psi=0.99976, the log law lands within 5% while a cubic iota undershoots by 18%. `edge_q_law` returns nothing when the diverging model does not describe the equilibrium, judged by knot count, sign of A, and the log fit beating a plain linear one. That test now gates the edge density floor, which was previously applied unconditionally: a limited plasma with finite edge q is no longer packed as if q blew up. Marked `!` because that changes the auto grid, and so results, for limited decks (Solovev, LAR); diverted decks are unaffected. Also adds an optional `psihigh` to `refined_psi_grid`, which truncates the measured knot density so a pass-1 equilibrium can supply the density for a reduced domain. --- src/Equilibrium/GridRefinement.jl | 123 +++++++++++++++++++++++++++++- 1 file changed, 121 insertions(+), 2 deletions(-) diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..87aee692c 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -57,6 +57,86 @@ const CORE_MODEL_PSI_MAX = 0.03 const EDGE_MODEL_PSI_MIN = 0.9 # θ-lines subsample stride for the 2D geometry channels const THETA_STRIDE = 8 +# --- shared separatrix edge q-law ------------------------------------------------------------- +# Minimum knots in the edge band before a fit is attempted. +const EDGE_FIT_MIN_KNOTS = 4 +# The diverging model must explain the edge q this well in absolute terms. Measured over the +# shipped decks: DIII-D-like 0.9957 and 0.9989 (diverted, q -> inf at the separatrix) against +# 0.972 for the a10 fixed-boundary case, 0.9035 for LAR and 0.7600 for Solovev (all limited, +# finite edge q). Rejecting is the safe direction -- it only means no extrapolation. +const EDGE_FIT_MIN_R2 = 0.99 + +# Least-squares slope and coefficient of determination for y = a + b*x. +function _linfit_r2(x::Vector{Float64}, y::Vector{Float64}) + x_bar = sum(x) / length(x) + y_bar = sum(y) / length(y) + sxx = sum((x .- x_bar) .^ 2) + sxx > 0 || return (NaN, NaN, NaN, NaN) + b = sum((x .- x_bar) .* (y .- y_bar)) / sxx + ss_res = sum((y .- (y_bar .+ b .* (x .- x_bar))) .^ 2) + ss_tot = sum((y .- y_bar) .^ 2) + r2 = ss_tot > 0 ? 1 - ss_res / ss_tot : NaN + return (b, r2, x_bar, y_bar) +end + +""" + edge_q_law(equil; psi_max, psi_min=EDGE_MODEL_PSI_MIN, min_knots=EDGE_FIT_MIN_KNOTS, + min_r2=EDGE_FIT_MIN_R2) -> nothing | (; A, q_bar, u_bar, n_knots, r2_log, r2_linear) + +Least-squares fit of the separatrix edge law `q = q̄ + A·(ln(1−ψ) − ū)` over the equilibrium's +outer knots — the single shared statement of that model, used both by the grid-refinement edge +density floor and by the resistive-layer overlap scan's out-of-grid surface search. + +Returns `nothing` when the diverging model does **not** describe this equilibrium's edge, so a +plasma with finite edge q is never extrapolated as if q blew up: + + - fewer than `min_knots` knots in the band; + - `A ≥ 0`, i.e. q not rising toward ψ = 1; + - `r2_log < min_r2` — the log law does not actually fit; + - `r2_log ≤ r2_linear` — a plain linear-in-ψ fit explains the edge q at least as well, which is + what a **limited** plasma looks like. This comparison carries no scale and is what separates + the shipped limited decks (Solovev 0.760 vs 0.9996 linear; LAR 0.904 vs 0.998) from the + diverted ones (DIII-D 0.996 vs 0.865, 0.999 vs 0.594). + +This is a test of the **model**, not a topology classification: it asks whether q diverges +logarithmically here, not whether an x-point exists. Geometric x-point detection is +[`classify_topology`](@ref), which is a separate concern. +""" +function edge_q_law(equil::PlasmaEquilibrium; + psi_max::Real=Float64(equil.profiles.xs[end]), + psi_min::Real=EDGE_MODEL_PSI_MIN, + min_knots::Int=EDGE_FIT_MIN_KNOTS, + min_r2::Real=EDGE_FIT_MIN_R2) + xs = collect(Float64, equil.profiles.xs) + band = findall(x -> x >= psi_min && x < psi_max, xs) + if length(band) < min_knots + n_tail = max(min_knots, length(xs) ÷ 10) + band = filter(i -> xs[i] < psi_max, collect(max(1, length(xs) - n_tail + 1):length(xs))) + end + length(band) >= min_knots || return nothing + + q = [Float64(equil.profiles.q_spline(xs[i])) for i in band] + u = [log(1.0 - xs[i]) for i in band] + all(isfinite, u) && all(isfinite, q) || return nothing + + A, r2_log, u_bar, q_bar = _linfit_r2(u, q) + _, r2_linear, _, _ = _linfit_r2([xs[i] for i in band], q) + (isfinite(A) && A < 0) || return nothing # q must rise toward the edge + (isfinite(r2_log) && r2_log >= min_r2) || return nothing + (isfinite(r2_linear) && r2_log > r2_linear) || return nothing + + return (A=A, q_bar=q_bar, u_bar=u_bar, n_knots=length(band), r2_log=r2_log, r2_linear=r2_linear) +end + +""" +ψ at which the edge law reaches `q_target`; closed form, no root-finding needed. +""" +edge_q_law_psi(fit, q_target::Real) = 1.0 - exp(fit.u_bar + (q_target - fit.q_bar) / fit.A) + +""" +dq/dψ from the edge law: q = q̄ + A·ln(1−ψ) + const ⇒ dq/dψ = −A/(1−ψ). +""" +edge_q_law_dqdpsi(fit, psi::Real) = -fit.A / (1.0 - psi) # Rational-surface bracketing (Δ′ robustness). The ideal-MHD Δ′ asymptotic matching samples the # cubic equilibrium splines' 2nd/3rd derivatives across each rational ψ_s over the matching stencil # [ψ_s − dpsi, ψ_s + dpsi], dpsi = singfac_min/|n·q′|. A cubic 3rd derivative is piecewise constant @@ -248,10 +328,15 @@ function _knot_density(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothin # nodal data of the smallest flux surfaces is dominated by integration and axis # extrapolation error, so measured curvature is not trusted below the core split. dlog = (4.0 * tau)^(1 / 3) + # The edge floor encodes the DIVERGING edge law q ≈ -A·ln(1-ψ), so it is applied only where + # that model actually describes the equilibrium. A limited plasma has finite edge q and must + # not be packed as if q blew up. The density itself stays A-independent (that is the point of + # the form: uniform relative q′ error regardless of A) -- the fit supplies validity, not slope. + edge_diverges = edge_q_law(equil) !== nothing @inbounds for i in 1:n if xs[i] <= CORE_MODEL_PSI_MAX rho_s[i] = 1.0 / (dlog * xs[i]) - elseif xs[i] >= EDGE_MODEL_PSI_MIN + elseif edge_diverges && xs[i] >= EDGE_MODEL_PSI_MIN rho_s[i] = max(rho_s[i], 1.0 / (dlog * (1.0 - xs[i]))) end rho_s[i] = max(rho_s[i], 1.0 / H_TARGET_MAX) @@ -409,9 +494,33 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64} return merged end +""" + _truncate_density(xs, rho, psihigh) -> (xs_t, rho_t) + +Restrict a measured knot density to `[xs[1], psihigh]`, linearly interpolating `rho` at the new +outer endpoint so the density integral stays continuous in `psihigh`. Errors when `psihigh` lies +outside the sampled grid, where the density is unmeasured. +""" +function _truncate_density(xs, rho::Vector{Float64}, psihigh::Float64) + xs_v = collect(Float64, xs) + psihigh > xs_v[1] || + error("_truncate_density: psihigh=$psihigh must exceed the inner grid bound $(xs_v[1])") + psihigh <= xs_v[end] + 1e-12 || + error( + "_truncate_density: psihigh=$psihigh exceeds the pass-1 grid end $(xs_v[end]); " * + "the knot density is unmeasured there — form the enlarged domain first, then refine against it" + ) + psihigh >= xs_v[end] - 1e-12 && return (xs_v, rho) + + k = searchsortedlast(xs_v, psihigh) + frac = (psihigh - xs_v[k]) / (xs_v[k+1] - xs_v[k]) + rho_end = rho[k] + frac * (rho[k+1] - rho[k]) + return (vcat(xs_v[1:k], psihigh), vcat(rho[1:k], rho_end)) +end + """ refined_psi_grid(equil::PlasmaEquilibrium; tau, kin=nothing, mandatory=Float64[], - singfac_min=1e-4, n_min=1, bracket_coef=BRACKET_COEF, + psihigh=nothing, singfac_min=1e-4, n_min=1, bracket_coef=BRACKET_COEF, min_spacing=MIN_KNOT_SPACING, N_cap=1024) -> Vector{Float64} Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curvature knot @@ -423,11 +532,18 @@ whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ va `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the Δ′ extraction's cubic 3rd derivative jump mid-stencil (see `BRACKET_COEF`). + +`psihigh` builds the grid for a domain *smaller* than the one `equil` was formed on: the measured +density is truncated there and the last node lands exactly on it. This lets a pass-1 equilibrium +supply the density for a re-form on a reduced domain without an extra solve. Passing a `psihigh` +beyond the pass-1 grid is an error — the density out there is unmeasured, so an enlarged domain +must be formed first and refined against that. """ function refined_psi_grid(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, mandatory::Vector{Float64}=Float64[], + psihigh::Union{Nothing,Real}=nothing, singfac_min::Float64=1e-4, n_min::Int=1, bracket_coef::Float64=BRACKET_COEF, @@ -435,6 +551,9 @@ function refined_psi_grid(equil::PlasmaEquilibrium; N_cap::Int=REFINED_N_CAP) xs = equil.profiles.xs rho = _knot_density(equil; tau, kin) + if psihigh !== nothing + xs, rho = _truncate_density(xs, rho, Float64(psihigh)) + end # Floor the density to a fixed (τ-independent) locally-uniform fine patch around each rational # so Δ′ has the resolution to sample 3rd derivatives there at any accuracy target (see # RATIONAL_RES_SPACING). From 73c31399c62ae1a2a2825b59e3afd7e2800d9605 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 20 Aug 2026 20:38:55 -0400 Subject: [PATCH 3/6] Tearing - FEATURE - Add the resistive-layer overlap scan Locates the q = m/n rational surfaces, computes each one's resistive layer width, and reports the outermost psi at which adjacent layers are still separated -- the outermost domain the matched-asymptotic treatment supports, per Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9. Surfaces beyond the equilibrium grid are located on the shared separatrix edge q-law, so the scan can report an overlap point lying outside psihigh; that is recorded rather than clipped, since it is the bookkeeping a user needs to see whether the domain was constrained by layer physics at all. Widths are compared in normalized flux, converted from metres by the Jacobian of whichever radial label `rs_method` selects -- the same label that sets the shear, tau_R and d_beta/r, because Eq. (100) is not covariant (delta ~ J^(1/2) under a change of radial variable) and is only valid in the coordinate it was derived in. On the DIII-D-like deck the label is worth one rational surface: the cut sits at psi=0.99752 under `:midplane` and 0.99939 under `:flux`, Fitzpatrick's Eq. (30) label. Reports both the |delta_s| and Eq. (100) channels separately so the two can be compared; they disagree by about two surfaces, which is the subject of issue #411. --- src/Tearing/LayerOverlap.jl | 385 ++++++++++++++++++++++++++++++++++++ src/Tearing/Tearing.jl | 2 + 2 files changed, 387 insertions(+) create mode 100644 src/Tearing/LayerOverlap.jl diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl new file mode 100644 index 000000000..1e67f7e5e --- /dev/null +++ b/src/Tearing/LayerOverlap.jl @@ -0,0 +1,385 @@ +# LayerOverlap.jl +# +# Resistive-layer overlap scan: how far out in ψ the equilibrium domain needs +# to extend before adjacent rational surfaces' resistive layers run into each +# other. Once two neighbouring layers overlap, neither surface has a +# well-separated inner region, so the matched-asymptotic treatment stops being +# meaningful there and extending `psihigh` past that point buys nothing. +# Criterion: Fitzpatrick, Nucl. Fusion 2025 (doi 10.1088/1741-4326/ae4fdd) Sect. 5.9 — +# retain surfaces with psi < 1 - eps_c, where overlap is judged with the +# diffusive-resistive layer width of its Eq. (100), the regime strong shear forces +# near the separatrix. +# +# Lives in `Tearing` rather than `InnerLayer/SLAYER` for the same reason +# `build_ggj_inputs` does: it needs `ForceFreeStates._find_rational_surfaces`, +# and `InnerLayer` loads before `ForceFreeStates`. +# +# Surfaces beyond the equilibrium's own ψ grid are located on the separatrix edge +# law q ≈ -A·ln(1-ψ), via the shared `Equilibrium.edge_q_law` helper that also gates +# the grid-refinement edge density model -- one statement of the model, one fit, one +# limited-plasma guard. The log form is the only one that survives the +# extrapolation: a cubic in ψ (on q or on ι = 1/q) saturates instead of diverging. Measured on the DIII-D-like deck, extrapolating +# from a grid ending at ψ=0.97 to the true q=8 surface at ψ=0.99976, the edge law +# lands within 5% while a cubic ι undershoots by 18% — enough to miss the surface +# entirely and report no overlap where there is one. +# +# The extrapolated surfaces inform the *choice* of domain only; the final +# equilibrium is always re-formed and solved on the accepted domain. + +using ..Utilities: KineticProfiles +using ..ForceFreeStates: _find_rational_surfaces, SingType +using ..InnerLayer.SLAYER: SLAYERParameters, build_slayer_inputs, slayer_layer_thickness, + surface_minor_radius, surface_da_dpsi, radial_label +using ..Utilities.NeoclassicalResistivity: NeoResistivityModel, SauterNeoModel +using ..Equilibrium: edge_q_law, edge_q_law_psi, edge_q_law_dqdpsi, InverseIngest +using FastInterpolations: cubic_interp, ExtendExtrap, DerivOp +using Roots: find_zero, Brent +using Printf: @sprintf + +""" + LayerOverlapScan + +Per-surface resistive layer widths and the `psihigh` they imply, from +[`resistive_layer_overlap`](@ref). + +All widths are in **normalized flux**, converted from the metre-valued +`LayerWidths` scales by `Δψ = w / |da/dψ|` so they can be compared against +surface *separations* in ψ. Comparing the metre values directly against ψ is a +dimensional error. + +# Fields + + - `m`, `n` -- resonant mode numbers per surface, ordered by increasing ψ + - `psi` -- surface location in normalized flux + - `rs` -- minor radius at the surface in meters + - `delta_s_m` -- `|δ_s|`, the Riccati resistive layer width, in meters + - `width_delta_s` -- the same width in normalized flux, `|δ_s|/|da/dψ|` + - `width_visco` -- the visco-resistive comparison scale as a Δψ width + - `width_dr` -- the diffusive-resistive width of Fitzpatrick (2025) Eq. (100) as a Δψ + width. This is the criterion the recommendation uses: the paper's Sect. 5.6 shows that + strong shear near the separatrix forces every layer there into the DR regime + - `extrapolated` -- true when the surface was located outside the + equilibrium's ψ grid via the edge q-law + - `psihigh_delta_s`, `psihigh_visco` -- domain implied by each width channel, + `nothing` when that channel never overlaps within the scanned range + - `psihigh` -- the recommendation: the smaller of the two, or `nothing` + - `first_overlap` -- index of the first surface that overlaps its inner + neighbour, `nothing` when none do + - `notes` -- surfaces that were located but could not be scored, and why +""" +struct LayerOverlapScan + m::Vector{Int} + n::Vector{Int} + psi::Vector{Float64} + rs::Vector{Float64} + delta_s_m::Vector{Float64} + width_delta_s::Vector{Float64} + width_visco::Vector{Float64} + width_dr::Vector{Float64} + extrapolated::Vector{Bool} + psihigh_delta_s::Union{Float64,Nothing} + psihigh_visco::Union{Float64,Nothing} + psihigh_dr::Union{Float64,Nothing} + psihigh::Union{Float64,Nothing} + first_overlap::Union{Int,Nothing} + notes::Vector{String} +end + +# Inverse equilibria (CHEASE) are solved on a PRESCRIBED boundary: `InverseIngest.sq_xs` spans +# [0, 1] and `sq_fs[:, 3]` is the code's own q there, finite at ψ = 1 (6.90 on the shipped +# fixture). So beyond `psihigh` there is nothing to model -- the real q is already known, and the +# scan reads it instead of extrapolating. That also means the search runs all the way to ψ = 1 +# rather than stopping short of a separatrix that a fixed-boundary equilibrium does not have. +# +# The direct/EFIT path cannot do this: a g-file also carries q on [0, 1], but it is the +# reconstruction's q, which goes to a finite qa where a diverted plasma's q diverges. GPEC +# recomputes q by field-line tracing out to psihigh, and past that nothing has been traced. +_inverse_q_spline(ingest::InverseIngest) = + cubic_interp(collect(Float64, ingest.sq_xs), collect(Float64, ingest.sq_fs[:, 3]); + extrap=ExtendExtrap()) + +# Locate q = q_target on the real q profile, between `psi_lo` and `psi_hi`. Returns nothing when +# the target is outside the range this equilibrium actually reaches. +function _real_q_surface(qspl, q_target::Real, psi_lo::Float64, psi_hi::Float64) + q_lo, q_hi = Float64(qspl(psi_lo)), Float64(qspl(psi_hi)) + (q_lo - q_target) * (q_hi - q_target) <= 0 || return nothing + return find_zero(p -> Float64(qspl(p)) - q_target, (psi_lo, psi_hi), Brent()) +end + +""" + resistive_layer_overlap(equil, profiles::KineticProfiles; n_tor, kwargs...) + -> LayerOverlapScan + +Locate the q = m/`n_tor` rational surfaces, compute each one's resistive layer +width, and report the outermost `psihigh` at which adjacent layers are still +separated. + +Surfaces inside the equilibrium grid come from +`ForceFreeStates._find_rational_surfaces` (Brent root-finding segmented between +q-extrema, so reverse shear is handled). Surfaces beyond the grid come from the +separatrix edge law `q ≈ -A·ln(1-ψ)` fitted over the outer knots when +`extrapolate=true`, which is what lets the scan recommend a domain *larger* than +the one it was handed. See the file header for why a cubic extrapolation is not +usable here. + +Layer widths come from [`slayer_layer_thickness`](@ref) via +[`build_slayer_inputs`](@ref), so the scan uses the same plasma inputs and the +same resistivity closure as the SLAYER analysis itself. + +# Arguments + + - `equil` -- a formed `PlasmaEquilibrium` + - `profiles` -- `KineticProfiles` spanning the scanned ψ range + +# Keyword arguments + + - `n_tor` -- toroidal mode number; surfaces are q = m/`n_tor` + + - `m_max` -- runaway backstop on the poloidal mode number (default `2000`). This is **not** + the physics limit: the outward search terminates when the edge law places a surface at or + beyond `psi_cap`, i.e. when it is indistinguishable from the separatrix. On a diverted + equilibrium rational surfaces accumulate without bound toward ψ = 1, so some cap is needed, + but it should never be what stops the scan — if it is, the scan says so in `notes`, because + "no overlap found" would otherwise be reported for a non-physical reason + - `psihigh_safe` -- treat this as the outer edge of trusted equilibrium data + (default `equil.profiles.xs[end]`) + - `psi_cap` -- never place a surface beyond this ψ on the **direct** path (default + `1 - 1e-9`). Purely numerical: it is not derived from the contour clamp or from any topology + test, and it stops the outward search from crowding indefinitely onto ψ = 1 where surfaces + become indistinguishable. Loosening it is cheap (1e-6 → 1e-9 adds ~5 surfaces on the + DIII-D-like deck at no measurable cost and does not change the answer) and guards against + missing an overlap that lies further out than the cap. Inverse equilibria ignore it: they + have a prescribed boundary and the search runs to ψ = 1 on the real q. + + Note the accuracy limit here is **not** `psi_cap` but `surface_da_dpsi`, which clamps its + finite-difference stencil to `1 - 1e-4`; surfaces closer than that to the boundary get + geometry evaluated at `1 - 1e-4` rather than at their own ψ (da/dψ = 0.348729 identically + for every ψ beyond it on the DIII-D-like deck). Their widths are therefore approximate + - `extrapolate` -- search past `psihigh_safe` on the edge q-law (default `true`) + - `theta` -- poloidal angle for the minor-radius chord (default `0.0`) + - `mu_i`, `zeff`, `chi_perp`, `chi_tor`, `resistivity_model`, `lnLambda_form` + -- passed through to `build_slayer_inputs` + +The overlap criterion is that surface `k` overlaps its inner neighbour when +`ψ_k − w_k/2 < ψ_{k-1} + w_{k-1}/2`. When that happens **both** surfaces are +contaminated, so the recommendation cuts at the inner edge of surface `k-1`, +leaving `k-2` as the outermost trustworthy surface. +""" +function resistive_layer_overlap(equil, profiles::KineticProfiles; + n_tor::Integer, + m_max::Integer=2000, + max_layer_solves::Integer=64, + psihigh_safe::Real=Float64(equil.profiles.xs[end]), + psi_cap::Real=1.0 - 1e-9, + extrapolate::Bool=true, + theta::Real=0.0, + rs_method::Symbol=:midplane, + mu_i::Real=2.0, + zeff::Real=1.0, + chi_perp=1.0, + chi_tor=1.0, + resistivity_model::NeoResistivityModel=SauterNeoModel(), + lnLambda_form::Symbol=:nrl) + + n_tor > 0 || throw(ArgumentError("resistive_layer_overlap: n_tor must be positive, got $n_tor")) + # Widths come back in metres of whichever radial label `rs_method` selects, so the + # conversion into normalised flux must use that same label's Jacobian. + _, _dr_dpsi = radial_label(equil; rs_method=rs_method, theta=theta) + psihigh_safe = Float64(psihigh_safe) + psi_cap = Float64(psi_cap) + notes = String[] + + # In-grid surfaces for this n, ordered outward. + found = [(m=s.m, psi=s.psifac, extrap=false) + for s in _find_rational_surfaces(equil, Int(n_tor), Int(n_tor)) + if s.psifac <= psihigh_safe && s.m <= m_max] + sort!(found; by=s -> s.psi) + + # Surfaces past the grid. An inverse equilibrium already knows its real q out to ψ = 1, so + # it is read rather than modelled; only the direct path needs the edge law. + edge_fit = nothing + inv_q = (getfield(equil, :ingest) isa InverseIngest) ? _inverse_q_spline(equil.ingest) : nothing + if extrapolate && inv_q !== nothing + psi_top = 1.0 # a prescribed boundary, not a separatrix + m_start = isempty(found) ? 1 : maximum(s.m for s in found) + 1 + reached_end = false + for m in m_start:Int(m_max) + psi_m = _real_q_surface(inv_q, m / n_tor, psihigh_safe, psi_top) + if psi_m === nothing + reached_end = true # q never reaches m/n inside this equilibrium + break + end + psi_m > psihigh_safe || continue + push!(found, (m=m, psi=psi_m, extrap=true)) + end + !reached_end && m_start <= m_max && + push!(notes, + "outward search stopped at the m_max = $m_max backstop rather than at the boundary " * + "q; a \"no overlap\" result here is not conclusive -- raise m_max") + push!(notes, "inverse equilibrium: surfaces beyond ψ = $(@sprintf("%.6f", psihigh_safe)) " * + "read from the real q profile out to ψ = 1 (no edge-law extrapolation)") + elseif extrapolate && psihigh_safe < psi_cap + edge_fit = edge_q_law(equil; psi_max=psihigh_safe) + if edge_fit === nothing + push!( + notes, + "edge q-law rejected: the diverging model q ~ -A*ln(1-ψ) does not describe " * + "this edge (limited plasma with finite edge q, q not rising, or too few outer " * + "knots), so no surfaces are placed beyond ψ = $(@sprintf("%.6f", psihigh_safe))" + ) + else + m_start = isempty(found) ? 1 : maximum(s.m for s in found) + 1 + # Terminates on psi_cap -- surfaces indistinguishable from the separatrix -- not on + # m_max. Reaching m_max means the scan was cut short for a non-physical reason and + # any "no overlap" verdict from it is unsafe, so record that. + reached_cap = false + for m in m_start:Int(m_max) + psi_m = edge_q_law_psi(edge_fit, m / n_tor) + psi_m > psihigh_safe || continue # already inside the grid + if psi_m >= psi_cap + reached_cap = true # separatrix reached: the physical end of the scan + break + end + push!(found, (m=m, psi=psi_m, extrap=true)) + end + if !reached_cap && m_start <= m_max + push!(notes, + "outward search stopped at the m_max = $m_max backstop rather than at the " * + "separatrix; a \"no overlap\" result here is not conclusive -- raise m_max") + end + end + end + + isempty(found) && return LayerOverlapScan(Int[], Int[], Float64[], Float64[], Float64[], Float64[], + Float64[], Bool[], nothing, nothing, nothing, nothing, + push!(notes, "no q = m/$n_tor surfaces found")) + + ms, ns, psis, rss = Int[], Int[], Float64[], Float64[] + dels_m, w_dels, w_visc, w_dr, extraps = Float64[], Float64[], Float64[], Float64[], Bool[] + + # Overlap is decided by the first crossing walking outward, so surfaces far beyond it cannot + # change the answer -- but each one costs a Riccati layer solve, and this scan runs whenever + # kinetic profiles are readable, flag or no flag. Bound the work rather than only the m index, + # and say so when the bound bites, since a "no overlap" verdict from a truncated list is not + # conclusive. + if length(found) > max_layer_solves + push!(notes, + "surface list truncated from $(length(found)) to the innermost $max_layer_solves for " * + "the layer solve (max_layer_solves); a \"no overlap\" verdict here is not conclusive") + found = found[1:max_layer_solves] + end + + for s in found + q_val = s.m / n_tor + # dq/dψ from the equilibrium spline in-grid, from the edge law outside it. + q1_val = if !s.extrap + Float64(equil.profiles.q_deriv(s.psi)) + elseif inv_q !== nothing + Float64(inv_q(s.psi; deriv=DerivOp(1))) # real profile, not a model + else + edge_q_law_dqdpsi(edge_fit, s.psi) + end + da_dpsi = _dr_dpsi(s.psi) + if !isfinite(da_dpsi) || da_dpsi == 0.0 + push!(notes, "m=$(s.m): da/dψ = $da_dpsi at ψ=$(@sprintf("%.6f", s.psi)); cannot convert width to flux units") + continue + end + + sing = SingType(; m=[s.m], n=[Int(n_tor)], psifac=s.psi, rho=sqrt(s.psi), q=q_val, q1=q1_val) + # Built one surface at a time so a single degenerate surface (e.g. ω_*e == ω_*i, + # which makes iota_e singular) is recorded and skipped rather than aborting the scan. + params = try + build_slayer_inputs(equil, [sing], profiles; rs_method=rs_method, + mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, + dr_val=0.0, dc_type=:none, theta=theta, + resistivity_model=resistivity_model, lnLambda_form=lnLambda_form) + catch err + err isa ArgumentError || rethrow() + push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): $(err.msg)") + continue + end + + lw = slayer_layer_thickness(params[1]) + if !isfinite(lw.delta_s_m) + push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): del_s Riccati did not converge") + continue + end + # Both widths must be finite. A NaN delta_visco would compare false against every + # neighbour in _first_overlap_limit, so the viscous criterion would silently report "no + # overlap" and `recommended` would quietly lose its more conservative half -- biasing the + # domain outward, which is the unsafe direction. + if !isfinite(lw.delta_dr) + push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): Eq. (100) width is not finite; " * + "surface excluded so the DR criterion cannot fail open") + continue + end + if !isfinite(lw.delta_visco) + push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): viscous width is not finite; " * + "surface excluded so the viscous criterion cannot fail open") + continue + end + + # Metres → normalized flux. This conversion is the whole reason the scan + # can compare a layer width against a surface separation. + push!(ms, s.m) + push!(ns, Int(n_tor)) + push!(psis, s.psi) + push!(rss, params[1].rs) + push!(dels_m, lw.delta_s_m) + push!(w_dels, lw.delta_s_m / abs(da_dpsi)) + push!(w_visc, lw.delta_visco / abs(da_dpsi)) + push!(w_dr, lw.delta_dr / abs(da_dpsi)) + push!(extraps, s.extrap) + end + + ph_dels, k_dels = _first_overlap_limit(psis, w_dels) + ph_visc, k_visc = _first_overlap_limit(psis, w_visc) + ph_dr, k_dr = _first_overlap_limit(psis, w_dr) + # The DR width is the criterion: near the separatrix the shear is strong enough that + # Fitzpatrick (2025) Sect. 5.6 puts every layer in that regime. The delta_s ODE and the + # viscous scale stay reported so the three can be compared, but they do not set the domain. + recommended = ph_dr + + return LayerOverlapScan(ms, ns, psis, rss, dels_m, w_dels, w_visc, w_dr, extraps, + ph_dels, ph_visc, ph_dr, recommended, k_dr, notes) +end + +# Walk outward and return (recommended psihigh, index of first overlapping surface). +# Overlap at k means surfaces k and k-1 are both contaminated, so the domain is cut +# at the inner edge of k-1. +function _first_overlap_limit(psi::Vector{Float64}, w::Vector{Float64}) + for k in 2:length(psi) + if psi[k] - w[k] / 2 < psi[k-1] + w[k-1] / 2 + return (psi[k-1] - w[k-1] / 2, k) + end + end + return (nothing, nothing) +end + +function Base.show(io::IO, ::MIME"text/plain", s::LayerOverlapScan) + println(io, "LayerOverlapScan: $(length(s.psi)) surface(s), n = ", + isempty(s.n) ? "-" : string(s.n[1])) + if !isempty(s.psi) + println(io, rpad("q", 8), rpad("psi", 13), rpad("r_s [m]", 11), + rpad("dpsi(del_s)", 14), rpad("dpsi(visco)", 14), rpad("dpsi(DR)", 14), "source") + for i in eachindex(s.psi) + println(io, + rpad("$(s.m[i])/$(s.n[i])", 8), + rpad(@sprintf("%.7f", s.psi[i]), 13), + rpad(@sprintf("%.5f", s.rs[i]), 11), + rpad(@sprintf("%.4e", s.width_delta_s[i]), 14), + rpad(@sprintf("%.4e", s.width_visco[i]), 14), + rpad(@sprintf("%.4e", s.width_dr[i]), 14), + s.extrapolated[i] ? "extrapolated" : "in-grid") + end + end + _fmt(x) = x === nothing ? "none (no overlap in range)" : @sprintf("%.7f", x) + println(io, " psihigh from |del_s| : ", _fmt(s.psihigh_delta_s)) + println(io, " psihigh from visco : ", _fmt(s.psihigh_visco)) + println(io, " psihigh from Eq.(100): ", _fmt(s.psihigh_dr)) + println(io, " recommended psihigh : ", _fmt(s.psihigh)) + for note in s.notes + println(io, " note: ", note) + end + return nothing +end diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 745a30857..80bce7b23 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -22,6 +22,7 @@ using ..Utilities import ..InnerLayer as InnerLayer include("LayerInputs.jl") +include("LayerOverlap.jl") include("Dispersion/Dispersion.jl") include("Runner/Runner.jl") @@ -30,5 +31,6 @@ import .Runner as Runner export InnerLayer, Dispersion, Runner export build_ggj_inputs +export resistive_layer_overlap, LayerOverlapScan end # module Tearing From d4cb367c38f26c6dd3fc078f502b08b3b3ee03cc Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 20 Aug 2026 20:39:17 -0400 Subject: [PATCH 4/6] ForceFreeStates - FEATURE - Cap the integration domain at the resistive-layer overlap `sing_lim!` gains an optional `psilim_cap`. Where the resistive layers of adjacent rational surfaces overlap, no surface retains a well-separated inner region and matched asymptotics is not defined, so that location is an upper bound on the useful domain. It is applied as a cap on `qlim` rather than by assigning `psilim` directly, which keeps it a bound rather than a setpoint: the existing `dmlim` / `qhigh` truncation still selects the final surface from inside it, and a cap lying beyond `psihigh` is inert by construction with no special case -- it simply loses the `min`. The domain can therefore only shrink, never grow, so no equilibrium ever has to be re-formed. Measured on the shipped DIII-D decks with the Eq. (100) width in the `:flux` label, the cap changes no retained surface on any of them: it is beyond `psihigh` for the ideal (0.995), riccati (0.995) and gal_resistive (0.993) decks, and on the SLAYER deck (0.9995) it binds by 1.1e-4 in psi without crossing a rational. The scan needs kinetic profiles, which two of those decks do not carry, so it is inert there regardless. --- src/ForceFreeStates/Sing.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl index 960530f6a..f1246c1c8 100644 --- a/src/ForceFreeStates/Sing.jl +++ b/src/ForceFreeStates/Sing.jl @@ -108,7 +108,8 @@ performed to find the corresponding `psilim` to integrate to. Note that the Newton iteration will be triggered if either `set_psilim_via_dmlim` is true or `ctrl.qhigh < equil.params.qmax`. Otherwise, the equilibrium edge values are used. """ -function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) +function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium; + psilim_cap::Union{Nothing,Real}=nothing) profiles = equil.profiles @@ -117,6 +118,20 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, intr.q1lim = profiles.q_deriv(profiles.xs[end]; hint=Ref(profiles.npts_minus_1)) intr.psilim = equil.params.psihigh_resolved + # Resistive-layer overlap imposes an UPPER BOUND on the domain: past the first pair of + # overlapping layers no surface retains a well-separated inner region, so matched asymptotics + # is not defined out there. Applied as a cap on qlim rather than as psilim directly, so the + # dmlim / qhigh truncation below still selects the final surface from inside the bound. + # Never widens the domain: a cap beyond psihigh is inert by construction. + if psilim_cap !== nothing && psilim_cap < intr.psilim + q_cap = profiles.q_spline(Float64(psilim_cap)) + if q_cap < intr.qlim + @info "Resistive-layer overlap caps the domain: qlim $(@sprintf("%.3f", intr.qlim)) -> " * + "$(@sprintf("%.3f", q_cap)) (psi $(@sprintf("%.6f", intr.psilim)) -> $(@sprintf("%.6f", Float64(psilim_cap))))" + intr.qlim = q_cap + end + end + # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent). The cutoff reads # the *resolved* toroidal range on `intr`, so callers must assign intr.nlow / intr.nhigh # before calling; an unresolved range is an error rather than a silent change of truncation From 5fe6b42753b72b54b71165228f61be49ceea1100 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 21 Aug 2026 11:39:06 -0400 Subject: [PATCH 5/6] InnerLayer.SLAYER - TEST - Check the layer width against the published Fitzpatrick JET case Rebuilds the model JET equilibrium of Fitzpatrick, Nucl. Fusion 2025 Sect. 5.8 from the paper's own definitions -- the Eqs. (32)-(35) safety factor with alpha-/alpha+ solved self-consistently against Eq. (36), and mtanh edge profiles anchored to its Fig. 8 -- and checks where adjacent resistive layers first overlap against the values the paper reports. The widths are driven through `slayer_parameters` and `slayer_layer_thickness` rather than recomputed here, so the test covers the shipped tau_R / tau_A / d_beta chain and the `delta_dr` formula itself. That distinction matters: the bug this test would have caught -- carrying Eq. (100)'s explicit (n|s|)^(-1/2) on top of a `lu` that already contains it -- lives in exactly that chain, and a test that reimplemented the width would have passed while the shipped one was wrong. Measured against the paper: Psi = 0.9985111 for n = 1 (paper 0.9985, deviation 1.1e-5) and 0.9949738 for n = 4 (paper 0.9952, deviation 2.3e-4). Bounds are set at 1e-4 and 1e-3, leaving roughly 9x and 4x margin; the model is quadrature and root-finding with no BLAS, so platform spread is orders of magnitude below either. The test logs both deviations so a future failure reports how far it moved rather than only that it moved. This is the suite's second check pinned to an external published reference rather than to GPEC's own history. --- test/runtests.jl | 1 + test/runtests_layer_overlap_fitzpatrick.jl | 129 +++++++++++++++++++++ 2 files changed, 130 insertions(+) create mode 100644 test/runtests_layer_overlap_fitzpatrick.jl diff --git a/test/runtests.jl b/test/runtests.jl index 7845eea60..540d99a28 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,7 @@ else include("./runtests_slayer_params.jl") include("./runtests_slayer_riccati.jl") include("./runtests_slayer_inputs.jl") + include("./runtests_layer_overlap_fitzpatrick.jl") include("./runtests_dispersion_residual.jl") include("./runtests_dispersion_coupled.jl") include("./runtests_dispersion_coupled_full.jl") diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl new file mode 100644 index 000000000..21b1efb50 --- /dev/null +++ b/test/runtests_layer_overlap_fitzpatrick.jl @@ -0,0 +1,129 @@ +# External-reference check of the diffusive-resistive layer width against a published +# result: Fitzpatrick, Nucl. Fusion 2025 (doi 10.1088/1741-4326/ae4fdd), Sect. 5.8-5.9. +# +# The paper reports that for its model JET equilibrium the resistive layers of adjacent +# rational surfaces first overlap at Psi = 0.9985 for n = 1 (its Fig. 9) and Psi = 0.9952 +# for n = 4 (its Fig. 10). This rebuilds that equilibrium from the paper's own definitions +# and drives GPEC's `delta_dr` through `slayer_parameters`, so the width formula and the +# tau_R / tau_A / d_beta chain feeding it are checked against a number in print rather than +# against GPEC's own history. +@testset "Layer overlap vs Fitzpatrick (2025) JET model" begin + using GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER: slayer_parameters, + slayer_layer_thickness, SpitzerHarmModel + using QuadGK: quadgk + using Roots: find_zero, Bisection + + # Sect. 5.8 parameters. + B0, R0, A_MIN = 3.45, 2.96, 1.25 + ZEFF, MNUM, CHI = 10.0, 2.0, 1.0 + Q0, Q95, Q105 = 1.01, 3.5, 4.0 + + # Model safety factor, Eqs. (32)-(35). alpha- depends on rhat95, which depends on the + # profile through Eq. (36), so the pair is solved self-consistently. + q_in(r, am) = Q0 - am * log(1 - r^2) + q_out(r, ap) = -ap * log(r^2 - 1) + q_any(r, am, ap) = r < 1 ? q_in(r, am) : q_out(r, ap) + psi_raw(r, am, ap) = quadgk(x -> x / q_any(x, am, ap), 0.0, r; rtol=1e-11)[1] + + am, ap = 1.0, 1.0 + for _ in 1:200 + psi_sep = psi_raw(1.0 - 1e-12, am, ap) + r95 = find_zero(r -> psi_raw(r, am, ap) / psi_sep - 0.95, (0.5, 1 - 1e-9), Bisection()) + f105(r) = (psi_sep + quadgk(x -> x / q_out(x, ap), 1 + 1e-12, r; rtol=1e-10)[1]) / psi_sep - 1.05 + hi = 1.0 + 1e-9 + while hi < 1.40 && f105(hi) < 0 # q_out turns negative past r = sqrt(2) + hi += 0.005 + end + r105 = hi >= 1.40 ? 1.20 : find_zero(f105, (1 + 1e-9, hi), Bisection()) + am_new = -(Q95 - Q0) / log(1 - r95^2) + ap_new = -Q105 / log(r105^2 - 1) + converged = abs(am_new - am) < 1e-12 && abs(ap_new - ap) < 1e-12 + am, ap = am_new, ap_new + converged && break + end + psi_sep = psi_raw(1.0 - 1e-12, am, ap) + q_of(r) = q_any(r, am, ap) + Psi_of(r) = r < 1 ? psi_raw(r, am, ap) / psi_sep : + (psi_sep + quadgk(x -> x / q_out(x, ap), 1 + 1e-12, r; rtol=1e-10)[1]) / psi_sep + dq_dr(r) = r < 1 ? am * 2r / (1 - r^2) : -ap * 2r / (r^2 - 1) + s_of(r) = r * dq_dr(r) / q_of(r) + + # mtanh edge profiles read off the paper's Fig. 8, anchored at Psi = 0.96 and at the + # separatrix (Sect. 5.9 states Te ~ 100 eV there). + function make_tanh(f096, f100, f_ped, f_sol) + a = (f_ped - f_sol) / 2 + x1 = atanh(clamp(1 - (f096 - f_sol) / a, -0.999, 0.999)) + x2 = atanh(clamp(1 - (f100 - f_sol) / a, -0.999, 0.999)) + d = 0.04 / (x2 - x1) + p0 = 0.96 - x1 * d + return P -> f_sol + a * (1 - tanh((P - p0) / d)) + end + te_of = make_tanh(400.0, 100.0, 600.0, 20.0) + ne_of = make_tanh(3.5e19, 1.2e19, 5.0e19, 3.0e18) + + # Width in units of rhat, through the shipped code path. + function width_hat(r, n) + P = Psi_of(r) + te = te_of(P) + ne = ne_of(P) + p = slayer_parameters(; n_e=ne, t_e=te, t_i=te, + omega=0.0, omega_e=1.0e4, omega_i=5.0e3, + qval=q_of(r), sval_r=s_of(r), bt=B0, + rs=r * A_MIN, R0=R0, mu_i=MNUM, zeff=ZEFF, + chi_perp=CHI, chi_tor=CHI, m=1, n=n, + resistivity_model=SpitzerHarmModel(), lnLambda_form=:nrl) + return slayer_layer_thickness(p).delta_dr / A_MIN + end + + function surfaces(n; rmin=0.90, rmax=1.35) + out = Tuple{Int,Float64}[] + for m in 1:400 + qt = m / n + qt < q_of(rmin) && continue + r = if qt < q_of(1 - 1e-13) + find_zero(x -> q_of(x) - qt, (rmin, 1 - 1e-13), Bisection()) + elseif qt > q_of(1 + 1e-13) && qt < q_of(rmax) + find_zero(x -> q_of(x) - qt, (1 + 1e-13, rmax), Bisection()) + else + nothing + end + r === nothing || push!(out, (m, r)) + end + sort!(out; by=t -> t[2]) + return out + end + + # Inner boundary of the overlap region: the first surface whose layer runs into its + # inner neighbour, reported at that neighbour's location. + function overlap_psi(n) + S = surfaces(n) + for k in 2:length(S) + gap = S[k][2] - S[k-1][2] + (width_hat(S[k-1][2], n) + width_hat(S[k][2], n)) / 2 >= gap && return Psi_of(S[k-1][2]) + end + return nothing + end + + # The self-consistent profile must reproduce the paper's own inputs before the widths + # mean anything. + @test isapprox(q_of(0.0), Q0; atol=1e-10) + @test am > 0 && ap > 0 + + psi1 = overlap_psi(1) + @test psi1 !== nothing + psi1 === nothing || @info "Fitzpatrick 2025 overlap, n = 1" psi=psi1 paper=0.9985 deviation=abs(psi1 - 0.9985) atol=1e-4 + # Paper: Psi = 0.9985 (n = 1); measured deviation 1.1e-5. The bound covers the pedestal read + # off Fig. 8 and the resistivity closure (Spitzer-Harm here against the paper's Eqs. 70-72), + # and still leaves ~9x margin. Pure quadrature and root-finding, no BLAS, so platform spread + # sits far below this. + @test isapprox(psi1, 0.9985; atol=1e-4) + + psi4 = overlap_psi(4) + @test psi4 !== nothing + psi4 === nothing || @info "Fitzpatrick 2025 overlap, n = 4" psi=psi4 paper=0.9952 deviation=abs(psi4 - 0.9952) atol=1e-3 + # Paper: Psi = 0.9952 (n = 4); measured deviation 2.3e-4. Looser than n = 1 because these + # surfaces sit further in, where the digitised pedestal shape matters more. + @test isapprox(psi4, 0.9952; atol=1e-3) + # The overlap region must move inward with n, which is the paper's reported scaling. + @test psi4 < psi1 +end From b3d65ecf43b8e87d21a7b7903b140b5ef2887fc4 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 21 Aug 2026 15:26:16 -0400 Subject: [PATCH 6/6] =?UTF-8?q?Regression=20-=20MINOR=20-=20Set=20the=20SL?= =?UTF-8?q?AYER=20=CE=B3=20threshold=20above=20the=20root=20search's=20own?= =?UTF-8?q?=20reproducibility?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The `diiid_slayer_n1` γ_Hz threshold of 1e-1 sits below what the threaded dispersion root search reproduces: re-running one commit against itself with fresh samples moved γ by 0.145 Hz at surface 3/1, 0.122 Hz at 2/1 and 0.076 Hz at 4/1, in discrete steps, with every layer input bit-identical. At 1e-1 the case therefore flags a null change as CHANGED roughly half the time. 2.5e-1 clears the largest measured step by 1.7x. This is a stopgap: the underlying fix is a reproducible root search, after which the threshold can come back down. --- regression-harness/cases/diiid_slayer_n1.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 540ca7e50..57ca1c688 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -119,7 +119,7 @@ h5path = "Tearing/Roots/gamma" type = "real_vector" extract = "first_3" label = "SLAYER γ_Hz [2/1,3/1,4/1]" -noise_threshold = 1e-1 +noise_threshold = 2.5e-1 # Same-source fresh re-runs step by 0.076 / 0.122 / 0.145 Hz at 4/1, 2/1, 3/1 (threaded root search); 1.7× the largest order = 33 # no_root flag (1 = extraction failed). Pinned for the inner three surfaces;