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

docs/development/hdf5-conventions.md

@@ -46,7 +46,7 @@ Top level (11 groups):
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, `TorqueResponse/` (kinetic forward solves: ψ-resolved `T_xe`, `T_source`), control-surface spectra |
 | `KineticForces/` | `<method>/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `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/development/references.md

@@ -92,6 +92,15 @@ The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are
   - Published: Nuclear Fusion (2025), doi 10.1088/1741-4326/ae4fdd (open access: arXiv:2511.07666)
   - Describes: Diffusive-resistive layer width (Eq. 100) and the resistive-layer-overlap truncation `0 < Ψ < 1 − ε_c` (Sect. 5.9)
 
+- **Cole & Fitzpatrick (2006)**: "Drift-magnetohydrodynamical model of error-field penetration in tokamak plasmas"
+  - Location: not in `docs/resources/`
+  - Published: Physics of Plasmas 13, 032503 (2006), doi 10.1063/1.2178167
+  - Describes: Steady-state torque balance at a rational surface (Eq. 61) and the critical resonant field for error-field penetration (Eq. 62), used by `Tearing.CriticalResonantField`
+- **Fitzpatrick (1993)**: "Interaction of tearing modes with external structures in cylindrical geometry"
+  - Location: not in `docs/resources/`
+  - Published: Nuclear Fusion 33, 1049 (1993)
+  - Describes: Bifurcated rotating states and the fold at which an error field penetrates (Sec. 7.4), the threshold definition used by `Tearing.CriticalResonantField`
+
 ## Resistive MHD Stability Analysis (Future Work)
 
 GPEC will eventually implement resistive MHD stability analysis based on:

docs/src/citations.md

@@ -152,6 +152,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`.
+
+---
+
 > D. A. Burgess, N. C. Logan, J.-K. Park, C. J. Hansen, and C. Paz-Soldan, "Tearing stability prediction combining toroidal calculations with a two-fluid slab layer,"
 > *Plasma Physics and Controlled Fusion* **68**, 065017 (2026).
 > DOI: [10.1088/1361-6587/ae7156](https://doi.org/10.1088/1361-6587/ae7156)

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
@@ -201,6 +202,81 @@ res = solve_ray(p, GGJ.inner_Q(p, γ))
 res.Δ, res.breaks, res.bc_cond      # res.Δ = (Δ₊, Δ₋) = (interchange, tearing)
 ```
 
+## 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 Δ̂). Ω_E is the
+  kinetic file's `omega_E` or the `omega_E_kHz` override; the analysis needs it, and a surface
+  without one is skipped.
+- **Approximations.** As in Fortran `gslayer.f`, the viscosity integral is taken as 1/2
+  (κ̂ = 2/s²).
+- **Stability α.** Each surface uses its own α = S^(−1/3)(Δ_crit − r_s Δ'), so α + Δ̂ is the
+  dispersion relation's residual and the torque balance shares its poles with SLAYER's free
+  modes (Fortran fixes α = 10⁻²). Cole and Fitzpatrick have no Δ_crit; including it treats the
+  curvature offset as independent of frequency. It is zero unless `D_c_type` is set. Cole's
+  threshold needs a tearing-stable surface, so a surface with α ≤ 0 is skipped with a warning.
+- **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)).
+- **Threshold.** Penetration is associated with a saddle-node bifurcation (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4):
+  as the field grows, Q moves from Q0 toward the electron fluid at Q_e, where the torque
+  vanishes, and the threshold is the first resolved saddle-node bifurcation on that path, above
+  which the shielded branch followed from Q0 no longer exists. Both the rise to and fall from
+  an interior maximum must exceed 0.1 % of its value within a contiguous retained interval.
+  This excludes monotonic clipped windows and small numerical extrema; excluded samples do
+  not connect separate intervals. A later, higher maximum does not replace the first resolved
+  saddle-node bifurcation. Points where |Δ̂| ≤ α, the
+  unshielded zone next to Q_e, are not candidates, so the balance spike there is never taken.
+  With no resolved saddle-node bifurcation no discrete threshold is reported (NaN). Fortran takes the global maximum on a
+  window instead, which picks that spike. The search assumes Q_i ≤ 0 ≤ Q_e; a surface where
+  that fails is skipped with a warning.
+
+Outputs go to `Tearing/CriticalResonantField/` (`br_crit` = b_r/B_φ, `br_crit_tesla`, `Q_peak`, `Q0`, with P in
+`PerSurface/P_tor`, and with `store_scan` the per-surface `Scan/Surface_<k>/{Q, balance, Delta}`).
+
+`br_crit` is normalized by the toroidal field, as in Cole and Fitzpatrick (2006, Eqs. 61–68) and
+Fitzpatrick, Phys. Plasmas **30**, 092512 (2023), Eq. 120. `br_crit_tesla` = `br_crit`·|B0|, with
+the layer's B0, is for comparison with `SingularCoupling/resonant_area_weighted_field`, the total
+resonant field with the ideal plasma response, which Park et al., Phys. Rev. Lett. **99**, 195003
+(2007) identify as the field that opens islands. The comparison is approximate: the theory's b_r is
+a cylindrical field at the layer, and no reference maps it onto the area-weighted toroidal field.
+`penetrated_area_weighted_field` already includes the layer response and is not the matching quantity.
+
+### 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 and Fortran's fixed α = 10⁻².
+
+- **Torque balance.** Fed Fortran's own Δ(Q) table, the Julia routine reproduces Fortran's
+  balance to all printed digits on every surface, so the thresholds differ only through the
+  selection rule above.
+- **Layer model.** Fortran solves Park's drift-MHD layer and Julia Fitzpatrick's, so Δ(Q)
+  differs by a median 24–78 % between Q_i and Q_e on 2/1, 3/1, 4/1 and 6/1 (7 % on 5/1).
+
+The dots in the figure mark each code's balance maximum between Q_e and Q0 at that fixed α,
+not the fold threshold above.
+
+![Fortran vs Julia inner-layer Δ(Q) and torque balance](assets/crf_fortran_benchmark.png)
+
 ## API Reference
 
 ### InnerLayer
@@ -227,6 +303,12 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER]
 Modules = [GeneralizedPerturbedEquilibrium.Dispersion]
 ```
 
+## CriticalResonantField
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField]
+```
+
 ## Runner
 
 ```@autodocs

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -102,3 +102,12 @@ 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 and Q0 = τ_k·n·Ω_E, so it
+# needs the kinetic file's omega_E (or omega_E_kHz). b_crit is taken at the penetration fold met
+# walking from Q0 toward Q_e; Qmin/Qmax may be set to override the sampled range.
+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,40 @@ 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, and the natural rotation Q0 = τ_k·n·Ω_E it is computed from.
+[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_br_crit_tesla]
+h5path = "Tearing/CriticalResonantField/br_crit_tesla"
+type = "real_vector"
+extract = "first_3"
+label = "CRF b_r crit in tesla [2/1,3/1,4/1]"
+noise_threshold = 1e-9
+order = 43
+
+[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
+
+[quantities.crf_Q0]
+h5path = "Tearing/CriticalResonantField/Q0"
+type = "real_vector"
+extract = "first_3"
+label = "CRF Q0 [2/1,3/1,4/1]"
+noise_threshold = 1e-10
+order = 42
+
 # 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, fold_search_range, torque_balance_scan
+
+end # module CriticalResonantField

src/Tearing/CriticalResonantField/TorqueBalance.jl

@@ -0,0 +1,140 @@
+# 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))]).
+#
+# The viscosity integral is taken as 1/2 (κ̂ = 2/s²). α = S^(-1/3)·(−r_s Δ'_s) is each
+# surface's own value (Cole Sec. IV C), which needs a tearing-stable surface (α > 0).
+# Q is on Cole's axis: the electron and ion diamagnetic frequencies sit at Q = Q_e
+# and Q = Q_i (see `cole_delta`).
+
+
+"""
+    TorqueBalance{M<:InnerLayerModel,T}
+
+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; must carry `Q_e`, `Q_i`, `m`, `n`
+  - `Q0`        -- normalized natural E×B rotation of the m/n mode, τ_k·n·ω_E
+  - `prandtl`   -- magnetic Prandtl number τ_R/τ_V (the layer's P_φ)
+  - `lu`        -- Lundquist number S
+  - `kappa_hat` -- Cole's slab-to-tokamak factor κ̂ (Eq. 61)
+  - `alpha`     -- outer stability α = S^(-1/3)·(−r_s Δ'_s) (Eq. 61), positive on a tearing-stable surface
+"""
+struct TorqueBalance{M<:InnerLayerModel,T}
+    model::M
+    params::T
+    Q0::Float64
+    prandtl::Float64
+    lu::Float64
+    kappa_hat::Float64
+    alpha::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 / (Δ + tb.alpha))
+    return 2.0 * tb.prandtl * (tb.Q0 - Q) / jxb, Δ
+end
+
+# A resolved saddle-node bifurcation requires a rise and fall exceeding this fraction of the maximum.
+const TORQUE_BALANCE_FOLD_DROP = 1e-3
+
+"""
+    fold_search_range(Q0, Q_e, Q_i) -> (Qa, Qb)
+
+The segment the steady state follows as the error field grows: from the natural rotation `Q0`
+toward the electron-fluid point `Q_e`, where the electromagnetic torque vanishes, ending at `Q_i`
+instead if it lies in between (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4).
+"""
+function fold_search_range(Q0::Real, Q_e::Real, Q_i::Real)
+    return (Q0, min(Q0, Q_e) < Q_i < max(Q0, Q_e) ? Q_i : 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 first resolved saddle-node bifurcation on the
+path from `Q0` along `fold_search_range`. Both the rise to and fall from an interior maximum
+must exceed `TORQUE_BALANCE_FOLD_DROP` times its value within a contiguous retained interval.
+Above this threshold the shielded branch followed from `Q0` no longer exists as the field grows.
+Points where |Δ| ≤ α, the unshielded zone next to Q_e, are not on the path. With no resolved
+saddle-node bifurcation there is no discrete threshold. `Qmin`/`Qmax` override the sampled range,
+but the search stays between `Q0` and the end of `fold_search_range`.
+Returns NaN for `Qpeak`, `br_crit` and 0 for `idx_peak`, with a warning, when no saddle-node bifurcation is resolved,
+when `Q0` is not finite, when the surface is tearing-unstable (α ≤ 0), or when the diamagnetic
+ordering Q_i ≤ 0 ≤ Q_e does not hold.
+"""
+function torque_balance_scan(tb::TorqueBalance; Qmin=nothing, Qmax=nothing, n::Integer=2000)
+    p = tb.params
+    skip(why) = (@warn "CriticalResonantField: $(p.m)/$(p.n) surface skipped: $why"; (Float64[], Float64[], NaN, NaN, 0, ComplexF64[]))
+    isfinite(tb.Q0) || return skip("its natural rotation Q0 = $(tb.Q0) is not finite.")
+    tb.alpha > 0 || return skip("it is tearing-unstable or marginal (α = $(tb.alpha)), so there is no penetration threshold.")
+    p.Q_i <= 0 <= p.Q_e && p.Q_i < p.Q_e ||
+        return skip("the fold search assumes Q_i ≤ 0 ≤ Q_e, but Q_e = $(p.Q_e), Q_i = $(p.Q_i).")
+    Qa, Qb = fold_search_range(tb.Q0, p.Q_e, p.Q_i)
+    Qs = range(something(Qmin, Qa), something(Qmax, Qb); length=n + 2)[2:(end-1)]   # open at both ends
+    out = [torque_balance_value(tb, Q) for Q in Qs]
+    bal = first.(out)
+    Δs = last.(out)
+
+    # Follow the first resolved interior maximum without joining samples across excluded intervals.
+    path = [i for i in sort(eachindex(Qs); by=i -> abs(Qs[i] - Qa))
+                  if min(Qa, Qb) < Qs[i] < max(Qa, Qb) && isfinite(bal[i]) && bal[i] > 0 && abs(Δs[i]) > tb.alpha]
+    idx = 0
+    valley = peak = previous = 0
+    for i in path
+        if previous == 0 || abs(i - previous) != 1
+            valley = peak = i
+        elseif bal[i] > bal[peak]
+            peak = i
+        elseif bal[i] < (1 - TORQUE_BALANCE_FOLD_DROP) * bal[peak]
+            if bal[valley] < (1 - TORQUE_BALANCE_FOLD_DROP) * bal[peak]
+                idx = peak
+                break
+            end
+            valley = peak = i
+        elseif bal[i] < bal[valley]
+            valley = peak = i
+        end
+        previous = i
+    end
+    if idx == 0
+        @warn "CriticalResonantField: no penetration saddle-node bifurcation between Q0 = $(tb.Q0) and Q = $Qb at the $(p.m)/$(p.n) " *
+              "surface: no interior maximum with a rise and fall exceeding $(100 * TORQUE_BALANCE_FOLD_DROP) % was resolved " *
+              "on a contiguous sampled interval, so no discrete threshold is resolved."
+        return Qs, bal, NaN, NaN, 0, Δs
+    end
+    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