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

docs/src/citations.md

@@ -144,6 +144,21 @@ Provides an efficient algorithm for computing the full Δ' matrix (coupling betw
 
 ---
 
+> G. S. Cho, E. H. Choi, and D. I. Choi, "Compressional effects on drift-resistive instabilities in general toroidal plasma,"
+> *IEEE Transactions on Plasma Science* **PS-15**, 418 (1987).
+
+Derives the toroidal resistive-layer equations with electron diamagnetic drift, ion acoustic coupling, ion polarization drift and perpendicular resistivity by generalizing the GGJ scheme; `CCCModel` in the `InnerLayer` module solves them.
+
+---
+
+> J.-Y. Kim, G.-S. Cho, and D.-I. Choi, "Numerical study of drift-resistive modes in general toroidal geometry,"
+> *Physics of Fluids* **31**, 2659 (1988).
+> DOI: [10.1063/1.866543](https://doi.org/10.1063/1.866543)
+
+States those equations as Eqs. (1a)–(1c), with the electron diamagnetic frequency Q* and the perpendicular-transport weight M⊥/M, and their large-x exponents (Eq. 2): the form `CCCModel` implements. Its tearing dispersion curves (Figs. 1b, 2b) plot the Fourier-space Δ_o, which its Eq. (10b) relates to the real-space Δ that `CCCModel` returns.
+
+---
+
 > A. H. Glasser and Z. R. Wang, "Asymptotic solutions and convergence studies of the resistive inner region equations,"
 > *Physics of Plasmas* **27**, 012506 (2020).
 > DOI: [10.1063/1.5134999](https://doi.org/10.1063/1.5134999)

docs/src/inner_layer.md

@@ -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 |
@@ -138,6 +139,20 @@ 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.
+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
+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 +173,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).

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/`.

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
 

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

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