Tearing - FEATURE - Add critical resonant field calculations to the tearing workflow

Project.toml

@@ -1,7 +1,7 @@
 name = "GeneralizedPerturbedEquilibrium"
 uuid = "462872dd-e066-4d2e-b993-6468b5239634"
 license = "MIT"
-authors = ["Nikolas Logan <ncl2128@columbia.edu>", "Jong-Kyu Park <jkpark@snu.ac.kr>", "Matthew Pharr <m.pharr@protonmail.com>", "Jacob Halpern <jmh2363@columbia.edu>", "Rithik Banerjee <rb3736@columbia.edu>", "Jaebeom Cho <aspire1019@snu.ac.kr>", "Daniel Burgess <dab2245@columbia.edu>", "Min-Gu Yoo <yoom@fusion.gat.com>"]
+authors = ["Nikolas Logan <ncl2128@columbia.edu>", "Jong-Kyu Park <jkpark@snu.ac.kr>", "Matthew Pharr <m.pharr@protonmail.com>", "Jacob Halpern <jmh2363@columbia.edu>", "Rithik Banerjee <rb3736@columbia.edu>", "Jaebeom Cho <aspire1019@snu.ac.kr>", "Daniel Burgess <dab2245@columbia.edu>", "Min-Gu Yoo <yoom@fusion.gat.com>", "Evan Bursch <emb2333@columbia.edu>"]
 version = "0.1.0"
 
 [deps]

docs/development/hdf5-conventions.md

@@ -46,7 +46,7 @@ Top level (11 groups):
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |
 | `KineticForces/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) |
-| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_<k>/` |
+| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_<k>/`, `CriticalResonantField/` (+ `Scan/Surface_<k>/`) |
 | `SurfaceGeometries/` | `{Plasma,Wall}/{x,y,z}` point clouds |
 
 Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegration/` — the third integrator backend slot alongside `ForwardIntegration` and `GalerkinIntegration`.

docs/src/citations.md

@@ -124,6 +124,14 @@ The basis for the presently implemented in Fortran SLAYER code which computes th
 
 ---
 
+> A. Cole and R. Fitzpatrick, "Drift-magnetohydrodynamical model of error-field penetration in tokamak plasmas,"
+> *Physics of Plasmas* **13**, 032503 (2006).
+> DOI: [10.1063/1.2178167](https://doi.org/10.1063/1.2178167)
+
+The torque-balance error-field penetration threshold (Eqs. 61–62) behind the critical resonant field computed in `Tearing.CriticalResonantField`.
+
+---
+
 > A. Burgess et al., "Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer,"
 > Preprint (2026).
 

docs/src/inner_layer.md

@@ -3,7 +3,8 @@
 The `Tearing` module groups the resistive tearing-mode analysis stack:
 `InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and
 SLAYER models), `Dispersion` (physics-agnostic complex-plane scan and
-contour-intersection root extraction), and `Runner` (user-facing TOML
+contour-intersection root extraction), `CriticalResonantField` (torque-balance
+error-field penetration threshold), and `Runner` (user-facing TOML
 configuration, profile loading, and HDF5 output).
 
 ## Layer Inputs
@@ -278,6 +279,60 @@ res = solve_ray(p, GGJ.inner_Q(p, γ))
 res.Δ, res.resid, res.bc_cond
 ```
 
+## Critical resonant field
+
+With `[SLAYER.CriticalResonantField] enabled = true`, each SLAYER surface also gets the
+critical resonant field for error-field penetration, from the steady-state torque balance of
+Cole and Fitzpatrick, Phys. Plasmas **13**, 032503 (2006). The electromagnetic torque on the
+layer balances the viscous restoring torque when (Cole Eq. 61)
+
+```math
+\frac{\operatorname{Im}\hat\Delta(Q)}{|\alpha + \hat\Delta(Q)|^2}
+  = \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,(b_r/B_\phi)^2},
+\qquad
+\hat\kappa = \left[\frac{2}{s}\right]^2 \int_{r_s}^{a} \frac{\mu(r_s)}{\mu(r)}\,\frac{dr}{r},
+```
+
+and balance is lost above (Cole Eq. 62)
+
+```math
+\left(\frac{b_r}{B_\phi}\right)^2_\mathrm{crit}
+  = \max_Q \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,\operatorname{Im}[-1/(\alpha + \hat\Delta(Q))]} .
+```
+
+- **Inputs.** P is each layer's `P_tor` (so χ_φ from the kinetic file, or the scalar
+  `chi_tor`), S its Lundquist number, and Q0 = τ_k·n·ω_E the natural E×B rotation of the
+  m/n mode (Cole's ω0; the diamagnetic drifts enter through Q_e and Q_i in Δ̂).
+- **Approximations.** As in Fortran `gslayer.f`, the viscosity integral is taken as 1/2
+  (κ̂ = 2/s²) and α = 10⁻² stands in for S^(−1/3)(−r_s Δ'_s) ≪ 1 (b_crit moves about 4 %
+  over 10⁻⁴–10⁻¹).
+- **Q axis.** Q is on Cole's axis, with the diamagnetic poles at Q_e and Q_i; `solve_inner`
+  mirrors the real axis, so the scan uses conj(Δ(−Q)).
+- **Maximum.** The scan window follows `gslayer.f` (between Q0 and the Q_e pole), and the
+  largest positive interior local maximum is taken, skipping any within 0.02 of Q_e or Q_i
+  (Fortran SLAYER's pole-regularization radius).
+
+Outputs go to `Tearing/CriticalResonantField/` (`br_crit` = b_r/B_φ, `q_peak`, `q0`, `p_phi`,
+and with `store_scan` the per-surface `Scan/Surface_<k>/{Q, balance, Delta}`).
+
+### Benchmark against Fortran
+
+On the DIII-D-like example (2/1–6/1, n = 1, P = 1, κ̂ = 2/s², Q ∈ [−5, 5] with 20001
+points), both codes were given Fortran's normalized inputs.
+
+- **Torque-balance routine.** Fed Fortran's own Δ(Q) table with the pole exclusion narrowed to
+  one grid step, the Julia routine reproduces Fortran's b_crit and Q_peak to all printed digits
+  on every surface.
+- **Pole handling.** Fortran's global maximum sits within 0.003 of the Q_e pole on 2/1, 3/1,
+  4/1 and 6/1. Skipping the pole region moves the threshold to the torque-balance maximum
+  between Q_e and Q0 (blue dots below), lowering the Fortran-table b_crit by 6–88 %.
+- **Layer model.** Fortran solves Park's drift-MHD layer and Julia Fitzpatrick's, so Δ(Q)
+  differs by a median 7–78 % between Q_i and Q_e, and Julia's b_crit is 82–89 % of the
+  Fortran-table value on 2/1, 3/1, 4/1 and 6/1. The 5/1 edge surface is dominated by a
+  zero of the electromagnetic torque in both codes and is not a usable threshold.
+
+![Fortran vs Julia inner-layer Δ(Q) and torque balance](assets/crf_fortran_benchmark.png)
+
 ## API Reference
 
 ### InnerLayer
@@ -304,6 +359,12 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER]
 Modules = [GeneralizedPerturbedEquilibrium.Dispersion]
 ```
 
+## CriticalResonantField
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField]
+```
+
 ## Runner
 
 ```@autodocs

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -90,3 +90,11 @@ Q_re_range = [-2.0, 2.0]       # Scan box in the normalized Q plane, Re(Q) axis
 Q_im_range = [-0.5, 3.0]       # Scan box in the normalized Q plane, Im(Q) axis
 nre        = 41                # Grid resolution along the Re(Q) axis
 nim        = 31                # Grid resolution along the Im(Q) axis
+
+[SLAYER.CriticalResonantField]
+# Critical resonant field b_r/B_φ for error-field penetration at each SLAYER surface, from the
+# Cole-Fitzpatrick 2006 torque balance with P = P_tor of each layer. Qmin/Qmax may be set to
+# override the per-surface scan window derived from Q0, Q_e and Q_i.
+enabled    = true      # run the analysis after SLAYER
+n          = 2000      # number of Q samples per surface
+store_scan = false     # write the per-surface Q, torque balance and Δ samples to gpec.h5

regression-harness/cases/diiid_slayer_n1.toml

@@ -152,6 +152,24 @@ label = "SLAYER no_root flags [2/1,3/1,4/1]"
 noise_threshold = 0
 order = 34
 
+# Critical resonant field (torque balance) for the inner three surfaces, matching
+# the root pins above.
+[quantities.crf_br_crit]
+h5path = "Tearing/CriticalResonantField/br_crit"
+type = "real_vector"
+extract = "first_3"
+label = "CRF b_r/B_φ crit [2/1,3/1,4/1]"
+noise_threshold = 1e-9
+order = 40
+
+[quantities.crf_q_peak]
+h5path = "Tearing/CriticalResonantField/q_peak"
+type = "real_vector"
+extract = "first_3"
+label = "CRF Q_peak [2/1,3/1,4/1]"
+noise_threshold = 1e-6
+order = 41
+
 # Settings (catches accidental config drift)
 [quantities.slayer_enabled]
 h5path = "Tearing/enabled"

src/GeneralizedPerturbedEquilibrium.jl

@@ -48,6 +48,7 @@ export Tearing
 # Backward-compat top-level aliases so callers can still reach these
 # directly; the canonical nested path is `Tearing.{Dispersion,Runner}`.
 import .Tearing.Dispersion as Dispersion
+import .Tearing.CriticalResonantField as CriticalResonantField
 import .Tearing.Runner as Runner
 export Dispersion, Runner
 

src/InnerLayer/SLAYER/Riccati.jl

@@ -225,6 +225,7 @@ function solve_inner(::SLAYERModel{:fitzpatrick},
     # growth-rate roots are identical to those of Δ̂_s(ĝ); it only reflects
     # the off-root residual surface, fixing the integration-path orientation
     # so |Δ| contours are consistent across the scan plane.
+    # Real-Q callers see a mirrored axis; the critical resonant field maps back via conj(Δ(−Q)).
     Q_c = im * conj(ComplexF64(Q))
 
     # Boundary condition at p_start

src/Tearing/CriticalResonantField/CriticalResonantField.jl

@@ -0,0 +1,18 @@
+# CriticalResonantField.jl
+#
+# Uses the inner-layer delta from InnerLayer.jl to find the critical resonant field
+# required for mode penetration at each rational surface. This can then be compared
+# to the actual resonant field from the perturbed equilibrium to determine which modes
+# are predicted to penetrate.
+
+module CriticalResonantField
+
+using Printf
+
+using ..InnerLayer: InnerLayerModel, solve_inner
+
+include("TorqueBalance.jl")
+
+export TorqueBalance, cole_delta, torque_balance_value, torque_balance_window, torque_balance_scan
+
+end # module CriticalResonantField

src/Tearing/CriticalResonantField/TorqueBalance.jl

@@ -0,0 +1,113 @@
+# TorqueBalance.jl
+#
+# Error-field penetration threshold from the steady-state torque balance at a
+# rational surface (Cole and Fitzpatrick, Phys. Plasmas 13, 032503 (2006), Sec. IV).
+# The electromagnetic and viscous torques balance when (Cole Eq. 61)
+#
+#   Im[Δ(Q)] / |α + Δ(Q)|² = 2·P·(Q0 − Q) / (S·κ̂·(b_r/B_φ)²),
+#
+# with κ̂ = [2/s(r_s)]² ∫_{r_s}^{a} [μ(r_s)/μ(r)] dr/r. Balance is lost when no Q
+# satisfies this, so the critical field is (Cole Eq. 62)
+#
+#   (b_r/B_φ)²_crit = max_Q 2·P·(Q0 − Q) / (S·κ̂·Im[−1/(α + Δ(Q))]).
+#
+# As in Fortran gslayer.f, the viscosity integral is taken as 1/2 (κ̂ = 2/s²) and α
+# is kept finite at 1e-2. Q is on Cole's axis: the electron and ion diamagnetic
+# frequencies sit at Q = Q_e and Q = Q_i (see `cole_delta`).
+
+# Finite stand-in for α = S^(-1/3)·(−r_s Δ'_s) ≪ 1; b_crit moves ~4% over 1e-4 to 1e-1.
+const TORQUE_BALANCE_ALPHA = 1e-2
+# Half-width around Q_e, Q_i where the layer response is pole-dominated (Fortran SLAYER layfac).
+const TORQUE_BALANCE_POLE_WIDTH = 0.02
+
+"""
+    TorqueBalance{M<:InnerLayerModel,P}
+
+Torque-balance inputs at one rational surface.
+
+## Fields
+
+  - `model`     -- inner-layer model passed to `solve_inner`
+  - `params`    -- that model's layer parameters at the surface
+  - `Q0`        -- normalized natural E×B rotation of the m/n mode, τ_k·n·ω_E
+  - `P`         -- magnetic Prandtl number τ_R/τ_V (the layer's P_φ)
+  - `lu`        -- Lundquist number S
+  - `kappa_hat` -- Cole's slab-to-tokamak factor κ̂ (Eq. 61)
+"""
+struct TorqueBalance{M<:InnerLayerModel,P}
+    model::M
+    params::P
+    Q0::Float64
+    P::Float64
+    lu::Float64
+    kappa_hat::Float64
+end
+
+"""
+    cole_delta(model, params, Q::Real) -> ComplexF64
+
+Inner-layer Δ at real frequency `Q` on Cole's axis. `solve_inner` evaluates the layer
+at `i·conj(Q)`, which mirrors the real axis, so Cole's Δ(Q) is `conj(Δ_solve_inner(−Q))`.
+"""
+cole_delta(model::InnerLayerModel, params, Q::Real) =
+    conj(solve_inner(model, params, ComplexF64(-Q)).tearing)
+
+"""
+    torque_balance_value(tb::TorqueBalance, Q::Real) -> (bal, Δ)
+
+Right-hand side of Cole Eq. 62 before the maximum, `bal = 2·P·(Q0 − Q) / Im[−1/(α + Δ)]`,
+and the layer `Δ(Q)` it used.
+"""
+function torque_balance_value(tb::TorqueBalance, Q::Real)
+    Δ = cole_delta(tb.model, tb.params, Q)
+    jxb = -imag(1.0 / (Δ + TORQUE_BALANCE_ALPHA))
+    return 2.0 * tb.P * (tb.Q0 - Q) / jxb, Δ
+end
+
+"""
+    torque_balance_window(Q0, Q_e, Q_i) -> (Qmin, Qmax)
+
+Q range bracketing the torque-balance branch between the natural rotation `Q0` and the
+electron diamagnetic pole `Q_e` (the rule of Fortran gslayer.f).
+"""
+function torque_balance_window(Q0::Real, Q_e::Real, Q_i::Real)
+    Q0 > Q_e && return (1.05 * Q_e, 2.0 * Q0)
+    Qmin = Q0 > 0 ? 0.8 * Q_i : 1.5 * min(Q0, Q_i)
+    return (Qmin, 0.95 * Q_e)
+end
+
+"""
+    torque_balance_scan(tb::TorqueBalance; Qmin=nothing, Qmax=nothing, n=2000)
+        -> (Qs, bal, Qpeak, br_crit, idx_peak, Δs)
+
+Sample `torque_balance_value` on `n` uniform real Q points and return the critical
+normalized field `br_crit = b_r/B_φ` (Cole Eq. 62) at the largest positive interior local
+maximum of `bal`, located at `Qpeak = Qs[idx_peak]`. A maximum within
+`TORQUE_BALANCE_POLE_WIDTH` (or one grid step, if larger) of the diamagnetic poles `Q_e`,
+`Q_i` is discarded as a pole artifact. `Qmin`/`Qmax`
+default to `torque_balance_window`. Returns NaN for `Qpeak`, `br_crit` and 0 for
+`idx_peak` when no valid maximum is found.
+"""
+function torque_balance_scan(tb::TorqueBalance; Qmin=nothing, Qmax=nothing, n::Integer=2000)
+    p = tb.params
+    wmin, wmax = torque_balance_window(tb.Q0, p.Q_e, p.Q_i)
+    Qs = range(something(Qmin, wmin), something(Qmax, wmax); length=n)
+    out = [torque_balance_value(tb, Q) for Q in Qs]
+    bal = first.(out)
+    Δs = last.(out)
+
+    w = max(step(Qs), TORQUE_BALANCE_POLE_WIDTH)
+    near_pole(Q) = abs(Q - p.Q_e) <= w || abs(Q - p.Q_i) <= w
+    peaks = [i for i in 2:(n-1) if isfinite(bal[i]) && bal[i] > 0 &&
+             bal[i-1] < bal[i] && bal[i] >= bal[i+1] && !near_pole(Qs[i])]
+    if isempty(peaks)
+        @warn "CriticalResonantField: no positive torque-balance maximum on " *
+              "Q ∈ [$(first(Qs)), $(last(Qs))] at the $(p.m)/$(p.n) surface."
+        return Qs, bal, NaN, NaN, 0, Δs
+    end
+    idx = peaks[argmax(bal[peaks])]
+    br_crit = sqrt(bal[idx] / (tb.lu * tb.kappa_hat))
+    @info @sprintf("CriticalResonantField: %d/%d surface b_r/B_φ crit = %.3e at Q = %.3f",
+        p.m, p.n, br_crit, Qs[idx])
+    return Qs, bal, Qs[idx], br_crit, idx, Δs
+end

src/Tearing/Runner/Control.jl

@@ -5,6 +5,46 @@
 # constructor or by parsing the `[SLAYER]` (and nested `[SLAYER.*]`)
 # section(s) of a `gpec.toml`.
 
+"""
+    CriticalResonantFieldControl
+
+Configuration for the critical resonant field (torque-balance) analysis, read from the
+`[SLAYER.CriticalResonantField]` TOML section via `critical_resonant_field_control_from_toml`
+or built with the `@kwdef` keyword constructor. The magnetic Prandtl number is the per-surface
+`P_tor` of the SLAYER layer parameters.
+
+## Fields
+
+  - `enabled`    -- run the analysis after SLAYER
+  - `Qmin`, `Qmax` -- override the real-Q scan window; `nothing` derives it per surface from Q0, Q_e and Q_i
+  - `n`          -- number of Q samples per surface
+  - `store_scan` -- write the per-surface Q, torque balance and Δ samples to `gpec.h5`
+"""
+@kwdef struct CriticalResonantFieldControl
+    enabled::Bool = false
+    Qmin::Union{Nothing,Float64} = nothing
+    Qmax::Union{Nothing,Float64} = nothing
+    n::Int = 2000
+    store_scan::Bool = false
+end
+
+"""
+    critical_resonant_field_control_from_toml(section::AbstractDict) -> CriticalResonantFieldControl
+
+Parse a `[SLAYER.CriticalResonantField]` TOML section. Unknown keys raise an error.
+"""
+function critical_resonant_field_control_from_toml(section::AbstractDict)
+    field_names = Set(String.(fieldnames(CriticalResonantFieldControl)))
+    unknown = [k for k in keys(section) if !(k in field_names)]
+    isempty(unknown) ||
+        throw(
+            ArgumentError("critical_resonant_field_control_from_toml: unknown keys " *
+                          "$(unknown) in [SLAYER.CriticalResonantField]. Known: " *
+                          "$(sort(collect(field_names))).")
+        )
+    return CriticalResonantFieldControl(; (Symbol(k) => v for (k, v) in section)...)
+end
+
 """
     SLAYERControl
 
@@ -100,6 +140,11 @@ there is one consistent interface for resistive and kinetic profiles.
 
   - `store_scan`  -- write the full Q/Δ scan grid to HDF5. `false` by
     default to keep the output file small.
+
+# Critical resonant field
+
+  - `critical_resonant_field` -- `CriticalResonantFieldControl` from the
+    `[SLAYER.CriticalResonantField]` subsection; requires `inner_model = :slayer_fitzpatrick`
 """
 @kwdef struct SLAYERControl
     enabled::Bool = false
@@ -162,6 +207,8 @@ there is one consistent interface for resistive and kinetic profiles.
     profile_group::String = "/"
 
     store_scan::Bool = false
+
+    critical_resonant_field::CriticalResonantFieldControl = CriticalResonantFieldControl()
 end
 
 const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_shooting, :ggj_galerkin)
@@ -198,6 +245,10 @@ function validate(ctrl::SLAYERControl)
         throw(ArgumentError("SLAYERControl: nre and nim must both be ≥ 2"))
     ctrl.amr_passes >= 0 ||
         throw(ArgumentError("SLAYERControl: amr_passes must be ≥ 0"))
+    !ctrl.critical_resonant_field.enabled || ctrl.inner_model === :slayer_fitzpatrick ||
+        throw(ArgumentError("SLAYERControl: CriticalResonantField requires inner_model=:slayer_fitzpatrick"))
+    ctrl.critical_resonant_field.n >= 3 ||
+        throw(ArgumentError("SLAYERControl: CriticalResonantField n must be ≥ 3"))
     return ctrl
 end
 
@@ -233,6 +284,8 @@ function slayer_control_from_toml(section::AbstractDict)
             haskey(v, "pole_threshold") && (flat["pole_threshold"] = v["pole_threshold"])
             haskey(v, "filter_above_poles") && (flat["filter_above_poles"] = v["filter_above_poles"])
             haskey(v, "filter_outside_re") && (flat["filter_outside_re"] = v["filter_outside_re"])
+        elseif k == "CriticalResonantField" && v isa AbstractDict
+            flat["critical_resonant_field"] = critical_resonant_field_control_from_toml(v)
         else
             flat[k] = v
         end