Inner Layer Module

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.

solve_inner(model::InnerLayerModel, params, γ::ComplexF64; kwargs...) -> SVector{2,ComplexF64}

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

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

GGJModel{S} <: InnerLayerModel

Glasser–Greene–Johnson resistive inner-layer model. The type parameter S selects the solver: :galerkin (default) for the Hermite-cubic finite element solver and :shooting for the backward stable-shoot solver. Both implementations consume the same inps asymptotic-basis kernel and return the parity-projected matching data.

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 Δ.

Fields are the same as the Fortran resist_type:

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

GGJShootingSystem

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

InnerAsymptoticsCache

Frozen state of the inps Wasow asymptotic-basis construction for a single (GGJParameters, Q) pair. All matrices are stored as SMatrix/SVector so the evaluator can run allocation-free on a hot path.

Index convention: P[k+1] holds the k-th-order matrix P_k, B[k+1] holds B_k, etc., for k = 0, 1, …, the upper bound documented in each field.

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. 49).
• Y, Z — 2×2 series matrices, k = 0..kmax (Eq. 52).

_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 at coordinate x to the physical (W, N, Θ) and (W', N', Θ') representation used by deltac/inpso. Returns (ua, dua) each 3×2 complex, where columns are the two algebraic solutions and rows are (W, N, Θ).

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

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

Compute the residual D₆(x) of the asymptotic basis at x for each of the two algebraic columns. Mirrors inps_delta (inps.f lines 713–752): 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.

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

Construct the inps Wasow asymptotic basis for the given GGJ parameters and dimensionless growth rate Q. Truncates each power series at order kmax (default 8). The returned cache can be evaluated at any x > 0 via evaluate_asymptotics and queried for an adaptive X_max via pick_xmax.

Reference: Glasser & Wang, Phys. Plasmas 27, 012506 (2020), Eqs. 7–53.

evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real;
                    derivative::Bool=true, apply_T::Bool=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 inps_delta for residual checks). The default apply_T=true returns the solutions in the original 6-component first-order-system basis used by inpso_get_uv and the shooting / Galerkin solvers.

glasser_wang_2020_eq55() -> 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.

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

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

Dimensionless scaled inner-layer growth rate Q = γ / Q₀ used by the inps Wasow basis (deltac_run line 151). The argument γ may be real or complex; the result is always complex.

mercier_di(p::GGJParameters) -> Float64

Mercier interchange index D_I = E + F + H − 1/4 (deltac_run line 143).

mercier_dr(p::GGJParameters) -> Float64

Resistive interchange index D_R = E + F + H² (deltarrun derivation, deltacrun line 142).

p1(p::GGJParameters) -> Float64

p₁ = √(−D_I). The Mercier-stable branch requires D_I < 0; this function asserts it.

pick_xmax(params::GGJParameters, Q::ComplexF64;
          eps::Float64=1e-7, kmax::Int=8,
          xlogmin::Float64=-1.0, xlogmax::Float64=4.0,
          dxlog::Float64=0.01) -> (Float64, InnerAsymptoticsCache)

Sweep x log-uniformly upward from 10^xlogmin and return the smallest x at which max(asymptotic_residual(cache, x)) < eps. Also returns the InnerAsymptoticsCache it built so callers can reuse it. Mirrors inps_xmax (inps.f lines 639–705).

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

q0(p::GGJParameters) -> Float64

Inner-layer growth-rate scale Q₀ = X₀ / τ_A (deltac_run line 150).

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

Apply the Wang 2020 X₀^(2√(−D_I)) rescaling that maps the inner-layer matching data back to physical Δ at the rational surface (deltacrun line 192). Operates element-wise on a 2-vector of `(Δodd, Δ_even)`.

sfac(p::GGJParameters) -> Float64

Lundquist-like ratio S = τ_R / τ_A used to define the inner-layer length and time scales.

x0(p::GGJParameters) -> Float64

Inner-layer length scale X₀ = S^(−1/3) (deltac_run line 149).

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)
            -> SVector{2,ComplexF64}

Solve the GGJ inner-layer matching problem using the Hermite-cubic finite element (Galerkin) method. Direct port of rmatch/deltac.f in the "resonant + noexp + inps" configuration.

Returns (Δ₁, Δ₂) with rescaling applied. The ordering matches deltac.f's output convention (swapped relative to deltar.f).

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()) -> SVector{2,ComplexF64}

Solve the GGJ inner-layer matching problem by stable backward shooting in the origin-diagonalized 4×4 basis. Direct port of the rmatch deltar.f algorithm.

Returns the parity-projected matching data (Δ₁, Δ₂) (already rescaled back to physical units via rescale_delta). Index ordering matches the Fortran deltar output.

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