[d-burg]

Release note

• Audience: users
• Numerical impact: existing harness quantities unchanged; GGJ root extraction now stays at positive growth and rejects out-of-domain stable roots (harness @ 8a6ae1b)
• Migration: run_slayer_from_inputs now throws for a GGJ inner model instead of returning placeholder growth rates; use ggj_growth_rates (and ggj_inner_deltas for the layer's Δ).

GGJ tearing growth rates can now be computed in the package. ggj_growth_rates(multi_surface_coupling_full(...); gamma_range, omega_range) scans the full Pletzer-Dewar 4m × 4m determinant, which keeps the interchange channel and therefore the Glasser stabilization, and classifies every root by winding number. It works for coupled surfaces, and for a single surface through its own 2 × 2 block of the raw Δ′. Results come back as a GrowthRateResult, as for SLAYER.

Changes since the review

• 8a6ae1bdd: constrain GGJ classification circles, expanded retries, polishing trust regions and derivative probes to gamma > 0; other residuals remain unrestricted. Add regression tests for stable-root escape, nearby unstable roots and boundary rejection.
• 0ed7b88ac: the two-surface test against a direct solve, the γ-range validation, the :no_root/:winding_unresolved warnings, the cost and frame notes in the docstring, and the runner's GGJ path turned into an error with its dead branches removed.
• Tearing - FEATURE - Classify tearing roots by winding number #597's review rework merged in (81fe84749), so ggj_growth_rates uses the consolidated winding classifier.

Regression report

Latest full SDCC harness: 36d54d2fa (parent) versus 8a6ae1bdd, all 22 cases / 44 tasks on sirius, Julia 1.12.6, pinned manifest 4c9d54fa, 8 Julia threads / 4 BLAS threads. All 317 comparable quantities across 21 cases are unchanged, including uncoupled, alternate-transport, coupled and winding SLAYER cases. All tasks and report job 1300503 completed with exit 0. The existing eight gal_resistive_pe quantities remain missing/N/A on both refs and provide no numerical validation. Report: /home/ITER/pharrm/Projects/GPEC_dev/regress-fanout/runs/pr598_8a6ae1bdd_full/report.txt.

GGJ domain-fix validation: a separate SDCC study tested the exact source patch committed as 8a6ae1bdd, with the production classifier/polisher active and determinant evaluations instrumented to reject gamma <= 0. All eight study cases completed, with zero nonpositive-growth evaluations and results matching the earlier domain-limited prototype. In the DIII-D-like three-surface coupled case, the original implementation returned gamma = -876.82 s^-1 from a requested positive interval of 0.72–14410 s^-1; the fix reports no_root on both coarse and refined scans. This establishes containment, not proof that the search found every root. The other seven comparisons were unchanged. The study used actual raw DIII-D-like matching matrices, SLAYER-derived resistivity and rho = ion_massM_Pn_e, with no rotation shift; constructed-root benchmarks used manufactured outer matrices and real GGJ kernels. Both versions still miss the manufactured normalized-growth 0.0003 root in the tested scans, so search completeness remains a limitation.

All 175 focused SDCC checks passed: GGJ growth rates (28), full matching determinant (26), winding (56), polishing (65). Study array 1300477, report job 1300483; report: /home/ITER/pharrm/Projects/GPEC_dev/pr598-domain-fix-20261029/report.txt. The general harness has no coupled GGJ growth-extraction case, so this targeted study is reported separately.

Earlier harness report (before the domain fix):

regress --refs fbdd6eb07,0ed7b88ac --force (all cases, against #597's head) on feynman (SLURM, 8 threads, Julia 1.11.9, manifest pinned): every case is unchanged, including diiid_slayer_n1_winding (17) and diiid_slayer_n1_coupled (7); gal_resistive_pe is N/A on both refs. No harness case runs the coupled GGJ path, which is reachable only through the API; the tests below cover it. runtests_dispersion_ggj.jl (15), runtests_dispersion_coupled_full.jl (26), runtests_dispersion_winding.jl (56) and runtests_slayer_runner.jl (156) pass on feynman at 0ed7b88ac.

Notes for reviewers

What was missing. The package could evaluate the GGJ determinant (multi_surface_coupling_full) but had no root finder for it. The runner rejects ggj_* inner models in SLAYERControl. Its GGJ path, reachable through run_slayer_from_inputs, scans the tearing-only residual and warned that contour matching did not reliably find GGJ roots and that its growth rates are placeholders. GGJ γ's have therefore been computed by out-of-tree drivers.

What this adds (src/Tearing/Dispersion/GGJGrowthRates.jl, about 70 lines):

• GGJDispersion(mc) maps the determinant into the Q plane the SLAYER scans use, g(Q) = mc(−i·Q·Q₀), with Q₀ = S^(−1/3)/τ_A of the reference surface (GWP2016 Eq. A16). The determinant takes s = γ + iω, and the SLAYER convention is Q = τ·(−ω + iγ), so s = −iQ/τ with τ = 1/Q₀. As a result the AMR scan, the contour extraction and GrowthRateResult (gamma_Hz = γ in s⁻¹, omega_Hz = −ω in rad/s) carry over unchanged.
• ggj_growth_rates runs amr_scan over a physical box and find_growth_rates(...; classification = :winding, zero_winding = +1) from Tearing - FEATURE - Classify tearing roots by winding number #597. The map s ↦ Q is holomorphic, so the determinant's zeros wind +1, the opposite of the SLAYER residual.

The method is the one the out-of-tree driver used on the TJ circular β scans: AMR in the scaled variable, winding-number acceptance, polish. The difference is that candidates are the shared classifier's contour intersections rather than the driver's per-cell windings, so SLAYER and GGJ go through one classifier.

Tests.

• One surface: the Glasser & Wang (2020) Eq. 55 layer (:ray), with the raw Δ′ block built in the parity basis P = M⁻¹DᵀM with P₁₂ = P₂₁ = 0, which makes the determinant (P₁₁ − Δ_int)(P₂₂ − Δ_tear); P₂₂ = Δ_tear(s₀) puts a root exactly at s₀ = 0.1234·Q₀. Checks the plane convention, winding +1 and γ = s₀ to 1e-6.
• Two coupled surfaces: the same layer plus the q = 4 benchmark surface, with a symmetric coupling of their tearing channels (off-diagonal block c·[1 −1; −1 1], which reaches d₋ through the matching identity) that moves the root by 5.9e-5. ggj_growth_rates agrees with an independent complex-secant solve of the 8 × 8 determinant to 1.5e-10 relative, which also shows that the polish converges on the coupled GGJ determinant.
• The γ-range validation and the runner's new error.

Runner. run_slayer_from_inputs now throws for a GGJ inner model and points here and to ggj_inner_deltas, instead of scanning the tearing-only residual and returning placeholder growth rates; its GGJ-only branches are removed. Routing [SLAYER] decks with GGJ through ggj_growth_rates needs the raw 2m × 2m Δ′ in the runner (it receives the m × m matrix today) and would let the function reuse SLAYERControl's scan settings; that is left for a follow-up.

Domain. gamma_range must lie at γ > 0. The principal square-root cut is on the negative real s-axis. Classification circles, expanded retries, polishing steps and finite-difference probes now remain strictly in the supported positive-growth half-plane; stable-root searches and continuation across the cut are not implemented. A solver failure inside the layer stops the scan rather than being turned into NaN. :no_root and :winding_unresolved results warn. The docstring gives the cost: the default 40 × 40 grid alone is about 1,600 determinant evaluations, each solving every surface.

Caveats. 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 roots are in the plasma frame.

Stacked on #597 (feature/tearing-winding-classification), which is stacked on #577.

---

🤖 Generated with Claude Code

[matt-pharr]

@d-burg I haven't fully read through these 3 PRs yet but I can see that you did not implement this through the new API standard, and solve(TearingProblem(),GGJModel()) still errors... IMO the public syntax for this should look like this:

prob = TearingProblem(ffs; eta=eta, rho=rho, coupling_mode=:coupled)
model = GGJModel(; solver=:ray,
    gamma_range=(γmin, γmax), omega_range=(ωmin, ωmax))
result = solve(prob, model)
save_to_hdf5("gpec.h5", result)

If you agree, I will work on reviewing these as-is and get them merged, and then I can open a new PR on top of #602 and add you as a reviewer on that, where I just purely move the new functionality you have here into the planned public API above (I would prefer to do it myself to lessen the need for back-and-forths as I am still working on #602 which is the big cleanup PR).

[matt-pharr]

We made a small change to keep classification circles, expand retries, polishing steps, and derivative probes inside γ > 0. This fixes a demonstrated DIII-D-like coupled case that returned γ ≈ −877 s⁻¹ despite a positive-growth search. Other residuals, including SLAYER, remain unrestricted. Added regression tests and corrected the branch-cut documentation.

Work to port this to the new API spec remains outstanding and I will do that on a PR stacked on top of V.