diff --git a/docs/development/references.md b/docs/development/references.md index 43f3daf74..978e6826f 100644 --- a/docs/development/references.md +++ b/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) diff --git a/docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf b/docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf new file mode 100644 index 000000000..184539270 Binary files /dev/null and b/docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf differ diff --git a/docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf b/docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf new file mode 100644 index 000000000..014a7212c Binary files /dev/null and b/docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf differ diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index ad866cf06..ac81a81a3 100644 --- a/docs/src/kinetic_forces.md +++ b/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 diff --git a/docs/src/stability.md b/docs/src/stability.md index b056753ab..4f16debf7 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -326,28 +326,47 @@ 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, and GPEC's +only option before), from `chi_i`, from the toroidal Prandtl number, from its `c_β²` floor, or +from a whole-plasma confinement time (`"tau_E"`); see `slayer_parameters`. 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 | |---|---| diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index 61ce57700..db179dde0 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -72,17 +72,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_perp_i = 1.0 # fallback only; the kinetic file supplies χ⊥,i(ψ) 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) 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] diff --git a/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 b/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 index dab58dc3f..185329e82 100644 Binary files a/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 and b/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 differ diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 448ac5d07..109feeaa5 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -81,7 +81,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) diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml index eb2180e03..79f9a4809 100644 --- a/examples/DIIID-like_truncation_example/gpec.toml +++ b/examples/DIIID-like_truncation_example/gpec.toml @@ -45,6 +45,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 diff --git a/regression-harness/cases/diiid_slayer_n1_chi_perp_e.toml b/regression-harness/cases/diiid_slayer_n1_chi_perp_e.toml new file mode 100644 index 000000000..2968a1750 --- /dev/null +++ b/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 diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 3d4132d72..a73bb5d35 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/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,6 +82,7 @@ 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}}} @@ -88,11 +90,11 @@ struct KineticProfileData 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") @@ -100,7 +102,7 @@ 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. diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index c04836766..62afff79d 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -583,9 +583,10 @@ function _overlap_slayer_inputs(inputs, dir_path::AbstractString) # fallbacks -- the same precedence run_slayer uses. Every other plasma parameter the layer # widths depend on is the run's own [SLAYER] setting, so the scan and the SLAYER stage see # one plasma. + chi = Runner._slayer_chi_inputs(loaded, ctrl; warn=false) return (; profiles=loaded.profiles, - chi_perp=loaded.chi_perp === nothing ? ctrl.chi_perp : loaded.chi_perp, - chi_tor=loaded.chi_tor === nothing ? ctrl.chi_tor : loaded.chi_tor, + chi.chi_perp_e, chi.chi_perp_i, chi.chi_tor, + tau_E=something(ctrl.tau_E, NaN), P_perp_model=ctrl.P_perp_model, zeff=loaded.zeff === nothing ? ctrl.zeff : loaded.zeff, bt=ctrl.bt, mu_i=ctrl.ion_mass, theta=ctrl.theta_sample, resistivity_model=Runner._build_resistivity_model(ctrl.resistivity_model), @@ -612,7 +613,8 @@ function _layer_overlap_scan(overlap_inputs, intr::ForceFreeStatesInternal, equi # Eq. (30), so the scan is driven in that label regardless of the SLAYER default. scans = [ Tearing.resistive_layer_overlap(equil, overlap_inputs.profiles; n_tor=n, rs_method=:flux, - chi_perp=overlap_inputs.chi_perp, chi_tor=overlap_inputs.chi_tor, + chi_perp_e=overlap_inputs.chi_perp_e, chi_perp_i=overlap_inputs.chi_perp_i, chi_tor=overlap_inputs.chi_tor, + tau_E=overlap_inputs.tau_E, P_perp_model=overlap_inputs.P_perp_model, bt=overlap_inputs.bt, mu_i=overlap_inputs.mu_i, zeff=overlap_inputs.zeff, theta=overlap_inputs.theta, resistivity_model=overlap_inputs.resistivity_model, lnLambda_form=overlap_inputs.lnLambda_form) for n in intr.nlow:intr.nhigh diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 949f02963..5f642b25c 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -185,8 +185,17 @@ profiles, without an intermediate file round-trip. - `mu_i` -- ion mass in proton-mass units (default `2.0` for D). - `zeff` -- effective charge. Scalar or a callable of `psi`, evaluated per surface (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` (no default: NaN); used by `P_perp_model = :chi_perp_i` and `:D_perp`. + - `tau_E` -- whole-plasma energy confinement time τ̄_E [s], one scalar for every surface; + required by `P_perp_model=:tau_E`, which also takes the minor radius `a` of the boundary + surface (ψ = 1) from the `:flux` label `√(2ψ_t/B₀)` whatever `rs_method` is. Every P_perp + model is μ₀χ/η, so `r_s` cancels except through η's trapped fraction; in this mode `a` sets + χ_E, and the flux label is the area-equivalent radius of a shaped boundary. + - `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²` @@ -237,7 +246,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; mu_i::Real=2.0, zeff=1.0, z_i::Real=1.0, - chi_perp=1.0, + chi_perp_e=1.0, + chi_perp_i=NaN, + tau_E::Real=NaN, + P_perp_model::Symbol=:chi_perp_e, chi_tor=1.0, dr_val=nothing, dgeo_val=nothing, @@ -261,19 +273,19 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; end _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta, q_beyond=q_beyond) + # Boundary minor radius on the flux label, needed only to turn a whole-plasma τ̄_E into a uniform χ_E. + a_minor = P_perp_model === :tau_E ? first(radial_label(equil; rs_method=:flux, theta=theta, q_beyond=q_beyond))(1.0) : NaN # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline derivatives. # In flux coordinates ω_* = n·(dp/dψ)/(e·n_e), so the surface's toroidal mode # 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) @@ -296,6 +308,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. The tearing solve 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) @@ -383,8 +396,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=_eval(zeff, psi), - 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, a_minor=a_minor, + P_perp_model=P_perp_model, m=m_res, n=n_res, dr_val=dr_val_k, dgeo_val=dgeo_val_k, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index cb8a02345..58b7d903e 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/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, :fitzpatrick, :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 @@ -134,7 +137,7 @@ end # large-aspect-ratio offsets are Connor, Ham, Hastie & Liu, PPCF 57, 065001 (2015) Eqs. (59) and (61), # after Lütjens, Luciani & Garbet, Phys. Plasmas 8, 4267 (2001). 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 || @@ -151,7 +154,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 @@ -167,13 +170,13 @@ 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 === :fitzpatrick 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 @@ -181,8 +184,11 @@ 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=NaN, tau_E=NaN, a_minor=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(), @@ -212,9 +218,36 @@ 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` feeds both - `P_perp` (which wants the one-fluid energy diffusivity χ_E) and the critical-Δ `W_d` (which wants - the electron channel); one profile serves both for now, see the stability docs + - `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` and `:D_perp` + - `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 τ̄_E [s], required by `P_perp_model=:tau_E` + - `a_minor` -- plasma minor radius `a` [m], required by `P_perp_model=:tau_E` (`build_slayer_inputs` takes it from the flux label) + - `P_perp_model` -- how `P_perp` is built: + + `:chi_perp_e` (default of this function, which has no profile gradients) -- τ_R·χ⊥,e/r_s². + The electron channel is a proxy for the one-fluid energy diffusivity; the `[SLAYER]` deck + default is `:D_perp` + + `: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 the perpendicular particle diffusivity D_⊥ of + Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eq. 16, with η_⊥ = η_∥: + D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)[η_e τ_e/(1+η_e) χ⊥,e + η_i τ_i/(1+η_i) χ⊥,i], + where η_s = dlnT_s/dlnn, τ_e = τ/(1+τ), τ_i = 1/(1+τ) and τ = (T_e/T_i)(1+η_e)/(1+η_i) + (Eqs. 4-6; the inverse orientation of the `tau` = T_i/T_e field of `SLAYERParameters`). + The weights equal n T_s′/p′, so the bracket is (n T_e′ χ⊥,e + n T_i′ χ⊥,i)/p′. Where both + temperature profiles are flat the χ terms vanish and P_perp falls to the `:c_beta` value; + where the density gradient opposes the temperature gradients (hollow density) the weights exceed 1 + + `:c_beta` -- P_perp = C² with C² = c_β² + (1 − c_β²)K and K = 0 (J.-K. Park, Phys. Plasmas 29, + 072506 (2022), Eq. 9, C² defined below Eq. 10). This is the χ⊥ → 0 limit of `:D_perp`, so + it is a floor on P_perp; realistic anomalous transport makes D_⊥ much larger + + `:tau_E` -- P_perp = P_tor = τ_R/τ_E(r_s) (Fitzpatrick 2023, Eq. 131), with the local + τ_E(r) = r²/χ_E of a uniform χ_E (Eqs. 81-83). χ_E comes from the whole-plasma τ̄_E through + τ̄_E = κ_e a²/χ_E (Eqs. 37-38) with κ_e = 3/16, the value for T_e ∝ 1 − r̂². κ_e is fixed + rather than computed because Eq. 38 needs dT_e/dr at the edge, which a pedestal makes + ill-conditioned; it carries a shape uncertainty of order 2-3, and in H-mode the χ_E implied + by τ̄_E falls below the core power-balance χ. `chi_tor` is ignored in this mode + `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`, `:fitzpatrick`, `:toroidal` @@ -267,8 +300,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=NaN, tau_E::Real=NaN, a_minor::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(), @@ -370,10 +406,39 @@ 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 && isfinite(a_minor) && a_minor > 0 || + throw(ArgumentError("slayer_parameters: P_perp_model=:tau_E needs positive tau_E and a_minor, got $tau_E, $a_minor")) + # P_⊥ = P_φ = τ_R/τ_E(r_s) with τ_E(r_s) = r_s²/χ_E and χ_E = κ_e a²/τ̄_E, κ_e = 3/16 (Fitzpatrick 2023, Eqs. 37-38, 81-83, 131) + chi_E = (3 / 16) * a_minor^2 / tau_E + P_tor = tau_r / (rs^2 / chi_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 (η_⊥ = η_∥); the weights η_s τ_s/(1+η_s) equal n T_s'/p'. + dp = t_e * (dlnn + dlnTe) + t_i * (dlnn + dlnTi) + D_perp = c_beta^2 * eta / MU_0 + + (2 / 3) * (1 - c_beta^2) * (t_e * dlnTe * chi_perp_e + t_i * dlnTi * chi_perp_i) / dp + # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥ + tau_r / (rs^2 / D_perp) + elseif P_perp_model === :c_beta + # P_⊥ = C² = c_β² + (1 − c_β²)K with K = 0 (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 9) + c_beta^2 + else + throw(ArgumentError("slayer_parameters: unknown P_perp_model=$P_perp_model. " * + "Allowed: $(ALLOWED_P_PERP_MODELS)")) + end + isfinite(P_perp) && P_perp > 0 || + throw(ArgumentError("slayer_parameters: P_perp_model=:$P_perp_model gave P_perp=$P_perp; check its inputs " * + "(chi_perp_e=$chi_perp_e, chi_perp_i=$chi_perp_i, dlnn=$dlnn, dlnTe=$dlnTe, dlnTi=$dlnTi)")) # Normalized beta-related width and Δ-normalization d_beta = c_beta * d_i @@ -382,7 +447,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) diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl index a07d4ea48..9e55b178d 100644 --- a/src/InnerLayer/SLAYER/LayerThickness.jl +++ b/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 diff --git a/src/InnerLayer/SLAYER/Riccati.jl b/src/InnerLayer/SLAYER/Riccati.jl index 87781259b..9f630838b 100644 --- a/src/InnerLayer/SLAYER/Riccati.jl +++ b/src/InnerLayer/SLAYER/Riccati.jl @@ -2,7 +2,7 @@ # # Inner-layer Δ via the Fitzpatrick two-fluid drift-MHD layer Riccati ODE # (Fitzpatrick, Tearing Mode Dynamics in Tokamak Plasmas, IOP 2023; Park -# et al. 2022, Phys. Plasmas 29, 122505; Burgess et al. 2026), for the +# 2022, Phys. Plasmas 29, 072506; Burgess et al. 2026), for the # pressureless tearing channel (no parallel flow, pe = 0). The A, B, C # coefficients, the dW/dp Riccati equation, the large-p asymptotic boundary # conditions, the regime test D² > ι_e P_⊥/P_tor^(2/3), the small-p @@ -169,7 +169,7 @@ pressureless layer produces only the tearing channel. # Algorithm Implements the Fitzpatrick two-fluid drift-MHD layer formulation -(Fitzpatrick 2023; Park et al. 2022) for the pressureless tearing channel +(Fitzpatrick 2023; Park 2022) for the pressureless tearing channel (no parallel flow, pe = 0). Integrates `dW/dp = -(fA'/p)·W − W²/p + (fB/(fA·fC))·p³` from a large `p_start` (selected by `_riccati_f_initial` according to whether diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index 7f6c95c00..efc81cd31 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -1,8 +1,8 @@ # SLAYER.jl # # SLAYER (Slab Layer) two-fluid drift-MHD inner-layer model (Fitzpatrick, -# Tearing Mode Dynamics in Tokamak Plasmas, IOP 2023; Park et al. 2022, -# Phys. Plasmas 29, 122505; Burgess et al. 2026). The dispersion path uses +# Tearing Mode Dynamics in Tokamak Plasmas, IOP 2023; Park 2022, +# Phys. Plasmas 29, 072506; Burgess et al. 2026). The dispersion path uses # the Fitzpatrick layer formulation — P_perp / P_tor transport, c_beta # compressibility, D_norm normalized ion-skin scale, two-fluid drift # coupling via Q_e, Q_i, iota_e (see `Riccati.jl`). diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl index 64c187d2d..d9fd1534f 100644 --- a/src/Tearing/LayerOverlap.jl +++ b/src/Tearing/LayerOverlap.jl @@ -155,7 +155,7 @@ Riccati solve; `delta_s` is solved per surface as a reported comparison only. - `rs_method` -- radial label the widths and their ψ-conversion Jacobian are taken in (default `:midplane`; see `radial_label`). The caller wiring this scan into the domain truncation passes `:flux`, the label Fitzpatrick Eq. (100) is anchored to - - `bt`, `mu_i`, `zeff`, `chi_perp`, `chi_tor`, `resistivity_model`, `lnLambda_form` + - `bt`, `mu_i`, `zeff`, `chi_perp_e`, `chi_perp_i`, `chi_tor`, `tau_E`, `P_perp_model`, `resistivity_model`, `lnLambda_form` -- passed through to `build_slayer_inputs`. The defaults are placeholders: a caller scanning a real run must pass that run's own values, as the domain-truncation wiring does @@ -176,8 +176,11 @@ function resistive_layer_overlap(equil, profiles::KineticProfiles; bt=nothing, mu_i::Real=2.0, zeff=1.0, - chi_perp=1.0, + chi_perp_e=1.0, + chi_perp_i=NaN, chi_tor=1.0, + tau_E::Real=NaN, + P_perp_model::Symbol=:chi_perp_e, resistivity_model::NeoResistivityModel=SauterNeoModel(), lnLambda_form::Symbol=:nrl) @@ -316,7 +319,8 @@ function resistive_layer_overlap(equil, profiles::KineticProfiles; # which makes iota_e singular) is recorded and skipped rather than aborting the scan. params = try build_slayer_inputs(equil, [sing], profiles; rs_method=rs_method, q_beyond=q_beyond, - bt=bt, mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, + bt=bt, mu_i=mu_i, zeff=zeff, chi_perp_e=chi_perp_e, chi_perp_i=chi_perp_i, chi_tor=chi_tor, + tau_E=tau_E, P_perp_model=P_perp_model, dr_val=0.0, dc_type=:none, theta=theta, resistivity_model=resistivity_model, lnLambda_form=lnLambda_form) catch err diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index a82a2ce18..9434873c1 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -40,10 +40,19 @@ constructor. - `impurity_charge` -- charge of the single impurity species that closes quasi-neutrality when `Z_eff(ψ)` is derived from `n_i/n_e` (default 6.0, carbon); the main ion has charge 1. Unused when the file supplies a `zeff` profile - - `chi_perp`, `chi_tor` -- fallback perpendicular / toroidal heat - diffusivity [m²/s], used only when the kinetic file carries no usable - `chi_e`/`chi_phi` profile (dataset absent or all-zero); otherwise the - file's χ⊥(ψ)/χ_φ(ψ) take precedence + - `chi_perp_e`, `chi_perp_i`, `chi_tor` -- fallback electron / ion perpendicular heat and + toroidal momentum diffusivity [m²/s], used only when the kinetic file carries no usable + `chi_e`/`chi_i`/`chi_phi` profile (dataset absent or all-zero); otherwise the file's + χ⊥,e(ψ)/χ⊥,i(ψ)/χ_φ(ψ) take precedence. `chi_perp_i` has no default (NaN), so a model that + needs χ⊥,i fails loudly when neither the file nor the deck supplies it + - `tau_E` -- whole-plasma energy confinement time τ̄_E [s], required by + `P_perp_model = :tau_E` (default `nothing`) + - `P_perp_model` -- how the perpendicular Prandtl number is built: `:D_perp` (default; D⊥ from + χ⊥,e, χ⊥,i and the n, T gradients, Fitzpatrick 2022 Eq. 16, so χ⊥,i is required), `:chi_perp_e` + (the electron-channel proxy used before this option existed), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:c_beta` + (P_perp = c_β²), or `:tau_E` (a uniform χ_E = (3/16)a²/τ̄_E sets both P_perp and P_tor, + so `chi_phi`/`chi_tor` are ignored); see `InnerLayer.SLAYER.slayer_parameters`. + `chi_perp_e` still sets the critical-Δ width in every mode - `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 @@ -90,9 +99,9 @@ there is one consistent interface for resistive and kinetic profiles. - `profile_file` -- path to a kinetic-profile file (relative to the run dir), required when SLAYER is enabled. HDF5 (`.h5`) files use the GPEC - kinetic schema and may carry `chi_e` (χ⊥) and `chi_phi` (χ_φ); ASCII + kinetic schema and may carry `chi_e` (χ⊥,e), `chi_i` (χ⊥,i) and `chi_phi` (χ_φ); ASCII tables are also accepted but carry no χ. When a χ dataset is absent or - all-zero, the scalar `chi_perp`/`chi_tor` fallbacks below are used (so a + all-zero, the scalar `chi_perp_e`/`chi_perp_i`/`chi_tor` fallbacks below are used (so a file can keep the χ keys set to 0 to defer to the scalars). - `profile_group` -- group within the HDF5 file (default `"/"`) @@ -114,8 +123,11 @@ there is one consistent interface for resistive and kinetic profiles. ion_mass::Float64 = 2.0 zeff::Float64 = 1.0 impurity_charge::Float64 = 6.0 - chi_perp::Float64 = 1.0 + chi_perp_e::Float64 = 1.0 + chi_perp_i::Float64 = NaN chi_tor::Float64 = 1.0 + tau_E::Union{Float64,Nothing} = nothing + P_perp_model::Symbol = :D_perp D_R::Union{Float64,Nothing} = nothing D_geo::Union{Float64,Nothing} = nothing theta_sample::Float64 = 0.0 @@ -168,6 +180,7 @@ const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_ray, :ggj_shooting, :ggj_ const _VALID_SCAN_MODES = (:amr, :brute_force) const _VALID_COUPLING_MODES = (:uncoupled, :coupled) const _VALID_D_C_TYPES = (:none, :lar, :fitzpatrick, :toroidal) +const _VALID_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta, :tau_E) const _VALID_RESISTIVITY_MODELS = (:sauter, :redl, :spitzer, :spitzer_harm) const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson) @@ -184,6 +197,11 @@ function validate(ctrl::SLAYERControl) ctrl.D_c_type in _VALID_D_C_TYPES || throw(ArgumentError("SLAYERControl: D_c_type=$(ctrl.D_c_type) " * "not in $(_VALID_D_C_TYPES)")) + 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)")) + ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && isfinite(ctrl.tau_E) && ctrl.tau_E > 0) || + throw(ArgumentError("SLAYERControl: P_perp_model=tau_E requires a positive, finite tau_E [s]")) ctrl.resistivity_model in _VALID_RESISTIVITY_MODELS || throw(ArgumentError("SLAYERControl: resistivity_model=$(ctrl.resistivity_model) " * "not in $(_VALID_RESISTIVITY_MODELS)")) @@ -241,9 +259,10 @@ function slayer_control_from_toml(section::AbstractDict) # Validate keys against the struct fields field_names = Set(String.(fieldnames(SLAYERControl))) unknown = [k for k in keys(flat) if !(k in field_names)] + hint = "chi_perp" in unknown ? " (`chi_perp` is now `chi_perp_e`.)" : "" isempty(unknown) || throw(ArgumentError("slayer_control_from_toml: unknown keys " * - "$(unknown) in [SLAYER] section. Known: " * + "$(unknown) in [SLAYER] section.$hint Known: " * "$(sort(collect(field_names))).")) # Coerce types where needed @@ -251,11 +270,11 @@ function slayer_control_from_toml(section::AbstractDict) for (k, v) in flat sym = Symbol(k) if sym in (:inner_model, :scan_mode, :coupling_mode, :D_c_type, - :resistivity_model, :lnLambda_form) + :resistivity_model, :lnLambda_form, :P_perp_model) 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) + elseif sym in (:bt, :D_R, :D_geo, :tau_E) # Allow explicit nothing (auto-derive) or a number (override) kwargs[sym] = v === nothing ? nothing : Float64(v) elseif sym === :boxes diff --git a/src/Tearing/Runner/TearingProblem.jl b/src/Tearing/Runner/TearingProblem.jl index 0ba447eb4..2b203f279 100644 --- a/src/Tearing/Runner/TearingProblem.jl +++ b/src/Tearing/Runner/TearingProblem.jl @@ -32,12 +32,8 @@ function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control) # physical B_T = F(psi)/(2*pi*R_0) per surface from the equilibrium's F-spline, which is # what its docstring already prescribes. bt = control.bt - # χ⊥/χ_φ from the kinetic file when present, else the scalar fallbacks. - chi_perp = loaded.chi_perp === nothing ? control.chi_perp : loaded.chi_perp - chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor - (loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn( - "SLAYER: no usable chi_e/chi_phi profile(s) (dataset absent or all-zero); " * - "using the scalar chi_perp/chi_tor fallback for the missing one(s).") + # χ⊥,e/χ⊥,i/χ_φ from the kinetic file when present, else the scalar fallbacks. + chi = _slayer_chi_inputs(loaded, control) # Z_eff(ψ) from the kinetic file when it carries n_i, else the scalar fallback. zeff = loaded.zeff === nothing ? control.zeff : loaded.zeff if loaded.zeff_source === :scalar @@ -51,8 +47,11 @@ function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control) bt=bt, mu_i=control.ion_mass, zeff=zeff, - chi_perp=chi_perp, - chi_tor=chi_tor, + chi_perp_e=chi.chi_perp_e, + chi_perp_i=chi.chi_perp_i, + chi_tor=chi.chi_tor, + tau_E=something(control.tau_E, NaN), + P_perp_model=control.P_perp_model, dr_val=control.D_R, dgeo_val=control.D_geo, dc_type=control.D_c_type, diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index f50f4101e..dfbd647b4 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -8,9 +8,9 @@ # --------------------------------------------------------------------- # Read kinetic profiles through the shared `Equilibrium.read_kinetic_file` # reader and adapt them to the SLAYER inner layer. Returns the spline-based -# `KineticProfiles` plus χ⊥(ψ)/χ_φ(ψ) callables built from the file's -# `chi_e`/`chi_phi` (or `nothing` when the file carries no usable χ, in which -# case the caller falls back to the scalar `control.chi_perp`/`chi_tor`), and +# `KineticProfiles` plus χ⊥,e(ψ)/χ⊥,i(ψ)/χ_φ(ψ) callables built from the file's +# `chi_e`/`chi_i`/`chi_phi` (or `nothing` when the file carries no usable χ, in which +# case the caller falls back to the scalar `control.chi_perp_e`/`chi_perp_i`/`chi_tor`), and # a Z_eff(ψ) callable: the file's own `zeff` profile when it carries one, else # derived from its `n_i/n_e` (or `nothing` when the file carries neither, in which # case the caller falls back to the scalar `control.zeff`). @@ -39,12 +39,13 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString) # χ⊥(ψ)/χ_φ(ψ) splines from the file. A χ array that is absent OR all-zero # is treated as "not provided" — χ must be positive (χ=0 ⇒ τ_⊥→∞), and the # all-zero sentinel lets a file keep the chi_e/chi_phi keys while deferring - # to the scalar control.chi_perp/chi_tor fallback. Build the spline only + # to the scalar control.chi_perp_e/chi_tor fallback. Build the spline only # from a usable (present, not all-zero) array. psi_xs = collect(Float64, data.psi) _chi_spline(v) = (v === nothing || all(iszero, v)) ? nothing : cubic_interp(psi_xs, collect(Float64, v)) - chi_perp = _chi_spline(data.chi_e) + chi_perp_e = _chi_spline(data.chi_e) + chi_perp_i = _chi_spline(data.chi_i) chi_tor = _chi_spline(data.chi_phi) # Z_eff(ψ): a profile supplied in the file wins, since a real impurity mix is not recoverable # from n_i/n_e. Otherwise quasi-neutrality with one impurity species of charge `impurity_charge` and a @@ -64,7 +65,34 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString) cubic_interp(psi_xs, z) end zeff_source = _usable(data.zeff) ? :file : zeff === nothing ? :scalar : :closure - return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source) + return (profiles=profiles, chi_perp_e=chi_perp_e, chi_perp_i=chi_perp_i, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source) +end + +# χ⊥,e/χ⊥,i/χ_φ for the layer: the kinetic file's profiles where present, else the control's +# scalar fallbacks. With `warn`, a fallback is reported only for a χ the selected P_perp model +# (and, for χ⊥,e, the critical Δ) actually uses. +function _slayer_chi_inputs(loaded, control::SLAYERControl; warn::Bool=true) + model = control.P_perp_model + pick(f, s) = f === nothing ? s : f + model in (:chi_perp_i, :D_perp) && loaded.chi_perp_i === nothing && !(isfinite(control.chi_perp_i) && control.chi_perp_i > 0) && + throw(ArgumentError("SLAYER: P_perp_model = $model needs the ion diffusivity χ⊥,i: add a chi_i profile to the " * + "kinetic file, set [SLAYER] chi_perp_i, or choose P_perp_model = \"chi_perp_e\".")) + if warn + uses_e = model in (:chi_perp_e, :D_perp) || control.D_c_type !== :none + uses_i = model in (:chi_perp_i, :D_perp) + uses_tor = model !== :tau_E + fallbacks = String[] + uses_e && loaded.chi_perp_e === nothing && push!(fallbacks, "chi_e → chi_perp_e = $(control.chi_perp_e)") + uses_i && loaded.chi_perp_i === nothing && push!(fallbacks, "chi_i → chi_perp_i = $(control.chi_perp_i)") + 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." + 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_tor=pick(loaded.chi_tor, control.chi_tor)) end # --------------------------------------------------------------------- diff --git a/test/runtests_dispersion_coupled.jl b/test/runtests_dispersion_coupled.jl index 1b950016e..30efddbd4 100644 --- a/test/runtests_dispersion_coupled.jl +++ b/test/runtests_dispersion_coupled.jl @@ -24,7 +24,7 @@ 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) + chi_perp_e=1.0, chi_tor=1.0, m=2, n=1) end @testset "Constructor validation" begin diff --git a/test/runtests_dispersion_residual.jl b/test/runtests_dispersion_residual.jl index 30ebf24ff..dc5bbc420 100644 --- a/test/runtests_dispersion_residual.jl +++ b/test/runtests_dispersion_residual.jl @@ -24,7 +24,7 @@ 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) + chi_perp_e=1.0, chi_tor=1.0, m=2, n=1) end @testset "Constructor scale defaults" begin diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index 4d141f626..6410bfba3 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -679,8 +679,9 @@ da = E.read_kinetic_file(gpeckf) db = E.read_kinetic_file(h5) # legacy ASCII carries no optional fields; the example .h5 carries chi - @test da.chi_e === nothing && da.chi_phi === nothing - @test db.chi_e !== nothing && db.chi_phi !== nothing + @test da.chi_e === nothing && da.chi_i === nothing && da.chi_phi === nothing + @test db.chi_e !== nothing && db.chi_i !== nothing && db.chi_phi !== nothing + @test all(>(0), db.chi_i) # core columns identical (the .h5 was generated from the .gpeckf) @test da.psi == db.psi @test da.n_e == db.n_e @@ -708,6 +709,9 @@ @test dc.chi_phi ≈ db.chi_phi @test dc.n_e ≈ db.n_e @test dc.provenance == "roundtrip" + @test dc.chi_i ≈ db.chi_i + E.write_kinetic_h5(tmp, E.KineticProfileData(; psi=db.psi, chi_e=db.chi_e)) + @test E.read_kinetic_file(tmp).chi_i === nothing rm(tmp; force=true) end diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl index 0cedfb9e9..a81dfdd70 100644 --- a/test/runtests_layer_overlap_fitzpatrick.jl +++ b/test/runtests_layer_overlap_fitzpatrick.jl @@ -74,7 +74,7 @@ omega=0.0, omega_e=1.0e4, omega_i=5.0e3, qval=q_of(r), sval_r=s_of(r), bt=B0, rs=r * A_MIN, R0=R0, mu_i=MNUM, zeff=ZEFF, - chi_perp=CHI, chi_tor=CHI, m=1, n=n, + chi_perp_e=CHI, chi_tor=CHI, m=1, n=n, resistivity_model=SpitzerHarmModel(), lnLambda_form=:nrl) return p end diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index ccb881ecb..bf3ef08d2 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -2,6 +2,7 @@ using GeneralizedPerturbedEquilibrium using GeneralizedPerturbedEquilibrium.Equilibrium using GeneralizedPerturbedEquilibrium.Utilities + using GeneralizedPerturbedEquilibrium.Utilities: MU_0 using GeneralizedPerturbedEquilibrium.InnerLayer using GeneralizedPerturbedEquilibrium.ForceFreeStates: SingType using TOML @@ -142,30 +143,48 @@ @test sl_double[1].lu != sl_default[1].lu end - @testset "build_slayer_inputs: chi_perp/chi_tor as scalars and callables" begin + @testset "build_slayer_inputs: chi_perp_e/chi_tor as scalars and callables" begin sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] # Scalar (dr_val=0.0 bypasses the sing.restype requirement; see comment above) sl_s = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=2.0, chi_tor=1.5, dr_val=0.0) + bt=2.0, chi_perp_e=2.0, chi_tor=1.5, dr_val=0.0) # Callable with matching value chi_p(psi) = 2.0 + 0.0*psi chi_t(psi) = 1.5 + 0.0*psi sl_c = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=chi_p, chi_tor=chi_t, dr_val=0.0) + bt=2.0, chi_perp_e=chi_p, chi_tor=chi_t, dr_val=0.0) @test sl_s[1].P_perp ≈ sl_c[1].P_perp @test sl_s[1].P_tor ≈ sl_c[1].P_tor # Callable with ψ-dependence changes the result chi_p_var(psi) = 1.0 + 10.0 * psi # χ⊥(0.5) = 6.0 > 2.0 sl_var = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=chi_p_var, chi_tor=1.5, dr_val=0.0) + bt=2.0, chi_perp_e=chi_p_var, chi_tor=1.5, dr_val=0.0) # P_perp = τ_r · χ⊥ / r² grows with χ⊥, so the varying-χ case at # ψ=0.5 (χ⊥=6) gives a *larger* P_perp than the scalar χ⊥=2. @test sl_var[1].P_perp > sl_s[1].P_perp @test sl_var[1].P_perp ≈ sl_s[1].P_perp * 6.0 / 2.0 rtol = 1e-10 end + @testset "build_slayer_inputs: P_perp_model forwards the profile gradients and minor radius" begin + sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] + kw = (; bt=2.0, chi_perp_e=2.0, chi_perp_i=1.0, chi_tor=1.5, dr_val=0.0) + ref = build_slayer_inputs(equil, sings, profiles; kw...)[1] + # Flat n_e and linear T_e, T_i: the D_⊥ weights are T_s′/(T_e′ + T_i′) = 7/13 and 6/13 + dp = build_slayer_inputs(equil, sings, profiles; kw..., P_perp_model=:D_perp)[1] + D_hand = dp.c_beta^2 * dp.eta / MU_0 + (2 / 3) * (1 - dp.c_beta^2) * (2.0 * 7 / 13 + 1.0 * 6 / 13) + @test dp.P_perp ≈ ref.P_perp * D_hand / 2.0 rtol = 1e-8 + # τ̄_E mode takes a from the boundary surface on the flux label, whatever rs_method is + te = build_slayer_inputs(equil, sings, profiles; kw..., P_perp_model=:tau_E, tau_E=0.1)[1] + a = first(radial_label(equil; rs_method=:flux))(1.0) + ref_fsa = build_slayer_inputs(equil, sings, profiles; kw..., rs_method=:fsa)[1] + te_fsa = build_slayer_inputs(equil, sings, profiles; kw..., P_perp_model=:tau_E, tau_E=0.1, rs_method=:fsa)[1] + @test te_fsa.P_perp ≈ ref_fsa.P_perp * ((3 / 16) * a^2 / 0.1) / 2.0 rtol = 1e-12 + @test te.P_perp ≈ ref.P_perp * ((3 / 16) * a^2 / 0.1) / 2.0 rtol = 1e-12 + @test te.P_tor == te.P_perp + end + @testset "build_slayer_inputs: rs_method radial labels are self-consistent" begin sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] got = Dict{Symbol,Any}() diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index ac7c15c8a..f69a33b72 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -12,7 +12,7 @@ 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, + chi_perp_e=1.0, chi_tor=1.0, m=2, n=1, dr_val=dr_val, dgeo_val=0.5, dc_type=dc_type, ising=3 @@ -79,7 +79,7 @@ # Compressibility is in (0,1) for finite β @test 0.0 < p.c_beta < 1.0 - # Prandtl-like ratios are positive and equal here (chi_perp=chi_tor=1) + # Prandtl-like ratios are positive and equal here (chi_perp_e=chi_tor=1) @test p.P_perp ≈ p.P_tor @test p.P_perp > 0 @@ -199,4 +199,43 @@ @test neg.lu > 0 @test neg.tauk > 0 end + + @testset "Test 4: P_perp_model options" begin + base = merge(_ref_kwargs(), (chi_tor=1.7,)) + ref = slayer_parameters(; base...) + @test ref.P_perp ≈ ref.tau_r * 1.0 / ref.rs^2 rtol = 1e-14 + @test slayer_parameters(; base..., P_perp_model=:chi_perp_i, chi_perp_i=2.0).P_perp ≈ 2 * ref.P_perp + @test slayer_parameters(; base..., P_perp_model=:P_phi).P_perp == ref.P_tor + # Eq. 16: with χ⊥,e = χ⊥,i = 0 only c_β²η/μ₀ survives, so P_perp = c_β² + grads = (dlnn=-1.0, dlnTe=-2.0, dlnTi=-3.0) + d0 = slayer_parameters(; base..., P_perp_model=:D_perp, grads..., chi_perp_e=0.0, chi_perp_i=0.0) + @test d0.P_perp ≈ d0.c_beta^2 rtol = 1e-12 + # η_e = 2, η_i = 3, T_e = T_i ⇒ τ = 3/4, τ_e = 3/7, τ_i = 4/7 + d1 = slayer_parameters(; base..., P_perp_model=:D_perp, grads..., chi_perp_i=2.0) + D_hand = d1.c_beta^2 * d1.eta / MU_0 + (2 / 3) * (1 - d1.c_beta^2) * ((2 / 3) * (3 / 7) * 1.0 + (3 / 4) * (4 / 7) * 2.0) + @test d1.P_perp ≈ d1.tau_r * D_hand / d1.rs^2 rtol = 1e-12 + # T_e = 2T_i ⇒ τ = 3/2, τ_e = 3/5, τ_i = 2/5: fixes the orientation of the temperature ratio + d2 = slayer_parameters(; base..., t_e=2000.0, P_perp_model=:D_perp, grads..., chi_perp_i=2.0) + D_hand2 = d2.c_beta^2 * d2.eta / MU_0 + (2 / 3) * (1 - d2.c_beta^2) * ((2 / 3) * (3 / 5) * 1.0 + (3 / 4) * (2 / 5) * 2.0) + @test d2.P_perp ≈ d2.tau_r * D_hand2 / d2.rs^2 rtol = 1e-12 + # Flat temperatures leave only the c_β² floor, and a flat ion profile is finite (no 0·∞) + @test slayer_parameters(; base..., P_perp_model=:D_perp, dlnn=-1.0, dlnTe=0.0, dlnTi=0.0, chi_perp_i=2.0).P_perp ≈ d0.P_perp rtol = 1e-12 + @test isfinite(slayer_parameters(; base..., P_perp_model=:D_perp, dlnn=0.0, dlnTe=-2.0, dlnTi=0.0, chi_perp_i=2.0).P_perp) + # A non-positive or missing χ⊥,i is an error in every model that uses it + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:chi_perp_i) + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:chi_perp_i, chi_perp_i=0.0) + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:D_perp) + @test slayer_parameters(; base..., P_perp_model=:c_beta).P_perp == ref.c_beta^2 + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:bogus) + # τ̄_E = 0.1 s and a = 0.6 m ⇒ uniform χ_E = (3/16)·0.36/0.1 m²/s sets both Prandtl numbers + te = slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1, a_minor=0.6) + @test te.P_tor ≈ te.tau_r * ((3 / 16) * 0.6^2 / 0.1) / te.rs^2 rtol = 1e-14 + @test te.P_perp == te.P_tor + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:tau_E) + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1) + # Only P_perp moves; the critical-Δ offset keeps using chi_perp_e. + cb = slayer_parameters(; merge(base, (dr_val=-1.0, dc_type=:lar))..., P_perp_model=:c_beta) + @test cb.dc_tmp == slayer_parameters(; merge(base, (dr_val=-1.0, dc_type=:lar))...).dc_tmp + @test cb.P_tor == ref.P_tor + end end diff --git a/test/runtests_slayer_riccati.jl b/test/runtests_slayer_riccati.jl index 9d4a95462..0755d51fd 100644 --- a/test/runtests_slayer_riccati.jl +++ b/test/runtests_slayer_riccati.jl @@ -18,7 +18,7 @@ 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, + chi_perp_e=1.0, chi_tor=1.0, m=2, n=1) end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 530b2d0ea..a79ca5c89 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -58,6 +58,7 @@ "bt" => 1.8, "ion_mass" => 2.0, "D_R" => 0.01, + "P_perp_model" => "D_perp", "scan_grid" => Dict{String,Any}( "Q_re_range" => [-5.0, 5.0], "Q_im_range" => [-1.0, 3.0], @@ -80,6 +81,7 @@ @test c.msing_max == 2 @test c.bt === 1.8 @test c.D_R == 0.01 + @test c.P_perp_model === :D_perp @test c.Q_re_range == (-5.0, 5.0) @test c.Q_im_range == (-1.0, 3.0) @test c.nre == 50 @@ -93,6 +95,19 @@ # Unknown keys should raise bad = merge(section, Dict{String,Any}("mistyped_key" => 42)) @test_throws ArgumentError slayer_control_from_toml(bad) + # The renamed χ⊥,e key and bad P_perp_model values are rejected + @test_throws "chi_perp` is now `chi_perp_e" slayer_control_from_toml(merge(section, Dict{String,Any}("chi_perp" => 1.0))) + @test_throws ArgumentError Runner.validate(slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "bogus")))) + @test_throws ArgumentError Runner.validate(slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "tau_E")))) + @test_throws ArgumentError Runner.validate(slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "tau_E", "tau_E" => Inf)))) + c_te = slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "tau_E", "tau_E" => 0.05))) + @test c_te.tau_E == 0.05 && Runner.validate(c_te) === c_te + # D_perp is the default, and it refuses to run without an ion diffusivity from the file or the deck. + @test SLAYERControl().P_perp_model === :D_perp + no_chi = (; chi_perp_e=nothing, chi_perp_i=nothing, chi_tor=nothing) + @test_throws "needs the ion diffusivity" Runner._slayer_chi_inputs(no_chi, SLAYERControl(); warn=false) + @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; chi_perp_i=0.7); warn=false).chi_perp_i == 0.7 + @test isnan(Runner._slayer_chi_inputs(no_chi, SLAYERControl(; P_perp_model=:chi_perp_e); warn=false).chi_perp_i) end @testset "TearingProblem solve forwards surfaces and Δ'" begin