InnerLayer.CCC - FEATURE - Add the Cho, Choi and Choi drift-resistive linear MHD inner layer

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)

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

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

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)

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