Tearing - BUGFIX - Stop the false chi_e fallback warning under chi_perp_Dc and test its forwarding

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 χ

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

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).

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

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)

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