diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a8a99f4c7..29b91a47f 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -75,6 +75,7 @@ mu_i = 2.0 # Ion mass in proton-mass units (2.0 = deut zeff = 1.0 # Effective charge chi_perp = 1.0 # fallback only; the kinetic file supplies χ⊥(ψ) chi_tor = 1.0 # fallback only; the kinetic file supplies χ_φ(ψ) +# omega_E_kHz = {"2/1" = 0.0, "3/1" = 3.0} # override of the E×B rotation Ω_E/2π per unit n [kHz], keyed by m/n; unlisted surfaces use the kinetic file 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 diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 1a12ff7e6..c42ec2903 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -32,6 +32,11 @@ constructor. - `bt` -- toroidal field `[T]`. `nothing` (default) resolves the physical `B_T = F(ψ)/(2π·R₀)` per surface from the equilibrium's F-spline; a scalar or a callable of `psi` overrides it + - `omega_E_kHz` -- optional override of the per-surface E×B rotation, keyed by the surface's + `"m/n"` (TOML: `omega_E_kHz = {"2/1" = 0.0, "3/1" = 3.0}`). Per unit toroidal mode number + like the kinetic file's `omega_E`, but given as Ω_E/2π in kHz rather than the file's rad/s. + Surfaces not listed (all of them by default) take Ω_E from the kinetic file; a key matching + no analysed surface is an error. Rotation enters only the coupled determinant - `mu_i` -- ion mass in proton-mass units (default 2.0 for D) - `zeff` -- effective charge - `chi_perp`, `chi_tor` -- fallback perpendicular / toroidal heat @@ -151,6 +156,8 @@ there is one consistent interface for resistive and kinetic profiles. # / failed-Δ'-BVP surface, not a real root. Flagged `:spurious`. validity_rtol::Float64 = 1e-3 + omega_E_kHz::Dict{String,Float64} = Dict{String,Float64}() + profile_file::String = "" profile_group::String = "/" @@ -185,6 +192,9 @@ function validate(ctrl::SLAYERControl) "not in $(_VALID_LNLAMBDA_FORMS)")) ctrl.msing_max >= 1 || throw(ArgumentError("SLAYERControl: msing_max=$(ctrl.msing_max) must be ≥ 1")) + all(isfinite, values(ctrl.omega_E_kHz)) || + throw(ArgumentError("SLAYERControl: omega_E_kHz contains a " * + "non-finite entry: $(ctrl.omega_E_kHz)")) ctrl.nre >= 2 && ctrl.nim >= 2 || throw(ArgumentError("SLAYERControl: nre and nim must both be ≥ 2")) ctrl.amr_passes >= 0 || @@ -249,6 +259,11 @@ function slayer_control_from_toml(section::AbstractDict) elseif sym in (:bt, :dr_val, :dgeo_val) # Allow explicit nothing (auto-derive) or a number (override) kwargs[sym] = v === nothing ? nothing : Float64(v) + elseif sym === :omega_E_kHz + v isa AbstractDict || + throw(ArgumentError("slayer_control_from_toml: omega_E_kHz must be a table keyed by m/n, " * + "e.g. omega_E_kHz = {\"2/1\" = 0.0, \"3/1\" = 3.0}; got $v")) + kwargs[sym] = Dict{String,Float64}(String(key) => Float64(x) for (key, x) in v) elseif sym === :boxes # `boxes` is a Vector{NTuple{4,Float64}}; from TOML this comes # in as a list of 4-element arrays. Coerce each. diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index b5395bfde..1dc6d9ee7 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -86,7 +86,7 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/tau_k" => (; long_name="Q-normalization time S^(1/3)·τ_H per surface (Q = τ_k·ω; diamagnetic inputs Q_e, Q_i carry the opposite sign by convention)", units="s", dims=("surface",)), "PerSurface/omega_E" => - (; long_name="E×B angular frequency Ω_E per unit n on each surface, from the kinetic file", units="rad/s", dims=("surface",)), + (; long_name="E×B angular frequency Ω_E per unit n on each surface (kinetic file or omega_E_kHz override)", units="rad/s", dims=("surface",)), "PerSurface/q_shift" => (; long_name="real E×B Doppler offset ΔRe(Q) = −τ_k·n·Ω_E applied to each surface's inner-layer Q in the coupled determinant (0 in uncoupled mode)", dims=("surface",)), "PerSurface/tau_R" => (; long_name="resistive diffusion time τ_R = μ₀r_s²/η per surface", units="s", dims=("surface",)), diff --git a/src/Tearing/Runner/Result.jl b/src/Tearing/Runner/Result.jl index 3a018a38d..1d2ddfb0b 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -38,8 +38,9 @@ downstream inspection and HDF5 output. - `layer_widths` -- `Vector{LayerWidths}`, one per surface: the resistive layer thickness (in meters) from the `del_s` Riccati solve plus FKR / visco-resistive sanity scales. Empty when disabled. - - `omega_E` -- E×B angular frequency Ω_E per unit n [rad/s] on each surface, - from the kinetic file. Recorded in both coupling modes + - `omega_E` -- E×B angular frequency Ω_E per unit n [rad/s] resolved on each + surface: the kinetic file's value, or the `control.omega_E_kHz` override for that m/n. + Recorded in both coupling modes - `q_shift` -- real Doppler offset `−τ_k·n·Ω_E` applied to each surface's inner-layer Q; zero outside `:coupled` mode, where rotation does not enter - `scan_data` -- scan results (per-surface in uncoupled, single diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 0556feeb2..ee4b89071 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -147,6 +147,32 @@ end _build_surface_coupling(model::GGJModel, params::GGJParameters, dp_diag, ::Real) = surface_coupling(model, params, dp_diag) +# "m/n" label keying `control.omega_E_kHz`; GGJ parameters carry no mode numbers. +_mn_label(p::SLAYERParameters) = "$(p.m)/$(p.n)" +_mn_label(::InnerLayerParameters) = nothing + +# Per-surface E×B angular frequency Ω_E per unit n [rad/s]: the kinetic file's `omega_E` at each +# rational surface, replaced on every surface whose m/n `control.omega_E_kHz` lists. +function _omega_E_per_surface(control::SLAYERControl, params::AbstractVector, omega_E_file::AbstractVector{<:Real}) + n = length(params) + isempty(omega_E_file) || length(omega_E_file) == n || + throw(ArgumentError("run_slayer: omega_E has $(length(omega_E_file)) entries but " * + "$n rational surfaces were analysed.")) + Ω_E = isempty(omega_E_file) ? zeros(Float64, n) : Float64.(omega_E_file) + isempty(control.omega_E_kHz) && return Ω_E + labels = [_mn_label(p) for p in params] + any(isnothing, labels) && + throw(ArgumentError("run_slayer: omega_E_kHz is keyed by m/n, which $(eltype(params)) surfaces do not carry.")) + unknown = setdiff(keys(control.omega_E_kHz), labels) + isempty(unknown) || + throw(ArgumentError("run_slayer: omega_E_kHz lists m/n $(sort(collect(unknown))) matching no analysed surface; " * + "the analysed surfaces are $(labels).")) + for (k, label) in enumerate(labels) + haskey(control.omega_E_kHz, label) && (Ω_E[k] = 2π * 1e3 * control.omega_E_kHz[label]) + end + return Ω_E +end + """ _q_shift(p::SLAYERParameters, Ω_E) -> Float64 @@ -225,8 +251,9 @@ Run the SLAYER tearing analysis given pre-built per-surface equilibrium-driven `build_slayer_inputs` step — use this when the parameters are already known (e.g. in unit tests or when rebuilding from cached HDF5 output). `omega_E` is the E×B angular frequency per unit n -at each surface [rad/s]; empty means no rotation. In coupled mode it Doppler-shifts -each surface's inner-layer Q; in both modes it is reported as `result.omega_E`. +at each surface [rad/s], replaced on any surface whose m/n `control.omega_E_kHz` +lists; empty means no rotation. In coupled mode it Doppler-shifts each surface's +inner-layer Q; in both modes the resolved value is reported as `result.omega_E`. """ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, dp_matrix::AbstractMatrix, @@ -297,9 +324,11 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, # E×B rotation only Doppler-shifts the coupled determinant: each uncoupled layer is already # solved in its own plasma frame, so a shift there would only relabel the reported frequency. coupled = control.coupling_mode === :coupled - isempty(omega_E) || length(omega_E) == n || - throw(ArgumentError("run_slayer: omega_E has $(length(omega_E)) entries but $n rational surfaces were analysed.")) - Ω_E = isempty(omega_E) ? zeros(Float64, n) : omega_E + !coupled && !isempty(control.omega_E_kHz) && + @warn( + "SLAYER: omega_E_kHz does not shift the layers with coupling_mode=:uncoupled; rotation only " * + "enters the coupled determinant, so it is just recorded in PerSurface/omega_E.") + Ω_E = _omega_E_per_surface(control, params, omega_E) coupled && _is_ggj(model) && any(!iszero, Ω_E) && @warn( "SLAYER: E×B rotation is not applied to GGJ surfaces, which carry no time normalization.") q_shifts = coupled ? Float64[_q_shift(params[k], Ω_E[k]) for k in 1:n] : zeros(Float64, n) diff --git a/test/runtests_dispersion_rotation.jl b/test/runtests_dispersion_rotation.jl index 488fa95b8..61e713b15 100644 --- a/test/runtests_dispersion_rotation.jl +++ b/test/runtests_dispersion_rotation.jl @@ -2,7 +2,8 @@ using GeneralizedPerturbedEquilibrium.InnerLayer using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel, solve_inner using GeneralizedPerturbedEquilibrium.Dispersion - using GeneralizedPerturbedEquilibrium.Tearing.Runner: _q_shift, _build_surface_coupling + using GeneralizedPerturbedEquilibrium.Tearing.Runner: SLAYERControl, + slayer_control_from_toml, validate, _q_shift, _build_surface_coupling using LinearAlgebra # Linear inner layer Δ(Q) = a + b·Q makes the applied Q offset readable @@ -84,4 +85,15 @@ @test sort(roots ./ tauk[ref]; by=real) ≈ expected rtol = 1e-10 end end + + @testset "omega_E_kHz parses and validates from TOML" begin + ctrl = slayer_control_from_toml(Dict("omega_E_kHz" => Dict("2/1" => 0, "3/1" => 3.0))) + @test ctrl.omega_E_kHz == Dict("2/1" => 0.0, "3/1" => 3.0) + # Default stays empty so existing decks are untouched. + @test isempty(slayer_control_from_toml(Dict{String,Any}()).omega_E_kHz) + @test_throws "table keyed by m/n" slayer_control_from_toml(Dict("omega_E_kHz" => [0.0, 3.0])) + # A non-finite shift would silently poison every Q evaluation; the + # validator rejects it (TOML cannot express NaN, so go through validate). + @test_throws ArgumentError validate(SLAYERControl(; omega_E_kHz=Dict("2/1" => NaN))) + end end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index de18d5623..30a52c89d 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -161,6 +161,14 @@ @test_throws "rational surfaces were analysed" run_slayer_from_inputs(params, dpm, SLAYERControl(; coupling_mode=:coupled, grid...); omega_E=[1.0]) + # omega_E_kHz takes precedence over the file on the surfaces it lists; the rest keep the file value. + r_ovr = run_slayer_from_inputs(params, dpm, + SLAYERControl(; coupling_mode=:coupled, omega_E_kHz=Dict("2/1" => 1.0), grid...); omega_E=Ω_E) + @test r_ovr.omega_E ≈ [2π * 1e3, Ω_E[2]] + @test r_ovr.q_shift ≈ [-1.0e-4 * 2π * 1e3, -1.2e-4 * Ω_E[2]] + @test_throws "matching no analysed surface" run_slayer_from_inputs(params, dpm, + SLAYERControl(; coupling_mode=:coupled, omega_E_kHz=Dict("5/1" => 1.0), grid...); omega_E=Ω_E) + # A Doppler offset beyond the scan box is flagged rather than silently losing the root: # τ_ref·n·Ω_E = 1e-4 · 3e4 = 3 lies outside Re(Q) ∈ [-1, 1]. @test_logs (:warn, r"outside Q_re_range") match_mode = :any run_slayer_from_inputs(params, dpm,