From 5931e7228e83c518f9f9a3464aec528feb31959d Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 5 Nov 2026 16:59:05 +0100 Subject: [PATCH] InnerLayer.CCC - FEATURE - Add the Cho, Choi and Choi drift-resistive layer as a ToroidalLayer model MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CCCModel solves the drift-resistive layer equations of Cho, Choi & Choi, IEEE Trans. Plasma Sci. PS-15, 418 (1987), in the form of Kim, Cho & Choi, Phys. Fluids 31, 2659 (1988), Eqs. (1a)-(1c): the GGJ layer with the electron diamagnetic frequency Q* in Ohm's law and vorticity and the perpendicular-transport weight M⊥/M on the pressure equation. CCCParameters holds a GGJParameters plus Qstar and mperp and takes GGJ's scales, Mercier exponents and Δ normalization. Its equations, its Wasow basis (a block- diagonalizing T equal to GW2020's at Q* = 0, m = 1, with a general Sylvester splitting solve) and the rates of its two exponential families (Kim Eq. 2) run on the shared rotated-ray solver; only :ray is available. The library gains two optional hooks for models whose exponential pairs have different rates, parabolic_resolution_rate and ray_angle_admissible, whose defaults leave GGJ unchanged. --- docs/src/citations.md | 15 ++ docs/src/inner_layer.md | 19 ++ src/InnerLayer/CCC/CCC.jl | 85 ++++++++ src/InnerLayer/CCC/CCCParameters.jl | 41 ++++ src/InnerLayer/CCC/Equations.jl | 195 ++++++++++++++++++ src/InnerLayer/InnerLayer.jl | 8 +- src/InnerLayer/ToroidalLayer/Ray.jl | 22 +- src/InnerLayer/ToroidalLayer/ToroidalLayer.jl | 20 +- test/runtests_innerlayer.jl | 74 +++++++ 9 files changed, 470 insertions(+), 9 deletions(-) create mode 100644 src/InnerLayer/CCC/CCC.jl create mode 100644 src/InnerLayer/CCC/CCCParameters.jl create mode 100644 src/InnerLayer/CCC/Equations.jl diff --git a/docs/src/citations.md b/docs/src/citations.md index 014a25839..5bb8432ea 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -144,6 +144,21 @@ Provides an efficient algorithm for computing the full Δ' matrix (coupling betw --- +> G. S. Cho, E. H. Choi, and D. I. Choi, "Compressional effects on drift-resistive instabilities in general toroidal plasma," +> *IEEE Transactions on Plasma Science* **PS-15**, 418 (1987). + +Derives the toroidal resistive-layer equations with electron diamagnetic drift, ion acoustic coupling, ion polarization drift and perpendicular resistivity by generalizing the GGJ scheme; `CCCModel` in the `InnerLayer` module solves them. + +--- + +> J.-Y. Kim, G.-S. Cho, and D.-I. Choi, "Numerical study of drift-resistive modes in general toroidal geometry," +> *Physics of Fluids* **31**, 2659 (1988). +> DOI: [10.1063/1.866543](https://doi.org/10.1063/1.866543) + +States those equations as Eqs. (1a)–(1c), with the electron diamagnetic frequency Q* and the perpendicular-transport weight M⊥/M, and their large-x exponents (Eq. 2): the form `CCCModel` implements. Its tearing dispersion curves (Figs. 1b, 2b) are reproduced by `CCCModel` up to a constant normalization of Δ. + +--- + > A. H. Glasser and Z. R. Wang, "Asymptotic solutions and convergence studies of the resistive inner region equations," > *Physics of Plasmas* **27**, 012506 (2020). > DOI: [10.1063/1.5134999](https://doi.org/10.1063/1.5134999) diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index bb98d6d80..d11d693e5 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -138,6 +138,19 @@ p = GGJ.q4_surface_benchmark() # GGJParameters for the q = 4 benchmark surf Δ = solve_inner(GGJModel(), p, γ) # Δ.interchange, Δ.tearing with the default :ray backend ``` +### CCC: diamagnetic drift and perpendicular transport + +[`CCCModel`](@ref InnerLayer.CCC.CCCModel) solves the drift-resistive layer equations of [Cho, Choi and Choi (1987)](citations.md#Tearing), in the form of [Kim, Cho and Choi (1988)](citations.md#Tearing), Eqs. (1a)–(1c): the GGJ equations with the electron diamagnetic frequency ``Q_* = \omega_{*e}/Q_0`` in Ohm's law and vorticity and the perpendicular-transport weight ``M_\perp/M`` on the pressure equation. +[`CCCParameters`](@ref InnerLayer.CCC.CCCParameters) holds a `GGJParameters` plus `Qstar` and `mperp`; ``Q_* = 0``, ``M_\perp/M = 1`` give the GGJ equations, and the GGJ scales, Mercier exponents and Δ normalization carry over unchanged. +Its far field splits GGJ's double exponential pair into two families with scales ``\sigma_j`` set by Kim's Eq. (2) ([`parabolic_scales`](@ref InnerLayer.CCC.parabolic_scales)), so the asymptotic basis uses a general Sylvester splitting solve and the ray is centred between the families, ``\theta = (\arg\sigma_1 + \arg\sigma_2)/4``; an angle at which either family stops decaying between the real axis and the ray is refused. +Only the `:ray` backend is available. +When the two families decay at very different rates (large ``K Q Q_*^2``) the march resolves the faster one across the slower one's whole band, and a solve can take tens of seconds. + +```julia +p = CCCParameters(GGJ.q4_surface_benchmark(), 0.5, 1.4) # GGJ coefficients, Q*, M⊥/M +Δ = solve_inner(CCCModel(), p, 5im * GGJ.q0(p.ggj)) +``` + ## API reference ### Tearing @@ -158,6 +171,12 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer] Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.GGJ] ``` +### CCC + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.CCC] +``` + ### 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). diff --git a/src/InnerLayer/CCC/CCC.jl b/src/InnerLayer/CCC/CCC.jl new file mode 100644 index 000000000..896bf6a4b --- /dev/null +++ b/src/InnerLayer/CCC/CCC.jl @@ -0,0 +1,85 @@ +# CCC.jl +# +# Drift-resistive inner-layer model of Cho, Choi & Choi, IEEE Trans. Plasma Sci. PS-15, 418 (1987), in the form of +# Kim, Cho & Choi, Phys. Fluids 31, 2659 (1988), Eqs. (1a)–(1c): the GGJ layer with electron diamagnetic drift (Q*) +# and perpendicular transport (M⊥/M). A ToroidalLayer model on the shared rotated-ray solver; Q* = 0 and M⊥/M = 1 +# give the GGJ equations. +# +# Its far field keeps GGJ's two Mercier power-like solutions, but the four parabolic exponentials split into two +# pairs with different rates and phases (Kim Eq. 2), so the Wasow basis needs a general Sylvester splitting solve +# and the ray angle is centred between the two families. + +module CCC + +using LinearAlgebra +using StaticArrays + +import ..InnerLayerResponse +import ..solve_inner, ..solve_inner_profile +using ..ToroidalLayer +import ..ToroidalLayer: inner_Q, rescale_delta, profile_conversions, power_exponents +import ..ToroidalLayer: + layer_rows, layer_weak_form, wasow_basis, split_block, default_ray_angle, parabolic_rate, + parabolic_resolution_rate, ray_angle_admissible, ohm_rate, fast_rates +using ..GGJ: GGJParameters + +""" + CCCModel{S,O} <: ToroidalLayerModel{S} + CCCModel(; solver=:ray, options...) + +Drift-resistive inner-layer model of Cho, Choi & Choi (1987) in the form of Kim, Cho & Choi (1988), solved on +[`CCCParameters`](@ref). `S` selects the +solver backend; only `:ray` (the rotated-contour solver of [`solve_ray`](@ref)) is available. Backend keywords given +at construction, e.g. `CCCModel(; cfl=0.1)`, apply to every solve; a keyword the backend does not take is an error +here, not at the first solve. + +## Fields + + - `options::NamedTuple` - Backend keywords forwarded to every solve; call-site keywords win. +""" +struct CCCModel{S,O<:NamedTuple} <: ToroidalLayerModel{S} + options::O +end + +function CCCModel(; solver::Symbol=:ray, options...) + solver === :ray || throw(ArgumentError("CCCModel: unknown solver :$solver (only :ray is available)")) + allowed = backend_options(solver) + unknown = setdiff(keys(options), allowed) + isempty(unknown) || throw(ArgumentError("CCCModel(solver=:$solver) does not take $(join(unknown, ", ")); its options are $(join(allowed, ", "))")) + return CCCModel{solver,typeof(values(options))}(values(options)) +end +CCCModel{S}() where {S} = CCCModel{S,@NamedTuple{}}((;)) + +# The backends implement the option-free model; a model carrying options forwards them. +const _BareCCC{S} = CCCModel{S,@NamedTuple{}} +function solve_inner(m::CCCModel{S}, params, γ::Number; kwargs...) where {S} + m isa _BareCCC && throw(MethodError(solve_inner, (m, params, γ))) + return solve_inner(CCCModel{S}(), params, γ; m.options..., kwargs...) +end +function solve_inner_profile(m::CCCModel{S}, params, γ::Number; kwargs...) where {S} + m isa _BareCCC && throw(MethodError(solve_inner_profile, (m, params, γ))) + return solve_inner_profile(CCCModel{S}(), params, γ; m.options..., kwargs...) +end + +include("CCCParameters.jl") +include("Equations.jl") + +""" + solve_inner(::CCCModel{:ray}, params::CCCParameters, γ::Number; kwargs...) -> InnerLayerResponse + +Rotated-ray matching data (Δ₊, Δ₋) as the `interchange` and `tearing` fields of `InnerLayerResponse` +([`ray_response`](@ref)); keywords forward to [`solve_ray`](@ref). `γ = s` is the growth rate in the surface's +E×B frame, in 1/s. +""" +solve_inner(::_BareCCC{:ray}, params::CCCParameters, γ::Number; kwargs...) = ray_response(params, γ; kwargs...) + +""" + solve_inner_profile(::CCCModel{:ray}, params::CCCParameters, γ::Number; kwargs...) + +Rotated-ray implementation of the `solve_inner_profile` interface ([`ray_profile`](@ref)). +""" +solve_inner_profile(::_BareCCC{:ray}, params::CCCParameters, γ::Number; kwargs...) = ray_profile(params, γ; kwargs...) + +export CCCModel, CCCParameters, parabolic_scales + +end # module CCC diff --git a/src/InnerLayer/CCC/CCCParameters.jl b/src/InnerLayer/CCC/CCCParameters.jl new file mode 100644 index 000000000..f035ea5f7 --- /dev/null +++ b/src/InnerLayer/CCC/CCCParameters.jl @@ -0,0 +1,41 @@ +# CCCParameters.jl +# +# Per-surface parameters of the CCC layer: the GGJ coefficients and local times, which set its normalization +# exactly as for GGJ, plus the electron diamagnetic frequency and the perpendicular-transport weight. + +""" + CCCParameters + +Parameters of the drift-resistive layer of Cho, Choi & Choi, IEEE Trans. Plasma Sci. PS-15, 418 (1987), in the form +of Kim, Cho & Choi, Phys. Fluids 31, 2659 (1988), Eqs. (1a)–(1c), at one rational surface. Its Γ is Υ here. Q* = 0 and M⊥/M = 1 give the GGJ equations. + +## Fields + + - `ggj::GGJParameters` - GGJ coefficients E, F, G, H, K, M and local times; Kim's coefficients are GGJ's, and + so are the time and length scales, the Mercier exponents of the far field and the Δ normalization. + - `Qstar::Float64` - Electron diamagnetic frequency scaled like Q, Q* = ω*ₑ/Q₀, in the sign convention of the + growth rate s = γ − iω (Kim's Ohm's law carries Q + iQ*). + - `mperp::Float64` - Perpendicular-transport weight M⊥/M of the QΥ'' term; must be positive. +""" +struct CCCParameters <: ToroidalLayerParameters + ggj::GGJParameters + Qstar::Float64 + mperp::Float64 + function CCCParameters(ggj::GGJParameters, Qstar::Real, mperp::Real) + mperp > 0 || throw(ArgumentError("CCCParameters: mperp = $mperp must be positive")) + return new(ggj, Float64(Qstar), Float64(mperp)) + end +end + +""" + CCCParameters(; Qstar=0.0, mperp=1.0, ggj_fields...) + +Keyword constructor: `ggj_fields` are the `GGJParameters` keywords (`E`, `F`, `G`, `H`, `K`, `M`, `taua`, `taur`, +`v1`, `ising`). +""" +CCCParameters(; Qstar::Real=0.0, mperp::Real=1.0, ggj_fields...) = CCCParameters(GGJParameters(; ggj_fields...), Qstar, mperp) + +inner_Q(p::CCCParameters, γ::Number) = inner_Q(p.ggj, γ) +rescale_delta(Δ::AbstractVector, p::CCCParameters) = rescale_delta(Δ, p.ggj) +profile_conversions(p::CCCParameters) = profile_conversions(p.ggj) +power_exponents(p::CCCParameters) = power_exponents(p.ggj) diff --git a/src/InnerLayer/CCC/Equations.jl b/src/InnerLayer/CCC/Equations.jl new file mode 100644 index 000000000..9ecad864c --- /dev/null +++ b/src/InnerLayer/CCC/Equations.jl @@ -0,0 +1,195 @@ +# Equations.jl +# +# The CCC layer as a ToroidalLayer model: Kim, Cho & Choi (1988) Eqs. (1a)–(1c) as coefficient rows and weak form, +# their Wasow basis with a general Sylvester splitting solve, and the rates of the two exponential families of +# their far field (Kim Eq. 2), which the ray solver sizes its contour and mesh by. + +""" + layer_rows([CT,] params::CCCParameters, Q, x) -> (Ũ, Ṽ) + +Kim, Cho & Choi (1988) Eqs. (1a)–(1c) solved for the second derivatives, `u'' = Ũ·u + Ṽ·u'` for u = (Ψ, Ξ, Υ): +GWP2016 Eq. 11 with the electron diamagnetic terms iQ*(Ψ − xΥ) in Ohm's law and vorticity, the weight M⊥/M on QΥ'' +and the coupling iQ*KQ³Ξ'' in the pressure equation, which enters the Υ'' row through Ξ''. +""" +function layer_rows(::Type{CT}, params::CCCParameters, Q::Number, x::Number) where {CT<:Complex} + g = params.ggj + e = CT(g.E) + f = CT(g.F) + gg = CT(g.G) + h = CT(g.H) + k = CT(g.K) + iqs = CT(0, params.Qstar) + m = real(CT)(params.mperp) + q = CT(Q) + q2 = q * q + x2 = x * x + # Ψ'' = (Q + iQ*) Ψ − Q x Ξ − iQ* x Υ + H Υ'; + # Ξ'' = −(x/Q)(1 + iQ*/Q) Ψ + (x²/Q) Ξ + (iQ* x² − E − F)/Q² Υ − (H/Q²) Ψ'; + # Υ'' = [−x Ψ − Q²(G−KE) Ξ + (x² + Q²(G+KF)) Υ + HKQ² Ψ']/(mQ) − c Ξ'', with c = iQ*KQ²/m + ξ1 = (-x / q) * (1 + iqs / q) + ξ2 = x2 / q + ξ3 = -(e + f) / q2 + iqs * x2 / q2 + ξ4 = -h / q2 + c = iqs * k * q2 / m + #! format: off + Ũ = SMatrix{3,3,CT}( + q + iqs, ξ1, -x / (m * q) - c * ξ1, + -q * x, ξ2, -(gg - k * e) * q / m - c * ξ2, + -iqs * x, ξ3, x2 / (m * q) + (gg + k * f) * q / m - c * ξ3) + Ṽ = SMatrix{3,3,CT}( + 0, ξ4, h * k * q / m - c * ξ4, + 0, 0, 0, + h, 0, 0) + #! format: on + return Ũ, Ṽ +end + +""" + layer_weak_form(params::CCCParameters, Q, x) -> (I, U, V) + +Kim (1988) Eqs. (1a)–(1c) as `I·u'' − V·u' − U·u = 0`: the rows of [`layer_rows`](@ref) multiplied by the u'' +coefficient matrix I = [1 0 0; 0 Q² 0; 0 iQ*KQ³ (M⊥/M)Q], which is lower triangular and independent of x. +""" +function layer_weak_form(params::CCCParameters, Q::ComplexF64, x::Number) + Ũ, Ṽ = layer_rows(params, Q, x) + a = ComplexF64(0, params.Qstar) * params.ggj.K * Q^3 # iQ*KQ³ Ξ'' in the Υ row + mq = params.mperp * Q + Imat = @SMatrix ComplexF64[1 0 0; 0 Q*Q 0; 0 a mq] + return Imat, Imat * Ũ, Imat * Ṽ +end + +# T, T⁻¹ and J_i = T⁻¹A_iT from A₀–A₂ (GW2020 Eqs. 3–10) with Kim's Q* and m = M⊥/M terms. T's columns 3–6 span the +# exponential subspace of A₀ in a basis regular in Q* and m that equals GW2020 Eq. 7 at Q* = 0, m = 1, so J₀ is +# block-diagonal everywhere and its exponential block is diagonal only at the GGJ point. +function _build_tjmat(p::CCCParameters, Q::ComplexF64) + g = p.ggj + e = ComplexF64(g.E) + f = ComplexF64(g.F) + h = ComplexF64(g.H) + gg = ComplexF64(g.G) + k = ComplexF64(g.K) + iqs = ComplexF64(0, p.Qstar) + qs2 = p.Qstar^2 + m = p.mperp + q = Q + q2 = q * q + λ = 1 / sqrt(q) + hqm = h * q * m + hkq3 = iqs * k * h * q * q2 + + #! format: off + T = @SMatrix ComplexF64[ + 1 0 hqm q2/λ+hkq3 hqm hkq3-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 -hkq3 -hqm + 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); Q* and m terms from Kim (1988) Eqs. (1a)–(1c). + A0 = zeros(ComplexF64, 6, 6) + A1 = zeros(ComplexF64, 6, 6) + A2 = zeros(ComplexF64, 6, 6) + + A0[4, 2] = -q + A0[4, 3] = -iqs + A0[5, 2] = 1 / q + A0[5, 3] = iqs / q2 + A0[6, 2] = -iqs * k * q / m + A0[6, 3] = 1 / (m * q) + k * qs2 / m + A0[1, 4] = 1 + A0[2, 5] = 1 + A0[3, 6] = 1 + A0[4, 6] = h + + A1[4, 1] = q + iqs + A1[5, 1] = -1 / q - iqs / q2 + A1[6, 1] = -1 / (m * q) + iqs * k * q / m - k * qs2 / m + A1[6, 2] = -(gg - k * e) * q / m + A1[5, 3] = -(e + f) / q2 + A1[6, 3] = (gg + k * f) * q / m + iqs * k * (e + f) / m + A1[4, 4] = 1 + A1[5, 4] = -h / q2 + A1[6, 4] = h * k * q / m + h * k * iqs / m + A1[5, 5] = -1 + A1[6, 6] = -1 + + A2[4, 1] = -2 + A2[5, 1] = h / q2 + A2[6, 1] = -(h * k * q / m + h * k * iqs / m) + + J0 = Tinv * SMatrix{6,6,ComplexF64}(A0) * T + J1 = Tinv * SMatrix{6,6,ComplexF64}(A1) * T + J2 = Tinv * SMatrix{6,6,ComplexF64}(A2) * T + return T, (J0, J1, J2) +end + +# The GW2020 Eq. (22) splitting equations for a general exponential block J22 = J₀[3:6, 3:6], kept as its LU. +struct SylvesterSplit{F} + J22lu::F +end + +function wasow_basis(p::CCCParameters, Q::ComplexF64) + T, J = _build_tjmat(p, Q) + return T, J, SylvesterSplit(lu(Matrix(J[1][3:6, 3:6]))) +end + +# With J011 = [0 1; 0 0]: J22·p21 − p21·J011 + K21 = 0 and J011·p12 − p12·J22 + K12 = 0, solved column by column +# and row by row; B is the block-diagonal part of K (GW2020 Eq. 21). +function split_block(s::SylvesterSplit, K::SMatrix{6,6,ComplexF64}) + Bm = zeros(ComplexF64, 6, 6) + Pm = zeros(ComplexF64, 6, 6) + Bm[1:2, 1:2] .= K[1:2, 1:2] + Bm[3:6, 3:6] .= K[3:6, 3:6] + p21_1 = -(s.J22lu \ Vector(K[3:6, 1])) + p21_2 = s.J22lu \ (p21_1 - Vector(K[3:6, 2])) + Pm[3:6, 1] .= p21_1 + Pm[3:6, 2] .= p21_2 + p12_2 = transpose(s.J22lu) \ Vector(K[2, 3:6]) + p12_1 = transpose(s.J22lu) \ (Vector(K[1, 3:6]) + p12_2) + Pm[2, 3:6] .= p12_2 + Pm[1, 3:6] .= p12_1 + return SMatrix{6,6,ComplexF64}(Bm), SMatrix{6,6,ComplexF64}(Pm) +end + +""" + parabolic_scales(params::CCCParameters, Q) -> SVector{2,ComplexF64} + +Scales σⱼ of the two exponential families exp(±x²/(2σⱼ)) at large x: the roots σ² = z of +z² − (1 + m + KQQ*²)Q z + mQ² = 0, Kim, Cho & Choi (1988) Eq. (2) with m = M⊥/M; (√Q, √m·√Q) at Q* = 0. +""" +function parabolic_scales(p::CCCParameters, Q::ComplexF64) + m = p.mperp + b = 1 + m + p.ggj.K * Q * p.Qstar^2 + d = sqrt(b * b - 4 * m) + return SVector(sqrt(Q * (b + d) / 2), sqrt(Q * (b - d) / 2)) +end + +# The ray centred between the two families, θ = (arg σ₁ + arg σ₂)/4, so both decay along it. +default_ray_angle(p::CCCParameters, Q::ComplexF64) = ((σ1, σ2) = parabolic_scales(p, Q); +(angle(σ1) + angle(σ2)) / 4) + +# Slowest decay rate min_j Re(e^{2iθ}/σⱼ) and fastest local rate max_j |e^{2iθ}/σⱼ| of the two families. +parabolic_rate(p::CCCParameters, Q::ComplexF64, θ::Real) = minimum(real(cis(2θ) / σ) for σ in parabolic_scales(p, Q)) +parabolic_resolution_rate(p::CCCParameters, Q::ComplexF64, θ::Real) = maximum(abs(cis(2θ) / σ) for σ in parabolic_scales(p, Q)) + +# Both families decay at every angle in [0, θ]; the condition is linear in the angle, so the end points decide it. +ray_angle_admissible(p::CCCParameters, Q::ComplexF64, θ::Real) = parabolic_rate(p, Q, 0.0) > 0 && parabolic_rate(p, Q, θ) > 0 + +# Fast Ohm pair |√(Q + iQ*)| and Υ-family pair √(|Q(G+KF) + iKQ*(E+F)|/m). +ohm_rate(p::CCCParameters, Q::ComplexF64) = abs(sqrt(Q + ComplexF64(0, p.Qstar))) +function fast_rates(p::CCCParameters, Q::ComplexF64) + g = p.ggj + return ohm_rate(p, Q) + sqrt(abs(Q * (g.G + g.K * g.F) + ComplexF64(0, p.Qstar) * g.K * (g.E + g.F)) / p.mperp) +end diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 48c68e319..9ea7f1520 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -2,8 +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. GGJ is built on the ToroidalLayer library of -# shared Wasow-basis, ray and Galerkin solvers. +# (Glasser–Greene–Johnson), CCC (Cho–Choi–Choi drift-resistive) and SLAYER submodules. GGJ and CCC are +# built on the ToroidalLayer library of shared Wasow-basis, ray and Galerkin solvers. module InnerLayer @@ -16,11 +16,13 @@ using ..Utilities include("InnerLayerInterface.jl") include("ToroidalLayer/ToroidalLayer.jl") include("GGJ/GGJ.jl") +include("CCC/CCC.jl") include("SLAYER/SLAYER.jl") 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 .CCC: CCCModel, CCCParameters import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths @@ -34,12 +36,14 @@ include("LayerPlasma.jl") The `inner_model` token the tearing output records for `model`, e.g. `"ggj_ray"` or `"slayer_fitzpatrick"`. """ model_name(::GGJModel{S}) where {S} = "ggj_$S" +model_name(::CCCModel{S}) where {S} = "ccc_$S" 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 CCC, CCCModel, CCCParameters 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/Ray.jl b/src/InnerLayer/ToroidalLayer/Ray.jl index 8fce660c5..6585cfb3e 100644 --- a/src/InnerLayer/ToroidalLayer/Ray.jl +++ b/src/InnerLayer/ToroidalLayer/Ray.jl @@ -33,7 +33,7 @@ end 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 +its [`parabolic_resolution_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). """ @@ -42,10 +42,11 @@ function ray_breaks(params::ToroidalLayerParameters, Q::ComplexF64, θ::Float64, ρ = 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) + ρres = parabolic_resolution_rate(params, Q, θ) hmax = S / min_cells breaks = [0.0] while breaks[end] < S - h = min(cfl / (μ + ρ * breaks[end]), hmax) + h = min(cfl / (μ + ρres * breaks[end]), hmax) snew = breaks[end] + h push!(breaks, snew >= S - 0.25 * h ? S : snew) end @@ -94,7 +95,7 @@ function decaying_pair(params::ToroidalLayerParameters, Q::ComplexF64, θ::Float for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * (S + ΔS)))) gmax_loc = max(gmax_loc, abs(real(λ))) end - rate = max(ρ * (S + ΔS) + ohm_rate(params, Q), gmax_loc) + rate = max(parabolic_resolution_rate(params, Q, θ) * (S + ΔS) + ohm_rate(params, Q), gmax_loc) nsteps = clamp(ceil(Int, dp_res * rate * ΔS), 200, 2_000_000) h = -ΔS / nsteps @@ -322,6 +323,8 @@ function solve_ray(params::ToroidalLayerParameters, Q::ComplexF64; ztrk::Float64=0.3, purge_budget::Float64=2.5, dp_res::Float64=20.0, handoff::Union{Nothing,NTuple{2,Vector{ComplexF64}}}=nothing, verbose::Bool=false) + ray_angle_admissible(params, Q, θ) || + throw(ArgumentError("contour angle θ=$θ lets an exponential family stop decaying between the real axis and the ray at Q=$Q; use default_ray_angle")) cache = build_asymptotics(params, Q; kmax=kmax) series_ok = true @@ -427,13 +430,20 @@ end Error bar for the matching data: re-solve with each numerical knob perturbed in turn (contour angle, mesh, series order and radius, march tolerance, BVP edge s_m, decaying-pair separation) and -report the relative change of (Δ₊, Δ₋); `spread` is the worst case per parity. θ moves outward -only, away from the pseudo-resonance x² = −Q²(G+KF). +report the relative change of (Δ₊, Δ₋); `spread` is the worst case per parity. θ moves outward, away from the +pseudo-resonance x² = −Q²(G+KF), or inward when the outward angle is not [`ray_angle_admissible`](@ref). """ +# The contour-angle knob of delta_convergence: outward by 20%, or inward by 20% when outward is inadmissible. +function θ_variation(params::ToroidalLayerParameters, Q::ComplexF64, θ::Float64) + θout = θ == 0 ? 0.1 : 1.2 * θ + ray_angle_admissible(params, Q, θout) && return ("θ → 1.2θ", (; θ=θout)) + return ("θ → 0.8θ", (; θ=0.8 * θ)) +end + 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.θ))), + θ_variation(params, Q, base.θ), ("mesh ×2", (; cfl=(base.cfl / 2), min_cells=2 * base.min_cells)), ("kmax+4, S↑", (; kmax=16, smax_tol=1e-10)), ("march_rtol/100", (; march_rtol=1e-11)), diff --git a/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl b/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl index 60ffed836..3f291c5a0 100644 --- a/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl +++ b/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl @@ -39,7 +39,8 @@ subtype only through this interface, so each model keeps whatever numbers its eq - [`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; + far-field rates that set the ray contour and mesh, with optional [`parabolic_resolution_rate`](@ref) and + [`ray_angle_admissible`](@ref) for models whose exponential pairs have different rates; - [`inner_Q`](@ref), [`rescale_delta`](@ref) and [`profile_conversions`](@ref) - the normalization between the physical growth rate, the outer region and the layer variables. """ @@ -124,6 +125,22 @@ Decay rate per unit s² of the slowest parabolic exponential pair along the ray """ function parabolic_rate end +""" + parabolic_resolution_rate(params, Q::ComplexF64, θ) -> Float64 + +Rate per unit s² of the fastest parabolic exponential pair along the ray, which every Hermite cell resolves; +[`parabolic_rate`](@ref) by default, when the pairs share one rate. +""" +parabolic_resolution_rate(params::ToroidalLayerParameters, Q::ComplexF64, θ::Real) = parabolic_rate(params, Q, θ) + +""" + ray_angle_admissible(params, Q::ComplexF64, θ) -> Bool + +Whether every exponential family decays at every angle between the real axis and θ, so the decaying pair on the +ray continues the real-axis one; `true` by default. +""" +ray_angle_admissible(::ToroidalLayerParameters, ::ComplexF64, ::Real) = true + """ ohm_rate(params, Q::ComplexF64) -> Float64 @@ -145,6 +162,7 @@ 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 parabolic_resolution_rate, ray_angle_admissible 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 diff --git a/test/runtests_innerlayer.jl b/test/runtests_innerlayer.jl index 810e115c0..a058c97c4 100644 --- a/test/runtests_innerlayer.jl +++ b/test/runtests_innerlayer.jl @@ -183,3 +183,77 @@ end @test full.Ψ[2, 2] ≈ c.Ψ[1, 2] rtol = 1e-2 # odd Ψ is even, so smooth through x = 0 @test full.certΔ < 1e-5 end + +@testset "InnerLayer CCC (Cho, Choi & Choi drift-resistive layer)" begin + q4 = IL.GGJ.q4_surface_benchmark() + ccc(; Qstar=0.0, mperp=1.0) = IL.CCCParameters(q4, Qstar, mperp) + TL = IL.ToroidalLayer + reldiff(a, b) = maximum(abs.(a .- b) ./ abs.(b)) + p = ccc(; Qstar=0.5, mperp=1.4) + + @testset "Wasow basis matches the layer equations" begin + Q = 0.07 + 0.11im + T, J = IL.CCC._build_tjmat(p, Q) + # J₀ is block-diagonal: power-like (1:2) and exponential (3:6) subspaces split. + @test norm(J[1][1:2, 3:6]) + norm(J[1][3:6, 1:2]) < 1e-12 + # x(A₀ + A₁/x² + A₂/x⁴), A = TJT⁻¹, is ode_matrix in the state u = (xΨ, Ξ, Υ, (xΨ)'/x, Ξ'/x, Υ'/x). + x = 1.7 + 0.4im + Sx = Matrix{ComplexF64}(Diagonal([x, 1, 1, 1, 1 / x, 1 / x])) + Sx[4, 1] = 1 / x + dSx = zeros(ComplexF64, 6, 6) + dSx[1, 1] = 1 + dSx[4, 1] = -1 / x^2 + dSx[5, 5] = dSx[6, 6] = -1 / x^2 + A = [T * J[i] / T for i in 1:3] + @test (Sx * TL.ode_matrix(p, Q, x) + dSx) / Sx ≈ x * (A[1] + A[2] / x^2 + A[3] / x^4) rtol = 1e-12 + # The series solves the equations: its residual (GW2020 Eq. 54) reaches round-off at moderate s. + c = TL.build_asymptotics(p, Q; kmax=12) + @test maximum(TL.asymptotic_residual(c, 10.0 * cis(TL.default_ray_angle(p, Q)))) < 1e-12 + end + + @testset "first-order Ξ coefficient matches the closed form" begin + # Q* enters the power-like solutions first through Ξ = x^{r−2}(1 + b₁/x² + …): + # b₁ − b₁^GGJ = −(i/8) K Q Q* (2s − 1)(1 + 4E + 4F + 4s + 4s²), s = ±√(−D_I) for the large/small column. + Q = 0.07 + 0.11im + cK = TL.build_asymptotics(ccc(; Qstar=0.3), Q; kmax=12) + cG = TL.build_asymptotics(q4, Q; kmax=12) + for (col, s) in ((1, GGJ.p1(q4)), (2, -GGJ.p1(q4))) + ratio(c, x) = (ua = TL.physical_ua_dua(c, x)[1]; x * ua[2, col] / ua[1, col]) + est(x) = x^2 * (ratio(cK, x) - ratio(cG, x)) + Δb1 = (4 * est(200.0) - est(100.0)) / 3 + @test Δb1 ≈ -(im / 8) * q4.K * Q * 0.3 * (2s - 1) * (1 + 4q4.E + 4q4.F + 4s + 4s^2) rtol = 1e-5 + end + end + + @testset "Δ is independent of the contour angle" begin + for Q in (0.5 + 0.5im, 5.0im) + r = TL.solve_ray(p, Q) + θ = TL.default_ray_angle(p, Q) + @test TL.ray_angle_admissible(p, Q, θ) + @test reldiff(TL.solve_ray(p, Q; θ=0.5 * θ).Δ, r.Δ) < 1e-5 + end + # Pin at Q = 5i (delta_convergence spread ~1e-7 there). + Δ = IL.solve_inner(IL.CCCModel(), p, 5.0im * GGJ.q0(q4)) + @test Δ.interchange ≈ 433.38897 + 2358.8069im rtol = 1e-5 + @test Δ.tearing ≈ 22.378564 + 127.277im rtol = 1e-5 + end + + @testset "reduces to GGJ" begin + for Q in (0.1 + 0.05im, 5.0im) + g = TL.solve_ray(q4, Q) + @test reldiff(TL.solve_ray(ccc(), Q).Δ, g.Δ) < 1e-8 # same equations, separate basis and split + @test reldiff(TL.solve_ray(ccc(; Qstar=1e-7, mperp=1 + 1e-7), Q).Δ, g.Δ) < 1e-5 + end + end + + @testset "model and guards" begin + @test IL.CCCModel() === IL.CCCModel{:ray}() + @test IL.model_name(IL.CCCModel()) == "ccc_ray" + @test_throws ArgumentError IL.CCCModel(; solver=:galerkin) + @test_throws ArgumentError IL.CCCModel(; nx=10) + @test_throws ArgumentError IL.CCCParameters(q4, 0.1, -1.0) + # An angle past the decaying sector of the slower family is refused. + @test_throws ArgumentError TL.solve_ray(p, 5.0im; θ=1.2) + @test IL.CCCParameters(; E=q4.E, F=q4.F, G=q4.G, H=q4.H, K=q4.K, taua=q4.taua, taur=q4.taur, Qstar=0.5).ggj.E == q4.E + end +end