Tearing - BUGFIX! - 🚨 Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation

regression-harness/cases/diiid_slayer_n1_coupled.toml

@@ -0,0 +1,97 @@
+# Regression case: DIII-D-like H-mode n=1 SLAYER coupled-determinant growth rate.
+# Reuses the n=1 SLAYER deck via [overrides] in coupled mode, the only harness coverage of the
+# multi-surface determinant: its inter-surface Q normalization and the per-surface E×B Doppler
+# shift taken from the kinetic file's omega_E. The Re(Q) window is widened because that shift
+# moves the coupled root to the lab frame, outside the deck's [-2, 2].
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "diiid_slayer_n1_coupled"
+description = "DIII-D-like H-mode equilibrium, n=1, SLAYER coupled multi-surface determinant with kinetic-file E×B rotation (reuses the n=1 deck via overrides)"
+example_dir = "examples/DIIID-like_SLAYER_example"
+
+# Run the shared SLAYER deck through the coupled determinant.
+[overrides]
+"SLAYER.coupling_mode" = "coupled"
+"SLAYER.Q_re_range" = [-4.0, 6.0]
+
+# Per-surface inputs to the determinant
+[quantities.slayer_tauk]
+h5path = "Tearing/PerSurface/tau_k"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER tauk"
+noise_threshold = 1e-12
+order = 10
+class = "physics_converged"
+
+# Applied E×B Doppler offset ΔRe(Q) = −τ_k·n·Ω_E per surface
+[quantities.slayer_q_shift]
+h5path = "Tearing/PerSurface/q_shift"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER q_shift"
+noise_threshold = 1e-10
+order = 11
+class = "physics_converged"
+
+# Coupled eigenvalue (one root for the whole determinant). The coupled path has no root
+# polishing, so the thresholds match the uncoupled case's absolute, modestly loose floors.
+# Pin Q, ω and γ with a measured tolerance, not the class default: re-extracting one set of scan
+# samples under different triangulations gave γ = 190.646 and 191.888 s⁻¹ (0.65 %). ω and Q
+# need the same measurement.
+[quantities.slayer_Q_coupled]
+h5path = "Tearing/Roots/Q_root"
+type = "complex_vector"
+extract = "all_complex"
+label = "SLAYER coupled Q_root"
+noise_threshold = 1e-4
+order = 30
+class = "physics_converged"
+
+[quantities.slayer_omega_coupled]
+h5path = "Tearing/Roots/omega"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER coupled ω"
+noise_threshold = 1.0
+order = 32
+class = "physics_converged"
+
+[quantities.slayer_gamma_coupled]
+h5path = "Tearing/Roots/gamma"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER coupled γ"
+noise_threshold = 1e-1
+order = 33
+class = "physics_converged"
+
+# no_root flag (1 = extraction failed); guards that the coupled root stays inside the window.
+[quantities.slayer_no_root_coupled]
+h5path = "Tearing/Roots/no_root"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER coupled no_root flag"
+noise_threshold = 0
+order = 34
+class = "physics_converged"
+
+# Settings (catches accidental config drift)
+[quantities.slayer_enabled]
+h5path = "Tearing/enabled"
+type = "int_scalar"
+extract = "value"
+label = "SLAYER enabled flag"
+noise_threshold = 0
+order = 90
+class = "topological"
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
+class = "diagnostic"

src/InnerLayer/SLAYER/LayerInputs.jl

@@ -238,7 +238,6 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     dgeo_val=nothing,
     dc_type::Symbol=:none,
     theta::Real=0.0,
-    compute_omega_star::Bool=true,
     resistivity_model::NeoResistivityModel=SauterNeoModel(),
     lnLambda_form::Symbol=:nrl)
     R0_use = R0 === nothing ? equil.ro : Float64(R0)
@@ -291,9 +290,8 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
         n_res = sing.n[1]
 
         prof = profiles(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)
+        # ω_*e, ω_*i from the density and temperature spline derivatives, at this surface's n.
+        ω_e_use, ω_i_use = _omega_star_at(psi, n_res)
 
         # Pull geometric trapped-fraction inputs from ResistGeometry when
         # available (populated by ForceFreeStates.resist_eval_all!); else
@@ -375,7 +373,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
 
         out[k] = slayer_parameters(;
             n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i,
-            omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use,
+            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),

src/InnerLayer/SLAYER/LayerParameters.jl

@@ -177,7 +177,7 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
 end
 
 """
-    slayer_parameters(; n_e, t_e, t_i, omega, omega_e, omega_i,
+    slayer_parameters(; n_e, t_e, t_i, omega_e, omega_i,
                         qval, sval_r, bt, rs, R0, mu_i, zeff,
                         chi_perp, chi_tor,
                         m, n,
@@ -199,7 +199,6 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl
   - `n_e` -- electron density [m⁻³]
   - `t_e` -- electron temperature [eV]
   - `t_i` -- ion temperature [eV]
-  - `omega`   -- toroidal rotation frequency at the surface [rad/s]
   - `omega_e` -- electron diamagnetic frequency [rad/s]
   - `omega_i` -- ion diamagnetic frequency [rad/s]
   - `qval`    -- safety factor q at the surface
@@ -260,7 +259,7 @@ dispersion relation.
 """
 function slayer_parameters(;
     n_e::Real, t_e::Real, t_i::Real,
-    omega::Real, omega_e::Real, omega_i::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,

src/Tearing/Dispersion/Coupled.jl

@@ -15,10 +15,10 @@
 #   det = mc(Q::ComplexF64)
 #
 # At each evaluation, for k = 1 .. msing_max, the inner-layer Δ is computed
-# at a Q rescaled by `tauk_ref / tauk_k`, then subtracted (with the dc
-# offset) from the diagonal of an `msing_max × msing_max` upper-left
-# submatrix of `dp_matrix`. The off-diagonal Δ' couplings are passed
-# through unchanged.
+# at a Q rescaled by `tauk_k / tauk_ref` and offset by that surface's real
+# Doppler shift `q_shift`, then subtracted (with the dc offset) from the
+# diagonal of an `msing_max × msing_max` upper-left submatrix of `dp_matrix`.
+# The off-diagonal Δ' couplings are passed through unchanged.
 
 """
     MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}}
@@ -29,11 +29,16 @@ normalization), and the truncation `msing_max` (number of surfaces actually
 participating in the determinant). Calling `mc(Q)` returns `det(M(Q))` where
 
 ```
-M[k,k] = dp_matrix[k,k] - scale_k · Δ_inner_k(Q · tauk_ref / tauk_k) - dc_k
+M[k,k] = dp_matrix[k,k] - scale_k · Δ_inner_k(Q · ratio_k + q_shift_k) - dc_k
 M[i,j] = dp_matrix[i,j]      for i ≠ j        (off-diagonal Δ' couplings)
 ```
 
 A root of `mc` in the complex `Q` plane is a coupled tearing eigenvalue.
+
+`ratio_k = tauk_k/tauk_ref`: `Q` is defined as `tauk·ω`, so one shared physical
+ω reaches surface `k` as `Q_k = tauk_k·ω`, on the same scaling as its `Q_e`/`Q_i`.
+
+`q_shift_k` is the real offset carried on each `SurfaceCoupling` (zero by default).
 """
 struct MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}}
     surfaces::V
@@ -91,7 +96,7 @@ function (mc::MultiSurfaceCoupling)(Q::Number)
     M = mc.dp_matrix[1:n, 1:n]
     @inbounds for k in 1:n
         sc = mc.surfaces[k]
-        Q_k = Qc * (ref_tauk / sc.tauk)
+        Q_k = Qc * (sc.tauk / ref_tauk) + sc.q_shift
         # m×m scalar coupling: use only the tearing channel. The
         # interchange (Glasser-stabilization) channel is carried in the
         # full 4m×4m dispersion in `CoupledFullMatch.jl`; this reduced

src/Tearing/Dispersion/CoupledFullMatch.jl

@@ -71,33 +71,27 @@ studies use the reduced m × m `MultiSurfaceCoupling` instead.
   - `dp_raw::Matrix{ComplexF64}` — 2m × 2m outer-region matrix (side-major).
   - `ref_idx::Int`              — reference surface for Q rescaling (1-based).
   - `msing_max::Int`            — number of surfaces to include (truncates).
-  - `rotation::Vector{Float64}` — per-surface rotation frequencies (s⁻¹).
-  - `ntor::Int`                 — toroidal mode number `n` (default 1).
 """
 struct MultiSurfaceCouplingFull{V<:AbstractVector{<:SurfaceCoupling},K<:NamedTuple}
     surfaces::V
     dp_raw::Matrix{ComplexF64}
     ref_idx::Int
     msing_max::Int
-    rotation::Vector{Float64}
-    ntor::Int
     inner_kwargs::K    # kwargs forwarded to solve_inner; e.g. (pfac=0.1, nx=128, nq=5)
 end
 
 """
     multi_surface_coupling_full(surfaces, dp_raw;
                                 ref_idx=1,
                                 msing_max=length(surfaces),
-                                rotation=zeros(length(surfaces)),
-                                ntor=1) -> MultiSurfaceCouplingFull
+                                inner_kwargs=NamedTuple()) -> MultiSurfaceCouplingFull
 
 Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the
 2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw`
 field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the
-Riccati path). `rotation[k]` is the per-surface rotation
-frequency; it shifts the per-surface inner Q argument by
-`i·ntor·rotation[k]`. Default zero rotation matches the static-equilibrium
-case.
+Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k + q_shift_k`, with
+`ratio_k` and the real E×B Doppler offset `q_shift_k` exactly as in the reduced
+`multi_surface_coupling`.
 
 # Keyword arguments
 
@@ -107,9 +101,6 @@ case.
     matching matrix becomes 4·msing_max × 4·msing_max, built from the
     corresponding 2·msing_max × 2·msing_max submatrix of `dp_raw`.
     Defaults to `length(surfaces)`.
-  - `rotation`  — per-surface rotation frequencies in s⁻¹ (length m).
-    Defaults to all zero.
-  - `ntor`      — toroidal mode number n. Defaults to 1.
   - `inner_kwargs` — NamedTuple of kwargs forwarded to `solve_inner` at
     every Q evaluation, e.g. `(pfac=0.1, xfac=10.0, nx=128, nq=5)` for
     Galerkin grid tuning. Defaults to `NamedTuple()`.
@@ -118,8 +109,6 @@ function multi_surface_coupling_full(surfaces::AbstractVector{<:SurfaceCoupling}
     dp_raw::AbstractMatrix;
     ref_idx::Integer=1,
     msing_max::Integer=length(surfaces),
-    rotation::AbstractVector{<:Real}=zeros(length(surfaces)),
-    ntor::Integer=1,
     inner_kwargs::NamedTuple=NamedTuple())
     m = length(surfaces)
     size(dp_raw) == (2m, 2m) ||
@@ -131,14 +120,9 @@ function multi_surface_coupling_full(surfaces::AbstractVector{<:SurfaceCoupling}
     1 <= msing_max <= m ||
         throw(ArgumentError("multi_surface_coupling_full: msing_max=$msing_max " *
                             "out of range 1:$m"))
-    length(rotation) == m ||
-        throw(ArgumentError("multi_surface_coupling_full: rotation length " *
-                            "$(length(rotation)) ≠ $m"))
     return MultiSurfaceCouplingFull(surfaces,
         Matrix{ComplexF64}(dp_raw),
         Int(ref_idx), Int(msing_max),
-        Float64.(collect(rotation)),
-        Int(ntor),
         inner_kwargs)
 end
 
@@ -165,10 +149,9 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number)
         idx3 = idx1 + s2       # d^k_+
         idx4 = idx2 + s2       # d^k_-
 
-        # Per-surface Q shift: guess_modify = Q + i·n·rotation[k].
-        # Also apply ref_tauk / sc.tauk rescaling (we keep the SurfaceCoupling
-        # tauk normalization that SLAYER needs; GGJ has tauk=1 so it's a no-op).
-        Q_k = Qc * (ref_tauk / sc.tauk) + 1im * mc.ntor * mc.rotation[k]
+        # Map the shared scanned Q onto this surface's normalization, then Doppler it into
+        # the surface's E×B frame (GGJ carries tauk = 1, so the ratio is a no-op there).
+        Q_k = Qc * (sc.tauk / ref_tauk) + sc.q_shift
         resp = solve_inner(sc.model, sc.params, Q_k; mc.inner_kwargs...)
 
         # delta1 = interchange (parity −), delta2 = tearing (parity +); named

src/Tearing/Dispersion/GrowthRateExtraction.jl

@@ -22,6 +22,7 @@
 
 using Contour
 using DelaunayTriangulation
+using Random: Xoshiro
 
 # ---------------------------------------------------------------------
 # Public result struct + main entry point.
@@ -955,7 +956,9 @@ function _extract_growth_rates_amr(Q::Vector{ComplexF64},
         throw(ArgumentError("_extract_growth_rates_amr: need ≥ 3 points to triangulate"))
 
     pts = [(real(q), imag(q)) for q in Q]
-    tri = triangulate(pts)
+    # AMR grids are full of cocircular points, where the Delaunay triangulation is not unique and
+    # the randomized insertion order picks the diagonals; a fixed seed makes the contours reproducible.
+    tri = triangulate(pts; rng=Xoshiro(0))
 
     # Segment types (carry complementary-field value at each endpoint)
     re_segs = NamedTuple{(:p1, :p2, :a1, :a2),

src/Tearing/Dispersion/SurfaceCoupling.jl

@@ -21,13 +21,15 @@
 """
     SurfaceCoupling{M<:InnerLayerModel, P}
 
-Per-surface dispersion data: `(model, params, dp_diag, dc, scale, tauk)`.
-Calling `sc(Q)` returns the complex residual
+Per-surface dispersion data: `(model, params, dp_diag, dc, scale, tauk, q_shift)`. Calling `sc(Q)` returns the complex residual
 
 ```
-r(Q) = dp_diag - scale * solve_inner(model, params, Q).tearing - dc
+r(Q) = dp_diag - scale * solve_inner(model, params, Q + q_shift).tearing - dc
 ```
 
+`q_shift` is a real offset on the layer's Q argument, zero by default. The Tearing runner
+sets it to the surface's E×B Doppler shift.
+
 A root of `sc` in the complex `Q` plane is a **tearing** eigenvalue at
 this surface in the *uncoupled* approximation (only the tearing channel
 of the inner-layer response appears — the interchange channel enters the
@@ -42,16 +44,18 @@ struct SurfaceCoupling{M<:InnerLayerModel,P}
     dc::Float64
     scale::Float64
     tauk::Float64
+    q_shift::Float64
 end
 
 function (sc::SurfaceCoupling)(Q::Number)
-    Δ = solve_inner(sc.model, sc.params, ComplexF64(Q)).tearing
+    Δ = solve_inner(sc.model, sc.params, ComplexF64(Q) + sc.q_shift).tearing
     return sc.dp_diag - sc.scale * Δ - sc.dc
 end
 
 """
     surface_coupling(model::SLAYERModel, params::SLAYERParameters,
-                     dp_diag::Number; dc::Real=0.0) -> SurfaceCoupling
+                     dp_diag::Number; dc::Real=0.0, q_shift::Real=0.0)
+        -> SurfaceCoupling
 
 SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)`, which
 maps the dimensionless inner-layer Δ from `riccati_f` to the r_s-referenced
@@ -63,9 +67,9 @@ before building couplings. `tauk` is taken from `params.tauk` for use by
 `MultiSurfaceCoupling` Q rescaling.
 """
 function surface_coupling(model::SLAYERModel, params::SLAYERParameters,
-    dp_diag::Number; dc::Real=0.0)
+    dp_diag::Number; dc::Real=0.0, q_shift::Real=0.0)
     return SurfaceCoupling(model, params, ComplexF64(dp_diag),
-        Float64(dc), params.lu^(1 / 3), params.tauk)
+        Float64(dc), params.lu^(1 / 3), params.tauk, Float64(q_shift))
 end
 
 """
@@ -83,25 +87,24 @@ interchange channel, which provides Glasser (Mercier) stabilization
 natively. A Δ_crit proxy (χ_parallel-matching offset on the diagonal) is
 meaningful only for tearing-only slab-layer approximations like SLAYER;
 for GGJ it would double-count the interchange physics. The `SurfaceCoupling`
-struct's `dc` field is hard-wired to 0 here.
+struct's `dc` field is hard-wired to 0 here, and so is `q_shift`.
 """
-function surface_coupling(model::GGJModel, params::GGJParameters,
-    dp_diag::Number)
+function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number)
     return SurfaceCoupling(model, params, ComplexF64(dp_diag),
-        0.0, 1.0, 1.0)
+        0.0, 1.0, 1.0, 0.0)
 end
 
 """
     surface_coupling(model::InnerLayerModel, params, dp_diag::Number;
-                     dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0)
-        -> SurfaceCoupling
+                     dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0,
+                     q_shift::Real=0.0) -> SurfaceCoupling
 
 Generic fallback constructor. Use this when wiring a new inner-layer model
 into the dispersion solver — pass the appropriate inner→outer-units `scale`
 and per-surface `tauk` explicitly.
 """
 function surface_coupling(model::InnerLayerModel, params, dp_diag::Number;
-    dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0)
+    dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0, q_shift::Real=0.0)
     return SurfaceCoupling(model, params, ComplexF64(dp_diag),
-        Float64(dc), Float64(scale), Float64(tauk))
+        Float64(dc), Float64(scale), Float64(tauk), Float64(q_shift))
 end