Utilities Module

Utilities

Shared mathematical and computational utilities for GPEC modules.

This module provides common functionality used across multiple GPEC modules, including efficient Fourier transforms, numerical integration, and other mathematical utilities.

Submodules

• FourierTransforms: Efficient Fourier transforms with pre-computed basis functions
• PhysicalConstants: SI physical constants matching Fortran GPEC/SLAYER values
• NeoclassicalResistivity: Spitzer/Sauter/Redl resistivity closures shared by the GGJ and SLAYER inner-layer models
• HDF5Annotations: self-describing metadata (long_name/units/dims attributes and HDF5 Dimension Scales) for the gpec.h5 output

FourierCoefficients

Lightweight container for Fourier coefficients without spline interpolation. Use this when you only need to access coefficients at original grid points.

Fields

• xs::Vector{Float64}: Radial coordinates
• mmax::Int: Highest retained poloidal Fourier mode (modes 0:mmax inclusive)
• nqty::Int: Number of quantities
• cos_coeffs::Array{Float64,3}: Cosine coefficients (npsi × nmodes × nqty)
• sin_coeffs::Array{Float64,3}: Sine coefficients (npsi × nmodes × nqty)

FourierCoefficients(xs, ys, fs, mmax)

Compute Fourier coefficients via FFT without creating splines.

Arguments

• xs::Vector{Float64}: Radial coordinates
• ys::Vector{Float64}: Poloidal coordinates (periodic domain)
• fs::Array{Float64,3}: Function values (npsi × ntheta × nqty)
• mmax::Int: Highest poloidal Fourier mode to retain

KineticProfiles

Radial kinetic-profile container. All six profiles are 1D cubic splines of the normalized poloidal flux ψ ∈ [0, 1].

Construct via the keyword constructor KineticProfiles(; psi, n_e, T_e, T_i, omega, omega_e, omega_i) with matched-length vectors. The SLAYER runner builds this object from a standardized kinetic-profile file via Equilibrium.read_kinetic_file.

Evaluate all profiles at a given ψ via the call operator:

vals = kp(0.5)    # NamedTuple(n_e=..., T_e=..., ..., omega_i=...)

(kp::KineticProfiles)(psi::Real) -> NamedTuple

Evaluate all profiles at psi and return them as a NamedTuple with fields (n_e, T_e, T_i, omega, omega_e, omega_i).

empty_FourierCoefficients()

Create an empty FourierCoefficients for initialization purposes.

get_complex_coeff(fc::FourierCoefficients, ipsi, mode, qty) -> ComplexF64

Get normalized complex FFT coefficient at grid point.

The FourierCoefficients internally stores:

• cos_coeffs = 2 * real(FFT) / ntheta (for m > 0)
• sin_coeffs = -2 * imag(FFT) / ntheta (for m > 0)

This function returns the normalized FFT coefficient (FFT/ntheta):

• c[0] = cos_coeffs[0] (DC component)
• c[m] = cos_coeffs[m]/2 - i*sin_coeffs[m]/2 (for m > 0)

get_complex_coeffs!(out, fc::FourierCoefficients, ipsi, qty)

Fill vector with normalized complex FFT coefficients for modes 0:mmax.

powspace(x_min, x_max, npower, nx, spacing)

Power-law spaced grid with analytic derivative calculation. Creates a grid with more points concentrated near boundaries.

Arguments

• x_min, x_max: Grid boundaries
• npower::Int: Power law exponent (1=linear, higher=more concentration at edges)
• nx::Int: Number of grid points
• spacing::String: "lower" (concentrate at xmin), "upper" (at xmax), or "both"

Returns

• Matrix{Float64}: (2, nx) where row 1 = positions, row 2 = dr/d(norm) derivatives

FourierTransforms

Pre-computed complex Fourier basis and functor interface for θ ↔ mode transforms.

basis[ℓ, i] = exp(-i(m_ℓ θ_i - n ν_i)) with shape (mpert, mtheta). Forward: basis * data / mtheta; inverse: adjoint(basis) * modes (Fortran iscdftf/iscdftb; see docs/src/conventions.md).

FourierTransform

Struct with precomputed complex Fourier basis for repeated θ ↔ mode transforms.

Fields

• mtheta: poloidal grid size
• mpert: number of poloidal modes
• mlow: lowest poloidal mode number
• basis: $\exp(-i(m\theta - n\nu))$, size (mpert, mtheta)

(ft::FourierTransform)(data)

Forward transform from θ-space to mode space.

Arguments

• data: real or complex samples on the poloidal grid — Vector{mtheta} or Matrix{mtheta, :}

Returns

• Mode coefficients — Vector{mpert} or Matrix{mpert, :}

FourierTransform(mtheta, mpert, mlow; n=0, ν=zeros(mtheta))

Construct a transform with precomputed basis for contiguous modes mlow:(mlow+mpert-1).

Arguments

• mtheta: number of poloidal grid points
• mpert: number of poloidal modes
• mlow: lowest poloidal mode number

Keyword Arguments

• n: toroidal mode number (default 0)
• ν: toroidal angle offset on the poloidal grid, length mtheta

Returns

• FourierTransform ready for forward and inverse transforms

compute_fourier_coefficients(mtheta, m_modes, nzeta, n_modes; nfp=1)

Build 3D basis for every (m, n) in m_modes × n_modes (rows m-fast, n-slow).

Arguments

• mtheta: poloidal grid points per toroidal plane
• m_modes: poloidal mode numbers
• nzeta: full-torus toroidal grid points (must be divisible by nfp)
• n_modes: toroidal mode numbers

Keyword Arguments

• nfp: number of field periods; when > 1, emits one period only (nzeta ÷ nfp planes)

Returns

• Basis matrix, size (length(m_modes) * length(n_modes), mtheta * (nzeta ÷ nfp))

compute_fourier_coefficients(mtheta, m_modes, n, ν)

Build complex basis $\exp(-i(m\theta - n\nu))$ on the uniform poloidal grid.

Arguments

• mtheta: number of poloidal grid points
• m_modes: poloidal mode numbers (one row per mode)
• n: toroidal mode number
• ν: toroidal angle offset on the poloidal grid, length mtheta

Returns

• Basis matrix, size (length(m_modes), mtheta)

inverse(ft, modes)

Inverse transform from mode space to θ-space.

Arguments

• ft: precomputed transform
• modes: complex mode coefficients — Vector{mpert} or Matrix{mpert, :}

Returns

• θ-space data — Vector{mtheta} or Matrix{mtheta, :}

inverse_transform!(output, ft, modes)

In-place inverse transform (mode → θ-space) writing into output: output .= adjoint(basis) * modes.

Allocation-free wrapper over inverse for reuse in hot loops. output must have first dimension mtheta and match the column count of modes (whose first dimension is mpert).

transform!(output, ft, data)

In-place forward transform (θ-space → mode) writing into output: output .= (basis * data) / mtheta.

Allocation-free wrapper over the FourierTransform functor for reuse in hot loops. output must have first dimension mpert and match the column count of data (whose first dimension is mtheta).

PhysicalConstants

Shared physical constants used across GPEC modules. Values match the Fortran GPEC/SLAYER conventions (sglobal_mod) so numerical results can be directly compared.

All quantities in SI units.

NeoclassicalResistivity

Spitzer + Sauter / Redl neoclassical resistivity closures, shared between the GGJ and SLAYER inner-layer models so both see identical plasma-input physics when the same NeoResistivityModel is selected.

Exports

Abstract tag for a neoclassical-resistivity closure.

Redl et al. 2021 F_33 neoclassical correction (Eqs. 17-18). Improved high-collisionality fit vs SauterNeoModel.

Sauter, Angioni & Lin-Liu 1999 F_33 neoclassical correction (Eqs. 13a,b).

Spitzer-Härm parallel resistivity 1/σ∥ as used by Fitzpatrick's TJ (LayerParameters.tex Eqs. 7-8) and legacy SLAYER for τR. No trapped-particle correction. Pair with lnLambda_form=:wesson to reproduce legacy SLAYER behaviour bit-identically.

Plain Spitzer resistivity — no trapped-particle correction.

coulomb_log_e(n_e, T_e; form=:nrl) -> Float64

Electron Coulomb logarithm. n_e in m⁻³, T_e in eV.

form=:nrl (default) uses the NRL Plasma Formulary 2009 expression, which OpenFUSIONToolkit's bootstrap.py also selects as the "more accurate" option. form=:sauter uses the simpler Sauter 1999 Eq. 18d form.

eta_neoclassical(model, n_e, T_e, Z_eff, f_t, nu_e_star;
                 lnLamb=nothing) -> Float64

Neoclassical resistivity η [Ω·m] under the chosen closure.

• SpitzerModel() – returns eta_spitzer(n_e, T_e, Z_eff; lnLamb) unchanged; f_t and nu_e_star are ignored.
• SpitzerHarmModel() – returns eta_spitzer_harm(n_e, T_e, Z_eff; lnLamb) (Fitzpatrick/TJ legacy); f_t and nu_e_star are ignored.
• SauterNeoModel() – Sauter 1999 Eq. 13: η = ηSp / F33(Sauter).
• RedlNeoModel() – Redl 2021 Eq. 17: η = ηSp / F33(Redl).

Note that σneo = σSp · F33, so ηneo = ηSp / F33. For a banana-regime plasma with ft ≈ 0.5 and ν*e ≪ 1, F33 ≈ 0.4–0.5, so ηneo is a factor of ~2 larger than η_Sp — this is the standard H-mode tearing correction.

eta_spitzer(n_e, T_e, Z_eff; lnLamb=nothing) -> Float64

Spitzer resistivity in Ω·m, using the Sauter 1999 Eq. 18a form

σ_Sp = 1.9012e4 · T_e^1.5 / (Z_eff · N(Z_eff) · lnΛ_e)
N(Z) = 0.58 + 0.74 / (0.76 + Z)
η_Sp = 1 / σ_Sp

n_e [m⁻³], T_e [eV]. lnLamb defaults to coulomb_log_e(n_e, T_e) (NRL).

eta_spitzer_harm(n_e, T_e, Z_eff; lnLamb=nothing) -> Float64

Spitzer-Härm parallel resistivity η∥ = 1/σ∥ in Ω·m, per Fitzpatrick TJ LayerParameters.tex Eq. 8:

σ_∥ = (√2 + 13 Z_eff/4) / (Z_eff (√2 + Z_eff)) · n_e e² τ_ee / m_e

This is the resistivity entering the legacy SLAYER / TJ resistive diffusion time τR = μ₀ rs² σ∥ (LayerParameters.tex Eq. 17). Agrees with `etaspitzer` to the fit accuracy of the two formulas (~1% at Z=1).

nu_star_e(n_e, T_e, R_major, eps, q, Z_eff; lnLamb=nothing) -> Float64

Electron collisionality ν*_e per Sauter 1999 Eq. 18b:

ν*_e = 6.921e-18 · |q| · R · n_e · Z_eff · lnΛ_e / (T_e² · ε^1.5)

n_e [m⁻³], T_e [eV], R_major [m]. Matches OFT bootstrap.py:640 and OMFIT utils_fusion.py:1278.

tau_ee_spitzer_harm(n_e, T_e; lnLamb=nothing) -> Float64

Electron-electron collision time per Fitzpatrick TJ LayerParameters.tex Eq. 7:

τ_ee = 6√2 π^1.5 ε₀² √m_e T_e^1.5 / (lnΛ e^2.5 n_e)

n_e in m⁻³, T_e in eV (the e^2.5 denominator absorbs the eV→J conversion). lnLamb defaults to coulomb_log_e(n_e, T_e) (NRL).

trapped_fraction(avg_B, avg_Bsq, B_min, B_max) -> Float64

Lin-Liu & Miller 1995, Phys. Plasmas 2, 1666, Eq. (6):

f_t = 1 − ⟨B⟩² / ⟨B²⟩ · (1 − √(1 − h) · (1 + h/2)),   h = B_min / B_max

Equivalent to the OMFIT f_t / f_c pair at full geometric accuracy (uses both the average-B ratio and the min/max extremes). Arguments are flux-surface averages computed from the θ-loop in the equilibrium.

trapped_fraction_eps(eps) -> Float64

Simple ε-only trapped-fraction approximation (OMFIT f_t):

f_c ≈ (1 − ε)² / (√(1 − ε²) · (1 + 1.46·√ε + 0.2·ε))
f_t = 1 − f_c

Used as a fallback when the full (⟨B⟩, ⟨B²⟩, Bmin, Bmax) moments are unavailable — e.g. when feeding SLAYER directly from minor-radius geometry without having evaluated ResistGeometry first.

HDF5Annotations

Self-describing metadata for gpec.h5 (the contract in docs/development/hdf5-conventions.md): every dataset carries a long_name and units attribute, array datasets carry a dims axis-name attribute, and coordinate datasets are marked as HDF5 Dimension Scales (netCDF-4 coordinate variables) attached to the arrays that share their axis, so h5py/xarray/HDFView read the file unaided.

Writers stay table-driven: each writer keeps a table of path => (; long_name, units, dims) entries next to it and calls annotate! once after its datasets are written. Paths absent from the file are skipped silently (many writes are conditional).

annotate!(parent, table)

Apply a metadata table to datasets under parent (an open HDF5.File or group). table iterates path => meta pairs where meta is a NamedTuple with fields long_name (required), units (default "1" = dimensionless), and optionally:

• dims — tuple of axis names in Julia (column-major) order, axis 1 first, stored as the greppable string attribute dims = "(psi, m)".
• scale — mark this dataset as an HDF5 Dimension Scale with the given name (netCDF-4 coordinate variable).
• attach — tuple of axis => scale_path pairs attaching declared scales to this dataset's Julia axes; scale_path is relative to parent and its dimension label comes from that entry's scale name (falling back to the path basename).

Scales are marked in a first pass so attachments within the same table resolve. Missing paths are skipped throughout.

attach_scale!(parent, path, julia_axis, scale_path, label)

Attach the Dimension Scale at scale_path to Julia axis julia_axis (axis 1 first) of the dataset at path, and label that dimension. The H5DS C API indexes file (row-major) dimensions, so Julia axis k of an N-d dataset is C index N - k. No-op when either path is absent or the axis lengths disagree.

make_scale!(parent, path, name)

Mark the dataset at path as an HDF5 Dimension Scale named name. No-op when the path is absent.

write_root_attrs!(file; title)

Stamp the file-level contract: schema_version, Conventions, references, title (run description), date_created (ISO 8601 UTC). The code version lives in Info/git_version.

write_imas(dd, result)

Write GPEC stability results into dd.mhd_linear. Creates one toroidal_mode entry per requested toroidal mode number, storing the least-stable (minimum real part) energy_perturbed for that n_tor. For multi-n runs the eigenvalue array et is sorted by stability across all n-blocks; n_tor_idx[i] identifies which n-block eigenvalue i belongs to, so each n_tor receives the correct least-stable δW regardless of how modes are interleaved in et.

The result argument is the named tuple returned by main; its ffs field carries the force-free-states result the energies are read from.