diff --git a/docs/src/stability.md b/docs/src/stability.md index 814a5207ae..0b07e6e775 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -349,7 +349,11 @@ total pressure; under the closure `T_s = T_s(n)` it is set by the perpendicular `D_⊥` of Fitzpatrick, Phys. Plasmas **29**, 032507 (2022), Eq. 16, which weights the electron and ion energy diffusivities by `n T_s′/p′`. The critical-Δ island width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature flattening against an electron parallel conductivity, so its χ_⊥ is the electron channel `chi_e` -in every case. `[SLAYER] P_perp_model` selects how `P_perp` is built: from `D_⊥` (`"D_perp"`, +in every case. `[SLAYER] chi_perp_Dc` (optional) sets that electron diffusivity for the critical Δ +alone, leaving `P_perp` on its own χ inputs. It is for a run whose `P_perp` uses a diffusivity +the critical Δ should not, for example a one-fluid χ_eff in both channels of `D_⊥` when an +interpretive power balance splits the conducted heat between species unreliably, while the +critical Δ takes the electron χ that balance allows. `[SLAYER] P_perp_model` selects how `P_perp` is built: from `D_⊥` (`"D_perp"`, the default), from `chi_e` alone (the single-χ convention of Fitzpatrick's TJ code), from `chi_i`, from the toroidal Prandtl number, from its `c_β²` floor, or from a whole-plasma confinement time (`"tau_E"`); see `slayer_parameters`. `"tau_E"` is for a plasma with no χ diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 0e72705b7b..bcb98b48a0 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -255,6 +255,8 @@ profiles, without an intermediate file round-trip. callable of `psi` (default `1.0`). - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], scalar or callable of `psi`; required by `P_perp_model = :chi_perp_i` and `:D_perp` (default `nothing`). + - `chi_perp_Dc` -- electron perpendicular heat diffusivity [m²/s] used only by the critical Δ, + scalar or callable of `psi`; `nothing` (default) uses `chi_perp_e`. See [`slayer_parameters`](@ref). - `tau_E` -- whole-plasma energy confinement time τ̄_E [s], required by `P_perp_model = :tau_E` (default `nothing`). That mode takes the minor radius `a` of the ψ = 1 surface from the `:flux` label, whatever `rs_method` is. @@ -323,6 +325,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; z_i::Real=1.0, chi_perp_e=1.0, chi_perp_i=nothing, + chi_perp_Dc=nothing, tau_E::Union{Real,Nothing}=nothing, P_perp_model::Symbol=:D_perp, chi_tor=1.0, @@ -484,6 +487,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; chi_perp_e=_eval(chi_perp_e, psi), chi_tor=_eval(chi_tor, psi), chi_perp_i=_eval(chi_perp_i, psi), + chi_perp_Dc=_eval(chi_perp_Dc, psi), dlnn=kin.dn_e / kin.n_e, dlnTe=kin.dT_e / kin.T_e, dlnTi=kin.dT_i / kin.T_i, tau_E=tau_E, a_minor=a_minor, P_perp_model=P_perp_model, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index df53b14af1..bcf812ead5 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -202,7 +202,7 @@ end m, n, chi_perp_i=nothing, tau_E=nothing, a_minor=nothing, dlnn=nothing, dlnTe=nothing, dlnTi=nothing, - P_perp_model=:D_perp, + P_perp_model=:D_perp, chi_perp_Dc=nothing, dr_val=0.0, dgeo_val=0.0, kpar_val=nothing, dc_type=:none, ising=0, resistivity_model=SauterNeoModel(), @@ -253,7 +253,13 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl a uniform χ_E = κ_e a²/τ̄_E with κ_e = 3/16, the value for T_e ∝ 1 − r̂² (Eq. 38). Other profile shapes change κ_e by a factor of order 2-3. Both Prandtl numbers come from χ_E, so `chi_tor` is ignored in this mode - `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode. + `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode, unless `chi_perp_Dc` is given. + - `chi_perp_Dc` -- electron perpendicular heat diffusivity [m²/s] used **only** by the critical Δ + (the W_d of the χ_∥ matching and the `(χ_∥/χ⊥)^(1/4)` factors of `:lar`/`:toroidal`). Default + `nothing` uses `chi_perp_e`, so the critical Δ and P_perp see one χ⊥,e as before. Set it when + the layer's P_perp is built from a different diffusivity than the electron channel the + critical Δ needs -- e.g. a one-fluid χ_eff in `:D_perp` when an interpretive power balance + splits the conducted heat between species unreliably. Must be positive and finite when given - `m`, `n` -- poloidal / toroidal mode numbers at the surface - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula: the resistive interchange index `D_R` and the Connor et al. 2015 Eq. 59 geometric factor @@ -314,7 +320,7 @@ function slayer_parameters(; m::Integer, n::Integer, chi_perp_i::Union{Real,Nothing}=nothing, tau_E::Union{Real,Nothing}=nothing, a_minor::Union{Real,Nothing}=nothing, dlnn::Union{Real,Nothing}=nothing, dlnTe::Union{Real,Nothing}=nothing, dlnTi::Union{Real,Nothing}=nothing, - P_perp_model::Symbol=:D_perp, + P_perp_model::Symbol=:D_perp, chi_perp_Dc::Union{Real,Nothing}=nothing, dr_val::Real=0.0, dgeo_val::Real=0.0, kpar_val::Union{Real,Nothing}=nothing, dc_type::Symbol=:none, ising::Integer=0, @@ -422,6 +428,8 @@ function slayer_parameters(; _need("chi_perp_e", chi_perp_e) P_perp_model === :tau_E || _need("chi_tor", chi_tor) P_perp_model in (:chi_perp_i, :D_perp) && _need("chi_perp_i", chi_perp_i) + chi_perp_Dc === nothing || (isfinite(chi_perp_Dc) && chi_perp_Dc > 0) || + throw(ArgumentError("slayer_parameters: chi_perp_Dc must be positive and finite when given; got $chi_perp_Dc")) tau_tor = rs^2 / chi_tor P_tor = tau_r / tau_tor if P_perp_model === :tau_E @@ -462,9 +470,9 @@ function slayer_parameters(; D_norm = (d_beta / rs) * lu^(1.0 / 3.0) * sqrt(iota_e) delta_n = lu^(1.0 / 3.0) / rs - # Critical-Δ offset from chi_parallel matching + # Critical-Δ offset from chi_parallel matching; its χ⊥,e is chi_perp_Dc when given dc_tmp = _solve_dc_tmp(; dc_type=dc_type, dr_val=dr_val, dgeo_val=dgeo_val, - chi_perp_e=chi_perp_e, t_e=t_e, zeff=zeff, + chi_perp_e=chi_perp_Dc === nothing ? chi_perp_e : chi_perp_Dc, t_e=t_e, zeff=zeff, tau_ee=tau_ee, rs=rs, R0=R0, sval_r=sval_r, n_tor=n, kpar_val=kpar_val) diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 0e0d6368a5..67db4626fb 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -95,6 +95,11 @@ constructor. χ⊥,i), `:chi_perp_e`, `:chi_perp_i`, `:P_phi`, `:c_beta` or `:tau_E`; see `InnerLayer.SLAYER.slayer_parameters`. `:tau_E` sets both P_perp and P_tor from τ̄_E, so it ignores `chi_phi`/`chi_tor`. `chi_perp_e` still sets the critical-Δ width in every mode + unless `chi_perp_Dc` is set + - `chi_perp_Dc` -- electron perpendicular heat diffusivity [m²/s] used only by the critical Δ + (default `nothing`: the critical Δ uses the same χ⊥,e as P_perp, from the kinetic file's + `chi_e` or `chi_perp_e`). When set, this scalar replaces χ⊥,e in the critical Δ only; + P_perp keeps the file's profiles or the `chi_perp_e`/`chi_perp_i` fallbacks. Must be positive - `D_R`, `D_geo` -- critical-Δ formula inputs. `nothing` (default) auto-derives them from the equilibrium: `D_R` from the resistive interchange index `D_R = E + F + H²` at each surface, `D_geo` from the @@ -178,6 +183,7 @@ file can keep the χ keys set to 0 to defer to the scalars). impurity_charge::Float64 = 6.0 chi_perp_e::Float64 = 1.0 chi_perp_i::Union{Float64,Nothing} = nothing + chi_perp_Dc::Union{Float64,Nothing} = nothing chi_tor::Float64 = 1.0 tau_E::Union{Float64,Nothing} = nothing P_perp_model::Symbol = :D_perp @@ -269,7 +275,8 @@ function validate(ctrl::SLAYERControl) ctrl.P_perp_model in _VALID_P_PERP_MODELS || throw(ArgumentError("SLAYERControl: P_perp_model=$(ctrl.P_perp_model) " * "not in $(_VALID_P_PERP_MODELS)")) - for (name, v) in (("chi_perp_e", ctrl.chi_perp_e), ("chi_perp_i", ctrl.chi_perp_i), ("chi_tor", ctrl.chi_tor)) + for (name, v) in (("chi_perp_e", ctrl.chi_perp_e), ("chi_perp_i", ctrl.chi_perp_i), ("chi_perp_Dc", ctrl.chi_perp_Dc), + ("chi_tor", ctrl.chi_tor)) v === nothing || (isfinite(v) && v > 0) || throw(ArgumentError("SLAYERControl: $name must be positive and finite [m²/s]; got $v")) end @@ -358,7 +365,7 @@ function slayer_control_from_toml(section::AbstractDict) kwargs[sym] = v isa Symbol ? v : Symbol(String(v)) elseif sym in (:Q_re_range, :Q_im_range) kwargs[sym] = _as_range(v) - elseif sym in (:bt, :D_R, :D_geo, :tau_E, :chi_perp_i) + elseif sym in (:bt, :D_R, :D_geo, :tau_E, :chi_perp_i, :chi_perp_Dc) # Allow explicit nothing (auto-derive) or a number (override) kwargs[sym] = v === nothing ? nothing : Float64(v) elseif sym === :omega_E_kHz diff --git a/src/Tearing/Runner/TearingProblem.jl b/src/Tearing/Runner/TearingProblem.jl index 4e4e1dd320..c23cd910fe 100644 --- a/src/Tearing/Runner/TearingProblem.jl +++ b/src/Tearing/Runner/TearingProblem.jl @@ -51,6 +51,7 @@ function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control) zeff=zeff, chi_perp_e=chi.chi_perp_e, chi_perp_i=chi.chi_perp_i, + chi_perp_Dc=chi.chi_perp_Dc, chi_tor=chi.chi_tor, tau_E=control.tau_E, P_perp_model=control.P_perp_model, diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 03bce89358..2e99129def 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -85,11 +85,14 @@ function _slayer_chi_inputs(loaded, control::SLAYERControl; warn::Bool=true) uses_tor && loaded.chi_tor === nothing && push!(fallbacks, "chi_phi → chi_tor = $(control.chi_tor)") isempty(fallbacks) || @warn "SLAYER: no usable kinetic-file profile (dataset absent or all-zero) for " * "$(join(fallbacks, ", ")); using the scalar fallback." + control.chi_perp_Dc !== nothing && control.D_c_type !== :none && + @info "SLAYER: the critical Δ uses chi_perp_Dc = $(control.chi_perp_Dc) m²/s for χ⊥,e; P_perp keeps its own χ inputs." model === :tau_E && @info "SLAYER: P_perp_model = tau_E sets P_perp and P_tor from a uniform " * "χ_E = (3/16)a²/τ̄_E with τ̄_E = $(control.tau_E) s; the χ_φ profile is not used." end return (; chi_perp_e=pick(loaded.chi_perp_e, control.chi_perp_e), chi_perp_i=pick(loaded.chi_perp_i, control.chi_perp_i), + chi_perp_Dc=control.chi_perp_Dc, chi_tor=pick(loaded.chi_tor, control.chi_tor)) end diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index a11fb1f3aa..adc45783fc 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -282,4 +282,23 @@ @test cb.dc_tmp == slayer_parameters(; P_perp_model=:chi_perp_e, merge(base, (dr_val=-1.0, dc_type=:lar))...).dc_tmp @test cb.P_tor == ref.P_tor end + + @testset "Test 4b: chi_perp_Dc sets only the critical Δ" begin + grads = (dlnn=-1.0, dlnTe=-2.0, dlnTi=-3.0) + for dc in (:lar, :fitzpatrick, :toroidal) + base = merge(_ref_kwargs(; dr_val=0.01, dc_type=dc), grads, (chi_perp_i=2.0, P_perp_model=:D_perp)) + a = slayer_parameters(; base...) + b = slayer_parameters(; base..., chi_perp_Dc=3.0) + # the critical Δ of b is the one chi_perp_e = 3 would give ... + @test b.dc_tmp == slayer_parameters(; merge(base, (chi_perp_e=3.0,))...).dc_tmp + @test b.dc_tmp != a.dc_tmp + # ... and every other layer quantity, P_perp included, is a's + for f in (:tau, :lu, :c_beta, :D_norm, :P_perp, :P_tor, :Q_e, :Q_i, :iota_e, :tauk, :tau_r, :delta_n, :eta, :d_beta) + @test getfield(b, f) == getfield(a, f) + end + end + base = merge(_ref_kwargs(; dr_val=0.01, dc_type=:fitzpatrick), (P_perp_model=:chi_perp_e,)) + @test_throws "chi_perp_Dc must be positive" slayer_parameters(; base..., chi_perp_Dc=0.0) + @test_throws "chi_perp_Dc must be positive" slayer_parameters(; base..., chi_perp_Dc=Inf) + end end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index ceb045863d..31b0a7f3c7 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -116,6 +116,12 @@ # Scalar diffusivities must be positive. @test_throws "chi_perp_i must be positive" Runner.validate(SLAYERControl(; chi_perp_i=-0.5)) @test_throws "chi_perp_e must be positive" Runner.validate(SLAYERControl(; chi_perp_e=0.0)) + # The critical-Δ diffusivity is optional, read from the deck, validated, and forwarded unchanged. + c_dc = slayer_control_from_toml(merge(section, Dict{String,Any}("chi_perp_Dc" => 1.4))) + @test c_dc.chi_perp_Dc == 1.4 && Runner.validate(c_dc) === c_dc + @test_throws "chi_perp_Dc must be positive" Runner.validate(SLAYERControl(; chi_perp_Dc=-1.0)) + @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; chi_perp_i=0.7, chi_perp_Dc=1.4); warn=false).chi_perp_Dc == 1.4 + @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; chi_perp_i=0.7); warn=false).chi_perp_Dc === nothing end @testset "TearingProblem solve forwards surfaces and Δ'" begin