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 Δ. 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 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. 48; exponents from 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 (the Wasow asymptotic basis U = TPQSY of GW2020 Eq. 53) at coordinate x to the physical (Ψ, Ξ, Υ) and (Ψ', Ξ', Υ') representation used by deltac/inpso. 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. Port of inpsogetuv. 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))

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). Mirrors inps_delta: 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).

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.

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 a Fortran rmatch/INPS run, not this paper alone.

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₀ (GWP2016 Eq. A15), the eigenvalue appearing in the inner-region equations (Eq. 11) and the inps Wasow basis. 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 (GWP2016 Eq. A9; also the quantity under the root in the GW2020 Eq. 49 power exponents). D_I > 0 indicates local ideal interchange instability.

mercier_dr(p::GGJParameters) -> Float64

Resistive interchange index D_R = E + F + H² = D_I + (H − 1/2)² (GWP2016 Eq. A10). D_R > 0 indicates local resistive interchange instability.

p1(p::GGJParameters) -> Float64

p₁ = √(−D_I), the Mercier power that sets the large-x Frobenius exponents r± = 3/2 ± √(−D_I) (GW2020 Eq. 49; μ = −1/2 ± √(−DI) in GWP2016 Eq. 26). The Mercier-stable branch requires `DI < 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 — 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. Mirrors inps_xmax.

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

q0(p::GGJParameters) -> Float64

Inner-layer growth-rate scale factor Q₀ = X₀ / τ_A (GWP2016 Eq. A15), relating the scaled growth rate Q to the physical rate s via s = Q₀ Q.

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

Map the scaled inner-layer matching data back to physical Δ at the rational surface by the X₀^(2√(−D_I)) rescaling implied by the power-like matching (GWP2016 Sec. IV; the inner solution ~ X^{r±} converts to physical x = X₀ X). Operates element-wise on a 2-vector of (Δ_odd, Δ_even).

sfac(p::GGJParameters) -> Float64

Lundquist number S = τ_R / τ_A (GWP2016 Sec. I) that sets the inner-layer length and time scales X₀ ∝ S^(−1/3), Q₀ ∝ S^(−1/3)/τ_A.

solve_inner_converged(::GGJModel{:galerkin}, params::GGJParameters, γ::Number;
                      rtol=1e-2, kmax0=12, xfac0=1.5, cells_per_unit=3.0,
                      nx_min=1024, nx_max=8192, kmax_step=2, kmax_max=28,
                      xfac_growth=1.5, max_levels=6, nq=6, pfac=1.0, cutoff=8, tol_res=1e-4)
    -> (; delta, converged, err, kmax, xfac, nx, nlevels)

Convergence-guarded GGJ inner-layer solve: only returns a Δ once it is stable under joint refinement of the three coupled accuracy knobs — series order kmax, asymptotic reach xfac, and grid resolution nx. Successively refines all three (raising kmax clears the high-|Q| series floor; raising xfac clears the reach floor; nx is scaled with xmax to hold cells-per-unit-x ≈ cells_per_unit, which prevents the grid-starvation breakup that otherwise corrupts Δ at large xfac/high |Q|) until the per-component relative change of (Δ₁, Δ₂) drops below rtol, or max_levels is hit (then converged=false).

The metric is per real/imag component with a significance floor (5% of |Δᵢ|), so a converged norm cannot mask a wrong reactive part Re(Δ₁) — the failure mode under-reach produces (Im dominates |Δ₁|).

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.

x0(p::GGJParameters) -> Float64

Inner-layer displacement scale factor X₀ = S^(−1/3) (GWP2016 Eq. A14), relating scaled X to physical x = ψ − ψ₀ via 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)
            -> SVector{2,ComplexF64}

Solve the GGJ inner-layer matching problem using the Hermite-cubic finite element (Galerkin) method (GWP2016 Sec. III). Direct port of rmatch/deltac.f 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. 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.