InnerLayer.SLAYER - BUGFIX - Use the shear magnitude in the layer timescales and widths

slayer_parameters took sval_r with its sign in five places that are magnitudes: the
parallel diffusivity chi_par_lmfp (twice), the island width Wd, the :lar critical-Delta,
and the Alfven time tau_h. On a reverse-shear surface tau_h came out negative, hence a
negative Lundquist number, and tauk = lu^(1/3)*tau_h threw a DomainError before any layer
quantity was produced. The two sqrt() sites would have thrown on the same surface.

Found while scanning DIII-D reconstructions: shot 153072_3415 has q dipping to 1.855
off-axis, giving two q = 2 surfaces, and the inner one is on the negative-shear branch.
Rational-surface finding is already reverse-shear-safe, so the layer build is handed
exactly those surfaces.

The layer timescales and widths depend on |dq/dr|, not its sign -- the module's own
diffusive-resistive width derivation states the convention as (n|s|). abs(s) == s wherever
the shear is positive, so no currently-working surface changes value; this only turns a
crash into a result.

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

Author: d-burg

src/InnerLayer/SLAYER/LayerParameters.jl

@@ -133,11 +133,11 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
     Wd = 0.1
     converged = false
     for _ in 1:max_iter
-        chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * sval_r * Wd)
+        chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * abs(sval_r) * Wd)
         chi_par = (chi_par_smfp * chi_par_lmfp) /
                   (chi_par_smfp + chi_par_lmfp)
         Wd_new = sqrt(8.0) * (chi_perp / chi_par)^0.25 *
-                 (1.0 / sqrt((rs / R0) * sval_r * n_tor))
+                 (1.0 / sqrt((rs / R0) * abs(sval_r) * n_tor))
         if abs(Wd_new - Wd) / max(abs(Wd), 1e-30) < tol
             Wd = Wd_new
             converged = true
@@ -147,13 +147,13 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
     end
     converged || error("SLAYERParameters: Wd iteration failed to converge")
 
-    chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * sval_r * Wd)
+    chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * abs(sval_r) * Wd)
     chi_par = (chi_par_smfp * chi_par_lmfp) / (chi_par_smfp + chi_par_lmfp)
 
     if dc_type === :lar
         return 0.5 * (-dr_val) * π^1.5 *
                (chi_par / chi_perp)^0.25 *
-               sqrt((n_tor * sval_r) / (R0 * rs))
+               sqrt((n_tor * abs(sval_r)) / (R0 * rs))
     elseif dc_type === :rfitzp
         return -(sqrt(2.0) * π^1.5 * dr_val) / Wd
     elseif dc_type === :toroidal
@@ -297,8 +297,10 @@ function slayer_parameters(;
 
     # Alfven time uses minor-radius shear directly (sval enters the
     # b_l = (n/m) r_s sval bt / R0 expression and cancels through to
-    # tau_h = R0 sqrt(mu0 rho) / (n sval bt)).
-    tau_h = R0 * sqrt(MU_0 * rho) / (n * sval_r * bt)
+    # tau_h = R0 sqrt(mu0 rho) / (n sval bt)). Magnitude only: the layer timescales and
+    # widths depend on |dq/dr|, not its sign, and a reverse-shear surface would otherwise
+    # give tau_h < 0, hence lu < 0 and a DomainError in lu^(1/3) below.
+    tau_h = R0 * sqrt(MU_0 * rho) / (n * abs(sval_r) * bt)
     # Resistive diffusion time τ_R = μ₀ r_s² / η (Fitzpatrick 2023), with
     # the selected η closure — neoclassical by default — setting the
     # Lundquist number.