From 5c5f6875374c298bc388d364847c2156fb8c0cda Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Wed, 4 Nov 2026 19:01:29 +0100 Subject: [PATCH 1/2] InnerLayer - REFACTOR - Move the GGJ Wasow, ray and Galerkin solvers into a shared ToroidalLayer library --- src/InnerLayer/{GGJ => ToroidalLayer}/Galerkin.jl | 0 src/InnerLayer/{GGJ => ToroidalLayer}/InnerAsymptotics.jl | 0 src/InnerLayer/{GGJ => ToroidalLayer}/Ray.jl | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename src/InnerLayer/{GGJ => ToroidalLayer}/Galerkin.jl (100%) rename src/InnerLayer/{GGJ => ToroidalLayer}/InnerAsymptotics.jl (100%) rename src/InnerLayer/{GGJ => ToroidalLayer}/Ray.jl (100%) diff --git a/src/InnerLayer/GGJ/Galerkin.jl b/src/InnerLayer/ToroidalLayer/Galerkin.jl similarity index 100% rename from src/InnerLayer/GGJ/Galerkin.jl rename to src/InnerLayer/ToroidalLayer/Galerkin.jl diff --git a/src/InnerLayer/GGJ/InnerAsymptotics.jl b/src/InnerLayer/ToroidalLayer/InnerAsymptotics.jl similarity index 100% rename from src/InnerLayer/GGJ/InnerAsymptotics.jl rename to src/InnerLayer/ToroidalLayer/InnerAsymptotics.jl diff --git a/src/InnerLayer/GGJ/Ray.jl b/src/InnerLayer/ToroidalLayer/Ray.jl similarity index 100% rename from src/InnerLayer/GGJ/Ray.jl rename to src/InnerLayer/ToroidalLayer/Ray.jl From 89c4807e1523b0b6f80a5bb4b43def2e1b590844 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 5 Nov 2026 15:54:08 +0100 Subject: [PATCH 2/2] InnerLayer - REFACTOR - Put the moved solvers behind a ToroidalLayer model interface with GGJ as its reference model --- docs/DOC_STANDARD.md | 2 +- .../inner_layer/make_backend_regime_map.jl | 2 +- docs/src/figures/manifest.toml | 2 +- docs/src/inner_layer.md | 8 + src/InnerLayer/GGJ/Equations.jl | 198 ++++++++++++++++++ src/InnerLayer/GGJ/GGJ.jl | 91 ++++++-- src/InnerLayer/GGJ/GGJParameters.jl | 11 +- src/InnerLayer/InnerLayer.jl | 10 +- src/InnerLayer/ToroidalLayer/Backends.jl | 91 ++++++++ src/InnerLayer/ToroidalLayer/Galerkin.jl | 137 ++---------- .../ToroidalLayer/InnerAsymptotics.jl | 162 ++------------ src/InnerLayer/ToroidalLayer/Ray.jl | 113 +++------- src/InnerLayer/ToroidalLayer/ToroidalLayer.jl | 152 ++++++++++++++ 13 files changed, 615 insertions(+), 364 deletions(-) create mode 100644 src/InnerLayer/GGJ/Equations.jl create mode 100644 src/InnerLayer/ToroidalLayer/Backends.jl create mode 100644 src/InnerLayer/ToroidalLayer/ToroidalLayer.jl diff --git a/docs/DOC_STANDARD.md b/docs/DOC_STANDARD.md index bafe45902f..6041fc5b80 100644 --- a/docs/DOC_STANDARD.md +++ b/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 diff --git a/docs/src/figures/inner_layer/make_backend_regime_map.jl b/docs/src/figures/inner_layer/make_backend_regime_map.jl index f0978ab303..516c6c6c8e 100644 --- a/docs/src/figures/inner_layer/make_backend_regime_map.jl +++ b/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"]) diff --git a/docs/src/figures/manifest.toml b/docs/src/figures/manifest.toml index fd45fad1e4..db6ebbe4ad 100644 --- a/docs/src/figures/manifest.toml +++ b/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" diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index e67b984ffd..bb98d6d803 100644 --- a/docs/src/inner_layer.md +++ b/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 diff --git a/src/InnerLayer/GGJ/Equations.jl b/src/InnerLayer/GGJ/Equations.jl new file mode 100644 index 0000000000..abd6a62549 --- /dev/null +++ b/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)) diff --git a/src/InnerLayer/GGJ/GGJ.jl b/src/InnerLayer/GGJ/GGJ.jl index 59082f6469..eacb5eaf86 100644 --- a/src/InnerLayer/GGJ/GGJ.jl +++ b/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. 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 diff --git a/src/InnerLayer/GGJ/GGJParameters.jl b/src/InnerLayer/GGJ/GGJParameters.jl index 3d5f59e8b4..9747d372a7 100644 --- a/src/InnerLayer/GGJ/GGJParameters.jl +++ b/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)) diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index a7205b6cde..48c68e3197 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -2,7 +2,8 @@ # # Resistive inner-layer models for matched-asymptotic resistive/extended MHD stability. # Provides an abstract `InnerLayerModel` interface and its implementations: the GGJ -# (Glasser–Greene–Johnson) and SLAYER submodules. +# (Glasser–Greene–Johnson) and SLAYER submodules. GGJ is built on the ToroidalLayer library of +# shared Wasow-basis, ray and Galerkin solvers. module InnerLayer @@ -13,11 +14,13 @@ using FastInterpolations: cubic_interp using ..Utilities include("InnerLayerInterface.jl") +include("ToroidalLayer/ToroidalLayer.jl") include("GGJ/GGJ.jl") include("SLAYER/SLAYER.jl") -import .GGJ: GGJModel, GGJParameters, mercier_di, mercier_dr, inner_Q, rescale_delta -import .GGJ: solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile +import .ToroidalLayer: ToroidalLayerModel, ToroidalLayerParameters, inner_Q, rescale_delta +import .ToroidalLayer: solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile +import .GGJ: GGJModel, GGJParameters, mercier_di, mercier_dr import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths @@ -35,6 +38,7 @@ model_name(::SLAYERModel{S}) where {S} = "slayer_$S" export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile, layer_response export layer_eta, layer_rho, layer_omega_E +export ToroidalLayer, ToroidalLayerModel, ToroidalLayerParameters export GGJ, GGJModel, GGJParameters export mercier_di, mercier_dr, inner_Q, rescale_delta export solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile diff --git a/src/InnerLayer/ToroidalLayer/Backends.jl b/src/InnerLayer/ToroidalLayer/Backends.jl new file mode 100644 index 0000000000..2cffce5450 --- /dev/null +++ b/src/InnerLayer/ToroidalLayer/Backends.jl @@ -0,0 +1,91 @@ +# Backends.jl +# +# The ray and Galerkin solves behind a model's `solve_inner` / `solve_inner_profile`, and the keywords each +# accepts, so a model type only forwards to them. + +""" + ray_response(params, γ::Number; kwargs...) -> InnerLayerResponse + +Rotated-ray matching data for a model's `solve_inner`: converts γ via `inner_Q`, forwards keywords to +[`solve_ray`](@ref), and returns (Δ₊, Δ₋) as the `interchange` and `tearing` fields of +[`InnerLayerResponse`](@ref). `real_profiles`, `npc` and `certify_rtol` are profile-only options +(see [`ray_profile`](@ref)) and are ignored here. +""" +function ray_response(params::ToroidalLayerParameters, γ::Number; + real_profiles=nothing, npc=nothing, certify_rtol=nothing, kwargs...) + res = solve_ray(params, inner_Q(params, γ); kwargs...) + return InnerLayerResponse(; interchange=res.Δ[1], tearing=res.Δ[2]) +end + +""" + ray_profile(params, γ::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. `Δ` and the +layer-center fields come from the solve at the model's default contour angle; 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 θ. +""" +function ray_profile(params::ToroidalLayerParameters, γ::Number; + real_profiles::Bool=false, npc::Int=8, certify_rtol::Float64=1e-3, kwargs...) + Q = inner_Q(params, γ) + rr = solve_ray(params, Q; kwargs...) # Δ and layer-center fields (optimal θ) + c = solution_profile(rr; npc=1) # first point is s = 0 + if !real_profiles + return (; Δ=rr.Δ, x=[0.0], Ψ=c.Ψ[1:1, :], Ξ=c.Ξ[1:1, :], profile_conversions(params)..., certΔ=NaN) + end + r0 = solve_ray(params, Q; θ=0.0, kwargs...) # real-axis profile solve + certΔ = maximum(abs.(r0.Δ .- rr.Δ) ./ max.(abs.(rr.Δ), 1e-300)) + certΔ > certify_rtol && @warn "Ray solver: θ=0 profile solve disagrees with certified Δ" Q certΔ + prof = solution_profile(r0; npc=npc) + prof.Ψ[1, :] .= c.Ψ[1, :] # x = 0 from the rotated solve, as without profiles + prof.Ξ[1, :] .= c.Ξ[1, :] + return (; Δ=rr.Δ, x=prof.s, Ψ=prof.Ψ, Ξ=prof.Ξ, profile_conversions(params)..., certΔ=certΔ) +end + +""" + galerkin_profile(params, γ::Number; kwargs...) -> (; Δ, x, Ψ, Ξ, dψdx, rescale) + +Hermite-FEM implementation of the [`solve_inner_profile`](@ref) interface: real-axis solve, so `Δ` and the +profiles come from the same solution. Same numerics and keywords as [`galerkin_response`](@ref). +""" +function galerkin_profile(params::ToroidalLayerParameters, γ::Number; kwargs...) + Δ, ws = _galerkin_solve(params, γ; kwargs...) + prof = _solution_profile(ws) + return (; Δ=Δ, x=prof.x, Ψ=prof.Ψ, Ξ=prof.Ξ, profile_conversions(params)...) +end + +""" + galerkin_response(params, γ::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 inner-layer matching problem with 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), 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. +""" +function galerkin_response(params::ToroidalLayerParameters, γ::Number; kwargs...) + Δ, _ = _galerkin_solve(params, γ; kwargs...) + return InnerLayerResponse(; interchange=Δ[1], tearing=Δ[2]) +end + +# 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), "...")]) + +""" + backend_options(solver::Symbol) -> Vector{Symbol} + +Keywords the shared `:ray` or `:galerkin` backend accepts, for validating a model's construction options. +""" +backend_options(solver::Symbol) = backend_options(Val(solver)) +backend_options(::Val{:ray}) = union(_kwarg_names(solve_ray), (:real_profiles, :npc, :certify_rtol)) +backend_options(::Val{:galerkin}) = _kwarg_names(_galerkin_solve) diff --git a/src/InnerLayer/ToroidalLayer/Galerkin.jl b/src/InnerLayer/ToroidalLayer/Galerkin.jl index 1fa75c2401..b2cdcc2f2a 100644 --- a/src/InnerLayer/ToroidalLayer/Galerkin.jl +++ b/src/InnerLayer/ToroidalLayer/Galerkin.jl @@ -1,78 +1,20 @@ # Galerkin.jl # -# Hermite-cubic finite-element (Galerkin) solver for the GGJ inner layer, porting the GWP2016 +# Hermite-cubic finite-element (Galerkin) solver for a toroidal inner layer, porting the GWP2016 # Sec. III singular-Galerkin method: the half-domain weak form (Eq. 32) on a packed grid # (Eq. 33), parity conditions at x = 0, and resonant/extension cells at x_max supplying the # matching data (Eqs. 34–35). Fortran deltac configuration: "resonant" + noexp + inps, -# Hermite basis, mpert = 3 (Ψ, Ξ, Υ). `_assemble_and_solve_point!` reuses the assembly on the -# rotated ray with a point far-field condition (Ray.jl). +# Hermite basis, mpert = 3 (Ψ, Ξ, Υ). The model's equations enter only through `layer_weak_form` and its +# Wasow basis. `_assemble_and_solve_point!` reuses the assembly on the rotated ray with a point far-field +# condition (Ray.jl). using LinearAlgebra using FastGaussQuadrature: gausslobatto using QuadGK: quadgk -""" - ggj_rows([CT,] params, 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) ([`_physical_uv`](@ref)); the ray march takes them as the lower half of its companion matrix -([`ode_matrix`](@ref)). -""" -ggj_rows(params::GGJParameters, Q::Number, x::Number) = ggj_rows(ComplexF64, params, Q, x) - -function ggj_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 - -""" - _physical_uv(params, 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 [`ggj_rows`](@ref) scaled by diag(1, Q², Q). -""" -function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Number) - Ũ, Ṽ = ggj_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 - # x_max and the two thin matching cells inside it, from where the series residual (GW2020 Eq. 54) # crosses three tolerances. Returns (xmax, dx1, dx2, cache). -function _xmax_3level(params::GGJParameters, Q::ComplexF64; +function _xmax_3level(params::ToroidalLayerParameters, Q::ComplexF64; kmax::Int=8, xfac::Float64=1.0, eps_vec::NTuple{3,Float64}=(1e-2, 5e-7, 1e-7)) cache = build_asymptotics(params, Q; kmax=kmax) @@ -316,7 +258,7 @@ end # I·(αᵢ' αⱼ') + V·(αᵢ αⱼ') + U·(αᵢ αⱼ), by Gauss–Lobatto quadrature. function _gauss_quad!(cell_mat::Array{ComplexF64,4}, cell::GalerkinCell, quad_nodes::Vector{Float64}, quad_weights::Vector{Float64}, - params::GGJParameters, Q::ComplexF64; θ::Float64=0.0) + params::ToroidalLayerParameters, Q::ComplexF64; θ::Float64=0.0) np_cell = cell.np nq = length(quad_nodes) x0 = (cell.x_right + cell.x_left) / 2 @@ -327,10 +269,10 @@ function _gauss_quad!(cell_mat::Array{ComplexF64,4}, cell::GalerkinCell, x = x0 + dx * quad_nodes[iq] w = dx * quad_weights[iq] if iszero(θ) - Imat, Umat, Vmat = _physical_uv(params, Q, x) + Imat, Umat, Vmat = layer_weak_form(params, Q, x) else # Rotated contour x = e^{iθ}s: d/dx = e^{-iθ} d/ds; times e^{2iθ}, I is unchanged, V gains e^{iθ}, U gains e^{2iθ}. - Imat, U0, V0 = _physical_uv(params, Q, cis(θ) * x) + Imat, U0, V0 = layer_weak_form(params, Q, cis(θ) * x) Umat, Vmat = cis(2θ) * U0, cis(θ) * V0 end pb, qb = hermite_basis(x, cell.x_left, cell.x_right) @@ -353,7 +295,7 @@ end function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64}, cell::GalerkinCell, quad_nodes::Vector{Float64}, quad_weights::Vector{Float64}, - params::GGJParameters, Q::ComplexF64, + params::ToroidalLayerParameters, Q::ComplexF64, cache::InnerAsymptoticsCache) np_cell = cell.np x0c = (cell.x_right + cell.x_left) / 2 @@ -369,7 +311,7 @@ function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64} x = x0c + dxc * quad_nodes[iq] w = dxc * quad_weights[iq] uax, duax = physical_ua_dua(cache, x) - Imat, Umat, Vmat = _physical_uv(params, Q, x) + Imat, Umat, Vmat = layer_weak_form(params, Q, x) pb, qb = hermite_basis(x, cell.x_left, cell.x_right) # Extra DOF φ = ua₁(x_r)·H₁₀ + ua₁'(x_r)·H₁₁, the large solution Hermite-blended at the right @@ -411,7 +353,7 @@ function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64} @inbounds for iq in 1:nq x = x0c + dxc * quad_nodes[iq] w = dxc * quad_weights[iq] - Imat, Umat, Vmat = _physical_uv(params, Q, x) + Imat, Umat, Vmat = layer_weak_form(params, Q, x) pb, qb = hermite_basis(x, cell.x_left, cell.x_right) if cell.etype == CT_EXT1 @@ -434,7 +376,7 @@ end # Resonant-cell scalar products ∫ ua₁ᵀ L ua_j dx (j = 1 stiffness, j = 2 RHS) with the asymptotic # solutions as test and trial functions; non-analytic, so adaptive QuadGK (GWP2016 Sec. III). -function _resonant_integral(cell::GalerkinCell, params::GGJParameters, +function _resonant_integral(cell::GalerkinCell, params::ToroidalLayerParameters, Q::ComplexF64, cache::InnerAsymptoticsCache; tol::Float64=1e-5) x_left = cell.x_left @@ -442,7 +384,7 @@ function _resonant_integral(cell::GalerkinCell, params::GGJParameters, function integrand_11(x) ua, dua = physical_ua_dua(cache, x) - Imat, Umat, Vmat = _physical_uv(params, Q, x) + Imat, Umat, Vmat = layer_weak_form(params, Q, x) ua1 = ua[:, 1] dua1 = dua[:, 1] return transpose(dua1) * Imat * dua1 + transpose(ua1) * Vmat * dua1 + transpose(ua1) * Umat * ua1 @@ -450,7 +392,7 @@ function _resonant_integral(cell::GalerkinCell, params::GGJParameters, function integrand_12(x) ua, dua = physical_ua_dua(cache, x) - Imat, Umat, Vmat = _physical_uv(params, Q, x) + Imat, Umat, Vmat = layer_weak_form(params, Q, x) dua1 = dua[:, 1] ua1 = ua[:, 1] ua2 = ua[:, 2] @@ -467,7 +409,7 @@ end # Global assembly and solve. function _assemble_and_solve!(ws::GalerkinWorkspace, - params::GGJParameters, Q::ComplexF64, + params::ToroidalLayerParameters, Q::ComplexF64, cache::InnerAsymptoticsCache; nq::Int=4, tol_res::Float64=1e-5) mpert = 3 @@ -540,7 +482,7 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, res_cell = ws.cells[ws.nx] x_lsode = res_cell.x_lsode ua_xmax, dua_xmax = physical_ua_dua(cache, x_lsode) - Imat_xmax, _, _ = _physical_uv(params, Q, x_lsode) + Imat_xmax, _, _ = layer_weak_form(params, Q, x_lsode) emap1 = res_cell.emap[1] if emap1 <= ws.ndim # Surface term: -ua_1^T I dua_1 at x_lsode (stiffness) @@ -575,7 +517,7 @@ end # Copy the base system to both parities, impose the x = 0 parity conditions, and solve the # leading n×n banded system of each parity by LAPACK banded LU. -function _parity_bcs_and_solve!(ws::GalerkinWorkspace, params::GGJParameters, Q::ComplexF64, n::Int) +function _parity_bcs_and_solve!(ws::GalerkinWorkspace, params::ToroidalLayerParameters, Q::ComplexF64, n::Int) mpert = 3 offset = ws.kl + ws.kl + 1 @@ -583,7 +525,7 @@ function _parity_bcs_and_solve!(ws::GalerkinWorkspace, params::GGJParameters, Q: ws.rhs[:, 2] .= ws.rhs[:, 1] cell1 = ws.cells[1] - Imat0, _, _ = _physical_uv(params, Q, cell1.x_left) + Imat0, _, _ = layer_weak_form(params, Q, cell1.x_left) # IBP surface term at x = 0 for ipert in 1:mpert, jpert in 1:mpert @@ -659,7 +601,7 @@ with Us, Ub, E given as (Ψ, Ξ, Υ, Ψ', Ξ', Υ') in x-derivatives; the last n replaced by (Δ, c₁, c₂). Returns the raw Δ and decaying-pair coefficients per parity (isol = 1 even/Δ₊, 2 odd/Δ₋); `ws.sol` holds the full nodal solution afterwards. """ -function _assemble_and_solve_point!(ws::GalerkinWorkspace, params::GGJParameters, Q::ComplexF64, +function _assemble_and_solve_point!(ws::GalerkinWorkspace, params::ToroidalLayerParameters, Q::ComplexF64, θ::Float64, Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}, E::Matrix{ComplexF64}; nq::Int=6) quad_nodes, quad_weights = gausslobatto(nq + 1) offset = ws.kl + ws.kl + 1 @@ -681,7 +623,7 @@ function _assemble_and_solve_point!(ws::GalerkinWorkspace, params::GGJParameters Kc[3ip+ipert, 3jp+jpert] = cell_mat[ipert, jpert, ip+1, jp+1] end # IBP boundary term at the far end, mirroring the x = 0 term of `_parity_bcs_and_solve!` (I is x-independent). - Imat, _, _ = _physical_uv(params, Q, 0.0) + Imat, _, _ = layer_weak_form(params, Q, 0.0) for ipert in 1:3, jpert in 1:3 Kc[6+ipert, 9+jpert] -= Imat[ipert, jpert] end @@ -770,7 +712,7 @@ end # Stock Galerkin solve shared by `solve_inner` and `solve_inner_profile`. The matching data are the # resonant cell's emap DOF, the coefficient of the large resonant solution (GWP2016 Eqs. 34–35), per # parity isol = 1 (even, Δ₊), 2 (odd, Δ₋); returns them rescaled, with the solved workspace. -function _galerkin_solve(params::GGJParameters, γ::Number; +function _galerkin_solve(params::ToroidalLayerParameters, γ::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) Q = inner_Q(params, γ) @@ -780,40 +722,3 @@ function _galerkin_solve(params::GGJParameters, γ::Number; emap1 = ws.cells[ws.nx].emap[1] return rescale_delta(SVector{2,ComplexF64}(ws.sol[emap1, 1], ws.sol[emap1, 2]), params), ws end - -""" - solve_inner_profile(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kwargs...) - -> (; Δ, x, Ψ, Ξ, dψdx, rescale) - -Hermite-FEM implementation of the [`solve_inner_profile`](@ref) interface: -real-axis solve, so `Δ` and the profiles come from the same solution. Same -numerics/kwargs as `solve_inner(GGJModel(; solver=:galerkin), ...)`. -""" -function solve_inner_profile(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) - Δ, ws = _galerkin_solve(params, γ; kwargs...) - prof = _solution_profile(ws) - return (; Δ=Δ, x=prof.x, Ψ=prof.Ψ, Ξ=prof.Ξ, _profile_conversions(params)...) -end - -""" - 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), 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. -""" -function solve_inner(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) - Δ, _ = _galerkin_solve(params, γ; kwargs...) - return InnerLayerResponse(; interchange=Δ[1], tearing=Δ[2]) -end diff --git a/src/InnerLayer/ToroidalLayer/InnerAsymptotics.jl b/src/InnerLayer/ToroidalLayer/InnerAsymptotics.jl index 8f30f4310d..b6ddca31c6 100644 --- a/src/InnerLayer/ToroidalLayer/InnerAsymptotics.jl +++ b/src/InnerLayer/ToroidalLayer/InnerAsymptotics.jl @@ -1,15 +1,17 @@ # InnerAsymptotics.jl # -# Wasow asymptotic basis ("inps" basis) for the GGJ inner-layer system. +# Wasow asymptotic basis ("inps" basis) of a toroidal inner-layer system. The model supplies A₀–A₂, the +# block-diagonalizing T and the splitting solve (`wasow_basis`, `split_block`); the recursion +# and evaluation below are shared. # # Reference: A. H. Glasser & Z. R. Wang, Phys. Plasmas 27, 012506 (2020) # ("GW2020"), Secs. II.A–II.G (Eqs. 1–53) and III (Eqs. 54–55). Notation map (all equation # numbers are GW2020): # # A0/A1/A2 -> Eqs. 3–5 (A = A₀+x⁻²A₁+x⁻⁴A₂ (Eq. 3); A₀, A₂ Eq. 4; A₁ Eq. 5) -# T, Tinv -> Eqs. 6–8 (Jordan basis of A₀; u = Tv, v' = xJv, λ = Q^{-1/2}) +# T, Tinv -> Eqs. 6–8 (block-diagonalizing basis of A₀; u = Tv, v' = xJv) # J0/J1/J2 -> Eqs. 9–10 (Jₖ = T^{-1} Aₖ T) -# P_k, B_k -> Eqs. 11–24 (splitting transformation V = PW; Lyapunov solve Eq. 22) +# P_k, B_k -> Eqs. 11–24 (splitting transformation V = PW; splitting solve Eq. 22) # Q_k, C_k -> Eqs. 25–39 (2×2 inner-block diagonalization W = QX) # D_k -> Eqs. 40–43 (shearing transformation X = SY) # Y0, R -> Eqs. 48–49 (lowest-order solution; Mercier Frobenius exponents r±) @@ -23,7 +25,7 @@ using StaticArrays """ InnerAsymptoticsCache -Frozen `inps` Wasow asymptotic-basis construction for one `(GGJParameters, Q)` +Frozen `inps` Wasow asymptotic-basis construction for one `(params, Q)` pair, stored as `SMatrix`/`SVector` for allocation-free evaluation. Index convention: `P[k+1]` holds the k-th-order matrix `P_k`, etc. @@ -37,8 +39,8 @@ convention: `P[k+1]` holds the k-th-order matrix `P_k`, etc. - `Qmat` — 2×2 inner-block transformation matrices, k = 0..kmax+2 (Eqs. 32–39). - `Y` — 2×2 series matrices Y_k = Y₀Z_k, k = 0..kmax (Eqs. 50–52). """ -struct InnerAsymptoticsCache - params::GGJParameters +struct InnerAsymptoticsCache{P<:ToroidalLayerParameters} + params::P Q::ComplexF64 kmax::Int R::SVector{2,Float64} @@ -49,126 +51,8 @@ struct InnerAsymptoticsCache Y::Vector{SMatrix{2,2,ComplexF64,4}} 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 - -# 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 - -# Splitting recurrence (GW2020 Eqs. 14–19): P_k, B_k for k = 0..kmax+2, each order by `_lyap_solve` of K_k. -function _split_recurrence(J::NTuple{3,SMatrix{6,6,ComplexF64}}, λ::ComplexF64, kmax::Int) +# Splitting recurrence (GW2020 Eqs. 14–19): P_k, B_k for k = 0..kmax+2, each order by the model's `split_block` of K_k. +function _split_recurrence(J::NTuple{3,SMatrix{6,6,ComplexF64}}, splitter, kmax::Int) N = kmax + 3 # store k = 0..kmax+2 P = Vector{SMatrix{6,6,ComplexF64,36}}(undef, N) B = Vector{SMatrix{6,6,ComplexF64,36}}(undef, N) @@ -193,7 +77,7 @@ function _split_recurrence(J::NTuple{3,SMatrix{6,6,ComplexF64}}, λ::ComplexF64, end Kacc = Kacc - P[l+1] * B[kml+1] end - Bk, Pk = _lyap_solve(Kacc, λ) + Bk, Pk = split_block(splitter, Kacc) B[kk] = Bk P[kk] = Pk end @@ -286,24 +170,22 @@ function _coefs(B::Vector{SMatrix{6,6,ComplexF64,36}}, end """ - build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) -> InnerAsymptoticsCache + build_asymptotics(params, 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`](@ref); pick a matching radius with [`pick_smax`](@ref). """ -function build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) - p1v = p1(params) - # Mercier Frobenius exponents r± = 3/2 ± √(−D_I) — GW2020 Eq. (49). - R = SVector{2,Float64}(p1v + 1.5, -p1v + 1.5) +function build_asymptotics(params::ToroidalLayerParameters, Q::ComplexF64; kmax::Int=8) + R = power_exponents(params) - T, J, λ = _build_tjmat(params, Q) - P, B = _split_recurrence(J, λ, kmax) + T, J, splitter = wasow_basis(params, Q) + P, B = _split_recurrence(J, splitter, kmax) Qm, Y = _coefs(B, R, kmax) return InnerAsymptoticsCache(params, Q, kmax, R, T, J, P, Qm, Y) end -build_asymptotics(params::GGJParameters, Q::Number; kmax::Int=8) = +build_asymptotics(params::ToroidalLayerParameters, Q::Number; kmax::Int=8) = build_asymptotics(params, ComplexF64(Q); kmax=kmax) # Horner sum y[i] = (Σ_{k=0..n} c[i,k+1] zᵏ)·z^rvec[i] and its z-derivative, z = 1/x² (GW2020 Eq. 14 form); @@ -410,7 +292,7 @@ end Evaluate the inps asymptotic basis at `x`, real (`x > 0`) or complex on a rotated contour (principal branch of `x^r`; the series is in 1/x², so valid for |arg x| < π/2). Returns the 6×2 complex matrix `U` whose two columns are the power-like (Mercier) solutions of -the GGJ system (column 1 large, column 2 small), and (if `derivative=true`) the +the layer system (column 1 large, column 2 small), 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 basis (used by `asymptotic_residual`). @@ -550,15 +432,15 @@ function asymptotic_residual(cache::InnerAsymptoticsCache, x::Number) end """ - pick_smax(params, Q; θ=angle(Q)/4, eps=1e-9, kmax=12, cache=nothing, + pick_smax(params, Q; θ=default_ray_angle(params, Q), eps=1e-9, kmax=12, cache=nothing, slogmin=-1.0, slogmax=6.5, dslog=0.01) -> (S, cache, achieved) Sweep `s` log-uniformly upward and return the smallest `s` at which the series residual (GW2020 Eq. 54) along `x = e^{iθ}s` drops below `eps`. If never reached, returns the residual-minimizing `s` with `achieved = false`. """ -function pick_smax(params::GGJParameters, Q::ComplexF64; - θ::Float64=angle(Q) / 4, eps::Float64=1e-9, kmax::Int=12, +function pick_smax(params::ToroidalLayerParameters, Q::ComplexF64; + θ::Float64=default_ray_angle(params, Q), eps::Float64=1e-9, kmax::Int=12, cache::Union{Nothing,InnerAsymptoticsCache}=nothing, slogmin::Float64=-1.0, slogmax::Float64=6.5, dslog::Float64=0.01) c = cache === nothing ? build_asymptotics(params, Q; kmax=kmax) : cache @@ -581,5 +463,5 @@ function pick_smax(params::GGJParameters, Q::ComplexF64; return best_s, c, false end -pick_smax(params::GGJParameters, Q::Number; kwargs...) = +pick_smax(params::ToroidalLayerParameters, Q::Number; kwargs...) = pick_smax(params, ComplexF64(Q); kwargs...) diff --git a/src/InnerLayer/ToroidalLayer/Ray.jl b/src/InnerLayer/ToroidalLayer/Ray.jl index 32136da918..8fce660c5d 100644 --- a/src/InnerLayer/ToroidalLayer/Ray.jl +++ b/src/InnerLayer/ToroidalLayer/Ray.jl @@ -1,25 +1,25 @@ # Ray.jl # -# Rotated-contour backend (`GGJModel{:ray}`): the GGJ equations continued onto x = e^{iθ}s, -# θ = arg(Q)/4. Far-field series data at the radius S are marched inward to s_m modulo the +# Rotated-contour solver: the layer equations continued onto x = e^{iθ}s, θ from the model's +# `default_ray_angle` (arg(Q)/4 for GGJ). Far-field series data at the radius S are marched inward to s_m modulo the # decaying exponential pair, then the Galerkin.jl Hermite assembly solves [0, s_m] with # v(s_m) = Us + Δ·Ub + c₁E₁ + c₂E₂. Method and rationale: docs/src/inner_layer.md. """ - ode_matrix([CT,] p::GGJParameters, Q, x) -> SMatrix{6,6,CT} + ode_matrix([CT,] p::ToroidalLayerParameters, Q, x) -> SMatrix{6,6,CT} -Coefficient matrix `M(x)` of GWP2016 Eq. 11 as `dv/dx = M v`, `v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')`, for -real or complex `x`. `M` is polynomial in `x`, so x = 0 is an ordinary point and the system -continues analytically to complex `x`. `CT` selects the element type (default `ComplexF64`; +Coefficient matrix `M(x)` of the layer equations (GWP2016 Eq. 11 for GGJ) as `dv/dx = M v`, +`v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')`, for real or complex `x`, built from the model's [`layer_rows`](@ref). `M` is +polynomial in `x`, so x = 0 is an ordinary point and the system continues analytically to complex `x`. `CT` selects the element type (default `ComplexF64`; `Complex{Double64}` for the extended-precision march). """ -ode_matrix(p::GGJParameters, Q::ComplexF64, x::Number) = +ode_matrix(p::ToroidalLayerParameters, Q::ComplexF64, x::Number) = ode_matrix(ComplexF64, p, Q, x) -function ode_matrix(::Type{CT}, p::GGJParameters, Q::Number, x::Number) where {CT<:Complex} - Ũ, Ṽ = ggj_rows(CT, p, Q, CT(x)) +function ode_matrix(::Type{CT}, p::ToroidalLayerParameters, Q::Number, x::Number) where {CT<:Complex} + Ũ, Ṽ = layer_rows(CT, p, Q, CT(x)) m = zeros(MMatrix{6,6,CT}) - # rows 1–3: (Ψ, Ξ, Υ)' = (Ψ', Ξ', Υ'); rows 4–6: GWP2016 Eq. 11 + # rows 1–3: (Ψ, Ξ, Υ)' = (Ψ', Ξ', Υ'); rows 4–6: the layer equations m[1, 4] = 1 m[2, 5] = 1 m[3, 6] = 1 @@ -31,18 +31,17 @@ end """ ray_breaks(params, Q, θ, S; cfl=0.2, min_cells=1024) -> Vector{Float64} -Cell boundaries `0 = s₀ < s₁ < … = S` for the Hermite BVP: each cell resolves the full -local WKB rate μ + ρs (fast Ohm and Υ-family pair μ = |√Q| + √(|Q|·|G+KF|), parabolic pair -ρ = Re(e^{2iθ}/√Q)) whether or not the mode is present, since unresolved dormant modes -alias into spurious discrete ones. `min_cells` floors the resolution where the local rates +Cell boundaries `0 = s₀ < s₁ < … = S` for the Hermite BVP: each cell resolves the full local WKB rate +μ + ρs, the model's [`fast_rates`](@ref) μ (for GGJ the Ohm and Υ-family pairs |√Q| + √(|Q|·|G+KF|)) and +its [`parabolic_rate`](@ref) ρ (Re(e^{2iθ}/√Q) for GGJ), whether or not the mode is present, since +unresolved dormant modes alias into spurious discrete ones. `min_cells` floors the resolution where the local rates are small (|Q| ≪ 1). """ -function ray_breaks(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Float64; +function ray_breaks(params::ToroidalLayerParameters, Q::ComplexF64, θ::Float64, S::Float64; cfl::Float64=0.2, min_cells::Int=1024) - sq = sqrt(Q) - ρ = real(cis(2θ) / sq) - ρ > 0 || throw(ArgumentError("contour angle θ=$θ gives non-decaying parabolic pair (ρ=$ρ ≤ 0); use θ ≈ arg(Q)/4")) - μ = abs(sq) + sqrt(abs(Q) * abs(params.G + params.K * params.F)) + ρ = parabolic_rate(params, Q, θ) + ρ > 0 || throw(ArgumentError("contour angle θ=$θ gives non-decaying parabolic pair (ρ=$ρ ≤ 0); use θ ≈ default_ray_angle(params, Q)")) + μ = fast_rates(params, Q) hmax = S / min_cells breaks = [0.0] while breaks[end] < S @@ -75,10 +74,9 @@ 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. """ -function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Float64; +function decaying_pair(params::ToroidalLayerParameters, Q::ComplexF64, θ::Float64, S::Float64; growth::Float64=30.0, seed::Int=20260707, dp_res::Float64=20.0) - sq = sqrt(Q) - ρ = real(cis(2θ) / sq) + ρ = parabolic_rate(params, Q, θ) ρ > 0 || throw(ArgumentError("θ=$θ gives ρ=$ρ ≤ 0; decaying pair undefined")) ph = cis(θ) # Window ΔS: the slowest backward grower separates from the power pair by `growth` e-folds. @@ -96,7 +94,7 @@ function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Flo for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * (S + ΔS)))) gmax_loc = max(gmax_loc, abs(real(λ))) end - rate = max(ρ * (S + ΔS) + abs(sq), gmax_loc) + rate = max(ρ * (S + ΔS) + ohm_rate(params, Q), gmax_loc) nsteps = clamp(ceil(Int, dp_res * rate * ΔS), 200, 2_000_000) h = -ΔS / nsteps @@ -132,7 +130,7 @@ which Δ is defined): 2-stage Radau IIA with purges against a carried decaying f `rtol` sets the step factor c = (rtol/100)^{1/3}; `ztrk` caps resolved-band steps at ztrk/g_max; `purge_budget` is the e-folds between purges; `growth`, `seed`, `dp_res` go to [`decaying_pair`](@ref). """ -function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, +function march_boundary(params::ToroidalLayerParameters, Q::ComplexF64, θ::Float64, S::Float64, s_m::Float64, Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}; rtol::Float64=1e-9, growth::Float64=30.0, seed::Int=20260707, ztrk::Float64=0.3, purge_budget::Float64=2.5, dp_res::Float64=20.0) @@ -255,7 +253,7 @@ Boundary data at s = S taken directly from the inps series, used when the BVP ed normalization, `E` the numerically generated decaying pair, `cond_est` the condition number of the column-normalized basis. """ -function boundary_basis(params::GGJParameters, Q::ComplexF64, +function boundary_basis(params::ToroidalLayerParameters, Q::ComplexF64, cache::InnerAsymptoticsCache, θ::Float64, S::Float64; growth::Float64=30.0, seed::Int=20260707, dp_res::Float64=20.0) Us, Ub = _power_states(cache, cis(θ) * S) @@ -300,10 +298,10 @@ struct RaySolveResult end """ - solve_ray(params::GGJParameters, Q; θ=angle(Q)/4, kmax=12, S=nothing, smax_tol=1e-9, + solve_ray(params, Q; θ=default_ray_angle(params, Q), kmax=12, S=nothing, smax_tol=1e-9, cfl=0.2, min_cells=1024, warn_cells=6000, nq=6, verbose=false, …) -> RaySolveResult -Solve the GGJ inner-layer problem on the rotated ray x = e^{iθ}s. `Δ` is the rescaled pair +Solve the inner-layer problem of `params`' model on the rotated ray x = e^{iθ}s. `Δ` is the rescaled pair (Δ₊, Δ₋) = (interchange, tearing), directly comparable to the `:galerkin` backend; θ = 0 is a real-axis solve. @@ -316,8 +314,8 @@ no error estimate is returned (at small |Q| Δ is good to about 1e-8–2e-7 rela `march_rtol` is [`march_boundary`](@ref)'s `rtol`; `growth`, `seed`, `ztrk`, `purge_budget`, `dp_res` go to `march_boundary`/[`decaying_pair`](@ref); `handoff` (diagnostic) supplies (Us, Ub) at S. """ -function solve_ray(params::GGJParameters, Q::ComplexF64; - θ::Float64=angle(Q) / 4, kmax::Int=12, +function solve_ray(params::ToroidalLayerParameters, Q::ComplexF64; + θ::Float64=default_ray_angle(params, Q), kmax::Int=12, S::Union{Nothing,Float64}=nothing, smax_tol::Float64=1e-9, growth::Float64=30.0, march_rtol::Float64=1e-9, sm_fac::Float64=1.0, cfl::Float64=0.2, min_cells::Int=1024, warn_cells::Int=6000, nq::Int=6, seed::Int=20260707, @@ -371,7 +369,7 @@ function solve_ray(params::GGJParameters, Q::ComplexF64; breaks = ray_breaks(params, Q, θ, Sbvp; cfl=cfl, min_cells=min_cells) length(breaks) - 1 > warn_cells && - @warn "GGJ ray: $(length(breaks) - 1) Hermite cells exceed warn_cells = $warn_cells, so the solve is slow; the grid follows cfl = $cfl, min_cells = $min_cells and the local rates at this Q" Q maxlog = + @warn "Ray solver: $(length(breaks) - 1) Hermite cells exceed warn_cells = $warn_cells, so the solve is slow; the grid follows cfl = $cfl, min_cells = $min_cells and the local rates at this Q" Q maxlog = 1 _id = (:ggj_ray_cells, warn_cells) ws = _galerkin_workspace(breaks, fill(CT_NONE, length(breaks) - 1)) Δraw, cexp = _assemble_and_solve_point!(ws, params, Q, θ, Us, Ub, E; nq=nq) @@ -383,56 +381,9 @@ function solve_ray(params::GGJParameters, Q::ComplexF64; return RaySolveResult(Δ, Δraw, cexp, Q, θ, Sv, series_ok, bc_cond, breaks, cfl, min_cells, ws) end -solve_ray(params::GGJParameters, Q::Number; kwargs...) = +solve_ray(params::ToroidalLayerParameters, Q::Number; kwargs...) = solve_ray(params, ComplexF64(Q); kwargs...) -""" - solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...) - -> InnerLayerResponse - -Rotated-ray backend: 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. -""" -function solve_inner(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; - real_profiles=nothing, npc=nothing, certify_rtol=nothing, kwargs...) - res = solve_ray(params, inner_Q(params, γ); kwargs...) - return InnerLayerResponse(; interchange=res.Δ[1], tearing=res.Δ[2]) -end - -""" - 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. `Δ` 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 θ. -""" -function solve_inner_profile(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; - real_profiles::Bool=false, npc::Int=8, certify_rtol::Float64=1e-3, kwargs...) - Q = inner_Q(params, γ) - rr = solve_ray(params, Q; kwargs...) # Δ and layer-center fields (optimal θ) - c = solution_profile(rr; npc=1) # first point is s = 0 - if !real_profiles - return (; Δ=rr.Δ, x=[0.0], Ψ=c.Ψ[1:1, :], Ξ=c.Ξ[1:1, :], _profile_conversions(params)..., certΔ=NaN) - end - r0 = solve_ray(params, Q; θ=0.0, kwargs...) # real-axis profile solve - certΔ = maximum(abs.(r0.Δ .- rr.Δ) ./ max.(abs.(rr.Δ), 1e-300)) - certΔ > certify_rtol && @warn "GGJ ray: θ=0 profile solve disagrees with certified Δ" ising = params.ising Q certΔ - prof = solution_profile(r0; npc=npc) - prof.Ψ[1, :] .= c.Ψ[1, :] # x = 0 from the rotated solve, as without profiles - prof.Ξ[1, :] .= c.Ξ[1, :] - return (; Δ=rr.Δ, x=prof.s, Ψ=prof.Ψ, Ξ=prof.Ξ, _profile_conversions(params)..., certΔ=certΔ) -end - """ solution_profile(res::RaySolveResult; npc=8) -> (; s, x, Ψ, Ξ, Υ) @@ -452,7 +403,7 @@ Evaluate the analytic representation `u_small + Δ·u_big` on the ray for `s ≥ res.S`; concatenated with [`solution_profile`](@ref) it should join the numeric solution seamlessly across the march zone. """ -function asymptotic_profile(params::GGJParameters, res::RaySolveResult, +function asymptotic_profile(params::ToroidalLayerParameters, res::RaySolveResult, srange::AbstractVector{<:Real}) cache = build_asymptotics(params, res.Q; kmax=12) n = length(srange) @@ -479,7 +430,7 @@ angle, mesh, series order and radius, march tolerance, BVP edge s_m, decaying-pa report the relative change of (Δ₊, Δ₋); `spread` is the worst case per parity. θ moves outward only, away from the pseudo-resonance x² = −Q²(G+KF). """ -function delta_convergence(params::GGJParameters, Q::ComplexF64; verbose::Bool=true, kwargs...) +function delta_convergence(params::ToroidalLayerParameters, Q::ComplexF64; verbose::Bool=true, kwargs...) base = solve_ray(params, Q; kwargs...) variations = [ ("θ → 1.2θ", (; θ=(base.θ == 0 ? 0.1 : 1.2 * base.θ))), @@ -507,5 +458,5 @@ function delta_convergence(params::GGJParameters, Q::ComplexF64; verbose::Bool=t return (; Δ=base.Δ, spread=SVector(spread), base=base, table=table) end -delta_convergence(params::GGJParameters, Q::Number; kwargs...) = +delta_convergence(params::ToroidalLayerParameters, Q::Number; kwargs...) = delta_convergence(params, ComplexF64(Q); kwargs...) diff --git a/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl b/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl new file mode 100644 index 0000000000..60ffed8361 --- /dev/null +++ b/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl @@ -0,0 +1,152 @@ +# ToroidalLayer.jl +# +# Shared library for toroidal resistive inner-layer models whose equations take the GGJ form: three fields +# (Ψ, Ξ, Υ), second order, coefficients polynomial in x, and a far field of two Mercier power-like solutions +# and four parabolic exponentials. It holds the Wasow far-field construction (GW2020), the rotated-ray +# solver and the real-axis Hermite (Galerkin) solver; a model supplies its equations through the interface +# below and forwards its `solve_inner` / `solve_inner_profile` to the backends. GGJ is the reference model. +# +# GWP2016 — A. H. Glasser, Z. R. Wang & J.-K. Park, Phys. Plasmas 23, 112506 (2016) +# GW2020 — A. H. Glasser & Z. R. Wang, Phys. Plasmas 27, 012506 (2020) + +module ToroidalLayer + +using LinearAlgebra +using StaticArrays +using Random +using Printf +using DoubleFloats: Double64 +using ...Utilities.HermiteFEM: hermite_basis, pack_grid + +import ..InnerLayerModel, ..InnerLayerResponse, ..InnerLayerParameters +import ..solve_inner_profile + +""" + ToroidalLayerModel{S} <: InnerLayerModel + +Abstract supertype of the inner-layer models built on this library; `S` is the solver backend symbol +(`:ray`, `:galerkin`, or a model's own). Concrete models forward `solve_inner` to [`ray_response`](@ref) or +[`galerkin_response`](@ref) and `solve_inner_profile` to [`ray_profile`](@ref) or [`galerkin_profile`](@ref). +""" +abstract type ToroidalLayerModel{S} <: InnerLayerModel end + +""" + ToroidalLayerParameters <: InnerLayerParameters + +Abstract supertype of the per-surface parameters of a [`ToroidalLayerModel`](@ref). The shared solvers read a +subtype only through this interface, so each model keeps whatever numbers its equations need: + + - [`layer_rows`](@ref) and [`layer_weak_form`](@ref) - the layer equations; + - [`wasow_basis`](@ref), [`split_block`](@ref) and [`power_exponents`](@ref) - the far-field Wasow construction; + - [`default_ray_angle`](@ref), [`parabolic_rate`](@ref), [`ohm_rate`](@ref) and [`fast_rates`](@ref) - the + far-field rates that set the ray contour and mesh; + - [`inner_Q`](@ref), [`rescale_delta`](@ref) and [`profile_conversions`](@ref) - the normalization between the + physical growth rate, the outer region and the layer variables. +""" +abstract type ToroidalLayerParameters <: InnerLayerParameters end + +""" + layer_rows([CT,] params, Q, x) -> (Ũ, Ṽ) + +The layer equations solved for the second derivatives, `u'' = Ũ·u + Ṽ·u'` for u = (Ψ, Ξ, Υ), as 3×3 matrices +of element type `CT` (default `ComplexF64`) at real or complex `x`. +""" +function layer_rows end +layer_rows(params::ToroidalLayerParameters, Q::Number, x::Number) = layer_rows(ComplexF64, params, Q, x) + +""" + layer_weak_form(params, Q::ComplexF64, x) -> (I, U, V) + +3×3 coefficients of `I·u'' − V·u' − U·u = 0` at real or complex `x`, the form the Hermite assembly integrates +by parts (GWP2016 Eqs. 12–15, 32); `I` must be independent of `x`. +""" +function layer_weak_form end + +""" + wasow_basis(params, Q::ComplexF64) -> (T, (J₀, J₁, J₂), splitter) + +The far-field system u' = x(A₀ + A₁/x² + A₂/x⁴)u in the state u = (xΨ, Ξ, Υ, (xΨ)'/x, Ξ'/x, Υ'/x) (GW2020 +Eqs. 2–5), as a basis `T` whose first two columns span the power-like subspace and last four the +exponential one, the rotated `Jₖ = T⁻¹AₖT` (GW2020 Eqs. 9–10), and the `splitter` passed to +[`split_block`](@ref). +""" +function wasow_basis end + +""" + split_block(splitter, K) -> (B, P) + +One order of the GW2020 Eq. (22) splitting equations: `B` the block-diagonal part of `K` and `P` its +off-diagonal blocks, for the exponential block of J₀ the `splitter` from [`wasow_basis`](@ref) encodes. +""" +function split_block end + +""" + power_exponents(params) -> SVector{2,Float64} + +Exponents (r₊, r₋) of the two power-like far-field solutions in the Wasow state, x^r for xΨ (GW2020 Eq. 49); +r₊ the large solution, r₋ the small one. +""" +function power_exponents end + +""" + inner_Q(params, γ::Number) -> ComplexF64 + +The layer's scaled growth rate Q for the physical rate `γ` in 1/s. +""" +function inner_Q end + +""" + rescale_delta(Δ, params) -> SVector{2,ComplexF64} + +Matching data (Δ₊, Δ₋) in the layer normalization, as the solvers return them, converted to the outer region's. +""" +function rescale_delta end + +""" + profile_conversions(params) -> (; dψdx, rescale) + +For `solve_inner_profile`: poloidal-flux distance per unit layer coordinate, `δψ = dψdx · x`, and the factor +taking the layer-normalized profiles to the outer δψ normalization. +""" +function profile_conversions end + +""" + default_ray_angle(params, Q::ComplexF64) -> Float64 + +Contour angle θ of the ray x = e^{iθ}s on which the exponential solutions decay monotonically. +""" +function default_ray_angle end + +""" + parabolic_rate(params, Q::ComplexF64, θ) -> Float64 + +Decay rate per unit s² of the slowest parabolic exponential pair along the ray at angle θ; must be positive. +""" +function parabolic_rate end + +""" + ohm_rate(params, Q::ComplexF64) -> Float64 + +Rate of the fast exponential pair from the resistive Ohm's-law balance. +""" +function ohm_rate end + +""" + fast_rates(params, Q::ComplexF64) -> Float64 + +Sum of the fast exponential rates every Hermite cell resolves, alongside the parabolic pair. +""" +function fast_rates end + +include("InnerAsymptotics.jl") +include("Galerkin.jl") +include("Ray.jl") +include("Backends.jl") + +export ToroidalLayerModel, ToroidalLayerParameters +export layer_rows, layer_weak_form, wasow_basis, split_block, power_exponents, default_ray_angle, parabolic_rate, ohm_rate, fast_rates +export inner_Q, rescale_delta, profile_conversions +export solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile +export ray_response, ray_profile, galerkin_response, galerkin_profile, backend_options + +end # module ToroidalLayer