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                     |
@@ -329,6 +329,19 @@ function slayer_parameters(;
             )
         )
     iota_e = Q_e / Q_e_minus_Q_i
+    # D below takes sqrt(iota_e). iota_e lies in (0,1) whenever ω_*e and ω_*i have
+    # opposite signs — the physical case and the convention stated above. Surface a
+    # negative value as a clear error rather than a bare DomainError.
+    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
@@ -342,7 +355,10 @@ 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 = S^(1/3)·iota_e^(1/2)·(d_beta/r_s), Fitzpatrick's ion-sound-radius normalization.
+    # iota_e is the electron share of the total diamagnetic frequency; it reduces to
+    # 1/(1+tau) only when the electron and ion pressure-gradient scale lengths match.
+    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

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

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,