InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -73,15 +73,17 @@ coupling_mode           = "uncoupled" # "uncoupled" (per-surface) or "coupled" (
 dc_type                 = "none"     # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal"
 mu_i                    = 2.0        # Ion mass in proton-mass units (2.0 = deuterium)
 zeff                    = 1.0        # Effective charge
-chi_perp                = 1.0       # fallback only; the kinetic file supplies χ⊥(ψ)
+chi_perp_e              = 1.0       # fallback only; the kinetic file supplies χ⊥,e(ψ)
+chi_perp_i              = 1.0       # fallback only; the kinetic file supplies χ⊥,i(ψ)
 chi_tor                 = 1.0       # fallback only; the kinetic file supplies χ_φ(ψ)
+P_perp_model            = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s])
 pole_threshold_adaptive = true       # SLAYER |Δ| spans many orders; adapt to 10·median(|Δ|)
 filter_above_poles      = true       # Discard roots above the highest pole γ
 filter_outside_re       = true       # Condition the above-pole filter on the +γ step exiting the Re(Δ)=0 contour loop
 polish_roots            = true       # refine each root to the true zero; enables the validity gate
 store_scan              = false      # per-surface scans omitted from HDF5 (uncoupled has one per surface)
 # Standardized kinetic-profile file (GPEC HDF5 kinetic schema), shared with the
-# sibling ideal example; carries n_e, T_e, T_i, omega_E, chi_e (χ⊥), chi_phi (χ_φ).
+# sibling ideal example; carries n_e, T_e, T_i, omega_E, chi_e (χ⊥,e), chi_i (χ⊥,i), chi_phi (χ_φ).
 profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled
 
 [SLAYER.scan_grid]

src/Equilibrium/KineticProfiles.jl

@@ -66,7 +66,8 @@ temperatures eV, frequencies rad/s, diffusivities m²/s.
 | `T_i`/`T_e`         | ion / electron temperature                                                                                        |
 | `omega_E`           | ExB rotation                                                                                                      |
 | `omega_tor`         | toroidal rotation (optional)                                                                                      |
-| `chi_e`             | perpendicular heat diffusivity χ⊥                                                                                 |
+| `chi_e`             | electron perpendicular heat diffusivity χ⊥,e                                                                      |
+| `chi_i`             | ion perpendicular heat diffusivity χ⊥,i                                                                           |
 | `chi_phi`           | toroidal momentum diffusivity χ_φ                                                                                 |
 | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras |
 | `provenance`        | short string recording the source file/format                                                                     |
@@ -80,25 +81,26 @@ struct KineticProfileData
     omega_E::Union{Nothing,Vector{Float64}}
     omega_tor::Union{Nothing,Vector{Float64}}
     chi_e::Union{Nothing,Vector{Float64}}
+    chi_i::Union{Nothing,Vector{Float64}}
     chi_phi::Union{Nothing,Vector{Float64}}
     species_densities::Union{Nothing,Dict{String,Vector{Float64}}}
     provenance::String
 end
 
 function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing,
-    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing,
+    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_i=nothing, chi_phi=nothing,
     species_densities=nothing, provenance="")
     _v(x) = x === nothing ? nothing : Float64.(collect(x))
     return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e),
-        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), species_densities, String(provenance))
+        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), species_densities, String(provenance))
 end
 
 const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5")
 
 # Dataset units written into / expected from the HDF5 kinetic schema.
 const _KINETIC_H5_UNITS = Dict("psi" => "normalized poloidal flux", "n_i" => "m^-3",
     "n_e" => "m^-3", "T_i" => "eV", "T_e" => "eV", "omega_E" => "rad/s",
-    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_phi" => "m^2/s")
+    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s")
 
 """
     read_kinetic_file(path; group="/") -> KineticProfileData
@@ -137,7 +139,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
         haskey(ats, "provenance") && (prov = string(read(ats["provenance"])))
         # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density
         # profiles for explicit multi-ion input; other names are reserved for future schema fields.
-        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi")
+        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_i", "chi_phi")
         extras = Dict{String,Vector{Float64}}()
         for k in keys(g)
             (k in standard || !startswith(k, "n_")) && continue
@@ -147,7 +149,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
         return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))),
             n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"),
             omega_E=rd("omega_E"), omega_tor=rd("omega_tor"),
-            chi_e=rd("chi_e"), chi_phi=rd("chi_phi"),
+            chi_e=rd("chi_e"), chi_i=rd("chi_i"), chi_phi=rd("chi_phi"),
             species_densities=(isempty(extras) ? nothing : extras), provenance=prov)
     end
 end
@@ -178,6 +180,7 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData;
         put("omega_E", data.omega_E)
         put("omega_tor", data.omega_tor)
         put("chi_e", data.chi_e)
+        put("chi_i", data.chi_i)
         put("chi_phi", data.chi_phi)
         # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip.
         if data.species_densities !== nothing

src/InnerLayer/SLAYER/LayerInputs.jl

@@ -181,8 +181,14 @@ profiles, without an intermediate file round-trip.
 
   - `mu_i`      -- ion mass in proton-mass units (default `2.0` for D).
   - `zeff`      -- effective charge (default `1.0`).
-  - `chi_perp`  -- perpendicular heat diffusivity [m²/s]. Scalar or a
+  - `chi_perp_e`  -- electron perpendicular heat diffusivity [m²/s]. Scalar or a
     callable of `psi` (default `1.0`).
+  - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], scalar or callable of
+    `psi` (default `1.0`); used by `P_perp_model = :chi_perp_i` and `:D_perp`.
+  - `tau_E` -- whole-plasma energy confinement time [s], one scalar for every surface;
+    required by `P_perp_model=:tau_E`.
+  - `P_perp_model` -- how `P_perp` is built; see [`slayer_parameters`](@ref)
+    (default `:chi_perp_e`).
   - `chi_tor`   -- toroidal heat diffusivity [m²/s]. Scalar or a callable
     of `psi` (default `1.0`).
   - `dr_val`    -- resistive interchange index `D_R = E + F + H²`
@@ -232,7 +238,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     mu_i::Real=2.0,
     zeff::Real=1.0,
     z_i::Real=1.0,
-    chi_perp=1.0,
+    chi_perp_e=1.0,
+    chi_perp_i=1.0,
+    tau_E::Real=NaN,
+    P_perp_model::Symbol=:chi_perp_e,
     chi_tor=1.0,
     dr_val=nothing,
     dgeo_val=nothing,
@@ -262,13 +271,11 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     # number enters here. Main-ion density is taken equal to the electron density
     # (quasi-neutrality, matching the staging step).
     chi1 = 2π * equil.psio
+    _kinetic_at(ψ) = (n_e=Float64(profiles.n_e(ψ)), dn_e=Float64(profiles.n_e(ψ; deriv=DerivOp(1))),
+        T_e=Float64(profiles.T_e(ψ)), dT_e=Float64(profiles.T_e(ψ; deriv=DerivOp(1))),
+        T_i=Float64(profiles.T_i(ψ)), dT_i=Float64(profiles.T_i(ψ; deriv=DerivOp(1))))
     _omega_star_at(ψ, n_tor) = begin
-        n_e = Float64(profiles.n_e(ψ))
-        dn_e = Float64(profiles.n_e(ψ; deriv=DerivOp(1)))
-        T_e = Float64(profiles.T_e(ψ))
-        dT_e = Float64(profiles.T_e(ψ; deriv=DerivOp(1)))
-        T_i = Float64(profiles.T_i(ψ))
-        dT_i = Float64(profiles.T_i(ψ; deriv=DerivOp(1)))
+        (; n_e, dn_e, T_e, dT_e, T_i, dT_i) = _kinetic_at(ψ)
         ω_star_e = n_tor * (2π / chi1) * (T_e * dn_e / n_e + dT_e)
         ω_star_i = -n_tor * (2π / (Float64(z_i) * chi1)) * (T_i * dn_e / n_e + dT_i)
         return (ω_star_e, ω_star_i)
@@ -291,6 +298,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
         n_res = sing.n[1]
 
         prof = profiles(psi)
+        kin = _kinetic_at(psi)
         # Take ω_*e, ω_*i from the spline derivatives, or from `profiles` when the caller
         # supplies them directly. `run_slayer` supplies zeros, so the latter is a library path.
         ω_e_use, ω_i_use = compute_omega_star ? _omega_star_at(psi, n_res) : (prof.omega_e, prof.omega_i)
@@ -378,8 +386,12 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use,
             qval=q, sval_r=sval_r, bt=_bt_at(psi),
             rs=rs, R0=R0_use, mu_i=mu_i, zeff=zeff,
-            chi_perp=_eval(chi_perp, psi),
+            chi_perp_e=_eval(chi_perp_e, psi),
             chi_tor=_eval(chi_tor, psi),
+            chi_perp_i=_eval(chi_perp_i, 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,
+            P_perp_model=P_perp_model,
             m=m_res, n=n_res,
             dr_val=dr_val_k,
             dgeo_val=dgeo_val_k,

src/InnerLayer/SLAYER/LayerParameters.jl

@@ -16,7 +16,7 @@
 
 Dimensionless layer-physics parameters at one rational surface for the
 Fitzpatrick two-fluid drift-MHD SLAYER inner-layer model (Fitzpatrick
-2023; Park et al. 2022), plus dimensional auxiliaries required for
+2023; Park 2022), plus dimensional auxiliaries required for
 de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 `D_norm` (not the older `pr`/`pe`/`ds` set).
 
@@ -28,7 +28,7 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 | `lu`       | Lundquist number S = τ_R / τ_H                                    |
 | `c_beta`   | Compressibility √(β_local / (1 + β_local))                        |
 | `D_norm`   | (d_β/r_s) · S^(1/3) · √ι_e  (Fitzpatrick normalized scale)        |
-| `P_perp`   | Perpendicular Prandtl number τ_R / τ_⊥                            |
+| `P_perp`   | Perpendicular Prandtl number τ_R / τ_⊥ (set by `P_perp_model`)    |
 | `P_tor`    | Toroidal-direction Prandtl number τ_R / τ_‖tor                    |
 | `Q_e`      | Normalized electron diamagnetic: −tauk · ω_*e                     |
 | `Q_i`      | Normalized ion diamagnetic:      −tauk · ω_*i                     |
@@ -104,6 +104,9 @@ end
 # `dc_tmp = 0` branch.
 const ALLOWED_DC_TYPES = (:none, :lar, :rfitzp, :toroidal)
 
+# Allowed P_perp_model values selecting how P_perp is built. `:chi_perp_e` is the default.
+const ALLOWED_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta, :tau_E)
+
 """
     r_based_shear(rs, q, dq_dpsi, da_dpsi) -> Float64
 
@@ -132,7 +135,7 @@ end
 # Internal: solve the Wd self-consistency loop for the chi_parallel-based
 # critical Δ (Connor-Hastie-Helander 2015). Returns dc_tmp as a Float64.
 function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
-    chi_perp::Real, t_e::Real, zeff::Real, tau_ee::Real,
+    chi_perp_e::Real, t_e::Real, zeff::Real, tau_ee::Real,
     rs::Real, R0::Real, sval_r::Real, n_tor::Integer,
     max_iter::Integer=100, tol::Real=1e-10)
     dc_type in ALLOWED_DC_TYPES ||
@@ -149,7 +152,7 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
         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 *
+        Wd_new = sqrt(8.0) * (chi_perp_e / chi_par)^0.25 *
                  (1.0 / sqrt((rs / R0) * abs(sval_r) * n_tor))
         if abs(Wd_new - Wd) / max(abs(Wd), 1e-30) < tol
             Wd = Wd_new
@@ -165,22 +168,25 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
 
     if dc_type === :lar
         return 0.5 * (-dr_val) * π^1.5 *
-               (chi_par / chi_perp)^0.25 *
+               (chi_par / chi_perp_e)^0.25 *
                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
         return 0.5 * (-dr_val) * π^1.5 *
-               (chi_par / chi_perp)^0.25 * dgeo_val
+               (chi_par / chi_perp_e)^0.25 * dgeo_val
     end
     return 0.0
 end
 
 """
     slayer_parameters(; n_e, t_e, t_i, omega, omega_e, omega_i,
                         qval, sval_r, bt, rs, R0, mu_i, zeff,
-                        chi_perp, chi_tor,
+                        chi_perp_e, chi_tor,
                         m, n,
+                        chi_perp_i=1.0, tau_E=NaN,
+                        dlnn=NaN, dlnTe=NaN, dlnTi=NaN,
+                        P_perp_model=:chi_perp_e,
                         dr_val=0.0, dgeo_val=0.0,
                         dc_type=:none, ising=0,
                         resistivity_model=SauterNeoModel(),
@@ -210,7 +216,23 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl
   - `R0`      -- major radius [m]
   - `mu_i`    -- ion mass in proton-mass units (e.g. 2.0 for D)
   - `zeff`    -- effective charge
-  - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s]
+  - `chi_perp_e`, `chi_tor` -- electron perpendicular heat / toroidal momentum diffusivity [m²/s]
+  - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s] (used by `P_perp_model=:chi_perp_i`)
+  - `dlnn`, `dlnTe`, `dlnTi` -- radial log-gradients of n, T_e, T_i in any common radial variable
+    (only their ratios enter; required by `P_perp_model=:D_perp`)
+  - `tau_E`   -- whole-plasma energy confinement time [s] (required by `P_perp_model=:tau_E`)
+  - `P_perp_model` -- how `P_perp` is built:
+      + `:chi_perp_e` (default) -- τ_R·χ⊥,e/r_s²
+      + `:chi_perp_i` -- τ_R·χ⊥,i/r_s²
+      + `:P_phi`  -- P_perp = P_tor (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143)
+      + `:D_perp` -- τ_R/τ_⊥ with τ_⊥ = r_s²/D_⊥ and D_⊥ from Eq. 16 of Fitzpatrick, Phys. Plasmas 29,
+        032507 (2022), with η_⊥ = η_∥:
+        D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)[η_e τ_e/(1+η_e) χ⊥,e + η_i τ_i/(1+η_i) χ⊥,i]
+      + `:c_beta` -- P_perp = C² with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, with the
+        thermal-conduction term K dropped). This is the χ⊥ → 0 limit of `:D_perp` (P_⊥ = μ₀D_⊥/η = c_β²), so it is
+        a floor on P_perp; realistic anomalous transport makes D_⊥ much larger.
+      + `:tau_E`  -- P_perp = P_tor = τ_R/τ_E, overriding χ_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131)
+    `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode.
   - `m`, `n`  -- poloidal / toroidal mode numbers at the surface
   - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula
   - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal`
@@ -263,8 +285,11 @@ function slayer_parameters(;
     omega::Real, omega_e::Real, omega_i::Real,
     qval::Real, sval_r::Real, bt::Real,
     rs::Real, R0::Real, mu_i::Real, zeff::Real,
-    chi_perp::Real, chi_tor::Real,
+    chi_perp_e::Real, chi_tor::Real,
     m::Integer, n::Integer,
+    chi_perp_i::Real=1.0, tau_E::Real=NaN,
+    dlnn::Real=NaN, dlnTe::Real=NaN, dlnTi::Real=NaN,
+    P_perp_model::Symbol=:chi_perp_e,
     dr_val::Real=0.0, dgeo_val::Real=0.0,
     dc_type::Symbol=:none, ising::Integer=0,
     resistivity_model::NeoResistivityModel=SauterNeoModel(),
@@ -366,10 +391,41 @@ function slayer_parameters(;
     c_beta = sqrt(lbeta / (1.0 + lbeta))
 
     # Effective Prandtl-like transport ratios
-    tau_perp = rs^2 / chi_perp
-    P_perp = tau_r / tau_perp
     tau_tor = rs^2 / chi_tor
     P_tor = tau_r / tau_tor
+    if P_perp_model === :tau_E
+        isfinite(tau_E) && tau_E > 0 ||
+            throw(ArgumentError("slayer_parameters: P_perp_model=:tau_E needs a positive tau_E, got $tau_E"))
+        # P_⊥ = P_φ = τ_R/τ_E (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131)
+        P_tor = tau_r / tau_E
+    end
+    P_perp = if P_perp_model === :chi_perp_e
+        tau_r / (rs^2 / chi_perp_e)
+    elseif P_perp_model === :chi_perp_i
+        tau_r / (rs^2 / chi_perp_i)
+    elseif P_perp_model === :P_phi || P_perp_model === :tau_E
+        # P_⊥ = P_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143)
+        P_tor
+    elseif P_perp_model === :D_perp
+        # D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥), with
+        # η_s/(1+η_s) = dlnT_s/(dlnn+dlnT_s) and τ = (T_e/T_i)(1+η_e)/(1+η_i) written flat-density safe.
+        frac_e = dlnTe / (dlnn + dlnTe)
+        frac_i = dlnTi / (dlnn + dlnTi)
+        tau_F = (t_e / t_i) * (dlnn + dlnTe) / (dlnn + dlnTi)
+        D_perp = c_beta^2 * eta / MU_0 +
+                 (2 / 3) * (1 - c_beta^2) * (frac_e * tau_F / (1 + tau_F) * chi_perp_e + frac_i / (1 + tau_F) * chi_perp_i)
+        isfinite(D_perp) && D_perp > 0 ||
+            throw(ArgumentError("slayer_parameters: P_perp_model=:D_perp gave D_perp=$D_perp; " *
+                                "check dlnn/dlnTe/dlnTi (got $dlnn, $dlnTe, $dlnTi)"))
+        # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥
+        tau_r / (rs^2 / D_perp)
+    elseif P_perp_model === :c_beta
+        # C² = P_⊥ with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, K = 0)
+        c_beta^2
+    else
+        throw(ArgumentError("slayer_parameters: unknown P_perp_model=$P_perp_model. " *
+                            "Allowed: $(ALLOWED_P_PERP_MODELS)"))
+    end
 
     # Normalized beta-related width and Δ-normalization
     d_beta = c_beta * d_i
@@ -378,7 +434,7 @@ function slayer_parameters(;
 
     # Critical-Δ offset from chi_parallel matching
     dc_tmp = _solve_dc_tmp(; dc_type=dc_type, dr_val=dr_val, dgeo_val=dgeo_val,
-        chi_perp=chi_perp, t_e=t_e, zeff=zeff,
+        chi_perp_e=chi_perp_e, t_e=t_e, zeff=zeff,
         tau_ee=tau_ee, rs=rs, R0=R0, sval_r=sval_r,
         n_tor=n)
 

src/InnerLayer/SLAYER/LayerThickness.jl

@@ -2,7 +2,7 @@
 #
 # Resistive inner-layer *width* (del_s) diagnostic for the SLAYER two-fluid
 # drift-MHD slab layer (Fitzpatrick, Tearing Mode Dynamics in Tokamak
-# Plasmas, IOP 2023; Park et al. 2022, Phys. Plasmas 29, 122505; Burgess
+# Plasmas, IOP 2023; Park 2022, Phys. Plasmas 29, 072506; Burgess
 # et al. 2026). This is a distinct quantity from the matching index
 # Δ = π/W' returned by the dispersion solver (Riccati.jl): it integrates a
 # separate reduced Riccati ODE evaluated at the electron diamagnetic