diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 403fb7ffc..eb7246d4d 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}` (+ `*Winding` with `root_classification = winding`), `Scan/Surface_/`, `CriticalResonantField/` (+ `Scan/Surface_/`) | +| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}` (+ `*Winding` with `root_classification = winding`), `Scan/Surface_/`, `CriticalResonantField/` (+ `Scan/Surface_/`), `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 diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index e67b984ff..71e852ff0 100644 --- a/docs/src/inner_layer.md +++ b/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_/` | 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/`, `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/` | 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 the lab-frame ``n\omega_E`` with `omega_model = "exb"`. + +### 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] +``` diff --git a/regression-harness/cases/diiid_n1_callen.toml b/regression-harness/cases/diiid_n1_callen.toml new file mode 100644 index 000000000..c738c14b6 --- /dev/null +++ b/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" diff --git a/regression-harness/cases/diiid_slayer_n1_mre.toml b/regression-harness/cases/diiid_slayer_n1_mre.toml new file mode 100644 index 000000000..7bd0d9a6d --- /dev/null +++ b/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" diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17c..a69568bb3 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -360,7 +360,9 @@ end """ function tearing_from_deck(inputs::AbstractDict) haskey(inputs, "Tearing") || return nothing - table = inputs["Tearing"] + # A [Tearing] holding only [Tearing.ModifiedRutherford] runs it without a layer solve. + table = _without(inputs["Tearing"], "ModifiedRutherford") + isempty(table) && haskey(inputs["Tearing"], "ModifiedRutherford") && return nothing get(table, "omega_E_kHz", Dict{String,Any}()) isa AbstractDict || throw(ArgumentError("[Tearing] omega_E_kHz must be a table keyed by m/n, e.g. omega_E_kHz = { \"2/1\" = 3.0 }")) fields = _struct_fields(TearingProblem, (:ffs, :critical_resonant_field)) @@ -377,6 +379,21 @@ function tearing_from_deck(inputs::AbstractDict) return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config) end +""" + modified_rutherford_from_deck(inputs) -> Union{Nothing,Dict{Symbol,Any}} + +The `ModifiedRutherfordProblem` keywords of `[Tearing.ModifiedRutherford]`; `nothing` without that section. The key +`preset = "callen"` expands to `CALLEN_PRESET`, under the section's other keys. +""" +function modified_rutherford_from_deck(inputs::AbstractDict) + haskey(inputs, "Tearing") && haskey(inputs["Tearing"], "ModifiedRutherford") || return nothing + table = inputs["Tearing"]["ModifiedRutherford"] + preset = get(table, "preset", "none") + preset in ("none", "callen") || throw(ArgumentError("[Tearing.ModifiedRutherford] preset must be \"none\" or \"callen\", got \"$preset\"")) + kwargs = section_keywords(_without(table, "preset"), "Tearing.ModifiedRutherford", _struct_fields(ModifiedRutherfordProblem, (:ffs, :tearing, :pe))) + return preset == "callen" ? merge(Dict{Symbol,Any}(pairs(ModifiedRutherford.CALLEN_PRESET)), kwargs) : kwargs +end + """ perturbed_equilibrium_from_deck(inputs) -> Union{Nothing,NamedTuple} @@ -471,8 +488,27 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) error("[Match] needs Δ′: use integrator = \"Galerkin\" with [Galerkin] rpec_flag = true, or \"Riccati\" with vac_flag = true") end tearing = tearing_from_deck(inputs) + modified_rutherford = modified_rutherford_from_deck(inputs) + # The ModifiedRutherfordProblem guards, checked before the force-free-states solve rather than after it. + if modified_rutherford !== nothing + has_ffs || error("[Tearing.ModifiedRutherford] needs a [ForceFreeStates] solve") + (; keywords, algorithm) = force_balance + first(keywords[:nn]) == last(keywords[:nn]) || error("[Tearing.ModifiedRutherford] is single-n: set [ForceFreeStates] nn to one toroidal mode number") + # Without a Δ′ of its own (a forward solve), the run adds a Riccati solve, which needs vac_flag. + get(modified_rutherford, :delta_prime_source, :computed) === :cylindrical || + (algorithm isa Galerkin && algorithm.rpec_flag) || get(keywords, :vac_flag, false) || + error( + "[Tearing.ModifiedRutherford] needs Δ′: set [ForceFreeStates] vac_flag = true (Riccati or Forward), use \"Galerkin\" with " * + "[Galerkin] rpec_flag = true, or set delta_prime_source = \"cylindrical\"" + ) + isempty(get(get(inputs, "Equilibrium", Dict{String,Any}()), "kinetic_file", "")) && + error("[Tearing.ModifiedRutherford] needs kinetic profiles: set kinetic_file in [Equilibrium]") + end perturbed_equilibrium = perturbed_equilibrium_from_deck(inputs) forcing = forcing_terms_control(inputs, path) + modified_rutherford === nothing || perturbed_equilibrium === nothing || get(modified_rutherford, :forcing, true) === false || + forcing.forcing_data_format == "coil" || + error("[Tearing.ModifiedRutherford] forcing needs a coil forcing: set [ForcingTerms] forcing_data_format = \"coil\", or forcing = false") if haskey(inputs, "ErrorFields") perturbed_equilibrium === nothing && error("[ErrorFields] needs a [PerturbedEquilibrium] section") forcing.forcing_data_format == "coil" || @@ -483,7 +519,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it") haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) && error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it") - return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, + return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, modified_rutherford, kinetic_forces=haskey(inputs, "KineticForces") ? kinetic_forces_control(inputs) : nothing, perturbed_equilibrium, forcing, error_fields=error_fields_from_deck(inputs, path)) end @@ -578,6 +614,8 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an merge!(inputs, _tearing_sections(filter(p -> p.first != :critical_resonant_field, problem), get(problem, :critical_resonant_field, nothing), model, config)) end + mre = deck.modified_rutherford + mre === nothing || (get!(inputs, "Tearing", Dict{String,Any}())["ModifiedRutherford"] = _deck_table(mre)) deck.kinetic_forces === nothing || merge!(inputs, deck_sections(deck.kinetic_forces)) pe = deck.perturbed_equilibrium pe === nothing || (inputs["PerturbedEquilibrium"] = _deck_table(merge(pe.response, pe.drive))) @@ -626,6 +664,7 @@ deck_sections(m::ForceFreeStates.MatchResult) = merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ")))) deck_sections(kf::KineticForcesResult) = deck_sections(kf.control) deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config) +deck_sections(r::ModifiedRutherfordResult) = Dict{String,Any}("Tearing" => Dict{String,Any}("ModifiedRutherford" => _deck_table(r.problem, (:ffs, :tearing, :pe)))) deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES) function deck_sections(ffs::ForceFreeStatesResult) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index bda2b7c8d..1eaf65168 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -50,6 +50,7 @@ include("Tearing/Tearing.jl") import .Tearing as Tearing export Tearing import .Tearing.CriticalResonantField as CriticalResonantField +import .Tearing.ModifiedRutherford as ModifiedRutherford include("PerturbedEquilibrium/PerturbedEquilibrium.jl") import .PerturbedEquilibrium as PerturbedEquilibrium @@ -74,7 +75,7 @@ include("Rerun.jl") using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult using .ForceFreeStates: MatchProblem using .InnerLayer: GGJModel, SLAYERModel -using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl +using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl, ModifiedRutherfordProblem, ModifiedRutherfordResult using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult using .Equilibrium: PlasmaEquilibrium, KineticProfileData, IonSpecies, attach_kinetic_profiles! using .ForcingTerms: RMPField @@ -140,7 +141,7 @@ equilibrium and writes nothing. The raw forcing file contents are captured in `Input/RawInputs/ForcingTerms/`, and `preloaded_coil_sets` lets the rerun path inject the coil geometry stored in `Input/RawInputs/Coils/`. -Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each +Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing, modified_rutherford)`, each `nothing` when its section is absent. `equil` is the equilibrium the stages used (the solve's, on the two-pass grid), `ffs` the result the later stages read (the matched one when the deck has `[Match]`), and `match` its `MatchResult`. @@ -158,7 +159,7 @@ function main_from_inputs( runtimes = Vector{Pair{String,Float64}}() deck = parse_deck(inputs, path) (; output_file, verbose) = deck.run - local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = nothing + local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = modified_rutherford = nothing equil = _timed_stage(runtimes, "equilibrium", "Equilibrium") do Equilibrium.PlasmaEquilibrium(eq_config, additional_input) @@ -224,12 +225,29 @@ function main_from_inputs( save_to_hdf5(output_file, tearing) end + # The perturbed equilibrium needs a forward solve, so a computed Δ′ for the MRE comes from its own Riccati solve. + has_delta_prime = ffs !== nothing && ffs.delta_prime !== nothing && !isempty(ffs.delta_prime.matrix) + ric = has_delta_prime ? ffs : nothing + if !has_delta_prime && deck.modified_rutherford !== nothing && get(deck.modified_rutherford, :delta_prime_source, :computed) === :computed + ric = _timed_stage(runtimes, "tearing_delta_prime", "Tearing Δ′ (Riccati)") do + solve(ForceBalanceProblem(equil; deck.force_balance.keywords...), Riccati(); verbose) + end + end + + if deck.modified_rutherford !== nothing + modified_rutherford = _timed_stage(runtimes, "modified_rutherford", "Modified Rutherford") do + # With [Tearing], its linear solve gates the roots; with [PerturbedEquilibrium], the coil field forces them. + solve(ModifiedRutherfordProblem(something(ric, ffs); tearing, pe, deck.modified_rutherford...)) + end + save_to_hdf5(output_file, modified_rutherford) + end + total_dt = time() - total_start push!(runtimes, "total" => total_dt) ffs === nothing || _write_runtimes!(output_file, runtimes) @info "\n$_BANNER\n GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER" - return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing) + return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing, modified_rutherford) end """ @@ -305,7 +323,7 @@ export main, write_imas export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields export KineticProfileData, IonSpecies, attach_kinetic_profiles! export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse -export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel +export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel, ModifiedRutherfordProblem, ModifiedRutherfordResult export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index a7205b6cd..91332d2cc 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -34,7 +34,7 @@ model_name(::GGJModel{S}) where {S} = "ggj_$S" model_name(::SLAYERModel{S}) where {S} = "slayer_$S" export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile, layer_response -export layer_eta, layer_rho, layer_omega_E +export layer_plasma, layer_eta, layer_rho, layer_omega_E export GGJ, GGJModel, GGJParameters export mercier_di, mercier_dr, inner_Q, rescale_delta export solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile diff --git a/src/InnerLayer/LayerPlasma.jl b/src/InnerLayer/LayerPlasma.jl index f8dd77105..aab90c9e5 100644 --- a/src/InnerLayer/LayerPlasma.jl +++ b/src/InnerLayer/LayerPlasma.jl @@ -15,6 +15,34 @@ function _require_layer_inputs(who::AbstractString, equil; given...) "(kinetic_file in [Equilibrium], or attach_kinetic_profiles!)")) end +""" + layer_plasma(equil; control=nothing, log=true) -> NamedTuple + +The plasma every inner layer sees, from the equilibrium's kinetic data `equil.kinetic`. + +## Fields + + - `data` - The kinetic data. + - `profiles` - `KineticProfiles` splines of n_e, T_e, T_i and Ω_E on the data's ψ grid. + - `zeff` - Z_eff(ψ) from the data's profile or the n_i/n_e closure, else the scalar `zeff_default`. + - `mu_i`, `z_i` - Main-ion mass in proton masses (fraction-weighted over the species) and charge. + - `mu_rho` - Mean ion mass per electron in proton masses at `zeff`, so ρ = μ m_p n_e (a number or a ψ-callable). + - `resistivity_model`, `lnLambda_form` - The data's resistivity closure and Coulomb-log form. + - `chi_perp_e`, `chi_perp_i`, `chi_tor` - χ(ψ) splines the `SLAYERModel` `control` reads (`nothing` when absent or + unread, so the caller uses `chi_*_default`; all `nothing` without `control`). + - `has_omega_E` - Whether the data carry Ω_E. + +`log = false` skips the note on where Z_eff came from. +""" +function layer_plasma(equil; control=nothing, log::Bool=true) + data = _layer_kinetic_data(equil, "plasma") + loaded = SLAYER._slayer_profiles(data, control) + zeff = log ? SLAYER._layer_zeff(loaded, data) : something(loaded.zeff, data.zeff_default) + return (; data, loaded.profiles, zeff, mu_i=SLAYER._layer_ion_mass(data), z_i=Float64(first(data.ion_species).z), + mu_rho=SLAYER._layer_mass_per_electron(data, zeff), resistivity_model=SLAYER._build_resistivity_model(data.resistivity_model), + data.lnLambda_form, loaded.chi_perp_e, loaded.chi_perp_i, loaded.chi_tor, loaded.zeff_source, loaded.has_omega_E) +end + """ layer_eta(equil, surfaces) -> Vector{Float64} @@ -22,10 +50,9 @@ Resistivity η in Ω·m at each rational surface of `surfaces` (core to edge), f its resistivity closure, the one the SLAYER layer uses. """ function layer_eta(equil, surfaces) - data = _layer_kinetic_data(equil, "resistivity η") - loaded = SLAYER._slayer_profiles(data, nothing) - return SLAYER.layer_resistivity(equil, surfaces, loaded.profiles; zeff=SLAYER._layer_zeff(loaded, data), - resistivity_model=SLAYER._build_resistivity_model(data.resistivity_model), lnLambda_form=data.lnLambda_form) + _layer_kinetic_data(equil, "resistivity η") + plasma = layer_plasma(equil) + return SLAYER.layer_resistivity(equil, surfaces, plasma.profiles; plasma.zeff, plasma.resistivity_model, plasma.lnLambda_form) end """ @@ -36,11 +63,10 @@ Ion mass density ρ = μ·m_p·n_e in kg/m³ at each rational surface of `surfac layer uses. """ function layer_rho(equil, surfaces) - data = _layer_kinetic_data(equil, "mass density ρ") - loaded = SLAYER._slayer_profiles(data, nothing) - mu = SLAYER._layer_mass_per_electron(data, SLAYER._layer_zeff(loaded, data)) + _layer_kinetic_data(equil, "mass density ρ") + plasma = layer_plasma(equil) _at(x, ψ) = x isa Real ? x : x(ψ) - return Float64[_at(mu, s.psifac) * Utilities.M_P * loaded.profiles.n_e(s.psifac) for s in surfaces] + return Float64[_at(plasma.mu_rho, s.psifac) * Utilities.M_P * plasma.profiles.n_e(s.psifac) for s in surfaces] end """ diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl index dbf1ed5b9..9da12ff53 100644 --- a/src/Tearing/HDF5Output.jl +++ b/src/Tearing/HDF5Output.jl @@ -1,5 +1,5 @@ import ..Utilities: save_to_hdf5 -import ..Utilities.HDF5Annotations: replace_group!, write_ragged! +import ..Utilities.HDF5Annotations: write_ragged! # HDF5Output.jl # @@ -19,11 +19,12 @@ using HDF5 """ save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::TearingResult) -Write `result` into a `Tearing/` subgroup of `parent`, replacing any existing one. +Write `result` into a `Tearing/` subgroup of `parent`, replacing its previous contents but keeping the +`ModifiedRutherford/` group that its result writes there. """ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::TearingResult) - g = replace_group!(parent, "Tearing") + g = _reset_tearing_group!(parent) g["enabled"] = Int(result.enabled) # Which inner-layer model produced PerSurface/: the SLAYER and GGJ branches write # disjoint field sets, so readers must not have to infer it from the schema. @@ -182,6 +183,22 @@ const TEARING_RAGGED_H5_ANNOTATIONS = [ "offsets" => (; long_name="ragged-array offsets: row k spans offsets[k]+1:offsets[k+1]") ] +# Groups under Tearing/ that other results own, so a TearingResult write keeps them. +const _TEARING_FOREIGN_GROUPS = ("ModifiedRutherford",) + +# Empty Tearing/ of everything a TearingResult writes, or create it. +function _reset_tearing_group!(parent) + haskey(parent, "Tearing") || return create_group(parent, "Tearing") + g = parent["Tearing"] + for k in keys(g) + k in _TEARING_FOREIGN_GROUPS || HDF5.delete_object(g, k) + end + for a in keys(attrs(g)) + HDF5.delete_attribute(g, a) + end + return g +end + const TEARING_RAGGED_WINDING_H5_ANNOTATIONS = [ "flat" => (; long_name="winding number of the dispersion residual around each entry of the matching root list (NaN: circle unresolved)"), "offsets" => (; long_name="ragged-array offsets: row k spans offsets[k]+1:offsets[k+1]") diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl index 3375731a4..67327e4da 100644 --- a/src/Tearing/LayerOverlap.jl +++ b/src/Tearing/LayerOverlap.jl @@ -403,14 +403,12 @@ Eq. (100) is anchored to. Caveats are logged as warnings and notes as info. function layer_overlap_scan(equil, model::SLAYERModel, ns) data = equil.kinetic data === nothing && error("psilim_from_layer_overlap needs kinetic profiles: set kinetic_file in [Equilibrium] or attach_kinetic_profiles!") - loaded = _slayer_profiles(data, model) - chi = _slayer_chi_inputs(loaded, model, data; warn=false) - zeff = loaded.zeff === nothing ? data.zeff_default : loaded.zeff + plasma = layer_plasma(equil; control=model, log=false) + chi = _slayer_chi_inputs(plasma, model, data; warn=false) scans = [ - resistive_layer_overlap(equil, loaded.profiles; n_tor=n, rs_method=:flux, chi..., tau_E=data.tau_E, - P_perp_model=model.P_perp_model, mu_i=_layer_ion_mass(data), mu_rho=_layer_mass_per_electron(data, zeff), zeff=zeff, - theta=model.theta_sample, - resistivity_model=_build_resistivity_model(data.resistivity_model), lnLambda_form=data.lnLambda_form) + resistive_layer_overlap(equil, plasma.profiles; n_tor=n, rs_method=:flux, chi..., tau_E=data.tau_E, + P_perp_model=model.P_perp_model, plasma.mu_i, plasma.mu_rho, plasma.zeff, theta=model.theta_sample, + plasma.resistivity_model, plasma.lnLambda_form) for n in ns ] for (n, scan) in zip(ns, scans) diff --git a/src/Tearing/ModifiedRutherford/Geometry.jl b/src/Tearing/ModifiedRutherford/Geometry.jl new file mode 100644 index 000000000..09374f851 --- /dev/null +++ b/src/Tearing/ModifiedRutherford/Geometry.jl @@ -0,0 +1,159 @@ +# Geometry.jl +# +# Flux-surface averages the modified Rutherford terms need beyond `ResistGeometry`: the Hegna 1999 R̄², +# poloidal-field averages, the mean minor radius and the Sauter 2002 trapped fraction, as in Fortran RDCON `mercier.f`. +# +# θ-integrands, each weighted by `jac / v1` before integration: +# 1: B²/|∇ψ|² 2: B² 3: |∇ψ|² 4: R² v1/jac 5: B_pol 6: B_pol² 7: r (minor radius) + +""" + MRESurfaceGeometry + +Per-surface geometry of the modified Rutherford terms. Averages ⟨·⟩ are normalized flux-surface averages (⟨1⟩ = 1) +and `ψ_N` is normalized poloidal flux, as in [`ResistGeometry`](@ref GeneralizedPerturbedEquilibrium.ForceFreeStates.ResistGeometry). + +## Fields + + - `q`, `q1` - Safety factor and dq/dψ_N. + - `p1` - μ₀ dp/dψ_N. + - `v1` - dV/dψ_N [m³]. + - `psio` - Poloidal flux per radian at the boundary [Wb]. + - `avg_bsq_over_dpsisq` - ⟨B²/‖∇ψ_N‖²⟩. + - `avg_bsq` - ⟨B²⟩ [T²]. + - `avg_dpsisq` - ⟨‖∇ψ_N‖²⟩. + - `overbar_Rsq` - Hegna 1999 R̄² = ⟨R² V′/𝒥⟩ [m²]. + - `avg_bpol`, `avg_bpolsq` - ⟨B_pol⟩ [T] and ⟨B_pol²⟩ [T²]. + - `r_minor` - ⟨r⟩, mean distance from the magnetic axis [m]. + - `R_major`, `eps` - (R_max + R_min)/2 [m] and (R_max − R_min)/(R_max + R_min). + - `delta` - Triangularity, mean of the upper and lower values. + - `f_trap` - Sauter 2002 trapped fraction with the triangularity-corrected ε_eff. + - `Wc_prefac` - Schlutt & Hegna 2012 w_d⁴ prefactor in ψ_N⁴, before χ⊥/χ∥ and 1/m². +""" +struct MRESurfaceGeometry + q::Float64 + q1::Float64 + p1::Float64 + v1::Float64 + psio::Float64 + avg_bsq_over_dpsisq::Float64 + avg_bsq::Float64 + avg_dpsisq::Float64 + overbar_Rsq::Float64 + avg_bpol::Float64 + avg_bpolsq::Float64 + r_minor::Float64 + R_major::Float64 + eps::Float64 + delta::Float64 + f_trap::Float64 + Wc_prefac::Float64 +end + +""" + sauter2002_trapped_fraction(eps, delta) -> Float64 + +Trapped fraction of Sauter 2002 (as Fortran RDCON uses it): ε_eff = max(0, 0.67(1 − 1.4δ|δ|)ε) and +`f_t = 1 − √((1 − ε)/(1 + ε)) (1 − ε_eff)/(1 + 2√ε_eff)`, capped at 1. +""" +function sauter2002_trapped_fraction(eps::Real, delta::Real) + eps_eff = max(0.0, 0.67 * (1 - 1.4 * delta * abs(delta)) * eps) + return min(1 - sqrt((1 - eps) / (1 + eps)) * (1 - eps_eff) / (1 + 2 * sqrt(eps_eff)), 1.0) +end + +""" + bootstrap_coefficient(g::MRESurfaceGeometry, zeff) -> (D_nc, mu_ratio) + +Bootstrap drive coefficient of Hegna 1999 as Fortran RDCON evaluates it, +`D_nc = −(q/q′)(μ₀p′/⟨B²⟩) R̄² ⟨B²/‖∇ψ_N‖²⟩ ψ_o⁻² · μ_e/(ν_e + μ_e)`, with the banana-limit viscosity ratio +`μ_e/(ν_e + μ_e) = f_t a/((1 − f_t) + f_t a)`, `a = 1 + 0.533/Z_eff`. +""" +function bootstrap_coefficient(g::MRESurfaceGeometry, zeff::Real) + a = 1 + 0.533 / zeff + mu_ratio = g.f_trap * a / ((1 - g.f_trap) + g.f_trap * a) + return -g.q * (g.p1 / (g.q1 * g.avg_bsq)) * g.overbar_Rsq * g.avg_bsq_over_dpsisq / g.psio^2 * mu_ratio, mu_ratio +end + +""" + mre_geometry(equil, psifac, q1) -> MRESurfaceGeometry + +Flux-surface geometry of the modified Rutherford terms at the interior surface `psifac` with dq/dψ_N `q1`. +""" +function mre_geometry(equil, psifac::Real, q1::Real) + profiles = equil.profiles + chi1 = 2π * equil.psio + psi = Float64(psifac) + twopif = profiles.F_spline(psi) + p1 = profiles.P_deriv(psi) + v1 = profiles.dVdpsi_spline(psi) + q = profiles.q_spline(psi) + + ys = equil.rzphi_ys + ff = zeros(Float64, length(ys), 7) + R_max, R_min, Z_max, Z_min = -Inf, Inf, -Inf, Inf + R_at_Zmax = R_at_Zmin = equil.ro + for (i, theta) in enumerate(ys) + rfac = sqrt(equil.rzphi_rsquared((psi, theta))) + jac = equil.rzphi_jac((psi, theta)) + fy1 = FastInterpolations.deriv_view(equil.rzphi_rsquared, (0, 1))((psi, theta)) + fy2 = FastInterpolations.deriv_view(equil.rzphi_offset, (0, 1))((psi, theta)) + fy3 = FastInterpolations.deriv_view(equil.rzphi_nu, (0, 1))((psi, theta)) + eta = 2π * (theta + equil.rzphi_offset((psi, theta))) + r = equil.ro + rfac * cos(eta) + z = equil.zo + rfac * sin(eta) + + v21 = fy1 / (2 * rfac * jac) + v22 = (1 + fy2) * 2π * rfac / jac + v23 = fy3 * r / jac + v33 = 2π * r / jac + bsq = chi1^2 * (v21^2 + v22^2 + (v23 + q * v33)^2) + dpsisq = (2π * r)^2 * (v21^2 + v22^2) + bpolsq = max(bsq - (twopif / (2π * r))^2, 0.0) + + R_max, R_min = max(R_max, r), min(R_min, r) + z > Z_max && ((Z_max, R_at_Zmax) = (z, r)) + z < Z_min && ((Z_min, R_at_Zmin) = (z, r)) + ff[i, :] .= (bsq / dpsisq, bsq, dpsisq, r^2 * v1 / jac, sqrt(bpolsq), bpolsq, rfac) + @views ff[i, :] .*= jac / v1 + end + @views ff[end, :] .= ff[1, :] + avg = FastInterpolations.integrate(cubic_interp(ys, Series(ff); bc=PeriodicBC())) + + R_major = 0.5 * (R_max + R_min) + a_minor = 0.5 * (R_max - R_min) + eps = a_minor / R_major + delta = 0.5 * ((R_major - R_at_Zmax) + (R_major - R_at_Zmin)) / a_minor + # Schlutt & Hegna 2012 ⟨𝒥B²⟩⟨𝒥|∇ψ|²⟩/(dι/dψ)² in toroidal flux, converted to ψ_N⁴. + Wc_prefac = (v1 * avg[2] / (q * equil.psio)) * v1 * avg[3] * (q * equil.psio) * (q^6 / (q1 / equil.psio)^2) / (q * equil.psio)^4 + return MRESurfaceGeometry(q, q1, p1, v1, equil.psio, avg[1], avg[2], avg[3], avg[4], avg[5], avg[6], avg[7], + R_major, eps, delta, sauter2002_trapped_fraction(eps, delta), Wc_prefac) +end + +""" + current_density_profile(equil) -> Function + +Dimensionless toroidal current density `J_z(ψ_N) = (R₀²/B₀) ⟨μ₀ j_φ/R⟩` of the cylindrical saturation theory, from the +Grad–Shafranov flux-surface identity `⟨μ₀ j_φ/R⟩ = μ₀ dp/dψ + F dF/dψ ⟨R⁻²⟩` (ψ the poloidal flux per radian), signed +positive in the core. It reduces to `R₀μ₀j_z/B_z` in the large-aspect-ratio cylinder. +""" +function current_density_profile(equil) + profiles = equil.profiles + R0, B0, psio = Float64(equil.ro), Float64(equil.params.b0), Float64(equil.psio) + ys = equil.rzphi_ys + function avg_inv_Rsq(psi) + ff = zeros(Float64, length(ys), 2) + for (i, theta) in enumerate(ys) + rfac = sqrt(equil.rzphi_rsquared((psi, theta))) + jac = equil.rzphi_jac((psi, theta)) + r = equil.ro + rfac * cos(2π * (theta + equil.rzphi_offset((psi, theta)))) + ff[i, 1], ff[i, 2] = jac / r^2, jac + end + @views ff[end, :] .= ff[1, :] + avg = FastInterpolations.integrate(cubic_interp(ys, Series(ff); bc=PeriodicBC())) + return avg[1] / avg[2] + end + raw(psi) = (profiles.P_deriv(psi) + profiles.F_spline(psi) * profiles.F_deriv(psi) / (4π^2) * avg_inv_Rsq(psi)) / psio + # Sign from a core surface inside the grid, so the core current is positive. + sgn = sign(raw(max(0.05, Float64(first(equil.rzphi_xs))))) + sgn == 0 && throw(ArgumentError("current_density_profile: zero current density at the core surface")) + return psi -> sgn * R0^2 / B0 * raw(Float64(psi)) +end diff --git a/src/Tearing/ModifiedRutherford/HDF5Output.jl b/src/Tearing/ModifiedRutherford/HDF5Output.jl new file mode 100644 index 000000000..276d036a1 --- /dev/null +++ b/src/Tearing/ModifiedRutherford/HDF5Output.jl @@ -0,0 +1,129 @@ +# HDF5Output.jl +# +# Writer for a `ModifiedRutherfordResult`, under the caller's parent group: +# +# Tearing/ +# └── ModifiedRutherford/ +# ├── w -- width grid (ψ_N) +# ├── w_m_star, b_m_star -- smallest seed threshold and its field over the surfaces +# ├── b_crit_star -- smallest critical field of the forcing over the surfaces +# ├── PerSurface/ -- surface identity, indices, term inputs, roots and flags +# └── Terms/ -- F(w), each term and w_d(w), surface × width + +""" + save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::ModifiedRutherfordResult) + +Write `result` into `Tearing/ModifiedRutherford/` of `parent`, replacing any existing one and keeping the rest of +`Tearing/`, and record its deck section. +""" +function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::ModifiedRutherfordResult) + tearing = haskey(parent, "Tearing") ? parent["Tearing"] : create_group(parent, "Tearing") + g = replace_group!(tearing, "ModifiedRutherford") + # One dataset per annotated field; flags as 0/1 integers, readable without an HDF5 boolean type. + for (key, _) in MRE_H5_ANNOTATIONS + value = getfield(result, Symbol(basename(key))) + g[key] = value isa AbstractVector{Bool} ? Int8.(value) : value + end + Utilities.HDF5Annotations.annotate!(g, MRE_H5_ANNOTATIONS) + Utilities.record_inputs!(parent, result) + return g +end + +# Julia's surface × width matrices land in HDF5 as (width, surface). +const _MRE_TERM_DIMS = ("w", "surface") +const _MRE_BRANCH_DIMS = ("branch", "surface") + +const MRE_H5_ANNOTATIONS = [ + "w" => (; long_name="full island width grid", units="psi_N", dims=("w",)), + "w_m_star" => (; long_name="smallest seed-threshold full width over the surfaces, NaN if none", units="psi_N"), + "b_m_star" => (; long_name="smallest b_m over the surfaces, NaN if none", units="T"), + "b_crit_star" => (; long_name="smallest b_crit over the surfaces, NaN if none", units="T"), + "PerSurface/rational_psi" => (; long_name="normalized poloidal flux of the rational surface", units="psi_N", dims=("surface",)), + "PerSurface/rational_q" => (; long_name="safety factor at the surface", dims=("surface",)), + "PerSurface/m" => (; long_name="resonant poloidal mode number", dims=("surface",)), + "PerSurface/n" => (; long_name="resonant toroidal mode number", dims=("surface",)), + "PerSurface/delta_prime" => (; + long_name="real part of the diagonal Delta' (PEST-3 convention; cylindrical -2m referred to psi_N with delta_prime_source = cylindrical), times delta_prime_factor", + dims=("surface",) + ), + "PerSurface/delta_prime_imag" => (; long_name="imaginary part of the diagonal Delta'", dims=("surface",)), + "PerSurface/D_I" => (; long_name="Mercier index D_I = E + F + H - 1/4", dims=("surface",)), + "PerSurface/D_R" => (; long_name="resistive interchange index D_R = E + F + H^2", dims=("surface",)), + "PerSurface/H" => (; long_name="Glasser-Greene-Johnson coefficient H", dims=("surface",)), + "PerSurface/D_nc" => (; long_name="bootstrap drive coefficient D_nc (Hegna 1999)", dims=("surface",)), + "PerSurface/f_trap" => (; long_name="Sauter 2002 trapped fraction", dims=("surface",)), + "PerSurface/mu_ratio" => (; long_name="banana-limit viscosity ratio mu_e/(nu_e + mu_e)", dims=("surface",)), + "PerSurface/zeff" => (; long_name="effective charge Z_eff", dims=("surface",)), + "PerSurface/Wc_prefac" => (; long_name="transport-width prefactor W_c (Schlutt & Hegna 2012), before chi_perp/chi_par and 1/m^2", units="psi_N^4", + dims=("surface",)), + "PerSurface/chi_perp" => (; long_name="perpendicular electron heat diffusivity used in w_d", units="m^2/s", dims=("surface",)), + "PerSurface/v_te" => (; long_name="electron thermal speed sqrt(2T_e/m_e)", units="m/s", dims=("surface",)), + "PerSurface/R0" => (; long_name="major radius (R_max + R_min)/2 of the surface, used in lambda_par", units="m", dims=("surface",)), + "PerSurface/omega_star_i" => (; long_name="ion diamagnetic frequency from the density gradient", units="rad/s", dims=("surface",)), + "PerSurface/omega_star_e" => (; long_name="electron diamagnetic frequency from the density gradient", units="rad/s", dims=("surface",)), + "PerSurface/omega_hat" => (; long_name="island frequency normalized to the ion diamagnetic frequency", dims=("surface",)), + "PerSurface/omega_hat_pi" => (; long_name="omega_*pi/omega_*i = 1 + eta_i", dims=("surface",)), + "PerSurface/beta_p" => (; long_name="poloidal ion beta 2 mu0 n_i T_i/", dims=("surface",)), + "PerSurface/rho_theta" => (; long_name="poloidal ion Larmor radius", units="m", dims=("surface",)), + "PerSurface/w_ib" => (; long_name="ion banana width", units="psi_N", dims=("surface",)), + "PerSurface/eps" => (; long_name="inverse aspect ratio of the surface", dims=("surface",)), + "PerSurface/gamma" => (; long_name="linear tearing growth rate of the gate, NaN without one", units="1/s", dims=("surface",)), + "PerSurface/delta_r" => (; long_name="resistive layer width of the gate, NaN without one", units="psi_N", dims=("surface",)), + "PerSurface/w_lo" => (; long_name="lower end of the root search max(delta_r, w_min)", units="psi_N", dims=("surface",)), + "PerSurface/D_c" => (; long_name="SLAYER toroidal critical offset used by curvature_model = lutjens (r_s reference), NaN otherwise", + dims=("surface",)), + "PerSurface/W_d_layer" => (; long_name="transport width sqrt(2) pi^1.5 (-D_R)/D_c in the (r - r_s)/r_s reference, NaN otherwise", + dims=("surface",)), + "PerSurface/alpha_s" => + (; long_name="saturation current-gradient coefficient alpha_s = -r q J_z'/s (Fitzpatrick 2023 Eq. 9.16), J_z = (R0^2/B0) ", dims=("surface",)), + "PerSurface/beta_s" => (; long_name="saturation current-curvature coefficient beta_s = -r^2 q J_z''/s", dims=("surface",)), + "PerSurface/C_sat" => (; long_name="saturation coefficient 0.8 alpha_s^2 - 0.27 beta_s - 0.09 alpha_s (Fitzpatrick 2023 Eq. 9.15)", + dims=("surface",)), + "PerSurface/w_m" => (; long_name="seed threshold: smallest F crossing from negative to positive, NaN if none", units="psi_N", dims=("surface",)), + "PerSurface/w_sat" => (; long_name="saturated width: first falling F crossing above the island start (w_lo or w_m), NaN if none", units="psi_N", + dims=("surface",)), + "PerSurface/fully_stabilized" => (; long_name="1 if F < 0 on the whole search range", dims=("surface",)), + "PerSurface/unstable_any_seed" => (; long_name="1 if F > 0 at w_lo", dims=("surface",)), + "PerSurface/unsaturated" => (; long_name="1 if F > 0 at w_max: no saturation in the search range", dims=("surface",)), + "PerSurface/outside_plasma" => (; long_name="1 if psi +- w_sat/2 leaves the equilibrium psi grid", dims=("surface",)), + "PerSurface/gate_disagrees" => (; long_name="1 if F(w_lo) > 0 with gamma < 0, or F(w_lo) < 0 with gamma > 0", dims=("surface",)), + "PerSurface/interchange_unstable" => (; long_name="1 if D_R > 0 (excluded)", dims=("surface",)), + "PerSurface/linear_unstable" => (; long_name="1 if the gate's growth rate is positive (no seed threshold, w_sat kept)", dims=("surface",)), + "PerSurface/below_w_ib" => (; long_name="1 if w_m < w_ib", dims=("surface",)), + "PerSurface/large_delta_prime_w" => (; long_name="1 if |Delta'(w)| w >= 1 at w_m or w_sat, Delta'(w) the classical term", dims=("surface",)), + "PerSurface/saturation_unconverged" => (; long_name="1 if C_sat moves more than 20% between half and double the fit window, or the window leaves the grid", + dims=("surface",)), + "PerSurface/w_sat_all" => (; long_name="every saturated width (falling F crossing) in increasing order, NaN-padded; w_sat is the first", + units="psi_N", dims=_MRE_BRANCH_DIMS), + "PerSurface/w_seed_all" => (; long_name="seed threshold reaching each w_sat_all entry (rising crossing below it), NaN if F > 0 from w_lo", + units="psi_N", dims=_MRE_BRANCH_DIMS), + "PerSurface/w_seed_unsaturated" => (; long_name="seed above which the island grows past w_max (F > 0 there), NaN otherwise", units="psi_N", + dims=("surface",)), + "PerSurface/b_m" => (; long_name="resonant area-weighted field whose vacuum island has full width w_m, B = pi |n q'| chi1 w^2/(8 A)", + units="T", dims=("surface",)), + "PerSurface/b_seed_all" => (; long_name="resonant area-weighted field whose vacuum island is as wide as each w_seed_all entry", units="T", + dims=_MRE_BRANCH_DIMS), + "PerSurface/b_seed_unsaturated" => (; long_name="resonant area-weighted field whose vacuum island is as wide as w_seed_unsaturated", + units="T", dims=("surface",)), + "PerSurface/n_sat" => (; long_name="number of saturated widths", dims=("surface",)), + "PerSurface/multiple_saturated" => (; long_name="1 if more than one saturated width (see w_sat_all)", dims=("surface",)), + "PerSurface/b_vac" => (; long_name="vacuum (coil) resonant area-weighted field of the perturbed equilibrium, NaN without one", units="T", dims=("surface",)), + "PerSurface/b_res" => (; long_name="total resonant area-weighted field of the perturbed equilibrium, NaN without one", units="T", dims=("surface",)), + "PerSurface/delta_rmp" => (; long_name="Callen Delta_RMP = |Delta|/(2 pi R0 q) B0/|B_res| in 1/m of rho, NaN without a perturbed equilibrium", + units="1/m", dims=("surface",)), + "PerSurface/w_v" => (; long_name="vacuum island full width at the applied field, NaN without a perturbed equilibrium", units="psi_N", dims=("surface",)), + "PerSurface/b_crit" => (; long_name="vacuum resonant field at which a forced window F > 0 opens below the unforced seed threshold", units="T", + dims=("surface",)), + "PerSurface/b_trigger" => (; long_name="vacuum resonant field at which the forced window reaches the unforced seed threshold", units="T", + dims=("surface",)), + "PerSurface/w_crit" => (; long_name="full width where the forced window opens at b_crit, NaN if none", units="psi_N", dims=("surface",)), + "PerSurface/complex_delta_prime" => (; long_name="1 if |Im Delta'| > 1e-3 |Re Delta'|", dims=("surface",)), + "Terms/F" => (; long_name="modified Rutherford F(w), the sum of the switched-on terms", dims=_MRE_TERM_DIMS), + "Terms/classical" => (; long_name="classical term Delta'(w)", dims=_MRE_TERM_DIMS), + "Terms/saturation" => (; long_name="finite-width saturation term Delta_sat(w)", dims=_MRE_TERM_DIMS), + "Terms/bootstrap" => (; long_name="bootstrap term Delta_nc(w)", dims=_MRE_TERM_DIMS), + "Terms/curvature" => (; long_name="curvature term Delta_GGJ(w) of curvature_model", dims=_MRE_TERM_DIMS), + "Terms/polarization" => (; long_name="ion polarization term Delta_pol(w), including pol_sign", dims=_MRE_TERM_DIMS), + "Terms/forcing" => (; long_name="coil forcing term Delta_RMP (drho/dpsi_N) (w_v/w)^2, 0 without a perturbed equilibrium", dims=_MRE_TERM_DIMS), + "Terms/w_d" => (; long_name="transport threshold width w_d(w)", units="psi_N", dims=_MRE_TERM_DIMS) +] diff --git a/src/Tearing/ModifiedRutherford/ModifiedRutherford.jl b/src/Tearing/ModifiedRutherford/ModifiedRutherford.jl new file mode 100644 index 000000000..561879741 --- /dev/null +++ b/src/Tearing/ModifiedRutherford/ModifiedRutherford.jl @@ -0,0 +1,32 @@ +# The steady modified Rutherford equation: switchable terms and their sum F(w) per rational surface, from the +# force-balance Δ′, the resistive-layer geometry and the equilibrium's kinetic data. + +module ModifiedRutherford + +using HDF5 +import CommonSolve +import FastInterpolations +using FastInterpolations: cubic_interp, Series, PeriodicBC, DerivOp + +import ...Utilities +import ...Equilibrium +import ...Utilities: save_to_hdf5 +using ...Utilities.HDF5Annotations: replace_group! +using ...Utilities.PhysicalConstants: M_E, M_P, E_CHG, MU_0 +using ...InnerLayer: layer_plasma +using ...InnerLayer.SLAYER: radial_label, build_slayer_inputs +using Roots: find_zero, Brent +using ..ResonantForcing: resonant_forcing, toroidal_flux_radius, callen_polarization_width, forcing_thresholds + +include("Geometry.jl") +include("Terms.jl") +include("Problem.jl") +include("HDF5Output.jl") + +export MRESurfaceGeometry, mre_geometry, current_density_profile, sauter2002_trapped_fraction, bootstrap_coefficient +export classical_term, saturation_coefficient, saturation_term, bootstrap_term, curvature_term, lutjens_curvature_term, layer_transport_width, polarization_term, + reduced_polarization_term, + transport_width +export mre_roots, vacuum_island_field, ModifiedRutherfordProblem, ModifiedRutherfordResult, CALLEN_PRESET + +end # module ModifiedRutherford diff --git a/src/Tearing/ModifiedRutherford/Problem.jl b/src/Tearing/ModifiedRutherford/Problem.jl new file mode 100644 index 000000000..289221d52 --- /dev/null +++ b/src/Tearing/ModifiedRutherford/Problem.jl @@ -0,0 +1,520 @@ +# Problem.jl +# +# `solve(ModifiedRutherfordProblem(ffs; ...))`: the terms of the steady modified Rutherford equation and their sum F(w) +# on a log-spaced grid of full island widths, per rational surface, and its roots w_m and w_sat. + +""" + ModifiedRutherfordProblem(ffs; kwargs...) + +The steady modified Rutherford equation (MRE) on every rational surface of a finished force-balance solve, + + F(w) = s_Δ Δ′(w) + s_sat Δ_sat(w) + s_nc Δ_nc(w) + s_GGJ Δ_GGJ(w) + s_pol Δ_pol(w) + s_RMP Δ_RMP(w), + +with `w` the full island width in normalized poloidal flux ψ_N. Each term is evaluated on the grid whether or not its +switch is on; the switches select which terms `F` sums. Δ′ is the real part of the diagonal of `ffs.delta_prime` +(`delta_prime_source = :computed`) or the cylindrical vacuum value r_s Δ′ = −2m referred to ψ_N as SLAYER refers Δ′ to +r_s (`:cylindrical`, which needs no Δ′ solve); +the Mercier, resistive-interchange and Pfirsch-Schlüter indices D_I, D_R = E + F + H², H come from the surfaces' +`ResistGeometry`, and the plasma (densities, temperatures, Z_eff, χ⊥,e, E×B rotation, ion species) is the +equilibrium's kinetic data, as every tearing layer reads it. A deck states it as `[Tearing.ModifiedRutherford]`. + +The roots of F are refined between grid points from `w_lo = max(δ_r, w_min)`: the seed threshold `w_m` is the +smallest crossing from F < 0 to F > 0 and the saturated width `w_sat` the first crossing from F > 0 to F < 0 above where +the island starts (`w_lo` when F(w_lo) > 0, else `w_m`). Every saturated width and the seed threshold that reaches it are +kept as well (`w_sat_all`, `w_seed_all`), so a surface with several saturated solutions shows all of them. Surfaces that +are resistive-interchange unstable (D_R > 0) get no roots; linearly tearing unstable ones (γ > 0 in `tearing`) grow +from noise, so they get no seed threshold but keep `w_sat`. + +With a coil-driven perturbed equilibrium `pe`, the applied field forces the island: Δ_RMP(w) = Δ_RMP (dρ/dψ_N) (w_v/w)² +(Callen's Eq. 17, [`resonant_forcing`](@ref)), with w_v the vacuum island width and the island locked in phase with the +field. It scales with the vacuum resonant field B_vac, so the roots and branches are those of the forced F, and the +critical fields `b_crit` and `b_trigger` follow from the unforced F ([`forcing_thresholds`](@ref)). + +`CALLEN_PRESET` (deck `preset = "callen"`) reduces F to Callen's steady state Δ_c + Δ_RMP (w_v/w)² − w_pol²/w³ in his +units (Fortran GPEC 1.6.0's w_pol, Δ_RMP and Δ_c), so its roots are the full widths of Callen's cubic +Δ_c w³ + Δ_RMP w_v² w − w_pol² = 0 and its `b_crit` is his critical field B_crit = |B_vac| (w_v,crit/w_v)². + +## Fields + + - `ffs` - Outer solution, single n, with a Δ′ matrix (Riccati with `vac_flag=true`, or Galerkin with `rpec_flag=true`) + unless `delta_prime_source = :cylindrical`. + - `tearing` - Linear tearing solve on the same Δ′ (the gate and the layer width δ_r), or `nothing`; set by the deck + when `[Tearing]` is present. + - `pe` - Single-n, coil-driven perturbed equilibrium on the same equilibrium and n, or `nothing`; set by the deck when + `[PerturbedEquilibrium]` is present. Its rows are matched to the surfaces by ψ. + - `classical`, `saturation`, `bootstrap`, `curvature`, `polarization`, `forcing` - Term switches s_Δ, s_sat, s_nc, s_GGJ, + s_pol, s_RMP. `saturation` is the finite-width correction to Δ′ ([`saturation_term`](@ref)), off by default; + `forcing` acts only with `pe`. + - `saturation_window` - Half-width in ψ_N of the quadratic fit that gives J_z′ and J_z″ in [`saturation_coefficient`](@ref). + - `delta_prime_source` - `:computed` (Δ′_mn of `ffs`) or `:cylindrical` (r_s Δ′_mn = −2m, a reference case). + - `delta_prime_model` - `:mercier_power` (Benjamin Eq. 10), `:constant` (Δ′(w) = Δ′_mn), or `:callen` (constant, + Callen's Δ_c = K^(2α) Δ′_mn/ρ, or −2m/ρ with `:cylindrical`, in ψ_N; K = ρ dψ_N/dρ on the toroidal-flux radius ρ). + - `delta_prime_factor`, `D_nc_factor`, `rho_theta_factor` - Multipliers on Δ′_mn, D_nc and the poloidal ion Larmor + radius ρ_θi (which also sets w_ib), for scans. + - `curvature_model` - `:lutjens` (default; Lütjens et al. 2001, [`lutjens_curvature_term`](@ref), whose small-w limit is + the SLAYER toroidal critical offset D_c of the linear gate) or `:benjamin` (Benjamin Eq. 12, [`curvature_term`](@ref), + for comparison with the tearing-physics suite). `:lutjens` takes W_d and D_c from SLAYER's toroidal layer parameters + built on each surface with this problem's χ⊥. + - `k1` - Bootstrap coefficient k₁, and curvature coefficient of `curvature_model = :benjamin`. + - `c_wd` - Prefactor of the transport threshold width w_d. + - `chi_perp_factor` - Multiplier on χ⊥,e (the kinetic `chi_e` profile, else its `chi_e_default`). + - `pol_model` - `:full` (Benjamin Eq. C.1), `:reduced` (−c_pol w_pol²/w³, w_pol = C_pol w_ib) or `:callen` + (−c_pol w_pol²/w³ with Callen's w_pol, [`callen_polarization_width`](@ref)). + - `c_pol`, `C_pol` - Polarization coefficients. + - `pol_sign` - ±1, multiplies the polarization term. + - `omega_model` - Island frequency ω in the E_r = 0 frame: `:electron_diamagnetic`, ω = ω_*e(1 + η_e/4) (Wilson + 1996), or `:exb`, the lab-frame ω = n Ω_E from the kinetic `omega_E`. + - `omega_factor` - Multiplier on the ω of `omega_model`; 0 is a stalled island. + - `omega_ratio` - ω/ω_*i on every surface, replacing `omega_model`; `nothing` (default) uses the model. + - `w_min`, `w_max`, `nw` - Log-spaced width grid in ψ_N. +""" +@kwdef struct ModifiedRutherfordProblem{R,T,S} + ffs::R + tearing::T = nothing + pe::S = nothing + classical::Bool = true + bootstrap::Bool = true + curvature::Bool = true + polarization::Bool = false + saturation::Bool = false + forcing::Bool = true + saturation_window::Float64 = 0.02 + curvature_model::Symbol = :lutjens + delta_prime_source::Symbol = :computed + delta_prime_model::Symbol = :mercier_power + delta_prime_factor::Float64 = 1.0 + D_nc_factor::Float64 = 1.0 + rho_theta_factor::Float64 = 1.0 + k1::Float64 = 1.7 + c_wd::Float64 = 1.0 + chi_perp_factor::Float64 = 1.0 + pol_model::Symbol = :full + c_pol::Float64 = 1.0 + C_pol::Float64 = 1.9 + pol_sign::Float64 = 1.0 + omega_model::Symbol = :electron_diamagnetic + omega_factor::Float64 = 1.0 + omega_ratio::Union{Nothing,Float64} = nothing + w_min::Float64 = 1e-5 + w_max::Float64 = 0.2 + nw::Int = 200 +end + +ModifiedRutherfordProblem(ffs; kwargs...) = ModifiedRutherfordProblem(; ffs=ffs, kwargs...) + +""" + CALLEN_PRESET + +`ModifiedRutherfordProblem` keywords that reduce F to Callen's steady state Δ_c + Δ_RMP (w_v/w)² − w_pol²/w³: a +constant classical term in Callen's units, his polarization width, the forcing, and no other term. Keywords given with +the preset override it, e.g. `ModifiedRutherfordProblem(ffs; pe, CALLEN_PRESET..., delta_prime_source=:cylindrical)`. +""" +const CALLEN_PRESET = (; classical=true, bootstrap=false, curvature=false, saturation=false, polarization=true, forcing=true, + delta_prime_model=:callen, pol_model=:callen, c_pol=1.0, pol_sign=1.0) + +""" + ModifiedRutherfordResult + +Output of `solve(::ModifiedRutherfordProblem)`, built by keyword. Per-surface vectors run core to edge; matrices are +surface × width. + +## Fields + + - `problem` - The producing problem. + - `w` - Width grid in ψ_N. + - `rational_psi`, `rational_q`, `m`, `n` - Surface identity. + - `delta_prime`, `delta_prime_imag`, `D_I`, `D_R`, `H` - Outer-region and layer indices (Δ′ after `delta_prime_factor`). + - `D_nc`, `f_trap`, `mu_ratio`, `zeff` - Bootstrap coefficient, Sauter 2002 trapped fraction, μ_e/(ν_e + μ_e), Z_eff. + - `Wc_prefac`, `chi_perp`, `v_te`, `R0` - Inputs of w_d. + - `omega_star_i`, `omega_star_e`, `omega_hat`, `omega_hat_pi` - Ion and electron diamagnetic frequencies [rad/s] from + the density gradients, ω/ω_*i and ω_*pi/ω_*i. + - `beta_p`, `rho_theta`, `w_ib`, `eps` - Polarization inputs: poloidal ion beta, poloidal ion Larmor radius [m], + banana width in ψ_N, inverse aspect ratio. + - `gamma`, `delta_r` - Linear growth rate [1/s] and resistive layer width in ψ_N from `tearing` (`NaN` without it). + - `w_lo` - Lower end of the root search. + - `D_c`, `W_d_layer` - SLAYER toroidal critical offset (r_s reference) and the transport width it implies, in the + (r − r_s)/r_s reference, used by `curvature_model = :lutjens` (`NaN` for `:benjamin`). + - `alpha_s`, `beta_s`, `C_sat` - Current-gradient coefficients of the saturation term and its coefficient. + - `w_m`, `w_sat` - Seed threshold and saturated full widths in ψ_N; `NaN` when absent or excluded, and `w_m` `NaN` when γ > 0. + - `w_sat_all`, `w_seed_all` - Every saturated width per surface (surface × branch, increasing, `NaN`-padded) and the + seed threshold that reaches each: seeds between `w_seed_all[s, k]` and the next threshold settle at `w_sat_all[s, k]`; + `w_seed_all` is `NaN` on a branch F > 0 reaches from `w_lo`. `w_sat` is the first column. + - `w_seed_unsaturated` - Seed above which the island grows past `w_max`, when F > 0 there; `NaN` otherwise. + - `n_sat` - Number of saturated widths; `multiple_saturated` when more than one. + - `w_m_star` - Smallest `w_m` over the surfaces. + - `b_m`, `b_seed_all`, `b_seed_unsaturated`, `b_m_star` - Resonant area-weighted field [T] whose vacuum island is as wide as + `w_m`, each `w_seed_all` and `w_seed_unsaturated` ([`vacuum_island_field`](@ref)), comparable with the perturbed + equilibrium's resonant field and with error-field thresholds; `b_m_star` is the smallest `b_m` over the surfaces (not + necessarily at the `w_m_star` surface). + - `b_vac`, `b_res`, `delta_rmp`, `w_v` - Vacuum and total resonant fields [T], Δ_RMP [1/m] and the full vacuum island + width in ψ_N from `pe` ([`resonant_forcing`](@ref)); `NaN` without it. + - `b_crit`, `b_trigger`, `w_crit` - Critical vacuum resonant fields [T] of the forcing, from the unforced F + ([`forcing_thresholds`](@ref)): a forced window opens at `b_crit` (at full width `w_crit`) and reaches the intrinsic + seed threshold at `b_trigger`; `NaN` without `pe`. `b_crit_star` is the smallest `b_crit` over the surfaces. + - `fully_stabilized` - F < 0 on the whole search range. + - `unstable_any_seed` - F > 0 at `w_lo`. + - `unsaturated` - F > 0 at `w_max`: an island past its last rising root does not saturate in the search range. + - `outside_plasma` - ψ ± `w_sat`/2 leaves the equilibrium's ψ grid. + - `gate_disagrees` - F(w_lo) > 0 with γ < 0, or F(w_lo) < 0 with γ > 0. + - `interchange_unstable`, `linear_unstable` - D_R > 0 (excluded), or γ > 0 (no seed threshold; `w_sat` kept). + - `below_w_ib` - `w_m` below the ion banana width, where the polarization ordering fails. + - `large_delta_prime_w` - |Δ′(w)| w ≥ 1 at `w_m` or `w_sat`, with Δ′(w) the classical term. + - `complex_delta_prime` - |Im Δ′_mn| > 10⁻³ |Re Δ′_mn|. + - `saturation_unconverged` - `C_sat` not converged in the fit window (see [`saturation_coefficient`](@ref)). + - `F`, `classical`, `saturation`, `bootstrap`, `curvature`, `polarization`, `forcing`, `w_d` - F(w), each term, and + w_d(w); `forcing` is 0 without `pe`. +""" +@kwdef struct ModifiedRutherfordResult{P} + problem::P + w::Vector{Float64} + rational_psi::Vector{Float64} + rational_q::Vector{Float64} + m::Vector{Int} + n::Vector{Int} + delta_prime::Vector{Float64} + delta_prime_imag::Vector{Float64} + D_I::Vector{Float64} + D_R::Vector{Float64} + H::Vector{Float64} + D_nc::Vector{Float64} + f_trap::Vector{Float64} + mu_ratio::Vector{Float64} + zeff::Vector{Float64} + Wc_prefac::Vector{Float64} + chi_perp::Vector{Float64} + v_te::Vector{Float64} + R0::Vector{Float64} + omega_star_i::Vector{Float64} + omega_star_e::Vector{Float64} + omega_hat::Vector{Float64} + omega_hat_pi::Vector{Float64} + beta_p::Vector{Float64} + rho_theta::Vector{Float64} + w_ib::Vector{Float64} + eps::Vector{Float64} + gamma::Vector{Float64} + delta_r::Vector{Float64} + w_lo::Vector{Float64} + D_c::Vector{Float64} + W_d_layer::Vector{Float64} + alpha_s::Vector{Float64} + beta_s::Vector{Float64} + C_sat::Vector{Float64} + w_m::Vector{Float64} + w_sat::Vector{Float64} + w_sat_all::Matrix{Float64} + w_seed_all::Matrix{Float64} + w_seed_unsaturated::Vector{Float64} + n_sat::Vector{Int} + b_m::Vector{Float64} + b_seed_all::Matrix{Float64} + b_seed_unsaturated::Vector{Float64} + w_m_star::Float64 + b_m_star::Float64 + b_vac::Vector{Float64} + b_res::Vector{Float64} + delta_rmp::Vector{Float64} + w_v::Vector{Float64} + b_crit::Vector{Float64} + b_trigger::Vector{Float64} + w_crit::Vector{Float64} + b_crit_star::Float64 + fully_stabilized::Vector{Bool} + unstable_any_seed::Vector{Bool} + unsaturated::Vector{Bool} + multiple_saturated::Vector{Bool} + outside_plasma::Vector{Bool} + gate_disagrees::Vector{Bool} + interchange_unstable::Vector{Bool} + linear_unstable::Vector{Bool} + below_w_ib::Vector{Bool} + large_delta_prime_w::Vector{Bool} + complex_delta_prime::Vector{Bool} + saturation_unconverged::Vector{Bool} + F::Matrix{Float64} + classical::Matrix{Float64} + saturation::Matrix{Float64} + bootstrap::Matrix{Float64} + curvature::Matrix{Float64} + polarization::Matrix{Float64} + forcing::Matrix{Float64} + w_d::Matrix{Float64} +end + +function validate(prob::ModifiedRutherfordProblem) + ffs = prob.ffs + ffs.npert == 1 || throw(ArgumentError("ModifiedRutherfordProblem is single-n; the solve spans n = $(ffs.nlow):$(ffs.nhigh)")) + prob.delta_prime_source in (:computed, :cylindrical) || throw(ArgumentError("delta_prime_source must be :computed or :cylindrical")) + if prob.delta_prime_source === :computed + dp = ffs.delta_prime + (dp === nothing || isempty(dp.matrix)) && + throw( + ArgumentError( + "ModifiedRutherfordProblem needs Δ′: solve with Riccati() and vac_flag=true, Galerkin(; rpec_flag=true), " * + "or set delta_prime_source = :cylindrical" + ) + ) + size(dp.matrix, 1) == length(ffs.surfaces) || + throw(ArgumentError("ModifiedRutherfordProblem: the Δ′ matrix is $(size(dp.matrix, 1))×$(size(dp.matrix, 2)) for $(length(ffs.surfaces)) surfaces")) + end + ffs.equil.kinetic === nothing && throw(ArgumentError("ModifiedRutherfordProblem needs kinetic profiles; set kinetic_file in [Equilibrium]")) + prob.delta_prime_model in (:mercier_power, :constant, :callen) || throw(ArgumentError("delta_prime_model must be :mercier_power, :constant or :callen")) + prob.curvature_model in (:lutjens, :benjamin) || throw(ArgumentError("curvature_model must be :lutjens or :benjamin")) + prob.pol_model in (:full, :reduced, :callen) || throw(ArgumentError("pol_model must be :full, :reduced or :callen")) + pe = prob.pe + if pe !== nothing + pffs = pe.response.ffs + pffs.npert == 1 && pffs.nlow == ffs.nlow || throw(ArgumentError("ModifiedRutherfordProblem: pe must be single-n with the n of ffs")) + end + prob.omega_model in (:electron_diamagnetic, :exb) || throw(ArgumentError("omega_model must be :electron_diamagnetic or :exb")) + abs(prob.pol_sign) == 1 || throw(ArgumentError("pol_sign must be +1 or -1")) + 0 < prob.w_min < prob.w_max && prob.nw >= 2 || throw(ArgumentError("the width grid needs 0 < w_min < w_max and nw ≥ 2")) + prob.k1 > 0 && prob.c_wd > 0 && prob.chi_perp_factor > 0 && prob.rho_theta_factor > 0 || + throw(ArgumentError("k1, c_wd, chi_perp_factor and rho_theta_factor must be positive")) + return nothing +end + +# Value and ψ-derivative of a kinetic profile column at ψ. +function _profile_at(psi_xs, values, psi) + itp = cubic_interp(psi_xs, collect(Float64, values)) + return Float64(itp(psi)), Float64(itp(psi; deriv=DerivOp(1))) +end + + +""" + mre_roots(F, w_lo, w_grid) -> (; w_m, w_sat, w_sat_all, w_seed_all, w_seed_unsaturated, fully_stabilized, unstable_any_seed, unsaturated) + +Roots of `F` on `[w_lo, maximum(w_grid)]`, bracketed on `w_lo` and the grid points above it and refined by Brent's +method: `w_m` is the smallest crossing from negative to positive, and `w_sat` the first crossing from positive to +negative above where an island starts, `w_lo` when F(w_lo) > 0, else `w_m` (`NaN` when absent). `unsaturated` is +F > 0 at the top of the grid. + +Every saturated branch is returned too: `w_sat_all` holds each falling crossing in increasing order, and `w_seed_all[k]` +the rising crossing just below `w_sat_all[k]` (`NaN` when F > 0 from `w_lo` up to it). An island seeded between +`w_seed_all[k]` and the next seed threshold settles at `w_sat_all[k]`. `w_seed_unsaturated` is the rising crossing +above the last falling one when F > 0 at the top of the grid: seeds above it grow past the grid (`NaN` otherwise). +""" +function mre_roots(F, w_lo::Real, w_grid::AbstractVector) + nodes = vcat(Float64(w_lo), filter(>(w_lo), w_grid)) + vals = F.(nodes) + rising, falling = Float64[], Float64[] + for j in 1:(length(nodes)-1) + a, b = vals[j], vals[j+1] + isfinite(a) && isfinite(b) && a != 0 && (a * b < 0 || b == 0) || continue + push!(a < 0 ? rising : falling, find_zero(F, (nodes[j], nodes[j+1]), Brent())) + end + w_m = isempty(rising) ? NaN : first(rising) + start = vals[1] > 0 ? nodes[1] : w_m + above = filter(>(start), falling) + w_sat = isempty(above) ? NaN : first(above) + # F is continuous, so rising and falling crossings alternate; each falling one is reached from the rising one below it. + seed_below(x) = (b = filter(<(x), rising); isempty(b) ? NaN : last(b)) + w_seed_all = seed_below.(above) + top = isempty(above) ? filter(>(-Inf), rising) : filter(>(last(above)), rising) + w_seed_unsaturated = vals[end] > 0 && !isempty(top) ? last(top) : NaN + return (; w_m, w_sat, w_sat_all=above, w_seed_all, w_seed_unsaturated, fully_stabilized=all(<(0), vals), unstable_any_seed=vals[1] > 0, + unsaturated=vals[end] > 0) +end + +""" + vacuum_island_field(equil, psi, n, q1, w) -> Float64 + +Resonant area-weighted field [T] whose vacuum island at the surface `psi` (toroidal mode number `n`, dq/dψ_N `q1`) has +full width `w` in ψ_N: the inverse of the perturbed equilibrium's island width, B = π |n q1| χ₁ w²/(8A), with χ₁ = 2πψ₀ +and A the flux-surface area. +""" +function vacuum_island_field(equil, psi::Real, n::Integer, q1::Real, w::Real) + area = Equilibrium.flux_surface_area(equil, Float64(psi), length(equil.rzphi_ys) - 1) + return π * abs(n * q1) * 2π * equil.psio * w^2 / (8 * area) +end + +# Per-surface scalars, flags and width rows of the result, in field order. +const _MRE_SCALARS = (:rational_psi, :rational_q, :m, :n, :delta_prime, :delta_prime_imag, :D_I, :D_R, :H, :D_nc, :f_trap, :mu_ratio, :zeff, + :Wc_prefac, :chi_perp, :v_te, :R0, :omega_star_i, :omega_star_e, :omega_hat, :omega_hat_pi, :beta_p, :rho_theta, :w_ib, :eps, :gamma, + :delta_r, :w_lo, :D_c, :W_d_layer, :alpha_s, :beta_s, :C_sat, :w_m, :w_sat, :w_seed_unsaturated, :n_sat, :b_m, :b_seed_unsaturated, :b_vac, :b_res, :delta_rmp, :w_v, + :b_crit, :b_trigger, :w_crit) +const _MRE_FLAGS = (:fully_stabilized, :unstable_any_seed, :unsaturated, :multiple_saturated, :outside_plasma, :gate_disagrees, :interchange_unstable, + :linear_unstable, :below_w_ib, :large_delta_prime_w, :complex_delta_prime, :saturation_unconverged) +const _MRE_TERMS = (:F, :classical, :saturation, :bootstrap, :curvature, :polarization, :forcing, :w_d) + +# SLAYER's toroidal critical offset D_c and the W_d it implies at one surface, with χ⊥ for every χ input. +function _layer_critical_offset(equil, s, plasma, chi_perp) + p = only( + build_slayer_inputs(equil, [s], plasma.profiles; plasma.mu_i, plasma.mu_rho, plasma.zeff, chi_perp_e=chi_perp, chi_perp_i=chi_perp, + chi_perp_Dc=chi_perp, chi_tor=chi_perp, tau_E=plasma.data.tau_E, P_perp_model=:chi_perp_e, dc_type=:toroidal, + plasma.resistivity_model, plasma.lnLambda_form) + ) + return p.dc_tmp, layer_transport_width(p.dr_val, p.dc_tmp) +end + +# The linear tearing row at ψ in `tearing`, or `nothing`. +_tearing_row(tearing, psi) = findfirst(p -> isapprox(p, psi; atol=1e-8), tearing.rational_psi) + +# What every surface shares: the plasma, the radial label, the current profile and the gate. +function _mre_context(prob::ModifiedRutherfordProblem) + equil = prob.ffs.equil + plasma = layer_plasma(equil; log=false) + data = plasma.data + usable(v) = v !== nothing && !all(iszero, v) + usable(data.chi_e) || data.chi_e_default !== nothing || + throw(ArgumentError("ModifiedRutherfordProblem needs χ⊥,e: a chi_e profile or chi_e_default in the kinetic data")) + gate = prob.tearing + # Metres map to ψ_N with the toroidal-flux radial label's dr/dψ, the label SLAYER's layer widths use. + r_at, dr_dpsi = radial_label(equil) + return (; equil, plasma, data, psi_xs=collect(Float64, data.psi), usable, gate, r_at, dr_dpsi, J_z=current_density_profile(equil), + forcing=prob.pe === nothing ? nothing : resonant_forcing(prob.pe), + psi_lim=(first(equil.profiles.xs), last(equil.profiles.xs)), + chi_gate=gate === nothing || !hasproperty(gate.model, :chi_perp_Dc) ? nothing : gate.model.chi_perp_Dc) +end + +# Geometry, indices, plasma and term inputs at the k-th rational surface. +function _surface_inputs(prob::ModifiedRutherfordProblem, ctx, k, s) + (; equil, plasma, data, psi_xs, usable, r_at, dr_dpsi) = ctx + rg = s.restype + rg === nothing && throw(ArgumentError("ModifiedRutherfordProblem: the surface at ψ = $(s.psifac) has no ResistGeometry")) + psi = s.psifac + m, n = s.m[1], s.n[1] + g = mre_geometry(equil, psi, s.q1) + D_I = rg.E + rg.F + rg.H - 0.25 + D_R = rg.E + rg.F + rg.H^2 + D_I < 0 || @warn "ModifiedRutherfordProblem: D_I = $D_I ≥ 0 at the $m/$n surface; its terms are NaN" + + # Plasma at the surface; the ion density falls back to n_e where the data have no n_i. + prof = plasma.profiles + n_e, dn_e = Float64(prof.n_e(psi)), Float64(prof.n_e(psi; deriv=DerivOp(1))) + T_e, dT_e = Float64(prof.T_e(psi)), Float64(prof.T_e(psi; deriv=DerivOp(1))) + T_i, dT_i = Float64(prof.T_i(psi)), Float64(prof.T_i(psi; deriv=DerivOp(1))) + n_i, dn_i = usable(data.n_i) ? _profile_at(psi_xs, data.n_i, psi) : (n_e, dn_e) + zeff = plasma.zeff isa Real ? Float64(plasma.zeff) : Float64(plasma.zeff(psi)) + chi_perp = prob.chi_perp_factor * (usable(data.chi_e) ? first(_profile_at(psi_xs, data.chi_e, psi)) : data.chi_e_default) + ctx.chi_gate === nothing || isapprox(ctx.chi_gate, chi_perp; rtol=1e-6) || + @warn "ModifiedRutherfordProblem: the gate's chi_perp_Dc = $(ctx.chi_gate) differs from the w_d χ⊥ = $chi_perp at the $m/$n surface" + D_nc, mu_ratio = bootstrap_coefficient(g, zeff) + + # Polarization inputs (Benjamin App. C); lengths in metres map to ψ_N with dr/dψ. + m_i, z_i = plasma.mu_i * M_P, plasma.z_i + v_te = sqrt(2 * E_CHG * T_e / M_E) + v_ti = sqrt(2 * E_CHG * T_i / m_i) + rho_theta = prob.rho_theta_factor * m_i * v_ti / (g.avg_bpol * E_CHG) + rho_theta_psi = rho_theta / abs(dr_dpsi(psi)) + # Diamagnetic frequencies from the density gradients, in SLAYER's sign convention (charge −1 for electrons). + omega_star_i = -n * T_i * (dn_i / n_i) / (z_i * g.psio) + omega_star_e = n * T_e * (dn_e / n_e) / g.psio + omega_model = + prob.omega_model === :electron_diamagnetic ? omega_star_e * (1 + 0.25 * (dT_e / T_e) / (dn_e / n_e)) : + n * (plasma.has_omega_E ? Float64(prof.omega_E(psi)) : 0.0) + + # Cylindrical: r_s Δ′ = −2m, referred back to ψ_N by K^(−2α), K = r_s dψ_N/dr, α = √(−D_I) (as delta_prime_to_rs_reference). + K = r_at(psi) / abs(dr_dpsi(psi)) + dp_complex = prob.delta_prime_source === :computed ? prob.ffs.delta_prime.matrix[k, k] : complex(-2 * m * K^(-2 * sqrt(max(-D_I, 0.0)))) + # Saturation coefficients from the equilibrium current in the same radial label. + sc = saturation_coefficient(ctx.J_z, r_at, x -> abs(dr_dpsi(x)), x -> Float64(equil.profiles.q_spline(x)), + x -> Float64(equil.profiles.q_deriv(x)), psi; window=prob.saturation_window, ctx.psi_lim) + # SLAYER's toroidal critical offset at this surface with this χ⊥, so the curvature term meets the gate at small w. + D_c, W_d_layer = prob.curvature_model === :lutjens ? _layer_critical_offset(equil, s, plasma, chi_perp) : (NaN, NaN) + # Callen's toroidal-flux radius (without the F/B₀R₀ factor of the layer label), his Δ_c and w_pol, in ψ_N. + rho, drho = toroidal_flux_radius(equil, psi, s.q) + delta_c = prob.delta_prime_source === :computed ? (rho / drho)^(2 * sqrt(max(-D_I, 0.0))) * real(dp_complex) / rho : -2 * m / rho + w_pol_callen = callen_polarization_width(equil, plasma, psi, s.q)[2] / drho + delta_prime = prob.delta_prime_factor * (prob.delta_prime_model === :callen ? delta_c * drho : real(dp_complex)) + # The coil forcing at this surface: Δ_RMP (dρ/dψ_N) w_v² is the coefficient of 1/w², proportional to B_vac. + rf = ctx.forcing + row = rf === nothing ? nothing : findfirst(p -> isapprox(p, psi; atol=1e-8), rf.rational_psi) + b_vac, b_res, delta_rmp, w_v = row === nothing ? (NaN, NaN, NaN, NaN) : (rf.b_vac[row], rf.b_res[row], rf.delta_rmp[row], 2 * rf.hw_v[row]) + forcing_amp = row === nothing ? 0.0 : delta_rmp * drho * w_v^2 + return (; rational_psi=psi, rational_q=s.q, m, n, g, b_per_w2=vacuum_island_field(equil, psi, n, s.q1, 1.0), D_I, D_R, H=rg.H, D_nc=prob.D_nc_factor * D_nc, mu_ratio, zeff, + chi_perp, v_te, + omega_star_i, omega_star_e, omega_hat=something(prob.omega_ratio, prob.omega_factor * omega_model / omega_star_i), + omega_hat_pi=1 + (dT_i / T_i) / (dn_i / n_i), beta_p=2 * MU_0 * n_i * E_CHG * T_i / g.avg_bpolsq, rho_theta, rho_theta_psi, + w_ib=sqrt(g.eps) * rho_theta_psi, gradient_ratio=(dn_i / n_i) / (g.q1 / g.q), K, dp_complex, + delta_prime, sc, D_c, W_d_layer, w_pol_callen, b_vac, b_res, delta_rmp, w_v, forcing_amp, has_forcing=row !== nothing) +end + +# Each term and their switched sum F at the full width w, from a surface's inputs. +function _mre_terms(prob::ModifiedRutherfordProblem, si, w) + (; g, D_I, D_R, n, m) = si + model = prob.delta_prime_model === :callen ? :constant : prob.delta_prime_model + wd = transport_width(w; prob.c_wd, si.chi_perp, si.v_te, R0=g.R_major, g.q, g.q1, n, m, g.Wc_prefac) + cl = D_I < 0 ? classical_term(si.delta_prime, D_I, w; model) : NaN + st = D_I < 0 ? saturation_term(si.sc.C_sat, si.K, D_I, w; model) : NaN + bs = bootstrap_term(si.D_nc, wd, w; prob.k1) + cv = + D_I >= 0 ? NaN : + prob.curvature_model === :lutjens ? lutjens_curvature_term(D_R, D_I, si.K, si.W_d_layer, w; model) : + curvature_term(D_R, si.H, D_I, wd, w; prob.k1) + pl = + prob.pol_sign * ( + prob.pol_model === :full ? + polarization_term(w; prob.c_pol, g.eps, si.rho_theta_psi, si.beta_p, si.gradient_ratio, si.omega_hat, si.omega_hat_pi, si.w_ib) : + prob.pol_model === :reduced ? reduced_polarization_term(w; prob.c_pol, prob.C_pol, si.w_ib) : + reduced_polarization_term(w; prob.c_pol, C_pol=1.0, w_ib=si.w_pol_callen) + ) + fc = si.forcing_amp / w^2 + F0 = prob.classical * cl + prob.saturation * st + prob.bootstrap * bs + prob.curvature * cv + prob.polarization * pl + return (; F=F0 + prob.forcing * fc, F0, classical=cl, saturation=st, bootstrap=bs, curvature=cv, polarization=pl, forcing=fc, w_d=wd) +end + +# The linear tearing growth rate and resistive layer width in ψ_N from `tearing`, and the bottom of the root search. +function _gate(prob::ModifiedRutherfordProblem, ctx, psi) + gate = ctx.gate + row = gate === nothing ? nothing : _tearing_row(gate, psi) + gamma = row === nothing ? NaN : gate.gamma[row] + delta_r = row === nothing || isempty(gate.layer_widths) ? NaN : gate.layer_widths[row].delta_s_m / abs(ctx.dr_dpsi(psi)) + return (; gamma, delta_r, w_lo=isnan(delta_r) ? prob.w_min : max(delta_r, prob.w_min)) +end + +# Roots and flags at one surface. +function _surface_roots(prob::ModifiedRutherfordProblem, ctx, si, gt, w) + F(x) = _mre_terms(prob, si, x).F + r = mre_roots(F, gt.w_lo, w) + interchange_unstable = si.D_R > 0 + linear_unstable = gt.gamma > 0 + excluded = interchange_unstable || si.D_I >= 0 + # A linearly unstable surface grows from noise, so it has no seed threshold; the MRE still sets its saturation. + w_m, w_sat = excluded ? (NaN, NaN) : (linear_unstable ? NaN : r.w_m, r.w_sat) + w_sat_all, w_seed_all = excluded ? (Float64[], Float64[]) : (r.w_sat_all, r.w_seed_all) + w_seed_unsaturated = excluded ? NaN : r.w_seed_unsaturated + psi = si.rational_psi + # Resonant field of the vacuum island as wide as each seed threshold; B ∝ w². + b(x) = si.b_per_w2 * x^2 + # Critical fields of the forcing from the unforced F; the forcing is (forcing_amp/b_vac) B_vac/w². + F0(x) = _mre_terms(prob, si, x).F0 + th = + excluded || !si.has_forcing ? (; b_crit=NaN, b_trigger=NaN, w_crit=NaN) : + forcing_thresholds(F0, si.forcing_amp / si.b_vac, gt.w_lo, w, mre_roots(F0, gt.w_lo, w).w_m) + return (; th..., w_m, w_sat, w_sat_all, w_seed_all, w_seed_unsaturated, n_sat=length(w_sat_all), b_m=b(w_m), b_seed_all=b.(w_seed_all), + b_seed_unsaturated=b(w_seed_unsaturated), + multiple_saturated=length(w_sat_all) > 1, fully_stabilized=!excluded && r.fully_stabilized, unstable_any_seed=!excluded && r.unstable_any_seed, + unsaturated=!excluded && r.unsaturated, outside_plasma=isfinite(w_sat) && (psi - w_sat / 2 < ctx.psi_lim[1] || psi + w_sat / 2 > ctx.psi_lim[2]), + gate_disagrees=!excluded && isfinite(gt.gamma) && (r.unstable_any_seed ? gt.gamma < 0 : gt.gamma > 0), interchange_unstable, + linear_unstable, below_w_ib=w_m < si.w_ib, + # Δ′(w) w at the root, dimensionless at any D_I. + large_delta_prime_w=any(x -> isfinite(x) && abs(_mre_terms(prob, si, x).classical) * x >= 1, (w_m, w_sat)), + complex_delta_prime=abs(imag(si.dp_complex)) > 1e-3 * abs(real(si.dp_complex)), saturation_unconverged=!si.sc.converged) +end + +function CommonSolve.solve(prob::ModifiedRutherfordProblem) + validate(prob) + ctx = _mre_context(prob) + w = exp.(range(log(prob.w_min), log(prob.w_max); length=prob.nw)) + rows = map(enumerate(prob.ffs.surfaces)) do (k, s) + si = _surface_inputs(prob, ctx, k, s) + gt = _gate(prob, ctx, si.rational_psi) + terms = [_mre_terms(prob, si, wj) for wj in w] + extra = (; si.g.f_trap, si.g.Wc_prefac, R0=si.g.R_major, si.g.eps, delta_prime_imag=imag(si.dp_complex), si.sc.alpha_s, si.sc.beta_s, si.sc.C_sat) + rowterms = NamedTuple{_MRE_TERMS}(Tuple(getfield.(terms, key) for key in _MRE_TERMS)) + merge(si, extra, gt, _surface_roots(prob, ctx, si, gt, w), rowterms) + end + column(key) = [getfield(r, key) for r in rows] + matrix(key) = isempty(rows) ? zeros(0, length(w)) : permutedims(reduce(hcat, column(key))) + # Branches padded with NaN to the most found on any surface. + nbranch = maximum((r.n_sat for r in rows); init=1) + branches(key) = [k <= length(getfield(r, key)) ? getfield(r, key)[k] : NaN for r in rows, k in 1:nbranch] + w_m = column(:w_m) + finite_min(v) = all(isnan, v) ? NaN : minimum(filter(!isnan, v)) + return ModifiedRutherfordResult(; problem=prob, w=collect(w), (key => column(key) for key in _MRE_SCALARS)..., b_crit_star=finite_min(column(:b_crit)), + w_sat_all=branches(:w_sat_all), w_seed_all=branches(:w_seed_all), b_seed_all=branches(:b_seed_all), + b_m_star=all(isnan, column(:b_m)) ? NaN : minimum(filter(!isnan, column(:b_m))), + w_m_star=all(isnan, w_m) ? NaN : minimum(filter(!isnan, w_m)), (key => Vector{Bool}(column(key)) for key in _MRE_FLAGS)..., + (key => matrix(key) for key in _MRE_TERMS)...) +end diff --git a/src/Tearing/ModifiedRutherford/Terms.jl b/src/Tearing/ModifiedRutherford/Terms.jl new file mode 100644 index 000000000..582e7e25f --- /dev/null +++ b/src/Tearing/ModifiedRutherford/Terms.jl @@ -0,0 +1,134 @@ +# Terms.jl +# +# The terms of the steady modified Rutherford equation F(w) as pure functions of the full island width w in ψ_N +# (Benjamin et al. 2026 Eqs. 7, 10-14 and C.1; Leuthold et al. 2026 Eq. 4.6). Surface quantities are arguments, +# so each term can be checked alone. + +""" + classical_term(delta_prime, D_I, w; model=:mercier_power) -> Float64 + +Classical drive with the finite-pressure correction, `Δ′(w) = Δ′_mn (w/2)^(−2α₋) (−4D_I)^(1/2)` with +`α₋ = ½ − (−D_I)^(1/2)` (Benjamin Eq. 10); it reduces to `Δ′_mn` at `D_I = −¼`. `model = :constant` returns `Δ′_mn`. +""" +function classical_term(delta_prime::Real, D_I::Real, w::Real; model::Symbol=:mercier_power) + model === :constant && return Float64(delta_prime) + alpha_minus = 0.5 - sqrt(-D_I) + return delta_prime * (w / 2)^(-2 * alpha_minus) * sqrt(-4 * D_I) +end + +""" + bootstrap_term(D_nc, w_d, w; k1=1.7) -> Float64 + +Bootstrap drive with the transport threshold, `k₁ D_nc w / (w² + (k₁/0.5) w_d²)` (Benjamin Eq. 11): tends to +`k₁ D_nc / w` for `w ≫ w_d` and to `0.5 D_nc w / w_d²` for `w ≪ w_d`. +""" +bootstrap_term(D_nc::Real, w_d::Real, w::Real; k1::Real=1.7) = k1 * D_nc * w / (w^2 + (k1 / 0.5) * w_d^2) + +""" + curvature_term(D_R, H, D_I, w_d, w; k1=1.7) -> Float64 + +Glasser-Greene-Johnson curvature term, `[k₁ D_R/(α₊ − H)] / [w + k₁ w_d/(0.3(1 + α₊))]` with +`α₊ = ½ + (−D_I)^(1/2)` (Benjamin Eq. 12); stabilizing for `D_R < 0`. +""" +function curvature_term(D_R::Real, H::Real, D_I::Real, w_d::Real, w::Real; k1::Real=1.7) + alpha_plus = 0.5 + sqrt(-D_I) + return (k1 * D_R / (alpha_plus - H)) / (w + k1 * w_d / (0.3 * (1 + alpha_plus))) +end + +""" + lutjens_curvature_term(D_R, D_I, K, W_d, w; model=:mercier_power) -> Float64 + +Curvature term of Lütjens, Luciani & Garbet (2001) in the form of Urso et al. (2010) Eq. 1, +`r_s Δ_GGJ = 6.35 D_R / √((W/r_s)² + 0.65 W_d²)`, with `W/r_s = w/K` and `W_d` the transport width in the same +`(r − r_s)/r_s` reference ([`layer_transport_width`](@ref)), referred to ψ_N as the outer Δ′ is (`K^(−2α)`, +`α = √(−D_I)`) and carried through [`classical_term`](@ref) as [`saturation_term`](@ref) is, so it has the units of +the other terms. Its `w → 0` limit is then `classical_term(−D_c K^(−2α), D_I, w)` with SLAYER's critical offset +`D_c = −√2 π^(3/2) D_R/W_d` (6.35/√0.65 = √2 π^(3/2)), so classical plus curvature change sign where SLAYER's gate does; +for `w ≫ W_d` its Δ′ form is the large-aspect-ratio `6.35 D_R/W` (La Haye 2006). The source has no `D_I ≠ −¼` form; +the factor of `classical_term` is the convention that keeps the seam at every `D_I`. Stabilizing for `D_R < 0`. +""" +function lutjens_curvature_term(D_R::Real, D_I::Real, K::Real, W_d::Real, w::Real; model::Symbol=:mercier_power) + return classical_term(K^(-2 * sqrt(-D_I)) * 6.35 * D_R / sqrt((w / K)^2 + 0.65 * W_d^2), D_I, w; model) +end + +""" + layer_transport_width(D_R, D_c) -> Float64 + +Transport width `W_d` in the `(r − r_s)/r_s` reference implied by a SLAYER critical offset, +`D_c = −√2 π^(3/2) D_R/W_d` (Connor et al. 2015 Eq. 59 with Fitzpatrick's W_d); `Inf` when `D_c = 0`. +""" +layer_transport_width(D_R::Real, D_c::Real) = D_c == 0 ? Inf : sqrt(2) * π^1.5 * (-D_R) / D_c + +""" + saturation_coefficient(J, r, dr_dpsi, q, dq_dpsi, psi; window=0.02, psi_lim=(0.0, 1.0), npts=11) -> NamedTuple + +Current-gradient coefficients of the finite-width saturation term (Fitzpatrick 2023 Eqs. 9.15–9.19, after Hastie, +Militello & Porcelli 2005): `α_s = −r q J_z′/s`, `β_s = −r² q J_z″/s` and `C_sat = 0.8 α_s² − 0.27 β_s − 0.09 α_s`, with +`s = r q′/q` and `′ = d/dr` in the radial label `r(ψ)` (derivative `dr_dpsi`), and `J(ψ)` the dimensionless current +density ([`current_density_profile`](@ref)). Eq. 9.17 prints `r q J_z″/s`; the `r²` form is the dimensionless one and +the one that reproduces Fitzpatrick's Fig. 9.2. `J_z′` and `J_z″` come from a least-squares quadratic in `r − r_s` over +`ψ ± window` (`npts` samples), the island-width scale of the theory's expansion in W. `converged` is false when the +doubled window leaves `psi_lim` or `C_sat` at half or double the window differs from it by more than 20%. +""" +function saturation_coefficient(J, r, dr_dpsi, q, dq_dpsi, psi::Real; window::Real=0.02, psi_lim=(0.0, 1.0), npts::Int=11) + rs, qs = r(psi), q(psi) + s = rs * dq_dpsi(psi) / (qs * dr_dpsi(psi)) + clamped = psi - 2window < psi_lim[1] || psi + 2window > psi_lim[2] + function coefficients(W) + W = min(W, 0.999 * (psi - psi_lim[1]), 0.999 * (psi_lim[2] - psi)) + xs = range(psi - W, psi + W; length=npts) + x = r.(xs) .- rs + c = hcat(ones(npts), x, x .^ 2) \ J.(xs) + alpha_s = -rs * qs * c[2] / s + beta_s = -rs^2 * qs * 2 * c[3] / s + return (; alpha_s, beta_s, C_sat=0.8 * alpha_s^2 - 0.27 * beta_s - 0.09 * alpha_s, J_s=c[1]) + end + mid = coefficients(window) + spread = maximum(abs(coefficients(f * window).C_sat - mid.C_sat) for f in (0.5, 2.0)) / max(abs(mid.C_sat), eps()) + return (; mid..., converged=!clamped && spread <= 0.2) +end + +""" + saturation_term(C_sat, K, D_I, w; model=:mercier_power) -> Float64 + +Finite-width saturation term: Fitzpatrick's replacement `r_sΔ′ → r_sΔ′(0) − C_sat W/r_s` (Eq. 9.15) referred to ψ_N +as the outer Δ′ is (`W/r_s = w/K`, `Δ′_ψ = K^(−2α) r_sΔ′`, `K = r_s dψ_N/dr`, `α = √(−D_I)`), then carried through +[`classical_term`](@ref), which is linear in Δ′: `−classical_term(C_sat K^(−1−2α) w, D_I, w)`. At `D_I = −¼` it is +`−C_sat w/K²`, so classical plus saturation vanish at Fitzpatrick's `W_pw/r_s = r_sΔ′(0)/C_sat` (Eq. 9.21). +""" +function saturation_term(C_sat::Real, K::Real, D_I::Real, w::Real; model::Symbol=:mercier_power) + return -classical_term(C_sat * K^(-1 - 2 * sqrt(-D_I)) * w, D_I, w; model) +end + +""" + polarization_term(w; c_pol, eps, rho_theta_psi, beta_p, gradient_ratio, omega_hat, omega_hat_pi, w_ib) -> Float64 + +Ion-polarization term of Benjamin Eq. C.1, +`−c_pol 3.28 ε^(3/2) I₁ (ρ_iθ,ψ/w)² (β_p/w) G² ω̂(ω̂ − ω̂_pi) w⁴/(w⁴ + w_ib⁴)`, with `I₁ = 0.71`, the poloidal Larmor +radius `ρ_iθ,ψ` in ψ_N (Benjamin's `(ρ_iθ/r_s)/(w/2ψ_s)` with its `dψ/dr = 2ψ_s/r_s` replaced by the layer label's), +`G = (n_i⁻¹ dn_i/dψ)/(q⁻¹ dq/dψ)` and the island frequency `ω̂ = ω/ω_*i` and `ω̂_pi = ω_*pi/ω_*i` normalized to the +ion diamagnetic frequency. Zero at `ω̂ = 0` and at `ω̂ = ω̂_pi`; stabilizing for `ω̂(ω̂ − ω̂_pi) > 0`. +""" +function polarization_term(w::Real; c_pol::Real, eps::Real, rho_theta_psi::Real, beta_p::Real, gradient_ratio::Real, omega_hat::Real, + omega_hat_pi::Real, w_ib::Real) + return -c_pol * 3.28 * eps^1.5 * 0.71 * (rho_theta_psi / w)^2 * (beta_p / w) * gradient_ratio^2 * + omega_hat * (omega_hat - omega_hat_pi) * w^4 / (w^4 + w_ib^4) +end + +""" + reduced_polarization_term(w; c_pol, C_pol, w_ib) -> Float64 + +Reduced polarization term `−c_pol w_pol²/w³` with `w_pol = C_pol w_ib`; it carries no island frequency. +""" +reduced_polarization_term(w::Real; c_pol::Real, C_pol::Real, w_ib::Real) = -c_pol * (C_pol * w_ib)^2 / w^3 + +""" + transport_width(w; c_wd, chi_perp, v_te, R0, q, q1, n, m, Wc_prefac) -> Float64 + +Transport threshold width `w_d = c_wd [(χ⊥/χ∥(w)) W_c/m²]^(1/4)` in ψ_N (Benjamin Eq. 14 in the Schlutt & Hegna flux +form, `W_c` from [`mre_geometry`](@ref)), with `χ∥(w) = v_Te λ∥(w)` and `λ∥(w) = R₀ q/(n |dq/dψ| w)` (Benjamin p. 4). +""" +function transport_width(w::Real; c_wd::Real, chi_perp::Real, v_te::Real, R0::Real, q::Real, q1::Real, n::Integer, m::Integer, Wc_prefac::Real) + chi_par = v_te * R0 * abs(q) / (abs(n * q1) * w) + return c_wd * ((chi_perp / chi_par) * abs(Wc_prefac) / m^2)^0.25 +end diff --git a/src/Tearing/ResonantForcing.jl b/src/Tearing/ResonantForcing.jl new file mode 100644 index 000000000..288a598f5 --- /dev/null +++ b/src/Tearing/ResonantForcing.jl @@ -0,0 +1,151 @@ +# The resonant forcing of a single-n coil-driven perturbed equilibrium, per rational surface, in the form a steady island +# equation takes it (Callen, UW-CPTC 16-4, 2016): the vacuum and total resonant fields, Δ_RMP and the vacuum island width, +# on the toroidal-flux radius ρ that Callen's widths use. The modified Rutherford equation's forcing term and Callen preset use it. + +module ResonantForcing + +import FastInterpolations + +import ...ForcingTerms +using ...Utilities.PhysicalConstants: M_P, E_CHG + +export toroidal_flux_radius, callen_polarization_width, resonant_forcing, forcing_thresholds + +""" + toroidal_flux_radius(equil, psi, q) -> (rho, drho) + +Toroidal-flux radius ρ = (Φ_tor/πB₀)^(1/2) [m] at `psi`, with Φ_tor = 2πψ₀ ∫q dψ_N and B₀ the vacuum field `bt0`, and +dρ/dψ_N = qψ₀/(B₀ρ) [m] with the surface's safety factor `q`. +""" +function toroidal_flux_radius(equil, psi::Real, q::Real) + B0 = equil.params.bt0 + rho = sqrt(FastInterpolations.integrate(equil.profiles.q_spline, first(equil.profiles.xs), psi) * 2 * equil.psio / B0) + return rho, q * equil.psio / (B0 * rho) +end + +""" + callen_polarization_width(equil, plasma, psi, q) -> (rho_i, w_pol) + +Ion gyroradius ρ_i = (2T_i/m_i)^(1/2)/(Z_i e B₀/m_i) and Callen's polarization width w_pol = ½ q ρ_i/(ρ/R₀)^(1/2), both in +metres, from a `layer_plasma` view of the kinetic data; R₀ is the axis major radius and B₀ the vacuum field `bt0`. +""" +function callen_polarization_width(equil, plasma, psi::Real, q::Real) + m_i, z_i, B0 = plasma.mu_i * M_P, plasma.z_i, equil.params.bt0 + rho_i = sqrt(2 * E_CHG * Float64(plasma.profiles.T_i(psi)) / m_i) / (z_i * E_CHG * B0 / m_i) + rho, _ = toroidal_flux_radius(equil, psi, q) + return rho_i, 0.5 * q * rho_i / sqrt(rho / equil.ro) +end + +# Vacuum resonant field |B_vac| [T] of the coil forcing on the surface ψ, harmonic m. +function _vacuum_resonant_field(pe, equil, psi, n, m, area) + rmp = pe.forcing + rmp isa ForcingTerms.RMPSource && rmp.ctrl.forcing_data_format == "coil" || + throw(ArgumentError("the resonant forcing needs a single coil forcing (forcing_data_format = \"coil\") for the vacuum resonant field")) + cfg = ForcingTerms.CoilConfig(rmp.ctrl) + grid = ForcingTerms.CoilForcingGrid(equil, cfg, n; psi) + amplitude = only(ForcingTerms.coil_forcing_modes(pe.coil_sets, grid, n, m, m)).amplitude + return abs(rmp.scale * amplitude) / area +end + +""" + resonant_forcing(pe) -> NamedTuple of vectors + +The resonant forcing on every resonant row of the single-n, coil-driven perturbed equilibrium `pe`: + + - `rational_psi`, `rational_q`, `m`, `n`, `surface_idx` - Surface identity; `surface_idx` indexes `pe.response.ffs.surfaces`. + - `rho`, `drho` - Toroidal-flux radius and dρ/dψ_N [m] ([`toroidal_flux_radius`](@ref)). + - `b_res`, `b_vac` - |total| and |vacuum| resonant area-weighted field [T]. + - `delta_rmp` - Δ_RMP = |Δ|/(2πR₀q) B₀/|B_res| [1/m of ρ], Δ the singular-current jump over 2πχ₁. + - `hw_isl`, `hw_v` - Total-field and vacuum island half-widths in ψ_N, hw_v = (4 B_vac A/(2π n q′ χ₁))^(1/2). + +In ψ_N the forcing of a full width w is then Δ_RMP (dρ/dψ_N) (2 hw_v/w)², which scales with B_vac. +""" +function resonant_forcing(pe) + resp = pe.response + ffs = resp.ffs + equil = ffs.equil + ffs.npert == 1 || throw(ArgumentError("the resonant forcing is single-n; the solve spans n = $(ffs.nlow):$(ffs.nhigh)")) + isempty(pe.island_half_width) && throw(ArgumentError("the resonant forcing needs the singular coupling of the perturbed equilibrium")) + R0, B0, chi1 = equil.ro, equil.params.bt0, 2π * equil.psio + nrow = length(resp.rational_psi) + rows = map(1:nrow) do row + psi, q, m, n, area = resp.rational_psi[row], resp.rational_q[row], resp.rational_m_res[row], resp.rational_n[row], resp.rational_area[row] + idx = resp.rational_surface_idx[row] + q1 = ffs.surfaces[idx].q1 + rho, drho = toroidal_flux_radius(equil, psi, q) + b_res = abs(pe.resonant_area_weighted_field[row]) + b_vac = _vacuum_resonant_field(pe, equil, psi, n, m, area) + return (; rational_psi=psi, rational_q=q, m, n, surface_idx=idx, rho, drho, b_res, b_vac, delta_rmp=abs(pe.delta_prime[row]) / (2π * R0 * q) * B0 / b_res, + hw_isl=pe.island_half_width[row], hw_v=sqrt(abs(4 * b_vac * area / (2π * n * q1 * chi1)))) + end + keys = (:rational_psi, :rational_q, :m, :n, :surface_idx, :rho, :drho, :b_res, :b_vac, :delta_rmp, :hw_isl, :hw_v) + return NamedTuple{keys}(Tuple([getfield(r, k) for r in rows] for k in keys)) +end + +# Golden-section extremum of f on [a, b]: the minimum, or with `sign = -1` the maximum. +function _golden(f, a, b; sign=1, tol=1e-10) + phi = (sqrt(5) - 1) / 2 + c, d = b - phi * (b - a), a + phi * (b - a) + fc, fd = sign * f(c), sign * f(d) + while b - a > tol * (abs(a) + abs(b)) + if fc < fd + b, d, fd = d, c, fc + c = b - phi * (b - a) + fc = sign * f(c) + else + a, c, fc = c, d, fd + d = a + phi * (b - a) + fd = sign * f(d) + end + end + x = (a + b) / 2 + return x, f(x) +end + +""" + forcing_thresholds(F0, k, w_lo, w_grid, w_m0) -> (; b_crit, b_trigger, w_crit) + +Critical fields of a forcing k B_vac/w² added to an unforced steady island equation `F0(w)`, with the island locked in +phase with the field (Callen's assumption). The forced F is positive where k B_vac > G(w) = −F0(w) w², so: + + - `b_crit` - The vacuum field at which a forced window F > 0 first opens below the intrinsic seed threshold `w_m0` + (anywhere on the grid when `w_m0` is `NaN`), apart from the intrinsic region: the smallest local minimum of G, + `w_lo` included, at width `w_crit`. For Callen's F0 it is his B_crit. It is 0 when F0 > 0 at `w_lo`, or when G + falls all the way to `w_m0` (any field lowers the seed threshold); `NaN` when it falls all the way to the top of + the grid without a `w_m0`. Without a term that beats 1/w² at small w (Callen's −w_pol²/w³ does), G rises from + `w_lo`, the window opens there (`w_crit == w_lo`) and `b_crit` depends on `w_lo`, the layer width when the gate + sets it. + - `b_trigger` - The field at which the window merges with the region above `w_m0`, so that an island in it grows + past `w_m0` towards the intrinsic saturated width: the maximum of G on `[w_crit, w_m0]` (on `[w_lo, w_m0]` + without a window, the field at which no seed is needed). It is 0 when F0 > 0 at `w_lo`, and `NaN` without `w_m0`. + +G is sampled on `w_lo` and the grid points above it, and each interior extremum is refined by golden-section search. +""" +function forcing_thresholds(F0, k::Real, w_lo::Real, w_grid::AbstractVector, w_m0::Real) + G(w) = -F0(w) * w^2 + none = (; b_crit=NaN, b_trigger=NaN, w_crit=NaN) + (k > 0 && isfinite(k)) || return none + top = isfinite(w_m0) ? w_m0 : maximum(w_grid) + nodes = vcat(Float64(w_lo), filter(w -> w_lo < w < top, w_grid), top) + vals = G.(nodes) + any(!isfinite, vals) && return none + vals[1] < 0 && return (; b_crit=0.0, b_trigger=0.0, w_crit=NaN) + n = length(nodes) + interior = [_golden(G, nodes[j-1], nodes[j+1]) for j in 2:(n-1) if vals[j] <= vals[j-1] && vals[j] <= vals[j+1]] + minima = n > 1 && vals[1] <= vals[2] ? vcat([(nodes[1], vals[1])], interior) : interior + if isempty(minima) + isfinite(w_m0) || return none + w_crit, b_crit = NaN, 0.0 + else + w_crit, g_crit = minima[argmin(last.(minima))] + b_crit = g_crit / k + isfinite(w_m0) || return (; b_crit, b_trigger=NaN, w_crit) + end + # The highest barrier of G between the window (or w_lo) and the intrinsic seed threshold. + start = isnan(w_crit) ? nodes[1] : w_crit + j = argmax([w >= start ? v : -Inf for (w, v) in zip(nodes, vals)]) + g_max = 1 < j < n ? last(_golden(G, max(nodes[j-1], start), nodes[j+1]; sign=-1)) : vals[j] + return (; b_crit, b_trigger=max(g_max, vals[j]) / k, w_crit) +end + +end # module ResonantForcing diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 0055cb273..02317fa08 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -6,6 +6,8 @@ # extraction (consumes any InnerLayerModel) # CriticalResonantField -- real-Q torque-balance scan for the critical # resonant field of error-field penetration +# ModifiedRutherford -- terms of the steady modified Rutherford equation and +# F(w) per surface, `solve(ModifiedRutherfordProblem(ffs; ...))` # TearingProblem / TearingConfig / TearingResult -- the public solve, # `solve(TearingProblem(ffs; ...), model; config)`, with one # method per inner-layer model (SLAYERModel, GGJModel) @@ -33,15 +35,19 @@ import ..InnerLayer as InnerLayer using .InnerLayer: InnerLayerParameters, InnerLayerResponse, solve_inner, layer_response, SLAYERModel, SLAYERParameters, build_slayer_inputs, GGJModel, GGJParameters, - LayerWidths, slayer_layer_thickness -using .InnerLayer.SLAYER: _slayer_profiles, _layer_zeff, _layer_ion_mass, _layer_mass_per_electron, _build_resistivity_model + LayerWidths, slayer_layer_thickness, layer_plasma include("LayerOverlap.jl") include("Dispersion/Dispersion.jl") include("CriticalResonantField/CriticalResonantField.jl") +include("ResonantForcing.jl") +include("ModifiedRutherford/ModifiedRutherford.jl") import .Dispersion as Dispersion import .CriticalResonantField as CriticalResonantField +import .ModifiedRutherford as ModifiedRutherford +import .ResonantForcing as ResonantForcing +using .ModifiedRutherford: ModifiedRutherfordProblem, ModifiedRutherfordResult using .Dispersion: SurfaceCoupling, surface_coupling, MultiSurfaceCoupling, multi_surface_coupling, ScanResult, brute_force_scan, @@ -57,7 +63,8 @@ include("SLAYERTearing.jl") include("GGJTearing.jl") include("HDF5Output.jl") -export InnerLayer, Dispersion, CriticalResonantField +export InnerLayer, Dispersion, CriticalResonantField, ModifiedRutherford +export ModifiedRutherfordProblem, ModifiedRutherfordResult export resistive_layer_overlap, layer_overlap_scan, LayerOverlapScan export TearingProblem, TearingConfig, TearingResult export run_slayer_from_inputs diff --git a/src/Tearing/TearingProblem.jl b/src/Tearing/TearingProblem.jl index 1d1b5a75d..dea53e673 100644 --- a/src/Tearing/TearingProblem.jl +++ b/src/Tearing/TearingProblem.jl @@ -39,16 +39,15 @@ end TearingProblem(ffs; eta=nothing, rho=nothing, kwargs...) = TearingProblem{typeof(ffs),typeof(eta),typeof(rho)}(; ffs=ffs, eta=eta, rho=rho, kwargs...) -function _tearing_params(model::SLAYERModel, equil, surfaces, loaded) - data = equil.kinetic +function _tearing_params(model::SLAYERModel, equil, surfaces, plasma) + data = plasma.data # χ⊥,e/χ⊥,i/χ_φ from the kinetic file when present, else the kinetic data's defaults. - chi = _slayer_chi_inputs(loaded, model, data) - zeff = _layer_zeff(loaded, data) + chi = _slayer_chi_inputs(plasma, model, data) # B_T = F(ψ)/(2π·R₀) per surface comes from the equilibrium's F-spline (bt=nothing). - return build_slayer_inputs(equil, surfaces, loaded.profiles; - mu_i=_layer_ion_mass(data), - mu_rho=_layer_mass_per_electron(data, zeff), - zeff=zeff, + return build_slayer_inputs(equil, surfaces, plasma.profiles; + plasma.mu_i, + plasma.mu_rho, + plasma.zeff, chi_perp_e=chi.chi_perp_e, chi_perp_i=chi.chi_perp_i, chi_perp_Dc=model.chi_perp_Dc, @@ -59,8 +58,8 @@ function _tearing_params(model::SLAYERModel, equil, surfaces, loaded) dgeo_val=model.D_geo, dc_type=model.D_c_type, theta=model.theta_sample, - resistivity_model=_build_resistivity_model(data.resistivity_model), - lnLambda_form=data.lnLambda_form) + plasma.resistivity_model, + plasma.lnLambda_form) end """ @@ -93,10 +92,10 @@ function CommonSolve.solve(prob::TearingProblem, model::SLAYERModel; config::Tea equil.kinetic === nothing && error("TearingProblem: the equilibrium carries no kinetic profiles; set kinetic_file in [Equilibrium] or attach_kinetic_profiles!") - loaded = _slayer_profiles(equil.kinetic, model) + plasma = layer_plasma(equil; control=model) _warn_past_layer_overlap(ffs, model) - params = _tearing_params(model, equil, surfaces, loaded) + params = _tearing_params(model, equil, surfaces, plasma) # Prefer the full inter-surface Δ' matrix; fall back to a diagonal of each surface's scalar delta_prime. dpm = ffs.delta_prime === nothing ? Matrix{ComplexF64}(undef, 0, 0) : ffs.delta_prime.matrix @@ -125,8 +124,8 @@ function CommonSolve.solve(prob::TearingProblem, model::SLAYERModel; config::Tea rational_q = Float64[surfaces[p.ising].q for p in params] # E×B rotation per unit n at each surface, for the Doppler shifts and the lab-frame roots; # empty when the kinetic file carries none. - omega_E = loaded.has_omega_E ? Float64[loaded.profiles.omega_E(ψ) for ψ in rational_psi] : Float64[] - loaded.has_omega_E || @warn "SLAYER: the kinetic file has no omega_E; surfaces without an omega_E_kHz override assume Ω_E = 0, so their lab-frame ω equals their E×B-frame ω." + omega_E = plasma.has_omega_E ? Float64[plasma.profiles.omega_E(ψ) for ψ in rational_psi] : Float64[] + plasma.has_omega_E || @warn "SLAYER: the kinetic file has no omega_E; surfaces without an omega_E_kHz override assume Ω_E = 0, so their lab-frame ω equals their E×B-frame ω." result = run_slayer_from_inputs(model, params, dp, prob, config; rational_psi=rational_psi, rational_q=rational_q, omega_E=omega_E) return result diff --git a/src/Utilities/Output.jl b/src/Utilities/Output.jl index bf75c28f8..f255e717d 100644 --- a/src/Utilities/Output.jl +++ b/src/Utilities/Output.jl @@ -28,13 +28,20 @@ function deck_sections end """ record_inputs!(target, result) -Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, replacing any section of -the same name, so every saved result adds the deck it was computed from. +Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, so every saved result adds +the deck it was computed from. A result replaces its section's scalars and the sub-tables it states, and keeps recorded +sub-tables it does not state (a section stated only as sub-tables keeps the recorded scalars too). """ function record_inputs!(target::Union{HDF5.File,HDF5.Group}, result) path = "Input/gpec_toml_raw" deck = haskey(target, path) ? TOML.parse(read(target[path])) : Dict{String,Any}() - merge!(deck, deck_sections(result)) + # A result owns its section's scalars and the sub-tables it states; recorded sub-tables it does not state stay. + for (k, v) in deck_sections(result) + old = get(deck, k, nothing) + deck[k] = + !(v isa AbstractDict && old isa AbstractDict) ? v : + all(x -> x isa AbstractDict, values(v)) ? merge(old, v) : merge(filter(p -> p.second isa AbstractDict, old), v) + end haskey(target, path) && HDF5.delete_object(target, path) target[path] = sprint(TOML.print, deck) return nothing diff --git a/test/runtests_modified_rutherford.jl b/test/runtests_modified_rutherford.jl new file mode 100644 index 000000000..24eedeffa --- /dev/null +++ b/test/runtests_modified_rutherford.jl @@ -0,0 +1,178 @@ +using TOML + +# The modified Rutherford terms: the limits and signs each term must show +# alone, and the surface geometry against Fortran RDCON (MRE_flag output) on the DIII-D-like example. +@testset "ModifiedRutherford terms" begin + MR = GeneralizedPerturbedEquilibrium.ModifiedRutherford + + @testset "classical term: Mercier power and its zero-pressure limit" begin + @test MR.classical_term(-3.0, -0.25, 0.01) ≈ -3.0 + @test MR.classical_term(-3.0, -0.3, 0.01; model=:constant) == -3.0 + # Slope of log Δ′(w) is −2α₋ = 2√(−D_I) − 1. + D_I = -0.3 + slope = log(MR.classical_term(2.0, D_I, 0.02) / MR.classical_term(2.0, D_I, 0.01)) / log(2) + @test slope ≈ 2 * sqrt(-D_I) - 1 rtol = 1e-12 + end + + @testset "bootstrap term: positive, peaked, k1·D_nc/w far above w_d" begin + D_nc, w_d, k1 = 0.4, 1e-3, 1.7 + w = exp.(range(log(1e-5), log(0.2); length=400)) + b = MR.bootstrap_term.(D_nc, w_d, w; k1) + @test all(>(0), b) + @test w[argmax(b)] ≈ sqrt(k1 / 0.5) * w_d rtol = 0.02 + @test MR.bootstrap_term(D_nc, w_d, 1e3 * w_d; k1) ≈ k1 * D_nc / (1e3 * w_d) rtol = 1e-5 + @test MR.bootstrap_term(D_nc, w_d, 1e-3 * w_d; k1) ≈ 0.5 * D_nc * 1e-3 * w_d / w_d^2 rtol = 1e-5 + end + + @testset "curvature term: stabilizing for D_R < 0" begin + @test MR.curvature_term(-0.02, 0.01, -0.27, 1e-3, 0.01) < 0 + @test MR.curvature_term(0.02, 0.01, -0.27, 1e-3, 0.01) > 0 + end + + @testset "polarization term: sign follows ω̂(ω̂ − ω̂_pi), zero at ω̂ = 0" begin + args = (; c_pol=1.0, eps=0.2, rho_theta_psi=0.012, beta_p=0.5, gradient_ratio=-0.8, omega_hat_pi=1.5, w_ib=1e-3) + pol(ω̂) = MR.polarization_term(0.02; args..., omega_hat=ω̂) + @test pol(0.0) == 0 + @test pol(1.5) == 0 + @test pol(0.7) > 0 # ω̂ inside (0, ω̂_pi): destabilizing + @test pol(2.0) < 0 && pol(-1.0) < 0 + @test MR.reduced_polarization_term(0.01; c_pol=1.0, C_pol=1.9, w_ib=1e-3) ≈ -(1.9e-3)^2 / 1e-6 + end + + @testset "Lütjens curvature: SLAYER's critical offset at small w, 6.35 D_R/W at large w" begin + D_R, D_c, K = -0.02, 1.5, 0.9 + W_d = MR.layer_transport_width(D_R, D_c) + @test MR.lutjens_curvature_term(D_R, -0.25, K, W_d, 1e-9) * K ≈ -D_c rtol = 1e-3 + @test MR.lutjens_curvature_term(D_R, -0.25, K, W_d, 1e3 * W_d * K) ≈ 6.35 * D_R / (1e3 * W_d * K) rtol = 1e-6 + # Away from D_I = −¼, classical plus curvature still change sign where SLAYER's gate does, at any small w. + for D_I in (-0.3, -0.4), w in (1e-6, 1e-3 * W_d * K), f in (0.99, 1.01) + dp = f * D_c * K^(-2 * sqrt(-D_I)) + @test sign(MR.classical_term(dp, D_I, w) + MR.lutjens_curvature_term(D_R, D_I, K, W_d, w)) == sign(f - 1) + end + @test MR.layer_transport_width(D_R, 0.0) == Inf + end + + @testset "saturation term: Wesson cylinder of Fitzpatrick 2023 Fig. 9.2" begin + # q0 = 0.8, q(a) = 4.4, r = √ψ (a = 1), J_z = (2/q0)(1 − r²)^ν; reference α_s, β_s from the independent shooting script. + q0, qa = 0.8, 4.4 + nu = qa / q0 - 1 + q(x) = qa * x / (1 - (1 - x)^(nu + 1)) + dq(x) = (q(x + 1e-7) - q(x - 1e-7)) / 2e-7 + lo, hi = 0.1, 0.9 + for _ in 1:60 + mid = (lo + hi) / 2 + q(mid) < 2 ? (lo = mid) : (hi = mid) + end + psi_s = (lo + hi) / 2 + sc = MR.saturation_coefficient(x -> (2 / q0) * (1 - x)^nu, sqrt, x -> 0.5 / sqrt(x), q, dq, psi_s; window=0.005) + @test sc.converged + @test sqrt(psi_s) ≈ 0.65942681 rtol = 1e-6 + @test sc.alpha_s ≈ 1.642673 rtol = 1e-4 + @test sc.beta_s ≈ -7.204691 rtol = 1e-3 + # Classical plus saturation vanish at W_pw/r_s = r_sΔ′(0)/C_sat (Eq. 9.21); r_sΔ′(0) = 1.64 from the cylindrical + # tearing equation, which puts W_pw/a at the 0.27 peak of the r_w = a curve. K = r_s dψ/dr = 2ψ_s. + E_ss, K = 1.64, 2 * psi_s + grid = exp.(range(log(1e-4), log(2.0); length=200)) + r = MR.mre_roots(w -> MR.classical_term(E_ss / K, -0.25, w) + MR.saturation_term(sc.C_sat, K, -0.25, w), 1e-4, grid) + @test r.w_sat / K ≈ E_ss / sc.C_sat rtol = 1e-8 + @test sqrt(psi_s) * r.w_sat / K ≈ 0.27 rtol = 0.02 + end + + @testset "transport width grows as w^(1/4)" begin + kw = (; c_wd=1.0, chi_perp=1.0, v_te=5e7, R0=1.7, q=2.0, q1=3.0, n=1, m=2, Wc_prefac=3e4) + @test MR.transport_width(16e-3; kw...) / MR.transport_width(1e-3; kw...) ≈ 2 rtol = 1e-12 + end + + @testset "roots: seed threshold rises through zero, saturation falls through it" begin + grid = exp.(range(log(1e-5), log(0.2); length=60)) + r = MR.mre_roots(w -> (w - 3e-3) * (0.05 - w) / w^2, 1e-5, grid) + @test r.w_m ≈ 3e-3 rtol = 1e-10 + @test r.w_sat ≈ 0.05 rtol = 1e-10 + @test !r.fully_stabilized && !r.unstable_any_seed + @test MR.mre_roots(w -> -1 / w, 1e-5, grid).fully_stabilized + s = MR.mre_roots(w -> (0.05 - w) / w, 1e-5, grid) + @test s.unstable_any_seed && isnan(s.w_m) + @test s.w_sat ≈ 0.05 rtol = 1e-10 + # Roots at 2e-3 (falling), 3e-3 (rising), 0.02 (falling), 0.05 (rising): from w_lo the island saturates at 2e-3; F > 0 at the top. + f4(w) = (2e-3 - w) * (w - 3e-3) * (0.02 - w) * (w - 0.05) / w^4 + u = MR.mre_roots(f4, 1e-5, grid) + @test u.unstable_any_seed + @test u.w_m ≈ 3e-3 rtol = 1e-10 + @test u.w_sat ≈ 2e-3 rtol = 1e-10 + @test u.unsaturated + # Every branch: 2e-3 from any seed below 3e-3, 0.02 from seeds in (3e-3, 0.05), growth past the grid above 0.05. + @test u.w_sat_all ≈ [2e-3, 0.02] rtol = 1e-10 + @test isnan(u.w_seed_all[1]) + @test u.w_seed_all[2] ≈ 3e-3 rtol = 1e-10 + @test u.w_seed_unsaturated ≈ 0.05 rtol = 1e-10 + # With F < 0 at w_lo the island starts at w_m = 2e-3 and saturates at the next falling root, 3e-3, not at 0.05. + v = MR.mre_roots(w -> -f4(w), 1e-5, grid) + @test v.w_m ≈ 2e-3 rtol = 1e-10 + @test v.w_sat ≈ 3e-3 rtol = 1e-10 + @test !v.unsaturated + @test v.w_sat_all ≈ [3e-3, 0.05] rtol = 1e-10 + @test v.w_seed_all ≈ [2e-3, 0.02] rtol = 1e-10 + @test isnan(v.w_seed_unsaturated) + # One rising root and F > 0 at the top: w_m is also the seed of unbounded growth; a stable branch has none. + @test MR.mre_roots(w -> (w - 3e-3) / w, 1e-5, grid).w_seed_unsaturated ≈ 3e-3 rtol = 1e-10 + @test isnan(s.w_seed_unsaturated) + @test s.w_sat_all ≈ [0.05] rtol = 1e-10 + end + + @testset "Callen's steady state: the root finder recovers his Fig. 8 widths" begin + # Callen 2016 Fig. 8 (DIII-D 158115, ρ = 71.88 cm), the cubic Δ_c w³ + Δ_RMP w_v² w − w_pol² = 0 over w³, in cm: + # full widths w_min ≈ 1.28 cm, w_sat ≈ 2.24 cm. + a, b, c = -16 / 71.88, 26.4 / 71.88 * 2.45^2, -1.5^2 + r = MR.mre_roots(w -> a + b / w^2 + c / w^3, 0.1, exp.(range(log(0.1), log(10.0); length=200))) + @test abs(a * r.w_m^3 + b * r.w_m + c) < 1e-10 + @test abs(a * r.w_sat^3 + b * r.w_sat + c) < 1e-10 + @test r.w_m ≈ 1.28 rtol = 0.1 + @test r.w_sat ≈ 2.24 rtol = 0.1 + end + + @testset "forcing thresholds: Callen's B_crit, and a forced window below an intrinsic seed" begin + RF = GeneralizedPerturbedEquilibrium.Tearing.ResonantForcing + grid = exp.(range(log(1e-5), log(0.2); length=200)) + # Callen's F0 = Δ_c + c/w³: the window opens where the cubic's roots merge, b = 3|a| (c/2a)^(2/3). + a, c = -40.0, -(1.5e-2)^2 + th = RF.forcing_thresholds(w -> a + c / w^3, 2.0, 1e-5, grid, NaN) + @test th.b_crit ≈ 3 * abs(a) * (c / (2a))^(2 / 3) / 2 rtol = 1e-8 + @test th.w_crit ≈ cbrt(c / (2a)) rtol = 1e-6 + @test isnan(th.b_trigger) + # G = −F0 w² with a local minimum near 0.01 and a barrier before the intrinsic seed threshold 0.03. + G(w) = (0.03 - w) * ((w - 0.01)^2 + 1e-5) + th = RF.forcing_thresholds(w -> -G(w) / w^2, 2.0, 1e-5, grid, 0.03) + fine = range(1e-5, 0.03; length=100_001) + @test th.b_crit ≈ minimum(G, filter(<(0.02), fine)) / 2 rtol = 1e-6 + @test th.b_trigger ≈ maximum(G, filter(>=(th.w_crit), fine)) / 2 rtol = 1e-6 + @test th.w_crit < 0.03 + # G rising from w_lo opens the window there; G falling all the way to the seed threshold has none, so any field + # lowers it and b_trigger is the field that needs no seed. F0 > 0 at w_lo: no seed is needed. + th = RF.forcing_thresholds(w -> (w - 0.03) / w, 2.0, 1e-5, grid, 0.03) + @test th.w_crit == 1e-5 + @test th.b_crit ≈ 1e-5 * (0.03 - 1e-5) / 2 rtol = 1e-12 + @test th.b_trigger ≈ 0.015^2 / 2 rtol = 1e-8 + th = RF.forcing_thresholds(w -> (w - 0.03) / w^3, 2.0, 1e-5, grid, 0.03) + @test th.b_crit == 0 && isnan(th.w_crit) + @test th.b_trigger ≈ (0.03 - 1e-5) / 1e-5 / 2 rtol = 1e-12 + @test RF.forcing_thresholds(w -> 1 / w, 2.0, 1e-5, grid, NaN).b_trigger == 0 + @test isnan(RF.forcing_thresholds(w -> -1 / w, 0.0, 1e-5, grid, NaN).b_crit) + end + + @testset "surface geometry matches Fortran RDCON on the DIII-D-like example" begin + # Fortran values at q = 2 (ψ_N = 0.5183), interpolated from rdcon_output_n1.nc with MRE_flag = t and Z_eff = 1.52. + dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_gal_resistive_example") + deck = TOML.parsefile(joinpath(dir, "gpec.toml")) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(GeneralizedPerturbedEquilibrium.equilibrium_config(deck, dir), + GeneralizedPerturbedEquilibrium.analytic_config(deck)) + psi = 0.51830565 + g = MR.mre_geometry(equil, psi, equil.profiles.q_deriv(psi)) + @test g.f_trap ≈ 0.61419 rtol = 2e-3 + @test g.overbar_Rsq ≈ 3.0113 rtol = 2e-3 + @test g.avg_bpol ≈ 0.28709 rtol = 2e-3 + @test g.Wc_prefac ≈ 33565 rtol = 2e-3 + D_nc, mu_ratio = MR.bootstrap_coefficient(g, 1.52) + @test mu_ratio ≈ 0.68255 rtol = 2e-3 + @test D_nc ≈ 0.38065 rtol = 2e-3 + end +end diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index f38ada7aa..5d71b69ee 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -45,8 +45,9 @@ include("display_check.jl") ret = GeneralizedPerturbedEquilibrium.main([dir]) # One field per stage, nothing for the sections the deck does not hold. - @test keys(ret) == (:equil, :local_stability, :ffs, :match, :resp, :pe, :kinetic, :error_fields, :tearing) - @test ret.equil === ret.ffs.equil && all(isnothing, (ret.match, ret.resp, ret.pe, ret.kinetic, ret.error_fields, ret.tearing)) + @test keys(ret) == (:equil, :local_stability, :ffs, :match, :resp, :pe, :kinetic, :error_fields, :tearing, :modified_rutherford) + @test ret.equil === ret.ffs.equil && + all(isnothing, (ret.match, ret.resp, ret.pe, ret.kinetic, ret.error_fields, ret.tearing, ret.modified_rutherford)) ffs = ret.ffs forward_result[] = ffs @test ffs isa FFS.ForceFreeStatesResult diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 89b27dfbd..ab38e7718 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -179,13 +179,12 @@ # The model key, through the same assembly TearingProblem uses. Tearing = GeneralizedPerturbedEquilibrium.Tearing - loaded = (; profiles, chi_perp_e=nothing, chi_perp_i=nothing, chi_tor=nothing, zeff=nothing, zeff_source=:scalar) model(; c...) = SLAYERModel(; P_perp_model=:chi_perp_e, D_R=0.01, D_c_type=:fitzpatrick, c...) attached = equil.kinetic equil.kinetic = Equilibrium.KineticProfileData(; psi=psi_pts, n_e=fill(5.0e19, length(psi_pts)), T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), T_i=1000.0 .* (1.0 .- 0.6 .* psi_pts), chi_e_default=2.0, chi_phi_default=1.5) - params(m) = Tearing._tearing_params(m, equil, sings, loaded)[1] + params(m) = Tearing._tearing_params(m, equil, sings, GeneralizedPerturbedEquilibrium.InnerLayer.layer_plasma(equil; control=m, log=false))[1] without, with = params(model()), params(model(; chi_perp_Dc=6.0)) equil.kinetic = attached @test with.dc_tmp != without.dc_tmp diff --git a/test/runtests_slayer_profile_transport.jl b/test/runtests_slayer_profile_transport.jl index 0e9f0b25a..8ec9c40c1 100644 --- a/test/runtests_slayer_profile_transport.jl +++ b/test/runtests_slayer_profile_transport.jl @@ -1,7 +1,7 @@ using Test using GeneralizedPerturbedEquilibrium.Equilibrium: KineticProfileData using GeneralizedPerturbedEquilibrium: SLAYERModel -using GeneralizedPerturbedEquilibrium.Tearing: _slayer_profiles +using GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER: _slayer_profiles @testset "SLAYER loads only transport profiles used by its model" begin psi = collect(range(0.0, 1.0; length=5)) diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 434047642..a6b4a1f02 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -117,7 +117,7 @@ @testset "SLAYER view of the kinetic profiles" begin kin_path = joinpath(dirname(@__DIR__), "examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode_kinetic.h5") data = GeneralizedPerturbedEquilibrium.Equilibrium.read_kinetic_file(kin_path) - view = Tearing._slayer_profiles(data, SLAYERModel()) + view = GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER._slayer_profiles(data, SLAYERModel()) j = length(data.psi) ÷ 2 @test view.profiles.n_e(data.psi[j]) ≈ data.n_e[j] @test view.profiles.omega_E(data.psi[j]) ≈ data.omega_E[j] @@ -126,16 +126,17 @@ @test view.has_omega_E # Without an E×B rotation profile the view records its absence and carries zero rotation. no_rotation = GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_e=data.T_e, T_i=data.T_i) - bare_view = Tearing._slayer_profiles(no_rotation, SLAYERModel()) + bare_view = GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER._slayer_profiles(no_rotation, SLAYERModel()) @test !bare_view.has_omega_E @test bare_view.profiles.omega_E(data.psi[j]) == 0.0 - @test_throws "missing required dataset 'T_i'" Tearing._slayer_profiles( + @test_throws "missing required dataset 'T_i'" GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER._slayer_profiles( GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_e=data.T_e), SLAYERModel()) # A non-positive chi_e stops the run only where it is used: chi_perp_Dc replaces it in the critical Δ. bad_chi_e = GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_e=data.T_e, T_i=data.T_i, chi_e=fill(-1.0, length(data.psi))) - @test_throws "chi_e" Tearing._slayer_profiles(bad_chi_e, SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal)) - @test Tearing._slayer_profiles(bad_chi_e, SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal, chi_perp_Dc=1.4)).chi_perp_e === nothing + @test_throws "chi_e" GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER._slayer_profiles(bad_chi_e, SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal)) + @test GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER._slayer_profiles(bad_chi_e, SLAYERModel(; P_perp_model=:tau_E, D_c_type=:toroidal, chi_perp_Dc=1.4)).chi_perp_e === + nothing end # Δ′ is unified across formalisms, so a Galerkin run feeds SLAYER exactly as a Riccati one diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 2c17a146b..aa9a17b49 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -243,6 +243,13 @@ include("display_check.jl") # The deck reads and normalizes its forcing file exactly as a fresh RMPField does. @test pe.resonant_area_weighted_field ≈ reference.resonant_area_weighted_field rtol = 1e-10 + # The MRE's seed field inverts the perturbed equilibrium's vacuum-island width at every surface. + MRE = GeneralizedPerturbedEquilibrium.Tearing.ModifiedRutherford + for (row, k) in enumerate(resp.rational_surface_idx) + sk = ffs.surfaces[k] + @test MRE.vacuum_island_field(equil, resp.rational_psi[row], resp.rational_n[row], sk.q1, 2 * pe.island_half_width[row]) ≈ + abs(pe.resonant_area_weighted_field[row]) rtol = 1e-8 + end @test pe.forcing_b ≈ reference.forcing_b rtol = 1e-10 # One response serves any number of fields: forcing it again reproduces the result. @@ -561,6 +568,192 @@ include("display_check.jl") end end + @testset "modified Rutherford under [Tearing.ModifiedRutherford]" begin + # The MRE alone: a [Tearing] holding only the sub-table runs no layer solve, and its deck replays. + settings = (; polarization=true, omega_model=:exb, omega_factor=0.5, delta_prime_factor=0.9, D_nc_factor=1.1, rho_theta_factor=1.2, nw=40, + curvature_model=:benjamin) + mktempdir() do dir + psi = collect(range(0.0, 1.0; length=101)) + h5open(joinpath(dir, "kinetic.h5"), "w") do f + f["psi"] = psi + f["n_e"] = f["n_i"] = 5e19 .* (1 .- 0.8 .* psi) + f["T_e"] = f["T_i"] = 2e3 .* (1 .- 0.9 .* psi) + f["omega_E"] = 1e4 .* (1 .- psi) + end + ric = solve(attach_kinetic_profiles!(deepcopy(equil), joinpath(dir, "kinetic.h5")), riccati; nn=1, wall=wall, ffs_kwargs...) + api = solve(ModifiedRutherfordProblem(ric; settings...)) + deck_dir = _stage_deck(dir; eq_overrides=Dict{String,Any}("kinetic_file" => "kinetic.h5"), + ffs_overrides=Dict{String,Any}("integrator" => "Riccati"), + extra_sections=Dict{String,Any}( + "Tearing" => Dict{String,Any}("ModifiedRutherford" => + Dict{String,Any}(String(k) => (v isa Symbol ? String(v) : v) for (k, v) in pairs(settings))) + )) + reference = GPEC.main([deck_dir]) + + @test reference.tearing === nothing + @test api.F ≈ reference.modified_rutherford.F rtol = 1e-10 + @test api.polarization ≈ reference.modified_rutherford.polarization rtol = 1e-10 + @test isequal(api.w_m, reference.modified_rutherford.w_m) && isequal(api.w_sat, reference.modified_rutherford.w_sat) + @test api.fully_stabilized == reference.modified_rutherford.fully_stabilized + @test all(<(0), api.D_R) && all(isfinite, api.F) + h5open(joinpath(dir, "gpec.h5")) do f + @test read(f["Tearing/ModifiedRutherford/Terms/F"]) ≈ api.F rtol = 1e-12 + recorded = TOML.parse(read(f["Input/gpec_toml_raw"])) + @test recorded["Tearing"] == Dict{String,Any}("ModifiedRutherford" => recorded["Tearing"]["ModifiedRutherford"]) + @test recorded["Tearing"]["ModifiedRutherford"]["omega_factor"] == 0.5 + @test recorded["Tearing"]["ModifiedRutherford"]["omega_model"] == "exb" + @test recorded["Tearing"]["ModifiedRutherford"]["curvature_model"] == "benjamin" + @test isequal(read(f["Tearing/ModifiedRutherford/PerSurface/w_m"]), api.w_m) + end + end + + # With [Tearing], main hands the linear solve to the MRE as its gate, as the API keyword does. + sl_dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_SLAYER_example") + mktempdir() do dir + deck = TOML.parsefile(joinpath(sl_dir, "gpec.toml")) + deck["Run"]["output_file"] = joinpath(dir, "gpec.h5") + deck["Tearing"]["ModifiedRutherford"] = Dict{String,Any}("nw" => 40) + for key in ("eq_filename", "kinetic_file") + deck["Equilibrium"][key] = abspath(joinpath(sl_dir, deck["Equilibrium"][key])) + end + open(io -> TOML.print(io, deck), joinpath(dir, "gpec.toml"), "w") + reference = GPEC.main([dir]) + api = solve(ModifiedRutherfordProblem(reference.ffs; tearing=reference.tearing, nw=40)) + @test reference.modified_rutherford.problem.tearing === reference.tearing + @test isequal(api.gamma, reference.modified_rutherford.gamma) && all(isfinite, api.gamma) + @test isequal(api.w_m, reference.modified_rutherford.w_m) + # The default Lütjens curvature takes SLAYER's toroidal critical offset, so it equals a toroidal gate's at the same χ⊥. + # χ⊥ varies by surface here, so the gate is built with the 3/1 surface's χ⊥ and compared there. + @test all(isfinite, api.D_c) && all(>(0), api.W_d_layer) + k = findfirst(==(3), api.m) + gate = solve(TearingProblem(reference.ffs), SLAYERModel(; D_c_type=:toroidal, chi_perp_Dc=api.chi_perp[k])) + row = findfirst(x -> isapprox(x, api.rational_psi[k]; atol=1e-8), gate.rational_psi) + @test gate.params[row].dc_tmp ≈ api.D_c[k] rtol = 1e-10 + ben = solve(ModifiedRutherfordProblem(reference.ffs; tearing=reference.tearing, nw=40, curvature_model=:benjamin)) + @test all(isnan, ben.D_c) && !(ben.curvature ≈ api.curvature) + @test_throws "curvature_model" solve(ModifiedRutherfordProblem(reference.ffs; curvature_model=:other)) + # A TearingResult saved after the MRE keeps Tearing/ModifiedRutherford and its recorded sub-table. + h5open(joinpath(dir, "reversed.h5"), "w") do f + save_to_hdf5(f, api) + save_to_hdf5(f, reference.tearing) + @test isequal(read(f["Tearing/ModifiedRutherford/PerSurface/w_m"]), api.w_m) && haskey(f, "Tearing/PerSurface") + recorded = TOML.parse(read(f["Input/gpec_toml_raw"]))["Tearing"] + @test recorded["ModifiedRutherford"]["nw"] == 40 && recorded["inner_model"] == "SLAYER" + end + end + + # A forward solve without vac_flag has no Δ′ and allows no Riccati solve, so the deck stops at parse time. + mktempdir() do dir + deck_dir = _stage_deck(dir; ffs_overrides=Dict{String,Any}("vac_flag" => false), + extra_sections=Dict{String,Any}("Tearing" => Dict{String,Any}("ModifiedRutherford" => Dict{String,Any}()))) + @test_throws "needs Δ′" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) + end + # omega_ratio is the one optional (nothing-default) key: a deck value reaches the problem as a Float64. + @test GPEC.modified_rutherford_from_deck(Dict{String,Any}("Tearing" => Dict{String,Any}("ModifiedRutherford" => + Dict{String,Any}("omega_ratio" => 1))))[:omega_ratio] == 1.0 + # preset = "callen" expands under the section's own keys. + preset = GPEC.modified_rutherford_from_deck(Dict{String,Any}("Tearing" => Dict{String,Any}("ModifiedRutherford" => + Dict{String,Any}("preset" => "callen", "nw" => 50)))) + @test preset[:pol_model] === :callen && preset[:delta_prime_model] === :callen && preset[:nw] == 50 && !haskey(preset, :preset) + @test_throws "preset" GPEC.modified_rutherford_from_deck(Dict{String,Any}("Tearing" => Dict{String,Any}("ModifiedRutherford" => + Dict{String,Any}("preset" => "other")))) + end + + @testset "MRE Callen preset against Callen's closed form" begin + # The coil-driven DIII-D-like example: the MRE's Callen preset with the coil forcing and its own Riccati Δ′ solve, + # against the API and against Callen's closed-form critical field and cubic roots. + MRE = GeneralizedPerturbedEquilibrium.Tearing.ModifiedRutherford + RF = GeneralizedPerturbedEquilibrium.Tearing.ResonantForcing + ex_dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") + mktempdir() do dir + deck = TOML.parsefile(joinpath(ex_dir, "gpec.toml")) + foreach(k -> delete!(deck, k), ("LocalStability", "Ballooning", "ErrorFields", "KineticForces")) + deck["Run"]["output_file"] = joinpath(dir, "gpec.h5") + deck["Tearing"] = Dict{String,Any}("ModifiedRutherford" => Dict{String,Any}("preset" => "callen", "nw" => 400)) + for key in ("eq_filename", "kinetic_file") + deck["Equilibrium"][key] = abspath(joinpath(ex_dir, deck["Equilibrium"][key])) + end + open(io -> TOML.print(io, deck), joinpath(dir, "gpec.toml"), "w") + reference = GPEC.main([dir]) + @test reference.tearing === nothing + @test reference.modified_rutherford isa ModifiedRutherfordResult + + ric = reference.modified_rutherford.problem.ffs + @test ric isa ForceFreeStatesResult + @test ric !== reference.ffs + @test !isempty(ric.delta_prime.matrix) + mre = solve(ModifiedRutherfordProblem(ric; pe=reference.pe, MRE.CALLEN_PRESET..., nw=400)) + @test mre.F ≈ reference.modified_rutherford.F rtol = 1e-10 + + # Callen's closed form (Fortran GPEC 1.6.0) with his cylindrical Δ_c = −2m/ρ, in metres of ρ: the critical + # vacuum half-width ½ w_pol^(2/3) ((27/4)|Δ_c|)^(1/6)/Δ_RMP^(1/2), B_crit = B_vac (h_v,crit/h_v)², and the + # positive roots of Δ_c x³ + Δ_RMP (2h_v)² x − w_pol² = 0 (trigonometric form), all widths back in ψ_N. + equil = reference.pe.response.ffs.equil + plasma = GeneralizedPerturbedEquilibrium.InnerLayer.layer_plasma(equil; log=false) + function callen(pe) + rf = RF.resonant_forcing(pe) + map(eachindex(rf.rational_psi)) do k + rho, drho, d_rmp, hw_v = rf.rho[k], rf.drho[k], rf.delta_rmp[k], rf.hw_v[k] + w_pol = RF.callen_polarization_width(equil, plasma, rf.rational_psi[k], rf.rational_q[k])[2] + dc = -2 * rf.m[k] / rho + hw_v_crit = 0.5 * w_pol^(2 / 3) * ((27 / 4) * abs(dc))^(1 / 6) / sqrt(d_rmp) / drho + p, q = d_rmp * (2 * hw_v * drho)^2 / dc, -w_pol^2 / dc + R = sqrt(-p / 3) + arg = -q / (2 * R^3) + x = abs(arg) <= 1 ? sort(filter(>(0), [2R * cos((acos(arg) + 2π * j) / 3) for j in 0:2])) ./ drho : [NaN, NaN] + (; psi=rf.rational_psi[k], b_vac=rf.b_vac[k], d_rmp, w_v=2 * hw_v, b_crit=rf.b_vac[k] * (hw_v_crit / hw_v)^2, + w_min=x[1], w_sat=x[end]) + end + end + row(res, psi) = findfirst(p -> isapprox(p, psi; atol=1e-6), res.rational_psi) + # On the forward solve's own surfaces, with the coils scaled past the largest critical field so every surface + # has both roots, the preset meets the closed form to the root tolerance. + cyl0 = callen(reference.pe) + strong = perturbed_equilibrium(reference.pe.response, 4 * maximum(c.b_crit / c.b_vac for c in cyl0) * reference.pe.forcing; + coil_sets=reference.pe.coil_sets, reference.pe.regularization_width) + cyl = callen(strong) + pre = solve(ModifiedRutherfordProblem(reference.ffs; pe=strong, MRE.CALLEN_PRESET..., delta_prime_source=:cylindrical, nw=400)) + @test [c.b_crit for c in cyl] ≈ [c.b_crit for c in cyl0] rtol = 1e-8 + @test all(c -> c.w_sat > c.w_min, cyl) + compared = 0 + for c in cyl + s = row(pre, c.psi) + @test pre.delta_rmp[s] ≈ c.d_rmp rtol = 1e-12 + @test pre.w_v[s] ≈ c.w_v rtol = 1e-12 + @test pre.b_crit[s] ≈ c.b_crit rtol = 1e-8 + c.w_sat < pre.w[end] && c.w_min > pre.w_lo[s] || continue + @test pre.w_m[s] ≈ c.w_min rtol = 1e-8 + @test pre.w_sat[s] ≈ c.w_sat rtol = 1e-8 + compared += 1 + end + @test compared > 0 + # The computed Δ′ moves Δ_c, and so the critical field, off the cylindrical one. + @test !(mre.b_crit ≈ solve(ModifiedRutherfordProblem(ric; pe=reference.pe, MRE.CALLEN_PRESET..., delta_prime_source=:cylindrical, nw=400)).b_crit) + # The general MRE takes the same forcing: F gains exactly the forcing term, which scales as (w_v/w)². + kw = (; delta_prime_source=:cylindrical, saturation=true, nw=40) + forced = solve(ModifiedRutherfordProblem(ric; pe=reference.pe, kw...)) + free = solve(ModifiedRutherfordProblem(ric; pe=reference.pe, forcing=false, kw...)) + @test forced.F ≈ free.F .+ forced.forcing rtol = 1e-12 + @test forced.forcing .* forced.w' .^ 2 ≈ repeat(forced.forcing[:, 1] .* forced.w[1]^2, 1, 40) rtol = 1e-10 + @test all(>(0), forced.forcing) + @test isequal(forced.b_crit, free.b_crit) + @test all(isnan, solve(ModifiedRutherfordProblem(ric; kw...)).b_crit) + h5open(joinpath(dir, "gpec.h5")) do f + @test !haskey(f["Tearing"], "Callen") + @test isequal(read(f["Tearing/ModifiedRutherford/PerSurface/b_crit"]), reference.modified_rutherford.b_crit) + @test read(f["Tearing/ModifiedRutherford/Terms/forcing"]) ≈ reference.modified_rutherford.forcing + recorded = TOML.parse(read(f["Input/gpec_toml_raw"]))["Tearing"] + @test recorded["ModifiedRutherford"]["pol_model"] == "callen" + end + end + + # The retired [Tearing.Callen] section stops the parse as an unknown [Tearing] key. + mktempdir() do dir + deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("Tearing" => Dict{String,Any}("Callen" => Dict{String,Any}()))) + @test_throws "Callen" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) + end + end + @testset "kinetic_format" begin # A rates-v1 kinetic file read through the deck's [Equilibrium] keys and through the API keyword. mktempdir() do dir diff --git a/test/suites.jl b/test/suites.jl index dfe79f9e4..f584bf4cf 100644 --- a/test/suites.jl +++ b/test/suites.jl @@ -47,6 +47,7 @@ const TEST_SUITES = [ "runtests_dispersion_winding.jl" => "errorfields", "runtests_dispersion_ggj.jl" => "errorfields", "runtests_critical_resonant_field.jl" => "errorfields", + "runtests_modified_rutherford.jl" => "errorfields", "runtests_slayer_runner.jl" => "errorfields", "runtests_kinetic.jl" => "core", "runtests_kinetic_validity.jl" => "core",