diff --git a/docs/src/stability.md b/docs/src/stability.md index a423ecb4c..ca8ccc191 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -354,13 +354,25 @@ q = 8, and the resulting domain (`psilim ≈ 0.99498`, last surface q = 6) repro 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. `SLAYERModel`'s `chi_perp_Dc` (`[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. `SLAYERModel`'s `P_perp_model` (`[SLAYER] P_perp_model` in a deck) selects how `P_perp` is built: from `D_⊥` (`"D_perp"`, +`n T_s′/p′`. The critical-Δ island width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` is the +width above which parallel transport, through an electron parallel conductivity, makes the electron +temperature a flux function (Fitzpatrick, Phys. Plasmas **32**, 062509 (2025), Eqs. 230–234). That +paper and the `"lar"` and `"toroidal"` forms of Connor et al., Plasma Phys. Control. Fusion **57**, +065001 (2015), Eqs. 59 and 61, are single-fluid: they have one χ_⊥, which in Connor et al. diffuses +the total pressure. GPEC reads it as the electron channel and uses `chi_e` for every `D_c_type`, +since the width is set by electron heat transport; that reading is a modeling choice, not a result of +those papers. `SLAYERModel`'s `chi_perp_Dc` (`[SLAYER] chi_perp_Dc`, optional) sets that electron diffusivity for the critical Δ +alone, leaving `P_perp` on its own χ inputs. The deck key is one value for every surface, over +any `chi_e` profile; `build_slayer_inputs` also takes a function of ψ. It moves the critical +resonant field through the critical Δ. It is for a run where one χ⊥,e cannot serve +both formulas. `D_⊥` needs only the total conducted heat, so a run can use a one-fluid χ_eff in both +channels when an interpretive power balance cannot separate the species; where T_i ≈ T_e the +electron-ion exchange term makes that split uncertain, and it can return χ⊥,i < 0. On the electron reading, +the critical Δ still needs an electron diffusivity, which χ_eff is not, and `chi_perp_Dc` makes the choice +explicit. Where the split is uncertain the electron χ is uncertain too, so treat the value as a +modeling choice and scan it. The same split can often be reached without this key, by choosing +`chi_e` and `chi_i` to keep the `D_⊥` total, but not when the wanted `chi_e` would leave `chi_i ≤ 0`. +`SLAYERModel`'s `P_perp_model` (`[SLAYER] P_perp_model` in a deck) 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/KineticView.jl b/src/InnerLayer/SLAYER/KineticView.jl index b8b646baa..b404b6b7d 100644 --- a/src/InnerLayer/SLAYER/KineticView.jl +++ b/src/InnerLayer/SLAYER/KineticView.jl @@ -27,7 +27,8 @@ function _slayer_profiles(data::KineticProfileData, control) end # `control === nothing` reads no diffusivity, for callers that need only the densities, temperatures and Z_eff. model = control === nothing ? nothing : control.P_perp_model - uses_e = model in (:chi_perp_e, :D_perp) || (control !== nothing && control.D_c_type !== :none) + # With chi_perp_Dc set, the critical Δ takes it in place of χ⊥,e, so only P_perp can need the file's chi_e. + uses_e = model in (:chi_perp_e, :D_perp) || (control !== nothing && control.D_c_type !== :none && control.chi_perp_Dc === nothing) chi_perp_e = uses_e ? _chi_spline("chi_e", data.chi_e) : nothing chi_perp_i = model in (:chi_perp_i, :D_perp) ? _chi_spline("chi_i", data.chi_i) : nothing chi_tor = model !== nothing && model !== :tau_E ? _chi_spline("chi_phi", data.chi_phi) : nothing diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index 9a6e78382..a244eb779 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -52,7 +52,9 @@ a dispersion `solve_inner` path. reference (consumed only by `:toroidal`, which also takes its parallel wavenumber from the equilibrium, `toroidal_kpar`). An explicit `D_R = 0.0` disables the offset for every `D_c_type`; `D_geo = 0.0` only for `:toroidal`. - `chi_perp_Dc` - Electron perpendicular heat diffusivity [m²/s] used only by the critical Δ; `nothing` (default) - uses the χ⊥,e that P_perp sees. Set it when P_perp is built from a diffusivity the critical Δ should not use. + uses the χ⊥,e that P_perp sees. Set it when P_perp is built from a diffusivity the critical Δ should not use. It + replaces χ⊥,e in the critical Δ at every surface, over any `chi_e` profile, reaches the critical resonant field through + the critical Δ, and has no effect with `D_c_type = :none`. - `theta_sample` - Poloidal angle (normalized to [0, 1)) of the minor-radius chord of the `:midplane` radial label, used only by API callers that pass `rs_method = :midplane`; the default toroidal-flux label does not depend on angle, so a nonzero value warns (default 0, the outboard midplane). diff --git a/src/Tearing/SLAYERTearing.jl b/src/Tearing/SLAYERTearing.jl index 27144eb91..8b522013c 100644 --- a/src/Tearing/SLAYERTearing.jl +++ b/src/Tearing/SLAYERTearing.jl @@ -19,7 +19,7 @@ function _slayer_chi_inputs(loaded, control, data::KineticProfileData; warn::Boo model === :tau_E && data.tau_E === nothing && throw(ArgumentError("SLAYER: P_perp_model = tau_E needs the confinement time tau_E on the kinetic data")) if warn - uses_e = model in (:chi_perp_e, :D_perp) || control.D_c_type !== :none + uses_e = model in (:chi_perp_e, :D_perp) || (control.D_c_type !== :none && control.chi_perp_Dc === nothing) uses_i = model in (:chi_perp_i, :D_perp) uses_tor = model !== :tau_E fallbacks = String[] @@ -28,8 +28,13 @@ function _slayer_chi_inputs(loaded, control, data::KineticProfileData; warn::Boo uses_tor && loaded.chi_tor === nothing && push!(fallbacks, "chi_phi → chi_phi_default = $(data.chi_phi_default)") isempty(fallbacks) || @warn "SLAYER: no usable kinetic-file profile (dataset absent or all-zero) for " * "$(join(fallbacks, ", ")); using the default." - 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." + if control.chi_perp_Dc !== nothing + if control.D_c_type === :none + @warn "SLAYER: chi_perp_Dc has no effect with D_c_type = \"none\"; there is no critical Δ to set." + else + @info "SLAYER: the critical Δ uses chi_perp_Dc = $(control.chi_perp_Dc) m²/s for χ⊥,e at every surface; P_perp keeps its own χ inputs." + end + end 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 = $(data.tau_E) s; the χ_φ profile is not used." end diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index a0cc7d79f..851d5dbda 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -167,6 +167,31 @@ @test sl_var[1].P_perp ≈ sl_s[1].P_perp * 6.0 / 2.0 rtol = 1e-10 end + @testset "chi_perp_Dc reaches the critical Δ alone, as a callable and through TearingProblem" begin + sings = [_mk_sing(; psi=0.5, q=2.4, q1=1.2, m=2, n=1)] + kw = (; P_perp_model=:chi_perp_e, bt=2.0, chi_perp_e=2.0, chi_tor=1.5, dr_val=0.01, dc_type=:fitzpatrick) + plain = build_slayer_inputs(equil, sings, profiles; kw...)[1] + scalar = build_slayer_inputs(equil, sings, profiles; kw..., chi_perp_Dc=6.0)[1] + varying = build_slayer_inputs(equil, sings, profiles; kw..., chi_perp_Dc=psi -> 1.0 + 10.0 * psi)[1] # 6.0 at ψ = 0.5 + @test varying.dc_tmp == scalar.dc_tmp + @test scalar.dc_tmp != plain.dc_tmp + @test scalar.P_perp == plain.P_perp + + # The model key, through the same assembly TearingProblem uses. + Tearing = GeneralizedPerturbedEquilibrium.Tearing + loaded = (; profiles, chi_perp_e=nothing, chi_perp_i=nothing, chi_tor=nothing, zeff=nothing, zeff_source=:scalar) + model(; c...) = SLAYERModel(; P_perp_model=:chi_perp_e, D_R=0.01, D_c_type=:fitzpatrick, c...) + attached = equil.kinetic + equil.kinetic = Equilibrium.KineticProfileData(; psi=psi_pts, n_e=fill(5.0e19, length(psi_pts)), + T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), T_i=1000.0 .* (1.0 .- 0.6 .* psi_pts), + chi_e_default=2.0, chi_phi_default=1.5) + params(m) = Tearing._tearing_params(m, equil, sings, loaded)[1] + without, with = params(model()), params(model(; chi_perp_Dc=6.0)) + equil.kinetic = attached + @test with.dc_tmp != without.dc_tmp + @test with.P_perp == without.P_perp + 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) diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index b37732bf3..947a2e7c1 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -90,6 +90,13 @@ @test GeneralizedPerturbedEquilibrium.slayer_model(Dict{String,Any}("SLAYER" => Dict{String,Any}("chi_perp_Dc" => 1.4))).chi_perp_Dc == 1.4 @test SLAYERModel().chi_perp_Dc === nothing @test_throws "chi_perp_Dc must be positive" SLAYERModel(; chi_perp_Dc=-1.0) + # With it set, the critical Δ needs no χ⊥,e default; without a critical Δ, it does nothing and says so. + tau_E_model(; kw...) = SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal, kw...) + tau_data = data(; tau_E=0.1) + @test_logs (:warn, r"chi_e → chi_e_default") min_level = Base.CoreLogging.Warn Tearing._slayer_chi_inputs(no_chi, tau_E_model(), tau_data) + @test_logs min_level = Base.CoreLogging.Warn Tearing._slayer_chi_inputs(no_chi, tau_E_model(; chi_perp_Dc=1.4), tau_data) + @test_logs (:warn, r"chi_perp_Dc has no effect") min_level = Base.CoreLogging.Warn Tearing._slayer_chi_inputs(no_chi, + SLAYERModel(; P_perp_model=:tau_E, chi_perp_Dc=1.4), tau_data) end @testset "TearingProblem solve forwards surfaces and Δ'" begin @@ -122,6 +129,11 @@ @test bare_view.profiles.omega_E(data.psi[j]) == 0.0 @test_throws "missing required dataset 'T_i'" Tearing._slayer_profiles( GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_e=data.T_e), SLAYERModel()) + # A non-positive chi_e stops the run only where it is used: chi_perp_Dc replaces it in the critical Δ. + bad_chi_e = GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_e=data.T_e, T_i=data.T_i, + chi_e=fill(-1.0, length(data.psi))) + @test_throws "chi_e" Tearing._slayer_profiles(bad_chi_e, SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal)) + @test Tearing._slayer_profiles(bad_chi_e, SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal, chi_perp_Dc=1.4)).chi_perp_e === nothing end # Δ′ is unified across formalisms, so a Galerkin run feeds SLAYER exactly as a Riccati one