Tearing - FEATURE - 🌱 Add a per-surface omega_E_kHz override of the kinetic-file E×B rotation

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

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.

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",)),

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

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)

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

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,