From b83e1936be727577fe82ea6a6ea98056ef71684c Mon Sep 17 00:00:00 2001 From: Evan Bursch Date: Wed, 30 Sep 2026 11:07:04 -0400 Subject: [PATCH 1/6] Equilibrium - FEATURE - Read optional chi_i ion heat diffusivity from kinetic HDF5 Add an optional chi_i (ion perpendicular heat diffusivity) dataset to the kinetic-profile reader and writer, alongside chi_e and chi_phi. Ship chi_i in the DIII-D-like example kinetic file, built like its chi_e: AOT CHI_I of DIII-D 147131 at 2300 ms, mapped rho_N -> psi_N with EFIT01 RHOVN, the spurious zero at psi_N = 1 dropped, Gaussian-smoothed with sigma = 12 of 201 rho_N points (the width that best reproduces the example chi_e from AOT CHI_E). Range 0.46-2.1 m^2/s. Co-Authored-By: Claude Opus 5.5 --- .../TkMkr_D3Dlike_Hmode_kinetic.h5 | Bin 57472 -> 67728 bytes src/Equilibrium/KineticProfiles.jl | 15 +++++++++------ test/runtests_kinetic.jl | 8 ++++++-- 3 files changed, 15 insertions(+), 8 deletions(-) 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 2f084f080315cfd157a6c6ed58ddd06cd72bf2e3..c1a305586b1b8c69213a9133ac9473a0b3be5c62 100644 GIT binary patch delta 4650 zcmY*dc{G&m|3)NRycMa?v@fy?i7wd{!Ye9eUy}$WUNxabNr)CymR4K#l%2FnX_FT6 z%oux&nK8yP(|6wGoZtD~|J=`KxvuN^e4fv_&$*vdO>>e{gp=3LNuDQ|qU9~E3`X|~<&Z%4_0?tcvab23rhEEFy%D3Yox94T)3U-VQ@N~uz6iSYWA z#^b3*RtL)|B5Vt(DIo!YQ|o}rurJBf%F0T`l%%1ptD&Vq@(3qsYN@M}g2G7aL;Zu@Y)L^LAt65A zr0`%;NLT>hjtB_%Cz+XBIjU=r94vP^*={0*g&YkH2tFD{it-5!{+r^TOhXNlcWCh6 z;7ue85AxBlumF$4Dq%sv!QuXV!hd+e0{ntJNQNX0&A(qolTUf{aDZ=cXb|aN@sQwv z!{K3@NPGkz;LTSK_viC_`h`Gk@FuFdEF*Jidsm2`NMBgrGw=O1)x z!iP2U{%y!haekNT#sBW0(5Zpc+8uU6$(+>gISwg}^HSS|jQ;qiRs3ARfAa<0Q%apv z)rIE@3a5S%&MB8!dR>H^OTS+C<7O9X8EUl`s9os%5HFsl(G6Ft2*2^O-T3{uSJIu* zjVzfrbsok)F|sVWOziGYxObgeswLb5+V<)tW*R-vDl=Ipx1$F^=4^pQu08m?JF%eT zNDsOlLzfY!dvIujL3`@89z?X7&1=l-frb5z-UBr~WPEi$99h%a1NX;LdTYM+0EN>U z5wsrc)}!n^Fy4dcg+((%vwY%x)q@;S0#)(S64w_HaAMxQ_+5d(A3Iad)~_XSapuZI zgDQdhYa~SDbO==DS;#O92?!kYOPOazAoko!-SM3ShQ^4@KnDUF>x@(K_mBy+ZvJ^8 z*OdV2ZnoDOF9OVJ%f{#a1oHAT1HFO>c;4zyryt=9O22)mb&SBHc`2XwMezl~EG2fH zB%rgUyNDA*pyNXDG4C@3ZmEdA%!(%h9}}<;QcW(Z=C}4>*tr+A1e$M&W%t(;Xsl}OSlmRw z-FTI=J`V=Vp5X9QkJRFD}j2<%N;*R%gMUufot|F<^;L{1$3z4sk~Y;ohn>GuRC zj9)ird?Fx5xT+lZLO>$yyYAL+1WIpdH7Rrvc)2NtJY~>Bpu765#sMk;dW*v0YkdSv z?b{^>2Kg0t#TDrNB4C(R;nzVY(5L7t9>pTy|6Yu_g+su`b8><1EP>G%9feu~6uc7c zTfaw`0!NB+LE=0LRMTg(nk6X+ZFVv^@COB%N<|hNSqeU!Zj>ooO2PhlAl&Sf{w`OsH%`?RXK@@~r9SFA$r$A?? zc;}hFDDWXatJyLYOF_D6)!f1h6qJ=rwja36w-u5|4jKHt;aO8yE(N!h{KK#2QxI_F zIx3v~c7!`25pt@{gqTJ>U7 zms?-Fw->AkM_NokFFN!6j!+_dp>>1#)Aa%$AMy2gR|el-BQn`n(2KvUa#*Hyxq-9Sv6BHZ=%()8KTUw{S9ohP*i=qGJg(Xy^&b|CLR{6Gd%Fn@SqyEIU!V;x!HW zC&;O8KWW(b+1p=_JV8V3z5Zp3B>JH6%$q5;q7N_CuH3BD?Sov%B4cmcK0N(qIwI}U z2fe5vrHtr4L@i$O#XhYMCBIt>B@6qI6I7gL)zXJGAC{4>+Mv^=ABO@RWXUPV`r#4NYDrG%$Kq~PW3z|Uz@e9Vv($b>&)R>M78t;G;=>JVg#q|pwfbJLWdQO5H@1m758!pcXpU$20Q!!& zP;OryK%G@|ky_~hP7cgDR`qTG2hMdh)Q=9pe0jBJ>ij_r6}|uBr8^QRx<`$%E+On3y+~4&sN*qT5X$2l0p7dnpt4Ak2y`9dD5z!rWyfg^aC3 zXemE%?1tMAPFZKR--;W8)T%<--TOm0y759+-J2oM98Bm-7(-Z6t2Z3Aco?fids!id z!&rZOd9{c4Fg6X6_mLeEhS5vwl-C6D0dWz!1Jyb(k^{eH6W;|TN{at7xrlH z&U8Ng3tM%!XNi>kLd<^l?V_K*(43+FqDE>I!(^sTkAdMR4&KhMtvWP{K?kLD%^Rc0 z(g@%%-;6?jfr*PKZxl%kh5DsxV~AO}(6!Zl4AD0;t~aNS;hnDO%JS!9c-_>xKYwZr zm4ZPo%8TiU)rob=(V$~j;5GZj)^r?WNZ8K$(h*fEH#a1fj*}rP>Zv#BAZK|Q|Kak4 z4$VJ1&enXUBQ=@cv4TOzv=fg=m_LqIDf+EI<#A*;cSN-*;|NF;e(xVQ z4vT76m4`XwSb6(;++@Qz1lk!t4Z6pn^tNN|TJAU~@~f4i6c`xT^F#B)Rt7W>qh;;F zz~-)9PgDM4kO6J8B|1PqNIm@fcVo+!T6fITH^K3CVpn zVq&iEG#lKt*`H6L2t9i2~~*+806e0 zvDGH9|J;J)6}u-Op|g*$kD9>CSH*g#vnF7qA$C{i*#sV4_xxjUXaX18cOA)=VMDtp z=-^cyHVz$jyj8iMjib%-pT5Sj@m~H#qG%o)gZF=(4SmPPyi8gGg*?H=%8im|Y!oJu z@kIU4LE}k?2AuOw517Q)b$2{WlP9r)muoLnHHrJhnw5iHlPJ+J3rI61)&p9ZrA{Txd;b2V?En?fkDX=z(bf#{Z z!X2NNrR1Y-Q}|-;ec9p66r|P*gsBxwA?N4E?9BEle9|pE|9xT#BFU#imo1xyaaTmA zq|r3&^K+gCcu&JUqV?hP_-Pb>E^=RAG>vO7#ohL{Pb0|DpfH{}4Wpi*sVcb{Xg;=T zkl!+cJ(Mgv_P!aMkYp1dqi0}vUVtoNbZZ7Lwj`AK*U!LVNOoUJ_Y7{^&gj?kX3)_5 z?2)gBl@(%LAirp`rV)Q)2z7F^tOA$Pud=*oqh%^dc!a4xv? za)F?DF3^9~z2*iNZuIa;`4TQj#iHE}WGvR9-00~g2dG##Z-xwz`Cb@1Id7eDR$ zH0yu!@lt!%6-e>WaDV$YmOKyoO?@joSMgv|Aa?()3J(<~w*?OC@Q{?WLQC0@2XXSk zRdM|$Jp8i1Tvlhr!>D3KWUe(2vgKEPp54WRp=E?wpaT#2_4>7k?0L8}Uu|-~9S=ub zS83Yt<70t}i&yXDL22(6h?(<+%y;{Z8u6ej^lCVH3l9!vnW0YVJmiX`sVT1G;kZ)P zEY%1QG;Zg%Z%suSd)^O$Vy6DAk_5+x^WDO{wsRmRiWxu{R7NsVscV$;{O zSXLevMdTklUDCKvTf17e`6L&XTdO2~`EW5P)N<;mH5c)(+b6i1T&%j(wxUapi+~^7 z8_!P7qIBaUSM?vWSiiPy`}^8iTquZ|Nh~xyq0NbG?OyouT4W{|9lJigVSjH9YZ|AdN70_m9`(pyGp_@Z@ zWE1kH9!vT$mIBSS8rGLe@{{~=JrM0=ZTlfqRdYDCkF9v)#rR`FhH zf(;XP?!_){D*X5c?tWb|+Lg<=se5&t_nhV=i+d zfr0DY;k&(#F>v-zsTR5q11r78ZjC!JkX047+R~DNQ=$jSsT+(KxR&WME~dvotINls zSWO1Ft>F#z8Vsl!$fZrIGq7+?)0w@S88}M2ap#Z$1HwELPL>G+y-O-yj5;uIUbHZ6 zt0x2bb>H&a!x`9eLTTffiwum>%P#8NVc_?Ot&7l820mR(ahahq;QD=J;@Mm#*l}e2 zPoLG9Sa5}+H@uIDdS3(WrVC6^yM=$I*D(<^ZppYd&V=X4PKAA{EbQ;MHCG5^p)+IB z<6996p38bUW-J!}^4T9PZ7_kGbLnZ%E=(Xpq;SsJuM_x?cu)MJDjOH)+r>AWWy7yI ziat!Rk@97+)H(hvRFFq1Hmtlqi4kp+h?6q>Gdj5Xx&wZOjt<9}P?Z9RL6T delta 125 zcmbO*k)`1v^8^j%h6fB2wUZecC-y5ca%^Ul`pzij;p!N`00a>r0th@H6vGLo$z2cq zHVZJTF*0&YHfO1ve1k<}vd5zZlN~&>Ha}wh=FZE`#K6#S0ivD}Xj%Z%X48P;CQhgc M3^z6}JanEH0QZ_8xc~qF diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 934256554..18037d7e7 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 χ⊥ | +| `chi_e` | electron perpendicular heat diffusivity χ⊥,e | +| `chi_i` | ion perpendicular heat diffusivity χ⊥,i | | `chi_phi` | toroidal momentum diffusivity χ_φ | | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras | | `provenance` | short string recording the source file/format | @@ -80,17 +81,18 @@ 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}} species_densities::Union{Nothing,Dict{String,Vector{Float64}}} provenance::String 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, + omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_i=nothing, chi_phi=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), species_densities, String(provenance)) + _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), species_densities, String(provenance)) end const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") @@ -98,7 +100,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") + "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s") """ read_kinetic_file(path; group="/") -> KineticProfileData @@ -137,7 +139,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") + standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_i", "chi_phi") extras = Dict{String,Vector{Float64}}() for k in keys(g) (k in standard || !startswith(k, "n_")) && continue @@ -147,7 +149,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"), + chi_e=rd("chi_e"), chi_i=rd("chi_i"), chi_phi=rd("chi_phi"), species_densities=(isempty(extras) ? nothing : extras), provenance=prov) end end @@ -178,6 +180,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) # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip. if data.species_densities !== nothing diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index d502385f7..a8a8a8399 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -605,8 +605,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 @@ -634,6 +635,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 From 3fd1885da8db0144f17b43fe54ae34631e4a590b Mon Sep 17 00:00:00 2001 From: Evan Bursch Date: Wed, 30 Sep 2026 11:07:04 -0400 Subject: [PATCH 2/6] InnerLayer.SLAYER - FEATURE! - Add P_perp_model options and rename chi_perp to chi_perp_e New [SLAYER] P_perp_model selects how the perpendicular Prandtl number is built: - chi_perp_e (default, unchanged numbers): tau_R*chi_perp_e/r_s^2 - chi_perp_i: tau_R*chi_perp_i/r_s^2 (kinetic-file chi_i, scalar fallback) - P_phi: P_perp = P_tor (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143) - D_perp: tau_perp = r_s^2/D_perp (Fitzpatrick, Phys. Plasmas 29, 032507 (2022)) - c_beta: C^2 = P_perp with C = c_beta (Park et al., Phys. Plasmas 29, 122505 (2022)) Breaking: the [SLAYER] chi_perp key is renamed chi_perp_e. chi_perp_e still sets the critical-Delta offset in every mode. Regression harness diiid_slayer_n1 develop vs local: 17 unchanged. Co-Authored-By: Claude Opus 5.5 --- examples/DIIID-like_SLAYER_example/gpec.toml | 2 +- src/InnerLayer/SLAYER/LayerInputs.jl | 17 +++++-- src/InnerLayer/SLAYER/LayerParameters.jl | 53 ++++++++++++++++---- src/Tearing/Runner/Control.jl | 28 ++++++++--- src/Tearing/Runner/run_slayer.jl | 30 +++++++---- test/runtests_dispersion_coupled.jl | 2 +- test/runtests_dispersion_residual.jl | 2 +- test/runtests_slayer_inputs.jl | 8 +-- test/runtests_slayer_params.jl | 19 ++++++- test/runtests_slayer_riccati.jl | 2 +- test/runtests_slayer_runner.jl | 7 +++ 11 files changed, 127 insertions(+), 43 deletions(-) diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index db8aec602..3dedd663c 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -73,7 +73,7 @@ coupling_mode = "uncoupled" # "uncoupled" (per-surface) or "coupled" ( dc_type = "none" # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal" mu_i = 2.0 # Ion mass in proton-mass units (2.0 = deuterium) zeff = 1.0 # Effective charge -chi_perp = 1.0 # fallback only; the kinetic file supplies χ⊥(ψ) +chi_perp_e = 1.0 # fallback only; the kinetic file supplies χ⊥,e(ψ) chi_tor = 1.0 # fallback only; the kinetic file supplies χ_φ(ψ) pole_threshold_adaptive = true # SLAYER |Δ| spans many orders; adapt to 10·median(|Δ|) filter_above_poles = true # Discard roots above the highest pole γ diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 252e81ae9..a59a1a4a1 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -181,8 +181,13 @@ profiles, without an intermediate file round-trip. - `mu_i` -- ion mass in proton-mass units (default `2.0` for D). - `zeff` -- effective charge (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`, `D_perp` -- ion perpendicular heat / perpendicular particle + diffusivity [m²/s], scalar or callable of `psi` (default `1.0`); used only by + the matching `P_perp_model`. + - `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²` @@ -232,7 +237,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; mu_i::Real=2.0, zeff::Real=1.0, z_i::Real=1.0, - chi_perp=1.0, + chi_perp_e=1.0, + chi_perp_i=1.0, + D_perp=1.0, + P_perp_model::Symbol=:chi_perp_e, chi_tor=1.0, dr_val=nothing, dgeo_val=nothing, @@ -378,8 +386,11 @@ 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=zeff, - 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), + D_perp=_eval(D_perp, psi), + 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 ae875dea4..c68e7f4cc 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -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, :rfitzp, :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) + """ r_based_shear(rs, q, dq_dpsi, da_dpsi) -> Float64 @@ -132,7 +135,7 @@ end # Internal: solve the Wd self-consistency loop for the chi_parallel-based # critical Δ (Connor-Hastie-Helander 2015). 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 || @@ -149,7 +152,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 @@ -165,13 +168,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 === :rfitzp 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 @@ -179,8 +182,10 @@ 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=1.0, D_perp=1.0, + P_perp_model=:chi_perp_e, dr_val=0.0, dgeo_val=0.0, dc_type=:none, ising=0, resistivity_model=SauterNeoModel(), @@ -210,7 +215,16 @@ 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_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`) + - `D_perp` -- perpendicular particle diffusivity [m²/s] (used by `P_perp_model=:D_perp`) + - `P_perp_model` -- how `P_perp` is built: + + `:chi_perp_e` (default) -- τ_R·χ⊥,e/r_s² + + `: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_⊥ (Fitzpatrick, Phys. Plasmas 29, 032507 (2022)) + + `:c_beta` -- P_perp = C² with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) + `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`, `:rfitzp`, `:toroidal` @@ -263,8 +277,10 @@ 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=1.0, D_perp::Real=1.0, + 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(), @@ -366,10 +382,25 @@ 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 + 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_⊥ = P_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143) + P_tor + elseif P_perp_model === :D_perp + # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥ (Fitzpatrick, Phys. Plasmas 29, 032507 (2022)) + tau_r / (rs^2 / D_perp) + elseif P_perp_model === :c_beta + # C² = P_⊥ with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) + c_beta^2 + else + throw(ArgumentError("slayer_parameters: unknown P_perp_model=$P_perp_model. " * + "Allowed: $(ALLOWED_P_PERP_MODELS)")) + end # Normalized beta-related width and Δ-normalization d_beta = c_beta * d_i @@ -378,7 +409,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/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 1a12ff7e6..728da81c0 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -34,10 +34,15 @@ constructor. callable of `psi` overrides it - `mu_i` -- ion mass in proton-mass units (default 2.0 for D) - `zeff` -- effective charge - - `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 + - `D_perp` -- perpendicular particle diffusivity [m²/s] (default 1.0), used + only by `P_perp_model = :D_perp` + - `P_perp_model` -- how the perpendicular Prandtl number is built: `:chi_perp_e` + (default), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:D_perp`, or `:c_beta` + (P_perp = c_β²); see `InnerLayer.SLAYER.slayer_parameters` - `dr_val`, `dgeo_val` -- critical-Δ formula inputs. `nothing` (default) auto-derives them from the equilibrium: `dr_val` from the resistive interchange index `D_R = E + F + H²` at each surface, `dgeo_val` from the @@ -84,9 +89,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 `"/"`) @@ -107,8 +112,11 @@ there is one consistent interface for resistive and kinetic profiles. bt::Union{Float64,Nothing} = nothing mu_i::Float64 = 2.0 zeff::Float64 = 1.0 - chi_perp::Float64 = 1.0 + chi_perp_e::Float64 = 1.0 chi_tor::Float64 = 1.0 + chi_perp_i::Float64 = 1.0 + D_perp::Float64 = 1.0 + P_perp_model::Symbol = :chi_perp_e dr_val::Union{Float64,Nothing} = nothing dgeo_val::Union{Float64,Nothing} = nothing theta_sample::Float64 = 0.0 @@ -161,6 +169,7 @@ const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_shooting, :ggj_galerkin) const _VALID_SCAN_MODES = (:amr, :brute_force) const _VALID_COUPLING_MODES = (:uncoupled, :coupled) const _VALID_DC_TYPES = (:none, :lar, :rfitzp, :toroidal) +const _VALID_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta) const _VALID_RESISTIVITY_MODELS = (:sauter, :redl, :spitzer, :spitzer_harm) const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson) @@ -177,6 +186,9 @@ function validate(ctrl::SLAYERControl) ctrl.dc_type in _VALID_DC_TYPES || throw(ArgumentError("SLAYERControl: dc_type=$(ctrl.dc_type) " * "not in $(_VALID_DC_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.resistivity_model in _VALID_RESISTIVITY_MODELS || throw(ArgumentError("SLAYERControl: resistivity_model=$(ctrl.resistivity_model) " * "not in $(_VALID_RESISTIVITY_MODELS)")) @@ -242,7 +254,7 @@ function slayer_control_from_toml(section::AbstractDict) for (k, v) in flat sym = Symbol(k) if sym in (:inner_model, :scan_mode, :coupling_mode, :dc_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) diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2459087cf..72c65a856 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -18,9 +18,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`). +# `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`). function _load_profiles(control::SLAYERControl, dir_path::AbstractString) isempty(control.profile_file) && error("run_slayer: [SLAYER] profile_file is empty — point it at a " * @@ -46,14 +46,15 @@ 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) - return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor) + return (profiles=profiles, chi_perp_e=chi_perp_e, chi_perp_i=chi_perp_i, chi_tor=chi_tor) end # --------------------------------------------------------------------- @@ -407,7 +408,7 @@ Orchestrate the full SLAYER analysis against a `ForceFreeStates.ForceFreeStatesR reading its equilibrium, singular surfaces and Δ' matrix. Kinetic profiles are read from `control.profile_file` (relative to `dir_path`) through the shared `Equilibrium.read_kinetic_file` reader; when the file carries `chi_e`/`chi_phi` -profiles they set χ⊥(ψ)/χ_φ(ψ), otherwise the scalar `control.chi_perp`/ +profiles they set χ⊥(ψ)/χ_φ(ψ), otherwise the scalar `control.chi_perp_e`/ `chi_tor` fallbacks are used. The toroidal field comes from `control.bt`; leaving it unset (the default) makes @@ -454,18 +455,25 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac # 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_perp_e = loaded.chi_perp_e === nothing ? control.chi_perp_e : loaded.chi_perp_e chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor - (loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn( + (loaded.chi_perp_e === nothing || loaded.chi_tor === nothing) && @warn( "SLAYER: kinetic file has no usable chi_e/chi_phi profile(s) " * "(dataset absent or all-zero); using the scalar " * - "control.chi_perp/chi_tor fallback for the missing one(s).") + "control.chi_perp_e/chi_tor fallback for the missing one(s).") + chi_perp_i = loaded.chi_perp_i === nothing ? control.chi_perp_i : loaded.chi_perp_i + (control.P_perp_model === :chi_perp_i && loaded.chi_perp_i === nothing) && @warn( + "SLAYER: kinetic file has no usable chi_i profile; using the scalar " * + "control.chi_perp_i fallback.") params = build_slayer_inputs(equil, surfaces, profiles; bt=bt, mu_i=control.mu_i, zeff=control.zeff, - chi_perp=chi_perp, + chi_perp_e=chi_perp_e, chi_tor=chi_tor, + chi_perp_i=chi_perp_i, + D_perp=control.D_perp, + P_perp_model=control.P_perp_model, dr_val=control.dr_val, dgeo_val=control.dgeo_val, dc_type=control.dc_type, 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 e8af06da4..2d68531a2 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_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index fce001036..26b51932c 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -142,24 +142,24 @@ @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 diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index 1f665fb1d..73b1a222d 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,19 @@ @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 + @test slayer_parameters(; base..., P_perp_model=:D_perp, D_perp=3.0).P_perp ≈ 3 * ref.P_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) + # 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 ee114e37f..96e706786 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -57,6 +57,8 @@ "bt" => 1.8, "mu_i" => 2.0, "dr_val" => 0.01, + "P_perp_model" => "D_perp", + "D_perp" => 0.5, "scan_grid" => Dict{String,Any}( "Q_re_range" => [-5.0, 5.0], "Q_im_range" => [-1.0, 3.0], @@ -79,6 +81,8 @@ @test c.msing_max == 2 @test c.bt === 1.8 @test c.dr_val == 0.01 + @test c.P_perp_model === :D_perp + @test c.D_perp == 0.5 @test c.Q_re_range == (-5.0, 5.0) @test c.Q_im_range == (-1.0, 3.0) @test c.nre == 50 @@ -92,6 +96,9 @@ # 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 ArgumentError 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")))) end @testset "run_slayer: result-facing form forwards surfaces and Δ'" begin From 80c6664aa1989d1d1367c66c551f8b6d9e1847e6 Mon Sep 17 00:00:00 2001 From: Evan Bursch Date: Wed, 30 Sep 2026 11:22:24 -0400 Subject: [PATCH 3/6] InnerLayer.SLAYER - FEATURE - Add tau_E P_perp_model (P_perp = P_phi = tau_R/tau_E) A single [SLAYER] tau_E [s] for the whole plasma sets P_perp = P_tor = tau_R/tau_E at every surface (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131). validate() requires a positive tau_E when this model is selected. Co-Authored-By: Claude Opus 5.5 --- src/InnerLayer/SLAYER/LayerInputs.jl | 4 ++++ src/InnerLayer/SLAYER/LayerParameters.jl | 16 ++++++++++++---- src/Tearing/Runner/Control.jl | 14 ++++++++++---- src/Tearing/Runner/run_slayer.jl | 1 + test/runtests_slayer_params.jl | 4 ++++ test/runtests_slayer_runner.jl | 3 +++ 6 files changed, 34 insertions(+), 8 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index a59a1a4a1..fb0ff7d46 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -186,6 +186,8 @@ profiles, without an intermediate file round-trip. - `chi_perp_i`, `D_perp` -- ion perpendicular heat / perpendicular particle diffusivity [m²/s], scalar or callable of `psi` (default `1.0`); used only by the matching `P_perp_model`. + - `tau_E` -- whole-plasma energy confinement time [s], one scalar for every surface; + required by `P_perp_model=:tau_E`. - `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 @@ -240,6 +242,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; chi_perp_e=1.0, chi_perp_i=1.0, D_perp=1.0, + tau_E::Real=NaN, P_perp_model::Symbol=:chi_perp_e, chi_tor=1.0, dr_val=nothing, @@ -390,6 +393,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; chi_tor=_eval(chi_tor, psi), chi_perp_i=_eval(chi_perp_i, psi), D_perp=_eval(D_perp, psi), + tau_E=tau_E, P_perp_model=P_perp_model, m=m_res, n=n_res, dr_val=dr_val_k, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index c68e7f4cc..a174e8c11 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -105,7 +105,7 @@ end const ALLOWED_DC_TYPES = (:none, :lar, :rfitzp, :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) +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 @@ -184,7 +184,7 @@ end qval, sval_r, bt, rs, R0, mu_i, zeff, chi_perp_e, chi_tor, m, n, - chi_perp_i=1.0, D_perp=1.0, + chi_perp_i=1.0, D_perp=1.0, tau_E=NaN, P_perp_model=:chi_perp_e, dr_val=0.0, dgeo_val=0.0, dc_type=:none, ising=0, @@ -218,12 +218,14 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `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`) - `D_perp` -- perpendicular particle diffusivity [m²/s] (used by `P_perp_model=:D_perp`) + - `tau_E` -- whole-plasma energy confinement time [s] (required by `P_perp_model=:tau_E`) - `P_perp_model` -- how `P_perp` is built: + `:chi_perp_e` (default) -- τ_R·χ⊥,e/r_s² + `: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_⊥ (Fitzpatrick, Phys. Plasmas 29, 032507 (2022)) + `:c_beta` -- P_perp = C² with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) + + `:tau_E` -- P_perp = P_tor = τ_R/τ_E, overriding χ_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131) `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 @@ -279,7 +281,7 @@ function slayer_parameters(; rs::Real, R0::Real, mu_i::Real, zeff::Real, chi_perp_e::Real, chi_tor::Real, m::Integer, n::Integer, - chi_perp_i::Real=1.0, D_perp::Real=1.0, + chi_perp_i::Real=1.0, D_perp::Real=1.0, tau_E::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, @@ -384,11 +386,17 @@ function slayer_parameters(; # Effective Prandtl-like transport ratios tau_tor = rs^2 / chi_tor P_tor = tau_r / tau_tor + if P_perp_model === :tau_E + isfinite(tau_E) && tau_E > 0 || + throw(ArgumentError("slayer_parameters: P_perp_model=:tau_E needs a positive tau_E, got $tau_E")) + # P_⊥ = P_φ = τ_R/τ_E (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131) + P_tor = tau_r / tau_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 + 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 diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 728da81c0..61aa879c7 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -40,9 +40,12 @@ constructor. all-zero); otherwise the file's χ⊥,e(ψ)/χ⊥,i(ψ)/χ_φ(ψ) take precedence - `D_perp` -- perpendicular particle diffusivity [m²/s] (default 1.0), used only by `P_perp_model = :D_perp` + - `tau_E` -- whole-plasma energy confinement time [s], required by + `P_perp_model = :tau_E` (default `nothing`) - `P_perp_model` -- how the perpendicular Prandtl number is built: `:chi_perp_e` - (default), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:D_perp`, or `:c_beta` - (P_perp = c_β²); see `InnerLayer.SLAYER.slayer_parameters` + (default), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:D_perp`, `:c_beta` + (P_perp = c_β²), or `:tau_E` (P_perp = P_tor = τ_R/τ_E); see + `InnerLayer.SLAYER.slayer_parameters` - `dr_val`, `dgeo_val` -- critical-Δ formula inputs. `nothing` (default) auto-derives them from the equilibrium: `dr_val` from the resistive interchange index `D_R = E + F + H²` at each surface, `dgeo_val` from the @@ -116,6 +119,7 @@ there is one consistent interface for resistive and kinetic profiles. chi_tor::Float64 = 1.0 chi_perp_i::Float64 = 1.0 D_perp::Float64 = 1.0 + tau_E::Union{Float64,Nothing} = nothing P_perp_model::Symbol = :chi_perp_e dr_val::Union{Float64,Nothing} = nothing dgeo_val::Union{Float64,Nothing} = nothing @@ -169,7 +173,7 @@ const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_shooting, :ggj_galerkin) const _VALID_SCAN_MODES = (:amr, :brute_force) const _VALID_COUPLING_MODES = (:uncoupled, :coupled) const _VALID_DC_TYPES = (:none, :lar, :rfitzp, :toroidal) -const _VALID_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta) +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) @@ -189,6 +193,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)")) + ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && ctrl.tau_E > 0) || + throw(ArgumentError("SLAYERControl: P_perp_model=tau_E requires a positive tau_E [s]")) ctrl.resistivity_model in _VALID_RESISTIVITY_MODELS || throw(ArgumentError("SLAYERControl: resistivity_model=$(ctrl.resistivity_model) " * "not in $(_VALID_RESISTIVITY_MODELS)")) @@ -258,7 +264,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, :dr_val, :dgeo_val) + elseif sym in (:bt, :dr_val, :dgeo_val, :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/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 72c65a856..c8fb91f3b 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -473,6 +473,7 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac chi_tor=chi_tor, chi_perp_i=chi_perp_i, D_perp=control.D_perp, + tau_E=something(control.tau_E, NaN), P_perp_model=control.P_perp_model, dr_val=control.dr_val, dgeo_val=control.dgeo_val, diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index 73b1a222d..d34b05759 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -209,6 +209,10 @@ @test slayer_parameters(; base..., P_perp_model=:D_perp, D_perp=3.0).P_perp ≈ 3 * ref.P_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) + te = slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1) + @test te.P_tor == te.tau_r / 0.1 + @test te.P_perp == te.P_tor + @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:tau_E) # 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 diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 96e706786..cf4bc5b87 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -99,6 +99,9 @@ # The renamed χ⊥,e key and bad P_perp_model values are rejected @test_throws ArgumentError 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")))) + 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 end @testset "run_slayer: result-facing form forwards surfaces and Δ'" begin From 915447fff7fe76e8cc95503a14c46771a2df178b Mon Sep 17 00:00:00 2001 From: Evan Bursch Date: Wed, 30 Sep 2026 11:28:41 -0400 Subject: [PATCH 4/6] InnerLayer.SLAYER - FEATURE - Compute D_perp from Fitzpatrick 2022 Eq. 16 P_perp_model = D_perp now builds D_perp per surface from Eq. 16 of Fitzpatrick, Phys. Plasmas 29, 032507 (2022), with eta_perp = eta_par: D_perp = c_beta^2 eta/mu0 + (2/3)(1 - c_beta^2)[eta_e tau_e/(1+eta_e) chi_perp_e + eta_i tau_i/(1+eta_i) chi_perp_i] using the n, T_e, T_i log-gradients at the surface (Eqs. 4-6, written flat-density safe). The scalar D_perp input is removed. Co-Authored-By: Claude Opus 5.5 --- src/InnerLayer/SLAYER/LayerInputs.jl | 19 ++++++++---------- src/InnerLayer/SLAYER/LayerParameters.jl | 25 +++++++++++++++++++----- src/Tearing/Runner/Control.jl | 5 +---- src/Tearing/Runner/run_slayer.jl | 3 +-- test/runtests_slayer_params.jl | 10 +++++++++- test/runtests_slayer_runner.jl | 2 -- 6 files changed, 39 insertions(+), 25 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index fb0ff7d46..a8530573f 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -183,9 +183,8 @@ profiles, without an intermediate file round-trip. - `zeff` -- effective charge (default `1.0`). - `chi_perp_e` -- electron perpendicular heat diffusivity [m²/s]. Scalar or a callable of `psi` (default `1.0`). - - `chi_perp_i`, `D_perp` -- ion perpendicular heat / perpendicular particle - diffusivity [m²/s], scalar or callable of `psi` (default `1.0`); used only by - the matching `P_perp_model`. + - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], scalar or callable of + `psi` (default `1.0`); used by `P_perp_model = :chi_perp_i` and `:D_perp`. - `tau_E` -- whole-plasma energy confinement time [s], one scalar for every surface; required by `P_perp_model=:tau_E`. - `P_perp_model` -- how `P_perp` is built; see [`slayer_parameters`](@ref) @@ -241,7 +240,6 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; z_i::Real=1.0, chi_perp_e=1.0, chi_perp_i=1.0, - D_perp=1.0, tau_E::Real=NaN, P_perp_model::Symbol=:chi_perp_e, chi_tor=1.0, @@ -273,13 +271,11 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; # 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) @@ -302,6 +298,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. `run_slayer` 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) @@ -392,7 +389,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), - D_perp=_eval(D_perp, 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, P_perp_model=P_perp_model, m=m_res, n=n_res, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index a174e8c11..d94de74d5 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -184,7 +184,8 @@ end qval, sval_r, bt, rs, R0, mu_i, zeff, chi_perp_e, chi_tor, m, n, - chi_perp_i=1.0, D_perp=1.0, tau_E=NaN, + chi_perp_i=1.0, tau_E=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, @@ -217,13 +218,16 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `zeff` -- effective charge - `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`) - - `D_perp` -- perpendicular particle diffusivity [m²/s] (used by `P_perp_model=: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 [s] (required by `P_perp_model=:tau_E`) - `P_perp_model` -- how `P_perp` is built: + `:chi_perp_e` (default) -- τ_R·χ⊥,e/r_s² + `: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_⊥ (Fitzpatrick, Phys. Plasmas 29, 032507 (2022)) + + `:D_perp` -- τ_R/τ_⊥ with τ_⊥ = r_s²/D_⊥ and D_⊥ from Eq. 16 of Fitzpatrick, Phys. Plasmas 29, + 032507 (2022), with η_⊥ = η_∥: + D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)[η_e τ_e/(1+η_e) χ⊥,e + η_i τ_i/(1+η_i) χ⊥,i] + `:c_beta` -- P_perp = C² with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) + `:tau_E` -- P_perp = P_tor = τ_R/τ_E, overriding χ_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131) `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode. @@ -281,7 +285,8 @@ function slayer_parameters(; rs::Real, R0::Real, mu_i::Real, zeff::Real, chi_perp_e::Real, chi_tor::Real, m::Integer, n::Integer, - chi_perp_i::Real=1.0, D_perp::Real=1.0, tau_E::Real=NaN, + chi_perp_i::Real=1.0, tau_E::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, @@ -400,7 +405,17 @@ function slayer_parameters(; # P_⊥ = P_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143) P_tor elseif P_perp_model === :D_perp - # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥ (Fitzpatrick, Phys. Plasmas 29, 032507 (2022)) + # D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥), with + # η_s/(1+η_s) = dlnT_s/(dlnn+dlnT_s) and τ = (T_e/T_i)(1+η_e)/(1+η_i) written flat-density safe. + frac_e = dlnTe / (dlnn + dlnTe) + frac_i = dlnTi / (dlnn + dlnTi) + tau_F = (t_e / t_i) * (dlnn + dlnTe) / (dlnn + dlnTi) + D_perp = c_beta^2 * eta / MU_0 + + (2 / 3) * (1 - c_beta^2) * (frac_e * tau_F / (1 + tau_F) * chi_perp_e + frac_i / (1 + tau_F) * chi_perp_i) + isfinite(D_perp) && D_perp > 0 || + throw(ArgumentError("slayer_parameters: P_perp_model=:D_perp gave D_perp=$D_perp; " * + "check dlnn/dlnTe/dlnTi (got $dlnn, $dlnTe, $dlnTi)")) + # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥ tau_r / (rs^2 / D_perp) elseif P_perp_model === :c_beta # C² = P_⊥ with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 61aa879c7..98943f1b1 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -38,12 +38,10 @@ constructor. 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 - - `D_perp` -- perpendicular particle diffusivity [m²/s] (default 1.0), used - only by `P_perp_model = :D_perp` - `tau_E` -- whole-plasma energy confinement time [s], required by `P_perp_model = :tau_E` (default `nothing`) - `P_perp_model` -- how the perpendicular Prandtl number is built: `:chi_perp_e` - (default), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:D_perp`, `:c_beta` + (default), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:D_perp` (D⊥ from χ⊥,e, χ⊥,i and the n, T gradients), `:c_beta` (P_perp = c_β²), or `:tau_E` (P_perp = P_tor = τ_R/τ_E); see `InnerLayer.SLAYER.slayer_parameters` - `dr_val`, `dgeo_val` -- critical-Δ formula inputs. `nothing` (default) @@ -118,7 +116,6 @@ there is one consistent interface for resistive and kinetic profiles. chi_perp_e::Float64 = 1.0 chi_tor::Float64 = 1.0 chi_perp_i::Float64 = 1.0 - D_perp::Float64 = 1.0 tau_E::Union{Float64,Nothing} = nothing P_perp_model::Symbol = :chi_perp_e dr_val::Union{Float64,Nothing} = nothing diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index c8fb91f3b..37cc6d1e0 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -462,7 +462,7 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac "(dataset absent or all-zero); using the scalar " * "control.chi_perp_e/chi_tor fallback for the missing one(s).") chi_perp_i = loaded.chi_perp_i === nothing ? control.chi_perp_i : loaded.chi_perp_i - (control.P_perp_model === :chi_perp_i && loaded.chi_perp_i === nothing) && @warn( + (control.P_perp_model in (:chi_perp_i, :D_perp) && loaded.chi_perp_i === nothing) && @warn( "SLAYER: kinetic file has no usable chi_i profile; using the scalar " * "control.chi_perp_i fallback.") params = build_slayer_inputs(equil, surfaces, profiles; @@ -472,7 +472,6 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac chi_perp_e=chi_perp_e, chi_tor=chi_tor, chi_perp_i=chi_perp_i, - D_perp=control.D_perp, tau_E=something(control.tau_E, NaN), P_perp_model=control.P_perp_model, dr_val=control.dr_val, diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index d34b05759..f3f8cd5df 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -206,7 +206,15 @@ @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 - @test slayer_parameters(; base..., P_perp_model=:D_perp, D_perp=3.0).P_perp ≈ 3 * ref.P_perp + # 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 + @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) te = slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1) diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index cf4bc5b87..f94e9dc8f 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -58,7 +58,6 @@ "mu_i" => 2.0, "dr_val" => 0.01, "P_perp_model" => "D_perp", - "D_perp" => 0.5, "scan_grid" => Dict{String,Any}( "Q_re_range" => [-5.0, 5.0], "Q_im_range" => [-1.0, 3.0], @@ -82,7 +81,6 @@ @test c.bt === 1.8 @test c.dr_val == 0.01 @test c.P_perp_model === :D_perp - @test c.D_perp == 0.5 @test c.Q_re_range == (-5.0, 5.0) @test c.Q_im_range == (-1.0, 3.0) @test c.nre == 50 From 27d07448c9aabbcdb061b0b7c1f6f162638ffa61 Mon Sep 17 00:00:00 2001 From: Evan Bursch Date: Mon, 5 Oct 2026 10:28:15 -0400 Subject: [PATCH 5/6] InnerLayer.SLAYER - DOCS - List chi_i and the P_perp_model keys in the DIII-D-like SLAYER example Co-Authored-By: Claude Opus 5.5 --- examples/DIIID-like_SLAYER_example/gpec.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index 3dedd663c..39d8130e3 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -74,14 +74,16 @@ dc_type = "none" # Critical-Δ offset selector: "none", "lar mu_i = 2.0 # Ion mass in proton-mass units (2.0 = deuterium) zeff = 1.0 # Effective charge 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 = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s]) 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] From 2083d02c8f66d7f22dfe22745ce2aa63731befd4 Mon Sep 17 00:00:00 2001 From: Evan Bursch Date: Mon, 5 Oct 2026 14:36:45 -0400 Subject: [PATCH 6/6] InnerLayer.SLAYER - DOCS - Cite Park 2022 as Phys. Plasmas 29, 072506 and describe c_beta as the no-transport limit of D_perp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The SLAYER file headers and the c_beta P_perp_model cited "Park et al., Phys. Plasmas 29, 122505 (2022)", which does not exist. The paper is J.-K. Park, Phys. Plasmas 29, 072506 (2022), doi:10.1063/5.0093079, as docs/src/citations.md already lists. Its Eq. 10 defines C² = c_β² + (1 − c_β²)K; C = c_β drops the thermal-conduction K, which makes P_perp = c_β² the χ⊥ → 0 limit of the D_perp model. Co-Authored-By: Claude Opus 5.5 --- src/InnerLayer/SLAYER/LayerParameters.jl | 8 +++++--- src/InnerLayer/SLAYER/LayerThickness.jl | 2 +- src/InnerLayer/SLAYER/Riccati.jl | 4 ++-- src/InnerLayer/SLAYER/SLAYER.jl | 4 ++-- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index d94de74d5..73b9e8c22 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). @@ -228,7 +228,9 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl + `:D_perp` -- τ_R/τ_⊥ with τ_⊥ = r_s²/D_⊥ and D_⊥ from Eq. 16 of Fitzpatrick, Phys. Plasmas 29, 032507 (2022), with η_⊥ = η_∥: D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)[η_e τ_e/(1+η_e) χ⊥,e + η_i τ_i/(1+η_i) χ⊥,i] - + `:c_beta` -- P_perp = C² with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) + + `:c_beta` -- P_perp = C² with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, with the + thermal-conduction term K dropped). This is the χ⊥ → 0 limit of `:D_perp` (P_⊥ = μ₀D_⊥/η = c_β²), so it is + a floor on P_perp; realistic anomalous transport makes D_⊥ much larger. + `:tau_E` -- P_perp = P_tor = τ_R/τ_E, overriding χ_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131) `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode. - `m`, `n` -- poloidal / toroidal mode numbers at the surface @@ -418,7 +420,7 @@ function slayer_parameters(; # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥ tau_r / (rs^2 / D_perp) elseif P_perp_model === :c_beta - # C² = P_⊥ with C = c_β (Park et al., Phys. Plasmas 29, 122505 (2022)) + # C² = P_⊥ with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, K = 0) c_beta^2 else throw(ArgumentError("slayer_parameters: unknown P_perp_model=$P_perp_model. " * diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl index 0b6db9761..431c3c806 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 033f71eca..9b4d0d9a9 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`).