InnerLayer - REFACTOR - Move the GGJ solvers into a shared ToroidalLayer library behind a layer-model interface

docs/DOC_STANDARD.md

@@ -87,7 +87,7 @@ and writes a machine-readable entry to `docs/src/figures/manifest.toml`:
 script  = "make_rotated_ray_contour.jl"
 commit  = "3a0837e3"
 date    = "2026-07-08"
-depends = ["src/InnerLayer/GGJ/Ray.jl"]
+depends = ["src/InnerLayer/ToroidalLayer/Ray.jl"]
 ```
 
 `depends` lists the repo-relative source files whose *numbers* the figure

docs/src/figures/inner_layer/make_backend_regime_map.jl

@@ -61,5 +61,5 @@ annotate!(plt, 8.0, 3e-6,
 
 save_doc_figure(plt, "inner_layer", "backend_regime_map";
     script=basename(@__FILE__),
-    depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl",
+    depends=["src/InnerLayer/ToroidalLayer/Ray.jl", "src/InnerLayer/ToroidalLayer/Galerkin.jl",
         "src/InnerLayer/GGJ/Shooting.jl"])

docs/src/figures/manifest.toml

@@ -14,5 +14,5 @@ script = "(not preserved)"
 [inner_layer.backend_regime_map]
 commit = "3a0837e3-dirty"
 date = "2026-07-08"
-depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl", "src/InnerLayer/GGJ/Shooting.jl"]
+depends = ["src/InnerLayer/ToroidalLayer/Ray.jl", "src/InnerLayer/ToroidalLayer/Galerkin.jl", "src/InnerLayer/GGJ/Shooting.jl"]
 script = "make_backend_regime_map.jl"

docs/src/inner_layer.md

@@ -158,6 +158,14 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer]
 Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.GGJ]
 ```
 
+### ToroidalLayer
+
+The shared library behind the GGJ `:ray` and `:galerkin` backends: the Wasow far-field construction, the rotated-ray solver and the Hermite assembly. A layer model supplies its equations through the interface on [`ToroidalLayerParameters`](@ref InnerLayer.ToroidalLayer.ToroidalLayerParameters).
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.ToroidalLayer]
+```
+
 ### SLAYER
 
 ```@autodocs

src/InnerLayer/GGJ/Equations.jl

@@ -0,0 +1,198 @@
+# Equations.jl
+#
+# The GGJ equations as a ToroidalLayer model: the coefficient rows and weak form of GWP2016 Eq. 11, the
+# Wasow basis of GW2020 Eqs. 3–10 with its closed-form splitting solve, and the far-field rates the ray
+# solver sizes its contour and mesh by.
+
+"""
+    layer_rows([CT,] params::GGJParameters, Q, x) -> (Ũ, Ṽ)
+
+GWP2016 Eq. 11 solved for the second derivatives, `u'' = Ũ·u + Ṽ·u'` for u = (Ψ, Ξ, Υ), at real or complex
+`x` (used as given). `CT` is the element type (default `ComplexF64`). The Galerkin weak form scales the rows by
+diag(1, Q², Q) ([`layer_weak_form`](@ref)); the ray march takes them as the lower half of its companion matrix
+([`ode_matrix`](@ref)).
+"""
+function layer_rows(::Type{CT}, params::GGJParameters, Q::Number, x::Number) where {CT<:Complex}
+    e = CT(params.E)
+    f = CT(params.F)
+    g = CT(params.G)
+    h = CT(params.H)
+    k = CT(params.K)
+    q = CT(Q)
+    q2 = q * q
+    x2 = x * x
+    # Ψ'' = Q Ψ − Q x Ξ + H Υ';  Ξ'' = −(x/Q) Ψ + (x²/Q) Ξ − ((E+F)/Q²) Υ − (H/Q²) Ψ';
+    # Υ'' = −(x/Q) Ψ − Q(G−KE) Ξ + (x²/Q + Q(G+KF)) Υ + HKQ Ψ'
+    #! format: off
+    Ũ = SMatrix{3,3,CT}(
+        q,          -x / q,             -x / q,
+        -q * x,     x2 / q,             -(g - k * e) * q,
+        0,          -(e + f) / q2,      x2 / q + (g + k * f) * q)
+    Ṽ = SMatrix{3,3,CT}(
+        0,          -h / q2,            h * k * q,
+        0,          0,                  0,
+        h,          0,                  0)
+    #! format: on
+    return Ũ, Ṽ
+end
+
+"""
+    layer_weak_form(params::GGJParameters, Q, x) -> (I, U, V)
+
+3×3 coefficients of `I·u'' − V·u' − U·u = 0` for u = (Ψ, Ξ, Υ) at real or complex `x`, with
+(I, V, U) = (A, −B, −C) of GWP2016 Eqs. 12–15: the rows of [`layer_rows`](@ref) scaled by diag(1, Q², Q).
+"""
+function layer_weak_form(params::GGJParameters, Q::ComplexF64, x::Number)
+    Ũ, Ṽ = layer_rows(params, Q, x)
+    q, q2 = Q, Q * Q
+    Imat = @SMatrix ComplexF64[1 0 0; 0 q2 0; 0 0 q]
+    #! format: off
+    U = @SMatrix ComplexF64[
+        Ũ[1, 1]         Ũ[1, 2]         0
+        Ũ[2, 1] * q2    Ũ[2, 2] * q2    Ũ[2, 3] * q2
+        Ũ[3, 1] * q     Ũ[3, 2] * q     Ũ[3, 3] * q
+    ]
+    V = @SMatrix ComplexF64[
+        0               0               Ṽ[1, 3]
+        Ṽ[2, 1] * q2    0               0
+        Ṽ[3, 1] * q     0               0
+    ]
+    #! format: on
+    return Imat, U, V
+end
+
+# T, T⁻¹, λ = Q^{-1/2} and J_i = T⁻¹A_iT from the coefficient matrices A₀–A₂ (GW2020 Eqs. 3–10).
+
+function _build_tjmat(p::GGJParameters, Q::ComplexF64)
+    e = ComplexF64(p.E)
+    f = ComplexF64(p.F)
+    h = ComplexF64(p.H)
+    g = ComplexF64(p.G)
+    k = ComplexF64(p.K)
+    q = Q
+    q2 = q * q
+    λ = 1 / sqrt(q)
+
+    # T (Eq. 7) and T⁻¹ (Eq. 8), listed row by row.
+    #! format: off
+    T = @SMatrix ComplexF64[
+        1 0 h*q q2/λ h*q -q2/λ
+        0 0 0 -1/λ 0 1/λ
+        0 0 -1/λ 0 1/λ 0
+        0 1 -h/λ -q2 h/λ -q2
+        0 0 0 1 0 1
+        0 0 1 0 1 0
+    ]
+
+    Tinv = @SMatrix ComplexF64[
+        1 q2 0 0 0 -h*q
+        0 0 -h 1 q2 0
+        0 0 -λ/2 0 0 1/2
+        0 -λ/2 0 0 1/2 0
+        0 0 λ/2 0 0 1/2
+        0 λ/2 0 0 1/2 0
+    ]
+    #! format: on
+
+    # A₀ and A₂: GW2020 Eq. (4); A₁: Eq. (5).
+    A0 = zeros(ComplexF64, 6, 6)
+    A1 = zeros(ComplexF64, 6, 6)
+    A2 = zeros(ComplexF64, 6, 6)
+
+    # A0
+    A0[4, 2] = -q
+    A0[5, 2] = 1 / q
+    A0[6, 3] = 1 / q
+    A0[1, 4] = 1
+    A0[2, 5] = 1
+    A0[3, 6] = 1
+    A0[4, 6] = h
+
+    # A1
+    A1[4, 1] = q
+    A1[5, 1] = -1 / q
+    A1[6, 1] = -1 / q
+    A1[6, 2] = -(g - k * e) * q
+    A1[5, 3] = -(e + f) / q2
+    A1[6, 3] = (g + k * f) * q
+    A1[4, 4] = 1
+    A1[5, 4] = -h / q2
+    A1[6, 4] = h * k * q
+    A1[5, 5] = -1
+    A1[6, 6] = -1
+
+    # A2
+    A2[4, 1] = -2
+    A2[5, 1] = h / q2
+    A2[6, 1] = -h * k * q
+
+    A0s = SMatrix{6,6,ComplexF64}(A0)
+    A1s = SMatrix{6,6,ComplexF64}(A1)
+    A2s = SMatrix{6,6,ComplexF64}(A2)
+
+    J0 = Tinv * A0s * T
+    J1 = Tinv * A1s * T
+    J2 = Tinv * A2s * T
+
+    return T, (J0, J1, J2), λ
+end
+
+# GW2020's closed-form split: J₀'s exponential block is diag(−λ, −λ, λ, λ) in the basis of Eq. 7.
+struct GW2020Split
+    λ::ComplexF64
+end
+
+wasow_basis(p::GGJParameters, Q::ComplexF64) = ((T, J, λ) = _build_tjmat(p, Q); (T, J, GW2020Split(λ)))
+split_block(s::GW2020Split, K::SMatrix{6,6,ComplexF64}) = _lyap_solve(K, s.λ)
+
+# Closed-form solve of the GW2020 Eq. (22) Lyapunov equations for the off-diagonal blocks of P_k,
+# using the structure of J₀'s lower-right block; B_k is the block-diagonal part of K_k (Eq. 21).
+function _lyap_solve(K::SMatrix{6,6,ComplexF64}, λ::ComplexF64)
+    Bm = zeros(ComplexF64, 6, 6)
+    Pm = zeros(ComplexF64, 6, 6)
+
+    # B is block-diagonal in the (r1, r2) split.
+    Bm[1, 1] = K[1, 1]
+    Bm[1, 2] = K[1, 2]
+    Bm[2, 1] = K[2, 1]
+    Bm[2, 2] = K[2, 2]
+    for i in 3:6, j in 3:6
+        Bm[i, j] = K[i, j]
+    end
+
+    # P[r1, r2] = P[1:2, 3:6] (top-right off-diagonal block)
+    @inbounds begin
+        Pm[2, 3] = -K[2, 3] / λ
+        Pm[2, 4] = -K[2, 4] / λ
+        Pm[2, 5] = K[2, 5] / λ
+        Pm[2, 6] = K[2, 6] / λ
+        Pm[1, 3] = -(K[1, 3] + Pm[2, 3]) / λ
+        Pm[1, 4] = -(K[1, 4] + Pm[2, 4]) / λ
+        Pm[1, 5] = (K[1, 5] + Pm[2, 5]) / λ
+        Pm[1, 6] = (K[1, 6] + Pm[2, 6]) / λ
+    end
+
+    # P[r2, r1] = P[3:6, 1:2] (bottom-left off-diagonal block)
+    @inbounds begin
+        Pm[3, 1] = K[3, 1] / λ
+        Pm[4, 1] = K[4, 1] / λ
+        Pm[5, 1] = -K[5, 1] / λ
+        Pm[6, 1] = -K[6, 1] / λ
+        Pm[3, 2] = (K[3, 2] - Pm[3, 1]) / λ
+        Pm[4, 2] = (K[4, 2] - Pm[4, 1]) / λ
+        Pm[5, 2] = -(K[5, 2] - Pm[5, 1]) / λ
+        Pm[6, 2] = -(K[6, 2] - Pm[6, 1]) / λ
+    end
+
+    return SMatrix{6,6,ComplexF64}(Bm), SMatrix{6,6,ComplexF64}(Pm)
+end
+
+# The exponent of the outer solutions, ±x²/(2√Q), is real on the ray at θ = arg(Q)/4.
+default_ray_angle(::GGJParameters, Q::ComplexF64) = angle(Q) / 4
+
+# Parabolic decay rate Re(e^{2iθ}/√Q) of the double exponential pair along the ray.
+parabolic_rate(::GGJParameters, Q::ComplexF64, θ::Real) = real(cis(2θ) / sqrt(Q))
+
+# Fast Ohm pair |√Q| and Υ-family pair √(|Q|·|G + KF|).
+ohm_rate(::GGJParameters, Q::ComplexF64) = abs(sqrt(Q))
+fast_rates(params::GGJParameters, Q::ComplexF64) = ohm_rate(params, Q) + sqrt(abs(Q) * abs(params.G + params.K * params.F))

src/InnerLayer/GGJ/GGJ.jl

@@ -1,12 +1,13 @@
 # GGJ.jl
 #
-# Glasser–Greene–Johnson resistive inner-layer model, with three solvers selected by the
-# `solver` type parameter of `GGJModel`:
+# Glasser–Greene–Johnson resistive inner-layer model, the reference ToroidalLayer model: its equations
+# and Wasow basis (Equations.jl) drive the shared ray and Galerkin solvers. Three solvers are selected by
+# the `solver` type parameter of `GGJModel`:
 #
 #   - `:shooting`  – backward shoot from X_max → 0; valid only for |Q| ≪ 1 (tests and comparisons)
 #   - `:galerkin`  – Hermite-cubic finite elements with resonant/extension matching cells
 #                    (GWP2016 Sec. III); accuracy depends on the surface and arg Q
-#   - `:ray`       – the same Hermite assembly on the rotated contour (Ray.jl), with the far-field
+#   - `:ray`       – the same Hermite assembly on the rotated contour, with the far-field
 #                    condition marched in from where the series converges; robust to |Q| ~ 500
 #
 # They return the matching data Δ± of GWP2016 Eqs. (34)–(35). Parity is that of the displacement Ξ:
@@ -23,22 +24,25 @@
 #   GW2020  — A. H. Glasser & Z. R. Wang, "Asymptotic solutions and convergence
 #             studies of the resistive inner region equations", Phys. Plasmas 27,
 #             012506 (2020). Wasow construction of the large-x asymptotic basis
-#             (Eqs. 1–55); implemented in InnerAsymptotics.jl.
+#             (Eqs. 1–55); recursion in ToroidalLayer/InnerAsymptotics.jl, GGJ matrices in Equations.jl.
 
 module GGJ
 
 using LinearAlgebra
 using StaticArrays
-using Random
-using Printf
 using DoubleFloats: Double64
-using ...Utilities.HermiteFEM: hermite_basis, pack_grid
 
-import ..InnerLayerModel, ..InnerLayerResponse, ..InnerLayerParameters
+import ..InnerLayerResponse
 import ..solve_inner, ..solve_inner_profile
+using ..ToroidalLayer
+import ..ToroidalLayer: inner_Q, rescale_delta, profile_conversions, power_exponents
+import ..ToroidalLayer: solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile
+import ..ToroidalLayer: layer_rows, layer_weak_form, wasow_basis, split_block, default_ray_angle, parabolic_rate, ohm_rate, fast_rates
+# Also reached as GGJ.<name> by callers and tests.
+import ..ToroidalLayer: ode_matrix, decaying_pair, RaySolveResult
 
 """
-    GGJModel{S,O} <: InnerLayerModel
+    GGJModel{S,O} <: ToroidalLayerModel{S}
     GGJModel(; solver=:ray, options...)
 
 Glasser–Greene–Johnson resistive inner-layer model. `S` selects the solver
@@ -52,7 +56,7 @@ take is an error here, not at the first solve.
 
   - `options::NamedTuple` - Backend keywords forwarded to every solve; call-site keywords win.
 """
-struct GGJModel{S,O<:NamedTuple} <: InnerLayerModel
+struct GGJModel{S,O<:NamedTuple} <: ToroidalLayerModel{S}
     options::O
 end
 
@@ -76,17 +80,70 @@ function solve_inner_profile(m::GGJModel{S}, params, γ::Number; kwargs...) wher
 end
 
 include("GGJParameters.jl")
-include("InnerAsymptotics.jl")
+include("Equations.jl")
 include("Reference.jl")
 include("Shooting.jl")
-include("Galerkin.jl")
-include("Ray.jl")
 
-# Each backend's option names: the keywords of the functions its solves forward them to.
-_kwarg_names(f) = unique!(Symbol[k for m in methods(f) for k in Base.kwarg_decl(m) if !endswith(string(k), "...")])
+"""
+    solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...) -> InnerLayerResponse
+
+Rotated-ray backend ([`ray_response`](@ref)): converts γ via `inner_Q`, forwards keywords to
+[`solve_ray`](@ref), and returns the parity-projected matching data (Δ₊, Δ₋) as the
+`interchange` and `tearing` fields of [`InnerLayerResponse`](@ref).
+Preferred for |Q| ≳ 1 and near the imaginary axis; use `solve_ray` directly
+when the full [`RaySolveResult`](@ref) is wanted. `real_profiles`, `npc` and `certify_rtol` are
+profile-only model options (see `solve_inner_profile`) and are ignored here.
+"""
+solve_inner(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; kwargs...) = ray_response(params, γ; kwargs...)
+
+"""
+    solve_inner_profile(::GGJModel{:ray}, params::GGJParameters, γ::Number;
+                        real_profiles=false, npc=8, certify_rtol=1e-3, kwargs...)
+        -> (; Δ, x, Ψ, Ξ, dψdx, rescale, certΔ)
+
+Rotated-ray implementation of the [`solve_inner_profile`](@ref) interface ([`ray_profile`](@ref)). `Δ` and the
+layer-center fields come from the solve at θ = arg(Q)/4; x = 0 lies on every contour, so
+Ψ(0), Ξ(0) are the real-axis values, and by default only they are returned (`x = [0.0]`,
+`certΔ = NaN`). With `real_profiles = true` the fields across the layer (x > 0) come from an
+extra θ = 0 solve (`npc` points per cell), a diagnostic; the x = 0 row stays the rotated
+solve's, and `certΔ` is the largest relative difference of the two solves' Δ₊, Δ₋, warned
+above `certify_rtol`. Other keywords forward to [`solve_ray`](@ref), except θ.
+"""
+solve_inner_profile(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; kwargs...) = ray_profile(params, γ; kwargs...)
+
+"""
+    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 ([`galerkin_response`](@ref)).
+
+Returns the parity-projected matching data (GWP2016 Eqs. 34–35), rescaled by
+[`rescale_delta`](@ref), as an `InnerLayerResponse`: `interchange` is the even-parity Δ₊
+(isol = 1), `tearing` the odd-parity Δ₋ (isol = 2).
+
+Knobs: `kmax` series order; `nx` cell count; `nq` Gauss–Lobatto order; `pfac` grid packing;
+`cutoff` extension-cell span (cells nx−cutoff−1 … nx−1); `xfac` scale on the two inner matching radii; `tol_res`
+resonant-cell quadrature tolerance.
+"""
+solve_inner(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) = galerkin_response(params, γ; kwargs...)
+
+"""
+    solve_inner_profile(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kwargs...)
+        -> (; Δ, x, Ψ, Ξ, dψdx, rescale)
+
+Hermite-FEM implementation of the [`solve_inner_profile`](@ref) interface ([`galerkin_profile`](@ref)):
+real-axis solve, so `Δ` and the profiles come from the same solution. Same
+numerics/kwargs as `solve_inner(GGJModel(; solver=:galerkin), ...)`.
+"""
+solve_inner_profile(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) = galerkin_profile(params, γ; kwargs...)
+
+# Each backend's option names: the shared backends' keywords, or the shooting solve's.
 function _ggj_options(solver::Symbol)
-    solver === :ray && return union(_kwarg_names(solve_ray), (:real_profiles, :npc, :certify_rtol))
-    solver === :galerkin && return _kwarg_names(_galerkin_solve)
+    solver in (:ray, :galerkin) && return backend_options(solver)
     solver === :shooting && return Base.kwarg_decl(which(solve_inner, (_BareGGJ{:shooting}, GGJParameters, ComplexF64)))
     throw(ArgumentError("GGJModel: unknown solver :$solver (expected :ray, :galerkin or :shooting)"))
 end

src/InnerLayer/GGJ/GGJParameters.jl

@@ -30,7 +30,7 @@ GWP2016 Eq. (A8); E, F, G, H, K enter the inner-region equations (Eq. 11).
 
 The growth rate `γ` is not stored here; it is a separate argument to `solve_inner`.
 """
-Base.@kwdef struct GGJParameters <: InnerLayerParameters
+Base.@kwdef struct GGJParameters <: ToroidalLayerParameters
     E::Float64
     F::Float64
     G::Float64
@@ -111,6 +111,9 @@ function rescale_delta(Δ::AbstractVector, p::GGJParameters)
     return SVector{2,ComplexF64}(Δ[1] * fac, Δ[2] * fac)
 end
 
-# Profile conversions for solve_inner_profile: δψ per unit inner coordinate (X = v₁·δψ/X₀) and the
-# large-solution amplitude rescale (v₁/X₀)^(1/2+p₁) to the outer δψ normalization.
-_profile_conversions(p::GGJParameters) = (; dψdx=x0(p) / p.v1, rescale=(p.v1 / x0(p))^(0.5 + p1(p)))
+# δψ per unit inner coordinate (X = v₁·δψ/X₀) and the large-solution amplitude rescale (v₁/X₀)^(1/2+p₁)
+# to the outer δψ normalization.
+profile_conversions(p::GGJParameters) = (; dψdx=x0(p) / p.v1, rescale=(p.v1 / x0(p))^(0.5 + p1(p)))
+
+# Mercier Frobenius exponents r± = 3/2 ± √(−D_I) — GW2020 Eq. (49).
+power_exponents(p::GGJParameters) = (p1v = p1(p); SVector{2,Float64}(p1v + 1.5, -p1v + 1.5))