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

[ebursch]

Release note

• Audience: users
• Numerical impact: none (harness @ c42d558)
• Migration: none

Allows users to calculate the critical resonant field required for tearing at a rational surface. The model is based on the simple torque balance model in Cole and Fitzpatrick PoP 2006. Resolves issue #371. Main user input not already present in a Tearing.jl run is supplying viscosity information. This can be done either as angular momentum diffusivity or a direct magnetic Prandtl number. Other improvements in progress: #416

Regression report

================================================================
Case: diiid_n1 — DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium
================================================================

Regression Report: diiid_n1
=================================================================================================
Ref 1: develop  @ c42d558e (2026-08-20)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 42ffd010 (pinned), 5 threads/5 BLAS
Ref 2: local  @ local (2026-08-20)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 42ffd010 (pinned), 5 threads/5 BLAS
-------------------------------------------------------------------------------------------------
Quantity                                      develop          local            Diff       Status
-------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.013943e-01     8.013943e-01     0.0e+00    OK    
total energy Im(et[1])                        1.233500e-04     1.233500e-04     0.0e+00    OK    
plasma energy Re(ep[1])                       -1.348114e+00    -1.348114e+00    0.0e+00    OK    
vacuum energy Re(ev[1])                       2.149508e+00     2.149508e+00     0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873975e-01     1.873975e-01     0.0e+00    OK    
plasma energy (all)                           [35 elem]        [35 elem]        0.0e+00    OK    
vacuum energy (all)                           [35 elem]        [35 elem]        0.0e+00    OK    
total energy (all)                            [35 elem]        [35 elem]        0.0e+00    OK    
ODE steps (saved)                             2655             2655             0.0e+00    OK    
ODE steps (total)                             4738             4738             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_t                                        1.327024e-02     1.327024e-02     0.0e+00    OK    
beta_n                                        1.372511e+00     1.372511e+00     0.0e+00    OK    
internal inductance li1                       8.842230e-01     8.842230e-01     0.0e+00    OK    
internal inductance li2                       7.080727e-01     7.080727e-01     0.0e+00    OK    
internal inductance li3                       7.304309e-01     7.304309e-01     0.0e+00    OK    
poloidal beta betap1                          6.680738e-01     6.680738e-01     0.0e+00    OK    
poloidal beta betap2                          5.349836e-01     5.349836e-01     0.0e+00    OK    
poloidal beta betap3                          5.518763e-01     5.518763e-01     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    
current beta betaj                            4.236479e-01     4.236479e-01     0.0e+00    OK    
plasma volume                                 1.829472e+01     1.829472e+01     0.0e+00    OK    
plasma current                                1.152130e+00     1.152130e+00     0.0e+00    OK    
mpert                                         35               35               0.0e+00    OK    
npert                                         1                1                0.0e+00    OK    
toroidal field bt0                            2.006573e+00     2.006573e+00     0.0e+00    OK    
wall field bwall                              3.880145e-01     3.880145e-01     0.0e+00    OK    
aspect ratio                                  2.845746e+00     2.845746e+00     0.0e+00    OK    
elongation kappa                              1.708322e+00     1.708322e+00     0.0e+00    OK    
q profile (checksum)                          ed7c21fd61df...  ed7c21fd61df...  identical  OK    
pressure profile (checksum)                   e15550827bf1...  e15550827bf1...  identical  OK    
Mercier D_I profile (checksum)                5a6fcb1c3a97...  5a6fcb1c3a97...  identical  OK    
resistive interchange D_R profile (checksum)  6284a4c9a75a...  6284a4c9a75a...  identical  OK    
ballooning Delta' profile (checksum)          44bf968c25d5...  44bf968c25d5...  identical  OK    
island half-widths                            [5 elem]         [5 elem]         0.0e+00    OK    
Chirikov parameter                            [5 elem]         [5 elem]         0.0e+00    OK    
||resonant area-weighted field||              5.207739e-04     5.207739e-04     0.0e+00    OK    
PE plasma energy                              3.422586e+00     3.422586e+00     0.0e+00    OK    
PE vacuum energy                              3.174509e+00     3.174509e+00     0.0e+00    OK    
PE surface energy                             5.841099e+00     5.841099e+00     0.0e+00    OK    
PE toroidal torque                            -5.062793e-02    -5.062793e-02    0.0e+00    OK    
NTV torque FGAR [N·m]                         5.587060e-01     5.587060e-01     0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                6.903492e-02     6.903492e-02     0.0e+00    OK    
Runtime (s)                                   242.4s           357.4s                      --    
resonant area-weighted field b^r              [5 elem]         [5 elem]         0.0e+00    OK    
=================================================================================================
Summary: 47 unchanged

Julia Torque Balance Output
<img width="688" height="377" alt="image" src="https://github.com/user-attachments/assets/7d1da4c2-2eb1-4d39-be61-3306359d029e" />
Similar Fortran Torque Balance Output
<img width="638" height="478" alt="image" src="https://github.com/user-attachments/assets/278ecfb8-c686-4823-9ffe-d144fc5dfe0d" />

---

[logan-nc]

Currently takes either an array or a scalar for all surfaces but will easily be updated to a true profile once I have angular momentum diffusivity profile data on hand to test with.

Sounds like this should be marked draft until this is done

This should perform as well as or better than the Fortran implementation.

"Should" is scary - benchmark it! There is a benchmark script for comparing to fortran - modify that for your needs (using scratch scripts - no need to commit run scripts)

Other improvements to the model and the numerical methods are in progress, but I wanted to have a working version available to users.

Mark this as draft until the IO and benchmarks are done then ping us for reviews. You can make a "Task" type issue describing future plans (and assign yourself) so we have it on record as work in progress.

[ebursch]

@logan-nc @d-burg Should be ready for review now. Added profile inputs and a Fortran comparison plot.

[d-burg]

Review (combining my findings with Claude's)

Structure is sound and Cole Eq. 62 is the right target. Physics claims below were checked numerically against a DIII-D-like 2/1 SLAYERParameters (n_e 5e19, T_e=T_i 1 keV, q 2, s 1, B_t 2 T, r_s 0.5, R_0 1.7, ω_*e −1e4, ω_*i 5e3, ω_E 3e4), not just read.

Checked, no action needed

• Q_e / Q_i signs check out. The scan axis is Q = +tauk·ω (Cole's convention), Q0_here matches it, and the odd-looking argmin(abs.(Qs .+ Q_e)) correctly locates Q = −Q_e (Im Δ changes sign at −0.860 vs predicted −0.8602).
• The positive-balance branch lands exactly on [−Q_e, Q0] = [−0.86, +2.58] — Cole's physical picture falling out on its own. Good self-consistency signal, and worth an assertion (see item 1).
• alpha = 1e-2 is not a free parameter: b_crit moves 4% across 1e-4 → 1e-1 with Qpeak unchanged. The max is interior, far from the branch edge. Defensible as-is; just state the measured insensitivity instead of leaving a bare literal.
• Grid resolution is benign for b_crit: n=200 → 4.5582e-3 vs n=2000 → 4.5593e-3.

Blocking

1. The conjugation in Riccati.jl is now significant but documented as cosmetic. solve_inner returns conj(Δ̂_s), justified as harmless because "conjugation maps zeros to zeros" — true for root-finding, which was the only consumer until now. Torque balance is the first consumer that depends on the sign of Im Δ. Undoing it: positive branch becomes the entire [−10, +10] window, Qpeak −2.246, b_crit 4.915e-3 — i.e. the max becomes an artifact of Qmin/Qmax. The code is correct as written, but nothing guards that. Please add a test pinning the Im Δ sign (or Qpeak ∈ (-Q_e, Q0)) and a note in Riccati.jl that the reflection is load-bearing for CriticalResonantField.

2. κ̂ is silently hardcoded. TorqueBalance.jl:131 does sqrt(maximum / lu * (sval^2/2)), i.e. 1/κ̂ → s²/2 (I think this is correct but double check for me), fixing Cole's viscosity integral ∫[μ(r_s)/μ(r)]dr/r at exactly 1/2. For flat μ that integral is ln(a/r_s): ≈0.51 at r_s/a≈0.6, but 1.20 at r_s/a=0.3 — κ̂ off 2.4×, b_crit off ~1.55× on inner surfaces. Also self-defeating: the viscosity profile is the headline new input, and κ̂ is the one place it's physically required.

3. Q0 is E×B only. run_slayer.jl:512 uses profiles(psi).omega, and _load_profiles sets omega_e/omega_i to zeros. Cole's ω₀ is the natural mode rotation frequency; his §III D specialises to ω_*=0, which isn't the regime here. Since the branch is [−Q_e, Q0], Q0 sets its upper edge and so resizes the domain of the max — this isn't just a shifted number.

4. Default viscous_input = 1.0 fabricates χ = 1.0 m²/s at every surface, silently, straight into b_crit ∝ √P. Remove the default and error if unset.

Should fix

• The "second local maximum" logic doesn't find local maxima (latent — never fired in my case). min_ΔQ = 0.001*|Qmax-Qmin| = 0.02, while grid spacing is 0.10 at the default n=200 (filter is a no-op, "second maximum" is the grid neighbour) and 0.010 at n=2000 (≥3 points away, still the same peak's flank). Broken at any resolution when it does fire. Sign change in the discrete derivative is a few lines. Related: pole detection uses exact index equality, and Q_e/Q_i are n=1 values so it degrades at n≥2.
• Scan window is user-fixed but the branch is physics-set. Only [−Q_e, Q0] contributes (~17% of samples here), so most of n is wasted and a surface with Q0 outside [−10,10] silently under-resolves. Derive the window from (−Q_e, Q0) with margin.
• File-header equation is wrong: it writes T_v = 2P(Q0−Q)/((S·κ̂)·(b_r/B_φ))² with S·κ̂ inside the square. Cole Eq. 61 is linear in S·κ̂, which is what the code implements.

Example is internally inconsistent

gpec.toml sets viscous_input_type = "magnetic_prandtl_number" with viscous_input = "test_ang_mom_diff" — a dataset named angular-momentum-diffusivity fed in as a Prandtl number. Values are 1.25–3.0: plausible as χ_φ (though it's a pointname inherited from ONETWO), implausible as P (the code's own P < 1 warning sits right next to it). Correct path gives P ≈ 33 at this case's η (7.5e-8 Ω·m) — about an order out, so ~3× in b_crit. The notebook confirms the round trip is broken: ang_mom_diff = P .* eta / (4π*1e-7) recovers ~0.1 m²/s, not the 1.25–3.0 fed in. If the validation figure in the PR body came from this deck, please re-run it before we treat it as validation.

Separately, test_ang_mom_diff was injected into TkMkr_D3Dlike_Hmode_kinetic.h5 — a test-prefixed dataset added to a file shared with the sibling ideal example, outside the documented kinetic schema.

Merge readiness

• PR is CONFLICTING with develop; branch is 37 commits behind.
• Regression report pinned at c42d558 is stale — re-run after conflicts are resolved.
• CI "Title and release note" failing (4 runs).
• No tests and no regression case for ~430 lines of new physics + workflow. Item 1 names the invariant most worth pinning.

Cleanup

• All four new docstrings are detached — blank line between the closing """ and the definition (TorqueBalance, torque_balance_scan, CriticalResonantFieldResult, run_critical_resonant_field). Verified: @doc returns nothing. Invisible to ? and Documenter. torque_balance_value has none at all, and CriticalResonantField isn't in any @autodocs block.
• store_scan is dead. Never read; scan_data is always populated and the writer gates only on !isempty, so the full Q/Δ/balance grid always lands in gpec.h5 (n=2000 in the example), contradicting both docstrings.
• ASCII profile files crash: the error text advertises "HDF5 or ASCII table" but the viscous loader does an unconditional HDF5.h5open.
• torque_balance_scan docstring mismatches code: documents an 8-tuple return (actual 6), a (model, params, Q0, P, lu, sval) signature (actual (tb;)), and n=200 (actual default 2000).
• Abstract field types: viscous_input::Any, params::AbstractVector{<:InnerLayerParameters}, scan_data::Vector{NamedTuple}; the last also stores the whole params struct per surface.
• combined_result = SLAYERResult(...) re-lists all 13 fields positionally; adding a field breaks it silently.
• μ₀ hardcoded as 4π*1e-7 at run_slayer.jl:528 and :535 despite PhysicalConstants.jl (touched by this PR).
• abs(chi_here) swallows negative/garbage input instead of erroring.
• viscous_input has two undocumented semantics: array must be length nsurfaces (unknowable in advance), string resolves to a per-ψ interpolant.
• Dead code in torque_balance_scan: positive recomputed as idx_maxima; maxima unused; isempty(selected) unreachable; maximum shadows Base; println should be @info.
• Unused imports in the CriticalResonantField module: LinearAlgebra, StaticArrays, GGJModel, GGJParameters, SLAYERModel, SLAYERParameters.
• critical_resonant_field_control_from_toml has a vestigial :inner_model branch for a nonexistent field; flat is a verbatim copy; unknown keys give an opaque kwarg MethodError.
• HDF5: Qpeak/Q0/P aren't snake_case, and no new dataset is in the annotation table — br_crit ships with no units="T".

Hygiene

• .gitignore adds src/Tearing/CriticalResonantField/CRF Dev/, a scratch dir inside the source tree.
• ~1500 lines of trailing-whitespace stripping across 7 fixed-format coil .dat files, plus formatter churn in ~90 unrelated files. First-touch normalisation is expected, but at this scale splitting it into its own commit would make the diff reviewable.
• Notebook: hardcoded Pkg.activate("/Users/bursche/Documents/GitHub/JPEC_BCRIT"); for i in 1:6 when the case has 5 surfaces (KeyError); using Plots/h5open without using HDF5; ad-hoc plotting when the last cell already uses Analysis.Equilibrium.plot_equilibrium_summary.
• Typos: "eqution", "Cole PopP 2006" (x2, should be PoP), "resoant", "sufficent", "Prandlt".

Summary

Core calculation holds up better than expected: conventions consistent, branch lands where it should, α is a 4% effect. Before sign-off I'd want the item-1 test, κ̂, the example's viscous_input_type mismatch plus a re-run figure, and conflicts resolved with a fresh harness run. Cleanup can ride along.

_{Drafted with assistance from Claude Code; physics verified numerically against origin/feature/b_crit_calculation.}

[logan-nc]

Chiming in without having read all of the more thorough review 😝

• ensure this has regression tests
• ensure it is thoroughly documented in docs and renders correctly in the web (https://openfusiontoolkit.github.io/GPEC/previews/PR415/)
• The biggest usability hurdle here is the need for a chi_phi profile (thus an external transport code run) → please allow a simple override for a user to specify a simple 0D tau_E ("intuitively" approximated) instead (P from τ_R/τ_E, with χ_φ = κ_E a²/τ_E and κ_E the profile-shape factor of Fitzpatrick Phys. Plasmas 2023, Eqs. 38 and 55)
• Argument for not using the less approximated equation was to keep it the same as fortran, so some sort of fortran comparison seems warranted. This PR needs PLOTS!
• Would love to see smart pole treatment per Tearing.CriticalResonantField - Improve pole avoidance, input options, and grid refinement #416 but I guess this could merge with that still open as well as it is well documented