Stencils

AbstractStencil

A neighbourhood shape in index space. Materialize it with offsets.

Concrete shapes: Axial, VonNeumann, Moore (alias Vertex), Diagonal, Anisotropic, Custom.

A shape supplies one method, offsets(s, ::Val{N}), returning its offsets as a tuple; foreach_offset, fold_offsets, nstencil and reach all follow from it, and their loops still unroll at compile time. Put whatever the offsets depend on in the type, so the tuple is inferable:

struct Upwind{R} <: AbstractStencil end
Stencils.offsets(::Upwind{R}, ::Val{N}) where {R,N} =
    ntuple(i -> ntuple(d -> d == cld(i, R) ? mod1(i, R) : 0, Val(N)), Val(N * R))

Anisotropic(radii)

A box with its own radius per direction: |δ_d| ≤ radii[d], origin excluded. A direction with radius 0 contributes no offset of its own, which is how a stencil is confined to a subset of the directions.

Axial(r = 1)

Axis-aligned neighbours out to r cells: ±k·ê_d for every direction d and k = 1…r. 2·N·r offsets. Axial(1) is the classic 4-in-2-D / 6-in-3-D face stencil; beyond radius 1 "face" would be a misnomer, since these are the axis lines rather than the touching faces.

CellRadius(r)

A stencil extent of r cells, as distinct from a physical MetricBall. Stencil constructors accept a bare Integer too; this exists to name the distinction where it matters.

Custom(offsets)

An explicit offset set, e.g. Custom(((1, 0), (0, 1))) for a forward-only two-point stencil. The origin is rejected: a cell is not its own neighbour.

Diagonal(r = 1)

Pure diagonals only: every offset whose components are all ±k for a single k = 1…r. 2^N·r offsets, and no axis-aligned neighbour among them.

MetricBall(r)

Every cell within physical distance r, measured through the geometry rather than in cells. Distinct from a CellRadius: on a stretched or spherical grid the number of cells within r varies across the grid, so this cannot be reduced to a fixed offset set.

Query it with Connectivity.neighbors_within! / nneighbors_within / neighbors_within, which take it as their ball keyword.

Moore(r = 1)
Vertex(r = 1)

Every cell in the surrounding box: 0 < |δ|_∞ ≤ r, i.e. (2r+1)^N - 1 offsets. Moore(1) is the 8-in-2-D / 26-in-3-D vertex stencil.

Alias for Moore, under the name the face/vertex pairing uses.

VonNeumann(r = 1)

Every cell within r steps measured in the $L^1$ (taxicab) metric: 0 < |δ|₁ ≤ r. Equal to Axial at r = 1 and strictly larger beyond it, since it admits diagonal combinations whose step count still fits.

fold_offsets(f, init, stencil, Val(N))

Fold f(acc, δ) over the stencil's offsets, unrolled at compile time.

The accumulator is threaded through as a value rather than mutated, so nothing is captured and nothing is boxed — which is what a counting or filling kernel needs, and what a closure over a mutated local would cost.

foreach_offset(f, stencil, Val(N))

Apply f to each of the stencil's offsets, with the loop unrolled at compile time.

Each offset reaches f as a register-sized NTuple{N,Int} and the offset set is never materialized on the heap — unlike offsets, whose returned tuple is allocated once it outgrows a register. Every bulk neighbour kernel goes through this.

nstencil(stencil, Val(N)) -> Int

How many offsets the stencil has in N dimensions — the buffer length a neighbors! call needs.

offsets(stencil, Val(N)) -> NTuple{K,NTuple{N,Int}}

The stencil's offsets in N dimensions, as a tuple built at compile time so a loop over it unrolls.

Offsets come out in column-major order of the enclosing box (direction 1 varying fastest), which puts Axial(1) and Moore(1) in the conventional order.

This is the one method a new shape defines; see AbstractStencil.

reach(stencil, Val(N)) -> NTuple{N,Int}

The stencil's extent in cells along each direction: maximum(|δ_d|) over its offsets. This is the halo width a traversal must leave, and the window a distance query has to scan.