From 5931e7228e83c518f9f9a3464aec528feb31959d Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 5 Nov 2026 16:59:05 +0100 Subject: [PATCH 1/3] 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 014a258393..5bb8432ea4 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 bb98d6d803..d11d693e55 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 0000000000..896bf6a4bd --- /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 0000000000..f035ea5f73 --- /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 0000000000..9ecad864cb --- /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 48c68e3197..9ea7f15207 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 8fce660c5d..6585cfb3ee 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 60ffed8361..3f291c5a04 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 810e115c01..a058c97c4f 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 From adb3c51fd60ddd15395502c88acaa7571f1d1dae Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 12 Nov 2026 21:09:07 +0100 Subject: [PATCH 2/3] =?UTF-8?q?InnerLayer.CCC=20-=20FEATURE=20-=20Carry=20?= =?UTF-8?q?Q*=20and=20M=E2=8A=A5/M=20on=20CCCModel=20and=20reach=20Tearing?= =?UTF-8?q?,=20Match=20and=20the=20deck=20through=20the=20ToroidalLayer=20?= =?UTF-8?q?interface?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- docs/src/citations.md | 2 +- docs/src/inner_layer.md | 12 +-- docs/src/workflow.md | 4 +- examples/DIIID-like_SLAYER_example/gpec.toml | 2 +- regression-harness/cases/ccc_ray_q5i.toml | 56 +++++++++++++ regression-harness/src/runner.jl | 25 ++++++ src/Deck.jl | 52 +++++++----- src/ForceFreeStates/Matching/MatchProblem.jl | 20 ++--- src/ForceFreeStates/Matching/ResonantMatch.jl | 32 +++---- src/GeneralizedPerturbedEquilibrium.jl | 4 +- src/InnerLayer/CCC/CCC.jl | 83 ++++++++++++------- src/InnerLayer/CCC/CCCParameters.jl | 8 ++ src/InnerLayer/GGJ/GGJ.jl | 3 + src/InnerLayer/GGJ/GGJParameters.jl | 10 +++ src/InnerLayer/InnerLayer.jl | 11 ++- src/InnerLayer/LayerPlasma.jl | 13 +++ src/InnerLayer/ToroidalLayer/Ray.jl | 14 ++-- src/InnerLayer/ToroidalLayer/ToroidalLayer.jl | 44 +++++++++- src/Tearing/Dispersion/Dispersion.jl | 2 +- src/Tearing/Dispersion/GGJGrowthRates.jl | 8 +- src/Tearing/Dispersion/SurfaceCoupling.jl | 8 +- src/Tearing/GGJTearing.jl | 19 +++-- src/Tearing/HDF5Output.jl | 40 ++++----- src/Tearing/SLAYERTearing.jl | 2 +- src/Tearing/Tearing.jl | 4 +- src/Tearing/TearingProblem.jl | 4 +- test/runtests_innerlayer.jl | 8 ++ test/runtests_slayer_runner.jl | 21 ++++- test/runtests_solve_api.jl | 56 +++++++++++++ 29 files changed, 415 insertions(+), 152 deletions(-) create mode 100644 regression-harness/cases/ccc_ray_q5i.toml diff --git a/docs/src/citations.md b/docs/src/citations.md index 5bb8432ea4..cc7aefb1c8 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -155,7 +155,7 @@ Derives the toroidal resistive-layer equations with electron diamagnetic drift, > *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 Δ. +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) plot the Fourier-space Δ_o, which its Eq. (10b) relates to the real-space Δ that `CCCModel` returns. --- diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index d11d693e55..34dfbea164 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -5,15 +5,16 @@ The outer region, ideal MHD everywhere except at the rational surfaces, gives th ## What it computes -Two layer models are available, selected by `[Tearing] inner_model`: +Three layer models are available, selected by `[Tearing] inner_model`: - `"SLAYER"`, the two-fluid drift-MHD slab layer equations of Cole and Fitzpatrick (2006), in the form Park's SLAYER solver takes them with a few terms dropped ([Park (2022)](citations.md#Tearing), [Burgess et al. (2026)](citations.md#Tearing)), with resistivity, perpendicular and toroidal transport, diamagnetic drifts and E×B rotation taken from the kinetic profiles of the equilibrium; its options are in `[SLAYER]`. - `"GGJ"`, the toroidal resistive-MHD layer of Glasser, Greene and Johnson ([Glasser, Wang and Park (2016)](citations.md#Tearing), [Glasser and Wang (2020)](citations.md#Tearing)), whose coefficients E, F, G, H, K are the flux-surface averages stored under `SingularSurfaces/`; its solver options are in `[GGJ]`. +- `"CCC"`, the GGJ layer with electron diamagnetic drift and perpendicular transport of [Cho, Choi and Choi (1987)](citations.md#Tearing) (see [CCC](@ref "CCC: diamagnetic drift and perpendicular transport") below); its solver options, Q* and M⊥/M are in `[CCC]`. For each surface the `Dispersion` search scans the complex frequency plane, extracts the roots of the dispersion relation and classifies them, and `[Tearing.search]` sets the scan. With `[Tearing.CriticalResonantField]` the stage also computes, for the `SLAYER` layer, the resonant field at which a static error field overcomes the layer's shielding: the error-field penetration threshold. -The same GGJ layer serves the `[Match]` stage of the [Workflow](workflow.md), which does the complementary job: instead of finding the free modes it imposes the layer response on the driven force-free states and writes the matched solution to `SingularSurfaces/Match/`. +The same GGJ layer, or the CCC layer, serves the `[Match]` stage of the [Workflow](workflow.md), which does the complementary job: instead of finding the free modes it imposes the layer response on the driven force-free states and writes the matched solution to `SingularSurfaces/Match/`. ## Governing equations @@ -54,7 +55,7 @@ From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` ```toml [Tearing] -inner_model = "SLAYER" # or "GGJ" +inner_model = "SLAYER" # or "GGJ", "CCC" coupling_mode = "uncoupled" # one dispersion relation per surface; "coupled" uses the full Δ′ matrix [Tearing.search] @@ -95,8 +96,8 @@ The DIII-D-like kinetic file of the examples carries analytic χ profiles of the | Dataset | Meaning | |---|---| | `Tearing/PerSurface/rational_psi`, `rational_q`, `m`, `n` | the analyzed surfaces | -| `Tearing/PerSurface/S`, `tau`, `c_beta`, `P_perp`, `P_tor`, `Q_e`, `Q_i`, `tau_k`, `tau_R`, `eta`, `rs`, … | the layer parameters of each surface (SLAYER); `E`, `F`, `G`, `H`, `K`, `M`, `tau_A` and `tau_R` for GGJ | -| `Tearing/PerSurface/Delta_prime_raw` and `Delta_prime_matrix_rs` (SLAYER) or `Delta_prime_matrix` (GGJ) | the outer matrix each layer was matched to | +| `Tearing/PerSurface/S`, `tau`, `c_beta`, `P_perp`, `P_tor`, `Q_e`, `Q_i`, `tau_k`, `tau_R`, `eta`, `rs`, … | the layer parameters of each surface (SLAYER); `E`, `F`, `G`, `H`, `K`, `M`, `tau_A` and `tau_R` for GGJ, with `Q_star` and `M_perp_ratio` for CCC | +| `Tearing/PerSurface/Delta_prime_raw` and `Delta_prime_matrix_rs` (SLAYER) or `Delta_prime_matrix` (GGJ, CCC) | the outer matrix each layer was matched to | | `Tearing/PerSurface/Delta_prime_norm`, `delta_prime_conversion` | the Δ′ normalization of the SLAYER layer | | `Tearing/Roots/Q_root`, `omega`, `gamma`, `no_root` | the dispersion root of each surface in the lab frame (`frequency_frame` attribute), its real frequency and growth rate, and the flag of a surface with no usable root | | `Tearing/LayerWidths/delta_s`, `delta_s_abs`, `d_beta` | the complex resistive layer thickness, its magnitude in metres, and the ion drift scale | @@ -144,6 +145,7 @@ p = GGJ.q4_surface_benchmark() # GGJParameters for the q = 4 benchmark surf [`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. +In a deck, `[Tearing]` or `[Match]` takes `inner_model = "CCC"`, and `[CCC]` holds the solver options with `Qstar` and `mperp`, each a scalar or one value per surface, as the `CCCModel` keywords of the same names; without `Qstar`, Q* = −n·ω*ₑ/Q₀ comes from the kinetic data ([`layer_omega_star_e`](@ref InnerLayer.layer_omega_star_e)), and `mperp` defaults to 1. 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 diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 38b4cdd8ba..a76507dd72 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -68,7 +68,7 @@ Page: [Local Stability](ballooning.md). ### Match -Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, the only model, with its options in an optional `[GGJ]` section) and rebuilds the force-free states with the resulting layer response, writing the matched solution and the penetrated resonant fields to `SingularSurfaces/Match/`; `ideal = true` keeps the shielded limit for comparison. +Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, the default, or the drift-resistive `"CCC"`, with its options in an optional section of the same name) and rebuilds the force-free states with the resulting layer response, writing the matched solution and the penetrated resonant fields to `SingularSurfaces/Match/`; `ideal = true` keeps the shielded limit for comparison. Outputs: `SingularSurfaces/Match/`. Page: [Tearing](inner_layer.md). @@ -98,7 +98,7 @@ Page: [Error Fields](error_fields.md). ### Tearing -Solves the tearing dispersion relation of each rational surface by matching the outer-region Δ′ to an inner-layer model: the two-fluid slab `SLAYER` (`[SLAYER]`) or the toroidal resistive-MHD `GGJ` (`[GGJ]`), selected by `[Tearing] inner_model`. +Solves the tearing dispersion relation of each rational surface by matching the outer-region Δ′ to an inner-layer model: the two-fluid slab `SLAYER` (`[SLAYER]`), the toroidal resistive-MHD `GGJ` (`[GGJ]`) or its drift-resistive extension `CCC` (`[CCC]`), selected by `[Tearing] inner_model`. `[Tearing.search]` sets the root search in the complex frequency plane, and `[Tearing.CriticalResonantField]` adds the resonant field at which the layer's shielding is overcome. Needs the Δ′ matrix of `"Riccati"` with `vac_flag = true` or `"Galerkin"` with `[Galerkin] rpec_flag = true`, and the kinetic profiles of `[Equilibrium] kinetic_file` (`SLAYER` reads its densities, temperatures, rotation and diffusivities there). Outputs: `Tearing/PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/`. diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index fa41a0dd50..3d9ad8877b 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -68,7 +68,7 @@ renorm_threshold = 1e3 # Column-norm threshold that triggers solution re # large and likely invalid, no Re/Im contour crossing polishes to a true zero, and the # validity gate correctly reports them as no_root rather than a spurious growth rate. [Tearing] -inner_model = "SLAYER" # Inner-layer model, stated in its own section: "SLAYER" or "GGJ" +inner_model = "SLAYER" # Inner-layer model, stated in its own section: "SLAYER", "GGJ" or "CCC" coupling_mode = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant) # omega_E_kHz = {"2/1" = 0.0, "3/1" = 3.0} # override of the E×B rotation Ω_E/2π per unit n [kHz], keyed by m/n; unlisted surfaces use the kinetic file diff --git a/regression-harness/cases/ccc_ray_q5i.toml b/regression-harness/cases/ccc_ray_q5i.toml new file mode 100644 index 0000000000..85f9a47487 --- /dev/null +++ b/regression-harness/cases/ccc_ray_q5i.toml @@ -0,0 +1,56 @@ +# Regression case: CCC (Cho, Choi & Choi drift-resistive) inner layer on the rotated ray. +# A "computed" case (kind = "computed", no example_dir) that calls solve_inner with CCCModel on the +# physical q=4 benchmark surface with Q* = 0.5 (electron diamagnetic drift) and M⊥/M = 1.4 +# (perpendicular transport) at Q = 5i, where the two exponential families of the far field separate, +# and tracks the interchange (Δ₊) / tearing (Δ₋) matching coefficients. Guards the CCC equations, its +# Wasow basis and splitting solve, and the centred contour angle against numerical drift across +# changes. Each [quantities.*] block names an HDF5 path in the computed output, how to extract it, +# and the noise floor for diffs. +[case] +name = "ccc_ray_q5i" +description = "CCC inner layer on the ray, q=4 surface, Q* = 0.5, M⊥/M = 1.4, Q = 5i" +kind = "computed" + +# interchange (Δ₊) / tearing (Δ₋) from solve_inner(CCCModel(), CCCParameters(q4_surface_benchmark(), 0.5, 1.4), 5i·Q₀). +# Complex-valued; both real and imaginary parts are tracked separately. The components are O(20–2400); +# noise_threshold 1e-3 (absolute) sits above the ~1e-7 relative delta_convergence spread at this point +# and the cross-machine roundoff of the extended-precision march. +[quantities.delta_interchange_real] +h5path = "ccc/delta_interchange_real" +type = "real_scalar" +extract = "value" +label = "Re(Δ₊ interchange)" +noise_threshold = 1e-3 +order = 10 + +[quantities.delta_interchange_imag] +h5path = "ccc/delta_interchange_imag" +type = "real_scalar" +extract = "value" +label = "Im(Δ₊ interchange)" +noise_threshold = 1e-3 +order = 11 + +[quantities.delta_tearing_real] +h5path = "ccc/delta_tearing_real" +type = "real_scalar" +extract = "value" +label = "Re(Δ₋ tearing)" +noise_threshold = 1e-3 +order = 12 + +[quantities.delta_tearing_imag] +h5path = "ccc/delta_tearing_imag" +type = "real_scalar" +extract = "value" +label = "Im(Δ₋ tearing)" +noise_threshold = 1e-3 +order = 13 + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 90 diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index cb1f5959e7..4cc465e359 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -157,6 +157,29 @@ end %RUNINFO% """ +# CCC (Cho, Choi & Choi drift-resistive) layer on the rotated ray: the q=4 benchmark surface with Q* = 0.5 and +# M⊥/M = 1.4 at Q = 5i, where its two exponential families separate. Writes the parity matching data; runtime and +# environment go to ARGS[2] via the run-info epilogue. +const COMPUTED_CCC_RAY_SCRIPT_TEMPLATE = """ +using Pkg +%INSTANTIATE% +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.InnerLayer +using HDF5 +p = CCCParameters(InnerLayer.GGJ.q4_surface_benchmark(), 0.5, 1.4) +γ = 5.0im * InnerLayer.GGJ.q0(p.ggj) +t_start = time() +Δ = solve_inner(CCCModel(), p, γ) +elapsed = time() - t_start +h5open(ARGS[1], "w") do fid + fid["ccc/delta_interchange_real"] = real(Δ.interchange) + fid["ccc/delta_interchange_imag"] = imag(Δ.interchange) + fid["ccc/delta_tearing_real"] = real(Δ.tearing) + fid["ccc/delta_tearing_imag"] = imag(Δ.tearing) +end +%RUNINFO% +""" + # Self-contained separatrix-finder regression (PR #296). Loads a fixed-boundary EFIT whose # computational box hugs the LCFS (eps=0.05 TokaMaker aspect-scan g-file): outside the prescribed # LCFS the coil-vacuum flux turns back above the boundary value before the grid edge, so the old @@ -346,6 +369,8 @@ function _computed_script_template(case_spec::CaseSpec) return COMPUTED_GGJ_SCRIPT_TEMPLATE elseif case_spec.name == "ggj_ray_q500i" return COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE + elseif case_spec.name == "ccc_ray_q5i" + return COMPUTED_CCC_RAY_SCRIPT_TEMPLATE elseif case_spec.name == "efit_fixedbdy_separatrix" return COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE end diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17ce..03821a95ad 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -76,9 +76,13 @@ const _RETIRED_SECTIONS = Dict("DEBUG" => "[DEBUG] is retired: output_benchmark_ "SOL_INPUT" => "[SOL_INPUT] is now [Solovev]" * _ANALYTIC_HOME, "LAR_INPUT" => "[LAR_INPUT] is now [LargeAspectRatio]" * _ANALYTIC_HOME, "TJ_ANALYTIC_INPUT" => "[TJ_ANALYTIC_INPUT] is now [TJAnalytic], whose solver = \"direct\" replaces eq_type = \"tj_analytic_direct\"" * _ANALYTIC_HOME) +# The toroidal layer models a deck can name as `inner_model`, each configured by the section of its token. +const TOROIDAL_LAYER_MODELS = ("GGJ" => GGJModel, "CCC" => CCCModel) +const _TOROIDAL_LAYER_TOKENS = first.(TOROIDAL_LAYER_MODELS) + # Every section a deck may hold, besides the analytic-equilibrium sections. const _DECK_SECTIONS = ("Run", "Equilibrium", "Wall", "LocalStability", "Ballooning", "ForceFreeStates", "Forward", "Riccati", - "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields") + "Galerkin", "Match", _TOROIDAL_LAYER_TOKENS..., "Tearing", "SLAYER", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields") const _INTEGRATORS = ("Forward" => Forward, "Riccati" => Riccati, "Galerkin" => Galerkin) # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated. @@ -250,17 +254,21 @@ function slayer_model(inputs::AbstractDict) end """ - ggj_model(inputs) -> GGJModel + toroidal_layer_model(inputs, token) -> InnerLayer.ToroidalLayerModel -The GGJ layer model of the `[GGJ]` section: `solver` and the backend options, passed to -`GGJModel` as they are. The default `GGJModel()` when absent. +The toroidal layer model `token` names in `TOROIDAL_LAYER_MODELS` (`"GGJ"`, `"CCC"`), from its `[token]` section: +`solver` and the model's keywords, passed to its constructor as they are. The model's defaults when absent. """ -function ggj_model(inputs::AbstractDict) - kwargs = Dict{Symbol,Any}(Symbol(k) => v for (k, v) in get(inputs, "GGJ", Dict{String,Any}())) +function toroidal_layer_model(inputs::AbstractDict, token::AbstractString) + kwargs = Dict{Symbol,Any}(Symbol(k) => v for (k, v) in get(inputs, token, Dict{String,Any}())) haskey(kwargs, :solver) && (kwargs[:solver] = Symbol(kwargs[:solver])) - return GGJModel(; kwargs...) + return Dict(TOROIDAL_LAYER_MODELS)[token](; kwargs...) end +# The `inner_model` token of a layer model, as a deck names it. +_inner_model_token(::SLAYERModel) = "SLAYER" +_inner_model_token(model::InnerLayer.ToroidalLayerModel) = only(token for (token, T) in TOROIDAL_LAYER_MODELS if model isa T) + """ force_balance_keywords(inputs, path) -> Dict{Symbol,Any} @@ -340,15 +348,15 @@ end match_from_deck(inputs) -> Union{Nothing,NamedTuple} `(; problem, model)`: the `MatchProblem` keywords of `[Match]` and the inner-layer model its -`inner_model` names (`"GGJ"`, from `[GGJ]`). `nothing` without `[Match]`. +`inner_model` names (`"GGJ"` by default, or `"CCC"`, from the section of that name). `nothing` without `[Match]`. """ function match_from_deck(inputs::AbstractDict) haskey(inputs, "Match") || return nothing fields = Dict{Symbol,Any}(:eta => Any, :rho => Any, :rotation => Any, :gamma => Any, :ideal => Bool, :inner_model => String) problem = section_keywords(inputs["Match"], "Match", fields) - _selector(inputs["Match"], "Match", "inner_model", ("GGJ",), "GGJ") + inner_model = _selector(inputs["Match"], "Match", "inner_model", _TOROIDAL_LAYER_TOKENS, "GGJ") delete!(problem, :inner_model) - return (; problem, model=ggj_model(inputs)) + return (; problem, model=toroidal_layer_model(inputs, inner_model)) end """ @@ -356,7 +364,7 @@ end `(; problem, model, config)`: the `TearingProblem` keywords of `[Tearing]` (with `[Tearing.CriticalResonantField]`), the layer model its `inner_model` names (`"SLAYER"` from -`[SLAYER]`, `"GGJ"` from `[GGJ]`) and the `TearingConfig` of `[Tearing.search]`. `nothing` without `[Tearing]`. +`[SLAYER]`, `"GGJ"` or `"CCC"` from the section of that name) and the `TearingConfig` of `[Tearing.search]`. `nothing` without `[Tearing]`. """ function tearing_from_deck(inputs::AbstractDict) haskey(inputs, "Tearing") || return nothing @@ -366,7 +374,7 @@ function tearing_from_deck(inputs::AbstractDict) fields = _struct_fields(TearingProblem, (:ffs, :critical_resonant_field)) fields[:inner_model] = String problem = section_keywords(_without(table, "search", "CriticalResonantField"), "Tearing", fields) - inner_model = _selector(table, "Tearing", "inner_model", ("SLAYER", "GGJ")) + inner_model = _selector(table, "Tearing", "inner_model", ("SLAYER", _TOROIDAL_LAYER_TOKENS...)) delete!(problem, :inner_model) crf = get(table, "CriticalResonantField", nothing) crf === nothing || ( @@ -374,7 +382,7 @@ function tearing_from_deck(inputs::AbstractDict) section_keywords(crf, "Tearing.CriticalResonantField", _struct_fields(CriticalResonantFieldControl))...) ) config = TearingConfig(; section_keywords(get(table, "search", Dict{String,Any}()), "Tearing.search", _struct_fields(TearingConfig))...) - return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config) + return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : toroidal_layer_model(inputs, inner_model), config) end """ @@ -447,7 +455,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) end has_ffs = haskey(inputs, "ForceFreeStates") if !has_ffs - for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "LocalStability", "Ballooning", + for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", _TOROIDAL_LAYER_TOKENS..., "Tearing", "SLAYER", "LocalStability", "Ballooning", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields") haskey(inputs, section) && error("[$section] is set but there is no [ForceFreeStates] section, so the run stops after the equilibrium") end @@ -479,8 +487,11 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) error("[ErrorFields] linearizes coil geometry, so [ForcingTerms] must use forcing_data_format = \"coil\" (got \"$(forcing.forcing_data_format)\")") end tearing_model = haskey(inputs, "Tearing") ? get(inputs["Tearing"], "inner_model", nothing) : nothing - haskey(inputs, "GGJ") && match === nothing && tearing_model != "GGJ" && - error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it") + match_model = match === nothing ? nothing : _inner_model_token(match.model) + for token in _TOROIDAL_LAYER_TOKENS + haskey(inputs, token) && match_model != token && tearing_model != token && + error("[$token] is set but neither [Match] nor [Tearing] inner_model = \"$token\" uses it") + end haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) && error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it") return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, @@ -521,7 +532,8 @@ end deck_sections(ctrl::KineticForces.KineticForcesControl) = Dict{String,Any}("KineticForces" => _deck_table(ctrl)) deck_sections(model::SLAYERModel) = Dict{String,Any}("SLAYER" => merge(_deck_table(model), Dict{String,Any}("variant" => String(typeof(model).parameters[1])))) -deck_sections(model::GGJModel) = Dict{String,Any}("GGJ" => merge(_deck_table(Dict(pairs(model.options))), Dict{String,Any}("solver" => String(typeof(model).parameters[1])))) +deck_sections(model::InnerLayer.ToroidalLayerModel{S}) where {S} = + Dict{String,Any}(_inner_model_token(model) => merge(_deck_table(model, (:options,)), _deck_table(Dict(pairs(model.options))), Dict{String,Any}("solver" => String(S)))) deck_sections(wall::Vacuum.WallSettings) = Dict{String,Any}("Wall" => merge(_deck_table(wall), Dict{String,Any}("shape" => only(k for (k, T) in Vacuum.WALL_SHAPES if wall isa T)))) """ @@ -571,7 +583,7 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an ls === nothing || merge!(inputs, Dict{String,Any}("LocalStability" => _deck_table(ls.problem), "Ballooning" => _deck_table(ls.algorithm))) if deck.match !== nothing merge!(inputs, deck_sections(deck.match.model)) - inputs["Match"] = merge(_deck_table(deck.match.problem), Dict{String,Any}("inner_model" => "GGJ")) + inputs["Match"] = merge(_deck_table(deck.match.problem), Dict{String,Any}("inner_model" => _inner_model_token(deck.match.model))) end if deck.tearing !== nothing (; problem, model, config) = deck.tearing @@ -590,7 +602,7 @@ end # [Tearing] with its search and critical-resonant-field sub-tables, and the section of its layer model. function _tearing_sections(problem, crf, model, config) table = _deck_table(problem, (:ffs, :critical_resonant_field)) - table["inner_model"] = model isa SLAYERModel ? "SLAYER" : "GGJ" + table["inner_model"] = _inner_model_token(model) table["search"] = _deck_table(config) crf === nothing || (table["CriticalResonantField"] = _deck_table(crf)) return merge(deck_sections(model), Dict{String,Any}("Tearing" => table)) @@ -623,7 +635,7 @@ const _RECORDED_TOLERANCES = "Input/RawInputs/ErrorFields/tolerance_toml_raw" deck_sections(equil::Equilibrium.PlasmaEquilibrium) = deck_sections(equil.config, equil.analytic) deck_sections(ls::LocalStabilityResult) = Dict{String,Any}("LocalStability" => _deck_table(ls.problem, (:equil,)), "Ballooning" => _deck_table(ls.algorithm)) deck_sections(m::ForceFreeStates.MatchResult) = - merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ")))) + merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => _inner_model_token(m.model))))) deck_sections(kf::KineticForcesResult) = deck_sections(kf.control) deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config) deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES) diff --git a/src/ForceFreeStates/Matching/MatchProblem.jl b/src/ForceFreeStates/Matching/MatchProblem.jl index f80f52eb46..d96e0904ae 100644 --- a/src/ForceFreeStates/Matching/MatchProblem.jl +++ b/src/ForceFreeStates/Matching/MatchProblem.jl @@ -7,7 +7,7 @@ Driven inner-layer matching on a finished force-free-states solve `ffs`: each rational surface is matched at the forced eigenvalue γ = 2πi·n·f, with f the plasma rotation frequency there. -`solve(prob, GGJModel())` returns a copy of `ffs` with the matched solution and the penetrated +`solve(prob, GGJModel())` (or another toroidal layer model, e.g. `CCCModel()`) returns a copy of `ffs` with the matched solution and the penetrated field `bpen`. `ffs` needs the Δ′ coil block, from `Galerkin(; rpec_flag=true)` or `Riccati()` with `vac_flag=true`. @@ -84,11 +84,11 @@ end """ closure_capable(model) -> Bool -Whether `model` can solve a [`MatchProblem`](@ref): `GGJModel` with the `:ray` or `:galerkin` -backend can (they reconstruct the layer profiles); `SLAYERModel` (slab, tearing parity only) +Whether `model` can solve a [`MatchProblem`](@ref): a toroidal layer model (`GGJModel`, `CCCModel`) with the +`:ray` or `:galerkin` backend can (they reconstruct the layer profiles); `SLAYERModel` (slab, tearing parity only) cannot. """ -closure_capable(::InnerLayer.GGJModel{S}) where {S} = S in (:ray, :galerkin) +closure_capable(::InnerLayer.ToroidalLayerModel{S}) where {S} = S in (:ray, :galerkin) closure_capable(::InnerLayer.SLAYERModel) = false """ @@ -107,7 +107,7 @@ function CommonSolve.solve(prob::MatchProblem, model::InnerLayer.InnerLayerModel end # Inner-layer Δ per surface, the matching system, then the matched outer and inner profiles. -function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) +function _compute_match(prob::MatchProblem, model::InnerLayer.ToroidalLayerModel) ffs = prob.ffs dp = ffs.delta_prime sings = ffs.surfaces @@ -128,7 +128,7 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) inner_psi = Vector{Float64}[] inner_xi = Matrix{ComplexF64}[] inner_b = Matrix{ComplexF64}[] - inner_params = InnerLayer.GGJParameters[] + inner_params = InnerLayer.ToroidalLayerParameters[] rpec_eig = zeros(ComplexF64, msing) residual = 0.0 else @@ -142,12 +142,12 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₊, even inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₋, even inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₊, odd - inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) + InnerLayer.check_layer_inputs(model, msing) + inner_params = [InnerLayer.layer_parameters(model, ggj_parameters(sings[i], equil; eta=prob.eta[i], rho=prob.rho[i], gamma=prob.gamma, ising=i), + equil, sings[i]) for i in 1:msing] profiles = true for i in 1:msing - params = ggj_parameters(sings[i], equil; eta=prob.eta[i], rho=prob.rho[i], - gamma=prob.gamma, ising=i) - inner_params[i] = params + params = inner_params[i] γ = 2π * im * nn * prob.rotation[i] # forced eigenvalue, f in Hz rpec_eig[i] = γ inner = InnerLayer.solve_inner_profile(model, params, γ) # parity columns (Δ₊, Δ₋) diff --git a/src/ForceFreeStates/Matching/ResonantMatch.jl b/src/ForceFreeStates/Matching/ResonantMatch.jl index 3c542f3327..a53a84e3b7 100644 --- a/src/ForceFreeStates/Matching/ResonantMatch.jl +++ b/src/ForceFreeStates/Matching/ResonantMatch.jl @@ -19,7 +19,7 @@ fields are empty when the solve kept no Galerkin basis; the inner fields are emp - `xi`, `xi_deriv`: `(mpert, ngrid, ncoil)` matched outer ξ and ξ′, one column per coil. - `inner_psi`, `inner_xi`, `inner_b`: per surface, the inner ψ grid and the composite ξ_ψ and b^ψ on it; empty unless the inner model returns real-ψ profiles (`GGJModel(; solver=:galerkin)` or `GGJModel(; real_profiles=true)`). - - `inner_params`: per-surface `GGJParameters` used. + - `inner_params`: per-surface layer parameters used (`GGJParameters`, `CCCParameters`). - `model`: the inner-layer model solved, with its backend options. - `problem`: the [`MatchProblem`](@ref) solved, with its per-surface η, ρ and rotation. """ @@ -36,7 +36,7 @@ struct MatchResult{P} inner_psi::Vector{Vector{Float64}} inner_xi::Vector{Matrix{ComplexF64}} inner_b::Vector{Matrix{ComplexF64}} - inner_params::Vector{InnerLayer.GGJParameters} + inner_params::Vector{<:InnerLayer.ToroidalLayerParameters} model::InnerLayer.InnerLayerModel problem::P end @@ -90,20 +90,17 @@ function save_to_hdf5(out_h5::Union{HDF5.File,HDF5.Group}, m::MatchResult) out_h5["$g/Inner/b_$i"] = m.inner_b[i] end out_h5["$g/residual"] = m.residual - if !isempty(m.inner_params) - for f in (:E, :F, :G, :H, :K, :M) - out_h5["$g/InnerParams/$(f)"] = [getfield(pp, f) for pp in m.inner_params] - end - # Literature names, matching the Tearing PerSurface mapping for the same fields. - out_h5["$g/InnerParams/tau_A"] = [pp.taua for pp in m.inner_params] - out_h5["$g/InnerParams/tau_R"] = [pp.taur for pp in m.inner_params] - out_h5["$g/InnerParams/dVdpsi"] = [pp.v1 for pp in m.inner_params] + # The layer model's per-surface parameters, the same columns as the Tearing PerSurface group. + columns = isempty(m.inner_params) ? [] : InnerLayer.layer_parameter_columns(m.inner_params) + for (name, column) in columns + out_h5["$g/InnerParams/$name"] = column.values end a = attrs(out_h5[g]) a["inner_model"] = InnerLayer.model_name(m.model) - m.model isa InnerLayer.GGJModel && (a["inner_model_options"] = string(m.model.options)) + m.model isa InnerLayer.ToroidalLayerModel && (a["inner_model_options"] = string(m.model.options)) inner = ["$g/Inner/$(v)_$i" => meta for i in eachindex(m.inner_psi) for (v, meta) in _MATCH_INNER_H5_ANNOTATIONS] - Utilities.HDF5Annotations.annotate!(out_h5, vcat(MATCH_H5_ANNOTATIONS, inner)) + params = ["$g/InnerParams/$name" => (; column.long_name, column.units, _MATCH_SURF...) for (name, column) in columns] + Utilities.HDF5Annotations.annotate!(out_h5, vcat(MATCH_H5_ANNOTATIONS, inner, params)) return nothing end @@ -121,16 +118,7 @@ const MATCH_H5_ANNOTATIONS = [ "SingularSurfaces/Match/reconnected_flux" => (; long_name="reconnected flux Δ_coil + Δ_outᵀ·cout per surface parity, per unit coil-response column", dims=("surface_parity", "coil")), "SingularSurfaces/Match/rpec_eig" => (; long_name="forced eigenvalue γ = 2πi·n·f at each surface", units="rad/s", _MATCH_SURF...), - "SingularSurfaces/Match/residual" => (; long_name="relative residual of the matching linear solve"), - "SingularSurfaces/Match/InnerParams/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/H" => (; long_name="Glasser-Greene-Johnson coefficient H per surface", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/K" => (; long_name="Glasser-Greene-Johnson coefficient K per surface", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/tau_A" => (; long_name="Alfvén time τ_A per surface", units="s", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/tau_R" => (; long_name="resistive diffusion time τ_R per surface", units="s", _MATCH_SURF...), - "SingularSurfaces/Match/InnerParams/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", _MATCH_SURF...) + "SingularSurfaces/Match/residual" => (; long_name="relative residual of the matching linear solve") ] const _MATCH_INNER_H5_ANNOTATIONS = [ "psi" => (; long_name="normalized poloidal flux ψ_N grid across the inner layer of this surface", dims=("inner_psi",)), diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 76c5515b57..fe53a6696a 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -73,7 +73,7 @@ include("Rerun.jl") # constructor and the forcing description, re-exported so a user needs one `using`. using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult using .ForceFreeStates: MatchProblem -using .InnerLayer: GGJModel, SLAYERModel +using .InnerLayer: GGJModel, CCCModel, SLAYERModel using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult using .Equilibrium: PlasmaEquilibrium, KineticProfileData, IonSpecies, attach_kinetic_profiles! @@ -304,7 +304,7 @@ export main, write_imas export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields export KineticProfileData, IonSpecies, attach_kinetic_profiles! export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse -export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel +export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, CCCModel, SLAYERModel export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult diff --git a/src/InnerLayer/CCC/CCC.jl b/src/InnerLayer/CCC/CCC.jl index 896bf6a4bd..f8b6adefbd 100644 --- a/src/InnerLayer/CCC/CCC.jl +++ b/src/InnerLayer/CCC/CCC.jl @@ -16,69 +16,96 @@ using StaticArrays import ..InnerLayerResponse import ..solve_inner, ..solve_inner_profile +import ..layer_omega_star_e 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 +import ..ToroidalLayer: layer_parameters, check_layer_inputs, mercier_di, growth_rate_scale, layer_parameter_columns +using ..GGJ: GGJParameters, q0 """ - CCCModel{S,O} <: ToroidalLayerModel{S} - CCCModel(; solver=:ray, options...) + CCCModel{S,O,Q,M} <: ToroidalLayerModel{S} + CCCModel(; solver=:ray, Qstar=nothing, mperp=1.0, 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. +[`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. `Qstar` and `mperp` set +the parameters the tearing and matching stages build ([`layer_parameters`](@ref)); a solve on given +`CCCParameters` uses theirs. ## Fields - `options::NamedTuple` - Backend keywords forwarded to every solve; call-site keywords win. + - `Qstar` - Scaled electron diamagnetic frequency Q* at every surface (a number), at each surface core to edge + (a vector), or `nothing` to take Q* = −n·ω*ₑ/Q₀ from the equilibrium's kinetic data + ([`layer_omega_star_e`](@ref)). + - `mperp` - Perpendicular-transport weight M⊥/M at every surface (a number) or at each surface core to edge (a + vector); positive. """ -struct CCCModel{S,O<:NamedTuple} <: ToroidalLayerModel{S} +struct CCCModel{S,O<:NamedTuple,Q<:Union{Nothing,Float64,Vector{Float64}},M<:Union{Float64,Vector{Float64}}} <: ToroidalLayerModel{S} options::O + Qstar::Q + mperp::M end -function CCCModel(; solver::Symbol=:ray, options...) +function CCCModel(; solver::Symbol=:ray, Qstar=nothing, mperp=1.0, 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...) + per_surface(v) = v isa AbstractVector ? collect(Float64, v) : Float64(v) + qstar = Qstar === nothing ? nothing : per_surface(Qstar) + qstar === nothing || all(isfinite, qstar) || throw(ArgumentError("CCCModel: Qstar must be finite")) + m = per_surface(mperp) + all(v -> isfinite(v) && v > 0, m) || throw(ArgumentError("CCCModel: mperp must be positive and finite")) + return CCCModel{solver,typeof(values(options)),typeof(qstar),typeof(m)}(values(options), qstar, m) end +CCCModel{S}() where {S} = CCCModel(; solver=S) include("CCCParameters.jl") include("Equations.jl") """ - solve_inner(::CCCModel{:ray}, params::CCCParameters, γ::Number; kwargs...) -> InnerLayerResponse + solve_inner(m::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. +([`ray_response`](@ref)); the model's options and `kwargs` 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(m::CCCModel{:ray}, params::CCCParameters, γ::Number; kwargs...) = ray_response(params, γ; m.options..., kwargs...) """ - solve_inner_profile(::CCCModel{:ray}, params::CCCParameters, γ::Number; kwargs...) + solve_inner_profile(m::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...) +solve_inner_profile(m::CCCModel{:ray}, params::CCCParameters, γ::Number; kwargs...) = ray_profile(params, γ; m.options..., kwargs...) + +# A per-surface model input at surface `k`: the scalar for every surface, or the vector's entry. +_at_surface(v::Float64, k::Integer) = v +_at_surface(v::Vector{Float64}, k::Integer) = v[k] + +function check_layer_inputs(model::CCCModel, nsurfaces::Integer) + for (name, v) in (("Qstar", model.Qstar), ("mperp", model.mperp)) + v isa Vector && length(v) != nsurfaces && + throw(ArgumentError("CCCModel: $name has length $(length(v)), expected $nsurfaces (one per surface, core to edge)")) + end + return nothing +end + +function layer_parameters(model::CCCModel, ggj::GGJParameters, equil, surface) + Qstar = if model.Qstar === nothing + equil.kinetic === nothing && + throw(ArgumentError("CCC layer: give Qstar, or load kinetic profiles (kinetic_file in [Equilibrium], or attach_kinetic_profiles!)")) + -only(surface.n) * only(layer_omega_star_e(equil, (surface,))) / q0(ggj) + else + _at_surface(model.Qstar, ggj.ising) + end + return CCCParameters(ggj, Qstar, _at_surface(model.mperp, ggj.ising)) +end export CCCModel, CCCParameters, parabolic_scales diff --git a/src/InnerLayer/CCC/CCCParameters.jl b/src/InnerLayer/CCC/CCCParameters.jl index f035ea5f73..3dd50f22ea 100644 --- a/src/InnerLayer/CCC/CCCParameters.jl +++ b/src/InnerLayer/CCC/CCCParameters.jl @@ -39,3 +39,11 @@ 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) +mercier_di(p::CCCParameters) = mercier_di(p.ggj) +growth_rate_scale(p::CCCParameters) = growth_rate_scale(p.ggj) + +function layer_parameter_columns(params::AbstractVector{CCCParameters}) + return [layer_parameter_columns([p.ggj for p in params]); + "Q_star" => (; values=Float64[p.Qstar for p in params], long_name="scaled electron diamagnetic frequency Q* = −n·ω*ₑ/Q₀ per surface", units="1"); + "M_perp_ratio" => (; values=Float64[p.mperp for p in params], long_name="perpendicular-transport weight M⊥/M per surface", units="1")] +end diff --git a/src/InnerLayer/GGJ/GGJ.jl b/src/InnerLayer/GGJ/GGJ.jl index eacb5eaf86..79b096f176 100644 --- a/src/InnerLayer/GGJ/GGJ.jl +++ b/src/InnerLayer/GGJ/GGJ.jl @@ -36,6 +36,7 @@ import ..InnerLayerResponse import ..solve_inner, ..solve_inner_profile using ..ToroidalLayer import ..ToroidalLayer: inner_Q, rescale_delta, profile_conversions, power_exponents +import ..ToroidalLayer: layer_parameters, mercier_di, growth_rate_scale, layer_parameter_columns 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. @@ -84,6 +85,8 @@ include("Equations.jl") include("Reference.jl") include("Shooting.jl") +layer_parameters(::GGJModel, ggj::GGJParameters, equil, surface) = ggj + """ solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...) -> InnerLayerResponse diff --git a/src/InnerLayer/GGJ/GGJParameters.jl b/src/InnerLayer/GGJ/GGJParameters.jl index 9747d372a7..c121f4b284 100644 --- a/src/InnerLayer/GGJ/GGJParameters.jl +++ b/src/InnerLayer/GGJ/GGJParameters.jl @@ -117,3 +117,13 @@ profile_conversions(p::GGJParameters) = (; dψdx=x0(p) / p.v1, rescale=(p.v1 / x # 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)) + +growth_rate_scale(p::GGJParameters) = q0(p) + +function layer_parameter_columns(params::AbstractVector{GGJParameters}) + coefficient(name) = name => (; values=Float64[getfield(p, Symbol(name)) for p in params], long_name="Glasser-Greene-Johnson coefficient $name per surface", units="1") + return [coefficient.(("E", "F", "G", "H", "K", "M"))..., + "tau_A" => (; values=Float64[p.taua for p in params], long_name="Alfvén time τ_A per surface", units="s"), + "tau_R" => (; values=Float64[p.taur for p in params], long_name="resistive diffusion time τ_R per surface", units="s"), + "dVdpsi" => (; values=Float64[p.v1 for p in params], long_name="dV/dψ_N at each surface", units="m^3")] +end diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 9ea7f15207..b034edcae0 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -9,26 +9,28 @@ module InnerLayer using LinearAlgebra using StaticArrays -using FastInterpolations: cubic_interp +using FastInterpolations: cubic_interp, DerivOp 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: layer_parameters, check_layer_inputs, growth_rate_scale, layer_parameter_columns 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 import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo, toroidal_kpar include("LayerPlasma.jl") +# CCC reads the electron diamagnetic frequency from LayerPlasma. +include("CCC/CCC.jl") +import .CCC: CCCModel, CCCParameters """ model_name(model) -> String @@ -40,11 +42,12 @@ 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 layer_eta, layer_rho, layer_omega_E, layer_omega_star_e export ToroidalLayer, ToroidalLayerModel, ToroidalLayerParameters export GGJ, GGJModel, GGJParameters export CCC, CCCModel, CCCParameters export mercier_di, mercier_dr, inner_Q, rescale_delta +export layer_parameters, check_layer_inputs, growth_rate_scale, layer_parameter_columns export solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear diff --git a/src/InnerLayer/LayerPlasma.jl b/src/InnerLayer/LayerPlasma.jl index f8dd77105d..dbcc371b28 100644 --- a/src/InnerLayer/LayerPlasma.jl +++ b/src/InnerLayer/LayerPlasma.jl @@ -55,3 +55,16 @@ function layer_omega_E(equil, surfaces) spline = SLAYER._slayer_profiles(data, nothing).profiles.omega_E return Float64[spline(s.psifac) for s in surfaces] end + +""" + layer_omega_star_e(equil, surfaces) -> Vector{Float64} + +Electron diamagnetic frequency per unit n, (2π/χ₁)(T_e n_e′/n_e + T_e′) in rad/s with T_e in eV and ′ = d/dψ_N, at each +rational surface of `surfaces`, from the equilibrium's kinetic data; the SLAYER layer's ω*ₑ divided by n. +""" +function layer_omega_star_e(equil, surfaces) + profiles = SLAYER._slayer_profiles(_layer_kinetic_data(equil, "electron diamagnetic frequency"), nothing).profiles + chi1 = 2π * equil.psio + return Float64[(2π / chi1) * (profiles.T_e(ψ) * profiles.n_e(ψ; deriv=DerivOp(1)) / profiles.n_e(ψ) + profiles.T_e(ψ; deriv=DerivOp(1))) + for ψ in (s.psifac for s in surfaces)] +end diff --git a/src/InnerLayer/ToroidalLayer/Ray.jl b/src/InnerLayer/ToroidalLayer/Ray.jl index 6585cfb3ee..4e60cf7afa 100644 --- a/src/InnerLayer/ToroidalLayer/Ray.jl +++ b/src/InnerLayer/ToroidalLayer/Ray.jl @@ -425,6 +425,13 @@ function asymptotic_profile(params::ToroidalLayerParameters, res::RaySolveResult return (; s=collect(float.(srange)), x=cis(res.θ) .* srange, Ψ=Ψ, Ξ=Ξ, Υ=Υ) end +# 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 + """ delta_convergence(params, Q; verbose=true, kwargs...) -> (; Δ, spread, base, table) @@ -433,13 +440,6 @@ 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, 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 = [ diff --git a/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl b/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl index 3f291c5a04..8171d115b4 100644 --- a/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl +++ b/src/InnerLayer/ToroidalLayer/ToroidalLayer.jl @@ -42,7 +42,12 @@ subtype only through this interface, so each model keeps whatever numbers its eq 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. + physical growth rate, the outer region and the layer variables; + - [`mercier_di`](@ref), [`growth_rate_scale`](@ref) and [`layer_parameter_columns`](@ref) - what the tearing and + matching stages read. + +A model builds its parameters with [`layer_parameters`](@ref) and checks its own inputs with +[`check_layer_inputs`](@ref). """ abstract type ToroidalLayerParameters <: InnerLayerParameters end @@ -155,6 +160,42 @@ Sum of the fast exponential rates every Hermite cell resolves, alongside the par """ function fast_rates end +""" + layer_parameters(model, ggj, equil, surface) -> ToroidalLayerParameters + +The per-surface parameters `model` solves on, from the surface's GGJ coefficients and times `ggj` (a +`GGJParameters`), the equilibrium `equil` and the singular surface `surface`. +""" +function layer_parameters end + +""" + check_layer_inputs(model, nsurfaces) -> nothing + +Throws if the inputs `model` carries do not fit a problem with `nsurfaces` rational surfaces; a no-op by default. +""" +check_layer_inputs(::ToroidalLayerModel, ::Integer) = nothing + +""" + mercier_di(params) -> Float64 + +Mercier interchange index D_I of the surface; the layer needs D_I < 0. +""" +function mercier_di end + +""" + growth_rate_scale(params) -> Float64 + +The rate Q₀ in 1/s that scales the physical growth rate to the layer's, Q = γ/Q₀. +""" +function growth_rate_scale end + +""" + layer_parameter_columns(params::AbstractVector) -> Vector{Pair{String,NamedTuple}} + +The per-surface parameters as output columns, `name => (; values, long_name, units)`, one value per surface. +""" +function layer_parameter_columns end + include("InnerAsymptotics.jl") include("Galerkin.jl") include("Ray.jl") @@ -164,6 +205,7 @@ 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 layer_parameters, check_layer_inputs, mercier_di, growth_rate_scale, layer_parameter_columns 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/src/Tearing/Dispersion/Dispersion.jl b/src/Tearing/Dispersion/Dispersion.jl index a8df34e52d..4f85c1ebd2 100644 --- a/src/Tearing/Dispersion/Dispersion.jl +++ b/src/Tearing/Dispersion/Dispersion.jl @@ -12,7 +12,7 @@ using LinearAlgebra using StaticArrays using ..InnerLayer -using ..InnerLayer: InnerLayerModel, layer_response, GGJModel, GGJParameters, +using ..InnerLayer: InnerLayerModel, layer_response, ToroidalLayerModel, ToroidalLayerParameters, SLAYERModel, SLAYERParameters include("SurfaceCoupling.jl") diff --git a/src/Tearing/Dispersion/GGJGrowthRates.jl b/src/Tearing/Dispersion/GGJGrowthRates.jl index 4d80a97608..37b1f4029e 100644 --- a/src/Tearing/Dispersion/GGJGrowthRates.jl +++ b/src/Tearing/Dispersion/GGJGrowthRates.jl @@ -9,7 +9,7 @@ GGJDispersion(mc::MultiSurfaceCouplingFull) The GGJ full-matching determinant on the scan plane, `g(Q) = mc(−i·Q·Q₀)`, where `Q₀ = S^(−1/3)/τ_A` -of the reference surface `mc.surfaces[mc.ref_idx]`. Every surface must use a `GGJModel`. +of the reference surface `mc.surfaces[mc.ref_idx]`. Every surface must use a `ToroidalLayerModel`. """ struct GGJDispersion{M<:MultiSurfaceCouplingFull} mc::M @@ -17,9 +17,9 @@ struct GGJDispersion{M<:MultiSurfaceCouplingFull} end function GGJDispersion(mc::MultiSurfaceCouplingFull) - all(sc -> sc.model isa GGJModel, mc.surfaces) || - throw(ArgumentError("GGJDispersion: every surface must use a GGJModel")) - return GGJDispersion(mc, InnerLayer.GGJ.q0(mc.surfaces[mc.ref_idx].params)) + all(sc -> sc.model isa ToroidalLayerModel, mc.surfaces) || + throw(ArgumentError("GGJDispersion: every surface must use a ToroidalLayerModel")) + return GGJDispersion(mc, InnerLayer.growth_rate_scale(mc.surfaces[mc.ref_idx].params)) end (g::GGJDispersion)(Q::Number) = g.mc(-im * ComplexF64(Q) * g.Q0) diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index ae076bd0ef..97d82a57a0 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -36,7 +36,7 @@ end # Doppler-shift a scanned value into the layer's own frequency variable: real for SLAYER's Q, imaginary for GGJ's s. layer_argument(sc::SurfaceCoupling, Q::Number) = ComplexF64(Q) + sc.q_shift -layer_argument(sc::SurfaceCoupling{<:GGJModel}, s::Number) = ComplexF64(s) + im * sc.q_shift +layer_argument(sc::SurfaceCoupling{<:ToroidalLayerModel}, s::Number) = ComplexF64(s) + im * sc.q_shift """ surface_coupling(model::SLAYERModel, params::SLAYERParameters, @@ -54,10 +54,10 @@ function surface_coupling(model::SLAYERModel, params::SLAYERParameters, end """ - surface_coupling(model::GGJModel, params::GGJParameters, + surface_coupling(model::ToroidalLayerModel, params::ToroidalLayerParameters, dp_diag::Number; omega_E::Real=0.0, n::Integer=0) -> SurfaceCoupling -GGJ convenience constructor. `scale = 1` (GGJ's `rescale_delta` already returns Δ in outer units) and `tauk = 1`. +Toroidal-layer (GGJ, CCC) convenience constructor. `scale = 1` (GGJ's `rescale_delta` already returns Δ in outer units) and `tauk = 1`. `omega_E` is the surface's E×B rotation per unit toroidal mode number, in rad/s, and `n` the toroidal mode number (required when `omega_E` is nonzero). GGJ is a plasma-frame layer (Johnson & Greene 1967, Plasma Phys. 9, 611, @@ -66,7 +66,7 @@ Appendix B), so a rotating surface sees the Doppler-shifted growth rate `s + i· There is no `dc` keyword: the GGJ residual already includes the interchange channel, so a Δ_crit offset would double-count it; `dc` is 0. """ -function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number; omega_E::Real=0.0, n::Integer=0) +function surface_coupling(model::ToroidalLayerModel, params::ToroidalLayerParameters, dp_diag::Number; omega_E::Real=0.0, n::Integer=0) iszero(omega_E) || n >= 1 || throw(ArgumentError("surface_coupling: E×B rotation omega_E=$(omega_E) rad/s needs the toroidal mode number n ≥ 1; got n=$(n)")) return SurfaceCoupling(model, params, ComplexF64(dp_diag), diff --git a/src/Tearing/GGJTearing.jl b/src/Tearing/GGJTearing.jl index c6780aac95..30b3678600 100644 --- a/src/Tearing/GGJTearing.jl +++ b/src/Tearing/GGJTearing.jl @@ -1,6 +1,6 @@ # GGJTearing.jl # -# `solve(TearingProblem, GGJModel)`: tearing growth rates from the full GGJ matching +# `solve(TearingProblem, model::ToroidalLayerModel)` (GGJ, CCC): tearing growth rates from the full matching # determinant (Dispersion.ggj_growth_rates), each surface on its own when uncoupled, # the leading `msing_max` surfaces together when coupled. @@ -51,10 +51,11 @@ end # The SLAYER-only physics fields must stay at their defaults, so a GGJ solve never silently ignores one. const _SLAYER_ONLY_FIELDS = (:critical_resonant_field,) -function _check_ggj_inputs(prob::TearingProblem, config::TearingConfig) +function _check_ggj_inputs(prob::TearingProblem, config::TearingConfig, model::ToroidalLayerModel) defaults = TearingProblem(prob.ffs) for key in _SLAYER_ONLY_FIELDS - getfield(prob, key) == getfield(defaults, key) || throw(ArgumentError("TearingProblem: `$key` is a SLAYER input; GGJModel does not use it")) + getfield(prob, key) == getfield(defaults, key) || + throw(ArgumentError("TearingProblem: `$key` is a SLAYER input; $(nameof(typeof(model))) does not use it")) end prob.gamma > 0 && isfinite(prob.gamma) || throw(ArgumentError("TearingProblem: gamma (ratio of specific heats) must be positive and finite")) config.gamma_range === nothing && throw(ArgumentError("GGJ tearing needs TearingConfig(; gamma_range=(γ_min, γ_max)) in 1/s")) @@ -64,19 +65,21 @@ function _check_ggj_inputs(prob::TearingProblem, config::TearingConfig) return nothing end -function CommonSolve.solve(prob::TearingProblem, model::GGJModel; config::TearingConfig=TearingConfig()) +function CommonSolve.solve(prob::TearingProblem, model::ToroidalLayerModel; config::TearingConfig=TearingConfig()) validate(prob, config) - _check_ggj_inputs(prob, config) + _check_ggj_inputs(prob, config, model) ffs = prob.ffs surfaces = ffs.surfaces isempty(surfaces) && return empty_tearing_result(prob, config, model) raw = _ggj_delta_prime_raw(ffs) count = length(surfaces) + InnerLayer.check_layer_inputs(model, count) eta, rho = _ggj_layer_inputs(prob) - params = [ggj_parameters(s, ffs.equil; eta=eta[k], rho=rho[k], gamma=prob.gamma, ising=k) for (k, s) in enumerate(surfaces)] + params = [InnerLayer.layer_parameters(model, ggj_parameters(s, ffs.equil; eta=eta[k], rho=rho[k], gamma=prob.gamma, ising=k), ffs.equil, s) + for (k, s) in enumerate(surfaces)] labels = ["$(only(s.m))/$(only(s.n))" for s in surfaces] # The layer needs D_I < 0 at every surface; check them all before any scan runs. - unstable = [labels[k] for (k, p) in enumerate(params) if InnerLayer.GGJ.mercier_di(p) >= 0] + unstable = [labels[k] for (k, p) in enumerate(params) if InnerLayer.mercier_di(p) >= 0] isempty(unstable) || throw(ArgumentError("GGJ tearing: the layer needs Mercier-stable surfaces (D_I < 0); $(join(unstable, ", ")) have D_I ≥ 0")) coupled = prob.coupling_mode === :coupled @@ -114,7 +117,7 @@ function CommonSolve.solve(prob::TearingProblem, model::GGJModel; config::Tearin for (k, e) in enumerate(extractions) :no_root in e.warning_flags && continue omega[k] += n * Ω_shift[k] - Q_root[k] += n * Ω_shift[k] / InnerLayer.GGJ.q0(params[k]) + Q_root[k] += n * Ω_shift[k] / InnerLayer.growth_rate_scale(params[k]) end end return TearingResult(prob, true, config, model, params, diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl index dbf1ed5b9f..37e27d4f9b 100644 --- a/src/Tearing/HDF5Output.jl +++ b/src/Tearing/HDF5Output.jl @@ -6,7 +6,7 @@ import ..Utilities.HDF5Annotations: replace_group!, write_ragged! # Writer for a `TearingResult`. Layout under the caller's parent group (inputs live only under Input/): # # Tearing/ -# ├── PerSurface/ -- layer parameters (SLAYER or GGJ set), Δ' matrix, omega_E, q_shift +# ├── PerSurface/ -- the layer model's parameters, Δ' matrix, omega_E, q_shift # ├── Roots/ -- Q_root, omega, gamma, no_root (lab frame) # ├── LayerWidths/ -- resistive layer thickness per surface # ├── Diagnostics/ -- ValidRoots, Poles, FilteredRoots (+ *Winding), flat-plus-offsets ragged; @@ -25,8 +25,8 @@ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::TearingResult) g = replace_group!(parent, "Tearing") g["enabled"] = Int(result.enabled) - # Which inner-layer model produced PerSurface/: the SLAYER and GGJ branches write - # disjoint field sets, so readers must not have to infer it from the schema. + # Which inner-layer model produced PerSurface/: each model writes its own field set, + # so readers must not have to infer it from the schema. attrs(g)["layer_model"] = result.enabled ? _layer_model_token(eltype(result.params)) : "none" attrs(g)["inner_model"] = result.model === nothing ? "none" : InnerLayer.model_name(result.model) @@ -37,7 +37,7 @@ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, end _write_per_surface!(g, result.params, result.dp_matrix) - if result.model isa GGJModel && result.problem !== nothing + if result.model isa ToroidalLayerModel && result.problem !== nothing g["PerSurface/Delta_prime_raw"] = result.problem.ffs.delta_prime.raw end # E×B rotation resolved per surface, and the Doppler offset it applies to the inner-layer Q. @@ -77,10 +77,9 @@ function _write_critical_resonant_field!(g, crf::CriticalResonantFieldResult) return nothing end -# Token recorded in the Tearing group's layer_model attribute; keyed by the -# per-surface parameter type since the SLAYER and GGJ branches write disjoint fields. -_layer_model_token(::Type{SLAYERParameters}) = "slayer" -_layer_model_token(::Type{GGJParameters}) = "ggj" +# Token recorded in the Tearing group's layer_model attribute, from the per-surface parameter type: +# SLAYERParameters → "slayer", GGJParameters → "ggj", CCCParameters → "ccc". +_layer_model_token(::Type{P}) where {P<:InnerLayerParameters} = lowercase(replace(String(nameof(P)), "Parameters" => "")) # Metadata table for the Tearing group (paths relative to it); ragged Diagnostics # subgroups and Scan/Surface_ groups are annotated by iteration below. @@ -123,14 +122,6 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/k_ref" => (; long_name="reference-length ratio K = r_s·(dψ_N/dr) at each surface", dims=("surface",)), "PerSurface/alpha_mercier" => (; long_name="Mercier Frobenius exponent α = √(−D_I) at each surface (Glasser-Greene-Johnson 1975 Eq. 48)", dims=("surface",)), "PerSurface/delta_prime_conversion" => (; long_name="ψ_N → r_s reference-length factor K^(2α) applied to each Δ' diagonal", dims=("surface",)), - "PerSurface/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",)), - "PerSurface/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)), - "PerSurface/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)), - "PerSurface/H" => (; long_name="Glasser-Greene-Johnson coefficient H per surface", dims=("surface",)), - "PerSurface/K" => (; long_name="Glasser-Greene-Johnson coefficient K per surface", dims=("surface",)), - "PerSurface/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)), - "PerSurface/tau_A" => (; long_name="Alfvén time τ_A per surface (GGJ layer parameters)", units="s", dims=("surface",)), - "PerSurface/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces, ψ_N-referenced (GGJ path; identical to SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col") @@ -265,18 +256,15 @@ function _write_per_surface!(g, params::AbstractVector{SLAYERParameters}, return nothing end -# GGJ per-surface parameters carry the geometric coefficients (E,F,G,H,K,M) -# and resistive/Alfvén times rather than the SLAYER dimensionless set. -function _write_per_surface!(g, params::AbstractVector{GGJParameters}, - dp_matrix::Matrix{ComplexF64}) +# A toroidal layer model's per-surface parameters are its own output columns (InnerLayer.layer_parameter_columns). +function _write_per_surface!(g, params::AbstractVector{<:ToroidalLayerParameters}, dp_matrix::Matrix{ComplexF64}) ps = create_group(g, "PerSurface") - ps["rational_index"] = Int[p.ising for p in params] - for fname in (:E, :F, :G, :H, :K, :M) - ps[String(fname)] = Float64[getfield(p, fname) for p in params] + ps["rational_index"] = collect(1:length(params)) + columns = InnerLayer.layer_parameter_columns(params) + for (name, column) in columns + ps[name] = column.values end - ps["tau_A"] = Float64[p.taua for p in params] - ps["tau_R"] = Float64[p.taur for p in params] - ps["dVdpsi"] = Float64[p.v1 for p in params] + Utilities.HDF5Annotations.annotate!(ps, [name => (; column.long_name, column.units, dims=("surface",)) for (name, column) in columns]) ps["Delta_prime_matrix"] = dp_matrix return nothing end diff --git a/src/Tearing/SLAYERTearing.jl b/src/Tearing/SLAYERTearing.jl index 7913ac1cb7..cb7b61283b 100644 --- a/src/Tearing/SLAYERTearing.jl +++ b/src/Tearing/SLAYERTearing.jl @@ -48,7 +48,7 @@ _surface_label(p::SLAYERParameters) = "m=$(p.m), n=$(p.n)" _surface_label(p::InnerLayerParameters) = "ising=$(getfield(p, :ising))" # Is this a GGJ inner-layer model? Its growth rates come from `ggj_growth_rates`, not this runner. -_is_ggj(::GGJModel) = true +_is_ggj(::ToroidalLayerModel) = true _is_ggj(::Any) = false # Scan dispatch diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 0055cb2739..86ef17759e 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -8,7 +8,7 @@ # resonant field of error-field penetration # TearingProblem / TearingConfig / TearingResult -- the public solve, # `solve(TearingProblem(ffs; ...), model; config)`, with one -# method per inner-layer model (SLAYERModel, GGJModel) +# method for SLAYERModel and one for every ToroidalLayerModel (GGJ, CCC) # # `InnerLayer` itself lives at the top level (`src/InnerLayer/`) and is loaded # before `ForceFreeStates`, which depends on it for the matched-Δ′ Galerkin @@ -32,7 +32,7 @@ import ..ForceFreeStates: layer_overlap_scan import ..InnerLayer as InnerLayer using .InnerLayer: InnerLayerParameters, InnerLayerResponse, solve_inner, layer_response, SLAYERModel, SLAYERParameters, build_slayer_inputs, - GGJModel, GGJParameters, + ToroidalLayerModel, ToroidalLayerParameters, LayerWidths, slayer_layer_thickness using .InnerLayer.SLAYER: _slayer_profiles, _layer_zeff, _layer_ion_mass, _layer_mass_per_electron, _build_resistivity_model diff --git a/src/Tearing/TearingProblem.jl b/src/Tearing/TearingProblem.jl index 1d1b5a75d0..eb0ca50b1b 100644 --- a/src/Tearing/TearingProblem.jl +++ b/src/Tearing/TearingProblem.jl @@ -10,8 +10,8 @@ GGJ layer inputs; the SLAYER layer closures belong to [`SLAYERModel`](@ref) and [`TearingConfig`](@ref). The plasma itself (profiles, species, Z_eff and transport defaults, confinement time and resistivity closure) is the equilibrium's kinetic data, `ffs.equil.kinetic` (see [`KineticProfileData`](@ref)), so every layer model sees one plasma. -Solve with `solve(problem, GGJModel(); config)` or `solve(problem, SLAYERModel(); config)`. -GGJ needs the raw Δ′ matrix (Riccati with `vac_flag=true`, or Galerkin with `rpec_flag=true`). +Solve with `solve(problem, model; config)` for a toroidal layer model (`GGJModel()`, `CCCModel()`) or `SLAYERModel()`. +A toroidal layer model needs the raw Δ′ matrix (Riccati with `vac_flag=true`, or Galerkin with `rpec_flag=true`). Both models report lab-frame roots from each surface's E×B rotation. ## Fields diff --git a/test/runtests_innerlayer.jl b/test/runtests_innerlayer.jl index a058c97c4f..2a4b410725 100644 --- a/test/runtests_innerlayer.jl +++ b/test/runtests_innerlayer.jl @@ -252,6 +252,14 @@ end @test_throws ArgumentError IL.CCCModel(; solver=:galerkin) @test_throws ArgumentError IL.CCCModel(; nx=10) @test_throws ArgumentError IL.CCCParameters(q4, 0.1, -1.0) + # Per-surface inputs live on the model: a number for every surface or one value per surface. + m = IL.CCCModel(; Qstar=[0.1, 0.2], mperp=1.2) + @test m.Qstar == [0.1, 0.2] && m.mperp === 1.2 && m.options == (;) + @test_throws ArgumentError IL.CCCModel(; mperp=0.0) + @test_throws "Qstar has length 2, expected 3" IL.check_layer_inputs(m, 3) + @test IL.check_layer_inputs(m, 2) === nothing && IL.check_layer_inputs(IL.GGJModel(), 3) === nothing + @test IL.layer_parameters(m, IL.GGJParameters(; E=q4.E, F=q4.F, G=q4.G, H=q4.H, K=q4.K, taua=q4.taua, taur=q4.taur, ising=2), nothing, nothing).Qstar == 0.2 + @test [c.first for c in IL.layer_parameter_columns([ccc(; Qstar=0.3)])] == ["E", "F", "G", "H", "K", "M", "tau_A", "tau_R", "dVdpsi", "Q_star", "M_perp_ratio"] # 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 diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index e3f205cbc7..32748bcc5d 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -1,7 +1,7 @@ @testset "Tearing: TearingProblem + TearingConfig + HDF5 output" begin using GeneralizedPerturbedEquilibrium using GeneralizedPerturbedEquilibrium.InnerLayer - using GeneralizedPerturbedEquilibrium: GGJModel, SLAYERModel, TearingProblem, TearingConfig, solve + using GeneralizedPerturbedEquilibrium: GGJModel, CCCModel, SLAYERModel, TearingProblem, TearingConfig, solve using GeneralizedPerturbedEquilibrium.Tearing using GeneralizedPerturbedEquilibrium.Tearing.Dispersion using HDF5 @@ -192,6 +192,25 @@ TearingProblem(ffs; eta=1e-7, rho=1e-7, critical_resonant_field=GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantFieldControl()), GGJModel(); config=TearingConfig(; gamma_range=(600.0, 2600.0), omega_range=(-1500.0, 1500.0))) + + # CCC tearing on the same solve: at Q* = 0, M⊥/M = 1 it is the GGJ layer, through its own Wasow basis and splitting solve. + box = TearingConfig(; gamma_range=(600.0, 2600.0), omega_range=(-1500.0, 1500.0), nre=4, nim=4, amr_passes=1) + ccc = solve(TearingProblem(ffs; eta=1e-7, rho=1e-7), CCCModel(; Qstar=0.0); config=box) + @test ccc.model === CCCModel(; Qstar=0.0) && eltype(ccc.params) <: CCCParameters + @test ccc.gamma ≈ still.gamma rtol = 1e-6 + @test ccc.omega ≈ still.omega atol = 1e-3 * maximum(abs, still.gamma) + mktemp() do path, io + close(io) + h5open(f -> save_to_hdf5(f, ccc), path, "w") + h5open(path, "r") do f + @test read(attributes(f["Tearing"])["layer_model"]) == "ccc" + @test read(f["Tearing/PerSurface/Q_star"]) == zeros(length(ffs.surfaces)) + @test read(f["Tearing/PerSurface/M_perp_ratio"]) == ones(length(ffs.surfaces)) + end + end + # Without kinetic data CCC needs Q* given, and a per-surface Q* needs one value per surface. + @test_throws "give Qstar, or load kinetic profiles" solve(TearingProblem(ffs; eta=1e-7, rho=1e-7), CCCModel(); config=box) + @test_throws "Qstar has length" solve(TearingProblem(ffs; eta=1e-7, rho=1e-7), CCCModel(; Qstar=zeros(length(ffs.surfaces) + 1)); config=box) end end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 0b03b4c306..7b98785f91 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -404,6 +404,32 @@ using HDF5 @test MatchProblem(kgal; rotation=1.0).rotation == fill(1.0, n) end + @testset "CCC matching" begin + ric = solve(equil, riccati; nn=1, wall=wall, ffs_kwargs...) + n = length(ric.surfaces) + layer = (eta=1e-7, rho=1e-7, rotation=1.0) + ggj = solve(MatchProblem(ric; layer...), GGJModel()) + # Q* = 0, M⊥/M = 1 is the GGJ layer, through CCC's own Wasow basis and splitting solve. + @test solve(MatchProblem(ric; layer...), CCCModel(; Qstar=0.0)).bpen ≈ ggj.bpen rtol = 1e-6 + drift = solve(MatchProblem(ric; layer...), CCCModel(; Qstar=0.05, mperp=1.2)) + @test drift.closure === :matched && all(isfinite, drift.bpen) && drift.bpen != ggj.bpen + @test drift.match.model === CCCModel(; Qstar=0.05, mperp=1.2) && all(p -> p.Qstar == 0.05 && p.mperp == 1.2, drift.match.inner_params) + # Without kinetic data CCC needs Q* given. + @test_throws "give Qstar, or load kinetic profiles" solve(MatchProblem(ric; layer...), CCCModel()) + # From kinetic data, ω*ₑ = (2π/χ₁)(T_e n_e′/n_e + T_e′) per unit n, so a flat density and T_e = 1 keV·(1 − ψ/2) give + # (2π/χ₁)·(−500 eV); Q* = −n·ω*ₑ/Q₀ at each surface. + psi = collect(range(0.0, 1.0; length=64)) + keq = deepcopy(ric.equil) + attach_kinetic_profiles!(keq, KineticProfileData(; psi, n_e=fill(5e19, 64), n_i=fill(5e19, 64), T_e=1e3 .* (1 .- psi ./ 2), + T_i=fill(1e3, 64), omega_E=fill(2π * 1e3, 64))) + ωse = GPEC.InnerLayer.layer_omega_star_e(keq, ric.surfaces) + @test ωse ≈ fill(-500 / keq.psio, n) rtol = 1e-6 + ggj1 = FFS.ggj_parameters(ric.surfaces[1], keq; eta=1e-7, rho=1e-7, ising=1) + p = GPEC.InnerLayer.layer_parameters(CCCModel(), ggj1, keq, ric.surfaces[1]) + @test p.ggj === ggj1 && p.mperp == 1.0 + @test p.Qstar ≈ -ric.nlow * ωse[1] / GPEC.InnerLayer.growth_rate_scale(ggj1) rtol = 1e-12 + end + @testset "matched inner profiles: opt-in, and parity columns agree across inner models" begin mktempdir() do dir gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, ffs_kwargs...) @@ -510,6 +536,21 @@ using HDF5 end end + @testset "deck tearing under [Tearing] with the CCC layer" begin + t = GPEC.tearing_from_deck( + Dict{String,Any}("Tearing" => Dict{String,Any}("inner_model" => "CCC"), + "CCC" => Dict{String,Any}("cfl" => 0.1, "Qstar" => [0.05, 0.1], "mperp" => 1.2)) + ) + @test t.model isa CCCModel{:ray} && t.model.options == (; cfl=0.1) && t.model.Qstar == [0.05, 0.1] && t.model.mperp == 1.2 + # The model's deck section round-trips. + @test GPEC.toroidal_layer_model(GPEC.deck_sections(t.model), "CCC").Qstar == [0.05, 0.1] + # A [CCC] section no stage uses is refused. + mktempdir() do dir + deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("CCC" => Dict{String,Any}("cfl" => 0.1))) + @test_throws "[CCC] is set but neither" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) + end + end + @testset "deck matching under [Match]" begin layer = (eta=1e-7, rho=1e-7, rotation=1.0) # The deck's [Match] section for `n` surfaces and the [GGJ] solver. @@ -545,6 +586,21 @@ using HDF5 @test api.match.cout ≈ reference.match.cout rtol = 1e-12 end + # The CCC layer under [Match], named by inner_model with its options in [CCC]. + mktempdir() do dir + ric = solve(equil, riccati; nn=1, wall=wall, ffs_kwargs...) + n = length(ric.surfaces) + api = solve(MatchProblem(ric; layer...), CCCModel(; Qstar=fill(0.05, n), mperp=fill(1.2, n))) + sections = Dict{String,Any}("CCC" => Dict{String,Any}("solver" => "ray", "Qstar" => fill(0.05, n), "mperp" => fill(1.2, n)), + "Match" => Dict{String,Any}("inner_model" => "CCC", "eta" => fill(layer.eta, n), "rho" => fill(layer.rho, n), + "rotation" => fill(layer.rotation, n))) + reference = GPEC.main([_stage_deck(dir; ffs_overrides=Dict{String,Any}("integrator" => "Riccati"), extra_sections=sections)]).ffs + + @test reference.closure === :matched && reference.match.model isa CCCModel{:ray} && reference.match.model.Qstar == fill(0.05, n) + @test api.bpen ≈ reference.bpen rtol = 1e-12 + @test api.match.cout ≈ reference.match.cout rtol = 1e-12 + end + # A forward solve has no Δ′ to match against, so the deck stops at parse time instead of after the solve. mktempdir() do dir deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("Match" => Dict{String,Any}("ideal" => true))) From 979f829d895e5fad5915ab83f979d9c7866755aa Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 13 Nov 2026 10:35:45 +0100 Subject: [PATCH 3/3] =?UTF-8?q?InnerLayer.CCC=20-=20BUGFIX=20-=20Take=20Q*?= =?UTF-8?q?=20=3D=20n=C2=B7=CF=89*=E2=82=91/Q=E2=82=80=20so=20the=20CCC=20?= =?UTF-8?q?layer=20co-rotates=20with=20the=20electrons=20as=20SLAYER's=20d?= =?UTF-8?q?oes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With s = γ − iω, Kim's Ohm's law Q + iQ* vanishes at ω = Q*Q₀ in the E×B frame, so Q* = ω*ₑ/Q₀ with SLAYER's ω*ₑ = n·layer_omega_star_e. Q* from kinetic data had the opposite sign and put the CCC electron resonance at n·Ω_E − ω*ₑ instead of n·Ω_E + ω*ₑ. A test pins both layers' resonances at the same lab frequency, with and without E×B rotation. The SurfaceCoupling docstring now states s = γ − iω like the rest of the toroidal-layer path. --- docs/src/inner_layer.md | 2 +- src/InnerLayer/CCC/CCC.jl | 4 +-- src/InnerLayer/CCC/CCCParameters.jl | 7 ++-- src/Tearing/Dispersion/SurfaceCoupling.jl | 7 ++-- test/runtests_solve_api.jl | 44 ++++++++++++++++++++--- 5 files changed, 51 insertions(+), 13 deletions(-) diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index 34dfbea164..c95ad3e890 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -145,7 +145,7 @@ p = GGJ.q4_surface_benchmark() # GGJParameters for the q = 4 benchmark surf [`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. -In a deck, `[Tearing]` or `[Match]` takes `inner_model = "CCC"`, and `[CCC]` holds the solver options with `Qstar` and `mperp`, each a scalar or one value per surface, as the `CCCModel` keywords of the same names; without `Qstar`, Q* = −n·ω*ₑ/Q₀ comes from the kinetic data ([`layer_omega_star_e`](@ref InnerLayer.layer_omega_star_e)), and `mperp` defaults to 1. +In a deck, `[Tearing]` or `[Match]` takes `inner_model = "CCC"`, and `[CCC]` holds the solver options with `Qstar` and `mperp`, each a scalar or one value per surface, as the `CCCModel` keywords of the same names; without `Qstar`, Q* = n·ω*ₑ/Q₀ comes from the kinetic data ([`layer_omega_star_e`](@ref InnerLayer.layer_omega_star_e)), and `mperp` defaults to 1. 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 diff --git a/src/InnerLayer/CCC/CCC.jl b/src/InnerLayer/CCC/CCC.jl index f8b6adefbd..2c617ef54e 100644 --- a/src/InnerLayer/CCC/CCC.jl +++ b/src/InnerLayer/CCC/CCC.jl @@ -40,7 +40,7 @@ the parameters the tearing and matching stages build ([`layer_parameters`](@ref) - `options::NamedTuple` - Backend keywords forwarded to every solve; call-site keywords win. - `Qstar` - Scaled electron diamagnetic frequency Q* at every surface (a number), at each surface core to edge - (a vector), or `nothing` to take Q* = −n·ω*ₑ/Q₀ from the equilibrium's kinetic data + (a vector), or `nothing` to take Q* = n·ω*ₑ/Q₀ from the equilibrium's kinetic data ([`layer_omega_star_e`](@ref)). - `mperp` - Perpendicular-transport weight M⊥/M at every surface (a number) or at each surface core to edge (a vector); positive. @@ -100,7 +100,7 @@ function layer_parameters(model::CCCModel, ggj::GGJParameters, equil, surface) Qstar = if model.Qstar === nothing equil.kinetic === nothing && throw(ArgumentError("CCC layer: give Qstar, or load kinetic profiles (kinetic_file in [Equilibrium], or attach_kinetic_profiles!)")) - -only(surface.n) * only(layer_omega_star_e(equil, (surface,))) / q0(ggj) + only(surface.n) * only(layer_omega_star_e(equil, (surface,))) / q0(ggj) else _at_surface(model.Qstar, ggj.ising) end diff --git a/src/InnerLayer/CCC/CCCParameters.jl b/src/InnerLayer/CCC/CCCParameters.jl index 3dd50f22ea..65f2b3a6d5 100644 --- a/src/InnerLayer/CCC/CCCParameters.jl +++ b/src/InnerLayer/CCC/CCCParameters.jl @@ -13,8 +13,9 @@ of Kim, Cho & Choi, Phys. Fluids 31, 2659 (1988), Eqs. (1a)–(1c), at one ratio - `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*). + - `Qstar::Float64` - Electron diamagnetic frequency of the mode scaled like Q, Q* = ω*ₑ/Q₀, with ω*ₑ = n times + [`layer_omega_star_e`](@ref), SLAYER's. With the growth rate s = γ − iω, Kim's Ohm's law carries Q + iQ*, which + vanishes at ω = ω*ₑ in the E×B frame, where the layer co-rotates with the electrons. - `mperp::Float64` - Perpendicular-transport weight M⊥/M of the QΥ'' term; must be positive. """ struct CCCParameters <: ToroidalLayerParameters @@ -44,6 +45,6 @@ growth_rate_scale(p::CCCParameters) = growth_rate_scale(p.ggj) function layer_parameter_columns(params::AbstractVector{CCCParameters}) return [layer_parameter_columns([p.ggj for p in params]); - "Q_star" => (; values=Float64[p.Qstar for p in params], long_name="scaled electron diamagnetic frequency Q* = −n·ω*ₑ/Q₀ per surface", units="1"); + "Q_star" => (; values=Float64[p.Qstar for p in params], long_name="scaled electron diamagnetic frequency Q* = n·ω*ₑ/Q₀ per surface", units="1"); "M_perp_ratio" => (; values=Float64[p.mperp for p in params], long_name="perpendicular-transport weight M⊥/M per surface", units="1")] end diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index 97d82a57a0..7a7b6a40e7 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -13,8 +13,8 @@ r(Q) = dp_diag - scale * layer_response(model, params, Q + q_shift).tearing - dc `q_shift` is the surface's E×B Doppler shift, a real frequency offset in the layer's own frequency variable, zero by default. SLAYER's argument is the normalized Q, so the shift is added -to Q; GGJ's argument is the growth rate s = γ + iω, so the shift is added to ω and enters as -`s + i·q_shift`. The SLAYER and coupled GGJ Tearing runs set it. +to Q; a toroidal layer's argument is the growth rate s = γ − iω, so the layer frequency ω − q_shift enters as +`s + i·q_shift`. The SLAYER and coupled toroidal-layer Tearing runs set it. A root of `sc` in the complex `Q` plane is an uncoupled tearing eigenvalue at this surface (tearing channel only); coupled eigenvalues, including the interchange channel, come from `MultiSurfaceCoupling`. @@ -61,7 +61,8 @@ Toroidal-layer (GGJ, CCC) convenience constructor. `scale = 1` (GGJ's `rescale_d `omega_E` is the surface's E×B rotation per unit toroidal mode number, in rad/s, and `n` the toroidal mode number (required when `omega_E` is nonzero). GGJ is a plasma-frame layer (Johnson & Greene 1967, Plasma Phys. 9, 611, -Appendix B), so a rotating surface sees the Doppler-shifted growth rate `s + i·n·omega_E`, and `q_shift = n·omega_E`. +Appendix B), so with s = γ − iω a rotating surface sees the Doppler-shifted growth rate `s + i·n·omega_E`, the +frequency ω − n·omega_E, and `q_shift = n·omega_E`. There is no `dc` keyword: the GGJ residual already includes the interchange channel, so a Δ_crit offset would double-count it; `dc` is 0. diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 7b98785f91..7f05ac82d3 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -417,17 +417,53 @@ using HDF5 # Without kinetic data CCC needs Q* given. @test_throws "give Qstar, or load kinetic profiles" solve(MatchProblem(ric; layer...), CCCModel()) # From kinetic data, ω*ₑ = (2π/χ₁)(T_e n_e′/n_e + T_e′) per unit n, so a flat density and T_e = 1 keV·(1 − ψ/2) give - # (2π/χ₁)·(−500 eV); Q* = −n·ω*ₑ/Q₀ at each surface. + # (2π/χ₁)·(−500 eV); Q* = n·ω*ₑ/Q₀ at each surface. psi = collect(range(0.0, 1.0; length=64)) keq = deepcopy(ric.equil) - attach_kinetic_profiles!(keq, KineticProfileData(; psi, n_e=fill(5e19, 64), n_i=fill(5e19, 64), T_e=1e3 .* (1 .- psi ./ 2), - T_i=fill(1e3, 64), omega_E=fill(2π * 1e3, 64))) + attach_kinetic_profiles!( + keq, + KineticProfileData(; psi, n_e=fill(5e19, 64), n_i=fill(5e19, 64), T_e=1e3 .* (1 .- psi ./ 2), + T_i=1e3 .* (1 .- psi ./ 2), omega_E=fill(2π * 1e3, 64)) + ) ωse = GPEC.InnerLayer.layer_omega_star_e(keq, ric.surfaces) @test ωse ≈ fill(-500 / keq.psio, n) rtol = 1e-6 ggj1 = FFS.ggj_parameters(ric.surfaces[1], keq; eta=1e-7, rho=1e-7, ising=1) p = GPEC.InnerLayer.layer_parameters(CCCModel(), ggj1, keq, ric.surfaces[1]) @test p.ggj === ggj1 && p.mperp == 1.0 - @test p.Qstar ≈ -ric.nlow * ωse[1] / GPEC.InnerLayer.growth_rate_scale(ggj1) rtol = 1e-12 + @test p.Qstar ≈ ric.nlow * ωse[1] / GPEC.InnerLayer.growth_rate_scale(ggj1) rtol = 1e-12 + end + + @testset "CCC and SLAYER put the electron-frame resonance at the same lab frequency" begin + # With s = γ − iω, SLAYER's layer co-rotates with the electrons where ĝ − i·Q_e = 0, and CCC's where Q + i·Q* = 0 (Kim + # Eq. 1a). Both drift terms come from the same kinetic data, so at γ = 0 both zeros sit at ω_lab = n·Ω_E + ω*ₑ. + ric = solve(equil, riccati; nn=1, wall=wall, ffs_kwargs...) + psi = collect(range(0.0, 1.0; length=64)) + keq = deepcopy(ric.equil) + attach_kinetic_profiles!( + keq, + KineticProfileData(; psi, n_e=fill(5e19, 64), n_i=fill(5e19, 64), T_e=1e3 .* (1 .- psi ./ 2), + T_i=1e3 .* (1 .- psi ./ 2), omega_E=fill(2π * 1e3, 64)) + ) + TR, DS, IL = GPEC.Tearing, GPEC.Tearing.Dispersion, GPEC.InnerLayer + sing = ric.surfaces[1] + nn = only(sing.n) + slayer_model = SLAYERModel(; P_perp_model=:chi_perp_e) + sp = only(TR._tearing_params(slayer_model, keq, [sing], TR._slayer_profiles(keq.kinetic, slayer_model))) + cp = IL.layer_parameters(CCCModel(), FFS.ggj_parameters(sing, keq; eta=1e-7, rho=1e-7, ising=1), keq, sing) + ω_star_e = -sp.Q_e / sp.tauk # SLAYER's ω*ₑ (Q_e = −τ_k·ω*ₑ), the full mode's + @test ω_star_e ≈ nn * IL.layer_omega_star_e(keq, [sing])[1] rtol = 1e-10 + # A linear residual's zero in the real lab frequency ω. + zero_of(f) = (f0 = f(0.0); -real(f0 / (f(1.0) - f0))) + for Ω_E in (0.0, 2π * 3e3) + # SLAYER scans Q_lab = τ_k·ω_lab (γ = 0); its layer sees ĝ = −i·(Q_lab + q_shift). + slayer = DS.surface_coupling(slayer_model, sp, 0.0; q_shift=TR._q_shift(sp, Ω_E)) + ω_slayer = zero_of(ω -> -im * DS.layer_argument(slayer, sp.tauk * ω) - im * sp.Q_e) + # CCC's determinant takes the lab-frame s = −i·ω (γ = 0); its layer sees Q = inner_Q(s + i·n·Ω_E). + ccc = DS.surface_coupling(CCCModel(), cp, 0.0; omega_E=Ω_E, n=nn) + ω_ccc = zero_of(ω -> IL.inner_Q(cp, DS.layer_argument(ccc, -im * ω)) + im * cp.Qstar) + @test ω_slayer ≈ nn * Ω_E + ω_star_e rtol = 1e-10 + @test ω_ccc ≈ nn * Ω_E + ω_star_e rtol = 1e-10 + end end @testset "matched inner profiles: opt-in, and parity columns agree across inner models" begin