Diagnostics

AbstractSpectrumPolicy

What to do when a filtering spectral density is asked for with a kernel whose |Ĝ(k)|² is not monotone decreasing — StrictSpectrum, ForceSpectrum or NoSpectrum.

Sadek & Aluie (2018) eq. (21) guarantees Ẽ(k_ℓ) ≥ 0 only when d|Ĝ(k)|²/dk ≤ 0. That condition is sufficient, not necessary, and where it fails it tends to fail narrowly: the default TopHatKernel's |Ĝ|² falls to zero at kℓ ≈ 7.66 and climbs back to only 0.0175 at kℓ ≈ 10.27, so the violation sits in the far sub-filter tail at under 2% of the DC value while the rest of the curve is usable. Hence three settings rather than a veto: the safe reading stays the default, and the other one stays reachable.

Π and the cumulative energy carry no such condition and are unaffected by this choice.

EnergyWorkspace(grid; has_w = false)
EnergyWorkspace(grid, batch_size; has_w = false)

Scratch for the diagnostics that only need FILTERED VELOCITY — cumulative_energy! and, through it, filtering_spectrum. Those compute E(ℓ) = ½⟨|ū_ℓ|²⟩, which reads no stress, no strain and no quadratic product, so two buffers serve (three with a vertical component).

Interchangeable with ΠWorkspace wherever only the velocity buffers are used, so a sweep that already holds a flux workspace passes it straight through instead of allocating a second one.

EnstrophyFluxWorkspace(grid)

Scratch for enstrophy_flux!: the vorticity plus the tracer-flux scratch it is fed into.

Favre3DWorkspace(grid)

Scratch for the true-3-D compressible_flux!: the velocity triple, its unweighted filter, the mass flux and the Favre velocity, six product buffers, the six components of the Favre stress and of the resolved strain, the three pressure gradients and the three output fields.

FavreWorkspace(grid)

Scratch for compressible_flux!: the filtered density and pressure, the Favre velocities, the unweighted velocities, the three Favre stress components, the two unweighted mass-flux components, the four velocity gradients, the two pressure gradients, and the three output fields.

ForceSpectrum <: AbstractSpectrumPolicy

Compute the density regardless, warning once per session. The caller owns checking its sign — the right setting when the kernel's non-monotone band sits outside the range of scales being interpreted.

NoSpectrum <: AbstractSpectrumPolicy

Skip the density entirely; the field is filled with NaN rather than a number that would read as computed. Π and the cumulative energy are still produced.

PiDecomposed3DWorkspace(grid)

Scratch for the true-3D compute_Π_decomposed!: the three divergent components, the two filtered velocity triples, six product buffers, the six independent components of each of the three stress tensors and the two strain tensors, the four flux fields, and one derivative temporary.

Tensor components are stored in the (xx, xy, xz, yy, yz, zz) order shared with the spherical tau_decomposition!, so one loop over that order builds a whole tensor and feeds a single batched apply.

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.

PiStrainWorkspace(grid)

Scratch for compute_Π_strain_convergence!: the two filtered velocities, the three stress components, the four velocity-gradient components, the two rotation invariants and the three flux fields, plus the two shared product buffers.

SphericalAnalysis

Marker that ws.ux/uy/uz already hold the planetary-Cartesian velocity components for this sweep.

The spherical flux path works in planetary Cartesian coordinates, and the rotation into them is a function of (u, v, w, grid) alone — a coords lookup and a trigonometric basis change at every point, with no dependence on the filter scale. Those three buffers are not written again during a scale, so one rotation serves the whole sweep.

SphericalFavreWorkspace(grid)

Scratch for the spherical compressible_flux!: the planetary velocity, the mass flux and the Favre velocity in those coordinates, six product buffers, the six components of the Favre stress, the local velocity pairs, the resolved strain, the two pressure gradients and the three output fields.

SphericalPiDecomposed3DWorkspace(grid)

Scratch for the spherical volumetric compute_Π_decomposed!: the two planetary-Cartesian velocity triples and their filtered forms, six product buffers, the six components of each of the three stresses and the two strains, one local triple, and the four flux fields.

SphericalPiDecomposedWorkspace(grid)

Scratch for the spherical compute_Π_decomposed!: the two planetary-Cartesian velocity triples and their filtered forms, six product buffers, the six components of each of the three stresses, the local filtered velocity pair of each part, the two strains, and the four flux fields.

SphericalTracerFluxWorkspace(grid)

Scratch for the spherical tracer_variance_flux!: the planetary velocity, its filter, the three velocity–tracer products and their filter, the filtered tracer, the local subfilter flux pair and the resolved tracer gradient.

StrictSpectrum <: AbstractSpectrumPolicy

Refuse to produce a spectral density for a kernel that fails Kernels.transfer_monotone. The default: either a density guaranteed non-negative, or an error naming the alternatives.

Sym3TauWorkspace(grid)

Scratch for the three-component tau_decomposition!: a velocity triple and its residual, their filtered and double-filtered forms, six product buffers, and the six independent components of each of L, C and R.

Two paths take this shape. A spherical tangent split carries the planetary-Cartesian velocity and keeps the first three local components of each tensor; a true-3-D split carries (u, v, w) and keeps all six.

Where a rotation back to the local frame happens it is done in place into the leading slots of each tensor, since it reads a point's six planetary components and writes that point's local ones — so the returned components alias those slots and no separate output buffers exist.

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.

TracerFlux3DWorkspace(grid)

Scratch for the true-3-D tracer_variance_flux!: the velocity triple, its filter, the three velocity–tracer products and their filter, the local subfilter flux, the filtered tracer and its three-component resolved gradient.

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.

ΠWorkspace(grid; has_w = false)
ΠWorkspace(grid, batch_size; has_w = false)

Scratch for a flux computation. Passing batch_size sizes every buffer as (spatial..., batch...) so a whole batch of slices is held at once; the elementwise algebra then broadcasts over the trailing axes unchanged, and a filter apply can cover the batch in one pass instead of one per slice.

Only the buffers the requested configuration can reach are allocated, because which ones those are is fixed by the grid and by whether a vertical component is supplied — not discovered at run time:

has_w is what selects between the first two, and it must be given at construction because the buffers have to exist before the first call. A workspace built without them refuses a w rather than returning a wrong answer, with a message naming the fix.

Spherical and true-3-D grids always carry the vertical set: the spherical branch works in planetary Cartesian components, so it has three velocity components and six products whether or not the caller supplied a vertical velocity.

_area_mean(field, grid, total_area) -> T

Area-weighted mean of field over the ACTIVE cells, sharing active_area's denominator so every spatial average in this module is normalized the same way.

_fill_planetary!(ws, u, v, w, grid) -> nothing

Rotate the local (east, north, up) velocity into planetary Cartesian components in ws.ux/uy/uz. Inactive cells are zeroed rather than skipped, so the filter sees a genuine zero there.

_pair_moments!(Maa, Mab, Mbb, a, b, ā, b̄, pr, plan)

The three generalized second moments of a velocity pair, in whatever frame the pair is given:

M(a,a) ,   M(a,b) + M(b,a) ,   M(b,b) ,   M(x,y)_ij = (x_i y_j)‾ − x̄_i ȳ_j .

ā, b̄ are the already filtered triples and pr six scratch buffers, so each moment costs one batched apply of six products. The cross moment sums both orderings before filtering, which the linearity of the filter permits, so it costs one product per component.

_rotate_moments_to_local!(tensors, grid, geo, ::Val{NC})

Rotate each symmetric planetary-Cartesian tensor into the local frame, in place. NC is how many local components the caller keeps: three for a tangent split, six for a volume, both written over the leading slots of the tensor they came from — safe, since all six planetary values at a point are read before any local one is written.

_sph_pair_moments!(Maa, Mab, Mbb, a, b, ā, b̄, pr, plan, grid, geo)

The three generalized second moments of a planetary-Cartesian velocity pair,

M(a,a) ,   M(a,b) + M(b,a) ,   M(b,b) ,   M(x,y)_ij = (x_i y_j)‾ − x̄_i ȳ_j ,

each rotated from the planetary 3×3 into the local (east, north) frame. ā, b̄ are the already filtered triples; pr is six scratch buffers, so each moment costs one batched apply of six products. The cross moment sums both orderings before filtering, which the linearity of the filter permits.

Two decompositions of the flux use this with different pairs: the Germano split feeds it the filtered velocity and its residual, and the Helmholtz split feeds it the rotational and divergent parts.

A point's three local components are written over slots 1-3 of the same tensor, safe because all six of that point's values are read first.

active_area(grid) -> T

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

analyze_sweep(u, v, w, grid, ws, plan) -> analysis or nothing

Forward-transform the raw inputs of a flux computation once, for reuse across every scale of a sweep.

A spectral filter is analyze → multiply by Ĝ(|k|, ℓ) → synthesize, and only the multiply depends on the scale, while every field a flux computation filters is raw: the velocities and their pairwise products. So a sweep over S scales needs 5 + 5S transforms rather than 10S.

Returns nothing for an engine with no shareable analysis — a real-space filter does all its work per scale — and the caller then runs the ordinary per-scale path. plan may be any one of the sweep's per-scale plans; they share a forward transform.

band_energies(u, v, w, grid, kernel, scales; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (; bands, resolved, total, band_maps, resolved_map)

Split the kinetic energy into contributions from each scale band, using the repeated-filter generalization of the Germano identity rather than by band-passing the velocity.

With scales in ASCENDING order ℓ₁ < ℓ₂ < … < ℓ_N, define the repeatedly filtered fields

f₀ = u ,   f_n = G_{ℓ_n} * f_{n-1} ,

so f_n has had every scale below ℓ_n removed, successively. Band n holds the energy the n-th application removed, which is the generalized second moment at that level:

k_n = ½ τ_{ℓ_n}(f_{n-1}; f_{n-1}) = ½[ (|f_{n-1}|²)‾_{ℓ_n} − |f_n|² ] ,

and the decomposition is exact:

½⟨|u|²⟩ = Σ_{n=1}^N ⟨k_n⟩ + ½⟨|f_N|²⟩ .

The sum telescopes because each ⟨G * x⟩ = ⟨x⟩ — i.e. because the filter conserves the domain mean. Measured, that holds to round-off on a periodic, unmasked grid (relative error 5e-16 for the top-hat, 2e-15 for the Gaussian) and the identity is exact there.

Anywhere the footprint is truncated, it is not, and the identity carries a residual of order ℓ/L: measured on the same field, 1.1e-2 relative on a BOUNDED grid (the footprint runs off the domain edge) and 1.1e-2 on a masked periodic grid under ZeroFill (energy is smeared onto masked cells, which report zero). Deformable renormalizes that leakage away and does better on a masked domain — 4.8e-4 — at the cost of the commutation property ZeroFill is the default for. So: read the bands as exact on a periodic unmasked domain, and as carrying an O(ℓ/L) boundary residual otherwise.

Why not band-pass the velocity

The obvious alternative, u = ū₀ + Σ(ū_n − ū_{n-1}), gives ½⟨|u|²⟩ = ½⟨|ū₀|²⟩ + ½Σ_{n,m}⟨ū_n · ū_m⟩ — cross terms of indefinite sign, so there is no well-defined energy at a given scale at all (Aluie & Eyink 2009). The second-moment form above has no cross terms by construction, and k_n ≥ 0 pointwise iff the kernel is non-negative — so use a positive kernel here (TopHatKernel, GaussianKernel, SmoothHatKernel, HyperGaussianKernel); a signed one such as Kernels.HighOrderKernel can give negative band energies.

maps = true additionally returns the per-band and resolved MAPS; they are N+1 full fields and most callers reduce them straight to the scalars below, so they are not built by default and band_maps is then nothing. filter_plans accepts a prebuilt sweep family.

Returns the per-band domain-averaged energies bands (length N), the energy left in f_N (resolved), their sum total, and the corresponding pointwise maps.

References

• Germano, M. (1992). J. Fluid Mech. 238, eq. (33).
• Aluie, H., & Eyink, G. L. (2009). Localness of energy cascade in hydrodynamic turbulence. Phys. Fluids 21, 115108, Appendix 2.

compressible_flux!(ws::Favre3DWorkspace, u, v, w, ρ, P, grid, kernel, scale; ...)

In-place true-3-D compressible_flux. Two batched applies carry the budget: eight fields for the density, the pressure and both velocity forms, then the six mass-weighted products.

compressible_flux!(ws, u, v, ρ, P, grid, kernel, scale; filter_plan=nothing, deriv_plan=nothing, ...)
    -> (; Π, Λ, pressure_dilatation, ρ̄, P̄, ũ, ṽ)

In-place compressible_flux. 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.

compressible_flux!(ws::SphericalFavreWorkspace, u, v, ρ, P, grid, kernel, scale; ...)

In-place spherical compressible_flux. Two batched applies carry the whole budget: eight fields for the density, the pressure and both velocity forms, then the six mass-weighted products.

compressible_flux(u, v, w, ρ, P, grid::StructuredGrid{T,G,3}, kernel, scale; ...)
    -> (; Π, Λ, pressure_dilatation, ρ̄, P̄, ũ, ṽ, w̃)

True three-dimensional variable-density (Favre) budget. The three terms are the ones the 2-D methods compute, over all six independent components of the 3×3 Favre stress and the full 3×3 resolved strain — including the genuine vertical derivatives the layer-stack path drops:

Π = −ρ̄ S̃_ij τ̃(u_i,u_j) ,   Λ = (1/ρ̄) ∂_j P̄ · τ̄(ρ,u_j) ,   P̄ ∇·ū

On a Cartesian volume the components are (x, y, z) as given. On a spherical shell every quantity carrying a direction goes through planetary-Cartesian coordinates and comes back to local (east, north, radial), and the strain carries the shell's curvature terms — the convention the true-3-D compute_Π! uses. ∇·ū is taken as the trace of the unweighted strain, so it picks those terms up with it.

compressible_flux(u, v, ρ, P, grid::AbstractGrid{<:SphericalGeometry}, kernel, scale; ...)
    -> (; Π, Λ, pressure_dilatation, ρ̄, P̄, ũ, ṽ)

Spherical form of the variable-density (Favre) budget, on any grid resolving two tangent directions.

The three terms are the ones the Cartesian method computes; each piece that carries a direction is built in planetary-Cartesian coordinates and rotated back to the local frame, since a local (east, north) pair filtered component-wise is not a filtered vector (Aluie 2019). That applies to the Favre velocity ũ, to the Favre stress τ̃, and to the unweighted subscale mass flux τ̄(ρ,u_j) that baropycnal work contracts against.

Two spatial operators pick up the local frame's curvature: the resolved Favre strain that Π contracts, and the unweighted divergence in the pressure-dilatation term,

∇·ū = ∂ū_e/∂x + ∂v̄_n/∂y − v̄_n tanφ/R ,

taken here as the trace of the unweighted strain, so it carries that term by construction.

compressible_flux(u, v, ρ, P, grid, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (; Π, Λ, pressure_dilatation, ρ̄, P̄, ũ, ṽ)

The variable-density (Favre) cross-scale energy budget of Aluie (2013). Returns the three terms of that budget which act on the large-scale kinetic energy ρ̄|ũ|²/2, plus the filtered fields they are built from.

Favre filtering

f̃ ≡ (ρf)‾/ρ̄ is the density-weighted filter. It exists because it is the one that makes the filtered continuity equation close exactly, ∂_tρ̄ + ∂_i(ρ̄ũ_i) = 0; the unweighted filter does not. It is linear but does not commute with derivatives, so the two filters are not interchangeable and the budget genuinely needs both — which is the source of the trap below.

The three terms

Π = −ρ̄ ∂_j ũ_i τ̃(u_i,u_j) ,   τ̃(u_i,u_j) = (ρu_iu_j)‾/ρ̄ − ũ_iũ_j        deformation work
Λ = (1/ρ̄) ∂_j P̄ · τ̄(ρ,u_j) ,  τ̄(ρ,u_j)   = (ρu_j)‾ − ρ̄ū_j              baropycnal work
                                                                        (τ̄ is UNWEIGHTED)
P̄ ∇·ū                                                                   pressure dilatation

Π and Λ both pit a large-scale field against small-scale fluctuations, so both transfer energy across scales. P̄∇·ū involves only large scales and cannot — it is a conversion between kinetic and internal energy at the resolved scale, not a cascade term.

The trap

Λ is frequently absorbed into the pressure term by writing it as P̄∇·ũ (plus a transport term) and then dismissed as "large-scale pressure dilatation that needs no modelling". That is wrong: the budget term is P̄∇·ū with the unweighted divergence, and writing ∇·ũ silently destroys Λ — a genuine cross-scale transfer. This implementation keeps them separate and uses ū for the dilatation; the suite asserts that Λ is non-zero for a baroclinic configuration, so it cannot be quietly dropped.

Asymptotics

For a smooth field, Lees & Aluie (2019) give Λ ≈ (C₂ℓ²/ρ̄)·c_d·[∇P̄·S̄·∇ρ̄ + ½ ω̄·(∇ρ̄ × ∇P̄)] with C₂ the kernel's second moment — a strain-generation part plus a baroclinic part that survives even in pure solenoidal flow. That C₂ ≠ 0 requirement is another reason the flux framework wants a kernel with a NON-vanishing second moment; see Kernels.HighOrderKernel for the kernels that deliberately give it up.

References

• Aluie, H. (2013). Scale decomposition in compressible turbulence. Physica D 247, 54–65.
• Lees, A., & Aluie, H. (2019). Baropycnal work: a mechanism for energy transfer across scales. Fluids 4, 92. doi:10.3390/fluids4020092

compute_Π!(Π::AbstractArray{T,3}, u, v, w, grid::StructuredGrid{T,Cartesian,3}, kernel, scale; mask_strategy=ZeroFill(), 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{T,Spherical,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{T,Spherical,3}, kernel, scale; workspace=nothing, backend=AutoBackend(), mask_strategy=ZeroFill())

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, deriv_plan=nothing, backend=AutoBackend(), mask_strategy=ZeroFill(), method=RealSpace())

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 (Operators.gradient_plan + Operators.gradient!). method defaults to RealSpace() here. 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_Π!(Π::AbstractVector, u, grid::StructuredGrid{T,Cartesian,1}, kernel, scale; workspace=nothing, backend=AutoBackend(), mask_strategy=ZeroFill())

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=ZeroFill())

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 (Operators.gradient! over a Operators.gradient_plan, both components from one neighbour sweep) and real-space filtering uses the scattered per-point footprint. Pass a prebuilt deriv_plan = FG.Operators.gradient_plan(grid) (and a reusable workspace) to avoid rebuilding them per call across a scale sweep.

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

In-place true-3D 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.

Each stress tensor comes from one batched apply of its six symmetric products. The cross stress sums the two orderings before filtering, which the linearity of the filter permits:

τX_ij = M(uʳ_i, uᵈ_j) + M(uᵈ_i, uʳ_j) = (uʳ_i uᵈ_j + uᵈ_i uʳ_j)‾ − (ūʳ_i ūᵈ_j + ūᵈ_i ūʳ_j)

so it costs one product per component, and the diagonal i = j picks up its factor of two from the same expression. The two strains read the filtered triples the stresses already needed, so a whole evaluation is four batched applies over 24 fields.

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!(ws::SphericalPiDecomposed3DWorkspace, u, v, w, u_rot, v_rot, w_rot, grid, kernel, scale; ...)

In-place spherical volumetric compute_Π_decomposed. Four batched applies carry the whole split: six planetary components of the two parts, then six products per stress.

compute_Π_decomposed!(ws::SphericalPiDecomposedWorkspace, u, v, u_rot, v_rot, grid, kernel, scale; ...)
    -> (; total, rotational, cross, divergent)

In-place spherical compute_Π_decomposed. One batched apply carries the six planetary components of both parts; three more carry the three stresses.

compute_Π_decomposed(u, v, w, u_rot, v_rot, w_rot, grid::StructuredGrid{T,<:SphericalGeometry,3}, kernel, scale; ...)
    -> (; total, rotational, cross, divergent)

Spherical volumetric form of the rotational/divergent split: the same both-sides decomposition the Cartesian volume takes, with the shell's two conventions applied. The three stresses are formed as generalized second moments in planetary-Cartesian coordinates and rotated back to local (east, north, radial), and each part's strain carries the shell's curvature terms — so the channels still sum to the flux the volumetric compute_Π! computes.

compute_Π_decomposed(u, v, w, u_rot, v_rot, w_rot, grid::StructuredGrid{T,Cartesian,3}, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (; 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::AbstractGrid{<:SphericalGeometry}, kernel, scale; ...)
    -> (; total, rotational, cross, divergent)

Spherical form of the rotational/divergent split, on any grid resolving two tangent directions — the configuration of Buzzicotti et al. (2023), where the Helmholtz parts come from scalar potentials filtered as scalars.

Both sides of the bilinear contraction split exactly as they do on a plane; what the sphere changes is how each side is built. The stress is bilinear, so its three pieces τʳʳ, τᵈᵈ and τ_X are formed as generalized second moments in planetary-Cartesian coordinates and rotated back to the local frame (Aluie 2019) — the same moments the spherical tau_decomposition! takes, with the rotational and divergent parts in place of the filtered velocity and its residual. The strain is linear, so it splits with no cross term, and each part's strain carries the local frame's tanφ/R curvature terms.

compute_Π_decomposed(u, v, u_rot, v_rot, grid, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (; 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. Contracting the split stress against the full, undecomposed strain S̄ — a one-sided split — is only correct when S̄ᵈ ≡ 0, and silently wrong whenever the divergent part carries strain of its own.

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

compute_Π_strain_convergence!(ws, u, v, grid, kernel, scale; filter_plan=nothing, deriv_plan=nothing, ...)
    -> (; total, strain, convergence, divergence, strain_magnitude)

In-place compute_Π_strain_convergence. 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_Π_strain_convergence!(ws::ΠWorkspace, u, v, grid::AbstractGrid{<:SphericalGeometry}, kernel, scale; ...)
    -> (; total, strain, convergence, divergence, strain_magnitude)

In-place spherical compute_Π_strain_convergence, over the same workspace compute_Π! uses. Returns views of ws's buffers, valid until the next call on it.

compute_Π_strain_convergence(u, v, grid::AbstractGrid{<:SphericalGeometry}, kernel, scale; ...)
    -> (; total, strain, convergence, divergence, strain_magnitude)

Spherical form of the strain/convergence split, on any grid resolving two tangent directions.

The split is a statement about the filtered strain tensor, so it carries over from the plane unchanged once that tensor is the spherical one:

δ̄ = S̄_ee + S̄_nn        σ̄_n = S̄_ee − S̄_nn        σ̄_s = 2·S̄_en

Substituting S̄_xx = (δ̄+σ̄_n)/2, S̄_yy = (δ̄−σ̄_n)/2, S̄_xy = σ̄_s/2 into −S̄:τ̄ collapses to Π_α − Π_δ term by term, on any metric. What the sphere changes is the two tensors fed in: S̄ picks up the tanφ/R curvature terms of the local frame, and τ̄ is formed in planetary-Cartesian coordinates and rotated back (Aluie 2019), exactly as compute_Π! does — this shares that code, so the two agree by construction and the suite asserts total == compute_Π!.

Takes a ΠWorkspace: the spherical path needs the three planetary components and their six moments, which that workspace already carries.

compute_Π_strain_convergence(u, v, grid, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (; total, strain, convergence, divergence, strain_magnitude)

Split the 2D cross-scale flux into the two production terms of Srinivasan, Barkan & McWilliams (2023), eq. (10), by diagonalizing the filtered strain tensor. With

δ̄   = ū_x + v̄_y            divergence         (rotation invariant)
σ̄_n = ū_x − v̄_y            normal strain
σ̄_s = ū_y + v̄_x            shear strain
ᾱ   = √(σ̄_n² + σ̄_s²)       strain magnitude   (rotation invariant)

the flux separates into

Π = Π_α − Π_δ ,   Π_α = (τ_vv − τ_uu) σ̄_n/2 − τ_uv σ̄_s ,   Π_δ = (τ_vv + τ_uu) δ̄/2 ,

with Π_α the deformation/shear production — energy transferred by straining, present even in non-divergent flow — and Π_δ the convergence production, which vanishes identically for a non-divergent field and is the term that paper adds. Setting δ̄ = 0 recovers Polzin (2010); the equivalent Π = E′(γᵖ ᾱ − δ̄) form with E′ = (τ_vv + τ_uu)/2 recovers Jing et al. (2017), and |γᵖ| ≤ 1 gives the bound |Π_α| ≤ ᾱ E′.

Π_α − Π_δ is algebraically identical to the direct Π = −S̄:τ̄ that compute_Π! computes — expanding eq. (10) collapses to −τ_uu ū_x − τ_uv(ū_y + v̄_x) − τ_vv v̄_y. The two are therefore a genuine cross-check rather than a tautology: they contract different combinations of the same four derivatives, so a sign or an ordering error in either shows up as a disagreement. The suite asserts they match to round-off on masked and unmasked grids.

Returns flux maps in W m⁻³, plus the two rotation invariants, which are the natural axes to bin the flux against (divergence = δ̄, strain_magnitude = ᾱ).

References

• Srinivasan, K., Barkan, R., & McWilliams, J. C. (2023). A forward energy flux at submesoscales driven by frontogenesis. J. Phys. Oceanogr. 53(1), 287–305. doi:10.1175/JPO-D-22-0001.1

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

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=ZeroFill())

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 Diagnostics.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!(out, ws, grid, has_w, total_area) -> out

Per-slice E(ℓ) from a batched workspace: out holds one energy per slice, shaped like the workspace's trailing axes.

E(ℓ) is a mean over the domain, so unlike everything else on the batch path this reduces the spatial axes away while leaving the batch axes intact — it cannot simply broadcast. The mask and cell areas are spatial-only, so each slice reduces against the same geometry.

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.

enstrophy_flux!(Z, ws, u, v, grid, kernel, scale; filter_plan=nothing, deriv_plan=nothing, ...) -> Z

In-place enstrophy_flux. With ws and both plans supplied, a repeated evaluation allocates nothing.

enstrophy_flux(u, v, grid, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill()) -> Z

Cross-scale enstrophy flux (Rivera, Aluie & Ecke 2014, eq. 16),

Z_ℓ = −∂_j ω̄_ℓ · τ_ℓ(u_j, ω) ,   τ_j = (u_j ω)‾ − ū_j ω̄ ,   ω = ∂v/∂x − ∂u/∂y ,

the enstrophy analogue of compute_Π!: positive means enstrophy moving to smaller scales. In 2-D turbulence this is the quantity with a forward cascade while Π cascades inverse, so the two are usually read together.

Gauge

This is the deformation (subtracted) form, the same gauge Π uses: the resolved product ū_j ω̄ is subtracted, which is what makes it pointwise Galilean-invariant. The unsubtracted alternative −∂_j ω̄ (u_j ω)‾ differs from it by a transport divergence, and while the two share a spatial mean on a homogeneous domain they "differ qualitatively as well as quantitatively" on an inhomogeneous or masked one (Aluie 2011; Aluie, Hecht & Vallis 2018). Mixing gauges between Π and Z would make the pair internally inconsistent, so only this one is provided.

Structurally Z is tracer_variance_flux with θ = ω, and that is how it is computed — the enstrophy is the "variance" of the vorticity. The separate entry point exists because the caller should not have to know to form ω with the matching derivative operator.

References

• Rivera, M. K., Aluie, H., & Ecke, R. E. (2014). The direct enstrophy cascade of two-dimensional soap film flows. Phys. Fluids 26, 055105. doi:10.1063/1.4873579

favre_filter!(out, tmp, f, ρ, ρ̄, plan) -> out

f̃ = (ρf)‾/ρ̄, given an already-filtered ρ̄ and a scratch array tmp. The building block of compressible_flux, exposed because a caller filtering their own tracer Favre-style should not have to reimplement it (and get the weighting backwards).

filtering_spectrum(u, v, w, grid, kernel, scales; L=1, backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (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. scales need not be uniform. Returns the filtering wavenumbers k_ℓ and the density Ẽ per scale.

The k_ℓ = C/ℓ convention, and why it must be stated

L is the region length, and k_ℓ = L/ℓ is the Sadek–Aluie convention: with their Fourier series f(x) = Σ_k f̂(k) e^{i(2π/L)k·x}, k is a dimensionless index, so L is the domain size. The default L = 1 instead gives k_ℓ = 1/ℓ, matching Storer et al. (2022, 2023) and FlowSieve. A third convention, k_ℓ = 2π/ℓ (Rivera, Aluie & Ecke 2014), is L = 2π.

The choice rescales the answer. Under k_ℓ = C/ℓ the density carries a Jacobian dℓ/dk_ℓ = -ℓ²/C, so Ẽ scales as 1/C while k_ℓ scales as C. Comparing amplitudes — or peak locations — against a Fourier spectrum or against another code is meaningless unless the conventions match. The cumulative energy cumulative_energy is convention-free; only the density is not.

Limits

• Slope ceiling. Sadek & Aluie eq. (18): a kernel with p vanishing moments recovers a true k^{-α} spectrum only for α < p + 2, and otherwise saturates at k^{-(p+2)}. Both TopHatKernel and GaussianKernel have p = 1, so the measured slope locks at k⁻³. This bites hardest in 2-D and QG work, where the enstrophy-range target slope is ≈ k⁻³. The flux Π is unaffected — this is a limitation of the spectrum diagnostic alone.
• Kernel admissibility. Ẽ(k_ℓ) ≥ 0 is guaranteed only when d|Ĝ(k)|²/dk ≤ 0. By default this function throws for a kernel that fails it; pass policy = ForceSpectrum() to compute it anyway. See AbstractSpectrumPolicy and Kernels.transfer_monotone.

References

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

gate_spectrum(kernel, policy) -> Bool

Apply policy to kernel, returning whether a spectral density should be computed. Throws under StrictSpectrum for a non-monotone kernel; warns once under ForceSpectrum.

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::Sym3TauWorkspace, u, v, w, grid::StructuredGrid{T,G,3}, kernel, scale; ...)

In-place true-3-D tau_decomposition. Four batched applies carry the whole split: the three velocity components, the six second-filtered fields, then six products per tensor.

tau_decomposition!(ws::Sym3TauWorkspace, u, v, grid, kernel, scale; filter_plan=nothing, ...)

In-place spherical tau_decomposition.

Like compute_Π!'s spherical branch, the Leonard/Cross/Reynolds moments are formed in PLANETARY-CARTESIAN coordinates — where filtering commutes with the moment and residual operations (Aluie 2019) — and each resulting symmetric 3x3 tensor is then rotated back to the local (east, north) frame. L + C + R = τ still holds exactly, the rotation being linear.

Every filter here goes through a batched apply: three velocity components, then six second-filtered fields, then six products per tensor. The cross term uses the identity

C_ij = M(ū_i, u'_j) + M(u'_i, ū_j) = (ū_i u'_j + u'_i ū_j)‾ − (ū̄_i ū'_j + ū'_i ū̄_j)

so it costs one product per component rather than two — the filter being linear, the two orderings can be summed before filtering instead of after.

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, w, grid::StructuredGrid{T,G,3}, kernel, scale; ...) -> (; L, C, R)

True three-dimensional Germano split: the same generalized central moments as the 2-D methods, over all six independent components of the 3×3 subfilter stress. Each tensor comes back as (; xx, xy, xz, yy, yz, zz), and L + C + R = τ holds componentwise.

On a Cartesian metric the components are (x, y, z) as given. On a spherical volumetric shell the moments are formed in planetary-Cartesian coordinates and rotated back, so the components are local (east, north, radial) — the convention the true-3-D compute_Π! uses.

tau_decomposition(u, v, grid, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> (; 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::AbstractGrid{<:SphericalGeometry}, kernel, scale; ...) -> (; L, C, R)

Spherical counterpart of the Cartesian method above, on any grid resolving two tangent directions — structured, curvilinear, a scattered node set or a sphere pixelization: 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::TracerFlux3DWorkspace, u, v, w, θ, grid, kernel, scale; ...) -> Πθ

In-place true-3-D tracer_variance_flux, on either metric. One batched apply carries the whole flux: the three velocity components, the tracer, and the three velocity–tracer products.

tracer_variance_flux!(Πθ, ws::SphericalTracerFluxWorkspace, u, v, θ, grid, kernel, scale; ...) -> Πθ

In-place spherical tracer_variance_flux. One batched apply carries the whole flux: the three planetary velocity components, the tracer, and the three velocity–tracer products.

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{T,Cartesian,3}, kernel, scale; backend=AutoBackend(), mask_strategy=ZeroFill())
    -> Πθ

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{T,<:SphericalGeometry,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=ZeroFill())
    -> Πθ

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::AbstractGrid{<:SphericalGeometry}, kernel, scale; ...) -> Πθ

Spherical form of the tracer-variance flux, on any grid resolving two tangent directions — structured, curvilinear, a scattered node set or a sphere pixelization. τ_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.

vorticity!(ω, u, v, grid[, deriv_plan]; scratch = nothing) -> ω

Vertical (radial) component of the relative vorticity on a grid resolving two tangent directions — structured, curvilinear, a scattered node set or a sphere pixelization. Masked cells are zeroed, as everywhere else. Uses the same gradient operator the flux diagnostics do, so ω and the gradients it is later contracted against are consistent to the last bit.

On a Cartesian metric this is ω = ∂v/∂x − ∂u/∂y. On a sphere the curl in orthogonal curvilinear coordinates carries the metric's own curvature term,

ζ = (1/(R cosφ))[∂v/∂λ − ∂(u cosφ)/∂φ] = ∂v/∂x − ∂u/∂y + u tanφ/R ,

with ∂/∂x, ∂/∂y the distance derivatives the gradient operator returns. The spherical strain compute_Π! contracts carries the same tanφ/R factor, so ω and Π share one gauge.

The two derivatives cannot be accumulated into one array, so a second full-size buffer is needed. Pass scratch to supply it: a scale sweep calls this once per scale, and the in-place form is meant to allocate nothing.

vorticity(u, v, grid[, deriv_plan]) -> ω

Allocating vorticity!.