From fa7ecedd3015db6efdb57f125624d73f2d12a55a Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 29 Oct 2026 12:35:33 -0400 Subject: [PATCH 1/2] =?UTF-8?q?InnerLayer.SLAYER=20-=20FEATURE!=20-=20Opti?= =?UTF-8?q?onal=20chi=5Fperp=5FDc=20sets=20the=20critical-=CE=94=20electro?= =?UTF-8?q?n=20diffusivity=20apart=20from=20P=5Fperp?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The critical Δ's W_d and (χ∥/χ⊥)^(1/4) factors take χ⊥,e from the new optional chi_perp_Dc (slayer_parameters, build_slayer_inputs, [SLAYER] chi_perp_Dc); P_perp keeps its own χ inputs. Default nothing: unchanged behavior. Co-Authored-By: Claude Opus 5.5 --- docs/src/stability.md | 6 +++++- src/InnerLayer/SLAYER/LayerInputs.jl | 4 ++++ src/InnerLayer/SLAYER/LayerParameters.jl | 18 ++++++++++++----- src/Tearing/Runner/Control.jl | 11 +++++++++-- src/Tearing/Runner/TearingProblem.jl | 1 + src/Tearing/Runner/run_slayer.jl | 3 +++ test/runtests_slayer_params.jl | 25 ++++++++++++++++++++++++ test/runtests_slayer_runner.jl | 7 +++++++ 8 files changed, 67 insertions(+), 8 deletions(-) diff --git a/docs/src/stability.md b/docs/src/stability.md index 6ef6ebb639..08eb221ef5 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 f1e37610f2..5a6cc6ebb0 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -190,6 +190,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. @@ -255,6 +257,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, @@ -405,6 +408,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 b990f5d22a..07078977c6 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -188,7 +188,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, dc_type=:none, ising=0, resistivity_model=SauterNeoModel(), @@ -239,7 +239,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 - `dc_type` -- one of `:none`, `:lar`, `:fitzpatrick`, `:toroidal` @@ -296,7 +302,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, dc_type::Symbol=:none, ising::Integer=0, resistivity_model::NeoResistivityModel=SauterNeoModel(), @@ -403,6 +409,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 @@ -443,9 +451,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) diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 5dd2c5ccc1..5c5ff5545f 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 @@ there is one consistent interface for resistive and kinetic profiles. 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 @@ -266,7 +272,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 @@ -355,7 +362,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 a5c9aeb107..3501f22790 100644 --- a/src/Tearing/Runner/TearingProblem.jl +++ b/src/Tearing/Runner/TearingProblem.jl @@ -50,6 +50,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 26ba0c18cd..fdbb75ffc6 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -96,11 +96,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 f8f5a60810..554bcd57d5 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -240,4 +240,29 @@ @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 + # chi_perp_Dc equal to chi_perp_e reproduces the default exactly + same = slayer_parameters(; base..., chi_perp_Dc=1.0) + for f in (:P_perp, :D_norm, :dc_tmp) + @test getfield(same, 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=-1.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 3c07ff7dea..684d5122a4 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -118,6 +118,13 @@ # 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. + @test SLAYERControl().chi_perp_Dc === nothing + 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 From b6af97f48912b143466a121a722af300ae9793b4 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 30 Oct 2026 01:06:47 -0400 Subject: [PATCH 2/2] InnerLayer.SLAYER - TEST - Drop the chi_perp_Dc test cases the remaining checks already cover The equal-to-chi_perp_e block repeats what the default path and the harness pin, the negative value takes the same rejection branch as zero, and the default-nothing check is the last forwarding assertion's own case. Co-Authored-By: Claude Opus 5.5 --- test/runtests_slayer_params.jl | 6 ------ test/runtests_slayer_runner.jl | 1 - 2 files changed, 7 deletions(-) diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index e3677ae8fb..adc45783fc 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -296,15 +296,9 @@ 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 - # chi_perp_Dc equal to chi_perp_e reproduces the default exactly - same = slayer_parameters(; base..., chi_perp_Dc=1.0) - for f in (:P_perp, :D_norm, :dc_tmp) - @test getfield(same, 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=-1.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 0c6fbbbb71..31b0a7f3c7 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -117,7 +117,6 @@ @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. - @test SLAYERControl().chi_perp_Dc === nothing 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))