diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 0b94c9cb1..40f381ebc 100644 --- a/docs/development/hdf5-conventions.md +++ b/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/` | `/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `KineticMatrices/`); multi-ion runs add `PerSpecies///` 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}`, `Scan/Surface_/` | +| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_/`, `CriticalResonantField/` (+ `Scan/Surface_/`) | | `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`. diff --git a/docs/development/references.md b/docs/development/references.md index 8eb10b9e5..3ba308b7b 100644 --- a/docs/development/references.md +++ b/docs/development/references.md @@ -92,6 +92,15 @@ The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are - Published: Nuclear Fusion (2025), doi 10.1088/1741-4326/ae4fdd (open access: arXiv:2511.07666) - Describes: Diffusive-resistive layer width (Eq. 100) and the resistive-layer-overlap truncation `0 < Ψ < 1 − ε_c` (Sect. 5.9) +- **Cole & Fitzpatrick (2006)**: "Drift-magnetohydrodynamical model of error-field penetration in tokamak plasmas" + - Location: not in `docs/resources/` + - Published: Physics of Plasmas 13, 032503 (2006), doi 10.1063/1.2178167 + - Describes: Steady-state torque balance at a rational surface (Eq. 61) and the critical resonant field for error-field penetration (Eq. 62), used by `Tearing.CriticalResonantField` +- **Fitzpatrick (1993)**: "Interaction of tearing modes with external structures in cylindrical geometry" + - Location: not in `docs/resources/` + - Published: Nuclear Fusion 33, 1049 (1993) + - Describes: Bifurcated rotating states and the fold at which an error field penetrates (Sec. 7.4), the threshold definition used by `Tearing.CriticalResonantField` + ## Resistive MHD Stability Analysis (Future Work) GPEC will eventually implement resistive MHD stability analysis based on: diff --git a/docs/src/assets/crf_fortran_benchmark.png b/docs/src/assets/crf_fortran_benchmark.png new file mode 100644 index 000000000..80ce504eb Binary files /dev/null and b/docs/src/assets/crf_fortran_benchmark.png differ diff --git a/docs/src/citations.md b/docs/src/citations.md index f89debbd1..9bf7f5a03 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -152,6 +152,14 @@ The basis for the presently implemented in Fortran SLAYER code which computes th --- +> A. Cole and R. Fitzpatrick, "Drift-magnetohydrodynamical model of error-field penetration in tokamak plasmas," +> *Physics of Plasmas* **13**, 032503 (2006). +> DOI: [10.1063/1.2178167](https://doi.org/10.1063/1.2178167) + +The torque-balance error-field penetration threshold (Eqs. 61–62) behind the critical resonant field computed in `Tearing.CriticalResonantField`. + +--- + > D. A. Burgess, N. C. Logan, J.-K. Park, C. J. Hansen, and C. Paz-Soldan, "Tearing stability prediction combining toroidal calculations with a two-fluid slab layer," > *Plasma Physics and Controlled Fusion* **68**, 065017 (2026). > DOI: [10.1088/1361-6587/ae7156](https://doi.org/10.1088/1361-6587/ae7156) diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index ff3033c12..fa459f4a0 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -3,7 +3,8 @@ The `Tearing` module groups the resistive tearing-mode analysis stack: `InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and SLAYER models), `Dispersion` (physics-agnostic complex-plane scan and -contour-intersection root extraction), and `Runner` (user-facing TOML +contour-intersection root extraction), `CriticalResonantField` (torque-balance +error-field penetration threshold), and `Runner` (user-facing TOML configuration, profile loading, and HDF5 output). ## Layer Inputs @@ -201,6 +202,81 @@ res = solve_ray(p, GGJ.inner_Q(p, γ)) res.Δ, res.breaks, res.bc_cond # res.Δ = (Δ₊, Δ₋) = (interchange, tearing) ``` +## Critical resonant field + +With `[SLAYER.CriticalResonantField] enabled = true`, each SLAYER surface also gets the +critical resonant field for error-field penetration, from the steady-state torque balance of +Cole and Fitzpatrick, Phys. Plasmas **13**, 032503 (2006). The electromagnetic torque on the +layer balances the viscous restoring torque when (Cole Eq. 61) + +```math +\frac{\operatorname{Im}\hat\Delta(Q)}{|\alpha + \hat\Delta(Q)|^2} + = \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,(b_r/B_\phi)^2}, +\qquad +\hat\kappa = \left[\frac{2}{s}\right]^2 \int_{r_s}^{a} \frac{\mu(r_s)}{\mu(r)}\,\frac{dr}{r}, +``` + +and balance is lost above (Cole Eq. 62) + +```math +\left(\frac{b_r}{B_\phi}\right)^2_\mathrm{crit} + = \max_Q \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,\operatorname{Im}[-1/(\alpha + \hat\Delta(Q))]} . +``` + +- **Inputs.** P is each layer's `P_tor` (so χ_φ from the kinetic file, or the scalar + `chi_tor`), S its Lundquist number, and Q0 = τ_k·n·ω_E the natural E×B rotation of the + m/n mode (Cole's ω0; the diamagnetic drifts enter through Q_e and Q_i in Δ̂). Ω_E is the + kinetic file's `omega_E` or the `omega_E_kHz` override; the analysis needs it, and a surface + without one is skipped. +- **Approximations.** As in Fortran `gslayer.f`, the viscosity integral is taken as 1/2 + (κ̂ = 2/s²). +- **Stability α.** Each surface uses its own α = S^(−1/3)(Δ_crit − r_s Δ'), so α + Δ̂ is the + dispersion relation's residual and the torque balance shares its poles with SLAYER's free + modes (Fortran fixes α = 10⁻²). Cole and Fitzpatrick have no Δ_crit; including it treats the + curvature offset as independent of frequency. It is zero unless `D_c_type` is set. Cole's + threshold needs a tearing-stable surface, so a surface with α ≤ 0 is skipped with a warning. +- **Q axis.** Q is on Cole's axis, with the diamagnetic poles at Q_e and Q_i; `solve_inner` + mirrors the real axis, so the scan uses conj(Δ(−Q)). +- **Threshold.** Penetration is associated with a saddle-node bifurcation (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4): + as the field grows, Q moves from Q0 toward the electron fluid at Q_e, where the torque + vanishes, and the threshold is the first resolved saddle-node bifurcation on that path, above + which the shielded branch followed from Q0 no longer exists. Both the rise to and fall from + an interior maximum must exceed 0.1 % of its value within a contiguous retained interval. + This excludes monotonic clipped windows and small numerical extrema; excluded samples do + not connect separate intervals. A later, higher maximum does not replace the first resolved + saddle-node bifurcation. Points where |Δ̂| ≤ α, the + unshielded zone next to Q_e, are not candidates, so the balance spike there is never taken. + With no resolved saddle-node bifurcation no discrete threshold is reported (NaN). Fortran takes the global maximum on a + window instead, which picks that spike. The search assumes Q_i ≤ 0 ≤ Q_e; a surface where + that fails is skipped with a warning. + +Outputs go to `Tearing/CriticalResonantField/` (`br_crit` = b_r/B_φ, `br_crit_tesla`, `Q_peak`, `Q0`, with P in +`PerSurface/P_tor`, and with `store_scan` the per-surface `Scan/Surface_/{Q, balance, Delta}`). + +`br_crit` is normalized by the toroidal field, as in Cole and Fitzpatrick (2006, Eqs. 61–68) and +Fitzpatrick, Phys. Plasmas **30**, 092512 (2023), Eq. 120. `br_crit_tesla` = `br_crit`·|B0|, with +the layer's B0, is for comparison with `SingularCoupling/resonant_area_weighted_field`, the total +resonant field with the ideal plasma response, which Park et al., Phys. Rev. Lett. **99**, 195003 +(2007) identify as the field that opens islands. The comparison is approximate: the theory's b_r is +a cylindrical field at the layer, and no reference maps it onto the area-weighted toroidal field. +`penetrated_area_weighted_field` already includes the layer response and is not the matching quantity. + +### Benchmark against Fortran + +On the DIII-D-like example (2/1–6/1, n = 1, P = 1, κ̂ = 2/s², Q ∈ [−5, 5] with 20001 +points), both codes were given Fortran's normalized inputs and Fortran's fixed α = 10⁻². + +- **Torque balance.** Fed Fortran's own Δ(Q) table, the Julia routine reproduces Fortran's + balance to all printed digits on every surface, so the thresholds differ only through the + selection rule above. +- **Layer model.** 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 (7 % on 5/1). + +The dots in the figure mark each code's balance maximum between Q_e and Q0 at that fixed α, +not the fold threshold above. + +![Fortran vs Julia inner-layer Δ(Q) and torque balance](assets/crf_fortran_benchmark.png) + ## API Reference ### InnerLayer @@ -227,6 +303,12 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER] Modules = [GeneralizedPerturbedEquilibrium.Dispersion] ``` +## CriticalResonantField + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField] +``` + ## Runner ```@autodocs diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index d5f74e8da..3b033d2f8 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -102,3 +102,12 @@ Q_re_range = [-2.0, 2.0] # Scan box in the normalized Q plane, Re(Q) axis Q_im_range = [-0.5, 3.0] # Scan box in the normalized Q plane, Im(Q) axis nre = 41 # Grid resolution along the Re(Q) axis nim = 31 # Grid resolution along the Im(Q) axis + +[SLAYER.CriticalResonantField] +# Critical resonant field b_r/B_φ for error-field penetration at each SLAYER surface, from the +# Cole-Fitzpatrick 2006 torque balance with P = P_tor of each layer and Q0 = τ_k·n·Ω_E, so it +# needs the kinetic file's omega_E (or omega_E_kHz). b_crit is taken at the penetration fold met +# walking from Q0 toward Q_e; Qmin/Qmax may be set to override the sampled range. +enabled = true # run the analysis after SLAYER +n = 2000 # number of Q samples per surface +store_scan = false # write the per-surface Q, torque balance and Δ samples to gpec.h5 diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 90edf2e09..7a565c713 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -152,6 +152,40 @@ label = "SLAYER no_root flags [2/1,3/1,4/1]" noise_threshold = 0 order = 34 +# Critical resonant field (torque balance) for the inner three surfaces, matching the root pins +# above, and the natural rotation Q0 = τ_k·n·Ω_E it is computed from. +[quantities.crf_br_crit] +h5path = "Tearing/CriticalResonantField/br_crit" +type = "real_vector" +extract = "first_3" +label = "CRF b_r/B_φ crit [2/1,3/1,4/1]" +noise_threshold = 1e-9 +order = 40 + +[quantities.crf_br_crit_tesla] +h5path = "Tearing/CriticalResonantField/br_crit_tesla" +type = "real_vector" +extract = "first_3" +label = "CRF b_r crit in tesla [2/1,3/1,4/1]" +noise_threshold = 1e-9 +order = 43 + +[quantities.crf_Q_peak] +h5path = "Tearing/CriticalResonantField/Q_peak" +type = "real_vector" +extract = "first_3" +label = "CRF Q_peak [2/1,3/1,4/1]" +noise_threshold = 1e-6 +order = 41 + +[quantities.crf_Q0] +h5path = "Tearing/CriticalResonantField/Q0" +type = "real_vector" +extract = "first_3" +label = "CRF Q0 [2/1,3/1,4/1]" +noise_threshold = 1e-10 +order = 42 + # Settings (catches accidental config drift) [quantities.slayer_enabled] h5path = "Tearing/enabled" diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 43f6bd980..1b4b19da3 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -48,6 +48,7 @@ export Tearing # Backward-compat top-level aliases so callers can still reach these # directly; the canonical nested path is `Tearing.{Dispersion,Runner}`. import .Tearing.Dispersion as Dispersion +import .Tearing.CriticalResonantField as CriticalResonantField import .Tearing.Runner as Runner export Dispersion, Runner diff --git a/src/InnerLayer/SLAYER/Riccati.jl b/src/InnerLayer/SLAYER/Riccati.jl index 22c57e51c..f9951099c 100644 --- a/src/InnerLayer/SLAYER/Riccati.jl +++ b/src/InnerLayer/SLAYER/Riccati.jl @@ -225,6 +225,7 @@ function solve_inner(::SLAYERModel{:fitzpatrick}, # growth-rate roots are identical to those of Δ̂_s(ĝ); it only reflects # the off-root residual surface, fixing the integration-path orientation # so |Δ| contours are consistent across the scan plane. + # Real-Q callers see a mirrored axis; the critical resonant field maps back via conj(Δ(−Q)). Q_c = im * conj(ComplexF64(Q)) # Boundary condition at p_start diff --git a/src/Tearing/CriticalResonantField/CriticalResonantField.jl b/src/Tearing/CriticalResonantField/CriticalResonantField.jl new file mode 100644 index 000000000..65ef027a0 --- /dev/null +++ b/src/Tearing/CriticalResonantField/CriticalResonantField.jl @@ -0,0 +1,18 @@ +# CriticalResonantField.jl +# +# Uses the inner-layer delta from InnerLayer.jl to find the critical resonant field +# required for mode penetration at each rational surface. This can then be compared +# to the actual resonant field from the perturbed equilibrium to determine which modes +# are predicted to penetrate. + +module CriticalResonantField + +using Printf + +using ..InnerLayer: InnerLayerModel, solve_inner + +include("TorqueBalance.jl") + +export TorqueBalance, cole_delta, torque_balance_value, fold_search_range, torque_balance_scan + +end # module CriticalResonantField diff --git a/src/Tearing/CriticalResonantField/TorqueBalance.jl b/src/Tearing/CriticalResonantField/TorqueBalance.jl new file mode 100644 index 000000000..b9e45447c --- /dev/null +++ b/src/Tearing/CriticalResonantField/TorqueBalance.jl @@ -0,0 +1,140 @@ +# TorqueBalance.jl +# +# Error-field penetration threshold from the steady-state torque balance at a +# rational surface (Cole and Fitzpatrick, Phys. Plasmas 13, 032503 (2006), Sec. IV). +# The electromagnetic and viscous torques balance when (Cole Eq. 61) +# +# Im[Δ(Q)] / |α + Δ(Q)|² = 2·P·(Q0 − Q) / (S·κ̂·(b_r/B_φ)²), +# +# with κ̂ = [2/s(r_s)]² ∫_{r_s}^{a} [μ(r_s)/μ(r)] dr/r. Balance is lost when no Q +# satisfies this, so the critical field is (Cole Eq. 62) +# +# (b_r/B_φ)²_crit = max_Q 2·P·(Q0 − Q) / (S·κ̂·Im[−1/(α + Δ(Q))]). +# +# The viscosity integral is taken as 1/2 (κ̂ = 2/s²). α = S^(-1/3)·(−r_s Δ'_s) is each +# surface's own value (Cole Sec. IV C), which needs a tearing-stable surface (α > 0). +# Q is on Cole's axis: the electron and ion diamagnetic frequencies sit at Q = Q_e +# and Q = Q_i (see `cole_delta`). + + +""" + TorqueBalance{M<:InnerLayerModel,T} + +Torque-balance inputs at one rational surface. + +## Fields + + - `model` -- inner-layer model passed to `solve_inner` + - `params` -- that model's layer parameters at the surface; must carry `Q_e`, `Q_i`, `m`, `n` + - `Q0` -- normalized natural E×B rotation of the m/n mode, τ_k·n·ω_E + - `prandtl` -- magnetic Prandtl number τ_R/τ_V (the layer's P_φ) + - `lu` -- Lundquist number S + - `kappa_hat` -- Cole's slab-to-tokamak factor κ̂ (Eq. 61) + - `alpha` -- outer stability α = S^(-1/3)·(−r_s Δ'_s) (Eq. 61), positive on a tearing-stable surface +""" +struct TorqueBalance{M<:InnerLayerModel,T} + model::M + params::T + Q0::Float64 + prandtl::Float64 + lu::Float64 + kappa_hat::Float64 + alpha::Float64 +end + +""" + cole_delta(model, params, Q::Real) -> ComplexF64 + +Inner-layer Δ at real frequency `Q` on Cole's axis. `solve_inner` evaluates the layer +at `i·conj(Q)`, which mirrors the real axis, so Cole's Δ(Q) is `conj(Δ_solve_inner(−Q))`. +""" +cole_delta(model::InnerLayerModel, params, Q::Real) = + conj(solve_inner(model, params, ComplexF64(-Q)).tearing) + +""" + torque_balance_value(tb::TorqueBalance, Q::Real) -> (bal, Δ) + +Right-hand side of Cole Eq. 62 before the maximum, `bal = 2·P·(Q0 − Q) / Im[−1/(α + Δ)]`, +and the layer `Δ(Q)` it used. +""" +function torque_balance_value(tb::TorqueBalance, Q::Real) + Δ = cole_delta(tb.model, tb.params, Q) + jxb = -imag(1.0 / (Δ + tb.alpha)) + return 2.0 * tb.prandtl * (tb.Q0 - Q) / jxb, Δ +end + +# A resolved saddle-node bifurcation requires a rise and fall exceeding this fraction of the maximum. +const TORQUE_BALANCE_FOLD_DROP = 1e-3 + +""" + fold_search_range(Q0, Q_e, Q_i) -> (Qa, Qb) + +The segment the steady state follows as the error field grows: from the natural rotation `Q0` +toward the electron-fluid point `Q_e`, where the electromagnetic torque vanishes, ending at `Q_i` +instead if it lies in between (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4). +""" +function fold_search_range(Q0::Real, Q_e::Real, Q_i::Real) + return (Q0, min(Q0, Q_e) < Q_i < max(Q0, Q_e) ? Q_i : Q_e) +end + +""" + torque_balance_scan(tb::TorqueBalance; Qmin=nothing, Qmax=nothing, n=2000) + -> (Qs, bal, Qpeak, br_crit, idx_peak, Δs) + +Sample `torque_balance_value` on `n` uniform real Q points and return the critical normalized +field `br_crit = b_r/B_φ` (Cole Eq. 62) at the first resolved saddle-node bifurcation on the +path from `Q0` along `fold_search_range`. Both the rise to and fall from an interior maximum +must exceed `TORQUE_BALANCE_FOLD_DROP` times its value within a contiguous retained interval. +Above this threshold the shielded branch followed from `Q0` no longer exists as the field grows. +Points where |Δ| ≤ α, the unshielded zone next to Q_e, are not on the path. With no resolved +saddle-node bifurcation there is no discrete threshold. `Qmin`/`Qmax` override the sampled range, +but the search stays between `Q0` and the end of `fold_search_range`. +Returns NaN for `Qpeak`, `br_crit` and 0 for `idx_peak`, with a warning, when no saddle-node bifurcation is resolved, +when `Q0` is not finite, when the surface is tearing-unstable (α ≤ 0), or when the diamagnetic +ordering Q_i ≤ 0 ≤ Q_e does not hold. +""" +function torque_balance_scan(tb::TorqueBalance; Qmin=nothing, Qmax=nothing, n::Integer=2000) + p = tb.params + skip(why) = (@warn "CriticalResonantField: $(p.m)/$(p.n) surface skipped: $why"; (Float64[], Float64[], NaN, NaN, 0, ComplexF64[])) + isfinite(tb.Q0) || return skip("its natural rotation Q0 = $(tb.Q0) is not finite.") + tb.alpha > 0 || return skip("it is tearing-unstable or marginal (α = $(tb.alpha)), so there is no penetration threshold.") + p.Q_i <= 0 <= p.Q_e && p.Q_i < p.Q_e || + return skip("the fold search assumes Q_i ≤ 0 ≤ Q_e, but Q_e = $(p.Q_e), Q_i = $(p.Q_i).") + Qa, Qb = fold_search_range(tb.Q0, p.Q_e, p.Q_i) + Qs = range(something(Qmin, Qa), something(Qmax, Qb); length=n + 2)[2:(end-1)] # open at both ends + out = [torque_balance_value(tb, Q) for Q in Qs] + bal = first.(out) + Δs = last.(out) + + # Follow the first resolved interior maximum without joining samples across excluded intervals. + path = [i for i in sort(eachindex(Qs); by=i -> abs(Qs[i] - Qa)) + if min(Qa, Qb) < Qs[i] < max(Qa, Qb) && isfinite(bal[i]) && bal[i] > 0 && abs(Δs[i]) > tb.alpha] + idx = 0 + valley = peak = previous = 0 + for i in path + if previous == 0 || abs(i - previous) != 1 + valley = peak = i + elseif bal[i] > bal[peak] + peak = i + elseif bal[i] < (1 - TORQUE_BALANCE_FOLD_DROP) * bal[peak] + if bal[valley] < (1 - TORQUE_BALANCE_FOLD_DROP) * bal[peak] + idx = peak + break + end + valley = peak = i + elseif bal[i] < bal[valley] + valley = peak = i + end + previous = i + end + if idx == 0 + @warn "CriticalResonantField: no penetration saddle-node bifurcation between Q0 = $(tb.Q0) and Q = $Qb at the $(p.m)/$(p.n) " * + "surface: no interior maximum with a rise and fall exceeding $(100 * TORQUE_BALANCE_FOLD_DROP) % was resolved " * + "on a contiguous sampled interval, so no discrete threshold is resolved." + return Qs, bal, NaN, NaN, 0, Δs + end + br_crit = sqrt(bal[idx] / (tb.lu * tb.kappa_hat)) + @info @sprintf("CriticalResonantField: %d/%d surface b_r/B_φ crit = %.3e at Q = %.3f", + p.m, p.n, br_crit, Qs[idx]) + return Qs, bal, Qs[idx], br_crit, idx, Δs +end diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 493383660..a35f3ef9c 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -5,6 +5,44 @@ # constructor or by parsing the `[SLAYER]` (and nested `[SLAYER.*]`) # section(s) of a `gpec.toml`. +""" + CriticalResonantFieldControl + +Configuration for the critical resonant field (torque-balance) analysis, read from the +`[SLAYER.CriticalResonantField]` TOML section via `critical_resonant_field_control_from_toml` +or built with the `@kwdef` keyword constructor. The magnetic Prandtl number is the per-surface +`P_tor` of the SLAYER layer parameters, and the natural rotation is each surface's resolved Ω_E +(the kinetic file's `omega_E`, or the `omega_E_kHz` override), which it therefore requires. + +## Fields + + - `enabled` -- run the analysis after SLAYER + - `Qmin`, `Qmax` -- override the sampled real-Q range; `nothing` samples from Q0 toward Q_e per surface + - `n` -- number of Q samples per surface + - `store_scan` -- write the per-surface Q, torque balance and Δ samples to `gpec.h5` +""" +@kwdef struct CriticalResonantFieldControl + enabled::Bool = false + Qmin::Union{Nothing,Float64} = nothing + Qmax::Union{Nothing,Float64} = nothing + n::Int = 2000 + store_scan::Bool = false +end + +""" + critical_resonant_field_control_from_toml(section::AbstractDict) -> CriticalResonantFieldControl + +Parse a `[SLAYER.CriticalResonantField]` TOML section. Unknown keys raise an error. +""" +function critical_resonant_field_control_from_toml(section::AbstractDict) + field_names = Set(String.(fieldnames(CriticalResonantFieldControl))) + unknown = [k for k in keys(section) if !(k in field_names)] + isempty(unknown) || + throw(ArgumentError("critical_resonant_field_control_from_toml: unknown keys $(unknown) in " * + "[SLAYER.CriticalResonantField]. Known: $(sort(collect(field_names))).")) + return CriticalResonantFieldControl(; (Symbol(k) => v for (k, v) in section)...) +end + """ SLAYERControl @@ -113,6 +151,12 @@ there is one consistent interface for resistive and kinetic profiles. - `store_scan` -- write the full Q/Δ scan grid to HDF5. `false` by default to keep the output file small. + +# Critical resonant field + + - `critical_resonant_field` -- `CriticalResonantFieldControl` from the + `[SLAYER.CriticalResonantField]` subsection (keyword form: + `critical_resonant_field=CriticalResonantFieldControl(; enabled=true)`) """ @kwdef struct SLAYERControl enabled::Bool = false @@ -180,6 +224,8 @@ there is one consistent interface for resistive and kinetic profiles. profile_group::String = "/" store_scan::Bool = false + + critical_resonant_field::CriticalResonantFieldControl = CriticalResonantFieldControl() end const _VALID_INNER_MODELS = (:slayer_fitzpatrick,) @@ -230,6 +276,8 @@ function validate(ctrl::SLAYERControl) all(isfinite, values(ctrl.omega_E_kHz)) || throw(ArgumentError("SLAYERControl: omega_E_kHz contains a " * "non-finite entry: $(ctrl.omega_E_kHz)")) + crf = ctrl.critical_resonant_field + crf.n >= 3 || throw(ArgumentError("SLAYERControl: CriticalResonantField n must be ≥ 3")) ctrl.nre >= 2 && ctrl.nim >= 2 || throw(ArgumentError("SLAYERControl: nre and nim must both be ≥ 2")) ctrl.amr_passes >= 0 || @@ -269,6 +317,8 @@ function slayer_control_from_toml(section::AbstractDict) haskey(v, "pole_threshold") && (flat["pole_threshold"] = v["pole_threshold"]) haskey(v, "filter_above_poles") && (flat["filter_above_poles"] = v["filter_above_poles"]) haskey(v, "filter_outside_re") && (flat["filter_outside_re"] = v["filter_outside_re"]) + elseif k == "CriticalResonantField" && v isa AbstractDict + flat["critical_resonant_field"] = critical_resonant_field_control_from_toml(v) else flat[k] = v end diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index 013118e32..ab7914191 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -15,7 +15,8 @@ # ├── Roots/ -- complex Q_root, omega_Hz, gamma_Hz # ├── Diagnostics/ -- ValidRoots, Poles, FilteredRoots (uncoupled: each surface's E×B frame) # │ (flat-plus-offsets ragged encoding) -# └── Scan/ -- optional: full Q/Δ scan data +# ├── Scan/ -- optional: full Q/Δ scan data +# └── CriticalResonantField/ -- optional: torque-balance b_r/B_φ crit per surface using HDF5 @@ -57,10 +58,30 @@ function write_slayer_hdf5!(parent::Union{HDF5.File,HDF5.Group}, if result.control.store_scan && !isempty(result.scan_data) _write_scan_data!(g, result) end + result.critical_resonant_field.enabled && _write_critical_resonant_field!(g, result.critical_resonant_field) _annotate_tearing!(g) return g end +# ---------- critical resonant field ---------- +function _write_critical_resonant_field!(g, crf::CriticalResonantFieldResult) + cg = create_group(g, "CriticalResonantField") + cg["rational_index"] = crf.rational_index + cg["Q_peak"] = crf.Q_peak + cg["br_crit"] = crf.br_crit + cg["br_crit_tesla"] = crf.br_crit_tesla + cg["Q0"] = crf.Q0 + isempty(crf.scan) && return nothing + scan = create_group(cg, "Scan") + for (k, sc) in enumerate(crf.scan) + sg = create_group(scan, "Surface_$k") + sg["Q"] = sc.Q + sg["balance"] = sc.balance + sg["Delta"] = sc.Delta + end + return nothing +end + # Token recorded in the Tearing group's layer_model attribute; keyed by the # per-surface parameter type since the SLAYER and GGJ branches write disjoint fields. _layer_model_token(::Type{SLAYERParameters}) = "slayer" @@ -138,7 +159,24 @@ const TEARING_H5_ANNOTATIONS = [ "LayerWidths/delta_s_over_d_beta" => (; long_name="complex dimensionless layer thickness δ_s/d_β", dims=("surface",)), "LayerWidths/delta_s" => (; long_name="complex resistive layer thickness δ_s (Riccati)", dims=("surface",)), "LayerWidths/delta_s_abs" => (; long_name="physical resistive layer thickness |δ_s|", units="m", dims=("surface",)), - "LayerWidths/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)) + "LayerWidths/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)), + "CriticalResonantField/rational_index" => (; long_name="rational-surface index of each row", dims=("surface",)), + "CriticalResonantField/Q_peak" => + (; long_name="normalized frequency Q (Cole axis) at the torque-balance maximum (NaN = no maximum or surface skipped)", units="1", dims=("surface",)), + "CriticalResonantField/br_crit" => + (; long_name="critical normalized resonant field b_r/B_φ for error-field penetration (Cole-Fitzpatrick 2006 Eq. 62; NaN = no maximum or surface skipped)", + units="1", dims=("surface",)), + "CriticalResonantField/br_crit_tesla" => + (; long_name="critical resonant field br_crit·|B0| with the layer's B0, approximately comparable to SingularCoupling/resonant_area_weighted_field (NaN = no maximum or surface skipped)", + units="T", dims=("surface",)), + "CriticalResonantField/Q0" => + (; long_name="normalized natural E×B rotation Q0 = τ_k·n·Ω_E of the m/n mode (NaN where Ω_E is unknown)", units="1", dims=("surface",)) +] + +const TEARING_CRF_SCAN_H5_ANNOTATIONS = [ + "Q" => (; long_name="sampled real normalized frequency Q (Cole axis)", units="1"), + "balance" => (; long_name="torque balance 2·P·(Q0 − Q)/Im[−1/(α + Δ)]", units="1"), + "Delta" => (; long_name="complex inner-layer Δ(Q) on the Cole axis", units="1") ] const TEARING_RAGGED_H5_ANNOTATIONS = [ @@ -165,6 +203,11 @@ function _annotate_tearing!(g) ann.annotate!(g["Diagnostics"][sub], TEARING_RAGGED_H5_ANNOTATIONS) end end + if haskey(g, "CriticalResonantField/Scan") + for sub in keys(g["CriticalResonantField/Scan"]) + ann.annotate!(g["CriticalResonantField/Scan"][sub], TEARING_CRF_SCAN_H5_ANNOTATIONS) + end + end if haskey(g, "Scan") for sub in keys(g["Scan"]) sg = g["Scan"][sub] diff --git a/src/Tearing/Runner/Result.jl b/src/Tearing/Runner/Result.jl index 50acacbd5..bbba3e61f 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -4,6 +4,52 @@ # per-surface layer parameters, the extracted tearing eigenvalues, and (if # `control.store_scan`) the full Q-plane scan data for plotting. +""" + CriticalResonantFieldScan + +Real-Q torque-balance samples at one surface. + +## Fields + + - `Q` -- sampled normalized frequencies (Cole's axis) + - `balance` -- torque-balance value 2·P·(Q0 − Q)/Im[−1/(α + Δ)] at each Q + - `Delta` -- inner-layer Δ(Q) at each Q +""" +struct CriticalResonantFieldScan + Q::Vector{Float64} + balance::Vector{Float64} + Delta::Vector{ComplexF64} +end + +""" + CriticalResonantFieldResult + +Output of `run_critical_resonant_field`, one entry per SLAYER surface. + +## Fields + + - `enabled` -- the analysis ran + - `rational_index` -- rational-surface index of each entry + - `Q_peak` -- normalized frequency Q at the torque-balance maximum (NaN if none or skipped) + - `br_crit` -- critical normalized resonant field b_r/B_φ (NaN if none or skipped) + - `br_crit_tesla` -- `br_crit`·|B0| in tesla, with the layer's B0, approximately comparable to the total resonant field + `resonant_area_weighted_field` (Park et al., PRL 99, 195003 (2007)) + - `Q0` -- normalized natural E×B rotation τ_k·n·Ω_E (NaN where Ω_E is unknown) + - `scan` -- per-surface Q scans; empty unless `store_scan` +""" +struct CriticalResonantFieldResult + enabled::Bool + rational_index::Vector{Int} + Q_peak::Vector{Float64} + br_crit::Vector{Float64} + br_crit_tesla::Vector{Float64} + Q0::Vector{Float64} + scan::Vector{CriticalResonantFieldScan} +end + +empty_critical_resonant_field_result() = + CriticalResonantFieldResult(false, Int[], Float64[], Float64[], Float64[], Float64[], CriticalResonantFieldScan[]) + """ SLAYERResult @@ -48,6 +94,8 @@ downstream inspection and HDF5 output. inner-layer Q; zero outside `:coupled` mode, where rotation does not enter - `scan_data` -- scan results (per-surface in uncoupled, single entry in coupled). Empty unless `control.store_scan == true`. + - `critical_resonant_field` -- `CriticalResonantFieldResult`; `enabled=false` + unless `control.critical_resonant_field.enabled` """ struct SLAYERResult enabled::Bool @@ -66,6 +114,7 @@ struct SLAYERResult omega_E::Vector{Float64} q_shift::Vector{Float64} scan_data::Vector{Union{ScanResult,AMRResult}} + critical_resonant_field::CriticalResonantFieldResult end # Empty result (enabled=false path) @@ -78,5 +127,6 @@ function empty_slayer_result(control::SLAYERControl) GrowthRateResult[], nothing, LayerWidths[], Float64[], Float64[], - Union{ScanResult,AMRResult}[]) + Union{ScanResult,AMRResult}[], + empty_critical_resonant_field_result()) end diff --git a/src/Tearing/Runner/Runner.jl b/src/Tearing/Runner/Runner.jl index 152bce0a2..f02de563e 100644 --- a/src/Tearing/Runner/Runner.jl +++ b/src/Tearing/Runner/Runner.jl @@ -45,6 +45,7 @@ using ..Dispersion: SurfaceCoupling, surface_coupling, AMRResult, amr_scan, MultiBoxAMRResult, multi_box_amr_scan, as_amr_result, GrowthRateResult, find_growth_rates +using ..CriticalResonantField: TorqueBalance, torque_balance_scan include("Control.jl") include("Result.jl") @@ -56,5 +57,8 @@ export SLAYERControl, slayer_control_from_toml, validate export SLAYERResult, empty_slayer_result export TearingProblem, run_slayer_from_inputs, ggj_inner_deltas export write_slayer_hdf5! +export CriticalResonantFieldControl, critical_resonant_field_control_from_toml +export CriticalResonantFieldScan, CriticalResonantFieldResult, empty_critical_resonant_field_result +export run_critical_resonant_field end # module Runner diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 72eed83f4..cf07e5954 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -480,7 +480,51 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel, return SLAYERResult(true, control, model, params, rational_psi, rational_q, dp, Q_root, omega_Hz, gamma_Hz, per_surface_extraction, coupled_extraction, - layer_widths, Ω_E, q_shifts, scan_data_list) + layer_widths, Ω_E, q_shifts, scan_data_list, + run_critical_resonant_field(model, params, Ω_E, real.(diag(dp)), control.critical_resonant_field)) +end + +# --------------------------------------------------------------------- +# Critical resonant field (torque balance) +# --------------------------------------------------------------------- +""" + run_critical_resonant_field(model, params, omega_E, dp_diag, ctrl) -> CriticalResonantFieldResult + +Critical normalized resonant field b_r/B_φ for error-field penetration at each SLAYER +surface, from the torque balance of Cole and Fitzpatrick, Phys. Plasmas 13, 032503 (2006), +Eq. 62. The magnetic Prandtl number is the layer's `P_tor`, and the natural rotation is the +E×B frequency of the m/n mode, Q0 = τ_k·n·Ω_E, with `omega_E` the resolved Ω_E per unit n at +each surface [rad/s]. The outer stability is α = S^(-1/3)·(Δ_crit − r_s Δ'), with `dp_diag` the +r_s-referenced Δ' diagonal and Δ_crit the layer's critical-Δ offset, as in the dispersion relation. +A surface whose Ω_E is unknown (NaN), or that is tearing-unstable (α ≤ 0), is skipped with a +warning; if no surface has a known Ω_E, the analysis cannot run and this throws. +""" +function run_critical_resonant_field(model::InnerLayer.InnerLayerModel, params::AbstractVector{<:InnerLayerParameters}, + omega_E::AbstractVector{<:Real}, dp_diag::AbstractVector{<:Real}, ctrl::CriticalResonantFieldControl) + ctrl.enabled || return empty_critical_resonant_field_result() + model isa SLAYERModel{:fitzpatrick} && eltype(params) <: SLAYERParameters || + throw(ArgumentError("run_critical_resonant_field: the critical resonant field needs the SLAYER " * + "fitzpatrick model and SLAYER parameters; got $(typeof(model)).")) + nsurf = length(params) + length(dp_diag) == nsurf || + throw(ArgumentError("run_critical_resonant_field: dp_diag has $(length(dp_diag)) entries for $nsurf surfaces.")) + length(omega_E) == nsurf || + throw(ArgumentError("run_critical_resonant_field: omega_E has $(length(omega_E)) entries for $nsurf surfaces.")) + nsurf > 0 && all(isnan, omega_E) && + throw(ArgumentError("run_critical_resonant_field: the critical resonant field needs the E×B rotation, but the " * + "kinetic file has no omega_E and no [SLAYER] omega_E_kHz entry covers the surfaces.")) + Q_peak, br_crit = fill(NaN, nsurf), fill(NaN, nsurf) + # Cole Sec. IV: ω0 is the mode frequency in the E×B frame; diamagnetic parts enter via Q_e, Q_i. + Q0 = Float64[p.tauk * p.n * Ω for (p, Ω) in zip(params, omega_E)] + scan = CriticalResonantFieldScan[] + for (k, p) in enumerate(params) + alpha = (p.dc_tmp - dp_diag[k]) / p.lu^(1 / 3) + tb = TorqueBalance(model, p, Q0[k], p.P_tor, p.lu, 2.0 / p.sval_r^2, alpha) + Qs, bal, Q_peak[k], br_crit[k], _, Δs = torque_balance_scan(tb; Qmin=ctrl.Qmin, Qmax=ctrl.Qmax, n=ctrl.n) + ctrl.store_scan && push!(scan, CriticalResonantFieldScan(collect(Float64, Qs), bal, Δs)) + end + br_crit_tesla = br_crit .* Float64[abs(p.bt) for p in params] + return CriticalResonantFieldResult(true, Int[p.ising for p in params], Q_peak, br_crit, br_crit_tesla, Q0, scan) end # --------------------------------------------------------------------- diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 41055f9e6..834e48253 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -6,6 +6,8 @@ # InnerLayer -- pure physics: Δ_inner(Q) for GGJ or SLAYER models # Dispersion -- physics-agnostic scan + contour-intersection root # extraction (consumes any InnerLayerModel) +# CriticalResonantField -- real-Q torque-balance scan for the critical +# resonant field of error-field penetration # Runner -- user-facing orchestration: TOML config, profile # loading, HDF5 output, workflow hooks # @@ -22,12 +24,14 @@ import ..InnerLayer as InnerLayer include("LayerOverlap.jl") include("Dispersion/Dispersion.jl") +include("CriticalResonantField/CriticalResonantField.jl") include("Runner/Runner.jl") import .Dispersion as Dispersion +import .CriticalResonantField as CriticalResonantField import .Runner as Runner -export InnerLayer, Dispersion, Runner +export InnerLayer, Dispersion, CriticalResonantField, Runner export resistive_layer_overlap, LayerOverlapScan end # module Tearing diff --git a/test/runtests.jl b/test/runtests.jl index 20ce7fd4c..c52bbf8ab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,7 @@ else include("./runtests_dispersion_polish.jl") include("./runtests_slayer_profile_transport.jl") include("./runtests_slayer_runner.jl") + include("./runtests_critical_resonant_field.jl") include("./runtests_kinetic.jl") include("./runtests_kinetic_validity.jl") include("./runtests_multiion.jl") diff --git a/test/runtests_critical_resonant_field.jl b/test/runtests_critical_resonant_field.jl new file mode 100644 index 000000000..a285fe3a4 --- /dev/null +++ b/test/runtests_critical_resonant_field.jl @@ -0,0 +1,209 @@ +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.InnerLayer +using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel + +# Layer stub returning a prescribed torque balance `target(Q)`, to test the maximum selection. +struct _CRFStubLayer{F} <: InnerLayerModel + target::F +end +const _CRF_Q0_STUB = 10.0 +_crf_target(Q) = 0.1 + exp(-(Q - 1)^2 / 0.01) + 2exp(-(Q - 3)^2 / 0.01) +_CRFStubLayer() = _CRFStubLayer(_crf_target) +# A fold at Q = 3, then a dip, then a rise past the fold's height toward Q_e = 0.5. +_crf_target_rising(Q) = 0.1 + exp(-(Q - 3)^2 / 0.01) + 0.6 * max(2.5 - Q, 0.0) +function InnerLayer.solve_inner(s::_CRFStubLayer, p, Q::Number) + Qc = -real(Q) # undo the axis mirror cole_delta applies + jxb = 2 * (_CRF_Q0_STUB - Qc) / s.target(Qc) + return (; tearing=conj(im / jxb - 1e-2)) +end + +# Cole and Fitzpatrick (2006) Table I visco-resistive layer on Cole's axis, Δ̂ = 2.104[i(Q−Q_i)]^(1/6)[i(Q−Q_e)]^(5/6)P^(1/6). +struct _VRiStub <: InnerLayerModel end +const _VRI_P = 4.0 +_vri(Q, Q_e, Q_i) = 2.104 * (im * (Q - Q_i))^(1 / 6) * (im * (Q - Q_e))^(5 / 6) * _VRI_P^(1 / 6) +InnerLayer.solve_inner(::_VRiStub, p, Q::Number) = (; tearing=conj(_vri(-real(Q), p.Q_e, p.Q_i))) + +@testset "CriticalResonantField: torque balance" begin + using GeneralizedPerturbedEquilibrium + using GeneralizedPerturbedEquilibrium.InnerLayer + using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel, SLAYERModel, SLAYERParameters, GGJModel + using GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField + using GeneralizedPerturbedEquilibrium.Runner + using HDF5 + + include("h5_metadata_check.jl") + + # DIII-D-like 2/1 layer (Fortran-normalized inputs of the shipped SLAYER example). + _mk(; n=1, P_tor=1.0, Q_e=1.38, Q_i=-2.15, ising=1) = SLAYERParameters(; + tau=1.2, lu=6.44e7, c_beta=0.112, D_norm=3.0, P_perp=1.0, P_tor=P_tor, + Q_e=Q_e, Q_i=Q_i, iota_e=Q_e / (Q_e - Q_i), tauk=1.31e-4, tau_r=21.1, + delta_n=885.0, rs=0.453, R0=1.74, bt=2.0, sval_r=1.09, eta=1.22e-8, + d_beta=3.6e-3, m=2, n=n, ising=ising) + model = SLAYERModel{:fitzpatrick}() + + @testset "Cole axis: branch lies between Q_e and Q0" begin + p = _mk() + @test cole_delta(model, p, 0.7) ≈ conj(InnerLayer.solve_inner(model, p, -0.7).tearing) + Q0 = 5.45 + tb = TorqueBalance(model, p, Q0, 1.0, p.lu, 2 / p.sval_r^2, 1e-2) + Qs, bal, Qpeak, br, idx, _ = torque_balance_scan(tb; n=401) + # Cole's picture: with Q0 above the electron pole, balance holds only for Q_e < Q < Q0. + @test p.Q_e < Qpeak < Q0 + @test bal[idx] > 0 && isfinite(br) && br > 0 + @test br ≈ sqrt(bal[idx] / (p.lu * 2 / p.sval_r^2)) + # Refining the grid moves b_crit by well under a percent. + _, _, _, br_fine, _, _ = torque_balance_scan(tb; n=801) + @test isapprox(br_fine, br; rtol=5e-3) + end + + @testset "saddle-node selection follows the path from Q0 toward Q_e" begin + # The stub's balance peaks at Q = 3 and Q = 1; walking down from Q0 = 10 the fold is the higher one at 3. + p = _mk(; Q_e=0.5, Q_i=-50.0) + tb = TorqueBalance(_CRFStubLayer(), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, bal, Qpeak, br, idx, _ = torque_balance_scan(tb; n=4001) + @test Qpeak ≈ 3.0 atol = 0.005 + @test bal[idx] ≈ _crf_target(Qpeak) rtol = 1e-8 + # A wider sampled range does not extend the search past Q_e. + tb_e2 = TorqueBalance(_CRFStubLayer(), _mk(; Q_e=2.0, Q_i=-50.0), _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + @test torque_balance_scan(tb_e2; Qmin=0.0, Qmax=4.0, n=401)[3] ≈ 3.0 atol = 0.011 + # With the electron fluid at Q = 3 the balance rises all the way along the path: no fold, no threshold. + tb_e3 = TorqueBalance(_CRFStubLayer(), _mk(; Q_e=3.0, Q_i=-50.0), _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, _, Qnone, brnone, inone, _ = @test_logs (:warn, r"no penetration") torque_balance_scan(tb_e3; n=401) + @test isnan(Qnone) && isnan(brnone) && inone == 0 + # The first fold from Q0 is the threshold even when the balance later climbs higher toward Q_e. + tb_r = TorqueBalance(_CRFStubLayer(_crf_target_rising), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, bal_r, Qpeak_r, _, idx_r, _ = torque_balance_scan(tb_r; n=4001) + @test Qpeak_r ≈ 3.0 atol = 0.005 + @test maximum(bal_r) > bal_r[idx_r] + # The path ends at Q_i when it lies between Q0 and Q_e. + @test fold_search_range(-2.0, 1.0, -0.5) == (-2.0, -0.5) + @test fold_search_range(5.0, 1.0, -0.5) == (5.0, 1.0) + end + + @testset "saddle-node selection requires a resolved interior maximum" begin + p = _mk(; Q_e=0.5, Q_i=-50.0) + targets = ( + Q -> Q, + Q -> 10.0 - Q, + Q -> 1.0, + Q -> 1.0 + 5e-4 * exp(-(Q - 3)^2 / 0.01) - 0.1 * max(2.8 - Q, 0.0) + ) + # Monotonic windows, flat windows and a subthreshold bump followed by a large fall have no resolved saddle-node bifurcation. + for target in targets, (Qmin, Qmax) in ((4.0, 2.0), (2.0, 4.0)) + tb = TorqueBalance(_CRFStubLayer(target), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, _, Qpeak, br, idx, _ = @test_logs (:warn, r"no penetration") torque_balance_scan(tb; Qmin=Qmin, Qmax=Qmax, n=401) + @test isnan(Qpeak) && isnan(br) && idx == 0 + end + + # An initial decline does not prevent resolving a later interior maximum. + target = Q -> 1.0 + 0.5 * exp(-(Q - 2)^2 / 0.04) + 0.2 * max(Q - 3.0, 0.0) + tb = TorqueBalance(_CRFStubLayer(target), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, _, Qpeak, br, _, _ = torque_balance_scan(tb; Qmin=4.0, Qmax=1.0, n=401) + @test Qpeak ≈ 2.0 atol = 0.01 + @test isfinite(br) && br > 0 + + # A missing interval must not turn a rising segment and a lower plateau into a maximum. + target = Q -> 2.9 < Q < 3.1 ? NaN : Q > 3.0 ? 1.0 + 0.1 * (4.0 - Q) : 1.0 + tb = TorqueBalance(_CRFStubLayer(target), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, _, Qpeak, br, idx, _ = @test_logs (:warn, r"no penetration") torque_balance_scan(tb; Qmin=4.0, Qmax=2.0, n=401) + @test isnan(Qpeak) && isnan(br) && idx == 0 + + # A later, higher maximum does not replace the first saddle-node bifurcation. + target = Q -> 0.1 + exp(-(Q - 3)^2 / 0.01) + 2 * exp(-(Q - 1)^2 / 0.01) + tb = TorqueBalance(_CRFStubLayer(target), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0, 1e-2) + _, bal, Qpeak, _, idx, _ = torque_balance_scan(tb; n=4001) + @test Qpeak ≈ 3.0 atol = 0.005 + @test maximum(bal) > bal[idx] + end + + @testset "fold reproduces Cole's visco-resistive thresholds" begin + # E×B rotation only (Cole Eq. 63): fold at Q0/2, where b_crit·√S/P^(7/12) = √(2.104/2) ≈ 1.03 with κ̂ = 1. + p = _mk(; Q_e=1e-9, Q_i=-1e-9) + tb = TorqueBalance(_VRiStub(), p, 1.0, _VRI_P, p.lu, 1.0, 1e-6) + _, _, Qpeak, br, _, _ = torque_balance_scan(tb; n=20001) + @test Qpeak ≈ 0.5 atol = 1e-3 + @test br * sqrt(p.lu) / _VRI_P^(7 / 12) ≈ sqrt(2.104 / 2) rtol = 1e-4 + # Diamagnetic rotation only (Cole Eq. 66): fold at 7Q_e/12, coefficient 1.11. + p = _mk(; Q_e=1.0, Q_i=-1e-9) + tb = TorqueBalance(_VRiStub(), p, 0.0, _VRI_P, p.lu, 1.0, 1e-6) + _, _, Qpeak, br, _, _ = torque_balance_scan(tb; n=20001) + @test Qpeak ≈ 7 / 12 atol = 1e-3 + @test br * sqrt(p.lu) / _VRI_P^(7 / 12) ≈ 1.11 rtol = 1e-2 + end + + @testset "surfaces without a usable fold search are skipped" begin + # Unknown rotation: no Q0, so no threshold. + tb_nan = TorqueBalance(model, _mk(), NaN, 1.0, 6.44e7, 1.0, 1e-2) + _, _, Qn, brn, in_, _ = @test_logs (:warn, r"not finite") torque_balance_scan(tb_nan; n=51) + @test isnan(Qn) && isnan(brn) && in_ == 0 + # Reversed diamagnetic ordering (Q_e < 0 < Q_i) is not guessed at. + tb_rev = TorqueBalance(model, _mk(; Q_e=-1.38, Q_i=2.15), -0.5, 1.0, 6.44e7, 1.0, 1e-2) + _, _, Qr, brr, ir, _ = @test_logs (:warn, r"Q_i ≤ 0 ≤ Q_e") torque_balance_scan(tb_rev; n=51) + @test isnan(Qr) && isnan(brr) && ir == 0 + end + + @testset "runner: P from P_tor, Q0 = τ_k·n·Ω_E" begin + ω_E = 4.0e4 + params = [_mk(; n=1, P_tor=3.0, ising=1), _mk(; n=2, P_tor=5.0, ising=2)] + ctrl = CriticalResonantFieldControl(; enabled=true, n=201, store_scan=true) + r = run_critical_resonant_field(model, params, [ω_E, ω_E], [-2.0, -3.0], ctrl) + @test r.enabled + @test r.Q0 ≈ [p.tauk * p.n * ω_E for p in params] + @test r.rational_index == [1, 2] + # α is each surface's own S^(-1/3)·(Δ_crit − r_s Δ'); a tearing-unstable surface has no threshold. + r_unst = @test_logs (:warn, r"2/2 surface skipped: it is tearing-unstable") match_mode = :any run_critical_resonant_field( + model, params, [ω_E, ω_E], [-2.0, 0.5], ctrl) + @test isequal(r_unst.br_crit[1], r.br_crit[1]) && isnan(r_unst.br_crit[2]) + # Zero E×B rotation is a valid input (Cole Sec. IV E: purely diamagnetic rotation) and gives a finite threshold. + r_zero = run_critical_resonant_field(model, params[1:1], [0.0], [-2.0], ctrl) + @test r_zero.Q0 == [0.0] && isfinite(r_zero.br_crit[1]) && r_zero.br_crit[1] > 0 + @test length(r.scan) == 2 && length(r.scan[1].Q) == 201 + @test !run_critical_resonant_field(model, params, [ω_E, ω_E], [-2.0, -3.0], CriticalResonantFieldControl()).enabled + # Only the SLAYER layer defines the Cole-axis mirror the scan relies on. + @test_throws "needs the SLAYER" run_critical_resonant_field(GGJModel(), params, [ω_E, ω_E], [-2.0, -3.0], ctrl) + + # The analysis runs inside run_slayer_from_inputs, so script callers get it too, and its + # HDF5 output honours the metadata contract. + c = SLAYERControl(; enabled=true, scan_mode=:brute_force, nre=4, nim=4, critical_resonant_field=ctrl) + dp = ComplexF64[-2.0 0.0; 0.0 -3.0] + res = run_slayer_from_inputs(model, params, dp, c; rational_psi=[0.5, 0.7], rational_q=[2.0, 1.5], omega_E=[ω_E, ω_E]) + @test res.critical_resonant_field.enabled + @test isequal(res.critical_resonant_field.br_crit, r.br_crit) + mktemp() do path, io + close(io) + h5open(path, "w") do f + write_slayer_hdf5!(f, res) + end + h5open(path, "r") do f + @test isempty(_collect_metadata_violations(f)) + g = f["Tearing/CriticalResonantField"] + @test read(g["Q0"]) ≈ r.Q0 + @test read(g["br_crit"]) ≈ r.br_crit nans = true + @test read(g["br_crit_tesla"]) ≈ r.br_crit .* [abs(p.bt) for p in params] nans = true + @test !haskey(g, "p_phi") + @test haskey(g, "Scan/Surface_2/Delta") + end + end + + # Without any rotation there is no threshold to compute, so the run stops rather than assume Ω_E = 0. + @test_throws "needs the E×B rotation" run_slayer_from_inputs(model, params, dp, c; rational_psi=[0.5, 0.7], rational_q=[2.0, 1.5]) + # A surface the override leaves without a rotation is skipped; the others are still computed. + c_part = SLAYERControl(; enabled=true, scan_mode=:brute_force, nre=4, nim=4, critical_resonant_field=ctrl, + omega_E_kHz=Dict("2/1" => ω_E / (2π * 1e3))) + r_part = (@test_logs (:warn, r"2/2 surface skipped") match_mode = :any run_slayer_from_inputs(model, params, dp, c_part; + rational_psi=[0.5, 0.7], rational_q=[2.0, 1.5])).critical_resonant_field + @test r_part.Q0[1] ≈ r.Q0[1] && isnan(r_part.Q0[2]) + @test isequal(r_part.br_crit[1], r.br_crit[1]) && isnan(r_part.br_crit[2]) + end + + @testset "control: TOML parsing and validation" begin + c = critical_resonant_field_control_from_toml(Dict("enabled" => true, "n" => 500)) + @test c.enabled && c.n == 500 && c.Qmin === nothing && !c.store_scan + @test critical_resonant_field_control_from_toml(Dict("Qmin" => -2.0)).Qmin == -2.0 + @test_throws ArgumentError critical_resonant_field_control_from_toml(Dict("pole_width" => 0.02)) + @test_throws ArgumentError critical_resonant_field_control_from_toml(Dict("viscous_input" => 1.0)) + s = slayer_control_from_toml(Dict("enabled" => true, "CriticalResonantField" => Dict("enabled" => true))) + @test s.critical_resonant_field.enabled + @test_throws "n must be ≥ 3" Runner.validate(SLAYERControl(; critical_resonant_field=CriticalResonantFieldControl(; n=2))) + end +end