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

[ebursch]

Release note

• Audience: users
• Numerical impact: none; the critical resonant field is off by default (harness @ acd8bca)
• Migration: none

Tearing can now compute the critical resonant field b_r/B_φ for error-field penetration at each SLAYER surface, using the torque balance of Cole and Fitzpatrick, Phys. Plasmas 13, 032503 (2006), Eqs. 61–62. Turn it on with [SLAYER.CriticalResonantField] enabled = true. The magnetic Prandtl number is each surface's SLAYER P_tor, so χ_φ from the kinetic file (or a future τ_E closure) carries through with no extra input. Results are written to Tearing/CriticalResonantField/ in gpec.h5. Resolves #371.

What changed since the reviews

• P is inherited from InnerLayer (IL). The viscous_input / viscous_input_type options, the second HDF5 read of the kinetic file, the default χ = 1 m²/s, and the test_ang_mom_diff dataset are all gone. P is each surface's P_tor = μ0 χ_φ/η, the same quantity the old code computed, now taken straight from SLAYER.
• Cole's Q axis. solve_inner mirrors the real-Q axis (conj Δ_Julia(−Q) = Δ_Fortran(Q)). The PR was pairing that mirrored Δ with an unmirrored Q0. The torque balance now evaluates Δ(Q) = conj(Δ_solve_inner(−Q)), so the diamagnetic poles sit at +Q_e and +Q_i and the solution branch lies between Q_e and Q0, as in Cole and in Fortran gslayer.f. Growth-rate root finding is untouched.
• Q0 = τ_k·n·ω_E. This follows Cole's ω0, the mode frequency in the E×B frame; the diamagnetic drifts enter separately through Q_e and Q_i. The factor n now matches Q_e and Q_i. ω_E is evaluated at each surface's own ψ.
• κ̂ kept at Fortran parity. The viscosity integral stays fixed at 1/2 (κ̂ = 2/s²) for 2.0, and the assumption is documented. Profile-based κ̂ is tracked in Tearing - Compute κ̂ from the viscosity profile in the critical resonant field #510 for after the release.
• Maximum selection.
  • The scan window follows gslayer.f: between Q0 and the Q_e pole. Qmin/Qmax are now optional overrides.
  • The code takes the largest positive interior local maximum, skipping any within 0.02 of Q_e or Q_i (Fortran SLAYER's pole-regularization radius, layfac).
  • This replaces the "two largest samples" logic and the exact-index pole test.
• Cleanup.
  • Docstrings attached and corrected, header equation fixed, @info instead of println, unused imports removed.
  • Concrete result types; store_scan is honoured; combined_result is rebuilt from fieldnames.
  • snake_case HDF5 datasets with metadata (br_crit, q_peak, q0, p_phi, rational_index).
  • The notebook, the Cole PDF (cited by DOI instead), the .gitignore entry and the renamed example directory are removed. The CRF block now lives in examples/DIIID-like_SLAYER_example/gpec.toml.
  • About 90 files of formatter churn were reverted to develop's version in the develop merge.
• Tests and docs.
  • test/runtests_critical_resonant_field.jl (28 tests) covers:
    • the Cole-axis branch, Q_e < Q_peak < Q0;
    • the window rule;
    • local-maximum and pole selection;
    • P = P_tor and Q0 = τ_k·n·ω_E;
    • TOML parsing;
    • the HDF5 metadata contract.
  • diiid_slayer_n1 now pins b_r/B_φ crit and Q_peak on the inner three surfaces.
  • A "Critical resonant field" section with the benchmark figure was added to the Tearing docs page.

Fortran benchmark

Both codes were run on Fortran's normalized inputs for the DIII-D-like example: 2/1–6/1, n = 1, P = 1, κ̂ = 2/s², Q ∈ [−5, 5], 20001 points.

• Fortran side. gpec with singthresh_slayer_flag, from bcrit_fine_grid @ 96d0cfbe (ebursch/GPEC). That commit fixes maxbal there: b_crit was always 0, and MAXLOC was off by one.
• Julia side. Shared-level SLAYERParameters built from Fortran's S, τ_k, Q_e, Q_i, τ and c_β.

m/n	Q0	Q_e	Q_i	Fortran gslayer.f (global max)	Julia routine on Fortran's Δ table	Julia layer	median rel. Δ diff, Q_i < Q < Q_e
2/1	5.448	1.383	−2.147	3.049e-4 @ Q = 1.381	2.516e-4 @ Q = 2.445	2.235e-4 @ Q = 2.511	36 %
3/1	2.855	0.898	−0.738	1.948e-4 @ Q = 0.901	1.827e-4 @ Q = 1.721	1.499e-4 @ Q = 1.722	78 %
4/1	0.570	0.757	−1.874	2.347e-4 @ Q = 0.755	2.833e-5 @ Q = 0.671	2.510e-5 @ Q = 0.667	24 %
5/1	−1.289	4.798	−3.319	1.057e-1 @ Q = −0.370	1.057e-1 @ Q = −0.370	1.042e-1 @ Q = −0.575	7 %
6/1	−0.092	0.253	−0.430	7.919e-4 @ Q = 0.252	4.340e-4 @ Q = 0.100	3.785e-4 @ Q = 0.082	62 %

• The torque-balance port is exact. Fed Fortran's own Δ 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 all five surfaces.
• Fortran's global maximum is the Q_e pole spike. On 2/1, 3/1, 4/1 and 6/1 it sits within 0.003 of Q_e. Skipping the pole region selects the torque-balance maximum between Q_e and Q0 instead, which lowers b_crit by 6–88 % on the same Δ table (third column).
• The layer models differ. 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. Julia's b_crit is 82–89 % of the Fortran-table value on those surfaces.
• 5/1 is not a usable threshold. That edge surface is dominated by a zero of the electromagnetic torque (bal ≈ 2.7e5) in both codes. These are the same surfaces the SLAYER example already flags as unreliable.

The figure is also in the docs preview under Tearing → Critical resonant field. Run scripts are not committed.

Regression report

regress --cases diiid_slayer_n1,diiid_n1 --refs 4b3a7014a,acd8bca78 (develop vs this branch, both fresh, manifest pinned).

Regression Report: diiid_slayer_n1
==========================================================================================
Ref 1: 4b3a7014a  @ 4b3a7014a (2026-10-06)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 078abfaa (pinned), 112 threads/1 BLAS
Ref 2: acd8bca78  @ acd8bca78 (2026-10-06)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 078abfaa (pinned), 112 threads/1 BLAS
------------------------------------------------------------------------------------------
Quantity                            4b3a7014a  acd8bca78  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_e                          [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER Q_i                          [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   1.5e-05            OK           
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   6.1e-02            OK           
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   2.331e-01 (0.02%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00            OK           
CRF b_r/B_φ crit [2/1,3/1,4/1]      N/A        [3 elem]   N/A                N/A          
CRF Q_peak [2/1,3/1,4/1]            N/A        [3 elem]   N/A                N/A          
SLAYER enabled flag                 1          1          0.0e+00            OK           
Runtime (s)                         269.5s     295.4s                        --           
==========================================================================================
Summary: 1 changed, 16 unchanged, 2 missing/N/A


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

Regression Report: diiid_n1
=================================================================================================
Ref 1: 4b3a7014a  @ 4b3a7014a (2026-10-06)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 078abfaa (pinned), 112 threads/1 BLAS
Ref 2: acd8bca78  @ acd8bca78 (2026-10-06)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 078abfaa (pinned), 112 threads/1 BLAS
-------------------------------------------------------------------------------------------------
Quantity                                      4b3a7014a        acd8bca78        Diff       Status
-------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.012318e-01     8.012318e-01     0.0e+00    OK    
total energy Im(et[1])                        4.142530e-05     4.142530e-05     0.0e+00    OK    
plasma energy Re(ep[1])                       -1.348486e+00    -1.348486e+00    0.0e+00    OK    
vacuum energy Re(ev[1])                       2.149718e+00     2.149718e+00     0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873976e-01     1.873976e-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)                             1429             1429             0.0e+00    OK    
ODE steps (total)                             1424             1424             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.842392e-01     8.842392e-01     0.0e+00    OK    
internal inductance li2                       7.080847e-01     7.080847e-01     0.0e+00    OK    
internal inductance li3                       7.304433e-01     7.304433e-01     0.0e+00    OK    
poloidal beta betap1                          6.680744e-01     6.680744e-01     0.0e+00    OK    
poloidal beta betap2                          5.349834e-01     5.349834e-01     0.0e+00    OK    
poloidal beta betap3                          5.518761e-01     5.518761e-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.236478e-01     4.236478e-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.708350e+00     1.708350e+00     0.0e+00    OK    
q profile (checksum)                          0cd285cea88d...  0cd285cea88d...  identical  OK    
pressure profile (checksum)                   a1c48b266622...  a1c48b266622...  identical  OK    
Mercier D_I profile (checksum)                eeb06744e795...  eeb06744e795...  identical  OK    
resistive interchange D_R profile (checksum)  fa37296851f6...  fa37296851f6...  identical  OK    
ballooning Delta' profile (checksum)          6eb0ea075ccf...  6eb0ea075ccf...  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.189042e-04     5.189042e-04     0.0e+00    OK    
dominant-coupling singular values             [3 elem]         [3 elem]         0.0e+00    OK    
|forcing overlap with dominant mode|          1.415706e-04     1.415706e-04     0.0e+00    OK    
|delta_nominal| of coil set 1                 7.055344e-05     7.055344e-05     0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.246238e-04     1.246238e-04     0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.180652e-07     4.180652e-07     0.0e+00    OK    
PE plasma energy                              3.422677e+00     3.422677e+00     0.0e+00    OK    
PE vacuum energy                              3.174510e+00     3.174510e+00     0.0e+00    OK    
PE surface energy                             5.826684e+00     5.826684e+00     0.0e+00    OK    
PE toroidal torque                            5.087465e-02     5.087465e-02     0.0e+00    OK    
NTV torque FGAR [N·m]                         5.296746e-01     5.296746e-01     0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                7.924988e-02     7.924988e-02     0.0e+00    OK    
Runtime (s)                                   311.4s           302.2s                      --    
||forcing b~|| (root-area-weighted)           4.683328e-04     4.683328e-04     0.0e+00    OK    
resonant area-weighted field b^r              [5 elem]         [5 elem]         0.0e+00    OK    
=================================================================================================
Summary: 53 unchanged

diiid_n1 is unchanged. In diiid_slayer_n1 every layer parameter is unchanged and the two new CRF pins are recorded. SLAYER γ on the inner three surfaces differs by 0.23 Hz (0.02 %), just over the case's 0.1 Hz threshold, with Q_root and ω inside theirs. Nothing on the root-finding path changed in this PR: the only edit there is a comment in Riccati.jl, and the CRF runs after the roots are found. This is the run-to-run spread of the threaded AMR root search tracked in #420. A single-threaded rerun (GPEC_REGRESS_THREADS=1 regress --cases diiid_slayer_n1 --refs 4b3a7014a,acd8bca78 --force) confirms it: all 17 existing quantities are unchanged, with γ differing by 1.3e-5 Hz, Q_root by 5.1e-10 and ω by 5.9e-7 Hz.

Notes for reviewers

• Depends on nothing; complements InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501. InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501's P_perp_model = "tau_E" sets P_tor = τ_R/τ_E, so once it merges the 0D τ_E route logan-nc asked for reaches the critical resonant field with no further change.
• Mirror note. A one-line comment in src/InnerLayer/SLAYER/Riccati.jl records that real-Q callers see a mirrored axis. The sign of SLAYER's root frequencies is a separate question, raised in the Julia-vs-Fortran SLAYER comparison, and is left alone here.

---

🤖 Generated with Claude Code

[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

[ebursch]

Just fyi, working on a clean-up and mostly Fortran parallel implementation of this for GPEC 2.0. One substantial change vs the above work and conversation is that now that I properly understand the diffusion inputs and metrics for the inner layer, the whole critical resonant field specific momentum input stuff is not necessary. Going forward it will just inherit the P_phi calculated for the inner layer, whether that is from chi_phi or tau_e.

[ebursch]

Modest pole avoidance will be implemented here but a more robust, database-ready version along with some other nice to have features like profile-dependent kappa, and the explicit Delta prime inlcusion will come post v 2.0.

[ebursch]

@d-burg thanks for the thorough review. Point by point:

Blocking

1. conj in solve_inner. This turned out to be more than a sign to pin. solve_inner mirrors the real-Q axis (conj Δ_Julia(−Q) = Δ_Fortran(Q)), and the torque balance was pairing that mirrored Δ with an unmirrored Q0. The branch you found, [−Q_e, Q0], was mixing the two conventions.
   • The torque balance now uses Δ(Q) = conj(Δ_solve_inner(−Q)) (cole_delta). The poles sit at +Q_e and +Q_i, and the branch lies between Q_e and Q0, as in Cole and in gslayer.f (4f713fc).
   • The test pins Q_e < Q_peak < Q0, and a one-line note in Riccati.jl records the mirror.
   • Root finding is untouched.
2. κ̂. You're right that sval^2/2 fixes Cole's ∫μ(r_s)/μ(r) dr/r at 1/2. We're keeping it there for 2.0 so the code matches Fortran gslayer.f exactly. The assumption is now stated in the header and the docs, and κ̂ is an explicit kappa_hat field. The profile integral, with your ln(a/r_s) numbers, is tracked in Tearing - Compute κ̂ from the viscosity profile in the critical resonant field #510 for after the release.
3. Q0. Cole §IV defines ω0 as "the (angular) oscillation frequency of an m,n mode comoving with the E×B frame … in the absence of an error field". The diamagnetic drifts enter separately through Q_e and Q_i in the layer, so adding ω_* to Q0 would count them twice. Fortran gslayer.f also uses Q0 = τ_k·ω_E. Two real problems were there, though, and both are fixed (baf49d7):
   • Q0 lacked the factor n that Q_e and Q_i carry, so it is now τ_k·n·ω_E.
   • ω_E was read through intr[isurf]; it is now read at each surface's own ψ.
4. Default χ = 1. Removed along with the whole viscous input. P is now each surface's P_tor from SLAYER, so χ_φ from the kinetic file (or the scalar chi_tor fallback, which already warns) sets it (baf49d7).

Should fix

• Local maxima. These now come from sign changes in the discrete slope. A maximum within 0.02 of Q_e or Q_i is skipped; 0.02 is Fortran SLAYER's layfac pole-regularization radius. The Fortran benchmark showed why a one-step exclusion is not enough: on a 20001-point grid, Fortran's own global maximum sits on the Q_e spike on four of five surfaces (f032554). The pole test uses the n-correct Q_e and Q_i from develop.
• Scan window. This is now derived per surface from Q0, Q_e and Q_i using the gslayer.f rule, with Qmin/Qmax as optional overrides.
• Header equation. Fixed: linear in S·κ̂, with κ̂ defined.

Example

• The viscous_input mismatch is gone with the input itself, and test_ang_mom_diff is gone from the kinetic h5.
• The figure in the PR body has been re-run from the new code.

Merge readiness

• Merged develop, and reverted all formatter churn to develop's version, so the diff is now only the CRF files.
• Fresh harness run; report in the PR body.
• The release-note block is fixed.
• 28 unit tests added, plus b_crit and Q_peak pins in diiid_slayer_n1.

Cleanup and hygiene

All done:

• Docstrings attached and corrected, and torque_balance_value now has one.
• store_scan is honoured.
• The ASCII-profile crash, abs(chi), the hard-coded μ0 and the dual-semantics viscous_input all went with the viscous input.
• Concrete types, and combined_result is rebuilt from fieldnames.
• Unused imports, the :inner_model branch and the flat copy are removed.
• HDF5 datasets are snake_case with long_name and units.
• The .gitignore line, the notebook and the PDF are removed.
• Typos fixed.

🤖 Generated with Claude Code

[ebursch]

@logan-nc point by point:

• Regression tests. diiid_slayer_n1 now enables the critical resonant field and pins b_r/B_φ crit and Q_peak on the inner three surfaces. A new test/runtests_critical_resonant_field.jl has 28 unit tests. The harness report is in the PR body.

• Docs. There is a "Critical resonant field" section on the Tearing page with:

  • the equations, inputs and approximations;
  • the HDF5 outputs;
  • the benchmark figure.

  The module is in the API reference. Preview: https://openfusiontoolkit.github.io/GPEC/previews/PR415/

• 0D τ_E instead of a χ_φ profile. The critical resonant field no longer takes its own viscosity input; it uses each surface's SLAYER P_tor. InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501 adds P_perp_model = "tau_E", which sets P_tor = τ_R/τ_E (Fitzpatrick 2023, Eq. 131). Once InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501 merges, that override reaches b_crit with no extra code here.

• Fortran comparison with plots. These are in the PR body and the docs:

  • Fed Fortran's own Δ(Q), the Julia torque balance reproduces gslayer.f's b_crit to all printed digits.
  • Fortran's global maximum turns out to sit on the Q_e pole spike on four of five surfaces. Both codes now skip the pole region (0.02, Fortran's layfac).
  • The remaining 11–18 % b_crit gap on the inner surfaces is the Park vs Fitzpatrick layer model.
  • On the Fortran side, b_crit had been stuck at 0 on bcrit_fine_grid; that is fixed in ebursch/GPEC@96d0cfbe.

• Smart poles (Tearing.CriticalResonantField - Improve pole avoidance, input options, and grid refinement #416). This is still open, and this PR doesn't depend on it. The 0.02 pole exclusion (from Fortran) is the stopgap.

🤖 Generated with Claude Code