ForceFreeStates - FEATURE! - Cap the integration domain where resistive layers overlap

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;

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,18 +532,28 @@ 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,
     min_spacing::Float64=MIN_KNOT_SPACING,
     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).

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

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