Tearing Module

build_ggj_inputs(equil, sings, profiles; mu_i=2.0, zeff=1.0,
                  v1_scale=1.0,
                  resistivity_model::NeoResistivityModel=SpitzerModel(),
                  lnLambda_form::Symbol=:nrl) -> Vector{GGJParameters}

Construct a GGJParameters for each rational surface in sings. Each surface's geometric coefficients (E, F, G, H, K, M) come from the sing.restype::ResistGeometry populated by resist_eval_all!. Kinetic timescales are derived from the KineticProfiles at sing.psifac:

ρ(ψ)   = μ_i · m_p · n_e(ψ)
η(ψ)   = eta_neoclassical(model, n_e, T_e, Z_eff, f_t, ν*_e)     [Ω·m]
τ_A    = √(ρ · M · μ_0) / |2π · n · q' · χ₁ / V'|                 [Alfvén time]
τ_R    = (⟨B²/|∇ψ|²⟩ / ⟨B²⟩) · μ_0 / η                             [resistive diffusion]

The mode number n is taken from sings[k].n[1] (first resonant mode at the surface). χ₁ = 2π · psio. The v1_scale kwarg is an optional multiplicative factor on V' in the τ_A denominator (a v1 / volume normalization option); default 1.0 means use the raw V'.

Resistivity model

resistivity_model selects the η closure:

• SpitzerModel() (default) — Sauter 1999 Eq. 18a (Zeff-aware Spitzer), with the NRL Coulomb log.
• SauterNeoModel() — multiplies by Sauter 1999 F33 using ft and ν*_e from the surface's ResistGeometry. Produces the physically-correct trapped-particle-corrected η for H-mode tearing stability.
• RedlNeoModel() — Redl 2021 F_33 (improved high-ν* fit).

lnLambda_form selects :nrl (default), :sauter, or :wesson.

Throws if any surface's restype is still nothing — call ForceFreeStates.resist_eval_all!(intr, equil) first.

InnerLayerModel

Abstract supertype for resistive inner-layer models. Each concrete model is a small, parameter-free type tag (often parameterized by a solver-choice symbol) that selects a solve_inner method.

Implementations live in submodules of InnerLayer, e.g. InnerLayer.GGJ.

InnerLayerParameters

Abstract supertype for the per-surface physical-parameter structs consumed by the inner-layer models (SLAYERParameters, GGJParameters). Lets the dispersion runner and HDF5 output dispatch generically over a vector of single-model parameters.

InnerLayerResponse

Parity-projected inner-layer matching data at one rational surface. The two components correspond to the homogeneous parity solutions of the half-domain inner-layer problem (parity boundary conditions imposed at X = 0). They are the Δ_{j,±}(γ) of Glasser, Wang & Park, Phys. Plasmas 23, 112506 (2016), Eqs. (34)–(35).

Fields

• tearing — the odd-parity matching coefficient (GWP Δ_+, the "odd mode"). Corresponds to a flux perturbation W that is EVEN in x and a velocity/temperature perturbation that is ODD — i.e., the reconnecting mode with a current sheet at the rational surface. This is the tearing drive that appears as Δ' in the classical constant-ψ tearing equation. Must be populated by every resistive inner-layer model.

• interchange — the even-parity matching coefficient (GWP Δ_−, the "even mode"). Corresponds to W odd, N and Θ even — i.e., the non-reconnecting interchange/ballooning channel. Its dissipative piece in toroidal geometry is the Glasser, Greene & Johnson stabilization term that opposes tearing growth (Glasser, Greene & Johnson 1975; Lütjens-Bondeson-Roy 1993). Pressureless inner-layer models (e.g. SLAYER's Fitzpatrick Riccati) set this identically zero.

The naming follows the physics channel rather than a mathematical parity label because odd/even carries different meanings across the literature depending on whether you label by the parity of W (GWP paper convention) or the parity of (N, Θ). Using tearing and interchange avoids ambiguity.

solve_inner(model::InnerLayerModel, params, γ::Number; kwargs...) -> InnerLayerResponse

Compute the parity-projected matching data (Δ_tearing, Δ_interchange) for the given inner-layer model, physical parameters params, and complex growth rate γ. Concrete models specialize this function.

See InnerLayerResponse for the physics-oriented field definitions. Pressureless models (SLAYER) populate only tearing and leave interchange at zero; two-fluid / finite-β models (GGJ) populate both.

solve_inner_profile(model::InnerLayerModel, params, γ::Number; kwargs...)
    -> (; Δ, x, Ψ, Ξ, dψdx, rescale, ...)

Compute the inner-layer matching data and the reconstructed layer field profiles for the given model — everything an outer↔inner matching driver needs from the layer, so drivers never touch model internals. Returns a named tuple with at least:

• Δ — the same (Δ_odd, Δ_even) matching data as solve_inner
• x — real ascending grid in the model's stretched inner coordinate, x ≥ 0 with the rational surface at x = 0
• Ψ, Ξ — length(x) × 2 profiles, columns (odd, even) parity, in the model's inner normalization: Ψ the normal-field (reconnected-flux) variable, Ξ the displacement
• dψdx — conversion to poloidal-flux distance, δψ = dψdx · x
• rescale — amplitude factor converting the inner-normalized profiles to the outer δψ-normalized convention (companion of the Δ rescale)

Concrete models may return additional diagnostic fields (e.g. a solve-quality certificate). Solver-knob keywords are model-specific.

GGJModel{S} <: InnerLayerModel

Glasser–Greene–Johnson resistive inner-layer model. S selects the solver backend: :ray (default; robust at large |Q| on/near the imaginary axis), :galerkin (real-axis Hermite FEM; degrades for |Q| ≳ 1), or :shooting (|Q| ≪ 1 only). The backends take different numerical-knob keywords.

GGJParameters

Dimensionless parameters of the Glasser–Greene–Johnson inner-layer model at a single rational surface, plus the local Alfvén/resistive timescales needed to scale the matching data back to physical Δ. The equilibrium coefficients E, F, G, H, K, M are the flux-surface averages defined in GWP2016 Eq. (A8); they enter the inner-region equations (Eq. 11).

Fields:

The growth rate γ is not stored here; it is a separate argument to solve_inner.

GGJShootingSystem

Precomputed origin- and infinity-side arrays for the backward shoot. Construct via _build_shooting_system.

InnerAsymptoticsCache

Frozen inps Wasow asymptotic-basis construction for one (GGJParameters, Q) pair, stored as SMatrix/SVector for allocation-free evaluation. Index convention: P[k+1] holds the k-th-order matrix P_k, etc.

Fields:

• params, Q, kmax — input parameters and series truncation order.
• λ = 1/√Q — complex scale factor used by the Wasow split.
• R = (r₊, r₋) — Mercier-shifted Frobenius exponents at infinity (Eq. 49).
• T, Tinv — 6×6 eigenvector basis of A₀ (Eq. 7–8).
• J — (J₀, J₁, J₂), the J-rotated coefficient matrices (Eq. 9–10).
• P, B — splitting matrices, k = 0..kmax+2 (Eqs. 16, 22).
• K2 — 2×2 inner working matrices, k = 0..kmax+2 (Eq. 32; entry k=0 unused).
• Qmat, Cmat — 2×2 inner-block transformation matrices, k = 0..kmax+2 (Eqs. 32–39).
• Dmat — 2×2 shearing matrices, k = 0..kmax (Eq. 43).
• Y0, Y0inv — lowest-order Y matrix and its inverse (Eq. 48; exponents from Eq. 49).
• Y, Z — 2×2 series matrices, k = 0..kmax (Eq. 52).

RaySolveResult

Output of solve_ray. Δ is the rescaled matching data in the deltac output ordering (index 1 ↔ even-BC solution, 2 ↔ odd-BC solution, i.e. the same swap deltacsolve applies); Δraw is (isol=1, isol=2) in raw inps normalization before `rescaledelta`.

_build_shooting_system(p::GGJParameters, Q::ComplexF64; nps=8, rtol=1e-6, fmax=1.0)

Construct a GGJShootingSystem for the given parameters p and complex frequency Q, precomputing the origin and infinity asymptotic arrays used by the forward/backward shoots.

Convert the inps 6×2 output U_inps (the Wasow asymptotic basis U = TPQSY of GW2020 Eq. 53) at coordinate x to the physical (Ψ, Ξ, Υ) and (Ψ', Ξ', Υ') representation. The 6-component first-order state packs (Ψ, Ξ, Υ, Ψ', Ξ', Υ') with the GW2020 Eq. (2) scaling 𝚿 ≡ (xΨ, Ξ, Υ), hence the /x and ·x factors below. Returns (ua, dua) each 3×2 complex, where columns are the two power-like (Mercier) solutions and rows are the components (Ψ, Ξ, Υ).

Build the (I, U, V) coefficient matrices of the second-order system I·u'' − V·u' − U·u = 0 for u = (Ψ, Ξ, Υ) at coordinate x. All matrices are 3×3 complex.

These are the matrices A, B, C of GWP2016 Eq. (12), A Ψ'' + B Ψ' + C Ψ = 0, with A given in Eq. (14), B in Eq. (14), and C in Eq. (15). The code's weak-form layout flips the off-diagonal signs: (I, V, U) = (A, −B, −C).

Each matrix row is one GGJ inner-region equation (GWP2016 Eq. 11 ≡ GW2020 Eq. 1); reading I u'' − V u' − U u across a row reproduces the corresponding equation:

row 1 (Ψ): Ψxx − H Υx − Q(Ψ − x Ξ) = 0 ⇒ I=(1,0,0) V=(0,0,H) U=(Q, −Qx, 0) row 2 (Ξ): Q² Ξxx − Q x² Ξ + Q x Ψ + (E+F) Υ + H Ψx = 0 ⇒ I=(0,Q²,0) V=(−H,0,0) U=(−Qx, Q x², −(E+F)) row 3 (Υ): Q Υxx − x² Υ + x Ψ + Q²[G(Ξ−Υ) − K(E Ξ + F Υ + H Ψx)] = 0 ⇒ I=(0,0,Q) V=(K Q² H,0,0) U=(−x, −Q²(G−KE), x²+Q²(G+KF))

_solve_parities(Ic, Jc, Vc, rhs, ndof, N, p) -> (sol_odd, sol_even)

Solve both parity systems from one factorization: the two matrices differ only in the 3 parity rows at s = 0 (a rank-3 update), so one sparse LU of the row-equilibrated odd system plus a Woodbury correction replaces the second factorization, which dominates the linear-algebra cost.

asymptotic_profile(params, res::RaySolveResult, srange) -> (; s, x, Ψ, Ξ, Υ)

Evaluate the analytic representation u_small + Δ·u_big on the ray for s ≥ res.S; concatenated with solution_profile it should join the numeric solution seamlessly across the march zone.

asymptotic_residual(cache, x::Complex) -> SVector{2,Float64}

GW2020 Eq. (54) residual of the two series columns at complex x.

asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) -> SVector{2,Float64}

Compute the convergence measure Δ± of the asymptotic basis at x for each of the two algebraic columns (GW2020 Eq. 54). Returns ‖dU − x·matrix·U‖∞ / max(‖dU‖∞, ‖x·matrix·U‖∞) per column, where matrix = J₀ + xfac·J₁ + xfac²·J₂ is the J-rotated coefficient matrix (the residual of v' = xJv, GW2020 Eq. 6).

barycentric_eval(t, w, vals, τ) -> value

Evaluate the interpolant of vals (matrix: nodes × components) at τ. Returns a vector over components.

barycentric_weights(t) -> w

Barycentric interpolation weights for nodes t (any distribution).

boundary_basis(params, Q, cache, θ, S; growth=30.0)
    -> (Us, Ub, E, cond_est)

Boundary data at s = S directly from the inps series (no march; used when the series is already valid at the BVP edge): Us, Ub are the small/large power-like solutions as 6-component states in the inps normalization, E the numerically-generated decaying pair, cond_est the condition number of the column-normalized basis.

build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) -> InnerAsymptoticsCache

Construct the inps Wasow asymptotic basis (GW2020 Eqs. 7–53), truncating each power series at order kmax. Evaluate with evaluate_asymptotics; pick a cutoff with pick_xmax.

cheblobatto(p) -> (t, D)

Chebyshev–Lobatto nodes on [-1, 1] in ASCENDING order and the spectral differentiation matrix D for those nodes (Trefethen's cheb, reflected).

decaying_pair(params, Q, θ, S; growth=30.0, seed=20260707, dp_res=20.0) -> E (6×2)

Orthonormal basis of the forward-decaying pair at s = S on the ray: random seeds integrated backward over [S, S+ΔS] with fully-resolved RK4, ΔS sized so the pair separates from everything else by ≈ growth e-folds.

delta_convergence(params, Q; verbose=true, refine_tol=1e-9,
                  max_rounds=16, max_cells=12000, kwargs...)
    -> (; Δ, spread, base, table)

Error bar for the matching data: solve at a strict baseline, then re-solve with each numerical knob perturbed on an independent axis (contour angle θ, element order p, series order/radius, refinement depth, march tolerance, handoff radius, purification e-folds — see variations below) and report the relative change of (Δ₁, Δ₂) per knob. θ is perturbed outward only: inward moves toward the x² ≈ −Q²(G+KF) pseudo-resonance and measures the perturbed solve, not the baseline. spread is the per-parity worst case over all knobs.

evaluate_asymptotics(cache, x::Complex; derivative=true, apply_T=true) -> (U, dU)

Complex-x version of evaluate_asymptotics (GW2020 Eq. 53; principal branch of x^r, valid for |arg x| < π). Columns of the 6×2 U are the two power-like solutions in the GW2020 Eq. (2) state convention u = (xΨ, Ξ, Υ, (xΨ)'/x, Ξ'/x, Υ'/x).

evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; derivative=true, apply_T=true) -> (U, dU)

Evaluate the inps asymptotic basis at x > 0. Returns the 6×2 complex matrix U whose two columns are the algebraically-decaying ("small") asymptotic solutions of the GGJ system, and (if derivative=true) the 6×2 matrix dU of their derivatives dU/dx.

If apply_T=false, the result is left in the J-rotated coordinate basis (used by asymptotic_residual for residual checks). The default apply_T=true returns the solutions in the original 6-component first-order-system basis used by _physical_uv and the shooting / Galerkin solvers.

evaluate_solution(res::RaySolveResult, s::Real; isol=1) -> SVector{6,ComplexF64}

Evaluate the solved state v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ') at ray parameter s (diagnostic; barycentric interpolation on the containing cell).

glasser_wang_2020_eq55(; taua=1.0, taur=1e6, v1=1.0) -> GGJParameters

D-shaped aspect-ratio-2, q = 2 surface from Glasser & Wang, Phys. Plasmas 27, 012506 (2020), Eq. 55. This is the primary benchmark case for validating the inps Wasow basis convergence (their Figs. 1–4). This is useful only for benchmarking the galerkin solver and comparing to published results.

The five coefficients below are transcribed verbatim from Eq. 55; the paper's companion operating point is the scaled growth rate Q = 1.234e-1 (their Fig. 1). Note Eq. 55 does not tabulate an inner-region matching Δ(Q) — its Δ_± (Eq. 54) is a convergence-error norm — so a quantitative Δ(Q) cross-check needs an independent inner-layer reference, not this paper alone.

Timescale parameters (taua, taur, v1) are set to canonical normalization; callers should override them for physical cases.

initial_breaks(params, Q, θ, S, p; cfl=0.4, drop=40.0, ratio=1.3,
               hmax_frac=0.08) -> Vector{Float64}

Build the initial cell-boundary vector 0 = s₀ < s₁ < … = S using the three-zone grading described in the header. cfl ≈ resolved-exponential rate per node; drop = e-folds after which decaying exponential content is considered machine-dead; ratio = geometric growth in the outer zone.

inner_Q(p::GGJParameters, γ::Number) -> ComplexF64

Scaled inner-layer growth rate Q = γ / Q₀ (GWP2016 Eq. A15).

march_boundary(params, Q, θ, S, s_m, Us, Ub; rtol=1e-9, growth=30.0, seed=20260707)
    -> (Us_m, Ub_m, F)

Transport the power-pair boundary data from the series radius S inward to s_m along the ray, modulo the decaying exponential pair (the quotient in which Δ is defined). 2-stage Radau IIA (L-stable) with Φ-budgeted purges against a carried decaying frame; explicit marchers cannot do this at large S (see file header).

mercier_di(p::GGJParameters) -> Float64

Mercier interchange index D_I = E + F + H − 1/4 (GWP2016 Eq. A9); D_I > 0 means local ideal interchange instability.

mercier_dr(p::GGJParameters) -> Float64

Resistive interchange index D_R = E + F + H² (GWP2016 Eq. A10); D_R > 0 means local resistive interchange instability.

ode_matrix([CT,] p::GGJParameters, Q, x) -> SMatrix{6,6,CT}

Coefficient matrix M(x) of dv/dx = M v, v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ'), for real or complex x. CT selects the element type (default ComplexF64; Complex{Double64} for the extended-precision march).

p1(p::GGJParameters) -> Float64

p₁ = √(−D_I), setting the large-x Frobenius exponents r± = 3/2 ± p₁ (GW2020 Eq. 49). Throws unless D_I < 0 (Mercier-stable required).

parity_rows(isol) -> SVector{3,Int}

Component indices of v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ') that vanish at s = 0 for each parity solution, mirroring deltacsetboundary:

• isol = 1 ("odd"): Ψ'(0) = 0, Ξ(0) = 0, Υ(0) = 0 → components (4, 2, 3)
• isol = 2 ("even"): Ψ(0) = 0, Ξ'(0) = 0, Υ'(0) = 0 → components (1, 5, 6)

physical_ua_dua(cache, x::Number) -> (ua, dua)

Convert the inps 6×2 basis at (possibly complex) x to physical (Ψ, Ξ, Υ) values and x-derivatives, each 3×2 (column 1 = large power solution, 2 = small). Same map as the deltac inpso_get_ua/dua convention.

pick_smax(params, Q; θ=angle(Q)/4, eps=1e-9, kmax=12, cache=nothing,
          slogmin=-1.0, slogmax=6.5, dslog=0.01) -> (S, cache, achieved)

Ray analog of pick_xmax: the smallest s at which the series residual along x = e^{iθ}s drops below eps. If never reached, returns the residual-minimizing s with achieved = false.

pick_xmax(params, Q; eps=1e-7, kmax=8, xlogmin=-1.0, xlogmax=4.0, dxlog=0.01) -> (x_max, cache)

Sweep x log-uniformly upward from 10^xlogmin and return the smallest x at which max(asymptotic_residual(cache, x)) < eps — the cutoff x_max where the GW2020 Eq. (54) convergence measure drops below tolerance (GW2020 Sec. III, Fig. 3). Also returns the InnerAsymptoticsCache it built so callers can reuse it.

Throws an ErrorException if no x in the sweep range achieves the target tolerance.

q0(p::GGJParameters) -> Float64

Growth-rate scale Q₀ = X₀ / τ_A (GWP2016 Eq. A15); physical rate = Q₀ Q.

q4_surface_benchmark() -> GGJParameters

Physical q = 4 rational-surface benchmark point (S = τR/τA ≈ 4.58×10⁶, DI ≈ −0.31166, α = √(−DI) ≈ 0.5583). Primary validation point for the rotated-ray backend on the imaginary-Q axis (pinned at Q = 100i, 500i in the test suite).

rescale_delta(Δ, p::GGJParameters) -> SVector{2,ComplexF64}

Rescale the matching data to physical Δ by S^(2√(−D_I)/3) · v₁^(2√(−D_I)) (GWP2016 Sec. IV), element-wise on (Δ_odd, Δ_even).

sfac(p::GGJParameters) -> Float64

Lundquist number S = τ_R / τ_A.

solution_profile(res::RaySolveResult; npc=8) -> (; s, x, Ψ, Ξ, Υ)

Reconstruct (Ψ, Ξ, Υ) on the solution contour, npc points per cell; each field is npts × 2 with columns (odd, even) and x = e^{iθ}s. For rotated solves the fields live on the ray — continuation back to real x is a separate (unstable) problem, by design.

solve_inner_converged(::GGJModel{:galerkin}, params::GGJParameters, γ::Number;
                      rtol=1e-2, max_levels=6, kwargs...)
    -> (; delta, converged, err, kmax, xfac, nx, nlevels)

Convergence-guarded Galerkin solve: jointly refines the three coupled accuracy knobs (series order kmax, asymptotic reach xfac, grid nx scaled to hold cells-per-unit-x ≈ cells_per_unit) until the per-component relative change of (Δ₁, Δ₂) drops below rtol, or max_levels is hit (converged=false). The metric is per real/imag component with a 5% significance floor, so a converged norm cannot mask a wrong small component.

solve_ray(params::GGJParameters, Q; θ=angle(Q)/4, kmax=12, p=12, S=nothing,
          smax_tol=1e-9, refine_tol=1e-8, max_rounds=8, max_cells=6000,
          verbose=false, kwargs...) -> RaySolveResult

Solve the GGJ inner-layer problem on the rotated ray x = e^{iθ}s. Δ follows the deltac output convention (swap + rescale_delta), directly comparable to the :galerkin backend. θ = 0 reproduces a real-axis solve; S defaults to the inps series-residual radius (pick_smax); p is the polynomial order per spectral element; remaining keywords are march/mesh knobs (see march_boundary, initial_breaks, decaying_pair).

x0(p::GGJParameters) -> Float64

Inner-layer length scale X₀ = S^(−1/3) (GWP2016 Eq. A14); physical x = X₀ X.

solve_inner(::GGJModel{:galerkin}, params::GGJParameters, γ::Number;
            kmax::Int=8, nx::Int=512, nq::Int=4, pfac::Float64=1.0,
            cutoff::Int=5, xfac::Float64=1.0, tol_res::Float64=1e-5)
            -> InnerLayerResponse

Solve the GGJ inner-layer matching problem using the Hermite-cubic finite element (Galerkin) method (GWP2016 Sec. III), in the "resonant + noexp + inps" configuration.

Returns the parity-projected matching data (GWP2016 Eqs. 34–35) with the X₀^{2√(−D_I)} physical rescaling applied, as an InnerLayerResponse whose tearing/interchange fields are the isol=1/isol=2 element solutions respectively — no parity swap (see the boundary-condition block above for the parity derivation).

solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...)
    -> InnerLayerResponse

Rotated-ray backend: converts γ via inner_Q, forwards keywords to solve_ray, and returns the parity-projected matching data. The solve_ray pair is ordered (Δ₁, Δ₂) = (interchange, tearing), so the two channels are swapped into the named fields of InnerLayerResponse. Preferred for |Q| ≳ 1 and near the imaginary axis; use solve_ray directly when the full RaySolveResult is wanted.

solve_inner(::GGJModel{:shooting}, params::GGJParameters, γ::Number;
            reltol::Float64=1e-6, abstol::Float64=1e-6,
            rtol_origin::Float64=1e-6, nps::Int=8,
            fmax::Float64=1.0, solver=Tsit5()) -> InnerLayerResponse

Solve the GGJ inner-layer matching problem by stable backward shooting in the origin-diagonalized 4×4 basis.

Returns an InnerLayerResponse(tearing, interchange) with rescaling applied. _delta_from_c0 returns a (Δ₁, Δ₂) pair where Δ₁ is the interchange (anti-symmetric / W-odd) channel and Δ₂ is the tearing (symmetric / W-even) channel. We therefore map Δ₂ → tearing and Δ₁ → interchange into the named fields, matching the physics channel labels used by the Galerkin solver and by the InnerLayerResponse docstring.

Tolerances reltol/abstol are the integrator tolerances; rtol_origin controls the truncation error of the origin Frobenius series and the choice of tmin.

solve_inner_profile(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kwargs...)
    -> (; Δ, x, Ψ, Ξ, dψdx, rescale)

Hermite-FEM implementation of the solve_inner_profile interface: real-axis solve, so Δ and the profiles come from the same solution. Same numerics/kwargs as solve_inner(GGJModel(; solver=:galerkin), ...).

solve_inner_profile(::GGJModel{:ray}, params::GGJParameters, γ::Number;
                    npc=8, certify_rtol=1e-3, kwargs...)
    -> (; Δ, x, Ψ, Ξ, dψdx, rescale, certΔ)

Rotated-ray implementation of the solve_inner_profile interface. The certified Δ comes from the optimal-contour solve at θ = arg(Q)/4 (robust for |Q| ≳ 1, where real-axis methods drift); the profiles come from a θ = 0 re-solve on the real axis, valid at physical (RPEC) |Q| since the on-axis pseudo-resonance is a regular point resolved by the BVP refinement. The θ=0-vs-optimal-θ Δ agreement doubles as a runtime certificate (two maximally different contours through the same entire solutions): the relative mismatch is returned as certΔ and warns above certify_rtol. npc sets the profile points per mesh cell; extra keywords forward to both solve_ray calls (θ is fixed by the method — do not pass it).

solve_inner_profile(params::GGJParameters, γ::Number; kwargs...)

Diagnostic variant of solve_inner(GGJModel(:galerkin), ...) that also returns the reconstructed inner-layer field profiles. Returns (Δ, Q, profile) where profile = (x, Ψ, Ξ, Υ) (Ψ = inner ξ, columns odd/even parity). Same numerics/kwargs as solve_inner; for investigating resonant-layer structure vs resistivity, not the production matching path.

LayerWidths

Resistive inner-layer length scales at one rational surface, in meters. The primary quantity is delta_s_m, the resistive layer thickness from the SLAYER del_s Riccati solve; d_beta is the β-weighted ion scale it is built from, retained as a drift-scale reference.

Fields

• ising, m, n – surface index and mode numbers (traceability)
• dels_db – dimensionless δ_s / d_β from riccati_del_s
• delta_s – complex layer thickness δ_s = dels_db · d_β [m]
• delta_s_m – |δ_s|, the resistive layer thickness in meters (primary)
• d_beta – β-weighted ion scale c_β·d_i in meters (drift reference)

delta_s_m should sit within a few orders of magnitude of d_beta for a well-posed surface (dels_db is O(1)); a large gap flags a normalisation or input problem.

SLAYERModel{S} <: InnerLayerModel

SLAYER inner-layer model selector. The type parameter S selects the Riccati formulation:

• :fitzpatrick – Pperp/Ptor Fitzpatrick formulation (default; Fitzpatrick 2023 two-fluid drift-MHD layer)

Future dispersion variants (e.g. :standard) may be added but are not currently implemented. The del_s formulation is exposed separately as the layer-thickness diagnostic slayer_layer_thickness, not as a dispersion solve_inner path.

SLAYERParameters

Dimensionless layer-physics parameters at one rational surface for the Fitzpatrick two-fluid drift-MHD SLAYER inner-layer model (Fitzpatrick 2023; Park et al. 2022), plus dimensional auxiliaries required for de-normalization. The parametrization uses P_perp, P_tor, and D_norm (not the older pr/pe/ds set).

The complex normalized growth rate Q = ω + iγ is not stored here; it is passed as a separate argument to solve_inner.

build_slayer_inputs(equil, sings, profiles; …) -> Vector{SLAYERParameters}

Build a SLAYERParameters for each rational surface in sings, pulling geometry (minor radius, r-based shear, q, dq/dψ, R₀) from the in-memory equil::PlasmaEquilibrium and kinetic data (ne, Te, Ti, ω, ω\*e, ω_*i) from profiles::KineticProfiles.

Layer inputs are assembled directly from the in-memory equilibrium and profiles, without an intermediate file round-trip.

Arguments

• equil – PlasmaEquilibrium
• sings – Vector{SingType} (one per resonant surface)
• profiles – KineticProfiles valid across all sings ψ values

Keyword arguments

• bt – toroidal field [T]. Scalar, callable of psi, or nothing (default). When nothing, the physical B_T = F(ψ) / (2π·R₀) is computed per surface from the equilibrium's F-spline. Note: equil.config.b0exp is a normalization (often just 1.0), not the physical field, so passing it as a scalar is almost always wrong.

• mu_i – ion mass in proton-mass units (default 2.0 for D).

• zeff – effective charge (default 1.0).

• chi_perp – perpendicular heat diffusivity [m²/s]. Scalar or a callable of psi (default 1.0).

• chi_tor – toroidal heat diffusivity [m²/s]. Scalar or a callable of psi (default 1.0).

• dr_val – resistive interchange index D_R = E + F + H² (Glasser-Greene-Johnson 1975) feeding the critical-Δ formulas (:lar, :rfitzp, :toroidal). When nothing (default), Julia derives it per-surface from the equilibrium as dr_val_k = D_R(ψ_k) = E_k + F_k + H_k², consistent with Connor-Hastie-Helander 2015 (PPCF 57 065001) Eq. 59 which uses (−D_R) in the χ_‖-matching critical-Δ. Pass a scalar / vector / callable to override.

  NOTE: the χ‖-matching critical-Δ requires the resistive interchange index `DR = E + F + H²(Glasser-Greene-Johnson 1975), NOT the Mercier indexDI = E + F + H − 1/4. The two differ by(H − 1/2)², which is non-trivial on shaped equilibria (~factor 3 on DIII-D); this code uses the physically correctDR`.

• dgeo_val – Connor 2015 (PPCF 57 065001) Eq. 59 geometric factor used by dc_type=:toroidal. When nothing (default), an error is raised if dc_type=:toroidal is also requested — the auto-derived formula additionally needs ⟨|∇ψ|²⟩ FSA which ResistGeometry doesn't currently expose. Pass a scalar / vector / callable to use a prescribed value. (For dc_type=:rfitzp and :lar, dgeo_val is not consulted.)

• dc_type – :none (default), :lar, :rfitzp, or :toroidal.

• theta – poloidal angle at which to measure minor radius (default 0.0, outboard midplane).

• resistivity_model – SauterNeoModel() (default), RedlNeoModel(), SpitzerModel(), or SpitzerHarmModel() (legacy Fitzpatrick σ∥). Sets the η entering τR = μ₀rs²/η. With a neoclassical model, `ftrapand ν*_e are taken from the surface'sResistGeometryif populated (viaForceFreeStates.resistevalall!), otherwise fall back to the ε-only Lin-Liu-Miller form andrs/R_0` aspect ratio.

• lnLambda_form – Coulomb-log form passed through to slayer_parameters (default :nrl; :wesson + SpitzerHarmModel() reproduces legacy SLAYER exactly).

r_based_shear(rs, q, dq_dpsi, da_dpsi) -> Float64

Convert a ψ-based shear to the r-based (Fitzpatrick) convention used throughout SLAYER:

s_r = r_s · (dq/dr) / q  =  r_s · (dq/dψ) / (q · da/dψ)

rs is the minor radius at the surface, q the safety factor, dq_dpsi the radial derivative of q with respect to ψ, and da_dpsi the derivative of the surface minor radius with respect to ψ. The two ψ derivatives must use the same ψ convention (i.e., both with respect to ψ_norm or both with respect to physical ψ — the conversion factor cancels in the ratio).

Equivalent to the conversion s_Fitz = s_psiN · r_s / (psi_N · da_dpsiN).

riccati_del_s(p::SLAYERParameters;
              q_start=5*p.D_norm, q_min=1e-5,
              reltol=1e-10, abstol=1e-10, maxiters=50_000,
              solver=Rodas5P(autodiff=false)) -> ComplexF64

Solve the SLAYER del_s inner-layer Riccati ODE and return the dimensionless layer-thickness ratio δ_s / d_β at one rational surface (see file header for references).

Integrates dW/dq = W/q − W²/q + q·E/F inward from q_start = 5·D_norm to q_min, with asymptotic boundary value W = −α·q_start² − 0.5, α = √(P̂_⊥ / (1 + 1/τ)). The result is −(π / √(1 + 1/τ)) · W′(q_min), evaluated from a single RHS call at the inner endpoint.

Returns NaN + NaN·im if the stiff integration does not converge (the caller treats this as a missing diagnostic rather than a hard error).

Multiply the result by p.d_beta to obtain the resistive layer thickness in meters (see slayer_layer_thickness).

slayer_layer_thickness(p::SLAYERParameters; kwargs...) -> LayerWidths

Compute the resistive inner-layer thickness in meters at one rational surface.

Runs riccati_del_s for the dimensionless δ_s / d_β and scales by p.d_beta to obtain δ_s in meters. Keyword arguments are forwarded to riccati_del_s.

slayer_parameters(; n_e, t_e, t_i, omega, omega_e, omega_i,
                    qval, sval_r, bt, rs, R0, mu_i, zeff,
                    chi_perp, chi_tor,
                    m, n,
                    dr_val=0.0, dgeo_val=0.0,
                    dc_type=:none, ising=0,
                    resistivity_model=SauterNeoModel(),
                    f_trap=nothing, nu_e_star=nothing,
                    R_major_eff=nothing,
                    lnLambda_form=:nrl)
    -> SLAYERParameters

Build a SLAYERParameters for one rational surface from dimensional equilibrium and kinetic-profile inputs, in the Fitzpatrick two-fluid layer parametrization (Pperp/Ptor/Dnorm; the older magnetic/electron Prandtl pr/pe and ρs-based ds parameters are not used).

Arguments

• 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
• sval_r – r-based magnetic shear r·(dq/dr)/q (Fitzpatrick). Use r_based_shear to convert from ψ-based shear.
• bt – toroidal field [T]
• rs – minor radius at the surface [m]
• R0 – major radius [m]
• mu_i – ion mass in proton-mass units (e.g. 2.0 for D)
• zeff – effective charge
• chi_perp, chi_tor – perpendicular / toroidal heat diffusivity [m²/s]
• m, n – poloidal / toroidal mode numbers at the surface
• dr_val, dgeo_val – inputs for the critical-Δ formula
• dc_type – one of :none, :lar, :rfitzp, :toroidal
• ising – singular-surface index for traceability

Resistivity kwargs

The selected resistivity sets the resistive diffusion time τ_R = μ₀ r_s² / η and hence the Lundquist number and every normalized layer parameter derived from it.

• resistivity_model – SauterNeoModel() (default: Sauter 1999 F33 trapped-particle correction, the physically-appropriate choice for H-mode tearing stability), RedlNeoModel() (improved high-ν* fit), SpitzerModel() (Sauter 18a fit, no trapped-particle correction), or SpitzerHarmModel() (Fitzpatrick Spitzer-Härm σ∥ — the legacy SLAYER closure; pair with lnLambda_form=:wesson to reproduce legacy τ_R exactly).
• f_trap – trapped-particle fraction at this surface. If not provided with a neoclassical model, falls back to Lin-Liu-Miller ε-only form with ε = rs / (R_major_eff or R0).
• nu_e_star – electron collisionality. If nothing with a neoclassical model, computed from Sauter 1999 Eq. 18b using the same ε.
• R_major_eff – ⟨R⟩ at the surface for the ν*_e formula (default R0).
• lnLambda_form – :nrl (default), :sauter, or :wesson (legacy form).

Sign convention for diamagnetic frequencies

The diamagnetic normalization uses

Q_e = -tauk · ω_*e
Q_i = -tauk · ω_*i

For the standard plasma-physics input where ω*e is tabulated negative and ω*i positive (electrons and ions drifting in opposite directions), this produces Q_e > 0, Q_i < 0, matching the opposite-drift expectation of the dispersion relation.

surface_da_dpsi(equil, psi; theta=0.0, h=1e-5) -> Float64

Central finite-difference approximation of d(minor radius)/dψ at psi. Falls back to one-sided differences near the flux-coordinate boundaries (0 or 1).

surface_minor_radius(equil, psi; theta=0.0) -> Float64

Minor radius at normalized flux psi and poloidal angle theta, computed from equil.rzphi_rsquared as √((R − R₀)² + (Z − Z₀)²). theta = 0.0 (outboard midplane) is the default; pass θ = π to measure the inboard side if you want an average.

solve_inner(::SLAYERModel{:fitzpatrick},
            p::SLAYERParameters, Q::Number;
            pmin=1e-6, p_floor=6.0,
            reltol=1e-10, abstol=1e-10,
            maxiters=50_000,
            solver=Rodas5P(autodiff=false)) -> InnerLayerResponse

Solve the Fitzpatrick SLAYER inner-layer Riccati ODE for the complex normalized growth rate Q = ω + iγ. Returns an InnerLayerResponse with tearing = Δ and interchange = 0, so the result is interface-compatible with GGJModel.solve_inner (which populates both channels); SLAYER's pressureless layer produces only the tearing channel.

Algorithm

Implements the Fitzpatrick two-fluid drift-MHD layer formulation (Fitzpatrick 2023; Park et al. 2022) for the pressureless tearing channel (no parallel flow, pe = 0). Integrates dW/dp = -(fA'/p)·W − W²/p + (fB/(fA·fC))·p³ from a large p_start (selected by _riccati_f_initial according to whether D_norm² ≷ iota_e·P_perp/P_tor^(2/3)) inward to pmin, then computes Δ = π / W'(pmin) from a single RHS evaluation at the inner endpoint.

Solver

Default Rodas5P(autodiff=false) (Rosenbrock, stiff-friendly). The analytic Jacobian wired via the ODEFunction(jac=...) field accelerates the Newton solves. AD is disabled because complex Dual propagation through the chained denominators incurs allocations in this regime; finite-difference fallback is fast enough for the 1-equation system.

Note on solver swaps: sub-percent floating-point differences between ODE solvers cascade through the outer AMR's cell-flagging decisions (ContourSearchAMR.jl::_crosses_zero) and produce structurally different AMR cell trees. Swapping the ODE solver has been observed to change the classified root/pole inventory substantially even when the most-unstable γ agrees to within ~1e-5 relative. So solver choice is not just a per-call optimization — it affects the downstream root/pole inventory. Future solver swaps need to be validated against the topology fields (n_valid_roots, n_poles), not just γ.

Keyword arguments

• pmin – inner-layer cutoff (default 1e-6)
• p_floor – floor on p_start (default 6.0)
• reltol,abstol,maxiters – stiff-solver tolerance/iteration limits
• solver – any OrdinaryDiffEq algorithm; pass Tsit5() for the non-stiff path (rarely needed here)

AMRCell

A single axis-aligned-rectangle cell of an AMR scan. The four corner Q values (q_bl, q_br, q_tl, q_tr) and corresponding residual values (d_bl, d_br, d_tl, d_tr) are sufficient for marching-squares contour extraction.

AMRResult

Output of amr_scan.

BoxActivity

Why a box was retained or skipped by multi_box_amr_scan. NoActivity means the pre-screen grid showed no zero-isoline crossings and no large-|Δ| corners; the box is excluded from refinement. The other variants record which criterion fired first.

GrowthRateResult

Output of find_growth_rates.

MultiBoxAMRResult

Output of multi_box_amr_scan. Per-box AMRResults plus the aggregated cells/Q/Δ across all active boxes. Pre-screen-inactive boxes have nothing for their AMRResult and contribute nothing to the aggregated arrays.

The aggregated (cells, Q, Δ) are suitable for direct consumption by find_growth_rates. Pre-screen evaluations are still included in Q/Δ even for skipped boxes, so any downstream pole-magnitude diagnostic that uses the flat residual list sees the full sample.

MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}}

Multi-surface dispersion data: a vector of SurfaceCoupling, the full Δ' matrix, the index of the reference surface (whose tauk defines the Q 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[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.

MultiSurfaceCouplingFull{V<:AbstractVector{<:SurfaceCoupling}}

Full 4m × 4m tearing+interchange Pletzer-Dewar dispersion matrix (Wang et al. 2020, Phys. Plasmas 27, 122503).

Given the raw 2m × 2m outer-region matrix dp_raw (side-major ordering [L_s1, R_s1, L_s2, R_s2, …], from intr.delta_prime_raw) and a vector of SurfaceCoupling (each containing the inner-layer model and parameters), calling mc(Q) assembles the 4m × 4m Pletzer-Dewar matching matrix and returns det(mat).

This is the correct Pletzer-Dewar dispersion relation for tearing+interchange coupling. A naive 2m × 2m det(D' − D(γ)) form is not equivalent: it subtracts the inner Δ (parity ± basis) from the outer D' (side-major L/R basis), mixing two different bases. The 4m × 4m matching system introduced here keeps the bases separate via the explicit C^j_{L,R} / d^j_± unknowns. For pure-tearing (pressureless SLAYER) studies use the reduced m × m MultiSurfaceCoupling instead.

Fields

• surfaces::V — per-surface SurfaceCoupling.
• 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).

ScanResult

Output of a brute-force or AMR Q-plane scan.

SurfaceCoupling{M<:InnerLayerModel, P}

Per-surface dispersion data: (model, params, dp_diag, dc, scale, tauk). Calling sc(Q) returns the complex residual

r(Q) = dp_diag - scale * solve_inner(model, params, Q).tearing - dc

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 full 2m×2m dispersion via MultiSurfaceCoupling, not this scalar form). Coupled multi-surface eigenvalues come from MultiSurfaceCoupling evaluating the determinant of the modified Δ' matrix.

amr_scan(f, Q_re_range, Q_im_range;
          nre0, nim0, passes,
          max_cells=10_000_000,
          max_cells_action=:error,
          snapshot_callback=nothing,
          parallel=Threads.nthreads() > 1) -> AMRResult

Adaptively refine a Q-plane scan of the residual f(Q). An initial nre0 × nim0 axis-aligned grid of cells is built over Q_re_range × Q_im_range and passes rounds of refinement are applied. Each pass:

1. flags any cell whose 4 corner residuals straddle zero in Re(Δ) or Im(Δ);
2. subdivides each flagged cell into 4 quadrant children, evaluating f at 5 new midpoints;
3. unflagged cells are kept unchanged.

All evaluations of f are deduplicated through a Dict{ComplexF64, ComplexF64} hash cache so that adjacent cells share a single evaluation per corner. The returned AMRResult carries both the final cell list (for marching-squares contour extraction) and the flat list of all unique Q/Δ evaluations.

Keyword arguments

• nre0, nim0 – initial coarse-grid cell counts along each axis
• passes – number of refinement passes
• max_cells – safety cap on total cells; behavior on hit is set by max_cells_action
• max_cells_action – :error (raises) or :warn_truncate (logs a warning and returns the partial result). The latter is useful for convergence-vs-resolution studies where we deliberately push max_cells and want graceful degradation. Default :error preserves the prior safety-rail behaviour.
• snapshot_callback – if not nothing, a function called after each pass (and once for the initial grid, pass=0) with arguments (pass::Int, cells::Vector{AMRCell}, cache::Dict{ComplexF64,ComplexF64}). The callback receives live references — copy if you need persistence. Used by convergence studies to extract intermediate γ at each pass count.
• parallel – evaluate f in parallel via Threads.@threads within each phase (initial grid + each refinement pass). Defaults to true when more than one Julia thread is available. Per-call evaluations of f must be thread-safe. The result is deterministic regardless of thread count: cache updates and cell-list construction stay serial, AND BLAS is pinned to one thread during the parallel region so the det()-based residual does not pick up non-reproducible multithreaded-LAPACK reductions (which otherwise flip the extracted root between near-degenerate solutions).

as_amr_result(mbres::MultiBoxAMRResult) -> AMRResult

Wrap the aggregated cells/Q/Δ from a multi-box scan as a plain AMRResult so it can be passed directly to find_growth_rates(::AMRResult, tauk; ...).

brute_force_scan(f, Q_re_range, Q_im_range; nre, nim,
                  threaded::Bool=true) -> ScanResult

Evaluate the Q-callable residual f on a regular nre × nim grid spanning the rectangle Q_re_range × Q_im_range in the complex Q plane. f must accept a single Complex argument and return a Complex value (typically a SurfaceCoupling or MultiSurfaceCoupling, but any callable works).

Use find_growth_rates(scan, tauk; ...) to extract growth-rate eigenvalues from the result.

Arguments

• f – Q-callable residual (e.g. SurfaceCoupling, MultiSurfaceCoupling)
• Q_re_range – (re_min, re_max) tuple
• Q_im_range – (im_min, im_max) tuple

Keyword arguments

• nre, nim – grid resolution along each axis
• threaded – distribute Q evaluations across Threads.@threads

find_growth_rates(amr::AMRResult, tauk::Real;
                   re_target=0.0, im_target=0.0,
                   pole_threshold=10.0,
                   filter_above_poles=true,
                   filter_outside_re=true) -> GrowthRateResult

Extract tearing growth-rate eigenvalues from an AMR AMRResult via Delaunay triangulation + marching triangles on the scattered evaluation points. The pipeline after contour extraction (segment intersection, pole classification, outside-Re filter, physical-Hz conversion) is identical to the brute-force grid path — only the contour extractor changes. Hanging-node issues from the quadtree's mixed refinement levels are resolved by the triangulation respecting every evaluated point uniformly.

find_growth_rates(scan::ScanResult, tauk::Real;
                   re_target=0.0, im_target=0.0,
                   pole_threshold=10.0,
                   filter_above_poles=true,
                   filter_outside_re=true,
                   gap_kHz_threshold=1.0) -> GrowthRateResult

Extract tearing growth-rate eigenvalues from a brute-force ScanResult by contour-intersection analysis. tauk is the per-surface time normalization used to convert Q back to physical (Hz) units (SurfaceCoupling.tauk for single-surface scans; mc.surfaces[mc.ref_idx].tauk for coupled scans).

Keyword arguments

• re_target, im_target – contour levels (zero for vanilla dispersion root-finding; nonzero values let the caller probe iso-residual contours)
• pole_threshold – intersection is classified as a pole when max(|Re(Δ)|) along the local arc of the nearest Im=0 contour exceeds this value
• filter_above_poles – discard roots whose γ exceeds the highest pole γ
• filter_outside_re – restrict the above-pole rejection to roots whose +γ step along the Im=0 contour exits the Re=0 contour loop. When true, roots that are above a pole but geometrically inside the Re=0 contour survive (matches the Python default). Note this gate fails when the Re=0 contour is OPEN (e.g., exits the Q box edge), letting spurious upper-branch roots through; the :geom and :gap flags below cover that case.
• gap_kHz_threshold – if the highest-γ root is unstable (γ > 0) AND its γ exceeds the next root by more than this many kHz, it is flagged as a :gap warning. Default 1.0 kHz.
• residual – optional dispersion-residual callable f(Q::Complex). When supplied, each contour-intersection root is POLISHED to the true zero of f by a bounded, neighbour-aware local solve (see _polish_root), making the extracted root resolution- and thread-independent. nothing (default) reports the raw marching-squares estimate, which is sensitive to grid and floating-point rounding.
• polish_maxit – max polish iterations per root (default 20). Each costs a handful of f evaluations; failures fall back to the unpolished point.
• validity_rtol – with residual supplied, a polished root is kept only if |residual| < validity_rtol · median(|Δ|) (default 1e-3). Intersections that don't polish to a true zero (contour near-misses, or surfaces whose Δ' BVP failed → |Δ'| huge and uncancellable) are dropped from valid_roots, parked in filtered_roots, and flagged :spurious. If every candidate is spurious, the result carries :no_root.

Spurious-root recursion

After the per-intersection pole / above-pole filters, the remaining roots are sorted by descending γ. The selection loop walks down this list and at each candidate evaluates two flags:

• :geom — Re(Δ)=0 contour is locally a downward-concave "hill" at the candidate (clean polyline-following quadratic fit).
• :gap — candidate is unstable AND its γ exceeds the next root's by more than gap_kHz_threshold kHz.

If BOTH fire, the candidate is discarded as spurious and the next-most- unstable root is tried. If exactly ONE fires, the candidate is accepted as primary with that warning recorded, and the next root is exposed as Q_root_secondary so downstream tools can plot or reanalyse it. If neither fires, the candidate is accepted cleanly.

multi_box_amr_scan(f, boxes;
                   pole_magnitude_threshold,
                   prescreen_nre=25, prescreen_nim=25,
                   nre0=25, nim0=25, passes=4,
                   max_cells=10_000_000,
                   max_cells_action=:error,
                   parallel=Threads.nthreads() > 1) -> MultiBoxAMRResult

Run amr_scan over multiple Q-plane boxes with a coarse pre-screen step that skips inactive boxes entirely. The typical use case is the three-stripe ω-axis scan for SLAYER coupled tearing dispersion:

ω ∈ [-75, -25],  γ ∈ [-25, +25]   (left stripe)
ω ∈ [-25, +25],  γ ∈ [-25, +25]   (centre stripe)
ω ∈ [+25, +75],  γ ∈ [-25, +25]   (right stripe)

A single 150×50 box is wasteful when the dispersion is concentrated near a narrow ω band; splitting into stripes and pre-screening lets the AMR effort land on the active stripe.

Pre-screen logic

Each box is sampled on a prescreen_nre × prescreen_nim corner grid (default 25×25, matching the typical AMR initial-grid resolution). A box is ACTIVE if ANY pre-screen cell satisfies at least one criterion:

1. sign change of Re(Δ) across the cell's 4 corners (zero-isoline of Re(Δ) crosses the cell — root candidate);
2. sign change of Im(Δ) across the cell's 4 corners (zero-isoline of Im(Δ) crosses the cell — root candidate);
3. any corner with |Δ| ≥ pole_magnitude_threshold (likely pole — the sign-only criteria miss poles whose fringe doesn't straddle a corner).

Active boxes get the full amr_scan treatment. Inactive boxes are dropped (their AMRResult is nothing).

Arguments

• f: residual function Q::ComplexF64 → Δ::ComplexF64. Must be thread-safe if parallel=true.
• boxes: vector of (Q_re_range, Q_im_range) tuples, one per box. Boxes may overlap or share boundaries; the aggregator deduplicates Q values.

Required keyword

• pole_magnitude_threshold: activity threshold for |Δ|. A natural choice is ≈ |mean(Δ)| from a baseline (or the same value used for adaptive polethreshold in `findgrowth_rates`).

Optional keywords

• prescreen_nre, prescreen_nim (default 25 each): pre-screen grid resolution. Coarser misses small features; finer wastes evaluations on inactive boxes.
• nre0, nim0, passes, max_cells, max_cells_action, parallel: forwarded to each per-box amr_scan call. Defaults match amr_scan.

Returns

A MultiBoxAMRResult. The aggregated (cells, Q, Δ) can be wrapped in an AMRResult (helper as_amr_result below) for direct use with find_growth_rates.

Notes / TODO

• Each per-box amr_scan rebuilds its own cache, so the 25×25 pre-screen corners get re-evaluated by the AMR initial pass on active boxes (≈ 676 wasted evals per active box). A future refactor could thread a shared cache through amr_scan. For now the cost is small relative to the AMR refinement evals.
• Boxes that share a boundary line (e.g. the three ω-stripe layout above) duplicate ≈ prescreen_nim+1 corner evaluations per shared edge. Also small.

Example

boxes = [((-75.0, -25.0), (-25.0, 25.0)),
    ((-25.0, 25.0), (-25.0, 25.0)),
    ((25.0, 75.0), (-25.0, 25.0))]
result = multi_box_amr_scan(f_residual, boxes;
    pole_magnitude_threshold=1e-3,
    prescreen_nre=25, prescreen_nim=25,
    nre0=25, nim0=25, passes=4)
amr = AMRResult(result.cells, result.Q, result.Δ)
roots = find_growth_rates(amr, tauk; pole_threshold=1e-3)

multi_surface_coupling(surfaces, dp_matrix;
                        ref_idx=1,
                        msing_max=min(3, length(surfaces)))
    -> MultiSurfaceCoupling

Construct a multi-surface coupling from a vector of SurfaceCoupling and the full outer-region Δ' matrix. dp_matrix must be square with side length length(surfaces) (it is the same matrix returned by PerturbedEquilibrium.SingularCoupling's Δ' boundary-value problem).

Keyword arguments

• ref_idx – index of the reference surface whose tauk defines the Q normalization. Defaults to 1.
• msing_max – number of surfaces from the front of surfaces to include in the determinant. Defaults to min(3, length(surfaces)): Δ' off-diagonal couplings beyond the third surface tend to be erratic in practice, so the determinant is conservatively truncated to the upper-left msing_max × msing_max submatrix of dp_matrix. Set explicitly (up to length(surfaces)) to override.

multi_surface_coupling_full(surfaces, dp_raw;
                            ref_idx=1,
                            msing_max=length(surfaces),
                            rotation=zeros(length(surfaces)),
                            ntor=1) -> 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.

Keyword arguments

• ref_idx — index of the reference surface whose tauk defines the Q normalization (1 ≤ ref_idx ≤ m). Defaults to 1.
• msing_max — truncate to the leading msing_max surfaces; the matching matrix becomes 4·msingmax × 4·msingmax, built from the corresponding 2·msingmax × 2·msingmax 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().

surface_coupling(model::GGJModel, params::GGJParameters,
                 dp_diag::Number) -> 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 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).

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.

surface_coupling(model::InnerLayerModel, params, dp_diag::Number;
                 dc::Real=0.0, scale::Real=1.0, tauk::Real=1.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.

surface_coupling(model::SLAYERModel, params::SLAYERParameters,
                 dp_diag::Number; dc::Real=0.0) -> SurfaceCoupling

SLAYER convenience constructor. scale is set to params.lu^(1/3) so that the dimensionless Δ from riccati_f is mapped to outer ψ-units before subtraction from the Δ' diagonal. tauk is taken from params.tauk for use by MultiSurfaceCoupling Q rescaling.

SLAYERControl

Configuration for the SLAYER tearing-mode analysis. All fields are user-facing: read from the [SLAYER] TOML section of a gpec.toml via slayer_control_from_toml, or built directly via the @kwdef keyword constructor.

Core toggles

• enabled – run the analysis at all (default false)
• inner_model – :slayer_fitzpatrick (default), :ggj_shooting, or :ggj_galerkin
• scan_mode – :amr (default) or :brute_force
• coupling_mode – :uncoupled (default, per-surface) or :coupled (multi-surface determinant)
• dc_type – critical-Δ offset selector, one of :none, :lar, :rfitzp, :toroidal (χ_‖-matching critical-Δ formulas, Connor-Hastie-Helander 2015)
• msing_max – number of surfaces to include in the coupled determinant (default 3; capped at length(sings) at runtime)

Physics knobs

• bt – toroidal field [T]. nothing (default) resolves the physical B_T = F(ψ)/(2π·R₀) per surface from the equilibrium's F-spline; a scalar or a callable of psi overrides it
• mu_i – ion mass in proton-mass units (default 2.0 for D)
• zeff – effective charge
• chi_perp, chi_tor – fallback perpendicular / toroidal heat diffusivity [m²/s], used only when the kinetic file carries no usable chi_e/chi_phi profile (dataset absent or all-zero); otherwise the file's χ⊥(ψ)/χ_φ(ψ) take precedence
• dr_val, dgeo_val – critical-Δ formula inputs. nothing (default) auto-derives them from the equilibrium: dr_val from the resistive interchange index D_R = E + F + H² at each surface, dgeo_val from the toroidal geometric factor (required only by dc_type=:toroidal). Supply a scalar only to override the auto-derivation; an explicit 0.0 disables the critical-Δ offset (Δ_crit ≡ 0)
• theta_sample – poloidal angle at which to sample minor radius (default 0.0, outboard midplane)
• resistivity_model – η closure setting τR = μ₀rs²/η: :sauter (default, neoclassical F33), :redl, :spitzer (Sauter 18a, no trapped correction), or `:spitzerharm` (legacy Fitzpatrick σ_∥)
• lnLambda_form – Coulomb-log form: :nrl (default), :sauter, or :wesson (legacy; pair with :spitzer_harm to reproduce old SLAYER)

Scan grid (used for both brute-force and AMR initial mesh)

• Q_re_range, Q_im_range – box in the normalized Q plane
• nre, nim – grid resolution along each axis

AMR refinement

• amr_passes – max refinement levels
• amr_max_cells – hard safety cap

Growth-rate-extraction filters

• pole_threshold – threshold for pole classification (default 10)
• pole_threshold_adaptive – if true, pole_threshold is OVERRIDDEN per scan with 10 × median(|Δ|) over the scan grid. The median is robust where |mean(Δ)| is not: it resists inflation from the near-pole samples that dominate a mean, and avoids the phase-cancellation that can shrink |mean|. Useful when |Δ| spans 8+ orders of magnitude (e.g. SLAYER scans where the hardcoded 10.0 default is too restrictive and classifies all intersections as poles).
• filter_above_poles – discard roots above the highest pole γ
• filter_outside_re – condition the above-pole filter on the +γ step exiting the Re(Δ)=0 contour loop

Kinetic-profile source

SLAYER reads kinetic profiles through the shared Equilibrium.read_kinetic_file reader (the same standardized object used by the kinetic/NTV physics), so there is one consistent interface for resistive and kinetic profiles.

• profile_file – path to a kinetic-profile file (relative to the run dir), required when SLAYER is enabled. HDF5 (.h5) files use the GPEC kinetic schema and may carry chi_e (χ⊥) and chi_phi (χφ); ASCII tables are also accepted but carry no χ. When a χ dataset is absent or all-zero, the scalar `chiperp/chi_tor` fallbacks below are used (so a file can keep the χ keys set to 0 to defer to the scalars).
• profile_group – group within the HDF5 file (default "/")

Output control

• store_scan – write the full Q/Δ scan grid to HDF5. false by default to keep the output file small.

SLAYERResult

Output of run_slayer. Carries both summary eigenvalues (ωHz, γHz) and full diagnostic detail (valid roots, poles, filtered roots, contours) for downstream inspection and HDF5 output.

Fields

• enabled – true only when the analysis actually ran
• control – the SLAYERControl used (frozen snapshot)
• params – Vector{SLAYERParameters}, one per surface
• rational_psi, rational_q – normalized poloidal flux ψN and safety factor q of each analyzed surface, aligned with params. Empty when the analysis was built from bare parameters (`runslayerfrominputs` without the surface list), in which case the HDF5 writer skips them.
• dp_matrix – outer-region Δ' matrix used in the analysis
• Q_root – tearing eigenvalue(s) in normalized Q
  • length nsurfaces in :uncoupled mode
  • length 1 in :coupled mode (global eigenvalue normalized by params[1].tauk)
• omega_Hz, gamma_Hz – physical rotation frequency / growth rate
• per_surface_extraction – Vector{GrowthRateResult} of length nsurfaces in uncoupled mode (each includes polelines, pole list, valid roots, filtered roots). Empty in coupled mode.
• coupled_extraction – single GrowthRateResult in coupled mode. nothing otherwise.
• layer_widths – Vector{LayerWidths}, one per surface: the resistive layer thickness (in meters) from the del_s Riccati solve plus FKR / visco-resistive sanity scales. Empty when disabled.
• scan_data – scan results (per-surface in uncoupled, single entry in coupled). Empty unless control.store_scan == true.

ggj_inner_deltas(params::AbstractVector{GGJParameters}, Q::Number;
                 solver=:galerkin) -> Vector{NamedTuple}

Evaluate both parity channels of the GGJ inner-layer matching data at the complex normalized growth rate Q for each rational surface. Returns a vector of (ising, tearing, interchange) named tuples, where tearing (GWP Δ+) is the reconnecting channel and interchange (GWP Δ−) 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.

run_slayer(equil, surfaces, delta_prime_matrix, control; dir_path="./") -> SLAYERResult

Loose-argument form of run_slayer, taking the equilibrium, the singular-surface vector and the outer-region Δ' matrix directly. Per-surface parameters are built via build_slayer_inputs; an empty or wrong-sized delta_prime_matrix falls back to a diagonal built from the sing.delta_prime stubs.

run_slayer(result, control; dir_path="./") -> SLAYERResult

Orchestrate the full SLAYER analysis against a ForceFreeStates.ForceFreeStatesResult, reading its equilibrium, singular surfaces and Δ' matrix. Kinetic profiles are read from control.profile_file (relative to dir_path) through the shared Equilibrium.read_kinetic_file reader; when the file carries chi_e/chi_phi profiles they set χ⊥(ψ)/χφ(ψ), otherwise the scalar `control.chiperp/chi_tor` fallbacks are used.

The toroidal field comes from control.bt; leaving it unset (the default) makes build_slayer_inputs evaluate the physical B_T = F(ψ)/(2π·R₀) per surface.

Returns an enabled=false SLAYERResult when control.enabled is false.

run_slayer_from_inputs(params::Vector{SLAYERParameters},
                        dp_matrix::AbstractMatrix,
                        control::SLAYERControl) -> SLAYERResult

Run the SLAYER tearing analysis given pre-built per-surface SLAYERParameters and the outer-region Δ' matrix. Bypasses the equilibrium-driven build_slayer_inputs step — use this when the parameters are already known (e.g. in unit tests or when rebuilding from cached HDF5 output).

slayer_control_from_toml(section::AbstractDict) -> SLAYERControl

Parse a [SLAYER] TOML section into a SLAYERControl. Known nested subsections ([SLAYER.scan_grid], [SLAYER.amr], [SLAYER.growth_rate_filter]) are flattened into the top-level fields. Unknown keys raise an error so typos don't silently produce defaults.

write_slayer_hdf5!(parent::Union{HDF5.File,HDF5.Group},
                    result::SLAYERResult)

Write result into a Tearing/ subgroup of parent. The subgroup is created if missing and overwritten if it already exists (keeps the output file reproducible across reruns).