InnerLayer.SLAYER - FEATURE - Add the Fitzpatrick flux radial label and unify the label plumbing

Factor the rs_method radial-label options out of build_slayer_inputs into one
exported radial_label(equil; rs_method, theta) returning (r, dr/dpsi) closures,
and add a fifth label :flux — Fitzpatrick's toroidal-flux surface label
(Nucl. Fusion 2025, Eq. 30), r = sqrt(2 psi_t / B0) with the g = F/(B0 R0)
correction carried. All five labels now use analytic psi-derivatives from the
interpolants themselves; the finite-difference stencils (which silently clamped
within 1e-4 of the flux boundaries) are removed. The :fsa branch gains the
max(r^2, 0) guard needed on extrapolated surfaces.

k_ref and the K^(2mu) Delta-prime conversion consume the same closures, so the
label choice drives S, the r-based shear, W_d and the Delta-prime reference
length together for every option. Default :midplane behavior is unchanged up to
the analytic-derivative refinement. Shared infrastructure for the resistive
layer-overlap psihigh cap branch, which evaluates its criterion in :flux.

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

Author: d-burg

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, integrate
+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,19 +204,11 @@ 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
-    (S, r-based shear, W_d, and the Δ' reference-length factor `k_ref` all
-    follow it together, so every choice is self-consistent). `:midplane`
-    (default): outboard-midplane chord from the magnetic axis — natural for
-    comparison with midplane diagnostics. `:halfwidth`: midplane half-chord
-    `(R_out − R_in)/2` — shift-free, the closest stand-in for the circular
-    flux label the Fitzpatrick layer formulas are derived in (reproduces the
-    cylindrical-theory label to ~1% on circular benchmark equilibria).
-    `:fsa`: θ-mean surface radius. `:volume`: cylinder-equivalent label
-    `√(V(ψ)/(2π²R₀))`, the Rutherford-literature convention. On shaped
-    equilibria the labels agree at low-q surfaces (~±3% in growth rate at
-    q=2) and diverge strongly near the edge, where the slab-layer matching
-    is label-ambiguous regardless of choice. Not exposed via TOML —
+  - `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).
@@ -170,47 +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; `:halfwidth` = midplane half-chord
-    # (R_out − R_in)/2, the shift-free circular flux label; `:volume` = the
-    # cylinder-equivalent label r_V = √(V(ψ)/(2π²R₀)).
-    _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
-        elseif rs_method === :halfwidth
-            0.5 * (surface_minor_radius(equil, ψ; theta=0.0) +
-                   surface_minor_radius(equil, ψ; theta=0.5))
-        elseif rs_method === :volume
-            lo = 1e-4
-            V = integrate(equil.profiles.dVdpsi_spline, lo, Float64(ψ))
-            sqrt(max(V, 0.0) / (2π^2 * equil.ro))
-        else
-            surface_minor_radius(equil, ψ; theta=theta)
-        end
-    _da_dpsi_at(ψ) =
-        if rs_method !== :midplane
-            # 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

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

test/runtests_slayer_inputs.jl

@@ -42,8 +42,8 @@
         @test r1 > 0
     end
 
-    @testset "surface_da_dpsi: FD agrees with numerical derivative" begin
-        # Reference via a tighter FD
+    @testset "surface_da_dpsi: agrees with FD reference" begin
+        # The analytic ψ-derivative must reproduce a tight FD of the minor radius.
         for psi in (0.1, 0.4, 0.7)
             h_ref = 1e-4
             r_p = surface_minor_radius(equil, psi + h_ref)
@@ -53,16 +53,24 @@
         end
     end
 
-    @testset "surface_da_dpsi: one-sided near boundaries" begin
-        # Near ψ=0 and ψ=1, the function falls back to one-sided FD and
-        # should still produce a finite positive number (minor radius is
-        # still increasing).
+    @testset "surface_da_dpsi: finite near the boundaries" begin
+        # The analytic form needs no clamping: near ψ=0 and ψ=1 it still returns a
+        # finite positive number (large near the axis, where a ~ √ψ).
         d_near_axis = surface_da_dpsi(equil, 1e-6)
         d_near_edge = surface_da_dpsi(equil, 1.0 - 1e-6)
         @test isfinite(d_near_axis) && d_near_axis > 0
         @test isfinite(d_near_edge) && d_near_edge > 0
     end
 
+    @testset "radial_label: analytic derivatives match FD" begin
+        for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux)
+            r_at, dr_at = radial_label(equil; rs_method=rsm)
+            h = 1e-5
+            ref = (r_at(0.5 + h) - r_at(0.5 - h)) / (2h)
+            @test dr_at(0.5) ≈ ref rtol = 1e-3
+        end
+    end
+
     @testset "build_slayer_inputs: returns correct per-surface data" begin
         sings = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1),
             _mk_sing(psi=0.6, q=3.0, q1=2.5, m=3, n=1)]
@@ -161,7 +169,7 @@
     @testset "build_slayer_inputs: rs_method radial labels are self-consistent" begin
         sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)]
         got = Dict{Symbol,Any}()
-        for rsm in (:midplane, :halfwidth, :fsa, :volume)
+        for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux)
             sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, rs_method=rsm)
             got[rsm] = sl[1]
             @test isfinite(sl[1].rs) && sl[1].rs > 0
@@ -174,6 +182,7 @@
         # Labels genuinely differ (each self-consistent set has its own rs, S, k_ref)
         @test got[:fsa].rs != got[:midplane].rs
         @test got[:volume].lu != got[:midplane].lu
+        @test got[:flux].rs != got[:midplane].rs
     end
 
     @testset "build_slayer_inputs: dc_type propagates and dr_val activates offset" begin