Tearing.ModifiedRutherford - FEATURE! - Add a configurable steady modified Rutherford equation with seed thresholds, saturated widths and coil-forcing critical fields

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}` (+ `*Winding` with `root_classification = winding`), `Scan/Surface_<k>/`, `CriticalResonantField/` (+ `Scan/Surface_<k>/`) |
+| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}` (+ `*Winding` with `root_classification = winding`), `Scan/Surface_<k>/`, `CriticalResonantField/` (+ `Scan/Surface_<k>/`), `ModifiedRutherford/` (`w`, `w_m_star`, `b_m_star`, `b_crit_star`, `PerSurface/` with branch × surface `w_sat_all`-type arrays, `Terms/` surface × width) |
 | `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`.
@@ -62,7 +62,7 @@ Every dataset outside `Input/` (raw snapshot) must answer "what is this, in what
 
 Root-level file attributes: `schema_version` (currently `"2.0"`; bump on breaking schema changes — readers dispatch on it), `Conventions = "GPEC-HDF5-2.0"`, `references`, `title` (run description), `date_created` (ISO 8601 UTC). The code version stays in `Info/git_version`.
 
-Mechanism: writers stay table-driven — each writer keeps a `path => (; long_name, units, dims, scale, attach)` table next to it (`scale` marks a coordinate dataset as a dimension scale; `attach` binds axes to scales — see the `Utilities.HDF5Annotations.annotate!` docstring) (`src/Utilities/HDF5Schema.jl` for the main writer; alongside `write_galerkin!`, the PerturbedEquilibrium writer, `KineticForces/Output.jl`, and `Tearing/HDF5Output.jl` for the rest) and applies it post-write via `Utilities.HDF5Annotations.annotate!`. Entries for conditionally-written datasets are simply skipped when absent. When adding a dataset, add its table entry in the same commit — the schema test fails otherwise.
+Mechanism: writers stay table-driven — each writer keeps a `path => (; long_name, units, dims, scale, attach)` table next to it (`scale` marks a coordinate dataset as a dimension scale; `attach` binds axes to scales — see the `Utilities.HDF5Annotations.annotate!` docstring) (`src/Utilities/HDF5Schema.jl` for the main writer; alongside `write_galerkin!`, the PerturbedEquilibrium writer, `KineticForces/Output.jl`, `Tearing/HDF5Output.jl` and `Tearing/ModifiedRutherford/HDF5Output.jl` for the rest) and applies it post-write via `Utilities.HDF5Annotations.annotate!`. Entries for conditionally-written datasets are simply skipped when absent. When adding a dataset, add its table entry in the same commit — the schema test fails otherwise.
 
 ## File-wide conventions
 

docs/src/inner_layer.md

@@ -103,6 +103,44 @@ The DIII-D-like kinetic file of the examples carries analytic χ profiles of the
 | `Tearing/Diagnostics/ValidRoots`, `Poles`, `FilteredRoots` | the root classification of the scan |
 | `Tearing/Scan/Surface_<k>/` | the scanned Δ(Q) with `store_scan = true` |
 | `Tearing/CriticalResonantField/br_crit`, `br_crit_tesla`, `Q_peak`, `Q0` | the penetration threshold as ``b_r/B_\phi`` and in tesla, the Q of the balance maximum and the natural rotation |
+| `Tearing/ModifiedRutherford/w`, `Terms/F`, `Terms/<term>`, `Terms/w_d`, `PerSurface/` | the modified Rutherford width grid, ``F(w)`` and each switched term per surface, the transport width, and the surface inputs |
+| `Tearing/ModifiedRutherford/w_m_star`, `PerSurface/w_m`, `PerSurface/w_sat`, `PerSurface/<flag>` | the seed threshold and saturated widths per surface, their minimum over the surfaces, and the validity and gate flags (0/1) |
+| `PerSurface/w_sat_all`, `PerSurface/w_seed_all`, `PerSurface/w_seed_unsaturated`, `PerSurface/n_sat` | every saturated width per surface (branch × surface, NaN-padded) with the seed threshold that reaches it, the seed above which the island grows past `w_max`, and the number of saturated widths |
+| `PerSurface/b_m`, `PerSurface/b_seed_all`, `PerSurface/b_seed_unsaturated`, `b_m_star` | the resonant area-weighted field [T] whose vacuum island is as wide as each seed threshold, ``B = \pi |n q'| \chi_1 w^2/(8A)`` (the inverse of the perturbed equilibrium's island width), to compare with `PerturbedEquilibrium` resonant fields and error-field thresholds |
+| `PerSurface/b_vac`, `b_res`, `delta_rmp`, `w_v`, `Terms/forcing` | with a perturbed equilibrium, the vacuum and total resonant fields, Callen's ``\Delta'_\mathrm{RMP}``, the vacuum island width and the forcing term per width |
+| `PerSurface/b_crit`, `b_trigger`, `w_crit`, `b_crit_star` | the vacuum resonant fields [T] at which a forced window opens below the unforced seed threshold (at width `w_crit`) and at which it reaches that threshold, and the smallest `b_crit` over the surfaces |
+
+## Modified Rutherford diagnostics
+
+`[Tearing.ModifiedRutherford]` evaluates the steady modified Rutherford equation on a grid of full island widths ``w`` (in ``\psi_N``) for each rational surface of a single-``n`` Riccati or Galerkin run with Δ′:
+
+```math
+F(w) = \Delta'(w) + \Delta_\mathrm{nc}(w) + \Delta_\mathrm{GGJ}(w) + \Delta_\mathrm{pol}(w; \omega) + \Delta_\mathrm{RMP}(w),
+```
+
+the classical drive with its finite-pressure power law, the bootstrap drive with the transport threshold ``w_d``, the Glasser–Greene–Johnson curvature term and the ion-polarization term (Benjamin et al. 2026, Eqs. 7, 10–14 and C.1). Each term has a switch and a coefficient; the polarization term is off by default. A finite-width saturation term, `saturation = true` (off by default), replaces ``r_s\Delta'`` by ``r_s\Delta'(0) - C_\mathrm{sat}W/r_s`` with ``C_\mathrm{sat} = 0.8\alpha_s^2 - 0.27\beta_s - 0.09\alpha_s`` from the current-density gradients at the surface (Fitzpatrick 2023, Eq. 9.15), so that a surface with ``\Delta' > 0`` saturates. The curvature term defaults to the Lütjens et al. (2001) form ``r_s\Delta_\mathrm{GGJ} = 6.35 D_R/\sqrt{(W/r_s)^2 + 0.65 W_d^2}`` with ``W_d`` from SLAYER's toroidal critical offset, so that at small widths it equals the critical ``\Delta`` of the linear SLAYER layer; `curvature_model = "benjamin"` selects Benjamin et al. (2026) Eq. 12. The surface geometry follows the Fortran RDCON `mercier.f` averages, and the densities, temperatures, ``\omega_E``, ``Z_\mathrm{eff}`` and ``\chi_\perp`` come from the equilibrium's kinetic profiles. A `[Tearing]` section holding only this sub-table runs no layer solve. In a script, `solve(ModifiedRutherfordProblem(ffs; polarization=true))` returns the same [`ModifiedRutherfordResult`](@ref Tearing.ModifiedRutherford.ModifiedRutherfordResult).
+
+The roots of ``F`` are bracketed on the grid and refined by Brent's method: the seed threshold ``w_m`` is the smallest crossing from ``F<0`` to ``F>0`` and the saturated width ``w_\mathrm{sat}`` the first crossing from ``F>0`` to ``F<0`` above where the island starts, ``w_\mathrm{lo}`` when ``F(w_\mathrm{lo})>0``, else ``w_m``.
+Flags mark an island that does not saturate in the search range (`unsaturated`), one whose ``w_\mathrm{sat}`` reaches past the equilibrium grid (`outside_plasma`), and a sign of ``F(w_\mathrm{lo})`` that disagrees with the gate's growth rate (`gate_disagrees`).
+When ``F`` crosses zero more than once, every saturated width is kept: an island seeded between `w_seed_all[k]` and the next seed threshold settles at `w_sat_all[k]`, and `multiple_saturated` marks surfaces where ``w_\mathrm{sat}`` is the first of several.
+When the deck also has a `[Tearing]` layer solve, its uncoupled growth rates gate the roots (a surface with ``\gamma>0`` grows from noise, so it has no seed threshold ``w_m``, but keeps its saturated width) and its resistive layer width, converted to ``\psi_N`` with the layer's own radial label, sets the lower end of the search; in a script, pass it as `ModifiedRutherfordProblem(ffs; tearing)`.
+Surfaces with ``D_R>0`` are excluded as resistive-interchange unstable, and roots below the ion banana width or with ``|\Delta'|w\ge1`` are flagged.
+Lengths in metres (the layer width, the poloidal ion Larmor radius and the banana width) map to ``\psi_N`` with the toroidal-flux label of `radial_label`, as SLAYER's do.
+The island frequency in ``\Delta_\mathrm{pol}`` is ``\omega_{*e}(1+\eta_e/4)`` (Wilson 1996) by default, or ``-n\omega_E``, an island at rest in the laboratory, with `omega_model = "locked"`.
+
+### Applied-field forcing and critical fields
+
+With a coil-driven `[PerturbedEquilibrium]` in the same deck (or `ModifiedRutherfordProblem(ffs; pe)` in a script), the applied field forces the island through Callen's term ``\Delta_\mathrm{RMP}(w) = \Delta'_\mathrm{RMP}(d\rho/d\psi_N)(w_v/w)^2`` (Callen 2016, Eq. 17), with ``w_v`` the vacuum island width and the island locked in phase with the field; `forcing = false` switches it off. The roots and branches are then those of the forced ``F``.
+The perturbed equilibrium needs a forward solve, so `main` adds one Riccati solve for the MRE's ``\Delta'``; the forward solve then needs `vac_flag = true`.
+Since the forcing is proportional to the vacuum resonant field ``B_\mathrm{vac}``, the forced ``F`` is positive where ``k B_\mathrm{vac} > G(w) = -F_0(w)w^2``, with ``F_0`` the unforced sum and ``k = \Delta'_\mathrm{RMP}(d\rho/d\psi_N)w_v^2/B_\mathrm{vac}``.
+`b_crit` is the field at which a forced window first opens below the unforced seed threshold (the smallest local minimum of ``G``), and `b_trigger` the field at which that window reaches the unforced seed threshold, so that a forced island grows on to the intrinsic saturated width (the largest ``G`` between the window and that threshold).
+Both are 0 on a surface that is unstable at ``w_\mathrm{lo}`` without forcing.
+Callen's ``-w_\mathrm{pol}^2/w^3`` beats the forcing at small widths and so gives a window away from ``w_\mathrm{lo}``; without such a term ``G`` rises from ``w_\mathrm{lo}``, the window opens there (`w_crit` equals `w_lo`), and `b_crit` depends on ``w_\mathrm{lo}``, which the `[Tearing]` gate sets to the resistive layer width. `b_trigger` does not.
+`preset = "callen"` (`CALLEN_PRESET` in a script) reduces ``F`` to Callen's reduced equation for RMP-driven islands, the steady state ``\Delta_c + \Delta'_\mathrm{RMP}(w_v/w)^2 - w_\mathrm{pol}^2/w^3`` in his units and in the form of Fortran GPEC (`delta_prime_model = "callen"`, `pol_model = "callen"`, the other terms off).
+Its roots are then the two positive roots of Callen's cubic ``\Delta_c w^3 + \Delta'_\mathrm{RMP}w_v^2w - w_\mathrm{pol}^2 = 0`` (Callen 2016, Eq. 17) in the full width ``w``, and its `b_crit` is his critical vacuum resonant field ``B_\mathrm{crit} = B_\mathrm{vac}(w_{v,\mathrm{crit}}/w_v)^2``; the critical vacuum island width is ``w_v\sqrt{b_\mathrm{crit}/B_\mathrm{vac}}``.
+The vacuum resonant field is the coils' resonant harmonic on the rational surface over its area.
+The intrinsic index ``\Delta_c`` is the computed ``\Delta'_{mn}`` referred to ``\rho`` by default, or Callen's cylindrical ``-2m/\rho`` of Fortran GPEC 1.6.0 with `delta_prime_source = "cylindrical"`.
+The same switch without the preset replaces ``\Delta'_{mn}`` by the cylindrical ``r_s\Delta' = -2m``, so the MRE then runs on any solve, forward included, as a reference case.
 
 ## Numerical notes
 
@@ -175,3 +213,15 @@ Modules = [GeneralizedPerturbedEquilibrium.Tearing.Dispersion]
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField]
 ```
+
+### ModifiedRutherford
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing.ModifiedRutherford]
+```
+
+### ResonantForcing
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing.ResonantForcing]
+```

regression-harness/cases/diiid_n1_callen.toml

@@ -0,0 +1,50 @@
+# Regression case: DIII-D-like n=1 Callen vacuum-island threshold and cubic roots, from the modified Rutherford
+# equation's Callen preset on the coil-driven perturbed equilibrium.
+# The cylindrical intrinsic Δ′ (−2m/ρ) keeps the Fortran GPEC parity case and needs no extra Δ′ solve.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "diiid_n1_callen"
+description = "DIII-D-like, n=1, Callen vacuum-island threshold and cubic roots from the MRE Callen preset on the coil-driven perturbed equilibrium"
+example_dir = "examples/DIIID-like_ideal_example"
+remove_sections = ["LocalStability", "Ballooning", "ErrorFields", "KineticForces"]
+
+[overrides]
+"Tearing.ModifiedRutherford.preset" = "callen"   # Reduce F(w) to Callen's steady state with the coil forcing
+"Tearing.ModifiedRutherford.delta_prime_source" = "cylindrical"   # Intrinsic Delta' -2m/rho, as Fortran GPEC 1.6.0; no extra Delta' solve
+
+[quantities.mre_callen_w_m]
+h5path = "Tearing/ModifiedRutherford/PerSurface/w_m"
+type = "real_vector"
+extract = "all_real"
+label = "MRE Callen preset: forced seed threshold full width (psi_N)"
+noise_threshold = 1e-10
+order = 18
+class = "physics_converged"
+
+[quantities.mre_callen_w_sat]
+h5path = "Tearing/ModifiedRutherford/PerSurface/w_sat"
+type = "real_vector"
+extract = "all_real"
+label = "MRE Callen preset: forced saturated full width (psi_N)"
+noise_threshold = 1e-10
+order = 20
+class = "physics_converged"
+
+[quantities.mre_callen_b_crit]
+h5path = "Tearing/ModifiedRutherford/PerSurface/b_crit"
+type = "real_vector"
+extract = "all_real"
+label = "MRE Callen preset: critical vacuum resonant field (T)"
+noise_threshold = 1e-12
+order = 22
+class = "physics_converged"
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
+class = "diagnostic"

regression-harness/cases/diiid_slayer_n1_mre.toml

@@ -0,0 +1,88 @@
+# Regression case: DIII-D-like H-mode n=1 modified Rutherford equation (MRE) behind the toroidal SLAYER gate.
+# A variant of diiid_slayer_n1_toroidal (its deck, overrides and quantities) plus [Tearing.ModifiedRutherford] with
+# polarization and saturation on, so the classical, saturation, bootstrap, Lütjens curvature and polarization terms,
+# the gate and the roots are tracked.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "diiid_slayer_n1_mre"
+base = "diiid_slayer_n1_toroidal"   # Inherit example_dir, [quantities.*] and [overrides] from this case
+description = "DIII-D-like H-mode, n=1, modified Rutherford equation with polarization and saturation behind the toroidal SLAYER gate"
+
+[overrides]
+"Tearing.ModifiedRutherford.polarization" = true   # Ion polarization term on (off by default)
+"Tearing.ModifiedRutherford.saturation" = true     # Finite-width saturation term on (off by default)
+
+# Surface coefficients the terms are built from.
+[quantities.mre_D_nc]
+h5path = "Tearing/ModifiedRutherford/PerSurface/D_nc"
+type = "real_vector"
+extract = "all_real"
+label = "MRE bootstrap coefficient D_nc"
+noise_threshold = 1e-8
+order = 10
+class = "physics_converged"
+
+[quantities.mre_C_sat]
+h5path = "Tearing/ModifiedRutherford/PerSurface/C_sat"
+type = "real_vector"
+extract = "all_real"
+label = "MRE saturation coefficient C_sat"
+noise_threshold = 1e-6
+order = 12
+class = "physics_converged"
+
+[quantities.mre_D_c]
+h5path = "Tearing/ModifiedRutherford/PerSurface/D_c"
+type = "real_vector"
+extract = "all_real"
+label = "MRE Lütjens critical offset D_c"
+noise_threshold = 1e-8
+order = 14
+class = "physics_converged"
+
+[quantities.mre_w_ib]
+h5path = "Tearing/ModifiedRutherford/PerSurface/w_ib"
+type = "real_vector"
+extract = "all_real"
+label = "MRE ion banana width w_ib (psi_N)"
+noise_threshold = 1e-10
+order = 16
+class = "physics_converged"
+
+# Roots: seed thresholds and saturated widths per surface (NaN where absent), and their minimum.
+[quantities.mre_w_m]
+h5path = "Tearing/ModifiedRutherford/PerSurface/w_m"
+type = "real_vector"
+extract = "all_real"
+label = "MRE seed threshold w_m (psi_N)"
+noise_threshold = 1e-8
+order = 20
+class = "physics_converged"
+
+[quantities.mre_w_sat]
+h5path = "Tearing/ModifiedRutherford/PerSurface/w_sat"
+type = "real_vector"
+extract = "all_real"
+label = "MRE saturated width w_sat (psi_N)"
+noise_threshold = 1e-8
+order = 22
+class = "physics_converged"
+
+[quantities.mre_w_m_star]
+h5path = "Tearing/ModifiedRutherford/w_m_star"
+type = "real_scalar"
+extract = "value"
+label = "MRE smallest seed threshold w_m* (psi_N)"
+noise_threshold = 1e-8
+order = 24
+class = "physics_converged"
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
+class = "diagnostic"