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; 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). 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 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 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 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