Spatial filters

Oceanostics provides spatial filters that operate directly on Oceananigans fields. All filters are built on top of Oceananigans' KernelFunctionOperation, so they compose with the rest of the Oceananigans ecosystem (outputs, reductions, other operations, etc.) and run on both CPU and GPU.

Box filter

The BoxFilter computes a local running-mean (box average) of a field over one or more grid directions. The stencil size is controlled by the N keyword: each filtered direction uses an N-point symmetric average centred on the current cell — i.e. N cells contribute to one filtered output value. N must be an odd integer ≥ 3, and refers to the size of the filter stencil, not the grid.

Multi-direction filters are fused into a single kernel at compile time, so a 3D box filter performs one pass over the data, not three.

Basic usage

Build a reusable filter once with the keyword-only constructor, then apply it to any field by calling it — bf(ψ) returns a KernelFunctionOperation that you wrap in a Field:

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(32, 32), x=(0, 1), z=(0, 1),
                              topology=(Periodic, Flat, Bounded));

julia> c = CenterField(grid);

julia> set!(c, (x, z) -> sin(2π * x) * z);

julia> bf = BoxFilter(; dims=(1, 3), N=5)
BoxFilter(dims=(1, 3), N=5, boundary=:shrink)

julia> c̄ = Field(bf(c));

julia> size(c̄)
(32, 1, 32)

The same bf can be applied to as many fields as you like. A one-step shortcut BoxFilter(ψ; dims, N, boundary), which builds and applies the filter in a single call, is also accepted throughout.

Boundary handling

For Periodic directions, stencil offsets are always wrapped — the boundary keyword is silently ignored. For Bounded directions the boundary keyword selects how out-of-bounds offsets are treated:

`boundary`	Behaviour
`:shrink` (default)	Drop out-of-bounds offsets from sum and count (honest local average; stencil shrinks near walls).
`:edge`	Replicate the nearest interior cell (`ψ[1]` or `ψ[N]`).
`(left=a, right=b)`	Pad with constant `a` on the low side and `b` on the high side.

A single spec applies to every filtered dimension, or a tuple gives per-dimension control:

julia> bf_edge  = BoxFilter(; dims=(1, 3), N=3, boundary=:edge);

julia> bf_mixed = BoxFilter(; dims=(1, 3), N=3, boundary=(:shrink, (left=0.0, right=0.0)));

julia> size(Field(bf_edge(c))) == size(Field(bf_mixed(c))) == (32, 1, 32)
true

API reference

Oceanostics.SpatialFilters.BoxFilter — Type

BoxFilter(ψ; dims, N, boundary)

One-step form: apply a box filter to ψ directly, returning the KernelFunctionOperation that computes its local box-average. Equivalent to BoxFilter(; dims, N, boundary)(ψ).

Refer to BoxFilter(; dims, N, boundary) for the full description of the keyword arguments and boundary handling.

source

BoxFilter(; dims, N, boundary)

Build a reusable, field-less box filter that computes a local box-average over the directions listed in dims. The returned object is callable: applying it to a field, bf(ψ), returns a KernelFunctionOperation. Build the filter once and reuse it across many fields, or pass it to other diagnostics that accept a filter.

dims is a tuple of distinct integers drawn from (1, 2, 3) (where 1, 2, 3 correspond to x, y, z). A single integer (e.g. dims=2) is accepted as shorthand for filtering along that one direction.

N is the total number of grid points used by the filter stencil along each filtered direction — i.e. how many cells are averaged together to produce one filtered output value (e.g. N=3 is a 3-point running mean, N=5 is a 5-point running mean). N must be an odd integer ≥ 3 so the stencil is symmetric around the current cell. (This is the size of the filter stencil — not the size of the grid.)

A multi-directional filter is assembled as a single KernelFunctionOperation whose kernel function is a 1D BoxFilterKernel{d₁}, with the next dimension's BoxFilterKernel{d₂} (and so on) threaded into the argument list. The box average is separable, so when the operation is the operand of a Field (the standard Field(bf(ψ)) / compute! path) it is evaluated as d sequential 1D passes through intermediate fields — d × N reads per output cell instead of Nᵈ. If the filtered field is composed into another AbstractOperation (e.g. 2 * bf(c)) the original fused single-kernel evaluation runs instead.

See Performance notes in the documentation for what that costs and how to avoid it.

Boundary handling

Stencil offsets that leave the interior 1:Nd_grid of a direction (where Nd_grid is the number of cells along that direction) are handled per-direction. For Periodic directions offsets are always wrapped periodically, independent of the boundary keyword. For Bounded directions the boundary keyword picks the policy (default: :shrink):

:shrink — drop out-of-bounds offsets from both the sum and the count, so the filter is an honest local average whose effective stencil shrinks near a wall. This is the default for Bounded directions.
:edge — replicate the boundary-cell value (reads ψ[1] or ψ[Nd_grid] for offsets past either end).
(left=a, right=b) — pad with constant a on the low-index side and b on the high-index side (a and b are promoted to a common type).

boundary may be a single spec applied to every filtered dim, or a tuple with one entry per dim in dims (in the order the user passed them):

BoxFilter(; dims=(1, 2), N=7, boundary=:edge)
BoxFilter(; dims=(1, 2), N=7, boundary=(:shrink, :edge))
BoxFilter(; dims=(1,),   N=7, boundary=(left=0.0, right=1.0))

Because every policy wraps, clamps, or skips indices up front, halo_size(grid) does not constrain N: a small halo on a bounded direction is fine. The output location matches the location of the field the filter is applied to, which can be any input that supports the standard Oceananigans ψ[i, j, k] indexing contract (e.g. a Field or any AbstractOperation).

For Periodic directions the stencil must span at most one period: N ≤ 2*Nd_grid + 1, where Nd_grid is the number of cells along that direction. This is enforced when the filter is applied.

Examples

Build a box filter and apply it on a 2D (xz) grid that is periodic in x and bounded in z. The boundary=:edge keyword explicitly overrides the default :shrink policy for the bounded z-direction; x is Periodic so the boundary spec is silently ignored for that direction regardless.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(8, 8), x=(0, 1), z=(0, 1),
                              topology=(Periodic, Flat, Bounded));

julia> c = CenterField(grid);

julia> bf = BoxFilter(; dims=(1, 3), N=5, boundary=:edge)
BoxFilter(dims=(1, 3), N=5, boundary=:edge)

julia> bf(c) isa KernelFunctionOperation
true

A one-step shortcut BoxFilter(ψ; dims, N, boundary) is also accepted, which applies the filter to ψ immediately (equivalent to BoxFilter(; dims, N, boundary)(ψ)).

source

Gaussian filter

The GaussianFilter computes a Gaussian-weighted local average of a field over one or more grid directions. Along a single filtered direction the filtered value at cell $i$ is

\[\bar{\psi}_i = \frac{\sum_m \Delta_m \, e^{-(x_m - x_i)^2 / (2\sigma^2)} \, \psi_m} {\sum_m \Delta_m \, e^{-(x_m - x_i)^2 / (2\sigma^2)}} ,\]

where $m$ runs over the cells of the stencil centred on $i$; $x_m$ and $\Delta_m$ are the coordinate and width of cell $m$ along the filtered direction, $x_i$ is the current cell's coordinate, and $\sigma$ is the kernel's standard deviation in physical units (the same as the grid spacing). In words: each cell's Gaussian weight $e^{-(x_m - x_i)^2 / (2\sigma^2)}$ is multiplied by that cell's own width $\Delta_m$, and the same $\Delta_m$-weighted sum appears in the normalizing denominator. The normalization returns a constant field unchanged and keeps every boundary policy consistent, while the $\Delta_m$ factor is the quadrature weight (the $dx'$ of $\int G_\sigma(x-x')\,\psi(x')\,dx' \big/ \int G_\sigma(x-x')\,dx'$) that makes the discrete sum approximate the continuous Gaussian convolution. Multi-direction filters apply this to each filtered direction in turn.

σ is the only required parameter beyond dims — N is inferred per-direction from σ and the minimum cell spacing in that direction so the stencil extends at least 2σ on each side of the current cell. To override, pass N explicitly: a single odd integer (≥ 3) applies to every filtered dim, or a tuple of odd integers sets one count per dim.

dims and boundary work identically to BoxFilter.

For a worked end-to-end example — filtering a turbulent flow, a filter-width sweep, and a subfilter tracer flux — see the Spatial filtering example.

Variably spaced (stretched) grids

GaussianFilter works on stretched grids, but it is around 4 times slower than on regular grids:

Along a uniformly spaced direction every $\Delta_m$ is equal, so it factors out of the sum and cancels between numerator and denominator. The weights then reduce to a plain $e^{-(x_m - x_i)^2 / (2\sigma^2)}$ that is identical for every cell, so they are precomputed once and looked up — the fast path used on regular grids.
Along a variably spaced direction $x_m$ and $\Delta_m$ differ from cell to cell, so the weights cannot be precomputed and are evaluated on the fly. Keeping the per-cell $\Delta_m$ factor is what stops the average from being biased toward finely resolved regions; it makes the filter preserve constants exactly and linear fields to quadrature accuracy.

Directions are handled independently, so a grid that is uniform in x and y but stretched in z uses the fast path for the horizontal directions and the node-distance path for the vertical:

julia> stretched_grid = RectilinearGrid(size=(16, 16),
                                        x=(0, 1),
                                        z=k -> -1 + (k-1)/16 + 0.05sin(2π*(k-1)/16),
                                        topology=(Periodic, Flat, Bounded));

julia> cz = CenterField(stretched_grid); set!(cz, (x, z) -> sin(2π*x) * z);

julia> gf = GaussianFilter(; dims=(1, 3), σ=0.1)
GaussianFilter(dims=(1, 3), σ=0.1, N=nothing, boundary=:shrink)

julia> c̄z = Field(gf(cz));

julia> c̄z isa Field
true

Because the filter stores only its parameters (not a grid or field), the same gf applies to a field on any grid — here the uniform grid from the box-filter example above:

julia> c̄_gauss = Field(gf(c));

julia> c̄_gauss isa Field
true

Pass N to override the default stencil:

julia> gf_wide   = GaussianFilter(; dims=(1,), σ=0.05, N=11);

julia> gf_perdim = GaussianFilter(; dims=(1, 3), σ=0.05, N=(7, 11));

julia> (Field(gf_wide(c)) isa Field, Field(gf_perdim(c)) isa Field)
(true, true)

Performance notes

The staged (separable) evaluation only fires when the filter operation is the direct operand of a Field. Composing it into another AbstractOperation hides it from that dispatch, and the filter silently falls back to the fused, single-kernel path:

Field(gf(ψ))         # staged: d one-dimensional passes
Field(gf(ψ) - φ)     # fused:  one Nᵈ-point kernel

Both forms produce the same result, so the only symptom is speed. The penalty has two parts. First, the read count per output cell grows from d × N to Nᵈ, and, more importantly, the filtered operand is re-evaluated at every one of those Nᵈ stencil points instead of once per cell.

Second, filtering a stored Field is comparatively cheap either way, because reading an array element is cheap. Filtering an expensive KernelFunctionOperation is not, because the operation is recomputed from scratch at every stencil point.

A quick calculation of the sub-filter scale dissipation filter(ε) - ε̄ illustrates this, since KineticEnergyDissipationRate is a nine-term viscous-flux contraction rather than a stored array. We build a small LES and filter it horizontally with a width of ℓ = 8Δx, which is a 15-point stencil along each filtered direction:

using Oceananigans, Oceanostics

grid  = RectilinearGrid(size=(32, 32, 32), extent = (1, 1, 1))
model = NonhydrostaticModel(grid; closure = SmagorinskyLilly())

Δx = minimum_xspacing(grid)
gf = GaussianFilter(; dims=(1, 2), σ = 8Δx / (2 * sqrt(2 * log(2))))

ε  = KineticEnergyDissipationRate(model)                   # an expensive KernelFunctionOperation
εˡ = FilteredKineticEnergyDissipationRate(model, gf)  # dissipation of the filtered flow

Oceanostics.SpatialFilters is the staged path, Oceananigans.AbstractOperations the generic fused one. Timing both, after a warm-up call so that compilation is not counted:

julia> ∫εˢ_fused  = ∫dV(      gf(ε)  - εˡ);
julia> ∫εˢ_staged = ∫dV(Field(gf(ε)) - εˡ);
julia> compute!(∫εˢ_fused); compute!(∫εˢ_staged);   # warm up so compilation is not timed
julia> t_fused  = @elapsed compute!(∫εˢ_fused)11.837868437
julia> t_staged = @elapsed compute!(∫εˢ_staged)1.042157911

Absolute timings depend on the machine, and the ratio grows with the stencil width; the point is that materializing the filtered field first speeds up the calculation at the cost of one array.

The same reasoning applies to any filtered quantity assembled from operations rather than stored fields. Wrap the filtered velocities in Fields before forming ½(ū² + v̄² + w̄²), for instance. This is also why subfilter_stress_tensor materializes its filtered velocities and momentum fluxes internally.

To check whether a given expression is affected, pass it to check_filter_staging. It walks the operation tree and warns (returning false) if any multi-direction filter is positioned to run fused, so the fused path does not go unnoticed:

julia> check_filter_staging(Field(gf(ε)) - εˡ)         # staged: filtered field materialized firsttrue
julia> check_filter_staging(gf(ε) - εˡ; warn=false)    # fused: filter nested in the subtractionfalse

API reference

Oceanostics.SpatialFilters.GaussianFilter — Type

GaussianFilter(ψ; dims, σ, N, boundary)

One-step form: apply a Gaussian filter to ψ directly, returning the KernelFunctionOperation that computes its Gaussian-weighted local average. Equivalent to GaussianFilter(; dims, σ, N, boundary)(ψ).

Refer to GaussianFilter(; dims, σ, N, boundary) for the full description of the keyword arguments, the Gaussian weighting, stretched-grid handling, and boundary handling.

source

GaussianFilter(; dims, σ, N, boundary)

Build a reusable, field-less Gaussian filter that computes a Gaussian-weighted local average over the directions listed in dims. The returned object is callable: applying it to a field, gf(ψ), returns a KernelFunctionOperation. Build the filter once and reuse it across many fields, or pass it to other diagnostics that accept a filter.

σ is the standard deviation of the Gaussian kernel in physical units (the same units as the grid spacing). The filter approximates the continuous Gaussian convolution ∫ G_σ(x-x') ψ(x') dx' / ∫ G_σ(x-x') dx': a stencil cell of width Δₘ whose centre sits a physical distance r from the current cell centre contributes weight Δₘ · exp(-r² / 2σ²), and the result is normalized by the sum of the surviving weights, so all boundary policies behave consistently. On a uniform direction the Δₘ factor is constant and cancels, recovering the plain exp(-r² / 2σ²) weighting.

GaussianFilter supports both uniformly and variably spaced (stretched) directions, choosing the implementation per direction so the regular-grid case keeps its original speed:

Uniform direction — the weights are identical for every cell, so they are precomputed once in cell-offset units (σ_cells = σ / Δ, weight exp(-Δi² / 2σ_cells²) at cell offset Δi) and looked up in the hot loop. This is the fast path used on regular grids.
Stretched direction — the physical footprint of a fixed cell offset varies from cell to cell, so the weights cannot be precomputed; each is evaluated on the fly from the node coordinates and widths (Δₘ · exp(-(xₘ - xᵢ)² / 2σ²)). The cell-width factor Δₘ is the quadrature weight of the continuous convolution; it stops the average from being biased toward finely resolved regions and keeps constants (and, to quadrature accuracy, linear fields) preserved. The two paths agree exactly where they overlap (a uniform direction has xₘ - xᵢ = Δm·Δ and constant Δₘ). Directions are decided independently, so a grid that is uniform in x but stretched in z uses the fast path for x and the node-distance path for z.

N is the total number of grid points used by the filter stencil along each filtered direction — i.e. how many cells contribute to a single filtered output value. N must be an odd integer ≥ 3 so the stencil is symmetric around the current cell. (This is the size of the filter stencil — not the size of the grid.) If unspecified, N is inferred per-direction from σ and the smallest spacing along that direction so that the stencil extends at least 2σ to each side of the current cell everywhere (on a stretched direction it therefore reaches farther than 2σ where cells are large, which is harmless since the Gaussian weights there are tiny). To override, pass either a single odd integer (applied to every filtered dim) or a tuple with one odd-integer count per dim in dims (in the order the user passed them). For Periodic directions the stencil must span at most one period (N ≤ 2*Nd_grid + 1, where Nd_grid is the number of cells along that direction); this is enforced when the filter is applied.

See BoxFilter for the dims and boundary keyword documentation.

Performance notes

A multi-direction Gaussian filter is mathematically separable. Applying the filter returns a single composable KernelFunctionOperation, but when that operation is the operand of a Field (the standard Field(gf(ψ)) / compute! path), the implementation evaluates the filter as a sequence of 1D passes through intermediate fields. This reduces the per-output read count from N^d to d × N, which is the main reason multi-direction filters with wide stencils are competitive on GPUs. Mixed-spacing filters stage the same way — each direction's 1D pass uses its own (uniform or stretched) kernel. If the filtered field is composed into another AbstractOperation (e.g. 2 * gf(c)) it falls back to the fused, single-kernel evaluation.

See Performance notes in the documentation for what that costs and how to avoid it.

Examples

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(8, 8), x=(0, 1), z=(0, 1),
                              topology=(Periodic, Flat, Bounded));

julia> c = CenterField(grid);

julia> gf = GaussianFilter(; dims=(1, 3), σ=0.1)
GaussianFilter(dims=(1, 3), σ=0.1, N=nothing, boundary=:shrink)

julia> gf(c) isa KernelFunctionOperation
true

A one-step shortcut GaussianFilter(ψ; dims, σ, N, boundary) is also accepted, which applies the filter to ψ immediately (equivalent to GaussianFilter(; dims, σ, N, boundary)(ψ)).

source

Oceanostics.SpatialFilters.check_filter_staging — Function

check_filter_staging(op; warn=true)

Inspect the operation tree op (an AbstractOperation, a Field, or a reduction such as Integral(...)) for multi-direction Oceanostics filters — a BoxFilter or GaussianFilter over two or three directions — that will evaluate on the slow fused single-kernel path instead of the fast staged (separable) path.

A multi-direction filter is evaluated by its staged kernel only when it is the direct operand of a Field, e.g. Field(filter(ψ)). Composing the filter into any other operation — `Field(filter(ψ)

φ),2 * filter(ψ),Integral(filter(ψ))`, and so on — hides it from that dispatch and it silently

falls back to the fused path (Nᵈ reads per output cell instead of d × N, with the filtered operand recomputed at every stencil point). Both paths return the same values; only the speed differs. See Performance notes.

The fix is to materialize the filtered field first: wrap filter(ψ) in its own Field(...) before composing it, e.g. write Field(filter(ψ)) - φ rather than filter(ψ) - φ.

op is analyzed as if it were about to be wrapped in a Field and computed. Returns true when every multi-direction filter in op is positioned to run staged (nothing to fix), and false when at least one will run fused. When warn = true (the default), a false result also emits a @warn describing the problem. One-direction filters always use the single unrolled kernel (staged and fused coincide) and are never flagged.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1));

julia> c = CenterField(grid); φ = CenterField(grid);

julia> gf = GaussianFilter(; dims=(1, 2), σ=0.1);

julia> check_filter_staging(Field(gf(c)))          # staged: direct Field operand
true

julia> check_filter_staging(Field(gf(c)) - φ)      # staged: filtered field materialized first
true

julia> check_filter_staging(gf(c) - φ; warn=false) # fused: filter nested in a larger operation
false

source