InnerLayer.SLAYER - BUGFIX! - Normalize D by the electron diamagnetic share, not the temperature ratio

src/InnerLayer/SLAYER/LayerParameters.jl

@@ -27,7 +27,7 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 | `tau`      | T_i / T_e                                                         |
 | `lu`       | Lundquist number S = τ_R / τ_H                                    |
 | `c_beta`   | Compressibility √(β_local / (1 + β_local))                        |
-| `D_norm`   | (d_β/r_s) · S^(1/3) · √(τ/(1+τ))  (Fitzpatrick normalized scale)  |
+| `D_norm`   | (d_β/r_s) · S^(1/3) · √ι_e  (Fitzpatrick normalized scale)        |
 | `P_perp`   | Perpendicular Prandtl number τ_R / τ_⊥                            |
 | `P_tor`    | Toroidal-direction Prandtl number τ_R / τ_‖tor                    |
 | `Q_e`      | Normalized electron diamagnetic: −tauk · ω_*e                     |
@@ -350,6 +350,16 @@ function slayer_parameters(;
             )
         )
     iota_e = Q_e / Q_e_minus_Q_i
+    iota_e > 0 ||
+        throw(
+            ArgumentError(
+                "slayer_parameters: iota_e = Q_e/(Q_e - Q_i) = $iota_e is " *
+                "not positive, so the ion-sound-radius normalization " *
+                "D_norm = (d_beta/r_s)·S^(1/3)·sqrt(iota_e) is undefined. " *
+                "Check the diamagnetic-frequency inputs ω_*e, ω_*i — they " *
+                "are expected to carry opposite signs."
+            )
+        )
 
     # Plasma beta and compressibility
     lbeta = (5.0 / 3.0) * MU_0 * n_e * E_CHG * (t_e + t_i) / bt^2
@@ -363,7 +373,7 @@ function slayer_parameters(;
 
     # Normalized beta-related width and Δ-normalization
     d_beta = c_beta * d_i
-    D_norm = (d_beta / rs) * lu^(1.0 / 3.0) * sqrt(tau / (1.0 + tau))
+    D_norm = (d_beta / rs) * lu^(1.0 / 3.0) * sqrt(iota_e)
     delta_n = lu^(1.0 / 3.0) / rs
 
     # Critical-Δ offset from chi_parallel matching

src/Tearing/Runner/HDF5Output.jl

@@ -74,7 +74,7 @@ const TEARING_H5_ANNOTATIONS = [
     "PerSurface/tau" => (; long_name="temperature ratio τ = T_i/T_e per surface", dims=("surface",)),
     "PerSurface/lu" => (; long_name="Lundquist number S per surface", dims=("surface",)),
     "PerSurface/c_beta" => (; long_name="compressibility factor c_β = √(β_local/(1+β_local)) per surface", dims=("surface",)),
-    "PerSurface/D_norm" => (; long_name="Fitzpatrick normalized ion-sound/drift scale D = (d_β/r_s)·S^(1/3)·√(τ/(1+τ)) per surface", dims=("surface",)),
+    "PerSurface/D_norm" => (; long_name="Fitzpatrick normalized ion-sound/drift scale D = (d_β/r_s)·S^(1/3)·√ι_e per surface", dims=("surface",)),
     "PerSurface/P_perp" => (; long_name="perpendicular magnetic Prandtl number per surface", dims=("surface",)),
     "PerSurface/P_tor" => (; long_name="toroidal (momentum) magnetic Prandtl number per surface", dims=("surface",)),
     "PerSurface/Q_e" => (; long_name="normalized electron diamagnetic frequency Q_e per surface", dims=("surface",)),

test/runtests_dispersion_coupled.jl

@@ -21,7 +21,7 @@
     function _slayer_ref()
         return slayer_parameters(
             n_e=5.0e19, t_e=1000.0, t_i=1000.0,
-            omega=0.0, omega_e=1.0e4, omega_i=5.0e3,
+            omega=0.0, omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
             chi_perp=1.0, chi_tor=1.0, m=2, n=1)

test/runtests_dispersion_residual.jl

@@ -21,7 +21,7 @@
     function _slayer_ref()
         return slayer_parameters(
             n_e=5.0e19, t_e=1000.0, t_i=1000.0,
-            omega=0.0, omega_e=1.0e4, omega_i=5.0e3,
+            omega=0.0, omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
             chi_perp=1.0, chi_tor=1.0, m=2, n=1)

test/runtests_slayer_params.jl

@@ -9,7 +9,7 @@
     function _ref_kwargs(; dr_val=0.0, dc_type=:none)
         return (
             n_e=5.0e19, t_e=1000.0, t_i=1000.0,
-            omega=0.0, omega_e=1.0e4, omega_i=5.0e3,
+            omega=0.0, omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
             chi_perp=1.0, chi_tor=1.0,
@@ -36,13 +36,13 @@
 
         # Trivially exact ratios
         @test p.tau ≈ 1.0
-        # Q_e = −tauk·1e4 = negative; Q_i = −tauk·5e3 = negative
-        # Q_e − Q_i = −tauk·5e3 = Q_i (since Q_e = 2·Q_i) ⇒ iota_e = Q_e/Q_i = 2
-        @test p.iota_e ≈ 2.0
+        # Physical opposite-drift inputs: ω_*e < 0 < ω_*i, so Q = −tauk·ω gives
+        # Q_e > 0 > Q_i and iota_e = Q_e/(Q_e − Q_i) = 1e4/(1e4 + 5e3) = 2/3
+        @test p.iota_e ≈ 2 / 3
 
-        # Sign convention check (SLAYER layerinputs)
-        @test p.Q_e == -p.tauk * 1.0e4
-        @test p.Q_i == -p.tauk * 5.0e3    # SLAYER params convention: Q_i = −tauk·ω*i
+        # Sign convention check (SLAYER layerinputs): Q = −tauk·ω
+        @test p.Q_e == p.tauk * 1.0e4
+        @test p.Q_i == -p.tauk * 5.0e3
 
         # Default resistivity closure is neoclassical (Sauter F_33). The
         # trapped-particle correction raises η above plain Spitzer at the
@@ -83,8 +83,9 @@
         @test p.P_perp ≈ p.P_tor
         @test p.P_perp > 0
 
-        # D_norm = (d_β/r_s)·S^(1/3)·√(τ/(1+τ))
-        D_norm_expected = (p.d_beta / p.rs) * p.lu^(1 / 3) * sqrt(p.tau / (1 + p.tau))
+        # D_norm = (d_β/r_s)·S^(1/3)·√ι_e — the electron share of the total
+        # diamagnetic frequency, not the temperature ratio.
+        D_norm_expected = (p.d_beta / p.rs) * p.lu^(1 / 3) * sqrt(p.iota_e)
         @test p.D_norm ≈ D_norm_expected rtol = 1e-12
 
         # delta_n = S^(1/3)/r_s
@@ -139,6 +140,14 @@
         @test p.ising == 0
     end
 
+    @testset "Test 1d: reject diamagnetic frequencies that leave D undefined" begin
+        # iota_e = ω_*e/(ω_*e − ω_*i) must be finite and positive for D ∝ √iota_e.
+        _drifts(ωe, ωi) = merge(_ref_kwargs(), (omega_e=ωe, omega_i=ωi))
+        @test_throws ArgumentError slayer_parameters(; _drifts(-1.0e4, -1.0e4)...)   # degenerate: iota_e singular
+        @test_throws ArgumentError slayer_parameters(; _drifts(-1.0e4, -2.0e4)...)   # same sign, |ω_*i| > |ω_*e|: iota_e = -1
+        @test_throws ArgumentError slayer_parameters(; _drifts(0.0, 5.0e3)...)       # no electron drift: iota_e = 0
+    end
+
     @testset "Test 2: r-based shear conversion" begin
         # Direct application of r_s · (dq/dψ) / (q · da/dψ).
         @test r_based_shear(0.5, 2.0, 4.0, 0.5) ≈ 2.0

test/runtests_slayer_riccati.jl

@@ -15,7 +15,7 @@
     function _ref_params_large_D()
         return slayer_parameters(;
             n_e=5.0e19, t_e=3000.0, t_i=3000.0,
-            omega=0.0, omega_e=1.0e4, omega_i=5.0e3,
+            omega=0.0, omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
             chi_perp=1.0, chi_tor=1.0,