TEARING - BUG FIX - Convert Δ' to the r_s reference length before slab-layer matching

The SLAYER dispersion relation paired a ψ_N-referenced BVP Δ' with an
r_s-referenced layer side (Δ(Q), S=τ_R/τ_H on r_s, and the critical-Δ,
whose Ŵ_d is W/r_s). Apply the Frobenius reference-length transform
Δ̂_ij = K_i^(1/2+μ_i)·Δ'_ij·K_j^(μ_j−1/2), K = r_s·(dψ_N/dr)|_s,
μ = √(−D_I), at the matching point. Verified parameter-free against the
TJ circular benchmarks (median residual +1% over 23 points; absolute
2/1 agreement 17.5%→2.5% (β) and 17.6%→1.1% (ε)). BVP Δ' outputs are
unchanged; only the SLAYER matching (γ, Δ_eff) moves. GGJ untouched.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: d-burg

src/InnerLayer/SLAYER/LayerInputs.jl

@@ -289,6 +289,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 +318,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/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

src/Tearing/Runner/Result.jl

@@ -20,7 +20,9 @@ downstream inspection and HDF5 output.
     factor q of each analyzed surface, aligned with `params`. Empty when the
     analysis was built from bare parameters (`run_slayer_from_inputs` without
     the surface list), in which case the HDF5 writer skips them.
-  - `dp_matrix`           -- outer-region Δ' matrix used in the analysis
+  - `dp_matrix`           -- outer-region Δ' matrix used in the analysis.
+    SLAYER path: r_s-referenced (the ψ_N BVP matrix transformed by
+    `delta_prime_to_rs_reference`); GGJ path: the ψ_N matrix unchanged
   - `Q_root`              -- tearing eigenvalue(s) in normalized Q
     * length `nsurfaces` in `:uncoupled` mode
     * length `1` in `:coupled` mode (global eigenvalue normalized by

src/Tearing/Runner/run_slayer.jl

@@ -151,6 +151,39 @@ function _build_surface_coupling(model::GGJModel, params::GGJParameters,
     return surface_coupling(model, params, dp_diag)
 end
 
+# ---------------------------------------------------------------------
+# Reference-length conversion of the outer Δ' for the slab layer
+# ---------------------------------------------------------------------
+"""
+    delta_prime_to_rs_reference(dp_matrix, params) -> Matrix{ComplexF64}
+
+Convert the outer-region Δ' matrix from its ψ_N reference length (the
+STRIDE/BVP convention: Frobenius coefficients normalized per unit Δψ_N) to
+the r_s-based `x̂ = (r − r_s)/r_s` reference the slab layer works in
+(Fitzpatrick 2023 convention — the same reference used by `dc_tmp` and by
+the `S^(1/3)` Δ(Q) scale, both built on r_s).
+
+Near surface `k` the tearing-parity outer solution is
+`A_L·|x|^(1/2−μ) + A_S·|x|^(1/2+μ)` with `μ = √(−D_I)` (Glasser, Wang &
+Park 2016, Eq. 26). Rescaling the radial variable `x_ψ = K·x̂` with
+`K = r_s·(dψ_N/dr)|_s` maps the coefficients as `Â_L = A_L·K^(1/2−μ)` and
+`Â_S = A_S·K^(1/2+μ)`, so the response matrix (small coefficient at surface
+`i` per unit large coefficient at surface `j`) transforms as
+
+    Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2)
+
+whose diagonal is `K^(2μ)·Δ'_kk`; at D_I = −1/4 (μ = 1/2) this reduces to
+the textbook `Δ̂ = r_s·Δ'_phys`. `K` and `μ` are carried per surface in
+`SLAYERParameters.k_ref` / `.mu_mercier`; hand-built parameters default to
+`k_ref = 1`, making the conversion the identity.
+"""
+function delta_prime_to_rs_reference(dp_matrix::AbstractMatrix,
+    params::AbstractVector)
+    dl = [p.k_ref^(0.5 + p.mu_mercier) for p in params]
+    dr = [p.k_ref^(p.mu_mercier - 0.5) for p in params]
+    return Diagonal(dl) * Matrix{ComplexF64}(dp_matrix) * Diagonal(dr)
+end
+
 # ---------------------------------------------------------------------
 # Core analysis entry point that takes pre-built parameters.
 # ---------------------------------------------------------------------
@@ -194,6 +227,14 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters},
                 "(build_slayer_inputs for SLAYER, build_ggj_inputs for GGJ).")
         )
 
+    # Slab-layer path: convert Δ' from its ψ_N reference length to the
+    # r_s-based convention shared by the layer Δ(Q) and the critical-Δ (see
+    # `delta_prime_to_rs_reference`). GGJ is genuinely toroidal/ψ-based
+    # (its `rescale_delta` handles inner→outer units natively) — no conversion.
+    if !_is_ggj(model)
+        dp = delta_prime_to_rs_reference(dp, params)
+    end
+
     # The coupled determinant uses the reduced m×m (tearing-only) form, which
     # drops the interchange channel. For GGJ that channel carries the Glasser
     # interchange stabilization, so coupled-GGJ results omit real physics.

test/runtests_slayer_runner.jl

@@ -153,6 +153,41 @@
         @test_throws ArgumentError run_slayer_from_inputs(params, bad_dp, c)
     end
 
+    @testset "delta_prime_to_rs_reference: ψ_N → r_s conversion" begin
+        # Two surfaces with distinct K and μ; the transform is
+        # Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2).
+        K1, mu1 = 0.9, 0.54
+        K2, mu2 = 1.2, 0.60
+        p1 = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0,
+            P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3,
+            tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.4, R0=1.7, bt=2.0,
+            sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=2, n=1, ising=1,
+            k_ref=K1, mu_mercier=mu1)
+        p2 = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0,
+            P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3,
+            tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.5, R0=1.7, bt=2.0,
+            sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=3, n=1, ising=2,
+            k_ref=K2, mu_mercier=mu2)
+        dp = ComplexF64[10.0+1im 2.0-0.5im; 3.0+0im 1.5+2im]
+        out = Runner.delta_prime_to_rs_reference(dp, [p1, p2])
+        # Diagonal carries the scalar K^(2μ)
+        @test out[1, 1] ≈ K1^(2mu1) * dp[1, 1]
+        @test out[2, 2] ≈ K2^(2mu2) * dp[2, 2]
+        # Off-diagonals carry the split row/column factors
+        @test out[1, 2] ≈ K1^(0.5 + mu1) * K2^(mu2 - 0.5) * dp[1, 2]
+        @test out[2, 1] ≈ K2^(0.5 + mu2) * K1^(mu1 - 0.5) * dp[2, 1]
+        # At the slab point μ = 1/2 the diagonal factor is exactly K
+        pslab = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0,
+            P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3,
+            tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.4, R0=1.7, bt=2.0,
+            sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=2, n=1, ising=1,
+            k_ref=0.8, mu_mercier=0.5)
+        dp1 = ComplexF64[5.0+0im;;]
+        @test Runner.delta_prime_to_rs_reference(dp1, [pslab])[1, 1] ≈ 0.8 * 5.0
+        # Default parameters (k_ref = 1) give the identity regardless of μ
+        @test Runner.delta_prime_to_rs_reference(dp, [_mk_params(), _mk_params()]) ≈ dp
+    end
+
     @testset "run_slayer_from_inputs: coupled mode finds known root" begin
         # Build a 2-surface problem with a known coupled root by construction.
         p1 = _mk_params(; rs=0.5, lu=1.0e7, tauk=1.0e-4, Q_e=-1.0, Q_i=0.5,