API Reference

calculate_energy_transfer(method, velocity_data, coords_or_ks; kwargs...)

Unified entry point for all energy transfer computations.

Arguments

• method::AbstractEnergyTransferMethod: Which method to use:
  • SpectralFluxMethod(binning) — spectral flux Π(K)
  • CoarseGrainingFluxMethod(filter, ℓ) — coarse-graining flux Π_ℓ(x)
  • ShellToShellTransferMethod(binning) — shell-to-shell T(n,m)
• velocity_data: For spectral methods, a complex array of size (ns..., D) containing Fourier coefficients; for coarse-graining, a tuple of D real physical-space arrays.
• coords_or_ks: For spectral methods, a tuple of 1D wavenumber vectors; for coarse-graining, a tuple of 1D coordinate vectors.

Returns

Method-specific result container: SpectralFluxResult, CoarseGrainingFluxResult, or ShellToShellResult.

Examples

using FlowInvariantTransfer, FFTW

# Spectral flux on a 32×32 periodic domain
N = 32; L = 2π
x = range(0.0, L; length=N+1)[1:N]
y = range(0.0, L; length=N+1)[1:N]
u = [cos(x) for x in x, y in y]
v = [sin(y) for x in x, y in y]
û = cat(FFTW.fft(u), FFTW.fft(v); dims=3) ./ N^2  # (N,N,2)
ks = FlowInvariantTransfer.Utils.wavenumber_grid((N,N), (L,L))

result = calculate_energy_transfer(SpectralFluxMethod(LinearBinning(2π/L)), û, ks)

calculate_spectral_flux(velocity_hat, ks; binning, dealiasing=OrszagTwoThirds()) -> SpectralFluxResult

Compute the spectral energy transfer spectrum T(k) and the cumulative energy flux Π(K) from Fourier-space velocity data.

Arguments

• velocity_hat: Complex array of size (ns..., D) — Fourier coefficients of the D velocity components on an N-point periodic grid.
• ks: Tuple of 1D physical-wavenumber vectors (matching FFTW fftfreq convention).

Keyword Arguments

• binning::AbstractShellBinning: Shell binning strategy; default LinearBinning(1.0).
• dealiasing::AbstractDealiasing=OrszagTwoThirds(): Apply 2/3 dealiasing rule when computing (u·∇)u.
• spectral::SpectralBackends.AbstractSpectralBackend: transform backend — SpectralBackends.DirectSumSpectralBackend() (default, no deps) or SpectralBackends.FFTSpectralBackend() (requires FFTW extension).
• execution::ComputationalBackends.AbstractExecutionBackend=ComputationalBackends.SerialBackend(): how the transfer-density write and the mode→shell reduction are parallelised (orthogonal to spectral, which threads the FFT itself):
  • ComputationalBackends.SerialBackend() (default) — host scalar reduction.
  • ComputationalBackends.ThreadedBackend() (requires using OhMyThreads) — the mode→shell scatter is split into chunks with per-chunk partial shell sums (race-free), an O(Nᴰ) parallel pass.
  • ComputationalBackends.DistributedBackend() (requires using Distributed) — the scatter is partitioned across worker processes and +-reduced. Note: this distributes only the reduction, not the (dominant) nonlinear-term FFT; to distribute a grid too large for one node, use pencil_spectral_flux (PencilFFTs).
  • ComputationalBackends.GPUBackend(dev) (requires using KernelAbstractions) — the transfer density is written by a device kernel and reduced per shell on-device (no host round-trips), so a device-resident field stays on the GPU.

Returns

SpectralFluxResult with fields:

• k_shells: Representative wavenumber per shell.
• transfer_spectrum: T(k_n) — energy input to shell n per unit time.
• flux: Π(Kn) — cumulative upscale flux (energy transferred to k > Kn).

Physics

T(kn) = Σ{|k| ∈ shelln} Re{ û*(k) · N̂(k) } Π(Kn) = −Σ{m ≤ n} T(km)

Positive Π: forward (downscale) cascade; negative Π: inverse (upscale) cascade.

References

• Verma et al. (2002) [arXiv:nlin/0204027]
• Alexakis, Mininni & Pouquet (2005)

calculate_spectral_flux!(result, ws, velocity_hat, ks, shell_idx;
    dealiasing, invariant, spectral, execution, advecting_hat)

In-place version of calculate_spectral_flux. Writes into result using preallocated buffers from ws and a precomputed shell_idx array (from assign_shells). Zero heap allocations in the serial hot path. execution selects the mode→shell reduction backend (see calculate_spectral_flux).

calculate_scalar_flux(velocity_hat, scalar_hat, ks; binning, dealiasing=OrszagTwoThirds(), spectral) -> SpectralFluxResult

Compute the passive-scalar variance transfer spectrum T_θ(k) and flux Π_θ(K), for a scalar θ advected by the velocity u (∂_tθ + (u·∇)θ = κ∇²θ):

T_θ(k_n) = Σ_{|k|∈shell_n} Re{ θ̂*(k) N̂_θ(k) },   N̂_θ = FFT[(u·∇)θ],   Π_θ(K) = Σ_{k≤K} T_θ(k).

Scalar variance is conserved for incompressible u (Σ_k T_θ ≈ 0) and cascades forward in any dimension (Obukhov–Corrsin). Thin wrapper over calculate_spectral_flux with invariant = PassiveScalar() and advecting_hat = velocity_hat.

Arguments

• velocity_hat: complex (ns..., D) velocity Fourier coefficients (the advecting field).
• scalar_hat: complex scalar field, either (ns...) or (ns..., 1).
• ks: tuple of nd 1D wavenumber vectors.

calculate_shell_to_shell_transfer(velocity_hat, ks;
    binning, dealiasing=OrszagTwoThirds(), verify_antisymmetry=true,
    spectral=SpectralBackends.DirectSumSpectralBackend(), execution=ComputationalBackends.SerialBackend())
    -> ShellToShellResult

Compute the directed shell-to-shell kinetic energy transfer matrix T(n,m).

Arguments

• velocity_hat: Complex array of size (ns..., D) — Fourier coefficients of the D velocity components on a periodic uniform grid.
• ks: Tuple of D 1D physical-wavenumber vectors.

Keyword Arguments

• binning::AbstractShellBinning: Shell binning; default LinearBinning(1.0).
• dealiasing::AbstractDealiasing=OrszagTwoThirds(): Apply 2/3 rule dealiasing.
• verify_antisymmetry::Bool=true: Compute max|T(n,m)+T(m,n)| and store in result.
• spectral::SpectralBackends.AbstractSpectralBackend: transform — SpectralBackends.DirectSumSpectralBackend() (default) or SpectralBackends.FFTSpectralBackend() (FFTW).
• execution::ComputationalBackends.AbstractExecutionBackend: outer (mediator-loop) parallelism — ComputationalBackends.SerialBackend() (default), ComputationalBackends.ThreadedBackend() (OhMyThreads), ComputationalBackends.DistributedBackend(), or ComputationalBackends.GPUBackend(...).

Returns

ShellToShellResult with:

• transfer_matrix[n,m]: Energy transferred from shell m to shell n.
• net_transfer[n] = Σ_m T(n,m): net energy gain of shell n.
• max_antisymmetry_error: validation diagnostic.

Algorithm (Verma 2002 formulation)

For each pair of receiver shell n and mediator shell m: T(n,m) = Σ{k∈Sn} Re{ ûn*(k) · N̂m(k) } where N̂m(k) = FFT[(um · ∇)u], um = IFFT(û · χm).

This formulation uses the mediator velocity restricted to shell m, so: T(n,m) + T(m,n) = 0 exactly (antisymmetry).

Cost

O(Nshells² · N^D log N^D) with FFTW; O(Nshells² · N^{2D}) direct-sum.

References

• Verma et al. (2002), arXiv:nlin/0204027
• Alexakis, Mininni & Pouquet (2005)

calculate_shell_to_shell_transfer!(result, ws, velocity_hat, ks; kwargs...)

In-place version. Writes into result using preallocated buffers from ws. Zero heap allocations in the hot path.

calculate_scalar_shell_to_shell_transfer(velocity_hat, scalar_hat, ks; kwargs...) -> ShellToShellResult

Shell-to-shell transfer of passive-scalar variance T_θ(n,m): the rate at which scalar variance is transferred from scalar-shell m to scalar-shell n, mediated by the velocity:

T_θ(n,m) = Σ_{k∈S_n} Re{ θ̂*(k) · 𝒩̂_m(k) },   𝒩̂_m = FFT[(u·∇)θ_m],   θ_m = θ̂·χ_m.

The scalar field is band-filtered and carried; the velocity advects it. Thin wrapper over calculate_shell_to_shell_transfer with invariant = PassiveScalar() and advecting_hat = velocity_hat. scalar_hat may be (ns...) or (ns..., 1). As with energy, T_θ(n,m) is antisymmetric for incompressible u and reduces to T_θ(k) over mediators.

calculate_mode_to_mode_transfer(velocity_hat, ks; invariant=KineticEnergy(), dealiasing=OrszagTwoThirds(),
                                spectral=SpectralBackends.DirectSumSpectralBackend(), max_scales=1024, force=false)
    -> ModeToModeTriadResult

Fully mode-resolved triad transfer S(k|p) — the rate at which the chosen quadratic invariant is delivered to receiver scale k from giver scale p (mediated by q = k−p), the finest object in the reduction hierarchy S(k|p) → T(K,Q) (shell-to-shell) → T(k), Π(K) (spectral flux).

It is built from the validated pseudospectral nonlinear term — for each giver scale p, N̂_p = (u·∇)u_p (the full velocity advecting the single-scale field u_p), and

S(k|p) = Re{ û*(k) · N̂_p(k) }   (generalised per invariant via `transfer_density!`).

This construction is exact and inherits the right structural properties (verified by tests):

• reduces: Σ_p S(k|p) = T(k) = the spectral transfer (calculate_spectral_flux),
• antisymmetric: S(k|p) + S(p|k) = 0 (incompressible, since ∫(u·∇)(u_p·u_k)=0),
• conserves: Σ_k Σ_p S(k|p) = 0.

Cost & memory

Resolving every receiver/giver pair is O(N_scales) nonlinear-term evaluations → O(N_scales · Nᴰ log N) time with SpectralBackends.FFTSpectralBackend (strongly recommended) and an O(N_scales²) result tensor. For the aggregates prefer the cheaper, coarser diagnostics: calculate_spectral_flux (T(k), Π) or calculate_shell_to_shell_transfer (T(n,m)). A guard errors when N_scales = prod(size grid) > max_scales; pass force=true to override.

Keyword arguments

• invariant::AbstractInvariant: which quadratic invariant (default KineticEnergy()).
• dealiasing::AbstractDealiasing=OrszagTwoThirds(): 2/3-rule dealiasing of the nonlinear term.
• spectral::SpectralBackends.AbstractSpectralBackend: transform (SpectralBackends.DirectSumSpectralBackend() default, SpectralBackends.FFTSpectralBackend() fast).
• max_scales::Int=1024, force::Bool=false: resolved-tensor size guard.

Returns

ModeToModeTriadResult with net_transfer (T(k), shape ns) and transfer (the resolved S(k|p), shape (ns..., ns...)).

calculate_scalar_mode_to_mode_transfer(velocity_hat, scalar_hat, ks; kwargs...) -> ModeToModeTriadResult

Fully mode-resolved passive-scalar variance transfer S_θ(k|p) — variance delivered to scalar scale k from scalar scale p (mediated by the velocity, q = k−p). Thin wrapper over calculate_mode_to_mode_transfer with invariant = PassiveScalar() and advecting_hat = velocity_hat; the scalar may be (ns...) or (ns..., 1).

calculate_band_to_band_transfer(velocity_hat, ks; bands::SmoothBands, dealiasing=OrszagTwoThirds(),
    invariant=KineticEnergy(), spectral=SpectralBackends.DirectSumSpectralBackend(), advecting_hat=velocity_hat,
    geometry=IsotropicShells())
    -> (centers, transfer_matrix, net_transfer, max_antisymmetry_error)

Smooth band-to-band transfer T(n,m) of a quadratic invariant between the graded spectral bands (Eyink & Aluie 2009) — the smooth-filter analogue of calculate_shell_to_shell_transfer. For incompressible flow T is antisymmetric (T(n,m) = −T(m,n)), conserves (Σ T = 0), and Σ_m T(n,m) is the band-summed transfer spectrum. Accepts an advecting_hat (primary field can be a passive scalar advected by the velocity) and an anisotropic geometry. See calculate_band_to_band_transfer! for the in-place, workspace-reusing variant.

calculate_band_to_band_transfer!(T, net, bws::BandTransferWorkspace, velocity_hat, ks; kwargs...)

In-place smooth band-to-band transfer: writes the nb×nb matrix T and the nb-vector net, reusing bws (0 alloc beyond those). Returns the same (centers, transfer_matrix, net_transfer, max_antisymmetry_error) named tuple as the allocating version.

calculate_partial_fluxes(velocity_hat, ks; decomposition=HelicalDecomposition(), binning,
    dealiasing=OrszagTwoThirds(), spectral=SpectralBackends.DirectSumSpectralBackend(), geometry=IsotropicShells())
    -> (channels::Dict{NTuple{3,Symbol},SpectralFluxResult}, total::SpectralFluxResult, k_shells)

Decompose the kinetic-energy flux into per-component partial fluxes Π^{s_k s_p s_q}(K), where each of the three fields in a triad interaction is one component of a velocity decomposition u = Σ_s u_s (e.g. ±-helical via HelicalDecomposition, or rotational/divergent via HelmholtzDecomposition):

T^{s_k s_p s_q}(k) = Re{ û_{s_k}*(k) · [ (u_{s_p}·∇) u_{s_q} ](k) },   Π = Σ_{k≤K} T.

With an n-component decomposition this gives n³ channels that sum to the full energy flux. For helical components the homochiral channels (s_k=s_p=s_q) drive the inverse cascade and the heterochiral ones the forward cascade (Biferale–Musacchio–Toschi 2012); for the Helmholtz split the off-diagonal channels are the rotational↔divergent cross-flux (zero for incompressible flow, since u_div = 0). channels is keyed by the component-name triple (s_k, s_p, s_q). Built from the decomposition + generalized nonlinear term, so it inherits all backends/dealiasing.

calculate_helical_partial_fluxes(velocity_hat, ks; kwargs...)

The eight helical partial energy fluxes Π^{s_k s_p s_q}(K), s ∈ {positive, negative} — calculate_partial_fluxes with decomposition = HelicalDecomposition() (3D only). Homochiral channels drive the inverse cascade, heterochiral the forward (Waleffe 1992; Biferale–Musacchio–Toschi 2012; Alexakis 2017).

calculate_compressible_flux(velocity_hat, density_hat, ks; binning, pressure_hat=nothing,
    decompose=true, geometry=IsotropicShells(), spectral=SpectralBackends.FFTSpectralBackend()) -> CompressibleFluxResult

Compressible kinetic-energy spectral transfer T_u(k) and cumulative flux Π(K) (Singh–Tiwari– Sharma–Verma 2025): momentum v = ρu, E_u(k) = ½Re[v·u*]; the nonlinear transfer conserves total KE (Σ_k T_u ≈ 0), and the KE↔internal-energy pressure-dilatation is returned separately when pressure_hat is supplied. decompose=true also returns the Helmholtz rotational/compressive flux channels. In the incompressible limit T_u reduces to −ρ × the incompressible transfer spectrum. This allocates a CompressibleWorkspace and delegates to calculate_compressible_flux!; build the workspace once and use the in-place form to loop over snapshots allocation-free.

calculate_compressible_flux!(ws::CompressibleWorkspace, velocity_hat, density_hat, ks; kwargs...)

In-place momentum-weighted compressible transfer reusing ws (FFT plans + all field scratch). Returns a fresh CompressibleFluxResult whose arrays are the small per-shell vectors; the O(field) intermediates live in ws and are reused across calls (build the workspace once, loop over snapshots).

calculate_coarse_graining_flux(velocity_fields, coords_vecs, ℓ, filter;
                               decomposition=NoDecomposition(), return_diagnostics=false, kwargs...)
                               -> CoarseGrainingFluxResult

Compute the pointwise cross-scale kinetic energy flux Πℓ(x) = −τ̄ᵢⱼ S̄ᵢⱼ at filter scale ℓ. Delegates entirely to CoarseGrainingEnergyFluxes.jl (`computeΠ!,filter_field!`).

Requires CoarseGrainingEnergyFluxes to be loaded:

using CoarseGrainingEnergyFluxes
result = calculate_coarse_graining_flux((u, v), (x, y), ℓ, GaussianFilter())

Arguments

• velocity_fields: Tuple of D real arrays (u, v[, w]) — velocity components.
• coords_vecs: Tuple of D coordinate vectors (x, y[, z]).
• ℓ::Real: Filter length scale (same units as coordinates).
• filter::AbstractFilter: GaussianFilter(), SharpSpectralFilter(), or TopHatFilter().

Keyword Arguments

• decomposition::AbstractFieldDecomposition: NoDecomposition() (default), HelmholtzDecomposition(), RotationalDecomposition(), or DivergentDecomposition().
• return_diagnostics::Bool=false: Also return τ̄ᵢⱼ and S̄ᵢⱼ fields.
• mask::Union{Nothing,AbstractMatrix{Bool}}=nothing: boolean point mask (true = point included in the computation, false = excluded). If nothing, all points are included.
• Any additional kwargs are forwarded to CoarseGrainingEnergyFluxes.compute_Π!.

Returns

CoarseGrainingFluxResult or NamedTuple of CoarseGrainingFluxResult depending on decomposition.

nufft_coarse_graining_flux(velocity_fields, scatter_coords, ℓ, filter, ms; kwargs...)

Coarse-graining energy flux Π_ℓ(x) at scattered (non-uniform) Cartesian points via FINUFFT. Requires using FINUFFT (the extension supplies the method).

calculate_spherical_transfer(...)

Extension entry point for the spherical spectral transfer. Implemented in the FastSphericalHarmonics (regular grid) and NUFSHT (scattered) extensions; loading one of those packages provides the method. Prefer the unified calculate_energy_transfer(SphericalTransferMethod(), …).

calculate_spherical_transfer!(ws::SphericalTransferWorkspace, vorticity; kwargs...) -> SphericalTransferResult

In-place spherical spectral transfer reusing ws. Implemented in the FastSphericalHarmonics / NUFSHT extensions. Reuses the workspace's spectral work arrays, per-mode reduction buffers, and the result vectors, so a repeat call on the same grid re-allocates only what the underlying transform library itself allocates (FastSphericalHarmonics has no in-place transform API — that portion is an irreducible floor). Build the workspace once with SphericalTransferWorkspace.

calculate_divergent_spherical_transfer(...)

Extension entry point for the divergent horizontal-KE spherical spectral transfer. Implemented in the FastSphericalHarmonics (regular grid) and NUFSHT (scattered) extensions; loading one provides the method. Prefer the unified calculate_energy_transfer(DivergentSphericalTransferMethod(), …).

calculate_divergent_spherical_transfer!(ws, u_θ, u_φ; kwargs...) -> DivergentSphericalTransferResult

In-place divergent spherical KE transfer reusing ws. Implemented in the FastSphericalHarmonics / NUFSHT extensions. Build the workspace once with DivergentSphericalTransferWorkspace (regular grid) or ScatteredDivergentSphericalTransferWorkspace (scattered points).

mpi_batch_map(f, items; comm=MPI.COMM_WORLD, reduce=:gather, root=0)

Distribute an embarrassingly-parallel batch of independent inputs across MPI ranks: each rank applies f to a round-robin subset of items (e.g. snapshots of a time series), then the per-item outputs are combined. With reduce=:gather (default) the results are collated into one Vector in the original order of items, returned on every rank; reduce=:sum/:mean returns the element-wise reduction (the outputs of f must support +, and / for :mean); a callable reduce is applied as a binary combiner. This is the "batch axis" of distribution — orthogonal to the pencil axis (pencil_spectral_flux), which splits a single grid.

Requires using MPI to load the extension.

pencil_spectral_flux(u_phys, ks; comm=MPI.COMM_WORLD, binning, dealiasing, invariant) -> (centers, transfer_spectrum, flux)

Distributed spectral transfer/flux for a single grid split across MPI ranks along the pencil axis: u_phys is the physical-space velocity held as a PencilArray (each rank owns a pencil of the global grid). The pseudospectral nonlinear term is evaluated with a transpose-based distributed FFT (PencilFFTs), the transfer density is shell-binned locally, and the per-shell spectrum is MPI.Allreduced to a global result identical on every rank (matching the serial calculate_spectral_flux on the same field). Use this when one snapshot's grid is too large for a single node; for many independent snapshots use mpi_batch_map instead.

Requires using MPI, PencilFFTs, PencilArrays to load the extension.

build_pencil_plan(ns, comm=MPI.COMM_WORLD; T=Float64) -> PencilFFTPlan

Convenience constructor for the distributed complex-to-complex FFT plan used by pencil_spectral_flux, with an auto-balanced MPI process grid. Requires using MPI, PencilFFTs, PencilArrays.

triadic_orthogonal_decomposition(X; kwargs...)
    -> TriadicOrthogonalDecompositionResult

Compute the Triadic Orthogonal Decomposition of data array X.

TOD decomposes triadic (three-wave) nonlinear interactions in time-series data, identifying coherent flow structures that optimally capture spectral momentum transfer. It produces a mode bispectrum (singular values quantifying coupling strength per frequency triad), convective/recipient modes, and a modal energy budget.

Arguments

• X::AbstractArray: Data of size (nt, nvar, spatial_dims...). First dimension is time, second is variable indices, remaining are spatial.

Keyword Arguments

• window: Temporal window. Vector → used directly. Integer → Hamming of that length. nothing → auto Hamming (length = 2^floor(log2(nt/5)), capped at 256).
• weight: Spatial inner-product weight (same spatial dims as X). nothing → uniform.
• noverlap: Block overlap in snapshots. nothing → 50% of window length.
• dt: Time step between snapshots. nothing → 1/nDFT (frequency index output).
• Q: Quadratic nonlinearity function (q1, q2) -> product. Default: element-wise product with permutation matching the MATLAB reference.
• LHS: Left-hand side operator q -> Lq. Default: identity.
• nmode: Modes per triad to store. nothing → nBlks.
• nfreq: Restrict to |l|, |k|, |n| ≤ nfreq. nothing → all.
• isreal_data: Whether data is real (restricts bispectrum to fn ≥ 0). nothing → auto.
• mean_type: :zero (default), :blockwise, or an array (long-time mean to subtract).
• return_coefficients::Bool=false: Also compute expansion coefficients.
• return_auxiliary_modes::Bool=false: Also compute donor/catalyst modes.
• spectral::SpectralBackends.AbstractSpectralBackend=SpectralBackends.DirectSumSpectralBackend(): temporal-DFT transform. SpectralBackends.FFTSpectralBackend() uses FFTW (much faster; requires using FFTW).
• execution::ComputationalBackends.AbstractExecutionBackend=ComputationalBackends.SerialBackend(): triad-loop parallelism. ComputationalBackends.ThreadedBackend() parallelises the triad loop (requires OhMyThreads).

Returns

TriadicOrthogonalDecompositionResult containing:

• frequencies: Frequency vector.
• mode_bispectrum: Singular values per triad per mode.
• modes: Dict of convective/recipient mode pairs.
• modal_energy_budget: Energy transfer per triad per mode.
• expansion_coefficients: Expansion coefficients (or nothing).
• auxiliary_modes: Donor/catalyst modes (or nothing).

References

• Yeung, Chu & Schmidt (2026), J. Fluid Mech. 1031, A34. DOI 10.1017/jfm.2026.11183

hamming_window(N, ::Type{T} = Float64) -> Vector{T}

Standard Hamming window of length N: w[n] = 0.54 − 0.46·cos(2πn/(N−1)), in element type T.

hann_window(N, ::Type{T} = Float64) -> Vector{T}

Hann (raised-cosine) window of length N: w[n] = ½(1 − cos(2πn/(N−1))), in element type T. Tapers to zero at both ends — lower spectral leakage than Hamming. Pass to triadic_orthogonal_decomposition via window.

tukey_window(N, ::Type{T} = Float64; α=0.5) -> Vector{T}

Tukey (tapered-cosine) window of length N in element type T: a flat middle with cosine tapers over a fraction α of the length at each end. α = 0 is rectangular (no taper), α = 1 is the Hann window; intermediate α trades main-lobe width against leakage.

compute_nonlinear_term(advected_hat, ks; dealiasing=OrszagTwoThirds(),
                       spectral=SpectralBackends.DirectSumSpectralBackend(), advecting_hat=advected_hat)

Compute the pseudospectral nonlinear term 𝒩̂ᵢ(k) = F̂[(uⱼ ∂fᵢ/∂xⱼ)] — the advection of an M-component field f (advected_hat) by a velocity u (advecting_hat). For the momentum self-advection term pass the velocity as both (the default), giving N̂ᵢ = F̂[(u·∇)uᵢ].

Arguments

• advected_hat: Array of size (ns..., M) — Fourier coefficients of the advected field f (M = D for momentum, M = 1 for a passive scalar / vector potential).
• ks: Tuple of 1D wavenumber vectors (length nd), one per spatial dimension.

Keyword Arguments

• dealiasing::AbstractDealiasing=OrszagTwoThirds(): dealiasing strategy (NoDealiasing / OrszagTwoThirds / PaddedThreeHalves).
• spectral::SpectralBackends.AbstractSpectralBackend: SpectralBackends.DirectSumSpectralBackend() (default, no deps) or SpectralBackends.FFTSpectralBackend() (requires the FFTW extension) for the O(N log N) path.
• advecting_hat: the advecting velocity u (shape (ns..., D), D ≥ nd); defaults to advected_hat (self-advection). Only the nd spatial components participate in (u·∇).

Returns

Array of size (ns..., M) containing 𝒩̂ᵢ(k).

compute_nonlinear_term!(ws, advected_hat, ks; dealiasing=OrszagTwoThirds(),
                        spectral=SpectralBackends.DirectSumSpectralBackend(), advecting_hat=advected_hat)

In-place version of compute_nonlinear_term. Writes result into ws.N̂. Pass a NonlinearTermWorkspace (sized for advected_hat) to avoid any allocations in the hot path.

transfer_density(invariant, velocity_hat, N̂, ks) -> Array

Allocating version of transfer_density!: returns a real array of shape ns with the per-mode transfer density.

transfer_density!(t, invariant, velocity_hat, N̂, ks) -> t

Write the real per-mode transfer density for invariant into t (shape ns), given Fourier-space velocity velocity_hat and nonlinear term N̂ (both shape (ns..., D)) and wavenumber vectors ks (length D). No allocations.

decompose_field(decomp::AbstractFieldDecomposition, fields::Tuple, coords::Tuple; kwargs...)

Decompose a physical-space velocity field fields (e.g. (u, v)) using the coordinate vectors coords and the specified decomposition strategy.

decompose_field(decomp::AbstractFieldDecomposition, velocity_hat::AbstractArray{<:Complex}, ks)

Decompose a spectral-space velocity field velocity_hat along the wavenumbers ks.

helmholtz_project_spectral!(args...; kwargs...)

In-place spectral Helmholtz projection of a velocity field into its rotational/divergent parts. Provided by the HelmholtzDecomposition.jl extension; this core stub errors until it is loaded.

AbstractEnergyTransferMethod

Abstract supertype for all energy transfer computation methods. Concrete subtypes dispatch calculate_energy_transfer to specific algorithms.

SpectralFluxMethod{B<:AbstractShellBinning} <: AbstractEnergyTransferMethod

Compute the spectral energy flux Π(K) and transfer spectrum T(k) using the pseudospectral method on a periodic uniform grid.

Fields

• binning::B: Shell binning strategy for grouping wavenumbers.

Notes

Requires Fourier-space velocity data on a uniform periodic grid. When FFTW is loaded, all transforms run in O(N log N); without it, falls back to an O(N²) direct-sum reference implementation.

ShellToShellTransferMethod{B<:AbstractShellBinning} <: AbstractEnergyTransferMethod

Compute the directed shell-to-shell transfer matrix T(n,m), where T(n,m) is the rate of energy transfer from shell Sm into shell Sn mediated by the nonlinear advection term.

Fields

• binning::B: Shell binning strategy.

Notes

The antisymmetry property T(n,m) = −T(m,n) holds when the mediator velocity is the full field u (Verma et al. 2002). This is automatically verified by default.

ModeToModeTransferMethod{B, I<:AbstractInvariant} <: AbstractEnergyTransferMethod

Compute the exact mode-to-mode triad transfer S(k|p|q) — energy (or other invariant) given to receiver mode k from giver p, mediated by q, with triad closure k = p + q:

S(k|p|q) = −Im{ [k · û(q)] [û*(k) · û(p)] }.

This is the most fundamental (delta-in-k) mode-to-mode object; it reduces to the shell-to-shell matrix and the spectral transfer T(k) under summation.

Fields

• binning::B: Optional shell binning for reductions to the magnitude-to-magnitude transfer T(K,Q). Use nothing to return the raw per-receiver transfer only.
• invariant::I: Which quadratic invariant to accumulate (default KineticEnergy()).

Cost

O(N^D) per receiver mode; O(N^{2D}) for the full tensor — exact but slow. Guard with a mode-count limit unless force=true.

References

• Dar, Verma & Eswaran (2001); Verma (2004 review, 2019 book).

CoarseGrainingFluxMethod{F<:AbstractFilter, S} <: AbstractEnergyTransferMethod

Compute the pointwise cross-scale energy flux Π_ℓ(x) = −τ̄ᵢⱼ S̄ᵢⱼ at filter scale ℓ.

Fields

• filter::F: Filter kernel (Gaussian, sharp-spectral, or top-hat).
• scale::S: Filter length scale ℓ (same units as the coordinate arrays).

Notes

Physical-space output; suitable for detecting spatial intermittency in the cascade.

TriadicOrthogonalDecompositionMethod{N, O, M} <: AbstractEnergyTransferMethod

Triadic Orthogonal Decomposition (Yeung, Chu & Schmidt 2026).

Operates on temporal snapshots, decomposing triadic (three-wave) nonlinear interactions in the temporal-frequency domain via the mode bispectrum.

Fields

• nfft: DFT block length. nothing for auto-selection.
• noverlap: Block overlap in snapshots. nothing for 50% of window.
• nmode: Number of modes per triad to retain. nothing for nblocks.

References

• Yeung, Chu & Schmidt (2026), J. Fluid Mech. 1031, A34. DOI 10.1017/jfm.2026.11183

SphericalTransferMethod{T<:Real} <: AbstractEnergyTransferMethod

Spectral energy/enstrophy transfer for 2D non-divergent (barotropic) flow on the sphere, in the spherical-harmonic degree spectrum l. Given the vorticity field ζ, with streamfunction ψ = ∇⁻²ζ (so ζ̂_lm = -l(l+1)/a² ψ̂_lm) and advection A = J(ψ,ζ) = u·∇ζ, u = k̂×∇ψ, the transfers are

T_E(l) = -Σ_m Re{ψ̂*_lm Â_lm},   T_Z(l) = Σ_m Re{ζ̂*_lm Â_lm},

both conserving (Σ_l T = 0). See THEORY.md §"Spherical spectral transfer".

Dispatched through calculate_energy_transfer: a regular colatitude–longitude grid (an AbstractMatrix vorticity field) routes to the FastSphericalHarmonics extension; scattered points (a vorticity vector + (θ, φ) coordinates) route to the NUFSHT extension.

Fields

• radius::T: sphere radius a (default 1.0).

DivergentSphericalTransferMethod{T<:Real} <: AbstractEnergyTransferMethod

Spectral kinetic-energy transfer for the full horizontal flow on the sphere — rotational and divergent — in the spherical-harmonic degree spectrum l. Generalises SphericalTransferMethod (which assumes non-divergent/barotropic flow, ∇·u = 0) to a velocity field carrying divergence, and reduces to it exactly in the non-divergent limit.

The input is the horizontal velocity u = (u_θ, u_φ) (colatitude, longitude components), Helmholtz- decomposed as u = k̂×∇ψ + ∇χ (rotational streamfunction ψ, divergent velocity potential χ). Writing the advection in Lamb (rotational) form (u·∇)u = ∇(½|u|²) + ζ (k̂×u) with vorticity ζ = k̂·(∇×u), the nonlinear KE transfer into degree l is the vector-harmonic projection

T(l) = Σ_m Re{ û*_lm · Â_lm},   Â = [(u·∇)u]^  (spin-1 vector-harmonic coefficients),

split by the toroidal (rotational) / spheroidal (divergent) parts of û into T = T_rot + T_div. Total KE is advectively conserved: Σ_l T(l) ≈ 0 (the rotational and divergent channels are not individually conserved — they exchange energy). The Lamb form needs only spin-0/spin-1 transforms (no spin-2). See THEORY.md §"Divergent spherical spectral transfer" (Augier–Lindborg 2013; Burgess–Erler–Shepherd 2013).

Dispatched through calculate_energy_transfer: a regular colatitude–longitude grid (two AbstractMatrix velocity components) routes to the FastSphericalHarmonics extension; scattered points (velocity-component vectors + (θ, φ) coordinates) route to the NUFSHT extension.

Fields

• radius::T: sphere radius a (default 1.0).

AbstractInvariant

Trait supertype selecting which quadratic inviscid invariant a transfer diagnostic accumulates. The same nonlinear-term machinery serves every invariant; only the per-mode transfer-density weighting changes (see Invariants.transfer_density!).

Concrete subtypes: KineticEnergy (default), Helicity (3D), Enstrophy (2D), PassiveScalar (any D).

Advected vs. carrier field

Every transfer diagnostic forms T(k) = Re{ ĉ*(k) · 𝒩̂(k) }, where 𝒩̂ = FFT[(u·∇)f] is the nonlinear term of the advected field f and ĉ is the carrier. For the momentum invariants (KE/helicity/enstrophy) both are the velocity (f = c = u, with vorticity weighting folded into the carrier for helicity/enstrophy). For PassiveScalar the advected and carrier field is the scalar θ, advected by the velocity u — handled by passing the scalar as the primary field and the velocity as advecting_hat.

KineticEnergy <: AbstractInvariant

Kinetic energy E = ½∫|u|². The default invariant; transfer density is Re{ û*(k) · N̂(k) }. Forward cascade in 3D, inverse in 2D.

Helicity <: AbstractInvariant

Helicity H = ∫ u·ω, ω = ∇×u (3D only). Transfer density is Re{ ω̂*(k) · N̂(k) } with ω̂ = i k × û. Co-directional (forward) with energy.

Enstrophy <: AbstractInvariant

Enstrophy Ω = ½∫|ω|², transfer density Re{ ω̂*(k) · N̂_ω(k) } with ω̂ = i k×û and N̂_ω = i k×N̂.

• 2D (scalar vorticity ω̂ = i(k_x û_y − k_y û_x)): enstrophy is an inviscid invariant — conserved (Σ_k T_Ω = 0), counter-directional forward cascade dual to the inverse energy cascade (Kraichnan–Batchelor).
• 3D (vector vorticity): N̂_ω = curl[(u·∇)u] = (u·∇)ω − (ω·∇)u includes vortex stretching, so enstrophy is not conserved (Σ_k T_Ω ≠ 0: net production). This is a valid transfer/budget diagnostic, not a conservative cascade.

Available in 2D and 3D across every diagnostic — spectral flux, shell-to-shell, and the resolved mode-to-mode triad form (the invariant weighting rides the generic transfer_density!).

PassiveScalar <: AbstractInvariant

Passive-scalar variance E_θ = ½∫θ² (Obukhov–Corrsin), advected by the velocity: ∂_tθ + (u·∇)θ = κ∇²θ. The transfer density is Re{ θ̂*(k) N̂_θ(k) } with N̂_θ = FFT[(u·∇)θ].

The scalar is the advected and carrier field; the velocity only advects it. Pass the scalar (shape (ns..., 1)) as the primary field and the velocity as advecting_hat (the convenience entry points calculate_scalar_* do this for you).

Scalar variance is an inviscid invariant for incompressible flow (∫θ(u·∇)θ = −½∫θ²∇·u = 0), so it is conserved (Σ_k T_θ ≈ 0) and cascades forward (to small scales) in both 2D and 3D — unlike kinetic energy there is no inverse-cascade dimension.

A family of canonical invariants

Other quadratic invariants are advected by the velocity exactly like a passive scalar, so their cross-scale transfer is computed by this same path (pass the field as the "scalar"):

• Buoyancy / available-potential-energy variance ½⟨b²⟩ (APE = ½⟨b²⟩/N²) in the Boussinesq system; the −N²w term is a KE↔APE conversion (a separate source, not a triad transfer), so the variance cascade is exactly the scalar transfer of b.
• QG potential enstrophy ½⟨q²⟩ with PV q = ∇²ψ + βy advected by the geostrophic velocity.

References

• Obukhov (1949); Corrsin (1951); Batchelor (1959); QG: Charney (1971); stratified APE: Lindborg (2006). See THEORY.md §0.5.

AbstractFieldDecomposition

Abstract supertype specifying the field decomposition/projection strategy (e.g., Helmholtz rotational/divergent decomposition).

NoDecomposition <: AbstractFieldDecomposition

No decomposition or projection is applied; use the full velocity field.

HelmholtzDecomposition <: AbstractFieldDecomposition

Decompose the velocity field into rotational (solenoidal) and divergent (dilatational) components, computing transfer results for both.

RotationalDecomposition <: AbstractFieldDecomposition

Only compute or retain the rotational (solenoidal/divergence-free) component.

DivergentDecomposition <: AbstractFieldDecomposition

Only compute or retain the divergent (dilatational/curl-free) component.

HelicalDecomposition <: AbstractFieldDecomposition

Decompose a 3D velocity field into its positive- and negative-helicity components via the Craya–Herring/helical basis (Waleffe 1992; Alexakis 2017). For each k ≠ 0 the plane ⊥ k is spanned by the orthonormal helical eigenvectors of the curl, h_±(k) = (e₁ ± i e₂)/√2 with i k̂ × h_± = ± h_± and the Alexakis √2 unit-norm convention (h_± · h_±* = 1, h_+ · h_-* = 0). The velocity projects as û = u_+ h_+ + u_- h_- (u_± = û · h_±*), so

E(k) = E⁺(k) + E⁻(k),   E^±(k) = ½|u_±|²,   H(k) = |k|(|u_+|² − |u_-|²) = 2|k|(E⁺ − E⁻),

recovering the realizability bound |H(k)| ≤ 2|k| E(k). Returns the two vector components (positive = u_+ h_+, negative = u_- h_-); for an incompressible field they sum back to û. 3D only. Used as the decomposition argument to calculate_spectral_flux to get helicity-resolved energy fluxes Π^±(K).

ToroidalPoloidalDecomposition <: AbstractFieldDecomposition

Split a 3D solenoidal velocity into toroidal (horizontal/vortical) and poloidal (vertical/wave) components in the Craya–Herring frame (Craya 1958; Herring 1974; Bartello 1995). For each k with horizontal part k_⊥ ≠ 0, the plane ⊥ k is spanned by e⁽¹⁾ = (k × ẑ)/|k × ẑ| (horizontal, the toroidal direction) and e⁽²⁾ = (k × e⁽¹⁾)/|k| (the poloidal direction); û = u₁ e⁽¹⁾ + u₂ e⁽²⁾. The toroidal part carries the vertical vorticity and has zero vertical velocity; the poloidal part carries the vertical velocity (the linear gravity-wave mode in stratified flow). For purely vertical k (k_⊥ = 0) the split is degenerate and (x̂, ŷ) are used as an arbitrary horizontal orthonormal pair.

Returns (toroidal = u₁ e⁽¹⁾, poloidal = u₂ e⁽²⁾); both are divergence-free and they sum back to the solenoidal part of û. 3D only.

AbstractDealiasing

Strategy for removing aliasing from the pseudospectral quadratic product, passed as the dealiasing keyword. Subtypes: NoDealiasing, OrszagTwoThirds (the default), and PaddedThreeHalves (exact 3/2 zero-padding).

NoDealiasing <: AbstractDealiasing

No dealiasing — the raw pseudospectral product, aliasing included.

OrszagTwoThirds <: AbstractDealiasing

Orszag 2/3-rule truncation: zero modes with |k_d| ≥ N_d/3 in the inputs and output. Exact on the retained band |k| < N/3; the default dealiasing.

PaddedThreeHalves <: AbstractDealiasing

Exact 3/2 zero-padding: form the quadratic product on a (3N/2)-point grid so no aliasing reaches the resolved band, then truncate back to N. Exact for the quadratic nonlinear term over every mode up to Nyquist (nothing discarded), at ~(3/2)^D higher transform cost. Requires FFTW for the fast path.

SpectralFluxResult{KS, V}

Result of a spectral energy flux computation.

Fields

• k_shells::KS: Representative wavenumber for each shell (midpoint of bin edges).
• transfer_spectrum::V: T(k) — energy transfer rate per shell.
• flux::V: Π(K) = +cumsum(T(k)) — cumulative energy flux (Π>0 forward/down-scale cascade, Π<0 inverse; Alexakis–Biferale 2018, THEORY.md §0.5).

Parametric with no element-type bound (works with Float32/Float64/Dual/Unitful, etc.). k_shells (host-side shell wavenumbers) is parametrised separately from the transfer_spectrum/flux data so a device computation can return device data while k_shells stays a host vector.

CompressibleFluxResult{KS, TS, FL, CH, PD}

Result of a compressible kinetic-energy spectral-transfer computation (Singh–Tiwari–Sharma–Verma 2025; see THEORY.md §0.5). The transfer is momentum-weighted (v = ρu), so unlike the incompressible diagnostics it needs the density field and — for the KE↔internal-energy exchange — the pressure.

Fields

• k_shells::KS: representative wavenumber per shell.
• transfer_spectrum::TS: T_u(k) — net momentum-weighted KE transfer into shell k (energy gain rate; sign is opposite the incompressible loss convention). Conserves total KE: Σ_k T_u(k) ≈ 0.
• flux::FL: Π(K) — cumulative flux, Π(K) = Σ_{k>K} T_u(k).
• channels::CH: rotational/compressive (Helmholtz u = u_R + u_C) flux channels as a NamedTuple (rotational, compressive, rot_to_comp, comp_to_rot), or nothing if not requested.
• pressure_dilatation::PD: KE↔IE conversion (rotational = Q_{I,R}(k), compressive = Q_{I,C}(k)) as a NamedTuple, or nothing if no pressure field was supplied.

Each array field carries its own type parameter (element-type/container generic — no shared or <:AbstractVector-bounded param); optional CH/PD resolve to Nothing or a concrete NamedTuple.

ShellToShellResult{V, M, E}

Result of a shell-to-shell energy transfer computation.

Fields

• shell_centers::V: Representative wavenumber for each shell.
• shell_edges::V: Shell boundary wavenumbers (length = N_shells + 1).
• transfer_matrix::M: T(n,m) — Nshells × Nshells matrix; T[n,m] is energy from shell m to shell n.
• net_transfer::V: Σ_m T(n,m) for each receiver shell n (net energy gain of shell n).
• max_antisymmetry_error::E: max |T(n,m) + T(m,n)| — antisymmetry validation metric.

Parametric on vector type V, matrix type M, and scalar type E = eltype(M).

ModeToModeTriadResult{I, KS, A, S}

Result of a fully mode-resolved triad transfer computation.

Fields

• invariant::I: The invariant that was accumulated (e.g. KineticEnergy()).
• ks::KS: The wavenumber vectors (kx, ky[, kz]) defining the spectral grid.
• net_transfer::A: T(k) = Σ_p S(k|p) — net per-mode transfer (shape ns); equals the spectral transfer from calculate_spectral_flux.
• transfer::S: the resolved S(k|p) — energy delivered to receiver mode k from giver mode p (mediated by q=k−p), shape (ns..., ns...) (receiver indices then giver indices). Antisymmetric (S(k|p)=−S(p|k)); summed over p gives net_transfer; summed over shells gives the shell-to-shell matrix.

Parametric on all array/field types — GPU-array friendly.

CoarseGrainingFluxResult{S, A}

Result of a coarse-graining energy flux computation (flux field only).

Fields

• filter_scale::S: Filter scale ℓ used.
• flux_field::A: Π_ℓ(x) pointwise energy flux field (same shape as input velocity).
• mean_flux::S: Area-weighted spatial mean ⟨Π_ℓ⟩.

See also CoarseGrainingFluxResultWithDiagnostics for stress/strain output.

CoarseGrainingFluxResultWithDiagnostics{S, A}

Result of a coarse-graining energy flux computation including stress/strain diagnostics.

Fields

• filter_scale::S: Filter scale ℓ used.
• flux_field::A: Π_ℓ(x) pointwise energy flux field.
• mean_flux::S: Area-weighted spatial mean ⟨Π_ℓ⟩.
• stress_tensor::T: τ̄ᵢʲ (component-indexed array, e.g. (Nx,Ny,2,2)).
• strain_rate::T: S̄ᵢʲ (same array type as stress_tensor).

Returned instead of CoarseGrainingFluxResult when return_diagnostics=true. The tensor diagnostics carry their own type parameter T (higher-rank than the scalar flux_field::A); all fields are always present — no Union{Nothing,...} type instability.

TriadicOrthogonalDecompositionResult{V, A3, PM, EC, XM}

Result container for Triadic Orthogonal Decomposition.

Fields

• frequencies::V: Frequency vector (length nFreq).
• mode_bispectrum::A3: Singular values λ(fl, fn, mode) — array of size (nFreq, nFreq, nmode).
• modes::PM: Dict mapping (l, n) index tuples to mode arrays. Each value contains convective modes (index 1 along first dim) and recipient modes (index 2) with spatial/variable dimensions.
• modal_energy_budget::A3: Energy transfer T(fl, fn, mode) per triad per mode. Same shape as mode_bispectrum.
• expansion_coefficients::EC: Expansion coefficients, or nothing if not requested.
• auxiliary_modes::XM: Dict mapping (l, n) to donor/catalyst modes, or nothing.

All fields are typed — the optional EC/XM parameters resolve to Nothing or the concrete container type at construction, so the struct is type-stable (no untyped Any fields).

SphericalTransferResult{V<:AbstractVector}

Result of a spherical spectral energy/enstrophy transfer (SphericalTransferMethod), indexed by spherical-harmonic degree l = 0…lmax.

Fields

• degrees::V: the degrees l.
• energy_transfer::V: T_E(l) — nonlinear kinetic-energy transfer into degree l; Σ_l T_E ≈ 0.
• enstrophy_transfer::V: T_Z(l) — enstrophy transfer into degree l; Σ_l T_Z ≈ 0.
• energy_flux::V: Π_E(L) = -Σ_{l≤L} T_E(l) — cumulative up-degree energy flux.
• enstrophy_flux::V: Π_Z(L) = -Σ_{l≤L} T_Z(l).

DivergentSphericalTransferResult{V<:AbstractVector}

Result of the divergent horizontal kinetic-energy spectral transfer (DivergentSphericalTransferMethod), indexed by spherical-harmonic degree l = 0…lmax.

Fields

• degrees::V: the degrees l.
• energy_transfer::V: total horizontal-KE transfer T(l) = T_rot(l) + T_div(l) into degree l; the skew-symmetric (energy-conserving) advection makes Σ_l T ≈ 0.
• energy_flux::V: Π(L) = -Σ_{l≤L} T(l) — cumulative up-degree KE flux.
• rotational_transfer::V: rotational-channel transfer T_rot(l) — projection of the advection onto the toroidal (streamfunction ψ) part of the velocity.
• divergent_transfer::V: divergent-channel transfer T_div(l) — projection onto the spheroidal (velocity-potential χ) part.
• rotational_flux::V, divergent_flux::V: cumulative fluxes -Σ_{l≤L} T_rot, -Σ_{l≤L} T_div.

Only the total is conserved (Σ_l T ≈ 0); the two channels exchange energy, so Σ_l T_rot and Σ_l T_div are individually nonzero (equal and opposite up to the total).

NonlinearTermWorkspace{CA, RA, GA, P}

Preallocated buffers for the generalized pseudospectral nonlinear term 𝒩(k) = FFT[(u·∇)f], where the advecting velocity u has nd advecting (spatial) components and the advected field f has M components. For the momentum term f = u (M = D); for passive-scalar / vector-potential advection f = θ/a (M = 1).

Fields

• u_phys::RA: (ns..., nd) real physical-space advecting velocity (rank nd+1); only the nd spatial directions of the velocity participate in (u·∇), so this never depends on D.
• grad_phys::GA: (ns..., M, nd) real physical-space gradients ∂fi/∂xj (rank nd+2).
• N_phys::RA: (ns..., M) real physical-space nonlinear term (rank nd+1).
• N̂::CA: (ns..., M) complex spectral output buffer (rank nd+1).
• plans::P: FFT plan/scratch bundle (set by the FFTW extension) or nothing.

Parametric on the concrete array types CA (complex), RA (real, rank nd+1), GA (real gradient buffer, rank nd+2), and the plan-bundle type P — no element-type bounds, and each field is concretely typed (grad_phys has a separate parameter because its rank differs from the others; u_phys and N_phys share RA — same rank/eltype, possibly different trailing extent).

SpectralFluxWorkspace{NW, V, A}

Preallocated buffers for calculate_spectral_flux!.

Fields

• nonlinear::NW: NonlinearTermWorkspace for computing N̂(k).
• T_spec::V: Shell transfer spectrum buffer (length N_sh).
• flux::V: Cumulative flux buffer (length N_sh).
• transfer_density::A: Per-mode transfer density buffer.

ShellToShellWorkspace{NW, CA, M, V, IA}

Preallocated buffers for calculate_shell_to_shell_transfer!.

Fields

• nonlinear::NW: NonlinearTermWorkspace (owns N̂_m and all physical-space temps).
• û_m::CA: Band-filtered mediator velocity buffer (reused each shell m).
• T_mat::M: Output transfer matrix (Nsh × Nsh), written in-place.
• net_transfer::V: Net per-shell transfer buffer (length N_sh).
• shell_idx::IA: Integer shell-index array (same shape as k_mag).

CompressibleWorkspace(velocity_hat, ks; spectral=SpectralBackends.FFTSpectralBackend())

Reusable field-scratch + transform context for calculate_compressible_flux! — holds the FFT plans and every (ns...)-sized intermediate of the momentum-weighted budget (velocity/density/ momentum physical & spectral fields, gradients, nonlinear terms, and the R/C-channel + pressure- dilatation scratch), so repeated per-snapshot calls allocate ~0 field memory (only the small shell vectors of each result). Build once for a given grid/precision and reuse across snapshots.

SphericalTransferWorkspace(lmax; ...)

Reusable buffers for calculate_spherical_transfer!. Constructed by the FastSphericalHarmonics (regular grid) or NUFSHT (scattered) extension; loading one provides the constructor. The buffer fields are typed via parameters so the core names no extension type.

ScatteredSphericalTransferWorkspace(coords, lmax; ...)

Reusable buffers for the SCATTERED-point spherical transfer !() (NUFSHT extension). Holds the three NUFSHT spin plans (spin-0 at lmax, spin-1 at lmax, spin-0 at the dealiased lwork = 2·lmax) with the scattered points preset, plus every coefficient/gradient/reduction buffer and the result vectors. The plans are the dominant cost (FINUFFT planning + the CG least-squares setup); building them once lets a snapshot time series on the same points reuse them. NUFSHT plans self-finalize their FINUFFT resources, so this struct needs no finalizer. Fields are typed via parameters so the core names no NUFSHT type. Requires using NUFSHT.

DivergentSphericalTransferWorkspace(lmax; ...)

Reusable buffers for the regular-grid divergent KE transfer !() (FastSphericalHarmonics extension; loading it provides the constructor). FastSphericalHarmonics has no in-place transform API, so the spin-weighted transforms/eth allocate internally on every call (an irreducible floor); the workspace therefore just carries the reused DivergentSphericalTransferResult and the resolution parameters. Fields are typed via parameters so the core names no extension type. Requires using FastSphericalHarmonics.

ScatteredDivergentSphericalTransferWorkspace(coords, lmax; ...)

Reusable buffers for the SCATTERED-point divergent KE transfer !() (NUFSHT extension; loading it provides the constructor). Holds the five NUFSHT spin plans (spin ±1 and spin 0 at lmax, spin 0 and spin +1 at the dealiased lwork = 2·lmax) with the points preset, plus every coefficient/field/reduction buffer and the result. The plans are the dominant reusable cost. Fields are typed via parameters so the core names no NUFSHT type. Requires using NUFSHT.

wavenumber_grid(ns, Ls) -> NTuple{D, Vector{Float64}}

Return a tuple of 1D physical-wavenumber vectors matching the FFTW fftfreq convention (centered at zero after fftshift).

Arguments

• ns::NTuple{D,Int}: Number of grid points along each dimension.
• Ls::NTuple{D,Float64}: Physical domain size along each dimension.

Returns

Tuple of length D; element d is the range kd ∈ [−⌊Nd/2⌋, ⌊(Nd−1)/2⌋] × (2π / Ld).

Example

ks = wavenumber_grid((16, 16), (2π, 2π))
# ks[1] and ks[2] are both [-8, -7, ..., 7] * (2π/2π)

wavenumber_magnitude_grid(ks) -> Array{Float64, D}

Compute the isotropic wavenumber magnitude |k| at every grid point.

Arguments

• ks::NTuple{D, AbstractVector}: Tuple of 1D wavenumber vectors (e.g., from wavenumber_grid).

Returns

D-dimensional array of the same size as the full grid, with entry [i₁,…,i_D] = sqrt(ks[1][i₁]² + … + ks[D][i_D]²).

dealiasing_mask(ns; rule=:twothirds) -> BitArray{D}

Build a spectral dealiasing mask.

Arguments

• ns::NTuple{D,Int}: Grid sizes.
• rule::Symbol: :twothirds (2/3 rule) or :half.

Returns

BitArray of the same shape as the spectral grid; true where the mode is kept (i.e., |kd| < Nd/2 * threshold for all d).

Notes

For the 2/3 rule, modes with |kd| ≥ Nd/3 along any dimension are zeroed.

dealiasing_mask!(mask, ns; rule=:twothirds) -> mask

In-place version of dealiasing_mask: write the keep/discard Bool mask for grid shape ns into mask (rule = :twothirds for the Orszag 2/3 cutoff, otherwise the Nyquist half).

assign_shells(k_mag, edges) -> Array{Int}

Return an integer array (same shape as k_mag) where [I] = n if edges[n] <= k_mag[I] < edges[n+1], and 0 if the mode falls outside all shells.

One integer per mode (single allocation, cache-friendly): the canonical shell-membership representation used by every transfer accumulation kernel.

shell_edges(binning, k_max) -> Vector{Float64}

Return the monotonically increasing shell boundary vector for binning up to k_max.

The resulting vector has length n_shells(binning, k_max) + 1; shell n covers wavenumbers in [edges[n], edges[n+1]).

shell_centers(binning, k_max) -> Vector{Float64}

Return the geometric midpoint of each shell. For logarithmic binnings, this is the geometric mean of the edge pair; for linear binnings, the arithmetic mean.

n_shells(binning, k_max) -> Int

Return the number of shells for binning up to k_max.

shell_coordinate(geometry, ks) -> Array

Return an array (shape ns) giving the wavenumber coordinate each mode is binned by under geometry. For ShellMagnitude this is √(Σ_{d∈dims} k_d²) — |k| when dims covers all dimensions (isotropic), k_⊥/k_∥ for an anisotropic projection.

AbstractShellBinning

Supertype for shell spacing strategies (how wavenumber space is partitioned into shells): LinearBinning, LogarithmicBinning, DyadicBinning, CustomBinning. Orthogonal to the shell coordinate (AbstractShellGeometry).

LinearBinning(Δk) <: AbstractShellBinning

Uniform shell spacing: k_n = n · Δk.

Fields

• Δk: Shell width in physical wavenumber units.

LogarithmicBinning(k₀, λ) <: AbstractShellBinning

Geometrically-spaced shells: k_n = k₀ · λⁿ.

Fields

• k₀: First shell lower edge (> 0).
• λ: Ratio between consecutive shell edges (> 1); λ = 2 gives dyadic.

DyadicBinning(k₀) <: AbstractShellBinning

Dyadic (octave) shells: k_n = k₀ · 2ⁿ. Equivalent to LogarithmicBinning(k₀, 2.0).

Fields

• k₀: First shell lower edge (> 0).

CustomBinning(edges) <: AbstractShellBinning

User-specified shell edges. Shell n covers wavenumbers in [edges[n], edges[n+1]).

Fields

• edges: Monotonically increasing edge values (length = N_shells + 1).

SmoothBands(centers; logwidth=0.6)

Graded (smooth) spectral bands for band-to-band transfer T(K,Q) (Eyink & Aluie 2009), as an alternative to the sharp shells of AbstractShellBinning. Each band n weights a mode at coordinate κ by a log-Gaussian exp(−(ln(κ/centers[n]))² / (2·logwidth²)), renormalized across bands to a partition of unity (Σ_n w_n(κ) = 1) so the smooth bands conserve and reduce to the band-summed transfer spectrum. Smaller logwidth → sharper, more shell-like bands.

Fields

• centers: band-center wavenumbers (monotonically increasing, all > 0).
• logwidth: Gaussian width in ln κ (dimensionless); default 0.6 (≈ one octave overlap).

AbstractShellGeometry

Abstract supertype selecting the wavenumber coordinate the shells partition (isotropic |k|, or an anisotropic projection like k_⊥/k_∥). Orthogonal to the binning spacing (AbstractShellBinning).

ShellMagnitude(dims) <: AbstractShellGeometry

Bin modes by the Euclidean magnitude of the wavenumber components in dims: κ(k) = √(Σ_{d∈dims} k_d²). dims === nothing uses all spatial dimensions (isotropic |k|).

Use the constructors IsotropicShells, PerpendicularShells, ParallelShells for the common cases.

IsotropicShells() -> ShellMagnitude

Isotropic shells over |k| (all dimensions) — the default geometry.

PerpendicularShells(dims=(1, 2)) -> ShellMagnitude

Cylindrical shells over k_⊥ = √(Σ_{d∈dims} k_d²) (the horizontal plane by default), giving the anisotropic perpendicular flux Π(k_⊥) for rotating/stratified flows.

ParallelShells(dims=(3,)) -> ShellMagnitude

Plane shells over k_∥ = √(Σ_{d∈dims} k_d²) (the vertical axis by default), giving the anisotropic parallel flux Π(k_∥).

AbstractFilter

Supertype for spectral filter kernels used by coarse-graining: SharpSpectralFilter, GaussianFilter, TopHatFilter.

SharpSpectralFilter <: AbstractFilter

Ideal low-pass (brick-wall) filter in spectral space: Ĝ(k, ℓ) = 1 if |k| < π/ℓ, else 0.

Provides exact scale separation but produces Gibbs ringing in physical space.

GaussianFilter <: AbstractFilter

Gaussian filter in spectral space: Ĝ(k, ℓ) = exp(−k² ℓ² / 24).

Excellent physical-space locality; widely used in LES and coarse-graining studies. The normalisation factor 24 follows the convention of Aluie et al. (2018).

TopHatFilter <: AbstractFilter

Top-hat (box) filter in physical space; sinc response in spectral space: Ĝ(k, ℓ) = sinc(k ℓ / (2π)).

Compact support in physical space; standard in LES.

filter_response(filter, k, ℓ) -> Real

Evaluate the spectral transfer function Ĝ(k, ℓ) of filter at wavenumber magnitude k and filter scale ℓ.

The filter scale ℓ is defined so that the filter retains scales larger than ℓ.

apply_filter_spectral(û_in, k_mag, filter, ℓ) -> Array

Non-mutating version of apply_filter_spectral!.

apply_filter_spectral!(û_out, û_in, k_mag, filter, ℓ)

Apply filter at scale ℓ in spectral space by pointwise multiplication: ûout[I] = Ĝ(|k[I]|, ℓ) * ûin[I].

Arguments

• û_out: Output spectral array (same shape as û_in).
• û_in: Input spectral array (complex or real).
• k_mag: Array of wavenumber magnitudes (same shape as û_in).
• filter::AbstractFilter: Filter kernel.
• ℓ::Real: Filter scale.

Modifies û_out in-place and returns it.