[logan-nc]

Release note

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

GPEC can now solve the steady modified Rutherford equation (MRE) per rational surface, with every term switchable: classical Δ′(w), finite-width saturation, bootstrap, curvature, ion polarization and the applied-field forcing of a coil-driven perturbed equilibrium. It reports seed thresholds, every saturated island width with the seed that reaches it, the resonant field equivalent to each seed island, and the critical applied fields at which a forced island opens and at which it triggers growth past the seed threshold. Callen's reduced threshold for coil-driven islands is a preset of the MRE.

What is added

• Tearing.ModifiedRutherford ([Tearing.ModifiedRutherford], solve(ModifiedRutherfordProblem(ffs; tearing, ...))).
  • Terms: F(w) on a log-spaced grid of full widths, built from the surfaces' Δ′, ResistGeometry (D_I, D_R, H) and the kinetic profiles.
  • Curvature: curvature_model = "lutjens" (default) matches SLAYER's toroidal critical offset D_c at small w; "benjamin" gives the Benjamin form.
  • Saturation: the finite-width term of Fitzpatrick 2023 Eq. 9.15 (off by default), with the current profile from Grad–Shafranov.
  • Gate: the SLAYER TearingProblem supplies γ, the layer width δ_r that starts the root search, and D_c.
  • Roots:
    • the seed threshold w_m and the first saturated width w_sat;
    • every saturated branch w_sat_all, with the seed that reaches it w_seed_all, and w_seed_unsaturated;
    • the resonant field of the vacuum island as wide as each seed (b_m, b_seed_all, b_seed_unsaturated, b_m_star, via the new public vacuum_island_field);
    • validity flags: unsaturated, multiple_saturated, outside_plasma, gate_disagrees, saturation_unconverged, among others.
  • Forcing: with pe (a coil-driven perturbed equilibrium), the term Δ_RMP (dρ/dψ_N)(w_v/w)² of Callen's Eq. 17 enters F, so the roots and branches include the applied field (forcing = false switches it off). The critical fields b_crit (a forced window opens below the unforced seed threshold) and b_trigger (the window reaches that threshold) follow from the unforced F, since the forcing scales with the vacuum resonant field.
  • Callen preset: preset = "callen" (CALLEN_PRESET) reduces F to Callen's steady state in his units. Its roots and b_crit equal the positive roots of Callen's cubic and his closed-form critical field (API test against the closed form evaluated outside the MRE, 10⁻⁸). The critical vacuum island width is w_v (b_crit/B_vac)^½.
  • Output: Tearing/ModifiedRutherford/{PerSurface,Terms} in gpec.h5, written from one annotation table.
• Tearing.ResonantForcing: the resonant forcing the MRE uses (vacuum and total resonant fields, Δ_RMP, vacuum island width on the toroidal-flux radius ρ) and forcing_thresholds. When the force-free-states solve has no Δ′ (the perturbed equilibrium needs a forward solve), main adds one Riccati solve for the MRE, timed as tearing_delta_prime.
• InnerLayer.layer_plasma: one public view of the kinetic plasma (profiles, Z_eff, ion mass and charge, resistivity model). The SLAYER tearing solve, the layer-overlap scan, and the MRE now use it in place of private helpers. Values are unchanged.
• Docs: the MRE section of docs/src/inner_layer.md, with the forcing and the Callen preset; its group in docs/development/hdf5-conventions.md.
• Harness cases diiid_slayer_n1_mre (a base = variant of diiid_slayer_n1_toroidal) and diiid_n1_callen (the MRE Callen preset on the coil forcing).

Part of #630 and #629. Still open after this PR:

• #630 (general MRE): a flux-surface-averaged bootstrap current ⟨j_bs·B⟩ (Sauter/Redl), so the default can use one normalization family for the bootstrap and curvature terms.
• #629 (Callen island features): pressure flattening across saturated islands.

Regression report

All 29 cases, origin/develop (a1149c5) against this branch at 862b2f6:

tj_n1 | 25 unchanged
diiid_slayer_n1_toroidal | 5 unchanged
gal_resistive_diiid_galinner | 11 unchanged
ggj_ray_q500i | 4 unchanged
solovev_kinetic_nuzero | 14 unchanged
diiid_error_field | 10 unchanged
diiid_slayer_n1 | 22 unchanged
diiid_truncation_n1 | 12 unchanged
diiid_slayer_n1_mre | 12 missing/N/A (new case)
tj_kinetic_calculated | 20 unchanged
diiid_slayer_n1_coupled | 7 unchanged
gal_resistive_diiid | 11 unchanged
solovev_n1 | 22 unchanged
ggj_reference | 4 unchanged
tj_kinetic_ntv | 7 unchanged
solovev_multi_n | 15 unchanged
diiid_n1 | 53 unchanged
tj_kinetic_ntv_rates | 2 unchanged
diiid_slayer_n1_chi_perp_e | 22 unchanged
diiid_n1_riccati | 17 unchanged
diiid_n1_callen | 3 missing/N/A (new case)
lar_ggj_tearing_n1 | 3 unchanged
diiid_multi_n | 19 unchanged
diiid_slayer_n1_winding | 22 unchanged
efit_fixedbdy_separatrix | 5 unchanged
tj_kinetic_multiion | 6 unchanged
diiid_kinetic_multi_n | 8 unchanged
diiid_n1_riccati_precompiled | 17 unchanged
lar_reversed_shear_n1 | 22 unchanged

The two new cases have no develop baseline, since the features are new on this branch. Against the previous head of this branch (02de80f), the three diiid_n1_callen preset quantities and all twelve diiid_slayer_n1_mre quantities are unchanged, so retiring CallenProblem moves no MRE values. The four Callen quantities the case used to track went with it.

Notes for reviewers

• The default F(w) is provisional. It pairs the Lütjens curvature with Benjamin's bootstrap term, which come from different normalization families. The follow-up adds ⟨j_bs·B⟩ so that both can take Fitzpatrick's 4I₂ coefficient (Fitzpatrick 2023 Eq. 14.196).
• Carrying Lütjens' term through M(w). Lütjens' form is large-aspect-ratio. Carrying it through the classical factor M(w) = (w/2)^(2α−1)√(−4D_I) is a convention that keeps the small-w match with SLAYER's gate at any D_I. The docstrings say so.
• Linearly unstable surfaces. SLAYER γ > 0 marks a surface linearly unstable. It then gets no seed threshold but keeps its saturated widths.
• Branch output. Branch arrays are surface × branch, padded with NaN to the most branches found on any surface. HDF5 dims are (branch, surface), like the term matrices.
• Critical fields. With the island locked in phase with the field, the forced F is positive where k B_vac > G(w) = −F₀(w)w². b_crit is the smallest local minimum of G below the unforced seed threshold, the bottom of the search range included. b_trigger is the largest G between there and the threshold. For Callen's F₀ the first is his critical field exactly; the unit test checks both against brute force on an analytic G with a window and a barrier.
• No separate Callen solver. An earlier commit on this branch had a closed-form CallenProblem with [Tearing.Callen]; the preset replaces it, and a deck with [Tearing.Callen] stops at parse time. The closed form survives in the API test as the preset's reference, and Callen's Fig. 8 widths are checked through the MRE root finder. Its extra outputs (ρ_i, w_pol, Δ_c, the root count) are not written. It is in the branch history if reviewers want it back.
• Seed field. vacuum_island_field is the exact inverse of the perturbed equilibrium's island width. The API test checks it against island_half_width and resonant_area_weighted_field at every surface.
• Test run. On the full suite earlier on this branch, runtests_equil.jl (TokaMaker re-form) failed identically on develop and every other file passed. At this head, the files the last commit touches pass: runtests_modified_rutherford.jl 78/78, runtests_result_struct.jl 241/241 and runtests_solve_api.jl 256/256.
• Left for follow-ups. Sharing the θ-metric kernel with ResistEval is left for a follow-up pull request.