diff --git a/Project.toml b/Project.toml index a758c3b5a..9b8e5098a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeneralizedPerturbedEquilibrium" uuid = "462872dd-e066-4d2e-b993-6468b5239634" license = "MIT" -authors = ["Nikolas Logan ", "Jong-Kyu Park ", "Matthew Pharr ", "Jacob Halpern ", "Rithik Banerjee ", "Jaebeom Cho ", "Daniel Burgess ", "Min-Gu Yoo "] +authors = ["Nikolas Logan ", "Jong-Kyu Park ", "Matthew Pharr ", "Jacob Halpern ", "Rithik Banerjee ", "Jaebeom Cho ", "Daniel Burgess ", "Min-Gu Yoo ", "Evan Bursch "] version = "0.1.0" [deps] diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index c4631c984..692407c96 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/`, control-surface spectra | | `KineticForces/` | `/` (torque/energy profiles, `EnergyIntegrals/`, `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/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 e78370553..076daa527 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -124,6 +124,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`. + +--- + > A. Burgess et al., "Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer," > Preprint (2026). diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index 7ddfecb95..20fc32bb9 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 @@ -278,6 +279,60 @@ res = solve_ray(p, GGJ.inner_Q(p, γ)) res.Δ, res.resid, res.bc_cond ``` +## 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 Δ̂). +- **Approximations.** As in Fortran `gslayer.f`, the viscosity integral is taken as 1/2 + (κ̂ = 2/s²) and α = 10⁻² stands in for S^(−1/3)(−r_s Δ'_s) ≪ 1 (b_crit moves about 4 % + over 10⁻⁴–10⁻¹). +- **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)). +- **Maximum.** The scan window follows `gslayer.f` (between Q0 and the Q_e pole), and the + largest positive interior local maximum is taken, skipping any within 0.02 of Q_e or Q_i + (Fortran SLAYER's pole-regularization radius). + +Outputs go to `Tearing/CriticalResonantField/` (`br_crit` = b_r/B_φ, `q_peak`, `q0`, `p_phi`, +and with `store_scan` the per-surface `Scan/Surface_/{Q, balance, Delta}`). + +### 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. + +- **Torque-balance routine.** Fed Fortran's own Δ(Q) 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 every surface. +- **Pole handling.** Fortran's global maximum sits within 0.003 of the Q_e pole on 2/1, 3/1, + 4/1 and 6/1. Skipping the pole region moves the threshold to the torque-balance maximum + between Q_e and Q0 (blue dots below), lowering the Fortran-table b_crit by 6–88 %. +- **Layer model.** Fortran solves Park's drift-MHD layer and Julia Fitzpatrick's, so Δ(Q) + differs by a median 7–78 % between Q_i and Q_e, and Julia's b_crit is 82–89 % of the + Fortran-table value on 2/1, 3/1, 4/1 and 6/1. The 5/1 edge surface is dominated by a + zero of the electromagnetic torque in both codes and is not a usable threshold. + +![Fortran vs Julia inner-layer Δ(Q) and torque balance](assets/crf_fortran_benchmark.png) + ## API Reference ### InnerLayer @@ -304,6 +359,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 a8a99f4c7..b3efc474c 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -89,3 +89,11 @@ 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. Qmin/Qmax may be set to +# override the per-surface scan window derived from Q0, Q_e and Q_i. +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..ce620d53d 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -152,6 +152,24 @@ 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. +[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_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 + # Settings (catches accidental config drift) [quantities.slayer_enabled] h5path = "Tearing/enabled" diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 23de76dad..291d2b6a0 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 87781259b..373180e24 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..0cc09daf2 --- /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, torque_balance_window, 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..9cffa549f --- /dev/null +++ b/src/Tearing/CriticalResonantField/TorqueBalance.jl @@ -0,0 +1,113 @@ +# 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))]). +# +# As in Fortran gslayer.f, the viscosity integral is taken as 1/2 (κ̂ = 2/s²) and α +# is kept finite at 1e-2. Q is on Cole's axis: the electron and ion diamagnetic +# frequencies sit at Q = Q_e and Q = Q_i (see `cole_delta`). + +# Finite stand-in for α = S^(-1/3)·(−r_s Δ'_s) ≪ 1; b_crit moves ~4% over 1e-4 to 1e-1. +const TORQUE_BALANCE_ALPHA = 1e-2 +# Half-width around Q_e, Q_i where the layer response is pole-dominated (Fortran SLAYER layfac). +const TORQUE_BALANCE_POLE_WIDTH = 0.02 + +""" + TorqueBalance{M<:InnerLayerModel,P} + +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 + - `Q0` -- normalized natural E×B rotation of the m/n mode, τ_k·n·ω_E + - `P` -- magnetic Prandtl number τ_R/τ_V (the layer's P_φ) + - `lu` -- Lundquist number S + - `kappa_hat` -- Cole's slab-to-tokamak factor κ̂ (Eq. 61) +""" +struct TorqueBalance{M<:InnerLayerModel,P} + model::M + params::P + Q0::Float64 + P::Float64 + lu::Float64 + kappa_hat::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 / (Δ + TORQUE_BALANCE_ALPHA)) + return 2.0 * tb.P * (tb.Q0 - Q) / jxb, Δ +end + +""" + torque_balance_window(Q0, Q_e, Q_i) -> (Qmin, Qmax) + +Q range bracketing the torque-balance branch between the natural rotation `Q0` and the +electron diamagnetic pole `Q_e` (the rule of Fortran gslayer.f). +""" +function torque_balance_window(Q0::Real, Q_e::Real, Q_i::Real) + Q0 > Q_e && return (1.05 * Q_e, 2.0 * Q0) + Qmin = Q0 > 0 ? 0.8 * Q_i : 1.5 * min(Q0, Q_i) + return (Qmin, 0.95 * 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 largest positive interior local +maximum of `bal`, located at `Qpeak = Qs[idx_peak]`. A maximum within +`TORQUE_BALANCE_POLE_WIDTH` (or one grid step, if larger) of the diamagnetic poles `Q_e`, +`Q_i` is discarded as a pole artifact. `Qmin`/`Qmax` +default to `torque_balance_window`. Returns NaN for `Qpeak`, `br_crit` and 0 for +`idx_peak` when no valid maximum is found. +""" +function torque_balance_scan(tb::TorqueBalance; Qmin=nothing, Qmax=nothing, n::Integer=2000) + p = tb.params + wmin, wmax = torque_balance_window(tb.Q0, p.Q_e, p.Q_i) + Qs = range(something(Qmin, wmin), something(Qmax, wmax); length=n) + out = [torque_balance_value(tb, Q) for Q in Qs] + bal = first.(out) + Δs = last.(out) + + w = max(step(Qs), TORQUE_BALANCE_POLE_WIDTH) + near_pole(Q) = abs(Q - p.Q_e) <= w || abs(Q - p.Q_i) <= w + peaks = [i for i in 2:(n-1) if isfinite(bal[i]) && bal[i] > 0 && + bal[i-1] < bal[i] && bal[i] >= bal[i+1] && !near_pole(Qs[i])] + if isempty(peaks) + @warn "CriticalResonantField: no positive torque-balance maximum on " * + "Q ∈ [$(first(Qs)), $(last(Qs))] at the $(p.m)/$(p.n) surface." + return Qs, bal, NaN, NaN, 0, Δs + end + idx = peaks[argmax(bal[peaks])] + 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 1a12ff7e6..370ef04c6 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -5,6 +5,46 @@ # 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. + +## Fields + + - `enabled` -- run the analysis after SLAYER + - `Qmin`, `Qmax` -- override the real-Q scan window; `nothing` derives it per surface from Q0, Q_e and Q_i + - `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 @@ -94,6 +134,11 @@ 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; requires `inner_model = :slayer_fitzpatrick` """ @kwdef struct SLAYERControl enabled::Bool = false @@ -155,6 +200,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, :ggj_shooting, :ggj_galerkin) @@ -189,6 +236,10 @@ function validate(ctrl::SLAYERControl) throw(ArgumentError("SLAYERControl: nre and nim must both be ≥ 2")) ctrl.amr_passes >= 0 || throw(ArgumentError("SLAYERControl: amr_passes must be ≥ 0")) + !ctrl.critical_resonant_field.enabled || ctrl.inner_model === :slayer_fitzpatrick || + throw(ArgumentError("SLAYERControl: CriticalResonantField requires inner_model=:slayer_fitzpatrick")) + ctrl.critical_resonant_field.n >= 3 || + throw(ArgumentError("SLAYERControl: CriticalResonantField n must be ≥ 3")) return ctrl end @@ -224,6 +275,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 b0fc8d2e4..93ddde447 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 # │ (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 @@ -53,10 +54,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["q0"] = crf.q0 + cg["p_phi"] = crf.p_phi + 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" @@ -106,8 +127,14 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)), "PerSurface/tau_A" => (; long_name="Alfvén time τ_A per surface (GGJ layer parameters)", units="s", dims=("surface",)), "PerSurface/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), - "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces, ψ_N-referenced (GGJ path; identical to SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")), - "PerSurface/Delta_prime_matrix_rs" => (; long_name="full complex Δ' matrix as used in the slab-layer matching, converted to the r_s reference length (K^(2α) on the diagonal; the ψ_N-referenced BVP matrix is SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")), + "PerSurface/Delta_prime_matrix" => (; + long_name="full complex Δ' matrix coupling the rational surfaces, ψ_N-referenced (GGJ path; identical to SingularSurfaces/Delta_prime_matrix)", + dims=("surface_row", "surface_col") + ), + "PerSurface/Delta_prime_matrix_rs" => (; + long_name="full complex Δ' matrix as used in the slab-layer matching, converted to the r_s reference length (K^(2α) on the diagonal; the ψ_N-referenced BVP matrix is SingularSurfaces/Delta_prime_matrix)", + dims=("surface_row", "surface_col") + ), "Roots/Q_root" => (; long_name="complex dispersion-root normalized frequency Q (NaN = no root)", dims=("surface",)), "Roots/omega" => (; long_name="mode rotation angular frequency ω = Re(Q)/τ_k of each root", units="rad/s", dims=("surface",)), @@ -120,7 +147,21 @@ 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)", 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)", + units="1", dims=("surface",)), + "CriticalResonantField/q0" => (; long_name="normalized natural E×B rotation Q0 = τ_k·n·ω_E", units="1", dims=("surface",)), + "CriticalResonantField/p_phi" => (; long_name="magnetic Prandtl number used in the torque balance (the layer P_tor)", 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 = [ @@ -147,6 +188,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 60a728869..5bd6861f3 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -3,6 +3,54 @@ # `SLAYERResult` packages the output of a full SLAYER analysis run: # per-surface layer parameters, the extracted tearing eigenvalues, and (if # `control.store_scan`) the full Q-plane scan data for plotting. +# +# `CriticalResonantFieldResult` packages the per-surface critical resonant field +# from the torque-balance analysis (and, if requested, its Q scans). + +""" + 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) + - `br_crit` -- critical normalized resonant field b_r/B_φ (NaN if none) + - `q0` -- normalized natural E×B rotation τ_k·n·ω_E + - `p_phi` -- magnetic Prandtl number used (the layer `P_tor`) + - `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} + q0::Vector{Float64} + p_phi::Vector{Float64} + scan::Vector{CriticalResonantFieldScan} +end + +empty_critical_resonant_field_result() = + CriticalResonantFieldResult(false, Int[], Float64[], Float64[], Float64[], Float64[], CriticalResonantFieldScan[]) """ SLAYERResult @@ -25,9 +73,10 @@ downstream inspection and HDF5 output. `delta_prime_to_rs_reference`, written as `PerSurface/Delta_prime_matrix_rs`); GGJ path: the ψ_N matrix unchanged (written as `PerSurface/Delta_prime_matrix`) - `Q_root` -- tearing eigenvalue(s) in normalized Q - * length `nsurfaces` in `:uncoupled` mode - * length `1` in `:coupled` mode (global eigenvalue normalized by - `params[1].tauk`) + + + length `nsurfaces` in `:uncoupled` mode + + length `1` in `:coupled` mode (global eigenvalue normalized by + `params[1].tauk`) - `omega_Hz`, `gamma_Hz` -- physical rotation frequency / growth rate - `per_surface_extraction` -- `Vector{GrowthRateResult}` of length `nsurfaces` in uncoupled mode (each includes polelines, pole list, @@ -39,6 +88,8 @@ downstream inspection and HDF5 output. plus FKR / visco-resistive sanity scales. Empty when disabled. - `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 @@ -54,16 +105,18 @@ struct SLAYERResult coupled_extraction::Union{Nothing,GrowthRateResult} layer_widths::Vector{LayerWidths} scan_data::Vector{Union{ScanResult,AMRResult}} + critical_resonant_field::CriticalResonantFieldResult end # Empty result (enabled=false path) function empty_slayer_result(control::SLAYERControl) return SLAYERResult(false, control, - SLAYERParameters[], - Float64[], Float64[], - zeros(ComplexF64, 0, 0), - ComplexF64[], Float64[], Float64[], - GrowthRateResult[], nothing, - LayerWidths[], - Union{ScanResult,AMRResult}[]) + SLAYERParameters[], + Float64[], Float64[], + zeros(ComplexF64, 0, 0), + ComplexF64[], Float64[], Float64[], + GrowthRateResult[], nothing, + LayerWidths[], + Union{ScanResult,AMRResult}[], + empty_critical_resonant_field_result()) end diff --git a/src/Tearing/Runner/Runner.jl b/src/Tearing/Runner/Runner.jl index 919088065..be4b74391 100644 --- a/src/Tearing/Runner/Runner.jl +++ b/src/Tearing/Runner/Runner.jl @@ -32,7 +32,8 @@ using ..Utilities using ..Utilities: KineticProfiles using ...Equilibrium: read_kinetic_file, KineticProfileData using ..InnerLayer -using ..InnerLayer: InnerLayerParameters, InnerLayerResponse, solve_inner, +using ..InnerLayer: + InnerLayerParameters, InnerLayerResponse, solve_inner, SLAYERModel, SLAYERParameters, build_slayer_inputs, GGJModel, GGJParameters, LayerWidths, slayer_layer_thickness @@ -44,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") @@ -54,5 +56,8 @@ export SLAYERControl, slayer_control_from_toml, validate export SLAYERResult, empty_slayer_result export run_slayer, 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 2459087cf..17403308f 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -364,7 +364,7 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, return SLAYERResult(true, control, params, rational_psi, rational_q, dp, Q_root, omega_Hz, gamma_Hz, per_surface_extraction, coupled_extraction, - layer_widths, scan_data_list) + layer_widths, scan_data_list, empty_critical_resonant_field_result()) end # --------------------------------------------------------------------- @@ -397,6 +397,36 @@ function ggj_inner_deltas(params::AbstractVector{GGJParameters}, Q::Number; return out end +# --------------------------------------------------------------------- +# Critical resonant field (torque balance) +# --------------------------------------------------------------------- +""" + run_critical_resonant_field(params, rational_psi, profiles, 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 ω_E read from `profiles` at +`rational_psi`. +""" +function run_critical_resonant_field(params::AbstractVector{SLAYERParameters}, rational_psi::AbstractVector{<:Real}, + profiles::KineticProfiles, ctrl::CriticalResonantFieldControl) + ctrl.enabled || return empty_critical_resonant_field_result() + model = SLAYERModel{:fitzpatrick}() + nsurf = length(params) + q_peak, br_crit, q0 = fill(NaN, nsurf), fill(NaN, nsurf), zeros(nsurf) + scan = CriticalResonantFieldScan[] + for (k, p) in enumerate(params) + # Cole Sec. IV: ω0 is the mode frequency in the E×B frame; diamagnetic parts enter via Q_e, Q_i. + q0[k] = p.tauk * p.n * profiles(rational_psi[k]).omega + tb = TorqueBalance(model, p, q0[k], p.P_tor, p.lu, 2.0 / p.sval_r^2) + 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(Qs), bal, Δs)) + end + return CriticalResonantFieldResult(true, Int[p.ising for p in params], q_peak, br_crit, q0, + Float64[p.P_tor for p in params], scan) +end + # --------------------------------------------------------------------- # Full pipeline: equilibrium + ForceFreeStates → parameters → analysis # --------------------------------------------------------------------- @@ -503,5 +533,8 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac rational_psi = Float64[surfaces[p.ising].psifac for p in params] rational_q = Float64[surfaces[p.ising].q for p in params] - return run_slayer_from_inputs(params, dp, control; rational_psi=rational_psi, rational_q=rational_q) + result = run_slayer_from_inputs(params, dp, control; rational_psi=rational_psi, rational_q=rational_q) + control.critical_resonant_field.enabled || return result + crf = run_critical_resonant_field(params, rational_psi, profiles, control.critical_resonant_field) + return SLAYERResult((f === :critical_resonant_field ? crf : getfield(result, f) for f in fieldnames(SLAYERResult))...) end diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 745a30857..d4c225246 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 # @@ -23,12 +25,14 @@ import ..InnerLayer as InnerLayer include("LayerInputs.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 build_ggj_inputs end # module Tearing diff --git a/test/runtests.jl b/test/runtests.jl index e1256e56d..bd96d26d2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,6 +56,7 @@ else include("./runtests_dispersion_amr.jl") include("./runtests_dispersion_polish.jl") include("./runtests_slayer_runner.jl") + include("./runtests_critical_resonant_field.jl") include("./runtests_kinetic.jl") include("./runtests_multiion.jl") include("./runtests_fullruns.jl") diff --git a/test/runtests_critical_resonant_field.jl b/test/runtests_critical_resonant_field.jl new file mode 100644 index 000000000..c8933f411 --- /dev/null +++ b/test/runtests_critical_resonant_field.jl @@ -0,0 +1,115 @@ +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.InnerLayer +using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel + +# Layer stub returning a prescribed torque balance, to test the maximum selection. +struct _CRFStubLayer <: InnerLayerModel end +const _CRF_Q0_STUB = 10.0 +_crf_target(Q) = 0.1 + exp(-(Q - 1)^2 / 0.01) + 2exp(-(Q - 3)^2 / 0.01) +function InnerLayer.solve_inner(::_CRFStubLayer, p, Q::Number) + Qc = -real(Q) # undo the axis mirror cole_delta applies + jxb = 2 * (_CRF_Q0_STUB - Qc) / _crf_target(Qc) + return (; tearing=conj(im / jxb - 1e-2)) +end + +@testset "CriticalResonantField: torque balance" begin + using GeneralizedPerturbedEquilibrium + using GeneralizedPerturbedEquilibrium.InnerLayer + using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel, SLAYERModel, SLAYERParameters + using GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField + using GeneralizedPerturbedEquilibrium.Runner + using GeneralizedPerturbedEquilibrium.Utilities: KineticProfiles + 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) + 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 "Fortran scan window" begin + @test torque_balance_window(5.0, 1.0, -2.0) == (1.05, 10.0) + @test torque_balance_window(0.5, 1.0, -2.0) == (0.8 * -2.0, 0.95) + @test torque_balance_window(-0.5, 1.0, -2.0) == (1.5 * -2.0, 0.95) + end + + @testset "local-maximum selection and pole rejection" begin + p = _mk(; Q_e=50.0, Q_i=-50.0) + tb = TorqueBalance(_CRFStubLayer(), p, _CRF_Q0_STUB, 1.0, p.lu, 1.0) + _, bal, Qpeak, br, idx, _ = torque_balance_scan(tb; Qmin=0.0, Qmax=4.0, n=401) + @test Qpeak ≈ 3.0 atol = 0.011 + @test bal[idx] ≈ _crf_target(Qpeak) rtol = 1e-8 + # A maximum on the electron pole is skipped for the next one. + p_pole = _mk(; Q_e=3.0, Q_i=-50.0) + tb_pole = TorqueBalance(_CRFStubLayer(), p_pole, _CRF_Q0_STUB, 1.0, p.lu, 1.0) + _, _, Qpeak_pole, _, _, _ = torque_balance_scan(tb_pole; Qmin=0.0, Qmax=4.0, n=401) + @test Qpeak_pole ≈ 1.0 atol = 0.011 + # No interior maximum (monotone balance) returns NaN. + _, _, Qnone, brnone, inone, _ = @test_logs (:warn,) torque_balance_scan(tb; Qmin=2.7, Qmax=2.95, n=51) + @test isnan(Qnone) && isnan(brnone) && inone == 0 + end + + @testset "runner: P from P_tor, Q0 = τ_k·n·ω_E" begin + psi = collect(range(0.0, 1.0; length=11)) + ω_E = 4.0e4 + prof = KineticProfiles(; psi=psi, n_e=fill(4e19, 11), T_e=fill(1.7e3, 11), T_i=fill(2e3, 11), + omega=fill(ω_E, 11), omega_e=fill(-1e4, 11), omega_i=fill(1.6e4, 11)) + 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(params, [0.5, 0.7], prof, ctrl) + @test r.enabled + @test r.p_phi == [3.0, 5.0] + @test r.q0 ≈ [p.tauk * p.n * ω_E for p in params] + @test r.rational_index == [1, 2] + @test length(r.scan) == 2 && length(r.scan[1].Q) == 201 + @test !run_critical_resonant_field(params, [0.5, 0.7], prof, CriticalResonantFieldControl()).enabled + + # HDF5 output honours the metadata contract. + c = SLAYERControl(; enabled=true, scan_mode=:brute_force, nre=4, nim=4, critical_resonant_field=ctrl) + base = run_slayer_from_inputs(params, ComplexF64[1.0 0.0; 0.0 1.0], c; rational_psi=[0.5, 0.7], rational_q=[2.0, 1.5]) + res = SLAYERResult((f === :critical_resonant_field ? r : getfield(base, f) for f in fieldnames(SLAYERResult))...) + 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["p_phi"]) == [3.0, 5.0] + @test read(g["br_crit"]) ≈ r.br_crit nans = true + @test haskey(g, "Scan/Surface_2/Delta") + end + end + 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("viscous_input" => 1.0)) + s = slayer_control_from_toml(Dict("enabled" => true, "CriticalResonantField" => Dict("enabled" => true))) + @test s.critical_resonant_field.enabled + @test_throws ArgumentError Runner.validate(SLAYERControl(; inner_model=:ggj_shooting, + critical_resonant_field=CriticalResonantFieldControl(; enabled=true))) + end +end