API Reference

CoarseGrainResult(scales, Π, cumulative_energy, wavenumber, filtering_spectrum)

Container for results of a complete coarse-graining multiscale analysis. Every field's container type is a type parameter, inferred by the constructor above, so nothing here is stored behind an abstract annotation.

Fields

• scales::AbstractVector{T}: filter scales ℓ in meters
• Π::A: energy-flux maps stacked into ONE contiguous (spatial dims..., Nscales) array (W/m³) — not a Vector of separately-allocated per-scale matrices, so the whole sweep is a single allocation and each scale's map is a zero-copy view (compute_Π! writes directly into its slice).
• cumulative_energy::AbstractArray{T}: cumulative coarse specific KE ½⟨|ū_ℓ|²⟩ per scale (Sadek–Aluie Eq.
  1. — a Vector (per scale) for coarse_grain, or a (Nlevels, Nscales) Matrix for
  coarse_grain_profile (per vertical level AND scale — deliberately not summed across levels, since that would need volume/thickness weighting this function doesn't have).
• wavenumber::AbstractVector{T}: filtering wavenumber k_ℓ = L/ℓ per scale (level-independent)
• filtering_spectrum::AbstractArray{T}: filtering spectral density Ẽ(k_ℓ) per scale (Eq. 14), same shape convention as cumulative_energy

Examples

res = coarse_grain(u, v, grid; scales=[10e3, 20e3, 30e3], kernel=TopHatKernel())
# Access results:
res.scales[1]              # First scale (10 km)
res.Π[:, :, 1]              # Energy-flux map at 10 km (a view; use `@view` to avoid copying)
res.cumulative_energy[1]   # cumulative coarse KE at 10 km
res.filtering_spectrum[1]  # filtering spectral density at k_ℓ = res.wavenumber[1]

coarse_grain!(result, u, v, w, grid; scales, kernel, workspace, filter_plans, deriv_plan, backend, mask_strategy, method, L)
coarse_grain!(result, u, v, grid; scales, ...)  # 2D convenience wrapper

In-place coarse_grain: refills an existing CoarseGrainResult's buffers — scales, the stacked Π array, cumulative energy, wavenumber, spectrum — instead of allocating fresh ones. Supplying workspace and filter_plans reuses the scratch arrays and the per-scale plans too, which is the zero-reallocation entry point for a sweep repeated across timesteps.

result must already be sized for length(scales) scales over grid's shape; a mismatch throws DimensionMismatch.

One driver for every grid architecture: the scale loop, the plan reuse and the spectrum are the same for all of them. Only the derivative plan and the trailing dimension of result.Π vary.

coarse_grain(u, v, w, grid; scales, kernel=TopHatKernel(), backend=AutoBackend(), mask_strategy=Deformable(), method=nothing, L=1)
coarse_grain(u, v, grid; scales, ...)  # 2D convenience wrapper

Perform complete coarse-graining analysis across multiple filter scales, allocating a fresh CoarseGrainResult and workspace. This is a thin wrapper around coarse_grain!; for repeated sweeps over the same grid/scales (e.g. successive timesteps), allocate the result once and call coarse_grain! directly to reuse its buffers.

Arguments

• u::AbstractMatrix: Eastward/zonal velocity component (m/s)
• v::AbstractMatrix: Northward/meridional velocity component (m/s)
• w::Union{Nothing,AbstractMatrix}: Vertical velocity (nothing for 2D)
• grid::StructuredGrid: Grid geometry and active-cell mask

Keyword Arguments

• scales::AbstractVector: Vector of filter scales ℓ in meters (e.g., 10e3:10e3:100e3)
• kernel::AbstractFilterKernel=TopHatKernel(): Filter kernel
• backend::AbstractExecutionBackend=AutoBackend(): Execution backend
• mask_strategy::AbstractMaskStrategy=Deformable(): Masking strategy (ZeroFill() or Deformable())
• method::Union{Nothing,AbstractFilterMethod}=nothing: filtering engine; nothing takes plan_filter's per-grid default (real space where a grid has that engine)
• L::Real=1: reference length setting the wavenumber normalization k_ℓ = L/ℓ

Returns

• CoarseGrainResult: Container with scales, Π maps, and spectrum

Examples

geom = SphericalGeometry(6371000.0)
grid = StructuredGrid(geom, lon_rad, lat_rad, mask)
scales = collect(10e3:10e3:100e3)
res = coarse_grain(u, v, grid; scales=scales, kernel=TopHatKernel())
plot(res.wavenumber, res.filtering_spectrum, xscale=:log10, yscale=:log10)
heatmap(res.Π[:, :, 3])  # 3rd scale is 30 km

coarse_grain_profile(u, v, w, grid; scales, kernel=TopHatKernel(), backend=AutoBackend(), mask_strategy=Deformable(), method=nothing, L=1)

Vertical-profile sweep: given 3D (x, y, z) velocity arrays, runs Diagnostics.compute_Π_profile! (the literature-standard independent-per-level 2D/2.5D method — see Diagnostics.compute_Π!'s thin-layer/QG regime note) at every scale, returning a CoarseGrainResult whose Π is one contiguous (Nx, Ny, Nlevels, Nscales) array. The workspace is built once and reused across the whole level × scale sweep.

PiDecomposedWorkspace(grid)

Scratch for compute_Π_decomposed!: the divergent components, the four filtered means, the three stress tensors, the two strain tensors, the four flux fields and three shared temporaries.

TauWorkspace(grid)

Scratch for tau_decomposition!: the nine output components, the filtered fields and residuals, and two product buffers. Allocated once and reused, so a repeated decomposition — over timesteps, or over scales — costs no allocation after the first.

TracerFluxWorkspace(grid)

Scratch for tracer_variance_flux!: the filtered fields, the two products, the subfilter flux components and the resolved tracer gradient. Allocated once and reused.

ΠWorkspace{T, A}

Pre-allocated arrays for computing cross-scale energy flux Π to avoid heap allocations in scale loops.

active_area(grid) -> T

Total area of the active cells: the denominator of every spatial average here.

compute_Π!(Π::AbstractArray{T,3}, u, v, w, grid::StructuredGrid{Cartesian,T,3}, kernel, scale; mask_strategy=Deformable(), backend=AutoBackend())

Full three-dimensional Cartesian cross-scale energy flux Π = -S̄ij τij with all nine strain components (the diagonal S_zz = ∂w̄/∂z and the off-diagonals S_xz, S_yz carry genuine vertical derivatives, unlike the 2.5D layer-by-layer path). The 3D grid carries a 3D mask, so masked cells are handled per-cell in all three directions.

The contraction is the symmetric six-term sum S̄:τ = S_xx τ_xx + S_yy τ_yy + S_zz τ_zz + 2(S_xy τ_xy + S_xz τ_xz + S_yz τ_yz).

Dispatched on a 3D output array + 3D Cartesian grid (the 2D method takes an AbstractMatrix); see the separate StructuredGrid{Spherical,T,3} method below for the spherical volumetric case (genuine radius axis, real ∂/∂r, full curvature-corrected strain). Pass a reusable workspace (a ΠWorkspace, dimension-generic) to avoid reallocating temporaries on every call — the same "build once, reuse many" pattern the 2D driver uses, now that ΠWorkspace infers its array type from the grid's actual shape instead of hardcoding Matrix.

compute_Π!(Π::AbstractArray{T,3}, u, v, w, grid::StructuredGrid{Spherical,T,3}, kernel, scale; workspace=nothing, backend=AutoBackend(), mask_strategy=Deformable())

Full three-dimensional spherical cross-scale energy flux Π = -S̄ij τij: a genuine radius axis r[k] (absolute distance from the planet center — see FlowGeometries.Grids.StructuredGrid's 3D constructor) and real vertical derivatives ∂/∂r, unlike the 2.5D layer-by-layer path (which drops the u_r/r curvature terms in S_ee/S_nn and the S_er/S_nr/S_rr radial strain entirely, since it has no radial axis to differentiate against).

Velocities are rotated to planetary Cartesian for filtering (Aluie 2019 commutativity), then rotated back to local (east, north, radial), through the same _rotate_stress_to_local_enr/_sfs_contraction kernels the 2D spherical driver uses — that rotation is fully 3×3-general, and the 2.5D caller simply discards its radial components. What differs here is the strain: the spherical strain-rate tensor in orthogonal curvilinear coordinates (scale factors h_λ = r cosφ, h_φ = r, h_r = 1),

S_ee = (1/(r cosφ))∂ū_e/∂λ - v̄_n·tanφ/r + w̄_r/r
S_nn = (1/r)∂v̄_n/∂φ + w̄_r/r
S_rr = ∂w̄_r/∂r
S_en = ½[(1/(r cosφ))∂v̄_n/∂λ + (1/r)∂ū_e/∂φ + ū_e·tanφ/r]
S_er = ½[(1/(r cosφ))∂w̄_r/∂λ + ∂ū_e/∂r - ū_e/r]
S_nr = ½[(1/r)∂w̄_r/∂φ + ∂v̄_n/∂r - v̄_n/r]

where ∂/∂λ/∂/∂φ/∂/∂r are Derivatives.ddx!/Derivatives.ddy!/Derivatives.ddz! (already metric-scaled using the LOCAL r[k], not the fixed reference radius). Pass a reusable workspace to avoid reallocating temporaries on every call, exactly as the Cartesian 3D method does.

compute_Π!(Π, u, v, w, grid, kernel, scale; workspace=nothing, backend=AutoBackend(), mask_strategy=Deformable())

Compute the cross-scale kinetic energy flux Π = -S̄ij τij at filter scale ℓ.

This implements the coarse-graining framework of Aluie et al. (2018) for computing energy transfer across scales in turbulent flows. Positive Π indicates forward cascade (energy from large to small scales), negative Π indicates inverse cascade.

Arguments

• Π::AbstractMatrix{T}: Output array for energy flux (modified in-place)
• u::AbstractMatrix: Eastward/zonal velocity component
• v::AbstractMatrix: Northward/meridional velocity component
• w::Union{Nothing,AbstractMatrix}: Vertical velocity (nothing for 2D calculations)
• grid::StructuredGrid: Grid geometry and coordinates
• kernel::AbstractFilterKernel: Filter kernel
• scale::T: Filter scale ℓ in meters

Keyword Arguments

• workspace=nothing: Pre-allocated ΠWorkspace for intermediate arrays
• backend::AbstractExecutionBackend=AutoBackend(): Execution backend
• mask_strategy::AbstractMaskStrategy=Deformable(): Masking strategy (ZeroFill() or Deformable())

Physics

The cross-scale energy flux is computed as:

Π = -S̄_ij * τ_ij

where:

• S̄_ij = 0.5 * (∂ū_i/∂x_j + ∂ū_j/∂x_i) is the resolved strain rate tensor
• τ_ij = [u_i*u_j]̄ - ū_i*ū_j is the subfilter-scale (SFS) stress tensor
• Overbar denotes filtered quantities

For spherical geometry, velocity components are transformed to planetary Cartesian coordinates before filtering to ensure commutativity with derivatives (Aluie 2019).

Physics regime: 2.5D thin-layer/quasi-geostrophic approximation when w is supplied

When w !== nothing, this method still computes only a SINGLE 2D layer's tensor: it includes the cross terms S_xz = ½∂ū/∂x, S_yz = ½∂v̄/∂y in the strain contraction, but sets S_zz = ∂w̄/∂z ≡ 0 and never differentiates u/v/w in the vertical — there is no 3rd spatial dimension in the input arrays for it to differentiate against. This is not a shortcut; it is the standard thin-layer (small aspect ratio δ = H/L) / quasi-geostrophic scaling used throughout large-scale ocean and atmosphere dynamics (Vallis, Atmospheric and Oceanic Fluid Dynamics, §5; Pedlosky, Geophysical Fluid Dynamics, ch. 6), under which vertical shear terms are genuinely subdominant to horizontal gradients — valid for the normal large-scale, stratified, rotating-flow regime this package targets, NOT for homogeneous/isotropic 3D turbulence (e.g. boundary-layer or Rayleigh–Taylor studies), where filtering genuinely blends all three directions and vertical derivatives are real, not assumed away. The literature on "vertical structure via coarse-graining" (Aluie, Hecht & Vallis 2018, JPO; Buzzicotti, Storer, Khatri, Griffies & Aluie 2023, JAMES) analyzes vertical structure by running this SAME 2D/2.5D method independently at each z level of a multi-level dataset and comparing/stacking the resulting profiles — not by computing a coupled 3D tensor — so Pipeline.coarse_grain_profile/ compute_Π_profile! (looping this method per level) is the literature-matching way to get a vertical-structure result. A genuinely coupled, all-nine-strain-component 3D method exists separately for the true-3D Cartesian case (see the AbstractArray{T,3} compute_Π! method).

Returns

• Π: Energy flux array (same as input), units of W/m³

Examples

Π = zeros(100, 100)
compute_Π!(Π, u, v, nothing, grid, TopHatKernel(), 30000.0)
# Π now contains energy flux at 30 km scale

References

• Aluie et al. (2018): https://doi.org/10.1175/JPO-D-17-0100.1
• Aluie (2019): https://doi.org/10.1007/s13137-019-0123-9
• Vallis, G.K., Atmospheric and Oceanic Fluid Dynamics, 2nd ed., Cambridge University Press, 2017.
• Pedlosky, J., Geophysical Fluid Dynamics, 2nd ed., Springer, 1987.
• Aluie, Hecht & Vallis (2018), J. Phys. Oceanogr. 48(2): https://doi.org/10.1175/JPO-D-17-0100.1
• Buzzicotti, Storer, Khatri, Griffies & Aluie (2023), J. Adv. Model. Earth Syst.: https://doi.org/10.1029/2021MS002583

compute_Π!(Π::AbstractVector, u, grid::StructuredGrid{Cartesian,T,1}, kernel, scale; workspace=nothing, backend=AutoBackend(), mask_strategy=Deformable())

1D cross-scale energy flux Π = -S̄xx τxx on a genuinely 1D StructuredGrid (a single scalar velocity component u along one axis — the 1D analog of the 2D tensor contraction, which reduces to a single term since there's only one strain/stress component). Not the 2D-with-singleton-dimension case (which reuses the 2D methods directly).

compute_Π!(Π, u, v, w, grid::CurvilinearGrid, kernel, scale; workspace=nothing, deriv_plan=nothing, backend=AutoBackend(), mask_strategy=Deformable())

Cross-scale kinetic energy flux Π = -S̄ij τij on a FlowGeometries.Grids.CurvilinearGrid. Identical physics to the StructuredGrid 2D method — it shares the same _compute_Π! tensor kernel — but the resolved strain uses the least-squares tangent-plane gradient (Discretization.gradient! over a Connectivity.gradient_plan, both components from one neighbour sweep) and real-space filtering uses the scattered per-point footprint. Pass a prebuilt deriv_plan = FG.Connectivity.gradient_plan(grid) (and a reusable workspace) to avoid rebuilding them per call across a scale sweep.

compute_Π!(Π, u, v, w, grid::UnstructuredGrid, kernel, scale; workspace=nothing, deriv_plan=nothing, backend=AutoBackend(), mask_strategy=Deformable(), method=Spectral())

Cross-scale kinetic energy flux Π = -S̄ij τij on a FlowGeometries.Grids.UnstructuredGrid (scattered points, node-indexed) — the same physics as the 2D methods (planetary-Cartesian rotation for spherical geometry), via _compute_Π!. The resolved strain uses the node-indexed WLSQ gradient (Connectivity.gradient_plan + Discretization.gradient!). method defaults to Spectral() here, unlike the other grid types' RealSpace() default: the transform is exact for a band-limited field and its per-apply cost does not grow with the filter scale. RealSpace() applies the kernel as written, with compact support; a transform's support is global.

compute_Π_decomposed!(ws, u, v, u_rot, v_rot, grid, kernel, scale; filter_plan=nothing, deriv_plan=nothing, ...)
    -> (; total, rotational, cross, divergent)

In-place compute_Π_decomposed. Returns views of ws's buffers, valid until the next call on the same workspace. With ws and both plans supplied, a repeated evaluation allocates nothing.

compute_Π_decomposed(u, v, w, u_rot, v_rot, w_rot, grid::StructuredGrid{Cartesian,T,3}, kernel, scale; backend=AutoBackend(), mask_strategy=Deformable())
    -> (; total, rotational, cross, divergent)

True three-dimensional analog of the 2D compute_Π_decomposed above: the same both-sides (strain AND stress) rotational/divergent split — see that method's docstring for the derivation — generalized to all six independent strain/stress tensor components, contracted the same way the true-3D compute_Π! method does (nine-term symmetric contraction).

compute_Π_decomposed(u, v, u_rot, v_rot, grid, kernel, scale; backend=AutoBackend(), mask_strategy=Deformable())
    -> (; total, rotational, cross, divergent)

Split the 2D Cartesian cross-scale KE flux Π = -S̄ij τij into rotational-rotational (ΠRR), divergent-divergent (ΠDD), and cross/interaction (Π_X — the "stimulated cascade" channel of Barkan, Srinivasan & McWilliams 2024, JPO) parts, by decomposing both sides of the bilinear contraction, not just the stress.

The Helmholtz decomposition itself is NOT recomputed here — pass the rotational (solenoidal, divergence-free) part (u_rot, v_rot) from a Helmholtz solver (e.g. HelmholtzDecomposition.jl); the divergent (irrotational) part is taken as the complement (u, v) - (u_rot, v_rot). Writing u = uʳ + uᵈ:

• The strain S̄ is LINEAR in velocity, so it splits with no cross term: S̄ = S̄ʳ + S̄ᵈ.
• The stress τ is BILINEAR (quadratic in velocity), so it splits into three pieces: τ(u,u) = τ(uʳ,uʳ) + τ(uᵈ,uᵈ) + [τ(uʳ,uᵈ) + τ(uᵈ,uʳ)] = τʳʳ + τᵈᵈ + τ_X.

Substituting both splits into Π = -S̄:τ = -(S̄ʳ+S̄ᵈ):(τʳʳ+τᵈᵈ+τ_X) and expanding the six resulting terms into three physically named channels:

Π_RR = -S̄ʳ:τʳʳ                                        (pure rotational-to-rotational cascade)
Π_DD = -S̄ᵈ:τᵈᵈ                                        (pure divergent-to-divergent cascade)
Π_X  = -(S̄ʳ:τᵈᵈ + S̄ᵈ:τʳʳ + S̄ʳ:τ_X + S̄ᵈ:τ_X)          (all rotational/divergent interaction terms)

so the channels sum exactly to the total flux, Π = ΠRR + ΠX + Π_DD — each piece constructed directly (not as a residual), yet the identity holds by the same bilinearity/linearity argument. An earlier version of this function computed all three channels by contracting the split stress against the full, undecomposed strain S̄ (a one-sided split); that's only correct in the special case S̄ᵈ ≡ 0, and silently wrong whenever the divergent part itself has nonzero strain.

Returns a named tuple of flux maps (W m⁻³): rotational = ΠRR, divergent = ΠDD, cross = Π_X.

compute_Π_profile!(Π, u, v, w, grid, kernel, scale; workspace=nothing, backend=AutoBackend(), mask_strategy=Deformable())

Vertical-profile energy flux: given 3D (x, y, z) velocity arrays, runs the 2D/2.5D compute_Π! INDEPENDENTLY at each z level — the literature-standard way to obtain vertical structure via coarse-graining (Aluie, Hecht & Vallis 2018; Buzzicotti et al. 2023; see the thin-layer/QG regime note on compute_Π!'s docstring) — writing each level's 2D result into the matching slice of the 3D Π output. This is NOT a coupled 3D tensor computation (no vertical derivatives are taken across levels); for that, see the true-3D compute_Π! method instead.

cumulative_energy!(spectrum, u, v, w, grid, kernel, scales; workspace=nothing, backend=AutoBackend(), mask_strategy=Deformable())

In-place cumulative_energy: writes into the caller-supplied spectrum vector and, when workspace (a ΠWorkspace) is supplied, reuses its u_filt/v_filt/w_filt scratch arrays instead of allocating fresh ones — the same buffers compute_Π! already fills at each scale, so a coarse_grain! sweep pays for this filtered-velocity scratch space once, not twice.

cumulative_energy(u, v, w, grid, kernel, scales; backend=AutoBackend(), mask_strategy=Deformable())

Cumulative coarse-grained kinetic energy E(ℓ) = 0.5 ⟨|ū_ℓ|²⟩ at each filter scale (Sadek & Aluie 2018, PRF, Eq. 15). This is the CUMULATIVE quantity; the filtering spectral DENSITY (comparable to a Fourier energy spectrum) is its derivative w.r.t. filtering wavenumber — see filtering_spectrum. Allocates a fresh spectrum vector each call; for a repeated sweep (e.g. inside coarse_grain!), call cumulative_energy! directly with a reused buffer.

Examples

scales = collect(10000.0:10000.0:100000.0)  # 10-100 km
E = cumulative_energy(u, v, nothing, grid, TopHatKernel(), scales)
# E[i] is the cumulative coarse KE at scale scales[i]

References

• Sadek & Aluie (2018), Phys. Rev. Fluids 3, 124610 — extracting the spectrum by filtering.

energy_from_filtered(ws, grid, has_w, total_area) -> E(ℓ)

E(ℓ) = ½⟨|ū_ℓ|²⟩ read from the filtered velocities already in ws. compute_Π! leaves exactly those there, so calling this straight after it filters u/v once per scale rather than twice.

filtering_spectrum(u, v, w, grid, kernel, scales; L=1, backend=AutoBackend(), mask_strategy=Deformable())
    -> (k_ℓ, Ẽ)

Filtering spectral DENSITY (Sadek & Aluie 2018, PRF, Eq. 14): the derivative of the cumulative coarse-grained KE w.r.t. the filtering wavenumber k_ℓ = L/ℓ,

Ẽ(k_ℓ) = d/dk_ℓ [ ½⟨|ū_ℓ|²⟩ ] = -(ℓ²/L) d/dℓ[ ½⟨|ū_ℓ|²⟩ ].

Unlike cumulative_energy (the cumulative quantity, Eq. 15), this is the spectral density comparable to a Fourier energy spectrum. L is the region length: pass the domain size for the Sadek–Aluie convention k_ℓ = L/ℓ; the default L = 1 gives the FlowSieve convention k_ℓ = 1/ℓ. scales need not be uniform. Returns the filtering wavenumbers k_ℓ and the density Ẽ per scale.

References

• Sadek & Aluie (2018), Phys. Rev. Fluids 3, 124610.

spectral_density!(g, C, k) -> g

In-place spectral_density: writes the non-uniform finite-difference derivative of C w.r.t. k into the caller-supplied g (central in the interior, one-sided at the ends). Fills zeros for fewer than two points.

spectral_density(C, k) -> dC/dk

Non-uniform finite-difference derivative of cumulative values C w.r.t. k (central in the interior, one-sided at the ends). Returns zeros for fewer than two points.

tau_decomposition!(ws::TauWorkspace, u, v, grid, kernel, scale; filter_plan=nothing, ...) -> (; L, C, R)

In-place tau_decomposition. Writes into ws and returns views of its component buffers, so the result is valid until the next call on the same workspace. Supplying filter_plan as well makes a repeated decomposition allocation-free.

tau_decomposition(u, v, grid, kernel, scale; backend=AutoBackend(), mask_strategy=Deformable())
    -> (; L, C, R)

Split the 2D subfilter-scale stress τ_ij = ⟨u_i u_j⟩ - ū_i ū_j into Leonard, Cross, and Reynolds contributions (Germano 1992, JFM 238, using generalized central moments so each piece is individually Galilean-invariant). With ū = G * u the filtered velocity and u' = u - ū the residual, and the generalized second moment M(f, g) = (fg)‾ - f̄ ḡ:

• Leonard L_ij = M(ū_i, ū_j) (resolved–resolved),
• Cross C_ij = M(ū_i, u'_j) + M(u'_i, ū_j),
• Reynolds R_ij = M(u'_i, u'_j) (subfilter–subfilter; backscatter),

with L + C + R = τ exactly. Returns a named tuple of named tuples, each holding the symmetric 2D components (; xx, xy, yy) as arrays.

tau_decomposition(u, v, grid::StructuredGrid{<:SphericalGeometry}, kernel, scale; ...) -> (; L, C, R)

Spherical counterpart of the Cartesian method above: like compute_Π!'s spherical branch, the Leonard/Cross/Reynolds moments are formed in PLANETARY-CARTESIAN coordinates (so filtering commutes with the moment/residual operations, Aluie 2019), then each of L, C, R's resulting 3×3 symmetric tensor is rotated back to the local (east, north) frame at every grid point. L+C+R = τ still holds exactly (the rotation is linear). Returns the same (; L, C, R) shape as the Cartesian method — local (; xx, xy, yy) (≡ east-east/east-north/north-north) components.

tracer_variance_flux!(Πθ, ws, u, v, θ, grid, kernel, scale; filter_plan=nothing, deriv_plan=nothing, ...) -> Πθ

In-place tracer_variance_flux. With ws, filter_plan and deriv_plan all supplied, a repeated evaluation — over timesteps or scales — allocates nothing.

tracer_variance_flux(u, v, w, θ, grid::StructuredGrid{Cartesian,T,3}, kernel, scale; backend=AutoBackend(), mask_strategy=Deformable())
    -> Πθ

True three-dimensional analog of the 2D tracer_variance_flux above: the subfilter tracer flux gets a genuine vertical component τ_z = ⟨wθ⟩ - w̄θ̄, contracted against the resolved 3D tracer gradient ∂_j θ̄ (all three components, including the real vertical derivative ∂θ̄/∂z).

tracer_variance_flux(u, v, w, θ, grid::StructuredGrid{<:SphericalGeometry,T,3}, kernel, scale; ...) -> Πθ

Volumetric spherical shell (lon, lat, radius): the 3D counterpart of the spherical 2D method, keeping the radial component of both the subfilter tracer flux and the resolved gradient. Velocities are rotated to planetary Cartesian for filtering and the filtered flux is rotated back to local (east, north, radial), the same convention the true-3D compute_Π! uses.

tracer_variance_flux(u, v, θ, grid, kernel, scale; backend=AutoBackend(), mask_strategy=Deformable())
    -> Πθ

Cross-scale flux of the tracer variance ½⟨θ'²⟩ at filter scale ℓ (the scalar analog of the kinetic energy flux Π; Aluie & Eyink):

Πθ(x) = -∂_j θ̄ · τ_j(u, θ),   τ_j = ⟨u_j θ⟩ - ū_j θ̄  (the subfilter tracer flux),

with the same sign convention as compute_Π!: Πθ > 0 is a forward cascade of tracer variance toward small scales, Πθ < 0 an inverse cascade.

Taking θ to be the buoyancy b = -g ρ'/ρ₀ makes this the cross-scale transfer of buoyancy variance (the available-potential-energy-related transfer). Unlike the full Lees & Aluie (2019) baropycnal work — which additionally requires the pressure field — this needs only (u, v, θ).

Cartesian and spherical, on a 2D grid; the true-3D Cartesian method is below. ddx!/ddy! supply the physical gradient in either geometry.

tracer_variance_flux(u, v, θ, grid::StructuredGrid{<:SphericalGeometry}, kernel, scale; ...) -> Πθ

Spherical form of the tracer-variance flux. τ_j = ⟨u_j θ⟩ - ū_j θ̄ is a vector, so — exactly as in compute_Π! and tau_decomposition — the velocity is rotated to planetary Cartesian before filtering (Aluie 2019 commutativity: component-wise filtering of a local east/north pair is not a filtered vector, since the local basis turns from point to point), and the filtered flux is rotated back to the local east/north frame to contract against ∂_j θ̄. The scalar θ needs no rotation. The radial component of τ is dropped, matching this 2-D shell's dropping of radial derivatives.

DEFAULT_CACHE_BYTE_BUDGET

Default byte budget for AutoCache's size check on a real-space footprint cache (256 MiB).

AbstractCacheStrategy

Whether a real-space footprint over a genuinely nonuniform axis (StructuredGrid with a Vector axis, CurvilinearGrid, or ND with a non-Range axis) stores its full per-point neighbour list, or recomputes it on the fly at apply time. The per-point neighbour/weight computation itself is always the same either way (there is no shared translation-invariant table to exploit on a nonuniform axis, unlike the FilterFootprint fast path) — this only controls whether that computation's RESULT is kept in memory for reuse across separate future filter_apply! calls, or re-derived each time.

AbstractFilterMethod

How the convolution is evaluated: RealSpace (physical-space footprint, any grid/mask) or Spectral (transform-space multiply — FFT for uniform periodic Cartesian, spherical harmonics for the uniform sphere, NUFFT/NUFSHT for scattered points; provided by extensions).

AbstractFilterPlan

A prebuilt filter (grid + kernel + scale + mask strategy + backend) that can be applied to many fields without redoing setup. Physical-space backends precompute a FilterFootprint; the spectral extensions (FFTW/FINUFFT/SHT) hold cached transform plans. Declared here rather than alongside PhysicalFilterPlan further down so that filter_field!'s filter_plan::Union{Nothing, AbstractFilterPlan} keyword annotation, just below, can name it.

AbstractMaskStrategy

How masked (inactive) cells enter the filter normalization.

AlwaysCache <: AbstractCacheStrategy

Force building the full per-point neighbour-list cache regardless of cache_byte_budget — for a caller who knows more memory is available than the conservative default budget assumes.

AutoCache <: AbstractCacheStrategy

Build and store the full per-point neighbour-list cache only if its estimated size is under cache_byte_budget (default DEFAULT_CACHE_BYTE_BUDGET); otherwise fall back to recomputing neighbours on the fly at apply time. This is the only cache-strategy knob most callers ever need — it caches whenever doing so is affordable, which is strictly better than not caching whenever a plan will be reused across more than one filter_apply! call.

AutoMethod <: AbstractFilterMethod

Pick the engine from real capability, the same contract AutoCache and AutoBackend follow: Spectral only where a transform is available AND exact for this grid, otherwise RealSpace.

Selects Spectral only where every axis is periodic and uniform and a transform backend is loaded — the case where a periodic transform is the exact filter. Otherwise RealSpace.

Deformable <: AbstractMaskStrategy

Masked cells are excluded from BOTH numerator and denominator, so the kernel is renormalized over the locally-included area only ("deformable kernel"). Excluded cells are genuinely dropped, but the kernel becomes inhomogeneous near a mask boundary (breaks the strict commutation theorems).

FilterFootprint{T}

Precomputed convolution footprint for a structured grid + kernel + scale. The in-support neighbour offsets (di, dj) and their geometric weights w = kernel_weight(distance) * cell_area are stored in a flat CSR-like layout, grouped into axis-2 (y) bands (ptr[b]:ptr[b+1]-1). For Cartesian grids the footprint is translation-invariant → a single band; for a spherical grid it is invariant in x (longitude) → one band per y (latitude) value. The weights are mask-independent (geometry only); masking is applied when the footprint is convolved with a field.

FilterFootprintND{N, T}

General N-dimensional footprint: in-support neighbour offsets (NTuple{N,Int}) and their geometric weights w = kernel_weight(distance) · cell_measure. Used for 1D and 3D (Cartesian, translation-invariant ⇒ a single offset set); the 2D path uses the optimized per-row FilterFootprint.

NDScatteredCache{N, T}

The full per-target-point neighbour list for an NDScatteredFilterPlan (absolute neighbour multi-indices + weights), built only when the plan's AbstractCacheStrategy calls for it.

NDScatteredFilterPlan{N, T, K}

N-D (1D or 3D) analog of ScatteredFilterPlan: compact scalar metadata (per-axis window limits, periodicity/period, geometry flag) for when at least one of the N axes is a plain AbstractVector (no type-level uniformity proof) — no translation-invariance assumption, correct for any spacing pattern. cache holds the materialized NDScatteredCache only when the plan's cache strategy decided to build it, nothing otherwise (apply-time recomputation).

NeverCache <: AbstractCacheStrategy

Force recomputing neighbours on the fly at every filter_apply!/filter_apply_batch! call, never storing the cache — for a genuine, known memory ceiling AutoCache's budget check doesn't already account for (e.g. a GPU's device memory budget, or deliberately running many large plans concurrently). Not a speed/memory preference: for any workflow that reuses a plan across more than one call, caching is strictly better whenever it fits in memory, so this should be reached for only when a specific external memory constraint is known, not by default.

NodeFilterPlan{T}

Real-space filter footprint for a node set: per-node CSR neighbour blocks with their geometric weights w = kernel_weight(d) · control volume.

A node set has no axes, so there is no index window to bound a search with and the neighbourhood cannot be re-derived per apply the way a structured grid's can. It is found once, at plan time, and stored — which also means the search is paid once per plan rather than once per field, and a single compute_Π! makes six to nine applies against one plan.

The node itself is included. Connectivity.neighbors_within excludes it, matching stencil semantics where a cell is not its own neighbour, but a filter's zero offset is a genuine contribution.

Physical-space plan: a precomputed footprint reused across all longitudes, fields, and layers — for EVERY backend, not just serial. kernel/scale are retained only so the cached-footprint path can still call each backend's row-parallel hook (which takes them positionally); they're not used to rebuild the footprint once footprint is already built.

PrefixSumTopHatPlan{T,VT,MT,DT,BT,WX,WY}

Exact O(N·dj_lim) top-hat footprint for a rectilinear 2D StructuredGrid — see the section comment above for the derivation. Holds the separable measure factors (wx,wy), the extended axis-1 coordinate array the two-pointer sweep walks (tripled when axis 1 is periodic, so a wrapped support interval is still one contiguous run), and preallocated scratch so filter_apply! allocates nothing.

prefix_den (Deformable masking) is the prefix sum of mask·wx — mask-only, so it is built ONCE at plan-build time, mirroring SeparableGaussianFootprint/FFTWFilterPlan's invrenorm convention. ZeroFill's denominator is mask-independent and needs only the 1-D prefix_wx.

Physical-space direct-sum convolution (works on any grid, mask, and geometry).

ScatteredCache{T}

The full per-target-point neighbour list for a ScatteredFilterPlan (absolute neighbour indices + weights), built only when the plan's AbstractCacheStrategy calls for it. Same CSR-style layout as before: target t = i + (j-1)*Nx, entries ptr[t]:ptr[t+1]-1.

ScatteredFilterPlan{T,K}

Real-space footprint for a genuinely nonuniform 2D StructuredGrid axis, or a CurvilinearGrid (exactly the periodic_x = periodic_y = false case of the same candidate-window/distance-gate computation). FilterFootprint is a translation-invariant cache — the SAME index offset (and its weight) is reused for every target point — which is only valid on a uniform axis; here that assumption is false (offset +3 means a different physical displacement depending on where you start), so there is no way to share one offset/weight set across points. That does NOT mean the result must be stored, though: since the search window (di_lim/dj_lim) is already a global scalar bound (not per-point), the exact same candidate enumeration + distance/kernel_weight gate that determines a point's neighbours can be re-run identically at apply time from these few scalars alone — cache holds the materialized ScatteredCache only when the plan's cache strategy decided to build it (see AbstractCacheStrategy), nothing otherwise (apply-time recomputation).

SeparableGaussianFootprint{T}

Precomputed 1D Gaussian weight vectors (gx,gy) plus preallocated scratch buffers for the row-pass/column-pass separable convolution. invrenorm (Deformable masking only) is the precomputed reciprocal local kernel mass over active cells — the SAME separable machinery run once, at plan-build time, on Float.(mask) instead of field, mirroring FFTWFilterPlan/SHTFilterPlan's established invrenorm pattern. Nx_profile/Ny_profile (ZeroFill masking) are the mask-INDEPENDENT denominator profiles: Σ w over valid (in-bounds/periodic) offsets is itself separable into one Nx-length and one Ny-length vector, since which offsets are valid depends only on i (resp. j) and periodicity, never on the mask.

Note: unlike the disk-truncated (d <= rad) non-separable footprint, this truncates each axis independently at the SAME per-axis rad (the Gaussian's own 1D marginal decays at the identical rate kernel_radius was derived from) — a square window, not a disk. The Gaussian has no true hard support (only a numerical truncation tolerance), so this is an equally valid truncation, just a different shape — matched against RealSpace within a measured tolerance, not asserted bit-identical.

SeparableGaussianFootprintND{N,T}

The separable Gaussian in N dimensions: one weight table per axis, applied as N successive 1-D passes.

This is the same factorization the 2-D path uses, and the reason it matters grows with N. A FilterFootprintND enumerates the whole ∏(2wᵈ+1) box per point; N passes cost ∑(2wᵈ+1). At w = 20 in 3-D that is 68,921 multiply-adds per point against 123.

Weight tables follow _sepw: a vector where the axis is uniform, an Nᵈ × (2wᵈ+1) matrix where it is stretched.

Spectral <: AbstractFilterMethod

Transform-space filtering (kernel applied as a multiply on the transformed field). Requires a spectral extension and a compatible grid (e.g. using FFTW for a uniform, periodic Cartesian grid).

ZeroFill <: AbstractMaskStrategy

Excluded cells are treated as zero-valued: they contribute to the denominator (kernel weight) but zero to the numerator. The kernel is homogeneous (same shape everywhere), which preserves domain averages and commutation with derivatives (the Storer 2022 / Aluie 2019 "fixed kernel" mode).

_build_footprint_curvilinear(grid, kernel, scale; cache_strategy=AutoCache(), cache_byte_budget=DEFAULT_CACHE_BYTE_BUDGET) -> ScatteredFilterPlan

Compact plan (and, if cache_strategy calls for it, the full per-point cache) for a FlowGeometries.Grids.CurvilinearGrid. Enumeration goes through the grid's own ball query, as it does for a StructuredGrid; what differs is the SIZE ESTIMATE the cache budget is checked against. Connectivity.metric_window bounds a window from per-axis spacing, and a curvilinear mesh has no separable axes to bound with, so the estimate comes from the smallest adjacent-node spacing in each index direction instead. It feeds no computation — only the AutoCache decision — so being loose costs a cache that would have fit, never a wrong answer.

_build_footprint_scattered(grid, kernel, scale; cache_strategy=AutoCache(), cache_byte_budget=DEFAULT_CACHE_BYTE_BUDGET) -> ScatteredFilterPlan

Build the compact nonuniform-axis plan (O(1) scalar metadata) and, only if cache_strategy calls for it, the full per-point ScatteredCache — correct for any spacing pattern (Cartesian or spherical, uniform or not), since it never assumes translation invariance. The search window comes from _scattered_window_bounds, i.e. the grid's own Connectivity.metric_window; the exact d <= rad check still gates inclusion, so a loose bound only costs iterations, never a missed cell.

_build_prefixsum_tophat(grid, kernel::TopHatKernel, scale; mask_strategy=Deformable(), kwargs...)

Build the exact O(N·dj_lim) prefix-sum top-hat plan. kwargs... absorbs (and ignores) cache_strategy/cache_byte_budget — there is no per-point neighbour list here to cache at all.

_build_separable_gaussian_footprint(grid, kernel::GaussianKernel, scale; mask_strategy=Deformable(), kwargs...) -> SeparableGaussianFootprint

Build the per-axis weight tables, the mask-dependent normalization data for mask_strategy, and the preallocated scratch buffers for the separable Gaussian path. kwargs... absorbs (and ignores) cache_strategy/cache_byte_budget — there is no per-point neighbour list to cache here at all, so those knobs (which only govern the scattered path) don't apply.

apply_footprint_row!(out, field, grid, fp, strategy, periodic_x, periodic_y, j)

Fill output row j (out[:, j]) from a precomputed footprint. Rows are independent (each writes a disjoint column of the column-major output), so this is the unit of parallelism for the threaded / distributed backends. Callers must fill!(out, 0) first (masked cells are left untouched here).

_nd_scattered_window_bounds(grid, rad) -> (lim, periodic, period, is_cartesian)

Compact per-axis scalar window-bound derivation for the nonuniform N-D (1D/3D) case — the same values an NDScatteredFilterPlan stores. The window itself is Connectivity.metric_window, read from the grid's cached axis statistics in O(1); nothing here scans the grid.

_rectilinear_measure_factors(grid) -> (wx, wy)

The grid's own per-axis measure factors, so that measure[i, j] == wx[i] * wy[j] exactly. Read from the stored SeparableMeasure rather than rebuilt here: rebuilding has to reproduce every convention the measure was constructed under, and a degenerate direction is where that fails — a single-latitude grid measures arc length R·Δλ and a single-longitude one R·Δφ, neither of which is the R²cosφ area form.

_scattered_window_bounds(grid::StructuredGrid{...,2}, rad) -> (di_lim, dj_lim, periodic_x, periodic_y, x_period, y_period, is_cartesian)

The widest window the grid's own ball query will scan, taken from Connectivity.metric_window rather than re-derived here. On a rectilinear grid the window depends on the row, not the column, so one evaluation per row covers the grid.

_sep_serial(f, indices)

The default pass driver: apply f to every index in order. Every output point of a pass is independent, so a backend can substitute a parallel driver of the same shape and get a bit-identical result — only the barrier BETWEEN passes is required.

_separable_axis_weights(x, lim, periodic, period, α, scale) -> AbstractVecOrMat

The Gaussian's per-axis weight table: exp(-α·(Δx/ℓ)²) over the stencil, in the layout _sepw reads.

Uniform axis: the displacement is ddi·Δ wherever the stencil sits, so one vector serves every position. Stretched axis: the displacement depends on the position too, so the table gains a position axis — (2·lim+1) × N, which is O(N·lim) against the O(N·lim²) of a per-point neighbour cache, and leaves the apply at O(N·lim) instead of O(N·lim²).

lim comes from the SMALLEST gap on the axis, so on a stretched axis a coarse region's stencil is wider than it needs to be; those slots hold exact zeros rather than being trimmed, which keeps the inner loop's bounds static. A periodic displacement carries the image offset, matching the tiling convention the scattered engine uses.

_sepw(g, i, k) -> T

Weight of stencil slot k at axis position i.

The Gaussian factorizes on ANY rectilinear grid — exp(-α(Δx²+Δy²)/ℓ²) = Gx(Δx)·Gy(Δy) needs no constant spacing — but on a uniform axis Gx depends on the OFFSET alone, while on a stretched one it depends on the position too. Both are the same convolution with a different weight table, so the two are one code path distinguished by the table's rank: a vector is shared across positions and a matrix is (2·lim+1) × N, column-major so each position's stencil is contiguous.

apply_footprint!(out, field, grid, fp::PrefixSumTopHatPlan, strategy, periodic_x, periodic_y) -> out

apply_footprint!-shaped entry point for the prefix-sum plan, so the generic whole-grid convolve name works uniformly across every footprint type. periodic_x/periodic_y are accepted for interface compatibility but IGNORED: unlike the offset-based footprints, this plan captured its periodicity (and the matching extended-coordinate layout) from the grid at build time, and cannot honour a different choice at apply time.

apply_footprint!(out, field, grid, fp::ScatteredFilterPlan, strategy, periodic_x, periodic_y)

Whole-grid convolve using a ScatteredFilterPlan (the nonuniform-axis/curvilinear fallback). periodic_x/periodic_y are accepted only for a uniform call signature with the FilterFootprint method above — periodicity for this footprint kind lives in fp itself, not these arguments.

apply_footprint!(out, field, grid, fp, strategy, periodic_x, periodic_y)

Convolve field with a precomputed fp into out, applying the mask strategy. out and field are 2D (a single layer). The masking branch specializes on the strategy type.

apply_footprint!(out, field, grid, fp::NodeFilterPlan, strategy) -> out

Weighted mean over each node's stored neighbourhood, with the same two mask conventions the structured engines use: ZeroFill keeps a masked neighbour in the denominator and contributes nothing for it, Deformable drops it from both.

apply_footprint_nd_batch!(outs, fields, grid, fp, strategy) -> outs

Batched point-indexed apply over the whole grid. The batch shares each point's neighbour enumeration across all K fields, so the enumeration is paid once rather than K times.

apply_footprint_nd_batch_over!(outs, fields, grid, fp, strategy, indices) -> outs

The batched apply restricted to indices. Each output point depends only on its own neighbourhood, so a parallel backend can hand disjoint index blocks to different tasks and get the serial answer. outs must already be zeroed — the caller owns that, since a block only writes its own points.

apply_footprint_row!(out, field, grid, fp::ScatteredFilterPlan, strategy, periodic_x, periodic_y, j)

Fill output row j from a ScatteredFilterPlan: if fp.cache !== nothing, read the precomputed per-point neighbour list (absolute (ii,jj) indices, periodic wrap already resolved at build time); otherwise recompute each point's neighbours/weights on the fly from fp's compact scalar metadata. Both branches enumerate candidates through _scattered_foldl, so they are bit-identical by construction rather than by convention. The accumulator is threaded through the fold's return value rather than captured and mutated, which is what keeps the streaming branch free of per-iteration allocation (verified by @allocated tests) without a second copy of the loop.

apply_prefixsum_tophat!(out, field, grid, fp, strategy) -> out

Whole-grid exact top-hat convolve: one O(N) prefix-sum pass over the field, then one O(N·dj_lim) two-pointer sweep.

apply_prefixsum_tophat_row!(out, grid, fp, strategy, j) -> out

Fill output row j from the (already current) prefix sums. For each row offset in the band, sweeps axis 1 with two MONOTONE pointers — O(1) amortized per point, not a per-point binary search — so the whole row costs O(Nx·dj_lim). Touches only row j of out/fp.den, so rows may run concurrently.

apply_separable_gaussian!(out, field, grid, fp::SeparableGaussianFootprint, strategy) -> out

Apply the separable Gaussian fast path: masked_input = mask .* field (the SAME numerator input for both mask strategies — see the struct docstring), one shared row-pass/column-pass convolution, then divide by whichever mask-strategy-specific denominator fp holds.

apply_separable_gaussian_nd!(out, field, grid, fp, strategy, driver = _sep_serial)

Run the N separable passes and the pointwise normalization. driver supplies the per-pass index sweep — see _sep_serial; a threaded backend passes its own and gets the same answer, since every point within a pass is independent and the passes themselves stay ordered.

apply_slice_serial!(out, field, plan) -> out

One slice, forced down the serial engine regardless of the backend recorded in plan. The slice loop owns the parallelism; see filter_slices!.

build_footprint(grid, kernel, scale; kwargs...) -> ScatteredFilterPlan

General path: at least one axis is a plain (non-Range) AbstractVector, which makes no type-level uniformity guarantee — its values might happen to be evenly spaced, but nothing proves it, so no assumption is made and the always-correct per-point plan is built instead. (Less specific than the method above, so Julia only reaches this one when the fast method's constraint doesn't match.)

build_footprint(grid::StructuredGrid{Cartesian,T,2}, kernel::GaussianKernel, scale; kwargs...) -> SeparableGaussianFootprint

Separability does not require constant spacing: exp(-α(Δx²+Δy²)/ℓ²) factorizes on any rectilinear grid, and a stretched axis only makes the per-axis weight depend on position as well as offset — see _separable_axis_weights. So a stretched Cartesian grid gets the same two-pass O(N·(wx+wy)) convolution rather than falling to the O(N·wx·wy) scattered engine, which for a Gaussian at w = 20 is a factor (2w+1)/2 in operations and a much larger one in per-operation cost.

The Range-axis method above is strictly more specific and resolves what would otherwise be an ambiguity with the generic Range-axis method; both build the same footprint.

build_footprint(grid::StructuredGrid{<:AbstractGeometry,T,2}, kernel::TopHatKernel, scale; kwargs...) -> PrefixSumTopHatPlan

Exact prefix-sum top-hat path for any rectilinear 2D StructuredGrid (see the section comment above). More specific than the generic 2D methods (constrained on kernel type), so Julia selects it whenever kernel isa TopHatKernel.

build_footprint(grid, kernel, scale) -> FilterFootprint

Fast path — real multiple dispatch, not a runtime check: both axes are AbstractRange, a compile-time proof of constant spacing, so the footprint is genuinely translation-invariant and can be shared via a single (Cartesian) or per-latitude-band (spherical) offset/weight cache. Spacing is read via step(...) directly from the axis that's already proven uniform by its type — not from the geometry's separately-stored dx/dy scalar, so there's no possibility of the two disagreeing.

build_footprint(grid::StructuredGrid{Cartesian,T,2,TP,<:Tuple{AbstractRange,AbstractRange}}, kernel::GaussianKernel, scale; kwargs...) -> SeparableGaussianFootprint

Fast path for a GaussianKernel on a uniform (Range-axis) Cartesian grid — see the "Separable Gaussian fast path" section above. More specific than the generic Range-axis method (constrained on kernel type too), so Julia picks this one whenever kernel isa GaussianKernel.

build_footprint(grid::CurvilinearGrid, kernel, scale; kwargs...) -> ScatteredFilterPlan

Real-space direct-sum footprint for a curvilinear grid (see _build_footprint_curvilinear).

build_footprint(grid::UnstructuredGrid, kernel, scale; kwargs...) -> NodeFilterPlan

Real-space footprint for a node set, from Connectivity.fold_within — the grid's own metric ball query, so the neighbourhood honours the geometry's distance and any periodic wrap without this package re-deriving either, and the fold hands back each neighbour's distance rather than making the weight recompute it.

The sweep visits every node, which is where a spatial index pays for itself: Connectivity.default_sweep_topology builds one when NearestNeighbors is loaded, taking the build from O(n²) to O(n log n), and returns the unindexed topology (same rows, linear scan) when it is not. Either way it is paid ONCE per plan and reused by every filter_apply!, including the six to nine a single compute_Π! makes.

filter_apply!(out, field, plan) -> out

Apply a prebuilt plan_filter to a single 2D field, dispatching to whichever backend the plan was built for — the footprint is ALWAYS the one cached in plan, never rebuilt here, for every backend (serial, threaded, distributed, GPU, MPI).

filter_apply_batch!(outs, fields, plan::AbstractFilterPlan) -> outs

Apply plan to every field in fields, writing into the matching entry of outs, deriving each target point's neighbour list/weight exactly ONCE and reusing it across the whole batch — not once per field. outs/fields must be equal-length, matching-shape collections of arrays: an NTuple{K,V} (single concrete array type V) for a compile-time-known batch size (fastest — see _batch_zeros), or an AbstractVector for a runtime-determined batch size.

filter_field!(out, field, grid, kernel, scale; mask_strategy=Deformable(), filter_plan=nothing, backend=AutoBackend())

Filter a field on a grid using kernel at characteristic full width scale (ℓ), writing the result to out (returned).

Keyword Arguments

• mask_strategy::AbstractMaskStrategy=Deformable(): masking strategy — ZeroFill() (excluded cells count in the denominator as zero; homogeneous kernel) or Deformable() (excluded cells dropped from numerator and denominator; renormalized over the locally-included area).
• filter_plan::Union{Nothing,AbstractFilterPlan}=nothing: a prebuilt plan_filter result to reuse instead of building one from scratch — the zero-(re)allocation entry point for a repeated sweep (many timesteps/fields over the same grid/kernel/scale). When supplied, mask_strategy/ backend/method are ignored (already baked into the plan); build it once with plan_filter and pass it here on every subsequent call.
• backend::AbstractExecutionBackend=AutoBackend(): execution backend (SerialBackend, ThreadedBackend, GPUBackend, …). Ignored when filter_plan is supplied.

For spherical grids the longitude footprint wraps only when the grid is periodic (isperiodic); distances use the great-circle (Haversine) metric.

Examples

geom = CartesianGeometry()
grid = StructuredGrid(geom, 0.0:1000.0:99_000.0, 0.0:1000.0:99_000.0, mask)
out = zeros(100, 100)
filter_field!(out, field, grid, TopHatKernel(), 5000.0; mask_strategy = Deformable())

filter_fields!(outs, fields, grid, kernel, scale; mask_strategy=Deformable(), backend=AutoBackend())

Filter several fields that share the same grid/kernel/scale, building the footprint/plan ONCE and applying it through filter_apply_batch!, so each target point's neighbours are enumerated once for the whole batch rather than once per field. outs and fields are indexable collections of matching arrays — a tuple of velocity components, or a vector of them.

filter_slices!(outs, fields, plans; backend = AutoBackend()) -> outs

Apply plans[t] to fields[t], writing outs[t], over a collection of independent slices.

This is a different parallel axis from filter_apply_batch!, which shares one grid across several fields: here each slice has its own grid, plan and point count, and slices share nothing, so there is no synchronization at all. Where a workload has many slices, this is the outermost race-free axis and the one that converts thread count into throughput — threading within one slice saturates once the slice is small enough that per-task overhead dominates its work.

Each slice runs serially inside, whatever backend its own plan carries: nesting a threaded apply under a threaded slice loop would have both levels claim the whole thread pool.

plan_filter(grid, kernel, scale; mask_strategy=Deformable(), backend=AutoBackend()) -> AbstractFilterPlan

Build a reusable filter plan: the footprint is precomputed ONCE regardless of backend (serial, threaded, distributed, GPU, or MPI) and reused across every subsequent filter_apply! call — no backend rebuilds it per call. Apply with filter_apply!(out, field, plan).

plan_filter(grid::UnstructuredGrid, kernel, scale; method = Spectral(), …)

A node set supports both methods, and both plan in O(n log n). Spectral() is the default: it is exact for a band-limited field and its per-apply cost is independent of the filter scale, where the real-space engine's grows with the neighbour count inside the ball. RealSpace() applies the kernel as written, with compact support; a transform's support is global.

prefixsum_fill_numerator!(fp, field, grid)

The single O(N) per-apply pass: refill fp.prefix_num (cumulative mask·field·wx per row) from the current field. Must run before any apply_prefixsum_tophat_row! call for that field — both the serial and threaded drivers guarantee this.

prefixsum_fill_numerator_row!(fp, field, grid, j)

Fill row j's numerator prefix scan. Writes only column j of fp.prefix_num, so rows are mutually independent and this may be run concurrently across j (the threaded backend does exactly that).

prepare_workspace(backend, grid, footprint) -> workspace

Backend hook run ONCE by plan_filter, whose result becomes the plan's stored workspace. The default returns the footprint unchanged; a backend that needs its own residency — the GPU's device buffers — returns something its apply step consumes directly, so no transfer is repeated per call.

slice_costs(plans) -> Vector{Int}

Per-slice work proxy: the number of target points each plan writes. A slice's cost grows at least linearly in this, so it is what a longest-first schedule should sort on.

AbstractFilterKernel

Abstract supertype for all filter kernels (TopHat, Gaussian, SharpSpectral, …).

GaussianKernel(; α = 6.0) <: AbstractFilterKernel

Real-space Gaussian filter G_ℓ(d) ∝ exp(-α (d/ℓ)²), with ℓ the full filter width.

• α = 6 (default) is the Pope/turbulence-literature convention: the Gaussian's second moment matches the top-hat box of width ℓ (σ² = ℓ²/12).
• α = 4 reproduces FlowSieve's default Gaussian (which also treats ℓ as a diameter), so GaussianKernel(; α = 4) is directly comparable to FlowSieve output.

SharpSpectralKernel <: AbstractFilterKernel

Sharp-spectral (brick-wall) filter: Ĝ_ℓ(k) = 1 for k ≤ k_c, else 0, with k_c = π/ℓ. Best applied in spectral space (FFTW / FINUFFT / spherical-harmonic extensions); the physical-space form below is a slowly-decaying sinc fallback.

TopHatKernel <: AbstractFilterKernel

Real-space top-hat (box) filter of full width ℓ: unit weight for d ≤ ℓ/2, zero otherwise.

kernel_radius(kernel::AbstractFilterKernel, ℓ::T) where {T<:AbstractFloat}

Distance beyond which the kernel weight is negligible, used to truncate physical-space convolution footprints.

kernel_weight(kernel::AbstractFilterKernel, d::T, ℓ::T) where {T<:AbstractFloat}

Evaluate the unnormalized kernel weight at distance d for filter width ℓ.

spectral_transfer(kernel, kmag::T, ℓ::T) where {T<:AbstractFloat}

Isotropic planar spectral transfer function Ĝ(|k|, ℓ): the factor by which a Fourier mode of physical wavenumber magnitude kmag (rad m⁻¹) is multiplied when filtering at width ℓ on a 2D Cartesian grid. Normalized so Ĝ(0, ℓ) = 1 (preserves the domain mean). Shared by the FFTW and FINUFFT backends (both 2D-Cartesian-only today). For the spherical-harmonic-degree analog used by the FastSphericalHarmonics/NUFSHT backends, see spectral_transfer_degree.

• GaussianKernel(α): Ĝ = exp(-k² ℓ² / (4α)) (the exact Fourier transform of exp(-α(r/ℓ)²)).
• SharpSpectralKernel: Ĝ = 1 for k ≤ π/ℓ, else 0.
• TopHatKernel: Ĝ = 2 J₁(kR)/(kR), R = ℓ/2 — the exact 2D (disk) Fourier transform of a top-hat (the "jinc" function, the circular-aperture analog of sinc). This oscillates and goes negative in k; that is the correct, exact behavior of a disk's Fourier transform, not an approximation error. This method is provided entirely by the SpecialFunctions weak dependency (CoarseGrainingEnergyFluxesSpecialFunctionsExt, for besselj1) — core has no method for TopHatKernel here (Julia disallows two modules defining the identical method signature, so a throwing core stub could never be replaced by the extension's real one); without using SpecialFunctions loaded, calling this is a MethodError with a registered hint pointing at the fix.

spectral_transfer_degree(::TopHatKernel, l::Integer, ℓ::T, R::T) where {T<:AbstractFloat}

Exact spherical-cap top-hat window function (Jekeli 1981's gravity-field averaging kernel; the sphere's analog of the planar top-hat's Bessel-J₁ transfer function): for a cap of angular radius θ0 = ℓ/(2R) (i.e. full physical width ℓ),

Ĝ_l = [P_{l-1}(x) - P_{l+1}(x)] / [(2l+1)(1 - x)]  ≡  (1 + x) P′_l(x) / (l(l+1)),   x = cosθ0

evaluated in the right-hand form, via the Legendre recurrences (n+1)P_{n+1}(x) = (2n+1)x Pₙ(x) - n P_{n-1}(x) and P′_{n+1}(x) = (2n+1)Pₙ(x) + P′_{n-1}(x) — no external dependency needed (unlike the planar case's Bessel J₁). Like the planar top-hat, this oscillates and goes negative in l; that is the exact, correct behavior of a spherical cap's harmonic content, not an artifact.

spectral_transfer_degree(kernel, l::Integer, ℓ::T, R::T) where {T<:AbstractFloat}

Spherical-harmonic-DEGREE-indexed transfer function Ĝ_l, used by the FastSphericalHarmonics/NUFSHT backends in place of spectral_transfer's continuous wavenumber kmag when a kernel's shape needs the discrete degree l itself, not just the Laplace–Beltrami eigenvalue k_l = √(l(l+1))/R. GaussianKernel/SharpSpectralKernel are smooth isotropic functions of k_l alone, so they simply delegate to spectral_transfer; TopHatKernel's spherical-cap window genuinely needs l.

StencilPlan(grid; order = 1, nodes = 3)

The finite-difference weights of every direction of grid, built once. Discretization.axis_stencils per axis; the derivative is then Discretization.derivative! reading a table it does not have to rebuild.

The weights depend only on the axis, its wrap period and the requested order — never on a field — so a caller taking many derivatives on one grid should build this once and pass it. Without it each call rebuilds an order-by-nodes table per axis sample, which is O(n) work and O(n·nodes) garbage against an O(nᴺ) apply: negligible on a large grid, several times the whole cost on a small one.

It also carries the degrade path's scratch, so a masked grid allocates nothing per call either. That scratch is written per cell, so a plan is one per task — the same contract as Connectivity.ball_scratch. Sharing one across concurrent tasks races; build one per task instead. (A threaded backend passed to apply_stencil! allocates its own set per chunk and ignores this one.)

compute_Π! builds one internally when its deriv_plan is nothing.

StencilPlan(axis::AbstractVector; order = 1, nodes = 3, period = nothing)

A one-direction plan over a bare axis, for the level-stack ddz! — there the vertical spacing is an argument rather than a grid axis, so there is no grid to take it from.

_dd!(∂f, f, grid, d) -> ∂f

Derivative of f with respect to distance along direction d. nodes = 3 for 2nd order on a stretched axis; ReduceInRun keeps the one-sided value at a mask edge, where the default writes zero.

ddx!(∂f∂x, f, grid[, plan]) -> ∂f∂x

Derivative of f with respect to distance along the Eastward/λ direction.

Pass a StencilPlan to reuse the weights across calls; without one they are rebuilt each time. The dimensionality is pinned per direction, so asking a grid for a derivative it has no axis for is a MethodError at the call rather than a bounds error inside the kernel.

ddy!(∂f∂y, f, grid[, plan]) -> ∂f∂y

Derivative of f with respect to distance along the Northward/φ direction — see ddx!.

ddz!(∂f∂z, f, grid[, plan]) -> ∂f∂z

Derivative of f with respect to distance along the third direction of a 3D grid, which supplies the axis — see ddx!. For a 3D field over a 2D grid, see the dz method below.

ddz!(∂f∂z, f, grid, dz[, plan]) -> ∂f∂z

Calculate spatial derivative of f in the vertical coordinate z, writing to ∂f∂z.

This is the level-stack case: a 3D field over a 2D grid (the same shape filter_field! accepts for a stack of levels). The grid describes only the horizontal, so it cannot supply the vertical spacing — dz is therefore an explicit argument rather than something read off the geometry. For a genuine 3D grid use the StructuredGrid{…,3} method, which takes its spacing from the grid's own third axis and handles nonuniform levels.

Since the vertical axis is not on the grid, its weights cannot come from a grid-built StencilPlan; build one over the axis instead — StencilPlan(range(0; step = dz, length = Nz)) — and pass it, or the table is rebuilt on every call at O(Nz).

plot_spectrum(res; which=:density) -> Figure

Plot the filtering spectrum from a CoarseGrainResult. which = :density plots the filtering spectral density Ẽ(k_ℓ) against filtering wavenumber k_ℓ (log x); which = :cumulative plots the cumulative coarse KE against scale ℓ (log–log). Provided by the CairoMakie package extension — run using CairoMakie to enable it.

plot_Π_map(res, scale_idx, grid; colormap=:balance, title=nothing) -> Figure

Heatmap of the cross-scale energy-flux map Π from a CoarseGrainResult at scale index scale_idx. Provided by the CairoMakie package extension — run using CairoMakie to enable it.