diff --git a/src/Tearing/Dispersion/Dispersion.jl b/src/Tearing/Dispersion/Dispersion.jl index 67731e5aa..2c503a452 100644 --- a/src/Tearing/Dispersion/Dispersion.jl +++ b/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 diff --git a/src/Tearing/Dispersion/GGJGrowthRates.jl b/src/Tearing/Dispersion/GGJGrowthRates.jl new file mode 100644 index 000000000..8d3927865 --- /dev/null +++ b/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 diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index c21408965..685fbea8a 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/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 diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index df9d76be7..26ba0c18c 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/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,51 +325,18 @@ 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. @@ -380,8 +344,6 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel, Ω_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) diff --git a/test/runtests.jl b/test/runtests.jl index 8f93b8c23..7e3253182 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/runtests_dispersion_ggj.jl b/test/runtests_dispersion_ggj.jl new file mode 100644 index 000000000..d5bae4159 --- /dev/null +++ b/test/runtests_dispersion_ggj.jl @@ -0,0 +1,125 @@ +struct GGJDomainTestResidual{G,F} + dispersion::G + residual::F + evaluated::Vector{ComplexF64} +end + +function (f::GGJDomainTestResidual)(Q) + push!(f.evaluated, ComplexF64(Q)) + return f.residual(Q) +end + +GeneralizedPerturbedEquilibrium.Dispersion._residual_domain_radius(f::GGJDomainTestResidual, Q::ComplexF64) = + GeneralizedPerturbedEquilibrium.Dispersion._residual_domain_radius(f.dispersion, Q) + +@testset "Dispersion GGJ growth rates (full matching, winding classification)" begin + using GeneralizedPerturbedEquilibrium.InnerLayer + using GeneralizedPerturbedEquilibrium.InnerLayer: GGJModel, solve_inner + using GeneralizedPerturbedEquilibrium.Dispersion + GGJ = GeneralizedPerturbedEquilibrium.InnerLayer.GGJ + + # Glasser & Wang (2020) Eq. 55 layer, one surface. The 2 × 2 raw Δ' block is + # built in the parity basis P = M⁻¹·Dᵀ·M (M = [1 −1; 1 1]) with P₁₂ = P₂₁ = 0, + # for which the 4 × 4 determinant vanishes where Δ_tear(s) = P₂₂. Setting + # P₂₂ = Δ_tear(s₀) puts a root at the real growth rate s₀. + p = GGJ.glasser_wang_2020_eq55() + Q0 = GGJ.q0(p) + model = GGJModel(; solver=:ray) + s0 = 0.1234 * Q0 + r0 = solve_inner(model, p, s0) + P = ComplexF64[real(r0.interchange)/2 0; 0 r0.tearing] + M = ComplexF64[1 -1; 1 1] + dp_raw = Matrix(transpose(M * P / M)) + mc = multi_surface_coupling_full([surface_coupling(model, p, r0.tearing)], dp_raw) + + @testset "GGJ classification and polishing stay at positive growth" begin + D = GeneralizedPerturbedEquilibrium.Dispersion + g = GGJDispersion(mc) + center = 0.01im + re_paths = [ComplexF64[-1 + center, 1 + center]] + im_paths = [ComplexF64[0.001im, 0.02im]] + # A coarse contour estimate must not pull a stable root into a positive-growth search. + f = GGJDomainTestResidual(g, Q -> Q + 0.1im, ComplexF64[]) + gr = D._run_analysis(re_paths, im_paths, [zeros(2)], 1.0; + pole_threshold=10.0, filter_above_poles=false, filter_outside_re=false, + residual=f, classification=:winding, zero_winding=1) + @test :no_root in gr.warning_flags + @test isempty(gr.valid_roots) + @test !isempty(f.evaluated) && all(Q -> imag(Q) > 0, f.evaluated) + + # Expansion may find a nearby unstable root, but its returned trust disk must remain in-domain. + nearby = GGJDomainTestResidual(g, Q -> Q - 0.014im, ComplexF64[]) + kind, W, radius = D._winding_classify(nearby, ComplexF64[center], 1, re_paths, im_paths, 2.0; + zero_winding=1, rfrac=0.001) + @test kind === :root + @test W ≈ 1 + @test 0 < radius < imag(center) + root, _, _ = D._polish_root(nearby, center, radius) + @test root ≈ 0.014im atol = 1e-10 + @test all(Q -> imag(Q) > 0, nearby.evaluated) + + # Even an oversized requested trust region cannot move the iterate or its probes below γ = 0. + empty!(f.evaluated) + root, _, _ = D._polish_root(f, center, 1.0) + @test imag(root) > 0 + @test !isempty(f.evaluated) && all(Q -> imag(Q) > 0, f.evaluated) + empty!(f.evaluated) + @test D._polish_root(f, 0.0im, 1.0) == (0.0im, 0, false) + kind, _, _ = D._winding_classify(f, ComplexF64[0.0im], 1, re_paths, im_paths, 2.0; + zero_winding=1, rfrac=0.5) + @test kind === :unresolved + @test isempty(f.evaluated) + end + + @testset "GGJDispersion: plane convention s = −i·Q·Q₀" begin + g = GGJDispersion(mc) + @test g.Q0 ≈ Q0 + Q = 0.03 + 0.11im + @test g(Q) == mc(-im * Q * Q0) + @test abs(g(0.1234im)) < 1e-8 * abs(g(0.03 + 0.11im)) # the constructed root + @test_throws ArgumentError GGJDispersion( + multi_surface_coupling_full( + [surface_coupling(GeneralizedPerturbedEquilibrium.InnerLayer.SLAYERModel(), nothing, 0.0im; scale=1.0, tauk=1.0)], dp_raw) + ) + end + + @testset "ggj_growth_rates finds the constructed root, winding +1" begin + gr, amr = ggj_growth_rates(mc; gamma_range=(0.08Q0, 0.17Q0), omega_range=(-0.031Q0, 0.047Q0), + nre0=6, nim0=6, passes=3, parallel=false) + @test gr.classification === :winding + @test !(:no_root in gr.warning_flags) + @test all(w -> w ≈ 1, gr.valid_roots_winding) + @test gr.gamma_Hz ≈ s0 rtol = 1e-6 + @test abs(gr.omega_Hz) < 1e-6 * Q0 + @test abs(gr.Q_root - 0.1234im) < 1e-6 + # The GGJ layer's branch cut runs along γ ≤ 0, so a box reaching there is refused. + @test_throws ArgumentError ggj_growth_rates(mc; gamma_range=(-0.01Q0, 0.17Q0), omega_range=(-0.031Q0, 0.047Q0)) + @test_throws ArgumentError ggj_growth_rates(mc; gamma_range=(0.08Q0, 0.17Q0), omega_range=(0.047Q0, -0.031Q0)) + end + + @testset "two coupled surfaces: the root agrees with a direct solve of the 8 × 8 determinant" begin + # A second surface with its own constructed root far outside the box, and a symmetric coupling of the two + # surfaces' tearing channels (the off-diagonal block c·[1 −1; −1 1] reaches d₋ through the matching identity), + # which moves surface 1's root off its uncoupled value. + p2 = GGJ.q4_surface_benchmark() + r02 = solve_inner(model, p2, 0.5 * Q0) + blk(r) = Matrix(transpose(M * ComplexF64[real(r.interchange)/2 0; 0 r.tearing] / M)) + B = 0.1 * abs(r0.tearing) * ComplexF64[1 -1; -1 1] + dp2 = [blk(r0) B; transpose(B) blk(r02)] + 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 + end + @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) + @test length(gr2.valid_roots) == 1 && all(w -> w ≈ 1, gr2.valid_roots_winding) + @test abs(gr2.Q_root - b) < 1e-8 * abs(b) + end +end