Kernels and derivatives

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: σ² = ℓ²/(2α) = ℓ²/12, which is the second moment of a top-hat of width ℓ in one dimension.
• α = 4 reproduces FlowSieve's default Gaussian (which also treats ℓ as a diameter), so GaussianKernel(; α = 4) is directly comparable to FlowSieve output.

Variance matching is dimension-dependent

Because the Gaussian is separable, its per-component variance is ℓ²/(2α) in every dimension. The top-hat here is not a separable box but the disk/ball of radius ℓ/2, whose per-component variance therefore shrinks with dimension. Measured (quadrature, units of ℓ²):

So α = 6 matches the top-hat exactly on a 1-D grid, and on a 2-D grid it is 15% wider in RMS than the top-hat of the same nominal ℓ. Use GaussianKernel(; α = 8) when the point is a like-for-like comparison against TopHatKernel on a 2-D grid, and α = 10 in 3-D.

Discretization behaves very differently between the two, also measured: the Gaussian's footprint reproduces its continuum ⟨x²⟩ to 2e-10 already at ℓ = 8Δx, while the top-hat's staircased disk boundary leaves it 2.0% low at ℓ = 8Δx, falling as O(Δx) to 0.24% at ℓ = 64Δx.

HighOrderKernel{P}(; b_over_ℓ = 1/8) <: AbstractFilterKernel

Piecewise-constant kernel with P vanishing moments (Sadek & Aluie 2018 §V): a positive body of half-width ℓ/2 flanked by alternating-sign limbs of width b = b_over_ℓ · ℓ, whose amplitudes are solved so that ∫xⁿG dx = 0 for n = 1 … P.

• P = 3 is the paper's M^I: one negative limb, support |x| < ℓ/2 + b.
• P = 5 is M^II: a negative limb then a positive one, support |x| < ℓ/2 + 2b.

At the paper's b = ℓ/8 the amplitudes reduce to (1, -64/61) and (1, -568/257, +200/257) respectively, relative to the body (the weights here are unnormalized, as everywhere in this module).

This is a SEPARABLE kernel, not a radial one

G(x, y) ≡ G(x) G(y) — the paper's own generalization, and a genuinely different object from a radially symmetric kernel: its support is a square, not a disk. kernel_weight therefore has no method for it and throws, and it can only be applied on a grid with separable axes (rectilinear StructuredGrid). A curvilinear or unstructured grid has no axes to factor over, so plan_filter refuses it there rather than silently applying a radial approximation.

Slope recovery, and the domain it needs

Measured on a synthetic k^-α field, N = 256, fit over ℓ = 8-32Δx: P = 3 recovers k⁻⁴ (-4.00 against a true -4) where TopHatKernel/GaussianKernel saturate at -2.83/-2.73. P = 5 does not reach k⁻⁷ in that box.

That last one is a property of the domain, not of the kernel. With 𝓔(ℓ) = ∫E(k)|Ĝ(kℓ)|²dk and x = kℓ,

Ẽ(k_ℓ) ∝ k_ℓ^{-α} · I ,     I = -∫₀^∞ x^{1-α} (|Ĝ|²)'(x) dx

so the slope is -α unless I diverges; since 1 - |Ĝ|² ∝ x^{p+1} the integrand is x^{p+1-α} and I diverges at small x = kℓ iff α ≥ p+2 (Sadek & Aluie eq. 18). The saturation is therefore an INFRARED effect, set by scales larger than the filter: seeing it needs many decades of spectrum below 1/ℓ, and the integrand sharpens with p. A 256² box leaves ~1.5 decades there — enough for p = 1, marginal for p = 3, not enough for p = 5. A finer grid does not help; it adds modes above 1/ℓ, which I does not see.

Two cautions

• It needs resolution. Every limb must span at least one cell, b ≥ Δx, i.e. ℓ ≥ 8Δx at the default b_over_ℓ = 1/8; plan_filter warns below that. Weights are the kernel's integral over each cell, not its value at the node (profile_cell_average): the discrete m₂ residual then falls as O(Δx²) rather than point sampling's O(Δx), 100× smaller at ℓ = 128Δx. Still nonzero at finite resolution, so the effective order is always below the nominal one.
• Do not use it for Π. The vanishing moments are bought with negative weights, so τ is not realizable, and |Ĝ|² reaches 0.282 (P = 3) and 1.0125 (P = 5) past its first zero — the latter amplifies. For a flux use TopHatKernel or GaussianKernel.

HyperGaussianKernel(; α = 1.0) <: AbstractFilterKernel

Super-Gaussian, G_ℓ(d) ∝ exp(-α (2d/ℓ)⁴). Flatter in the core and steeper in the skirt than a Gaussian of the same nominal width — closer to a box while staying smooth and strictly positive.

Real space only, for the same reason as SmoothHatKernel: exp(-r⁴) has no elementary Fourier transform. Its |Ĝ|² is not monotone (max sidelobe 0.056 measured), so it cannot produce a filtering spectrum either.

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.

SmoothHatKernel(; steepness = 10.0) <: AbstractFilterKernel

Tapered top-hat, G_ℓ(d) ∝ ½(1 − tanh(s(2d/ℓ − 1))), the kernel used in the Storer et al. papers. It is a top-hat whose rim is smoothed over a width ≈ ℓ/(2s), so it keeps the box's compactness while removing the discontinuity that makes the box's transfer function ring.

steepness = 10 is the published value, i.e. ½(1 − tanh((D − 1)/0.1)) with D = 2d/ℓ; larger values approach TopHatKernel, smaller ones approach a broad bell.

Real space only — there is no closed-form transfer function, so method = Spectral() throws rather than silently substituting a different kernel. Its |Ĝ|² is not monotone (max sidelobe 0.119 measured), so transfer_monotone is false and it cannot produce a filtering spectrum.

TopHatKernel <: AbstractFilterKernel

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

check_spectrum_kernel(kernel::AbstractFilterKernel)

Throw an ArgumentError unless kernel satisfies transfer_monotone. This is the strict reading of the admissibility condition, selected by Diagnostics.StrictSpectrum(); the cumulative energy needs no such check.

is_radial(kernel) -> Bool

Whether the kernel is a function of the scalar distance alone, i.e. whether kernel_weight has a method for it. HighOrderKernel is the one kernel here that is separable but NOT radial.

is_separable(kernel) -> Bool

Whether the kernel factors as G(x₁,…,x_N) = ∏ G(x_d), so it can be applied as successive 1-D passes and MUST be, if it has no radial form. See kernel_profile.

kernel_profile(kernel, δ::T, ℓ::T) -> T

One-dimensional profile of a SEPARABLE kernel at signed-or-unsigned axis offset δ. The full weight at a displacement (δ₁, …, δ_N) is the product ∏ kernel_profile(kernel, δ_d, ℓ).

Defined for HighOrderKernel, which is separable by construction, and for GaussianKernel, which is separable as an identity (exp(-α r²/ℓ²) factors), where it coincides with kernel_weight. A kernel with no method here is radial-only.

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 ℓ.

limb_amplitudes(k::HighOrderKernel{P}, ::Type{T}) -> NTuple{P == 3 ? 1 : 2, T}

Limb amplitudes relative to the body, from the vanishing-moment conditions. Dimensionless in β = b/ℓ, so ℓ does not enter.

profile_cell_average(kernel, δ::T, Δ::T, ℓ::T) -> T

The 1-D profile averaged over a cell of width Δ centred at offset δ, i.e. (1/Δ) ∫_{δ-Δ/2}^{δ+Δ/2} G₁(s) ds.

The default is the point sample kernel_profile(kernel, δ, ℓ) — for a smooth kernel the two agree to O(Δ²) and the point sample is what every existing engine and stored reference uses, so it is left exactly as it was. HighOrderKernel overrides it with the exact integral, which is what makes it well-posed on a non-uniform grid; see profile_integral.

profile_integral(kernel, t::T, ℓ::T) -> T

∫₀ᵗ G₁(s) ds for the 1-D profile, exact. Odd in t, since the profile is even.

This exists because HighOrderKernel is discontinuous with limbs only b = ℓ/8 wide. Point-sampling such a kernel onto a grid is a poor quadrature — the sampled mass of a limb depends on exactly where the nodes fall, which on a non-uniform axis can miss a limb entirely and drive the normalization denominator negative (measured: the order-5 profile reaches -1.77e3 on a ±3% jittered axis, and a filtered constant comes back sign-flipped). Integrating the kernel over each cell instead of sampling it at the node removes that failure completely and is exact for a piecewise-constant kernel; Sadek & Aluie 2018 §III.F point at the same remedy.

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.

transfer_monotone(kernel::AbstractFilterKernel) -> Bool

Whether d|Ĝ(k)|²/dk ≤ 0 holds on (0, ∞).

Sadek & Aluie (2018) eq. (21): this is the condition under which the filtering spectral density Ẽ(k_ℓ) is guaranteed non-negative. A kernel that fails it can produce a spectrum with negative values, which is why Diagnostics.StrictSpectrum — the default policy — refuses one.

The condition is sufficient, not necessary — Sadek & Aluie's own fallback argument concedes it is "not a rigorous proof" — so false means "not guaranteed", not "certainly negative". It is also narrower than it looks: it says nothing about the k^{-(p+2)} slope ceiling, which binds the Gaussian just as hard as the top-hat (see Diagnostics.filtering_spectrum).

There is deliberately no fallback method. A new kernel gets a MethodError here rather than a guessed answer, so its author has to establish which way it goes.

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 Operators.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).

gradient_plan(grid) -> Operators.GradientPlan

The least-squares tangent-plane gradient plan for grid, refusing a grid that carries no adjacency for it to fit.