Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching

src/InnerLayer/InnerLayer.jl

@@ -25,7 +25,7 @@ import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface
 
 import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear
 import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths
-import .SLAYER: surface_minor_radius, surface_da_dpsi, build_slayer_inputs
+import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs
 
 export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile
 export GGJ, GGJModel, GGJParameters
@@ -37,6 +37,6 @@ export delta_convergence, solution_profile, asymptotic_profile, q4_surface_bench
 
 export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear
 export riccati_del_s, slayer_layer_thickness, LayerWidths
-export surface_minor_radius, surface_da_dpsi, build_slayer_inputs
+export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs
 
 end # module InnerLayer

src/InnerLayer/SLAYER/LayerInputs.jl

@@ -9,14 +9,14 @@
 # Geometry extraction:
 #   - Minor radius at the outboard midplane (θ = 0) via
 #     `equil.rzphi_rsquared((ψ, 0.0))`.
-#   - `da/dψ` via central finite difference on the same bicubic.
+#   - `da/dψ` from the interpolant's own analytic ψ-derivative.
 #   - r-based magnetic shear via `r_based_shear(rs, q, q1, da/dψ)` (defined
 #     in LayerParameters.jl).
 
 using ..Utilities: KineticProfiles
 using ...Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel,
     coulomb_log_e, nu_star_e
-using FastInterpolations: DerivOp
+using FastInterpolations: DerivOp, integrate, cubic_interp, cumulative_integrate, ExtendExtrap
 
 """
     surface_minor_radius(equil, psi; theta=0.0) -> Float64
@@ -32,33 +32,123 @@ function surface_minor_radius(equil, psi::Real; theta::Real=0.0)
 end
 
 """
-    surface_da_dpsi(equil, psi; theta=0.0, h=1e-5) -> Float64
+    surface_da_dpsi(equil, psi; theta=0.0) -> Float64
 
-Central finite-difference approximation of `d(minor radius)/dψ` at `psi`.
-Falls back to one-sided differences near the flux-coordinate boundaries
-(0 or 1).
+Analytic ψ-derivative of the minor radius at `psi` and poloidal angle
+`theta`, taken from the `rzphi_rsquared` interpolant's own ψ-derivative as
+`da/dψ = (∂r²/∂ψ) / (2a)`. Valid wherever the interpolant is, including
+under extrapolation past the ψ grid. Diverges at the magnetic axis, where
+`a ~ √ψ`; callers evaluating near `ψ = 0` must check `isfinite`.
 """
-function surface_da_dpsi(equil, psi::Real; theta::Real=0.0, h::Real=1e-5)
-    psi_f = Float64(psi)
-    # Clamp to safe sampling range within (0, 1)
-    eps_edge = 10 * h
-    lo = psi_f - h
-    hi = psi_f + h
-    if lo < eps_edge
-        # one-sided forward
-        a0 = surface_minor_radius(equil, max(psi_f, eps_edge); theta=theta)
-        a1 = surface_minor_radius(equil, max(psi_f, eps_edge) + h; theta=theta)
-        return (a1 - a0) / h
-    elseif hi > 1.0 - eps_edge
-        # one-sided backward
-        a0 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge) - h; theta=theta)
-        a1 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge); theta=theta)
-        return (a1 - a0) / h
+function surface_da_dpsi(equil, psi::Real; theta::Real=0.0)
+    return _da_dpsi_at_theta(equil, Float64(psi), Float64(theta))
+end
+
+# d(√r²)/dψ at one (ψ, θ) from the interpolant's own ψ-derivative. Shared by all radial-label
+# conventions so none carries its own stencil.
+@inline function _da_dpsi_at_theta(equil, psi::Float64, theta::Float64)
+    r_sq = equil.rzphi_rsquared((psi, theta))
+    a = sqrt(max(r_sq, 0.0))
+    a > 0 || return Inf   # magnetic axis: a ~ √ψ, so da/dψ genuinely diverges
+    return equil.rzphi_rsquared((psi, theta); deriv=DerivOp(1, 0)) / (2a)
+end
+
+"""
+    radial_label(equil; rs_method=:midplane, theta=0.0) -> (r_at, dr_dpsi_at)
+
+Build the pair of closures `r_at(ψ)` and `dr_dpsi_at(ψ)` defining one radial
+label for the layer stack. Both closures must come from the same label,
+because the r-based shear `(r/q)dq/dr`, `τ_R = μ₀r²/η`, `τ_E = r²/χ`, the
+`d_β/r` normalization, the Δ' reference-length factor `k_ref = r_s/(da/dψ)`,
+and any metre-to-ψ width conversion all have to live in one coordinate; a
+mismatched `r` and `dr/dψ` silently corrupts every one of them. All
+derivatives are analytic — no label carries a finite-difference stencil.
+
+# Labels
+
+  - `:midplane` -- outboard-midplane chord from the magnetic axis at `theta`
+    (historical default), natural for comparison with midplane diagnostics.
+  - `:halfwidth` -- midplane half-chord, the mean of the outboard and inboard
+    chords at `θ = 0` and `θ = 0.5`; shift-free. Coincides with the flux label
+    on circular equilibria but is its own convention on shaped ones.
+  - `:fsa` -- θ-mean surface radius, a 128-point midpoint mean of the local
+    minor radius; the closest geometric approximation to `:flux` at interior
+    surfaces of shaped equilibria.
+  - `:volume` -- cylinder-equivalent label `√(V(ψ)/(2π²R₀))`, the
+    Rutherford-literature convention.
+  - `:flux` -- toroidal-flux label. Fitzpatrick, Nucl. Fusion (2025),
+    Eq. 30: `dψ_p/dr = B₀ r g/q` integrates to `ψ_t = B₀r²/2`, so
+    `r = √(2ψ_t/B₀)` with `ψ_t = psio·∫₀^ψ (q/g) dψ′` and `g = F/(B₀R₀)`.
+    Defined from flux alone, it carries no circular-cross-section assumption,
+    and its derivative `dr/dψ ∝ q` grows toward a separatrix where the
+    geometric labels' `da/dψ` collapses.
+
+On shaped equilibria the labels agree at low-q surfaces and diverge strongly
+near the edge, where the slab-layer matching is label-ambiguous regardless of
+choice. The label is selected programmatically; it is not exposed via TOML.
+"""
+function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0)
+    theta_f = Float64(theta)
+
+    _flux_r, _flux_dr = if rs_method === :flux
+        b0f = Float64(equil.params.b0)
+        R0f = Float64(equil.ro)
+        psiof = Float64(equil.psio)
+        xs_f = collect(Float64, equil.profiles.xs)
+        # g = F/(B0 R0) departs from 1 by ~3% on a DIII-D-like deck and ~0.8% on a circular
+        # one, so carry it rather than assuming g = 1: r² = 2∫(q/g)dψ_p/B0.
+        _g_at(x) = Float64(equil.profiles.F_spline(x)) / (2π * b0f * R0f)
+        qg = [Float64(equil.profiles.q_spline(x)) / _g_at(x) for x in xs_f]
+        Phi = collect(Float64, cumulative_integrate(cubic_interp(xs_f, qg)))
+        r_knots = sqrt.(max.(2 .* psiof .* Phi ./ b0f, 0.0))
+        rspl = cubic_interp(xs_f, r_knots; extrap=ExtendExtrap())
+        (ψ -> Float64(rspl(Float64(ψ))),
+            ψ -> psiof * Float64(equil.profiles.q_spline(Float64(ψ))) /
+                 (b0f * _g_at(Float64(ψ)) * max(Float64(rspl(Float64(ψ))), eps())))
     else
-        a_plus = surface_minor_radius(equil, psi_f + h; theta=theta)
-        a_minus = surface_minor_radius(equil, psi_f - h; theta=theta)
-        return (a_plus - a_minus) / (2h)
+        (nothing, nothing)
     end
+
+    _a_at(ψ, θ) = sqrt(max(equil.rzphi_rsquared((Float64(ψ), Float64(θ))), 0.0))
+
+    _rs_at(ψ) =
+        if rs_method === :fsa
+            N = 128
+            s = 0.0
+            @inbounds for k in 1:N
+                s += _a_at(ψ, (k - 0.5) / N)
+            end
+            s / N
+        elseif rs_method === :halfwidth
+            0.5 * (_a_at(ψ, 0.0) + _a_at(ψ, 0.5))
+        elseif rs_method === :volume
+            V = integrate(equil.profiles.dVdpsi_spline, 1e-4, Float64(ψ))
+            sqrt(max(V, 0.0) / (2π^2 * equil.ro))
+        elseif rs_method === :flux
+            _flux_r(ψ)
+        else
+            surface_minor_radius(equil, ψ; theta=theta_f)
+        end
+
+    _da_dpsi_at(ψ) =
+        if rs_method === :fsa
+            N = 128
+            s = 0.0
+            @inbounds for k in 1:N
+                s += _da_dpsi_at_theta(equil, Float64(ψ), (k - 0.5) / N)
+            end
+            s / N
+        elseif rs_method === :halfwidth
+            0.5 * (_da_dpsi_at_theta(equil, Float64(ψ), 0.0) + _da_dpsi_at_theta(equil, Float64(ψ), 0.5))
+        elseif rs_method === :volume
+            Float64(equil.profiles.dVdpsi_spline(ψ)) / (4π^2 * equil.ro * max(_rs_at(ψ), eps()))
+        elseif rs_method === :flux
+            _flux_dr(ψ)
+        else
+            _da_dpsi_at_theta(equil, Float64(ψ), theta_f)
+        end
+
+    return (_rs_at, _da_dpsi_at)
 end
 
 """
@@ -114,6 +204,12 @@ profiles, without an intermediate file round-trip.
     a prescribed value. (For `dc_type=:rfitzp` and `:lar`, dgeo_val is
     not consulted.)
   - `dc_type`   -- `:none` (default), `:lar`, `:rfitzp`, or `:toroidal`.
+  - `rs_method` -- radial label defining `r_s` for the whole layer stack:
+    `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See
+    [`radial_label`](@ref) for the definitions. S, the r-based shear, W_d,
+    and the Δ' reference-length factor `k_ref` all follow the choice
+    together, so every option is self-consistent. Not exposed via TOML —
+    programmatic use only.
   - `theta`     -- poloidal angle at which to measure minor radius (default
     `0.0`, outboard midplane).
   - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`,
@@ -156,38 +252,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             Float64(bt(ψ))
         end
 
-    # Minor-radius extractor: `:midplane` = outboard-midplane chord
-    # (original behavior); `:fsa` = θ-mean of √rzphi_rsquared, the
-    # flux-surface-averaged minor radius.
-    _rs_at(ψ) =
-        if rs_method === :fsa
-            integrand(θ) = sqrt(equil.rzphi_rsquared((Float64(ψ), Float64(θ))))
-            N = 128
-            s = 0.0
-            @inbounds for k in 1:N
-                s += integrand((k - 0.5) / N)
-            end
-            s / N
-        else
-            surface_minor_radius(equil, ψ; theta=theta)
-        end
-    _da_dpsi_at(ψ) =
-        if rs_method === :fsa
-            # central finite difference on _rs_at
-            h = 1e-5
-            lo = ψ - h
-            hi = ψ + h
-            eps_edge = 10h
-            if lo < eps_edge
-                (_rs_at(max(ψ, eps_edge) + h) - _rs_at(max(ψ, eps_edge))) / h
-            elseif hi > 1.0 - eps_edge
-                (_rs_at(min(ψ, 1.0 - eps_edge)) - _rs_at(min(ψ, 1.0 - eps_edge) - h)) / h
-            else
-                (_rs_at(ψ + h) - _rs_at(ψ - h)) / (2h)
-            end
-        else
-            surface_da_dpsi(equil, ψ; theta=theta)
-        end
+    _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta)
 
     # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline
     # derivatives. When `compute_omega_star=true` we override any ω_*e/ω_*i
@@ -289,6 +354,20 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             _eval(dgeo_val, psi)
         end
 
+        # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s
+        # and μ = √(−D_I) with D_I = E + F + H − 1/4 (Glasser-Greene-Johnson 1975).
+        # A Mercier-unstable surface (D_I ≥ 0) has no real exponent; μ → 0 there,
+        # which makes the conversion factor K^(2μ) → 1 continuously (Δ' left raw).
+        k_ref_k = da_dpsi == 0.0 ? 1.0 : rs / da_dpsi
+        mu_k = if rg === nothing
+            @warn("build_slayer_inputs: sing.restype not populated; using the " *
+                  "slab Mercier exponent μ = 1/2 for the Δ' reference-length " *
+                  "conversion at all such surfaces.", maxlog=1)
+            0.5
+        else
+            sqrt(max(-(rg.E + rg.F + rg.H - 0.25), 0.0))
+        end
+
         out[k] = slayer_parameters(;
             n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i,
             omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use,
@@ -304,7 +383,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             f_trap=f_trap_kw,
             nu_e_star=nu_e_star_kw,
             R_major_eff=R_major_eff,
-            lnLambda_form=lnLambda_form
+            lnLambda_form=lnLambda_form,
+            k_ref=k_ref_k,
+            mu_mercier=mu_k
         )
     end
     return out

src/InnerLayer/SLAYER/LayerParameters.jl

@@ -46,6 +46,15 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 | `d_beta`   | Beta-weighted ion length scale c_β · d_i [m]                      |
 | `dc_tmp`   | Critical-Δ offset from chi_parallel matching                      |
 | `dc_type`  | Selector for `dc_tmp` formula                                     |
+| `k_ref`    | Reference-length ratio K = r_s · (dψ_N/dr) at this surface (1 = no Δ' conversion) |
+| `mu_mercier` | Mercier exponent μ = √(−D_I) governing the Δ' reference-length conversion (1/2 = slab/cylindrical value) |
+
+`k_ref` and `mu_mercier` feed the ψ_N → r_s reference-length conversion of
+the outer Δ' matrix (see `delta_prime_to_rs_reference` in the Tearing
+runner): the slab layer, its `dc_tmp` critical-Δ, and the `S^(1/3)` Δ(Q)
+scale are all referenced to unit `x̂ = (r−r_s)/r_s`, while the outer BVP Δ'
+is referenced to unit Δψ_N. Hand-built parameters default to `k_ref = 1`,
+which makes the conversion the identity.
 
 The complex normalized growth rate `Q = ω + iγ` is **not** stored here;
 it is passed as a separate argument to `solve_inner`.
@@ -85,6 +94,10 @@ Base.@kwdef struct SLAYERParameters <: InnerLayerParameters
     # Critical-Δ offset
     dc_tmp::Float64 = 0.0
     dc_type::Symbol = :none
+
+    # Reference-length conversion inputs for the outer Δ' (ψ_N → r_s-based x̂)
+    k_ref::Float64 = 1.0
+    mu_mercier::Float64 = 0.5
 end
 
 # Allowed dc_type values for the critical-Δ offset. `:none` is the default
@@ -202,6 +215,12 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl
   - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula
   - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal`
   - `ising`   -- singular-surface index for traceability
+  - `k_ref`   -- reference-length ratio K = r_s·(dψ_N/dr) at the surface,
+    used by the Tearing runner to convert the ψ_N-referenced outer Δ' to
+    the r_s-referenced convention this layer works in (default `1.0`,
+    i.e. no conversion; `build_slayer_inputs` fills it from the equilibrium)
+  - `mu_mercier` -- Mercier exponent μ = √(−D_I) for the same conversion
+    (default `0.5`, the slab/cylindrical value at D_I = −1/4)
 
 # Resistivity kwargs
 
@@ -252,7 +271,9 @@ function slayer_parameters(;
     f_trap::Union{Real,Nothing}=nothing,
     nu_e_star::Union{Real,Nothing}=nothing,
     R_major_eff::Union{Real,Nothing}=nothing,
-    lnLambda_form::Symbol=:nrl)
+    lnLambda_form::Symbol=:nrl,
+    k_ref::Real=1.0,
+    mu_mercier::Real=0.5)
 
     # Coulomb logarithm shared by the resistivity closure and τ_ee.
     lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form)
@@ -358,6 +379,7 @@ function slayer_parameters(;
         rs=rs, R0=R0, bt=bt, sval_r=sval_r,
         dr_val=dr_val, dgeo_val=dgeo_val,
         eta=eta, d_beta=d_beta,
-        dc_tmp=dc_tmp, dc_type=dc_type
+        dc_tmp=dc_tmp, dc_type=dc_type,
+        k_ref=k_ref, mu_mercier=mu_mercier
     )
 end

src/InnerLayer/SLAYER/SLAYER.jl

@@ -56,7 +56,7 @@ include("LayerInputs.jl")
 export SLAYERModel, SLAYERParameters, slayer_parameters
 export r_based_shear
 export riccati_del_s, slayer_layer_thickness, LayerWidths
-export surface_minor_radius, surface_da_dpsi, build_slayer_inputs
+export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs
 export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel
 
 end # module SLAYER

src/Tearing/Dispersion/SurfaceCoupling.jl

@@ -53,10 +53,14 @@ end
     surface_coupling(model::SLAYERModel, params::SLAYERParameters,
                      dp_diag::Number; dc::Real=0.0) -> SurfaceCoupling
 
-SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)` so that
-the dimensionless Δ from `riccati_f` is mapped to outer ψ-units before
-subtraction from the Δ' diagonal. `tauk` is taken from `params.tauk` for use
-by `MultiSurfaceCoupling` Q rescaling.
+SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)`, which
+maps the dimensionless inner-layer Δ from `riccati_f` to the r_s-referenced
+outer convention (S^(1/3) with S = τ_R/τ_H built on r_s; the Lundquist
+number, `dc`, and the layer all share the `x̂ = (r−r_s)/r_s` reference
+length). `dp_diag` must already be in that same r_s reference — the Tearing
+runner converts the ψ_N-referenced BVP Δ' via `delta_prime_to_rs_reference`
+before building couplings. `tauk` is taken from `params.tauk` for use by
+`MultiSurfaceCoupling` Q rescaling.
 """
 function surface_coupling(model::SLAYERModel, params::SLAYERParameters,
     dp_diag::Number; dc::Real=0.0)

src/Tearing/Runner/HDF5Output.jl

@@ -95,6 +95,9 @@ const TEARING_H5_ANNOTATIONS = [
     "PerSurface/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)),
     "PerSurface/D_c_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor-Hastie-Helander 2015 Eq. 59)", dims=("surface",)),
     "PerSurface/D_c_type" => (; long_name="per-surface D_c prescription label", dims=("surface",)),
+    "PerSurface/k_ref" => (; long_name="reference-length ratio K = r_s·(dψ_N/dr) at each surface", dims=("surface",)),
+    "PerSurface/mu_mercier" => (; long_name="Mercier exponent μ = √(−D_I) at each surface", dims=("surface",)),
+    "PerSurface/delta_prime_conversion" => (; long_name="ψ_N → r_s reference-length factor K^(2μ) applied to each Δ' diagonal", dims=("surface",)),
     "PerSurface/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",)),
     "PerSurface/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)),
     "PerSurface/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)),
@@ -103,7 +106,7 @@ const TEARING_H5_ANNOTATIONS = [
     "PerSurface/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)),
     "PerSurface/tau_A" => (; long_name="Alfvén time τ_A per surface (GGJ layer parameters)", units="s", dims=("surface",)),
     "PerSurface/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)),
-    "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces", dims=("surface_row", "surface_col")),
+    "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces, as used in the matching (SLAYER path: r_s-referenced via K^(2μ); the ψ_N-referenced BVP matrix is SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")),
     "Roots/Q_root" => (; long_name="complex dispersion-root normalized frequency Q (NaN = no root)", dims=("surface",)),
     "Roots/omega" =>
         (; long_name="mode rotation angular frequency ω = Re(Q)/τ_k of each root", units="rad/s", dims=("surface",)),
@@ -184,6 +187,11 @@ function _write_per_surface!(g, params::AbstractVector{SLAYERParameters},
     ps["D_geo"] = Float64[p.dgeo_val for p in params]
     # Store dc_type per-surface as string array
     ps["D_c_type"] = String[String(p.dc_type) for p in params]
+    # Reference-length conversion applied to Δ' (ψ_N → r_s): K, μ, and the
+    # diagonal factor K^(2μ) actually multiplying each Δ'_kk.
+    ps["k_ref"] = Float64[p.k_ref for p in params]
+    ps["mu_mercier"] = Float64[p.mu_mercier for p in params]
+    ps["delta_prime_conversion"] = Float64[p.k_ref^(2 * p.mu_mercier) for p in params]
 
     ps["Delta_prime_matrix"] = dp_matrix
     return nothing