InnerLayer.SLAYER - FEATURE! - Derive the toroidal critical-Δ geometric factor in the r_s reference

[d-burg]

Release note

• Audience: users
• Numerical impact: none at defaults (dc_type = "none"); with dc_type = "toroidal" the critical-Δ offset is now computed from the equilibrium rather than from a user-supplied dgeo_val. The HDF5 dataset Tearing/PerSurface/D_geo carries derived values instead of zeros. (harness @ d433522)
• Migration: none. dc_type = "toroidal" no longer errors without dgeo_val; passing dgeo_val still overrides the derivation.

The toroidal critical-Δ of Connor, Ham, Hastie & Liu 2015 (PPCF 57 065001, Eq. 59) can now be used on any equilibrium: its geometric factor is derived per rational surface from the flux-surface averages GPEC already computes, expressed in the same r_s reference length as the rest of the SLAYER layer stack.

Regression report

regress --cases diiid_slayer_n1,diiid_n1_riccati --refs baf674659,d433522fa (#403 head vs this branch, both freshly run, same pinned manifest, julia 1.11.6):

Regression Report: diiid_slayer_n1
Ref 1: baf674659  @ baf674659 (2026-08-20)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
Ref 2: d433522fa  @ d433522fa (2026-09-03)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
------------------------------------------------------------------------------------------
Quantity                            baf674659  d433522fa  Diff               Status       
------------------------------------------------------------------------------------------
SLAYER surface indices              [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER poloidal m                   [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER toroidal n                   [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER minor radius rs              [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER r-based shear                [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER Lundquist S                  [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER D_norm                       [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER P_perp                       [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER tauk                         [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER iota_e                       [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   8.1e-06            OK           
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   3.0e-05            OK           
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   1.217e-01 (0.01%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00            OK           
SLAYER enabled flag                 1          1          0.0e+00            OK           
Runtime (s)                         134.1s     154.4s                        --           
Summary: 1 changed, 14 unchanged
Regression Report: diiid_n1_riccati
Ref 1: baf674659  @ baf674659 (2026-08-20)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
Ref 2: d433522fa  @ d433522fa (2026-09-03)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
----------------------------------------------------------------------------
Quantity                       baf674659      d433522fa      Diff     Status
----------------------------------------------------------------------------
delta prime (BVP diagonal)     [5 elem]       [5 elem]       0.0e+00  OK    
delta prime (raw side-major)   [10 elem]      [10 elem]      0.0e+00  OK    
edge coil response delta_coil  [10 elem]      [10 elem]      0.0e+00  OK    
total energy Re(et[1])         8.038607e-01   8.038607e-01   0.0e+00  OK    
plasma energy Re(ep[1])        -1.344628e+00  -1.344628e+00  0.0e+00  OK    
vacuum energy Re(ev[1])        2.148488e+00   2.148488e+00   0.0e+00  OK    
total energy (all)             [35 elem]      [35 elem]      0.0e+00  OK    
# singular surfaces            5              5              0.0e+00  OK    
singular psi locations         [5 elem]       [5 elem]       0.0e+00  OK    
singular q values              [5 elem]       [5 elem]       0.0e+00  OK    
ODE steps (saved)              51             51             0.0e+00  OK    
ODE steps (total)              1667           1667           0.0e+00  OK    
mpert                          35             35             0.0e+00  OK    
npert                          1              1              0.0e+00  OK    
q0                             1.204212e+00   1.204212e+00   0.0e+00  OK    
q95                            4.781723e+00   4.781723e+00   0.0e+00  OK    
beta_n                         1.372511e+00   1.372511e+00   0.0e+00  OK    
Runtime (s)                    131.7s         129.7s                  --    
Summary: 17 unchanged

Every SLAYER layer input (rs, shear, S, D_norm, P_perp, tauk, iota_e) and the whole Riccati case are bit-identical, as expected: with dc_type = "none" the new code only adds the D_geo diagnostic. The single flagged line, γ at 0.01 % (0.12 Hz) with Q_root within 8e-6 and ω within 3e-5 Hz, is the threaded AMR root-search sample dependence documented on #403 and #418; no input to the dispersion changed. No tolerance or acceptance criterion was altered.

Base branch

This PR is based on and targets bugfix/slayer-dprime-reference-length (#403), because #403 defines the r_s reference length K = r_s·dψ_N/dr that the conversion below uses and rewrites the same region of build_slayer_inputs. Retarget to develop once #403 merges.

Derivation

Connor's Eq. 59, in Hamada coordinates with ψ(V) the toroidal flux, χ(V) the poloidal flux, primes d/dV, ι/2π = χ'/ψ' = 1/q, Λ = ψ'² ι'/2π, α = 2πn/χ':

$$\Delta'_{crit} = \frac{\pi^{3/2}}{2}\left(\frac{\chi_\parallel}{\chi_\perp}\right)^{1/4} V_s \left(\frac{\alpha^2\Lambda^2}{\langle B^2\rangle\langle|\nabla V|^2\rangle}\right)^{1/4}(-D_R)$$

In GPEC quantities (ψ_N normalized poloidal flux, chi1 = 2π·psio = dχ/dψ_N with χ the full poloidal flux, v1 = dV/dψ_N, q1 = dq/dψ_N):

$$\iota' = 2\pi\frac{d(1/q)}{dV} = -\frac{2\pi}{q^2},\frac{q_1}{v_1},\qquad \psi_t' = \frac{d\psi_t}{d\psi_N}\frac{d\psi_N}{dV} = \frac{q,\texttt{chi1}}{v_1},\qquad \Lambda = \left(\frac{q,\texttt{chi1}}{v_1}\right)^2\left(-\frac{q_1}{q^2 v_1}\right),\qquad \alpha = \frac{2\pi n}{\texttt{chi1}/v_1}$$

$$\langle B^2\rangle = \oint B^2,\frac{\texttt{jac}}{V'},d\theta,\qquad \langle|\nabla V|^2\rangle = v_1^2\oint|\nabla\psi_N|^2,\frac{\texttt{jac}}{V'},d\theta$$

Two points where the paper's notation has to be read carefully, both verified against the paper's own large-aspect-ratio limit (Eq. 61):

1. Connor writes ⟨A⟩ = (1/2π)∮A dθ, but with Hamada θ of period 1 this must be the normalized flux-surface average (⟨1⟩ = 1), i.e. the jac/V' weighting above. A literal 1/2π leaves a (2π)^{1/2} mismatch with Eq. 61.
2. Eq. 59 is the jump in dΨ/dY with Y = (V−V_s)/V_s. GPEC's layer stack, the rfitzp critical-Δ, and the Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching #403-converted Δ' are all referenced to x̂ = (r−r_s)/r_s. Because Eq. 59 is derived at H = 0 (Frobenius exponent μ = 1/2), the conversion is linear: multiply by r_s(dV/dr)/V_s = K·v1/V_s. V_s cancels, so no volume integral is needed. Without this conversion Eq. 59 gives exactly half of rfitzp at large aspect ratio (V ∝ r²).

With the conversion the geometric factor reduces at large aspect ratio to √(n s r_s/R₀) with s = (r_s/q)dq/dr, which is exactly the factor inside rfitzp, so the two branches coincide there and the paper's Eq. 61 is recovered after dividing by r_s.

Implementation

• ResistGeometry (src/ForceFreeStates/ResistEval.jl) gains avg_dpsisq = ⟨|∇ψ_N|²⟩, an eighth column of the existing θ-loop.

• toroidal_dgeo (src/InnerLayer/SLAYER/LayerInputs.jl, exported) evaluates

  dgeo = k_ref · v1 · (α²Λ²/(⟨B²⟩ · v1² · ⟨|∇ψ_N|²⟩))^{1/4}
  

  with α and Λ as above. k_ref · v1 is the Y → x̂ reference conversion; this is the only place the coordinate transformation enters, and k_ref is the same K that delta_prime_to_rs_reference applies to Δ'.

• build_slayer_inputs derives dgeo_val whenever the surface carries a ResistGeometry (sing.restype, populated by resist_eval_all!), for every dc_type; only :toroidal consumes it, via the unchanged _solve_dc_tmp formula 0.5·(−D_R)·π^{3/2}·(χ∥/χ⊥)^{1/4}·dgeo. It throws only if dc_type=:toroidal is requested on a surface without a ResistGeometry. An explicit dgeo_val still overrides.

Verification

• Analytic chain Eq. 59 → Eq. 61 → Lutjens offset √2π^{3/2}(−D_R)/w_D (Eq. 65) → the rfitzp code formula: all consistent, pinning the prefactor.
• Large-aspect-ratio limit on TJ circular equilibria: dgeo/√(n s r_s/R₀) = 1.001–1.012 at ε = 0.05, ≤ 1.09 at ε = 0.3 near the edge (unit test at 1 %).
• Scale invariance: dgeo unchanged to 6 digits under B₀ → B₀/2 and (a, R₀) → 2(a, R₀) (unit test). The Fortran STRIDE form on the slayer_growthrate branch of the Fortran GPEC, which carries one power of ψ_t' in Λ instead of two, fails this by √2 and is off by ψ_t'^{-1/2} (2.3× at ε = 0.1, 3.2× at ε = 0.05). That branch also feeds D_I rather than D_R and writes log(q) as shear; none of that is used here.
• Independent metric: ⟨B²⟩, ⟨|∇ψ_N|²⟩ and V' recomputed from R(ψ,θ), Z(ψ,θ), F(ψ), psio alone agree with ResistGeometry to 1e-7–1e-4 on interior surfaces of the circular and DIII-D-like cases.
• Wiring: the radial label enters dgeo only through k_ref (dgeo/k_ref identical across :midplane, :flux, :volume to 1e-10; unit test).
• Tests: runtests_slayer_inputs.jl 70/70, runtests_tj_analytic.jl 20/20, runtests_slayer_params.jl 49/49.

Radial-label sensitivity (benchmarks)

Two benchmark scripts regenerate every figure: benchmarks/benchmark_toroidal_delta_crit.jl (convergence, scale invariance, independent metric, DIII-D Δ_crit per surface) and benchmarks/benchmark_delta_crit_radial_label.jl (threshold margin Δ'_rs/Δ_crit per radial label on the DIII-D-like SLAYER example with a fresh Riccati Δ', plus a TJ circular ε scan). Figures are archived in the maintainer's benchmark repo (maintainer access): figures + NOTES.md — lar_convergence.png, diiid_delta_crit_234.png, dgeo_diiid_delta_crit.png, diiid_threshold_margin_by_label.png, tj_epsilon_scan_by_label.png.

DIII-D-like SLAYER example, n = 1 (Δ' diagonal reproduces #403's table: 9.097, −6.196, −16.374):

	2/1	3/1	4/1
Δ_crit toroidal / rfitzp, :midplane (current default)	1.17	1.25	1.32
Δ_crit toroidal / rfitzp, :flux	0.96	0.94	0.91
Δ_crit toroidal / rfitzp, :volume	1.00	1.01	1.01
label spread of the margin Δ'_rs/Δ_crit, rfitzp	31 %	32 %	41 %
label spread of the margin Δ'_rs/Δ_crit, toroidal	7.7 %	10 %	15 %

Because dgeo ∝ K exactly and the converted Δ' ∝ K^{2μ} with μ ≈ 1/2, the toroidal threshold margin is nearly label-invariant; the cylindrical rfitzp formula (∝ r_s/√(da/dψ)) is not. The residual toroidal spread is ~1 % from K^{2μ−1} and ~6 % from χ∥ through the (cylindrical) W_d loop. With either branch the 2/1 sits at about half its threshold under every label, whereas with dc_type = "none" it is unstable.

TJ circular ε scan (fresh Δ' per point, surfaces 2/1 and 3/1): the label spread of the margin and the toroidal/rfitzp deviation both vanish roughly linearly in ε (rfitzp spread 1.3–1.7 % at ε = 0.05, toroidal 0.2 %; toroidal/rfitzp 1.006–1.012 midplane, 0.995–0.998 flux). No residual at small ε, so no implementation artefact; finite-ε differences are label physics.

Notes for reviewers

1. Reference-conversion exponent. The Y → x̂ conversion is linear, which is exact in the paper's H = 0 ordering. Using Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching #403's K^{2μ} exponent for the offset as well would multiply the toroidal Δ_crit by c^{2μ−1}: +5–12 % at the DIII-D 2/1–4/1 (c ≈ 2.2–2.6, μ ≈ 0.53–0.56), ×2 at the 5/1 (μ = 0.83, D_R = −0.67, outside the small-D_R ordering anyway). Documented in the toroidal_dgeo docstring; a reviewer decision.
2. Default radial label. This PR inherits :midplane from Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching #403. The evidence in InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417 and the tables above point to :flux (Fitzpatrick's toroidal-flux label, the one the layer formulas are derived in) as the default for shaped plasmas; :volume only coincides with rfitzp by accident of near-circularity. Flipping the default moves rs, S, shear and every SLAYER result and is deliberately left to InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417 as its own !-tagged PR.
3. D_geo is written for every dc_type once a ResistGeometry exists; it is a diagnostic and does not feed the dispersion unless dc_type = "toroidal".

---

NO MERGE WITHOUT HUMAN REVIEW. This PR is a draft and requires a named human reviewer and an assignee before it can be considered for merge. This is non-negotiable.

🤖 Generated with Claude Code