From 4f7f6657e207776f9a850b430eaf2bda4e0d76e7 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 29 Oct 2026 10:41:17 -0400 Subject: [PATCH 1/2] =?UTF-8?q?Tearing=20-=20BUGFIX=20-=20Doppler-shift=20?= =?UTF-8?q?each=20GGJ=20surface=20into=20its=20own=20E=C3=97B=20frame=20in?= =?UTF-8?q?=20the=20coupled=20growth=20rates?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit GGJ is a plasma-frame layer, so rotation enters only as a Doppler shift: surface k sees s + i·n·Ω_E,k. surface_coupling(::GGJModel, ...) takes omega_E and n, and one layer_argument helper applies each model's shift (real on SLAYER's Q, imaginary on GGJ's s) in the scalar, reduced and full determinants. Tests: a common rotation only shifts ω; a differential one moves the coupled γ, and only the relative rotation matters. Co-Authored-By: Claude Opus 5.5 --- src/Tearing/Dispersion/Coupled.jl | 7 ++--- src/Tearing/Dispersion/CoupledFullMatch.jl | 8 +++--- src/Tearing/Dispersion/GGJGrowthRates.jl | 9 ++++--- src/Tearing/Dispersion/SurfaceCoupling.jl | 31 +++++++++++++++++----- src/Tearing/Runner/run_slayer.jl | 2 +- test/runtests_dispersion_ggj.jl | 29 ++++++++++++++++++++ 6 files changed, 67 insertions(+), 19 deletions(-) diff --git a/src/Tearing/Dispersion/Coupled.jl b/src/Tearing/Dispersion/Coupled.jl index 43a4b7b87..38a8033e0 100644 --- a/src/Tearing/Dispersion/Coupled.jl +++ b/src/Tearing/Dispersion/Coupled.jl @@ -15,7 +15,7 @@ # det = mc(Q::ComplexF64) # # At each evaluation, for k = 1 .. msing_max, the inner-layer Δ is computed -# at a Q rescaled by `tauk_k / tauk_ref` and offset by that surface's real +# at a Q rescaled by `tauk_k / tauk_ref` and offset by that surface's # 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. @@ -38,7 +38,8 @@ 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). +`q_shift_k` is the E×B Doppler offset carried on each `SurfaceCoupling` (zero by default), +applied in the layer's own frequency variable (see `SurfaceCoupling`). """ struct MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}} surfaces::V @@ -96,7 +97,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 * (sc.tauk / ref_tauk) + sc.q_shift + Q_k = layer_argument(sc, Qc * (sc.tauk / ref_tauk)) # 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 diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index f1db6fd4f..72607f58d 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -89,9 +89,9 @@ end 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). 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`. +Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k` Doppler-shifted by its +E×B offset `q_shift_k`, exactly as in the reduced `multi_surface_coupling`: added to Q for a +SLAYER surface, and as `i·q_shift_k` to the growth rate s for a GGJ surface (see `SurfaceCoupling`). # Keyword arguments @@ -151,7 +151,7 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number) # 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 + Q_k = layer_argument(sc, Qc * (sc.tauk / ref_tauk)) resp = solve_inner(sc.model, sc.params, Q_k; mc.inner_kwargs...) # delta1 = interchange (parity +), delta2 = tearing (parity −); named diff --git a/src/Tearing/Dispersion/GGJGrowthRates.jl b/src/Tearing/Dispersion/GGJGrowthRates.jl index 8d3927865..bb4c284c3 100644 --- a/src/Tearing/Dispersion/GGJGrowthRates.jl +++ b/src/Tearing/Dispersion/GGJGrowthRates.jl @@ -53,10 +53,11 @@ its own 2 × 2 block of the raw Δ'). The box `gamma_range` × `omega_range` winding number; see `find_growth_rates`. The result is in the units of the SLAYER extraction: `gamma_Hz` is γ [s⁻¹], `omega_Hz` is −ω [rad/s], and `Q_root` is the normalized `(−ω + iγ)/Q₀`. The reported root is the least -stable zero in the box. Without drifts the GGJ roots come in ±ω pairs, so -either member of a pair can be reported. GGJ surfaces carry no time -normalization, so no E×B Doppler shift is applied and the roots are in the -plasma frame. +stable zero in the box. Each surface's layer is evaluated in its own E×B frame, +`s_k = s + i·n·Ω_E,k` from the `omega_E` and `n` given to `surface_coupling`; +the scanned s, and so ω, is then in the lab frame, and in the plasma frame when +no surface rotates. Without drifts or rotation the GGJ roots come in ±ω pairs, +so either member of a pair can be reported. `gamma_range` must lie at γ > 0: the GGJ layer is evaluated on the principal branch of a square root in s, with a cut on the negative real s-axis. The search diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index f243febec..95db72834 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -27,8 +27,10 @@ Per-surface dispersion data: `(model, params, dp_diag, dc, scale, tauk, q_shift) 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. +`q_shift` is the surface's E×B Doppler shift, a real frequency offset in the layer's own +frequency variable, zero by default. SLAYER's argument is the normalized Q, so the shift is added +to Q; GGJ's argument is the growth rate s = γ + iω, so the shift is added to ω and enters as +`s + i·q_shift`. The Tearing runner sets it for SLAYER surfaces. A root of `sc` in the complex `Q` plane is a **tearing** eigenvalue at this surface in the *uncoupled* approximation (only the tearing channel @@ -48,10 +50,15 @@ struct SurfaceCoupling{M<:InnerLayerModel,P} end function (sc::SurfaceCoupling)(Q::Number) - Δ = solve_inner(sc.model, sc.params, ComplexF64(Q) + sc.q_shift).tearing + Δ = solve_inner(sc.model, sc.params, layer_argument(sc, Q)).tearing return sc.dp_diag - sc.scale * Δ - sc.dc end +# The layer's own frequency argument for a scanned value already on this surface's normalization: +# a real Doppler shift of SLAYER's Q, an imaginary one of GGJ's s = γ + iω. +layer_argument(sc::SurfaceCoupling, Q::Number) = ComplexF64(Q) + sc.q_shift +layer_argument(sc::SurfaceCoupling{<:GGJModel}, s::Number) = ComplexF64(s) + im * sc.q_shift + """ surface_coupling(model::SLAYERModel, params::SLAYERParameters, dp_diag::Number; dc::Real=0.0, q_shift::Real=0.0) @@ -74,7 +81,7 @@ end """ surface_coupling(model::GGJModel, params::GGJParameters, - dp_diag::Number) -> SurfaceCoupling + dp_diag::Number; omega_E::Real=0.0, n::Integer=0) -> SurfaceCoupling GGJ convenience constructor. `scale` is `1.0` because GGJ's `solve_inner` applies its own `rescale_delta` (S^(2p₁/3)·v1^(2p₁)) internally, so the @@ -82,16 +89,26 @@ returned Δ is already in outer units. `tauk` defaults to `1.0` (GGJ has no direct analogue of SLAYER's per-surface time normalization, so multi-surface Q rescaling is a no-op for GGJ surfaces unless overridden). +`omega_E` is the surface's E×B rotation per unit toroidal mode number [rad/s] and `n` the +toroidal mode number. GGJ is a plasma-frame layer: rotation enters only beyond the order its +equations retain (Johnson & Greene 1967, Plasma Phys. 9, 611, Appendix B), so a rotating +surface sees the lab-frame growth rate Doppler-shifted, `s_k = s + i·n·omega_E`, and +`q_shift = n·omega_E`. With rotation set on every surface, the scanned `s` is in the lab +frame; differential rotation between coupled surfaces changes the coupled γ, while a common +rotation only shifts ω. + **No `dc` kwarg**: GGJ's 4m×4m Pletzer-Dewar residual already includes the 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, and so is `q_shift`. +struct's `dc` field is hard-wired to 0 here. """ -function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number) +function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number; omega_E::Real=0.0, n::Integer=0) + iszero(omega_E) || n >= 1 || + throw(ArgumentError("surface_coupling: E×B rotation omega_E=$(omega_E) rad/s needs the toroidal mode number n ≥ 1; got n=$(n)")) return SurfaceCoupling(model, params, ComplexF64(dp_diag), - 0.0, 1.0, 1.0, 0.0) + 0.0, 1.0, 1.0, Float64(n * omega_E)) end """ diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 26ba0c18c..835a6e492 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -238,7 +238,7 @@ function _q_shift(p::SLAYERParameters, Ω_E::Real) throw(ArgumentError("run_slayer: E×B rotation Ω_E=$(Ω_E) rad/s needs the toroidal mode number, but SLAYERParameters has n=$(p.n)")) return -p.tauk * p.n * Ω_E end -# GGJ carries no time normalization (tauk = 1), so a physical rotation has no Q-space image. +# The runner builds SLAYER surfaces only; a GGJ surface takes its rotation from `surface_coupling(...; omega_E, n)`. _q_shift(::InnerLayerParameters, ::Real) = 0.0 # --------------------------------------------------------------------- diff --git a/test/runtests_dispersion_ggj.jl b/test/runtests_dispersion_ggj.jl index d5bae4159..3dd5b8010 100644 --- a/test/runtests_dispersion_ggj.jl +++ b/test/runtests_dispersion_ggj.jl @@ -121,5 +121,34 @@ GeneralizedPerturbedEquilibrium.Dispersion._residual_domain_radius(f::GGJDomainT nre0=6, nim0=6, passes=3, parallel=false) @test length(gr2.valid_roots) == 1 && all(w -> w ≈ 1, gr2.valid_roots_winding) @test abs(gr2.Q_root - b) < 1e-8 * abs(b) + + @testset "E×B rotation Doppler-shifts each surface's s" begin + # Rotation enters the plasma-frame GGJ layer as s_k = s + i·n·Ω_E,k on each surface. + mc_rot(Ω1, Ω2) = multi_surface_coupling_full([surface_coupling(model, p, r0.tearing; omega_E=Ω1, n=1), + surface_coupling(model, p2, r02.tearing; omega_E=Ω2, n=1)], dp2) + Ω = 0.02 * Q0 + s = 0.11 * Q0 + 0.03im * Q0 + @test mc_rot(Ω, Ω)(s) == mc2(s + im * Ω) + @test_throws ArgumentError surface_coupling(model, p, r0.tearing; omega_E=Ω) + function secant_root(g, x0) + x1 = 1.01 * x0 + f0, f1 = g(x0), g(x1) + for _ in 1:60 + x0, f0, x1 = x1, f1, x1 - f1 * (x1 - x0) / (f1 - f0) + f1 = g(x1) + abs(x1 - x0) < 1e-14 && break + end + return x1 + end + # A common rotation is a change of frame: γ is unchanged and the lab ω shifts by n·Ω (Re Q by Ω/Q₀). + b_common = secant_root(GGJDispersion(mc_rot(Ω, Ω)), b + Ω / Q0) + @test abs(b_common - (b + Ω / Q0)) < 1e-8 * abs(b) + # Differential rotation moves the coupled growth rate; only the relative rotation matters. + Ω2 = 0.5 * Q0 + b_diff = secant_root(GGJDispersion(mc_rot(0.0, Ω2)), b) + @test abs(imag(b_diff) - imag(b)) > 1e-4 * imag(b) + b_frame = secant_root(GGJDispersion(mc_rot(-Ω2, 0.0)), b - Ω2 / Q0) + @test abs(b_frame - (b_diff - Ω2 / Q0)) < 1e-8 * abs(b) + end end end From c62678e11d1b0cc6cc1fa706340ccc359fd92517 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 29 Oct 2026 21:03:13 -0400 Subject: [PATCH 2/2] Tearing - MINOR - Drop the silent GGJ rotation fallback and consolidate the GGJ rotation tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Remove the unreachable _q_shift fallback, which would silently drop E×B rotation for a non-SLAYER surface; note in the ggj_growth_rates docstring that omega_range must cover the lab-frame ω when rotation is set; share one secant helper in the GGJ coupled tests and drop the assertions that restate the implementation or repeat the frame-change check. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/Tearing/Dispersion/GGJGrowthRates.jl | 3 +- src/Tearing/Dispersion/SurfaceCoupling.jl | 3 +- src/Tearing/Runner/run_slayer.jl | 2 -- test/runtests_dispersion_ggj.jl | 37 ++++++++--------------- 4 files changed, 16 insertions(+), 29 deletions(-) diff --git a/src/Tearing/Dispersion/GGJGrowthRates.jl b/src/Tearing/Dispersion/GGJGrowthRates.jl index bb4c284c3..59deb0179 100644 --- a/src/Tearing/Dispersion/GGJGrowthRates.jl +++ b/src/Tearing/Dispersion/GGJGrowthRates.jl @@ -56,7 +56,8 @@ SLAYER extraction: `gamma_Hz` is γ [s⁻¹], `omega_Hz` is −ω [rad/s], and stable zero in the box. Each surface's layer is evaluated in its own E×B frame, `s_k = s + i·n·Ω_E,k` from the `omega_E` and `n` given to `surface_coupling`; the scanned s, and so ω, is then in the lab frame, and in the plasma frame when -no surface rotates. Without drifts or rotation the GGJ roots come in ±ω pairs, +no surface rotates, so with rotation set `omega_range` must cover the lab-frame ω, +shifted by about n·Ω_E from the plasma frame. Without drifts or rotation the GGJ roots come in ±ω pairs, so either member of a pair can be reported. `gamma_range` must lie at γ > 0: the GGJ layer is evaluated on the principal diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index 95db72834..5e52b4183 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -54,8 +54,7 @@ function (sc::SurfaceCoupling)(Q::Number) return sc.dp_diag - sc.scale * Δ - sc.dc end -# The layer's own frequency argument for a scanned value already on this surface's normalization: -# a real Doppler shift of SLAYER's Q, an imaginary one of GGJ's s = γ + iω. +# Doppler-shift a scanned value into the layer's own frequency variable: real for SLAYER's Q, imaginary for GGJ's s. layer_argument(sc::SurfaceCoupling, Q::Number) = ComplexF64(Q) + sc.q_shift layer_argument(sc::SurfaceCoupling{<:GGJModel}, s::Number) = ComplexF64(s) + im * sc.q_shift diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 9b82db6a6..03bce8935 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -227,8 +227,6 @@ function _q_shift(p::SLAYERParameters, Ω_E::Real) throw(ArgumentError("run_slayer: E×B rotation Ω_E=$(Ω_E) rad/s needs the toroidal mode number, but SLAYERParameters has n=$(p.n)")) return -p.tauk * p.n * Ω_E end -# The runner builds SLAYER surfaces only; a GGJ surface takes its rotation from `surface_coupling(...; omega_E, n)`. -_q_shift(::InnerLayerParameters, ::Real) = 0.0 # --------------------------------------------------------------------- # Reference-length conversion of the outer Δ' for the slab layer diff --git a/test/runtests_dispersion_ggj.jl b/test/runtests_dispersion_ggj.jl index 3dd5b8010..04d40e7f8 100644 --- a/test/runtests_dispersion_ggj.jl +++ b/test/runtests_dispersion_ggj.jl @@ -109,13 +109,17 @@ GeneralizedPerturbedEquilibrium.Dispersion._residual_domain_radius(f::GGJDomainT mc2 = multi_surface_coupling_full([surface_coupling(model, p, r0.tearing), surface_coupling(model, p2, r02.tearing)], dp2) g2 = GGJDispersion(mc2) # Independent reference: a complex secant iteration on the determinant, started at the uncoupled root. - a, b = 0.1234im, 0.1234im * 1.01 - fa, fb = g2(a), g2(b) - for _ in 1:40 - a, fa, b = b, fb, b - fb * (b - a) / (fb - fa) - fb = g2(b) - abs(b - a) < 1e-14 && break + function secant_root(g, x0) + x1 = 1.01 * x0 + f0, f1 = g(x0), g(x1) + for _ in 1:60 + x0, f0, x1 = x1, f1, x1 - f1 * (x1 - x0) / (f1 - f0) + f1 = g(x1) + abs(x1 - x0) < 1e-14 && break + end + return x1 end + b = secant_root(g2, 0.1234im) @test abs(b - 0.1234im) > 1e-5 # the coupling moved the root gr2, _ = ggj_growth_rates(mc2; gamma_range=(0.08Q0, 0.17Q0), omega_range=(-0.031Q0, 0.047Q0), nre0=6, nim0=6, passes=3, parallel=false) @@ -127,28 +131,13 @@ GeneralizedPerturbedEquilibrium.Dispersion._residual_domain_radius(f::GGJDomainT mc_rot(Ω1, Ω2) = multi_surface_coupling_full([surface_coupling(model, p, r0.tearing; omega_E=Ω1, n=1), surface_coupling(model, p2, r02.tearing; omega_E=Ω2, n=1)], dp2) Ω = 0.02 * Q0 - s = 0.11 * Q0 + 0.03im * Q0 - @test mc_rot(Ω, Ω)(s) == mc2(s + im * Ω) @test_throws ArgumentError surface_coupling(model, p, r0.tearing; omega_E=Ω) - function secant_root(g, x0) - x1 = 1.01 * x0 - f0, f1 = g(x0), g(x1) - for _ in 1:60 - x0, f0, x1 = x1, f1, x1 - f1 * (x1 - x0) / (f1 - f0) - f1 = g(x1) - abs(x1 - x0) < 1e-14 && break - end - return x1 - end - # A common rotation is a change of frame: γ is unchanged and the lab ω shifts by n·Ω (Re Q by Ω/Q₀). + # A common rotation is a change of frame: γ is unchanged and the lab ω shifts by n·Ω (Re Q by Ω/Q₀, as in SLAYER). b_common = secant_root(GGJDispersion(mc_rot(Ω, Ω)), b + Ω / Q0) @test abs(b_common - (b + Ω / Q0)) < 1e-8 * abs(b) - # Differential rotation moves the coupled growth rate; only the relative rotation matters. - Ω2 = 0.5 * Q0 - b_diff = secant_root(GGJDispersion(mc_rot(0.0, Ω2)), b) + # Differential rotation moves the coupled growth rate. + b_diff = secant_root(GGJDispersion(mc_rot(0.0, 0.5 * Q0)), b) @test abs(imag(b_diff) - imag(b)) > 1e-4 * imag(b) - b_frame = secant_root(GGJDispersion(mc_rot(-Ω2, 0.0)), b - Ω2 / Q0) - @test abs(b_frame - (b_diff - Ω2 / Q0)) < 1e-8 * abs(b) end end end