Tearing - FEATURE - Extract GGJ growth rates from the full matching determinant

src/Tearing/Dispersion/Dispersion.jl

@@ -11,6 +11,8 @@
 #   - `brute_force_scan`    -- regular 2D Q-plane scan
 #   - `find_growth_rates`   -- contour-intersection root extraction (Re=0 ∩ Im=0)
 #   - `amr_scan`            -- adaptive Q-plane refinement
+#   - `ggj_growth_rates`    -- GGJ roots of the full 4m×4m determinant, classified
+#                              by winding number
 #
 # Roots are ISOLATED by 2D contour intersection on Nyquist-style Q-plane scans
 # (`find_growth_rates`, Re=0 ∩ Im=0), then optionally POLISHED to the true zero
@@ -43,6 +45,7 @@ include("CoupledFullMatch.jl")
 include("BruteForceScan.jl")
 include("ContourSearchAMR.jl")
 include("GrowthRateExtraction.jl")
+include("GGJGrowthRates.jl")
 
 export SurfaceCoupling, surface_coupling
 export MultiSurfaceCoupling, multi_surface_coupling
@@ -51,5 +54,6 @@ export ScanResult, brute_force_scan
 export AMRCell, AMRResult, amr_scan
 export BoxActivity, MultiBoxAMRResult, multi_box_amr_scan, as_amr_result
 export GrowthRateResult, find_growth_rates
+export GGJDispersion, ggj_growth_rates
 
 end # module Dispersion

src/Tearing/Dispersion/GGJGrowthRates.jl

@@ -0,0 +1,100 @@
+# GGJGrowthRates.jl
+#
+# Tearing growth rates from the GGJ inner layer. The dispersion function is the
+# full Pletzer-Dewar 4m × 4m determinant (CoupledFullMatch.jl), which keeps the
+# interchange channel and with it the Glasser stabilization. It is scanned in the
+# same normalized Q plane as the SLAYER residual and its roots are classified
+# by winding number (`find_growth_rates(...; classification=:winding)`).
+#
+# Plane convention. The determinant takes the complex growth rate
+# s = γ + iω [s⁻¹]. The SLAYER scans use Q = τ·(−ω + iγ), so s = −i·Q/τ; with
+# τ = 1/Q₀ of the reference surface (Q₀ = S^(−1/3)/τ_A, GWP2016 Eq. A16) the
+# scanned variable is the GGJ scaled growth rate rotated by a quarter turn,
+# Im Q = γ/Q₀ and Re Q = −ω/Q₀, and `GrowthRateResult.gamma_Hz` is γ in s⁻¹.
+# s ↦ Q is holomorphic, so a zero of the determinant winds +1 (the SLAYER
+# residual, a function of conj(Q), winds −1).
+
+"""
+    GGJDispersion(mc::MultiSurfaceCouplingFull)
+
+The GGJ full-matching determinant as a function of the normalized frequency
+`Q = (−ω + iγ)/Q₀` used by the SLAYER scans, `g(Q) = mc(−i·Q·Q₀)`, where `Q₀`
+is the growth-rate scale `S^(−1/3)/τ_A` of the reference surface
+`mc.surfaces[mc.ref_idx]`. Every surface must use a `GGJModel`.
+"""
+struct GGJDispersion{M<:MultiSurfaceCouplingFull}
+    mc::M
+    Q0::Float64
+end
+
+function GGJDispersion(mc::MultiSurfaceCouplingFull)
+    all(sc -> sc.model isa GGJModel, mc.surfaces) ||
+        throw(ArgumentError("GGJDispersion: every surface must use a GGJModel"))
+    return GGJDispersion(mc, InnerLayer.GGJ.q0(mc.surfaces[mc.ref_idx].params))
+end
+
+(g::GGJDispersion)(Q::Number) = g.mc(-im * ComplexF64(Q) * g.Q0)
+
+# Keep classifying circles, polishing steps and derivative probes strictly in γ > 0.
+_residual_domain_radius(::GGJDispersion, Q::ComplexF64) = 0.9 * max(imag(Q), 0.0)
+
+"""
+    ggj_growth_rates(mc::MultiSurfaceCouplingFull;
+                     gamma_range, omega_range,
+                     nre0=40, nim0=40, passes=4, max_cells=2_000_000,
+                     polish=true, winding_rfrac=0.5,
+                     parallel=Threads.nthreads() > 1) -> (GrowthRateResult, AMRResult)
+
+Tearing growth rates of the GGJ inner layer matched through the full 4m × 4m
+determinant `mc` (coupled for several surfaces, uncoupled for one surface with
+its own 2 × 2 block of the raw Δ'). The box `gamma_range` × `omega_range`
+(physical: γ in s⁻¹, ω in rad/s) is scanned adaptively in the plane of
+[`GGJDispersion`](@ref) and every contour intersection is classified by its
+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.
+
+`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
+supports only γ > 0; classification circles and polishing stay in that half-plane,
+including finite-difference probes. Stable-root searches and analytic continuation
+across the cut are not implemented. A solver failure inside the
+layer (for example a degenerate Frobenius pair near Mercier marginality) is
+not caught and stops the scan.
+
+`nre0` × `nim0` is the initial AMR grid and `passes` the number of refinement
+passes. Each determinant evaluation solves the inner layer at every surface:
+the default 40 × 40 grid alone is about 1,600 evaluations, and each contour
+intersection adds 64–4096 more for its classifying circles plus up to about
+120 for the polish. Size the grid for the cost of one inner-layer solve.
+`:no_root` and `:winding_unresolved` results warn; `gamma_Hz` is then not a
+growth rate.
+"""
+function ggj_growth_rates(mc::MultiSurfaceCouplingFull;
+    gamma_range::NTuple{2,<:Real}, omega_range::NTuple{2,<:Real},
+    nre0::Integer=40, nim0::Integer=40, passes::Integer=4,
+    max_cells::Integer=2_000_000,
+    polish::Bool=true, winding_rfrac::Real=0.5,
+    parallel::Bool=Threads.nthreads() > 1)
+    gamma_range[1] > 0 && gamma_range[2] > gamma_range[1] ||
+        throw(ArgumentError("ggj_growth_rates: gamma_range must satisfy 0 < γ_min < γ_max; got $(gamma_range). " *
+                            "The GGJ layer's branch cut runs along γ ≤ 0."))
+    omega_range[2] > omega_range[1] ||
+        throw(ArgumentError("ggj_growth_rates: omega_range must be increasing; got $(omega_range)"))
+    g = GGJDispersion(mc)
+    τ = 1 / g.Q0
+    re_range = (-Float64(omega_range[2]) * τ, -Float64(omega_range[1]) * τ)
+    im_range = (Float64(gamma_range[1]) * τ, Float64(gamma_range[2]) * τ)
+    amr = amr_scan(g, re_range, im_range; nre0=nre0, nim0=nim0, passes=passes,
+        max_cells=max_cells, parallel=parallel)
+    gr = find_growth_rates(amr, τ; residual=g, polish=polish,
+        classification=:winding, zero_winding=1, winding_rfrac=winding_rfrac)
+    :no_root in gr.warning_flags && @warn "ggj_growth_rates: no root found in the box; gamma_Hz = 0 is a placeholder"
+    :winding_unresolved in gr.warning_flags &&
+        @warn "ggj_growth_rates: some contour intersections could not be classified (the determinant failed or had a zero or pole on the circle)"
+    return gr, amr
+end

src/Tearing/Dispersion/GrowthRateExtraction.jl

@@ -599,8 +599,12 @@ end
 #
 # The local model is the real 2×2 forward-difference Jacobian of (Re f, Im f), exact for holomorphic
 # and anti-holomorphic f alike; SLAYER's residual is a function of conj(Q) (see Riccati.jl).
+# Radius of a disk contained in the residual's supported domain; unrestricted by default.
+_residual_domain_radius(f, Q::ComplexF64) = Inf
+
 function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8,
     tol_step::Float64=1e-12, tol_f::Float64=1e-8)
+    R = min(R, _residual_domain_radius(f, Q0))
     (R > 0) || return (Q0, 0, false)
     f0 = ComplexF64(f(Q0))
     nev = 1
@@ -613,6 +617,8 @@ function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8,
     Qbest, abest = Q0, a0
 
     for _ in 1:maxit
+        h = min(h, 0.5 * _residual_domain_radius(f, Qk))
+        h > 0 || break
         fx = ComplexF64(f(Qk + h))
         fy = ComplexF64(f(Qk + im * h))
         nev += 2
@@ -974,7 +980,9 @@ function _winding_classify(residual, pts::Vector{ComplexF64}, i::Int, re_paths,
     zero_winding::Int, rfrac::Float64)
     pt = pts[i]
     h_loc = _local_segment_length(re_paths, im_paths, pt, h_grid)
-    r = _polish_trust_radius(pts, i, h_loc; k_cell=rfrac, nn_frac=0.45)
+    domain_radius = _residual_domain_radius(residual, pt)
+    r = min(_polish_trust_radius(pts, i, h_loc; k_cell=rfrac, nn_frac=0.45), domain_radius)
+    r > 0 || return (:unresolved, NaN, r)
     W = NaN
     for _ in 1:3
         W, _ = _winding_number(residual, pt, r)
@@ -984,7 +992,7 @@ function _winding_classify(residual, pts::Vector{ComplexF64}, i::Int, re_paths,
     isfinite(W) || return (:unresolved, NaN, r)
     n = round(Int, W)
     if n == 0
-        r_wide = _polish_trust_radius(pts, i, h_loc; k_cell=8rfrac, nn_frac=0.9)
+        r_wide = min(_polish_trust_radius(pts, i, h_loc; k_cell=8rfrac, nn_frac=0.9), domain_radius)
         if r_wide > r
             W_wide, _ = _winding_number(residual, pt, r_wide)
             if isfinite(W_wide) && round(Int, W_wide) != 0

src/Tearing/Runner/run_slayer.jl

@@ -128,8 +128,7 @@ end
 _surface_label(p::SLAYERParameters) = "m=$(p.m), n=$(p.n)"
 _surface_label(p::InnerLayerParameters) = "ising=$(getfield(p, :ising))"
 
-# Is this an unvalidated GGJ inner-layer model? Used to gate the γ-extraction
-# future-work warning.
+# Is this a GGJ inner-layer model? Its growth rates come from `ggj_growth_rates`, not this runner.
 _is_ggj(::GGJModel) = true
 _is_ggj(::Any) = false
 
@@ -194,12 +193,6 @@ function _build_surface_coupling(model::SLAYERModel, params::SLAYERParameters,
     return surface_coupling(model, params, dp_diag; dc=params.dc_tmp, q_shift=q_shift)
 end
 
-# GGJ: scale = 1.0 (rescale_delta applied inside solve_inner), tauk = 1.0,
-# dc = 0 (the 4m×4m Pletzer-Dewar residual carries interchange stabilization
-# natively), and no E×B shift. See the GGJ `surface_coupling` method.
-_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
@@ -317,6 +310,10 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel,
     rational_q::Vector{Float64}=Float64[],
     omega_E::Vector{Float64}=Float64[])
     validate(control)
+    _is_ggj(model) && throw(ArgumentError(
+        "run_slayer_from_inputs: this runner scans SLAYER's tearing-only residual, which gives no GGJ growth rate. " *
+        "Use Dispersion.ggj_growth_rates on multi_surface_coupling_full, which matches both parity channels, and " *
+        "ggj_inner_deltas for the GGJ layer's Δ."))
     control.enabled || return empty_slayer_result(control)
     isempty(params) && return empty_slayer_result(control)
 
@@ -328,60 +325,25 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel,
 
     # Guard: the inner-layer model and the parameter eltype must match, or
     # `_build_surface_coupling` throws an opaque MethodError downstream.
-    expected_P = _is_ggj(model) ? GGJParameters : SLAYERParameters
+    expected_P = SLAYERParameters
     all(p -> p isa expected_P, params) ||
         throw(
             ArgumentError(
                 "run_slayer: a $(typeof(model)) inner model requires " *
                 "$(expected_P) per-surface parameters, but got eltype " *
-                "$(eltype(params)). Build inputs with the matching builder " *
-                "(build_slayer_inputs for SLAYER, ggj_parameters for GGJ).")
-        )
-
-    # Slab-layer path: convert Δ' from its ψ_N reference length to the
-    # r_s-based convention shared by the layer Δ(Q) and the critical-Δ (see
-    # `delta_prime_to_rs_reference`). GGJ is genuinely toroidal/ψ-based
-    # (its `rescale_delta` handles inner→outer units natively) — no conversion.
-    if !_is_ggj(model)
-        dp = delta_prime_to_rs_reference(dp, params)
-    end
-
-    # The coupled determinant uses the reduced m×m (tearing-only) form, which
-    # drops the interchange channel. For GGJ that channel carries the Glasser
-    # interchange stabilization, so coupled-GGJ results omit real physics.
-    if _is_ggj(model) && control.coupling_mode === :coupled
-        @warn(
-            "SLAYER: coupling_mode=:coupled with a GGJ inner model uses the " *
-            "reduced m×m tearing-only determinant, which DROPS the GGJ " *
-            "interchange (Glasser stabilization) channel — results are " *
-            "physically incomplete. Use the full Pletzer-Dewar matching " *
-            "(multi_surface_coupling_full) for coupled GGJ studies."
+                "$(eltype(params)). Build them with build_slayer_inputs.")
         )
-    end
 
-    # GGJ growth-rate extraction by Re/Im contour matching is not yet
-    # reliable (the contours do not robustly intersect at a dispersion-relation
-    # zero). The scan still runs so the user can inspect/plot the contours and
-    # access both parity Δ channels via `ggj_inner_deltas`, but the reported γ
-    # is unvalidated and should be treated as future work.
-    if _is_ggj(model)
-        @warn(
-            "SLAYER: GGJ γ-extraction by Re/Im contour matching is " *
-            "FUTURE WORK — contours do not reliably intersect at a " *
-            "dispersion-relation root. The scan grid and parity Δ channels " *
-            "are provided for inspection; reported growth rates are " *
-            "unvalidated placeholders."
-        )
-    end
+    # Convert Δ' from its ψ_N reference length to the r_s-based convention shared by the layer Δ(Q) and the
+    # critical-Δ (see `delta_prime_to_rs_reference`).
+    dp = delta_prime_to_rs_reference(dp, params)
 
     # 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
     Ω_E = _omega_E_per_surface(control, params, omega_E)
     # A surface with no known rotation applies no shift.
     Ω_shift = [isnan(Ω) ? 0.0 : Ω for Ω in Ω_E]
-    coupled && _is_ggj(model) && any(!iszero, Ω_shift) && @warn(
-        "SLAYER: E×B rotation is not applied to GGJ surfaces, which carry no time normalization.")
     q_shifts = coupled ? Float64[_q_shift(params[k], Ω_shift[k]) for k in 1:n] : zeros(Float64, n)
     scs = [_build_surface_coupling(model, params[k], dp[k, k], q_shifts[k]) for k in 1:n]
     # Coupled Q is scanned in the reference surface's normalization.
@@ -559,10 +521,9 @@ vector of `(ising, tearing, interchange)` named tuples, where `tearing`
 (GWP Δ₋, odd) is the reconnecting channel and `interchange` (GWP Δ₊, even) the
 non-reconnecting Glasser-stabilization channel (see `InnerLayerResponse`).
 
-This is the supported GGJ diagnostic: both parity Δ's are physical and
-directly accessible here. Contour-matching γ extraction from these (the
-Re/Im scan intersection) is not yet validated — see the warning in
-`run_slayer_from_inputs`.
+Both parity Δ's are physical and directly accessible here. For GGJ growth
+rates use `Dispersion.ggj_growth_rates`, which matches both channels through
+the full 4m×4m determinant.
 """
 function ggj_inner_deltas(params::AbstractVector{GGJParameters}, Q::Number;
     solver::Symbol=:galerkin)

test/runtests.jl

@@ -65,6 +65,7 @@ else
     include("./runtests_dispersion_polish.jl")
     include("./runtests_slayer_profile_transport.jl")
     include("./runtests_dispersion_winding.jl")
+    include("./runtests_dispersion_ggj.jl")
     include("./runtests_slayer_runner.jl")
     include("./runtests_critical_resonant_field.jl")
     include("./runtests_kinetic.jl")