InnerLayer.SLAYER - FEATURE! - Build P_perp from the perpendicular particle diffusivity by default and add an ion heat diffusivity profile

docs/development/references.md

@@ -76,6 +76,16 @@ The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are
   - Location: `docs/resources/2026-Burgess-Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer.pdf`
   - Describes: The two-fluid slab layer, its four regimes, and the `S^(1/3)` outer-inner matching
 
+- **Fitzpatrick (2022)**: "Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas"
+  - Location: `docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf`
+  - Published: Physics of Plasmas **29**, 032507 (2022)
+  - Describes: The perpendicular particle diffusivity D_⊥ built from the electron and ion energy diffusivities (Eqs. 4-6, 16), used by `P_perp_model = :D_perp`
+
+- **Fitzpatrick (2023)**: "Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas"
+  - Location: `docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf`
+  - Published: Physics of Plasmas **30**, 092512 (2023)
+  - Describes: P_⊥ = P_φ = τ_R/τ_E with the local τ_E = r²/χ_E (Eqs. 81-83, 131, 143) and the whole-plasma τ̄_E = κ_e a²/χ_E (Eqs. 37-38), used by `P_perp_model = :P_phi` and `:tau_E`
+
 - **Fitzpatrick (2025)**: "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation"
   - Location: `docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf`
   - Published: Nuclear Fusion (2025), doi 10.1088/1741-4326/ae4fdd (open access: arXiv:2511.07666)

docs/src/kinetic_forces.md

@@ -12,8 +12,9 @@ extension:
   - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets
     at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`,
     `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional
-    `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``),
-    `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``), and `zeff`
+    `omega_tor`, `chi_e` and `chi_i` (electron and ion perpendicular heat diffusivities
+    ``\chi_{\perp,e}``, ``\chi_{\perp,i}``), `chi_phi` (toroidal momentum diffusivity
+    ``\chi_\phi``), and `zeff`
     (effective charge ``Z_\mathrm{eff}``, dimensionless). Each dataset
     carries a `units` attribute; the root carries `schema_version` and
     `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
@@ -25,7 +26,7 @@ extension:
     `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
     compatibility. Header rows are skipped.
 
-The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`/`zeff`,
+The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_i`/`chi_phi`/`zeff`,
 when present, are carried for the resistive-layer (SLAYER) analysis and ignored
 here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file
 has one, otherwise from `n_i/n_e` with the `[SLAYER]` `impurity_charge`, otherwise

docs/src/stability.md

@@ -333,28 +333,49 @@ resolve it from, so the cap cannot be applied there and setting the flag warns i
 `psilim_from_layer_overlap = true`. It then caps `qlim` before the `rational_offset_fraction` step, so `rational_offset_fraction`
 still selects the final surface from inside the cap. A cap beyond `psihigh` has no effect.
 
-The flag is off by default, so every shipped deck is unchanged by it. With the shipped
-DIII-D-like profiles the overlap point is ψ ≈ 0.9985: outside the ideal and Riccati decks'
-`psihigh = 0.995`, where the cap could not bind even if enabled, but inside the SLAYER deck's
-0.9995, where enabling it would exclude the q = 7 surface. `examples/DIIID-like_truncation_example`
-raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and q = 8, and the resulting
-domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen `psihigh = 0.995` of the
-ideal deck. Whenever `psilim` lands past the overlap point, with the flag on or off, GPEC warns;
-on the shipped decks that includes the SLAYER deck.
-
-**Which perpendicular diffusivity.** The two formulas that take `chi_perp` want different
-quantities. The layer width goes through `P_perp = τ_R/τ_E` with `τ_E = r²/χ_E`, and χ_E is the
-one-fluid perpendicular *energy* diffusivity, the gradient-weighted combination of the electron and
-ion channels, `χ_E = (n_e χ_e ∇T_e + n_i χ_i ∇T_i)/(n_e ∇T_e + n_i ∇T_i)`. The critical-Δ island
-width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature flattening against an electron parallel
-conductivity, so its χ_⊥ is the electron channel. GPEC currently feeds both from one profile, the
-kinetic file's `chi_e`, which in the shipped DIII-D-like data is the electron power-balance
-diffusivity: the right channel for the critical Δ, and a proxy for χ_E in the layer width. This is
-the same single-χ convention as Fitzpatrick's TJ code, which also assumes `T_e = T_i`.
-Separating the two is deferred: it needs an ion `chi_i` in the kinetic-profile schema, and the
-power-balance ion diffusivity available upstream is not yet usable (it comes back negative across
-most of the radius on the cases checked). Since both widths scale as `χ^{1/4}`, a factor-2 error in
-χ moves them by about 19%, but that can be enough to move the overlap cut by a surface.
+The flag is off by default, so every shipped deck is unchanged by it. With the default
+`P_perp_model = "D_perp"` the layers on the shipped DIII-D-like profiles do not overlap inside
+ψ = 1. With the electron-channel `"chi_perp_e"` they are about 40% wider and the overlap point is
+ψ ≈ 0.9985: outside the ideal and Riccati decks' `psihigh = 0.995`, but inside the SLAYER deck's
+0.9995, where the cap would exclude the q = 7 surface. `examples/DIIID-like_truncation_example`
+pins `"chi_perp_e"` and raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and
+q = 8, and the resulting domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen
+`psihigh = 0.995` of the ideal deck. Whenever `psilim` lands past the overlap point, with the flag
+on or off, GPEC warns.
+
+**Which perpendicular diffusivity.** The two formulas that take a perpendicular diffusivity want
+different quantities. The layer width goes through `P_perp`, which multiplies the diffusion of the
+total pressure; under the closure `T_s = T_s(n)` it is set by the perpendicular particle diffusivity
+`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"`,
+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 χ
+profiles: its uniform χ_E sets both `P_perp` and the toroidal `P_tor`, so the file's `chi_phi`
+and `[SLAYER] chi_tor` are not used in that mode. On the DIII-D-like
+profiles `D_⊥` is 0.22–0.32 of `chi_e` at the outer rational surfaces, because the weights sum to
+less than one where density gradients carry part of `p′`. Since the widths scale as `χ^{1/4}`, a
+factor-2 change in χ moves them by about 19%, which can be enough to move the overlap cut by a
+surface.
+
+**The shipped diffusivities.** The DIII-D-like kinetic file carries analytic χ profiles in the
+normalized poloidal flux ψ, in m²/s:
+
+```
+chi_e   = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²)
+chi_i   = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²)
+chi_phi = 3/2 − (5/4) (1 − ψ)⁴
+```
+
+They are not fits to a measured discharge. Their magnitudes are of the same order as the TRANSP
+power-balance χ_e and χ_i reported for DIII-D H-modes by Hong *et al.*, Phys. Plasmas **30**,
+072512 (2023), Fig. 5; `chi_e` includes a factor-2 raise so that the default deck exercises the
+resistive-layer overlap cut. `chi_phi` is set to track `chi_i` (Prandtl number of order one, rising
+outward), as reported for DIII-D by Yoshida *et al.*, Nucl. Fusion **52**, 123005 (2012). The forms
+do not resolve the H-mode pedestal transport barrier, where power-balance χ falls to about
+0.1–0.3 m²/s (Callen *et al.*, Nucl. Fusion **50**, 064004 (2010)).
 
 | Dataset | Contents |
 |---|---|

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -82,18 +82,19 @@ scan_mode               = "amr"      # Q-plane scan strategy: "amr" (adaptive re
 coupling_mode           = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant)
 D_c_type                = "none"     # Critical-Δ offset selector: "none", "lar", "fitzpatrick", or "toroidal"
 ion_mass                = 2.0        # Ion mass in proton-mass units (2.0 = deuterium)
-zeff                    = 1.0        # fallback only; Z_eff(ψ) comes from the kinetic file (its zeff profile, else n_i/n_e with zimp)
+zeff                    = 1.0        # fallback only; Z_eff(ψ) comes from the kinetic file (its zeff profile, else n_i/n_e with impurity_charge)
 impurity_charge         = 6.0        # Impurity charge closing quasi-neutrality when Z_eff(ψ) is derived from n_i/n_e (6.0 = carbon)
-chi_perp                = 1.0       # fallback only; the kinetic file supplies χ⊥(ψ)
+chi_perp_e              = 1.0       # fallback only; the kinetic file supplies χ⊥,e(ψ)
 chi_tor                 = 1.0       # fallback only; the kinetic file supplies χ_φ(ψ)
+P_perp_model            = "D_perp"   # Perpendicular Prandtl number model: "D_perp" (default), "chi_perp_e", "chi_perp_i", "P_phi", "c_beta", or "tau_E" (needs tau_E [s]; also replaces the χ_φ-based P_tor)
 # omega_E_kHz           = {"2/1" = 0.0, "3/1" = 3.0} # override of the E×B rotation Ω_E/2π per unit n [kHz], keyed by m/n; unlisted surfaces use the kinetic file
 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]

examples/DIIID-like_ideal_example/gpec.toml

@@ -88,7 +88,7 @@ write_outputs_to_HDF5 = true            # Write ErrorFields/CoilSensitivities/ t
 verbose = true                          # Log per-coil-set progress and linearity diagnostics
 
 [KineticForces]
-kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5"  # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_phi). Legacy .gpeckf/.kin still readable.
+kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5"  # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_i,chi_phi). Legacy .gpeckf/.kin still readable.
 nn = 1                                  # Toroidal mode number (matches nn_low/nn_high)
 bounce_harmonic_max = 6                 # Bounce harmonics retained
 ion_charge = 1                          # Main-ion charge (deuterium)

examples/DIIID-like_truncation_example/gpec.toml

@@ -44,6 +44,8 @@ psilim_from_layer_overlap = true   # Cap psilim where adjacent resistive layers
 
 [SLAYER]
 # The growth-rate analysis is not run here; the section is present because the layer-overlap scan
-# takes its kinetic profiles, and its chi_perp(psi), from the SLAYER profile_file.
+# takes its kinetic profiles, and its chi_perp_e(psi), from the SLAYER profile_file. P_perp is pinned to the
+# electron channel so the cap binds: with the default D_perp the layers are ~30% narrower and do not overlap inside psi = 1.
 enabled = false                # Run the SLAYER tearing-mode analysis
+P_perp_model = "chi_perp_e"    # Perpendicular Prandtl number model (pinned; see above)
 profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled

regression-harness/cases/diiid_slayer_n1_chi_perp_e.toml

@@ -0,0 +1,10 @@
+# Regression case: the diiid_slayer_n1 deck with P_perp built from the electron diffusivity chi_e alone,
+# the model GPEC used before the D_perp default. It keeps a baseline for the non-default P_perp models
+# and for the electron-channel convention of Fitzpatrick's TJ code.
+[case]
+name = "diiid_slayer_n1_chi_perp_e"
+base = "diiid_slayer_n1"   # Inherit example_dir, [quantities.*] and [overrides] from this case
+description = "DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis with P_perp_model = chi_perp_e"
+
+[overrides]
+"SLAYER.P_perp_model" = "chi_perp_e"   # P_perp from chi_e alone, as before the D_perp default

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 χ⊥ (shipped data: electron power-balance χ_e; see the stability docs)              |
+| `chi_e`             | electron perpendicular heat diffusivity χ⊥,e                                                                      |
+| `chi_i`             | ion perpendicular heat diffusivity χ⊥,i                                                                           |
 | `chi_phi`           | toroidal momentum diffusivity χ_φ                                                                                 |
 | `zeff`              | effective charge Z_eff(ψ) supplied directly (optional); the resistive-layer analysis prefers it to its closure    |
 | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras |
@@ -81,26 +82,27 @@ 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}}
     zeff::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, zeff=nothing,
+    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_i=nothing, chi_phi=nothing, zeff=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), _v(zeff), species_densities, String(provenance))
+        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), _v(zeff), 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", "zeff" => "dimensionless")
+    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s", "zeff" => "dimensionless")
 
 """
     read_kinetic_file(path; group="/") -> KineticProfileData
@@ -139,7 +141,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", "zeff")
+        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_i", "chi_phi", "zeff")
         extras = Dict{String,Vector{Float64}}()
         for k in keys(g)
             (k in standard || !startswith(k, "n_")) && continue
@@ -149,7 +151,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"), zeff=rd("zeff"),
+            chi_e=rd("chi_e"), chi_i=rd("chi_i"), chi_phi=rd("chi_phi"), zeff=rd("zeff"),
             species_densities=(isempty(extras) ? nothing : extras), provenance=prov)
     end
 end
@@ -180,6 +182,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)
         put("zeff", data.zeff)
         # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip.