InnerLayer.SLAYER - FEATURE! - Optional chi_perp_Dc sets the critical-Δ electron diffusivity apart from P_perp

docs/src/stability.md

@@ -349,7 +349,11 @@ total pressure; under the closure `T_s = T_s(n)` it is set by the perpendicular
 `D_⊥` of Fitzpatrick, Phys. Plasmas **29**, 032507 (2022), Eq. 16, which weights the electron and ion energy diffusivities by
 `n T_s′/p′`. The critical-Δ island width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature
 flattening against an electron parallel conductivity, so its χ_⊥ is the electron channel `chi_e`
-in every case. `[SLAYER] P_perp_model` selects how `P_perp` is built: from `D_⊥` (`"D_perp"`,
+in every case. `[SLAYER] chi_perp_Dc` (optional) sets that electron diffusivity for the critical Δ
+alone, leaving `P_perp` on its own χ inputs. It is for a run whose `P_perp` uses a diffusivity
+the critical Δ should not, for example a one-fluid χ_eff in both channels of `D_⊥` when an
+interpretive power balance splits the conducted heat between species unreliably, while the
+critical Δ takes the electron χ that balance allows. `[SLAYER] P_perp_model` selects how `P_perp` is built: from `D_⊥` (`"D_perp"`,
 the default), from `chi_e` alone (the single-χ convention of Fitzpatrick's TJ code), from
 `chi_i`, from the toroidal Prandtl number, from its `c_β²` floor, or from a whole-plasma
 confinement time (`"tau_E"`); see `slayer_parameters`. `"tau_E"` is for a plasma with no χ

src/InnerLayer/SLAYER/LayerInputs.jl

@@ -255,6 +255,8 @@ profiles, without an intermediate file round-trip.
     callable of `psi` (default `1.0`).
   - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], scalar or callable of `psi`;
     required by `P_perp_model = :chi_perp_i` and `:D_perp` (default `nothing`).
+  - `chi_perp_Dc` -- electron perpendicular heat diffusivity [m²/s] used only by the critical Δ,
+    scalar or callable of `psi`; `nothing` (default) uses `chi_perp_e`. See [`slayer_parameters`](@ref).
   - `tau_E` -- whole-plasma energy confinement time τ̄_E [s], required by `P_perp_model = :tau_E`
     (default `nothing`). That mode takes the minor radius `a` of the ψ = 1 surface from the
     `:flux` label, whatever `rs_method` is.
@@ -323,6 +325,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     z_i::Real=1.0,
     chi_perp_e=1.0,
     chi_perp_i=nothing,
+    chi_perp_Dc=nothing,
     tau_E::Union{Real,Nothing}=nothing,
     P_perp_model::Symbol=:D_perp,
     chi_tor=1.0,
@@ -484,6 +487,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             chi_perp_e=_eval(chi_perp_e, psi),
             chi_tor=_eval(chi_tor, psi),
             chi_perp_i=_eval(chi_perp_i, psi),
+            chi_perp_Dc=_eval(chi_perp_Dc, psi),
             dlnn=kin.dn_e / kin.n_e, dlnTe=kin.dT_e / kin.T_e, dlnTi=kin.dT_i / kin.T_i,
             tau_E=tau_E, a_minor=a_minor,
             P_perp_model=P_perp_model,

src/InnerLayer/SLAYER/LayerParameters.jl

@@ -202,7 +202,7 @@ end
                         m, n,
                         chi_perp_i=nothing, tau_E=nothing, a_minor=nothing,
                         dlnn=nothing, dlnTe=nothing, dlnTi=nothing,
-                        P_perp_model=:D_perp,
+                        P_perp_model=:D_perp, chi_perp_Dc=nothing,
                         dr_val=0.0, dgeo_val=0.0, kpar_val=nothing,
                         dc_type=:none, ising=0,
                         resistivity_model=SauterNeoModel(),
@@ -253,7 +253,13 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl
         a uniform χ_E = κ_e a²/τ̄_E with κ_e = 3/16, the value for T_e ∝ 1 − r̂² (Eq. 38). Other
         profile shapes change κ_e by a factor of order 2-3. Both Prandtl numbers come from χ_E, so
         `chi_tor` is ignored in this mode
-    `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode.
+    `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode, unless `chi_perp_Dc` is given.
+  - `chi_perp_Dc` -- electron perpendicular heat diffusivity [m²/s] used **only** by the critical Δ
+    (the W_d of the χ_∥ matching and the `(χ_∥/χ⊥)^(1/4)` factors of `:lar`/`:toroidal`). Default
+    `nothing` uses `chi_perp_e`, so the critical Δ and P_perp see one χ⊥,e as before. Set it when
+    the layer's P_perp is built from a different diffusivity than the electron channel the
+    critical Δ needs -- e.g. a one-fluid χ_eff in `:D_perp` when an interpretive power balance
+    splits the conducted heat between species unreliably. Must be positive and finite when given
   - `m`, `n`  -- poloidal / toroidal mode numbers at the surface
   - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula: the resistive
     interchange index `D_R` and the Connor et al. 2015 Eq. 59 geometric factor
@@ -314,7 +320,7 @@ function slayer_parameters(;
     m::Integer, n::Integer,
     chi_perp_i::Union{Real,Nothing}=nothing, tau_E::Union{Real,Nothing}=nothing, a_minor::Union{Real,Nothing}=nothing,
     dlnn::Union{Real,Nothing}=nothing, dlnTe::Union{Real,Nothing}=nothing, dlnTi::Union{Real,Nothing}=nothing,
-    P_perp_model::Symbol=:D_perp,
+    P_perp_model::Symbol=:D_perp, chi_perp_Dc::Union{Real,Nothing}=nothing,
     dr_val::Real=0.0, dgeo_val::Real=0.0,
     kpar_val::Union{Real,Nothing}=nothing,
     dc_type::Symbol=:none, ising::Integer=0,
@@ -422,6 +428,8 @@ function slayer_parameters(;
     _need("chi_perp_e", chi_perp_e)
     P_perp_model === :tau_E || _need("chi_tor", chi_tor)
     P_perp_model in (:chi_perp_i, :D_perp) && _need("chi_perp_i", chi_perp_i)
+    chi_perp_Dc === nothing || (isfinite(chi_perp_Dc) && chi_perp_Dc > 0) ||
+        throw(ArgumentError("slayer_parameters: chi_perp_Dc must be positive and finite when given; got $chi_perp_Dc"))
     tau_tor = rs^2 / chi_tor
     P_tor = tau_r / tau_tor
     if P_perp_model === :tau_E
@@ -462,9 +470,9 @@ function slayer_parameters(;
     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
+    # Critical-Δ offset from chi_parallel matching; its χ⊥,e is chi_perp_Dc when given
     dc_tmp = _solve_dc_tmp(; dc_type=dc_type, dr_val=dr_val, dgeo_val=dgeo_val,
-        chi_perp_e=chi_perp_e, t_e=t_e, zeff=zeff,
+        chi_perp_e=chi_perp_Dc === nothing ? chi_perp_e : chi_perp_Dc, t_e=t_e, zeff=zeff,
         tau_ee=tau_ee, rs=rs, R0=R0, sval_r=sval_r,
         n_tor=n, kpar_val=kpar_val)
 

src/Tearing/Runner/Control.jl

@@ -95,6 +95,11 @@ constructor.
     χ⊥,i), `:chi_perp_e`, `:chi_perp_i`, `:P_phi`, `:c_beta` or `:tau_E`; see
     `InnerLayer.SLAYER.slayer_parameters`. `:tau_E` sets both P_perp and P_tor from τ̄_E, so it
     ignores `chi_phi`/`chi_tor`. `chi_perp_e` still sets the critical-Δ width in every mode
+    unless `chi_perp_Dc` is set
+  - `chi_perp_Dc` -- electron perpendicular heat diffusivity [m²/s] used only by the critical Δ
+    (default `nothing`: the critical Δ uses the same χ⊥,e as P_perp, from the kinetic file's
+    `chi_e` or `chi_perp_e`). When set, this scalar replaces χ⊥,e in the critical Δ only;
+    P_perp keeps the file's profiles or the `chi_perp_e`/`chi_perp_i` fallbacks. Must be positive
   - `D_R`, `D_geo`  -- critical-Δ formula inputs. `nothing` (default)
     auto-derives them from the equilibrium: `D_R` from the resistive
     interchange index `D_R = E + F + H²` at each surface, `D_geo` from the
@@ -178,6 +183,7 @@ file can keep the χ keys set to 0 to defer to the scalars).
     impurity_charge::Float64 = 6.0
     chi_perp_e::Float64 = 1.0
     chi_perp_i::Union{Float64,Nothing} = nothing
+    chi_perp_Dc::Union{Float64,Nothing} = nothing
     chi_tor::Float64 = 1.0
     tau_E::Union{Float64,Nothing} = nothing
     P_perp_model::Symbol = :D_perp
@@ -269,7 +275,8 @@ function validate(ctrl::SLAYERControl)
     ctrl.P_perp_model in _VALID_P_PERP_MODELS ||
         throw(ArgumentError("SLAYERControl: P_perp_model=$(ctrl.P_perp_model) " *
                             "not in $(_VALID_P_PERP_MODELS)"))
-    for (name, v) in (("chi_perp_e", ctrl.chi_perp_e), ("chi_perp_i", ctrl.chi_perp_i), ("chi_tor", ctrl.chi_tor))
+    for (name, v) in (("chi_perp_e", ctrl.chi_perp_e), ("chi_perp_i", ctrl.chi_perp_i), ("chi_perp_Dc", ctrl.chi_perp_Dc),
+        ("chi_tor", ctrl.chi_tor))
         v === nothing || (isfinite(v) && v > 0) ||
             throw(ArgumentError("SLAYERControl: $name must be positive and finite [m²/s]; got $v"))
     end
@@ -358,7 +365,7 @@ function slayer_control_from_toml(section::AbstractDict)
             kwargs[sym] = v isa Symbol ? v : Symbol(String(v))
         elseif sym in (:Q_re_range, :Q_im_range)
             kwargs[sym] = _as_range(v)
-        elseif sym in (:bt, :D_R, :D_geo, :tau_E, :chi_perp_i)
+        elseif sym in (:bt, :D_R, :D_geo, :tau_E, :chi_perp_i, :chi_perp_Dc)
             # Allow explicit nothing (auto-derive) or a number (override)
             kwargs[sym] = v === nothing ? nothing : Float64(v)
         elseif sym === :omega_E_kHz

src/Tearing/Runner/TearingProblem.jl

@@ -51,6 +51,7 @@ function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control)
         zeff=zeff,
         chi_perp_e=chi.chi_perp_e,
         chi_perp_i=chi.chi_perp_i,
+        chi_perp_Dc=chi.chi_perp_Dc,
         chi_tor=chi.chi_tor,
         tau_E=control.tau_E,
         P_perp_model=control.P_perp_model,

src/Tearing/Runner/run_slayer.jl

@@ -85,11 +85,14 @@ function _slayer_chi_inputs(loaded, control::SLAYERControl; warn::Bool=true)
         uses_tor && loaded.chi_tor === nothing && push!(fallbacks, "chi_phi → chi_tor = $(control.chi_tor)")
         isempty(fallbacks) || @warn "SLAYER: no usable kinetic-file profile (dataset absent or all-zero) for " *
                                     "$(join(fallbacks, ", ")); using the scalar fallback."
+        control.chi_perp_Dc !== nothing && control.D_c_type !== :none &&
+            @info "SLAYER: the critical Δ uses chi_perp_Dc = $(control.chi_perp_Dc) m²/s for χ⊥,e; P_perp keeps its own χ inputs."
         model === :tau_E && @info "SLAYER: P_perp_model = tau_E sets P_perp and P_tor from a uniform " *
               "χ_E = (3/16)a²/τ̄_E with τ̄_E = $(control.tau_E) s; the χ_φ profile is not used."
     end
     return (; chi_perp_e=pick(loaded.chi_perp_e, control.chi_perp_e),
         chi_perp_i=pick(loaded.chi_perp_i, control.chi_perp_i),
+        chi_perp_Dc=control.chi_perp_Dc,
         chi_tor=pick(loaded.chi_tor, control.chi_tor))
 end
 

test/runtests_slayer_params.jl

@@ -282,4 +282,23 @@
         @test cb.dc_tmp == slayer_parameters(; P_perp_model=:chi_perp_e, merge(base, (dr_val=-1.0, dc_type=:lar))...).dc_tmp
         @test cb.P_tor == ref.P_tor
     end
+
+    @testset "Test 4b: chi_perp_Dc sets only the critical Δ" begin
+        grads = (dlnn=-1.0, dlnTe=-2.0, dlnTi=-3.0)
+        for dc in (:lar, :fitzpatrick, :toroidal)
+            base = merge(_ref_kwargs(; dr_val=0.01, dc_type=dc), grads, (chi_perp_i=2.0, P_perp_model=:D_perp))
+            a = slayer_parameters(; base...)
+            b = slayer_parameters(; base..., chi_perp_Dc=3.0)
+            # the critical Δ of b is the one chi_perp_e = 3 would give ...
+            @test b.dc_tmp == slayer_parameters(; merge(base, (chi_perp_e=3.0,))...).dc_tmp
+            @test b.dc_tmp != a.dc_tmp
+            # ... and every other layer quantity, P_perp included, is a's
+            for f in (:tau, :lu, :c_beta, :D_norm, :P_perp, :P_tor, :Q_e, :Q_i, :iota_e, :tauk, :tau_r, :delta_n, :eta, :d_beta)
+                @test getfield(b, f) == getfield(a, f)
+            end
+        end
+        base = merge(_ref_kwargs(; dr_val=0.01, dc_type=:fitzpatrick), (P_perp_model=:chi_perp_e,))
+        @test_throws "chi_perp_Dc must be positive" slayer_parameters(; base..., chi_perp_Dc=0.0)
+        @test_throws "chi_perp_Dc must be positive" slayer_parameters(; base..., chi_perp_Dc=Inf)
+    end
 end

test/runtests_slayer_runner.jl

@@ -116,6 +116,12 @@
         # Scalar diffusivities must be positive.
         @test_throws "chi_perp_i must be positive" Runner.validate(SLAYERControl(; chi_perp_i=-0.5))
         @test_throws "chi_perp_e must be positive" Runner.validate(SLAYERControl(; chi_perp_e=0.0))
+        # The critical-Δ diffusivity is optional, read from the deck, validated, and forwarded unchanged.
+        c_dc = slayer_control_from_toml(merge(section, Dict{String,Any}("chi_perp_Dc" => 1.4)))
+        @test c_dc.chi_perp_Dc == 1.4 && Runner.validate(c_dc) === c_dc
+        @test_throws "chi_perp_Dc must be positive" Runner.validate(SLAYERControl(; chi_perp_Dc=-1.0))
+        @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; chi_perp_i=0.7, chi_perp_Dc=1.4); warn=false).chi_perp_Dc == 1.4
+        @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; chi_perp_i=0.7); warn=false).chi_perp_Dc === nothing
     end
 
     @testset "TearingProblem solve forwards surfaces and Δ'" begin