logan-nc · Oct 30, 2026 · Oct 27, 2026 · Oct 27, 2026 · Oct 29, 2026 · Oct 29, 2026
diff --git a/docs/src/stability.md b/docs/src/stability.md
@@ -379,7 +379,7 @@ do not resolve the H-mode pedestal transport barrier, where power-balance χ fal
 
 | Dataset | Contents |
 |---|---|
-| `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius |
+| `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius √(2ψ_t/B₀) |
 | `delta_s_abs`, `width_delta_s` | Riccati layer thickness, in metres and as a Δψ width |
 | `width_visco`, `width_dr` | visco-resistive and diffusive-resistive (Eq. 100) Δψ widths; `width_dr` sets the cut |
 | `extrapolated` | 1 where the surface lies beyond the equilibrium grid |

diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl
@@ -7,10 +7,9 @@
 # round-trip.
 #
 # Geometry extraction:
-#   - Minor radius at the outboard midplane (θ = 0) via
-#     `equil.rzphi_rsquared((ψ, 0.0))`.
-#   - `da/dψ` from the interpolant's own analytic ψ-derivative.
-#   - r-based magnetic shear via `r_based_shear(rs, q, q1, da/dψ)` (defined
+#   - r_s and dr/dψ from one radial label (`radial_label`; by default the
+#     toroidal-flux label r = √(2ψ_t/B₀)), both analytic.
+#   - r-based magnetic shear via `r_based_shear(rs, q, q1, dr/dψ)` (defined
 #     in LayerParameters.jl).
 
 using ..Utilities: KineticProfiles
@@ -54,7 +53,7 @@ end
 end
 
 """
-    radial_label(equil; rs_method=:midplane, theta=0.0) -> (r_at, dr_dpsi_at)
+    radial_label(equil; rs_method=:flux, 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,
@@ -91,7 +90,7 @@ 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, q_beyond=nothing)
+function radial_label(equil; rs_method::Symbol=:flux, theta::Real=0.0, q_beyond=nothing)
     theta_f = Float64(theta)
 
     _flux_r, _flux_dr = if rs_method === :flux
@@ -294,14 +293,14 @@ profiles, without an intermediate file round-trip.
     without one it falls back to the cylindrical `n·|s|/R₀` with a warning.
   - `dc_type`   -- `:none` (default), `:lar`, `:fitzpatrick`, or `:toroidal`.
   - `rs_method` -- radial label defining `r_s` for the whole layer stack:
-    `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See
+    `:flux` (default), `:midplane`, `:halfwidth`, `:fsa`, or `:volume`. 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).
+  - `theta`     -- poloidal angle of the minor-radius chord, used only by
+    `rs_method=:midplane` (default `0.0`, outboard midplane).
 
   - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`,
     `SpitzerModel()`, or `SpitzerHarmModel()` (legacy Fitzpatrick σ_∥).
@@ -317,7 +316,7 @@ profiles, without an intermediate file round-trip.
 function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     bt=nothing,
     R0=nothing,
-    rs_method::Symbol=:midplane,
+    rs_method::Symbol=:flux,
     q_beyond=nothing,
     mu_i::Real=2.0,
     zeff=1.0,

diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl
@@ -151,9 +151,9 @@ Riccati solve; `delta_s` is solved per surface as a reported comparison only.
     ψ = 1 on the real q. Geometry needs no clamp of its own -- `radial_label` derivatives are
     analytic and valid under extrapolation
   - `extrapolate` -- search past `psihigh_safe` on the separatrix q-law (default `true`)
-  - `theta` -- poloidal angle for the minor-radius chord (default `0.0`)
+  - `theta` -- poloidal angle for the minor-radius chord, used only by `rs_method=:midplane` (default `0.0`)
   - `rs_method` -- radial label the widths and their ψ-conversion Jacobian are taken in
-    (default `:midplane`; see `radial_label`). The caller wiring this scan into the domain
+    (default `:flux`; see `radial_label`). The caller wiring this scan into the domain
     truncation passes `:flux`, the label Fitzpatrick Eq. (100) is anchored to
   - `bt`, `mu_i`, `zeff`, `chi_perp_e`, `chi_perp_i`, `chi_tor`, `tau_E`, `P_perp_model`, `resistivity_model`, `lnLambda_form`
     -- passed through to `build_slayer_inputs`. The defaults are placeholders: a caller
@@ -172,7 +172,7 @@ function resistive_layer_overlap(equil, profiles::KineticProfiles;
     psi_gap::Real=1e-9,
     extrapolate::Bool=true,
     theta::Real=0.0,
-    rs_method::Symbol=:midplane,
+    rs_method::Symbol=:flux,
     bt=nothing,
     mu_i::Real=2.0,
     zeff=1.0,

diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl
@@ -104,8 +104,10 @@ constructor.
     offset (Δ_crit ≡ 0) for every `D_c_type`; `D_geo = 0.0` does so only for
     `:toroidal`. `:toroidal` also takes the χ∥ closure's parallel wavenumber
     from the equilibrium (`toroidal_kpar`), which has no override
-  - `theta_sample` -- poloidal angle at which to sample minor radius
-    (default 0.0, outboard midplane)
+  - `theta_sample` -- poloidal angle of the minor-radius chord of the
+    `:midplane` radial label, used only by API callers that pass
+    `rs_method = :midplane`; the default toroidal-flux label does not depend
+    on angle, so a nonzero value warns (default 0.0)
   - `resistivity_model` -- η closure setting τ_R = μ₀r_s²/η: `:sauter`
     (default, neoclassical F_33), `:redl`, `:spitzer` (Sauter 18a, no
     trapped correction), or `:spitzer_harm` (legacy Fitzpatrick σ_∥)
@@ -242,13 +244,19 @@ const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson)
 const _VALID_ROOT_CLASSIFICATIONS = (:threshold, :winding)
 
 function validate(ctrl::SLAYERControl)
+    iszero(ctrl.theta_sample) || @warn "SLAYERControl: theta_sample = $(ctrl.theta_sample) has no effect: the layer quantities " *
+          "use the toroidal-flux radial label, which does not depend on poloidal angle." maxlog = 1
     ctrl.inner_model in _VALID_INNER_MODELS || ctrl.inner_model in _UNSUPPORTED_INNER_MODELS ||
         throw(ArgumentError("SLAYERControl: inner_model=$(ctrl.inner_model) " *
                             "not in $(_VALID_INNER_MODELS)"))
     # A disabled section is never solved, so it may keep a GGJ model while serving only its profiles.
     !ctrl.enabled || ctrl.inner_model in _VALID_INNER_MODELS ||
-        throw(ArgumentError("SLAYERControl: inner_model=$(ctrl.inner_model) is not currently supported because GGJ tearing " *
-                            "growth rates are not implemented; use inner_model=slayer_fitzpatrick"))
+        throw(
+            ArgumentError(
+                "SLAYERControl: inner_model=$(ctrl.inner_model) is not currently supported because GGJ tearing " *
+                "growth rates are not implemented; use inner_model=slayer_fitzpatrick"
+            )
+        )
     ctrl.scan_mode in _VALID_SCAN_MODES ||
         throw(ArgumentError("SLAYERControl: scan_mode=$(ctrl.scan_mode) " *
                             "not in $(_VALID_SCAN_MODES)"))

diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl
@@ -114,7 +114,7 @@ const TEARING_H5_ANNOTATIONS = [
         (; long_name="real E×B Doppler offset ΔRe(Q) = −τ_k·n·Ω_E applied to each surface's inner-layer Q in the coupled determinant (0 in uncoupled mode)", dims=("surface",)),
     "PerSurface/tau_R" => (; long_name="resistive diffusion time τ_R = μ₀r_s²/η per surface", units="s", dims=("surface",)),
     "PerSurface/Delta_prime_norm" => (; long_name="Δ'-normalization factor S^(1/3)/r_s per surface", units="1/m", dims=("surface",)),
-    "PerSurface/rs" => (; long_name="minor radius of each rational surface", units="m", dims=("surface",)),
+    "PerSurface/rs" => (; long_name="minor radius √(2ψ_t/B₀) of each rational surface (toroidal-flux label)", units="m", dims=("surface",)),
     "PerSurface/R0" => (; long_name="major radius", units="m", dims=("surface",)),
     "PerSurface/bt" => (; long_name="toroidal field", units="T", dims=("surface",)),
     "PerSurface/sval_r" => (; long_name="r-based magnetic shear r_s·(dq/dr)/q (Fitzpatrick convention)", dims=("surface",)),

diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl
@@ -226,14 +226,14 @@ end
     @test binding([s1]) === s1
 
     # The scan must see the run's own plasma: every [SLAYER] setting the layer width depends on
-    # reaches it, and changing one changes the widths.
+    # reaches it, and changing one changes the widths. theta_sample is passed through but warns: the flux label ignores it.
     trunc_dir = joinpath(dirname(@__DIR__), "examples", "DIIID-like_truncation_example")
     trunc_deck = TOML.parsefile(joinpath(trunc_dir, "gpec.toml"))
     slayer_section = trunc_deck["SLAYER"]
     kin = Equilibrium.read_kinetic_file(joinpath(trunc_dir, trunc_deck["Equilibrium"]["kinetic_file"]))
     slayer_inputs(over) = Dict{String,Any}("SLAYER" => merge(slayer_section, over))
     overlap_inputs(over; kinetic=kin) = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(slayer_inputs(over), kinetic)
-    oi = overlap_inputs(Dict("ion_mass" => 3.0, "bt" => 1.7, "theta_sample" => 0.3, "lnLambda_form" => "sauter"))
+    oi = @test_logs (:warn, r"theta_sample") match_mode = :any overlap_inputs(Dict("ion_mass" => 3.0, "bt" => 1.7, "theta_sample" => 0.3, "lnLambda_form" => "sauter"))
     @test (oi.mu_i, oi.bt, oi.theta, oi.lnLambda_form) == (3.0, 1.7, 0.3, :sauter)
 
     # A disabled section feeds the scan whatever its inner model, and missing kinetic profiles only skip the scan;

diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl
@@ -94,9 +94,10 @@
         @test sl[1].R0 ≈ equil.ro
         @test sl[1].bt == 2.0
 
-        # Minor radius and r-based shear recovered from the equilibrium
-        rs1 = surface_minor_radius(equil, 0.3)
-        da1 = surface_da_dpsi(equil, 0.3)
+        # Minor radius and r-based shear recovered from the equilibrium in the default, toroidal-flux label
+        r_at, dr_at = radial_label(equil; rs_method=:flux)
+        rs1 = r_at(0.3)
+        da1 = dr_at(0.3)
         @test sl[1].rs ≈ rs1
         @test sl[1].sval_r ≈ rs1 * 1.5 / (2.0 * da1)