Stencils & Connectivity

Connectivity answers: which cells are neighbours?

Stencil shapes

A stencil is a neighbourhood shape in index space. It carries no dimensionality — the same shape applies to a 2-D and a 5-D grid — and its shape and radius live in its type, so the offset set is built at compile time and a loop over it unrolls.

using FlowGeometries: FlowGeometries as FG
S = FG.Stencils

S.offsets(S.Axial(1), Val(2))

((-1, 0), (1, 0), (0, -1), (0, 1))

S.offsets(S.Moore(1), Val(2))

((-1, -1), (0, -1), (1, -1), (-1, 0), (1, 0), (-1, 1), (0, 1), (1, 1))

Axial(1) and VonNeumann(1) are the same four-in-2-D set; beyond radius 1 they differ, because VonNeumann admits diagonal combinations whose step count still fits.

S.nstencil(S.Axial(2), Val(2)), S.nstencil(S.VonNeumann(2), Val(2)), S.nstencil(S.Moore(2), Val(2))

(8, 12, 24)

Any dimension:

S.nstencil(S.Axial(1), Val(5)), S.nstencil(S.Moore(1), Val(5))

(10, 242)

S.reach(S.Anisotropic((3, 1)), Val(2))     # halo width per direction

(3, 1)

Your own shape

A shape supplies Stencils.offsets and nothing else. Put whatever the offsets depend on in the type, so the tuple is inferable — then nstencil, reach, foreach_offset, fold_offsets and every neighbour query work on it, with the loop still unrolled and nothing allocated.

struct Upwind{R} <: S.AbstractStencil end
Upwind(r::Integer) = Upwind{Int(r)}()
S.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))

S.offsets(Upwind(2), Val(2)), S.nstencil(Upwind(2), Val(3)), S.reach(Upwind(2), Val(2))

(((1, 0), (2, 0), (0, 1), (0, 2)), 6, (2, 2))

gu = FG.Grids.StructuredGrid(FG.Geometry.CartesianGeometry(),
                             range(0.0; step = 1.0, length = 6),
                             range(0.0; step = 1.0, length = 5))
FG.Connectivity.nneighbors(gu, 2, 2; stencil = Upwind(1)),
FG.Connectivity.nneighbors(gu, 6, 5; stencil = Upwind(1))

(2, 0)

Two different things get called a radius, and they are separate types. Stencils.CellRadius counts cells in index space; Stencils.MetricBall is a physical distance measured through the geometry. On a stretched or spherical grid the number of cells within a given distance varies across the grid, so the two cannot be collapsed into one.

Neighbourhoods by distance

A Stencils.MetricBall is queried with Connectivity.neighbors_within and its buffer/counting forms — every cell whose centre lies within the given physical distance, under the geometry's own metric:

geo = FG.Geometry.SphericalGeometry()            # Earth radius, metres
λ = range(0, 2π; length = 25)[1:24]
φ = range(-π/2, π/2; length = 13)
g = FG.Grids.StructuredGrid(geo, λ, φ)
ball = S.MetricBall(2.0e6)                       # everything within 2000 km
FG.Connectivity.nneighbors_within(g, 5, 7; ball = ball)

4

The same call near a pole finds many more cells, because 2000 km spans every longitude there — which is exactly why this cannot be a stencil:

FG.Connectivity.nneighbors_within(g, 5, 13; ball = ball)

47

A bare number works as the ball, and the buffer form follows the neighbors! pattern — size with the count, fill in place:

n = FG.Connectivity.nneighbors_within(g, 5, 7; ball = 2.0e6)
buf = Vector{Int}(undef, n)
FG.Connectivity.neighbors_within!(buf, g, 5, 7; ball = 2.0e6)
buf

4-element Vector{Int64}:
 125
 148
 150
 173

The candidate window is Connectivity.metric_window — O(1) per direction on any separable axis, uniform or stretched, given the minimum steps Connectivity.MetricTopology carries — and each candidate is then kept or dropped by the geometry's distance, so the result is a genuine metric ball: great-circle on a sphere, Vincenty on a spheroid, the chord where a radial or height direction is present. Periodic directions wrap by minimum image, so a ball sitting on the seam is the same ball as anywhere else:

FG.Connectivity.neighbors_within(g, 1, 7; ball = 1.5e6) ==
    FG.Connectivity.neighbors_within(g, 1, 7; ball = S.MetricBall(1.5e6))

true

Distance between two cells

Geometry.distance also takes two cell indices, resolving them through the grid's own topology. Across a periodic seam it takes the short way round:

gper = FG.Grids.StructuredGrid(FG.Geometry.CartesianGeometry(),
                               range(0.0; step = 1.0, length = 10),
                               range(0.0; step = 1.0, length = 6);
                               periodic = true, period = 10.0)
FG.Geometry.distance(gper, (1, 1), (10, 1)), FG.Grids.displacement(gper, (1, 1), (10, 1))

(1.0, (-1.0, 0.0))

Grids.displacement is the signed per-direction offset that distance was taken from. It is a coordinate quantity, so it sits in Grids, while the distance extends Geometry.distance.

Convolutions over every image

A neighbour set visits each cell once, at its nearest image. A convolution on a torus cannot: with period L,

\[\\bar f(x) = \\int K(x-y) f(y)\\,dy = \\sum_k \\int_\\text{cell} K(x-y-kL) f(y)\\,dy\]

so once the kernel support exceeds L/2 one cell contributes through several images at different displacements. Connectivity.fold_within takes the convention as an argument — NearestImage or AllImages — and AllImages also widens the window past metric_window's one-turn cap, which excludes those images.

Nx, Δx = 32, 62.5; Lx = Nx * Δx
axx = range(0.0; step = Δx, length = Nx)
gt = FG.Grids.StructuredGrid(FG.Geometry.CartesianGeometry(), axx, axx;
                             periodic = (true, true), period = (Lx, Lx))
count_within(conv, rad) = FG.Connectivity.fold_within(
    (a, J, d) -> a + 1, 0, gt, 1, 1; ball = rad, images = conv)
count_within(FG.Connectivity.NearestImage(), 3500.0),
count_within(FG.Connectivity.AllImages(), 3500.0)

(1023, 9844)

Below L/2 the two agree exactly. Above it, filtering one Fourier mode with a Gaussian of width ℓ = L reproduces the analytic transfer exp(-k²ℓ²/4α) to roundoff under AllImages; under nearest-image the error stays at 40% however far the support is widened, since a dropped image is not recovered by searching further.

Summing images asserts that a periodic direction is a translation of the domain. On a sphere it is an identification — λ and λ+2π are the same point — so AllImages raises there, where summing images counts one cell repeatedly.

self = true folds the centre cell too, at distance zero. A neighbour set excludes it; a convolution needs it, and it carries the kernel's largest weight.

To query every cell, materialize the whole graph once instead — Connectivity.build_connectivity_within is the ball analogue of build_connectivity, row k holding exactly what the per-cell query returns for cell k, and symmetric because the metric is:

csr = FG.Connectivity.build_connectivity_within(g; ball = 2.0e6)
FG.Connectivity.nedges(csr), FG.Connectivity.is_symmetric_adjacency(csr)

(5184, true)

What a ball query reads about the grid

The observation the stencil side of the module rests on is that a neighbour computation never looks at a coordinate. It reads three things — extent per dimension, wrapping per dimension, and which cells are active — and nothing else. That triple is IndexTopology, under which a curvilinear grid needs no separate implementation from a structured one: it is the N = 2 case of the same algorithm.

A ball query is the one thing that cannot work that way. It needs coordinates, and the smallest step per direction, which bounds the candidate window on a stretched axis. That is Connectivity.MetricTopology, the same idea for the metric path. It is O(1) to build: it reads the per-axis reductions through Grids.minimum_spacing and Grids.bounds, which every layout answers without a scan:

mt = FG.Connectivity.MetricTopology(g)
FG.Connectivity.nneighbors_within(g, 5, 7; ball = ball, topology = mt)

4

So passing it changes nothing and omitting it costs nothing — the default is already free. There is no hoisting to remember here.

The spatial index is different: it is worth hoisting, and building one costs O(n), so it pays back over many queries. Curvilinear and node grids have no separable axes for a window to bound, so an unindexed query tests every cell.

Grids.cell_list is the one to reach for, and it needs no package at all. It bins the cell centres at the radius you mean to query at, and a query visits the bins its ball can reach:

gc = FG.Grids.HEALPixGrid(8)
top = FG.Connectivity.MetricTopology(gc; index = FG.Grids.cell_list(gc; ball = 2.0e6))
FG.Connectivity.nneighbors_within(gc, 1; ball = 2.0e6, topology = top)

18

Every cell lands in exactly one bin — periodicity wraps the bin coordinate, leaving the point where it is — so a query emits each cell once and needs no buffer to deduplicate into. It therefore runs inside a kernel, and the sweeps build one by default.

Connectivity.indexed is the alternative, a k-d tree, and needs NearestNeighbors:

using NearestNeighbors                                    # loads the extension
gu = FG.Grids.HEALPixGrid(8)
ixu = FG.Connectivity.indexed(gu)
scratch = FG.Connectivity.ball_scratch()
FG.Connectivity.nneighbors_within(gu, 1; ball = 2.0e6, topology = ixu, scratch = scratch)

18

The index only ever returns a superset of the ball — the exact distance ≤ r gate still decides membership — so an indexed query and a scan return the same cells, and loading the extension changes speed and nothing else.

A neighbour list is a set: the cells come back in whatever order enumerated them, and no query sorts, keeping an O(m) query O(m). Connectivity.sort_neighbors! sorts on request. Connectivity.ball_scratch is the candidate buffer, one per task; a query given none allocates its own.

Sweeping every cell

Better still, do not write the loop. Connectivity.foreach_within and Connectivity.mapreduce_within walk every cell's ball and build the topology and the index once for the whole sweep, making them O(n log n) where a per-cell loop that rebuilds the index is O(n²):

FG.Connectivity.mapreduce_within((I, J, d) -> 1, +, 0, g; ball = 1.0e6)   # total pairs within 1000 km

1392

FG.Connectivity.mapreduce_within((I, J, d) -> d, +, 0.0, g; ball = 1.0e6) # and the total distance

2.0322300096525317e8

foreach_within is the same traversal for an f that writes; under a threaded backend it runs on disjoint spans of cells, so what f writes must be determined by I. build_connectivity_within does the same hoisting internally.

The k nearest

A ball asks "everything within r". The other question is "the nearest k", and Connectivity.k_nearest answers it exactly under the same metric, on every architecture:

idx, dist = FG.Connectivity.k_nearest(g, 5, 7; k = 4)
dist

4-element Vector{Float64}:
 1.6679238996683806e6
 1.6679238996683809e6
 1.6679238996683809e6
 1.667923899668382e6

It widens a ball until k cells have been seen and keeps the k smallest in a bounded heap, so the result never depends on the starting radius and nothing is materialized along the way. Equal distances resolve by linear index, so an indexed query and a scan agree exactly. Connectivity.k_nearest! writes into your own buffers and allocates nothing.

Seeding from a coordinate

Every query above starts at a cell. Observational data does not arrive that way — a station, a ship track or a float has a coordinate, and the cell it belongs to is part of the question. Passing a coordinate in place of the cell indices asks the same questions about a point — written as a Tuple, NamedTuple, AbstractVector or SVector, as anywhere else a point is taken:

FG.Grids.locate(g, (0.4, 0.1))                    # the cell the point falls in

(3, 7)

length(FG.Connectivity.neighbors_within(g, (0.4, 0.1); ball = 1.0e6)),
FG.Connectivity.k_nearest(g, (0.4, 0.1); k = 3)[2]

(0, [1.0118794641175776e6, 1.0856083232992077e6, 1.2888499644977713e6])

Connectivity.fold_at is the fold behind them, as Connectivity.fold_within is for a cell. There is no seed cell to skip, so unlike the cell-seeded form every cell within the ball is visited — including the point's own.

locate is the containing cell on a rectilinear grid, per direction and wrapping a periodic one, and the nearest cell centre elsewhere. For a node set those coincide, its cells being the Voronoi regions of its nodes; on a strongly sheared curvilinear grid they can differ, so the contract is nearest-centre.

The traversal is the cell-seeded one with the window widened by how far the point sits from its cell's centre, keeping it O(1) per direction. Off the rectilinear grids the cost depends on the index: Grids.cell_list answers "which bins does this point's ball reach" directly, so a point query visits only those bins.

The ball, and the part of it you can get to

A ball is not a connected patch. With a mask, or a concave domain, it can contain cells that are close to the seed in space but reachable from it only by leaving the ball. Both sets are useful and they are different algorithms, so reach names which one you get:

wall = trues(9, 9); wall[:, 5] .= false                  # a barrier through the middle
gw = FG.Grids.StructuredGrid(FG.Geometry.CartesianGeometry(),
                             collect(0.0:8.0), collect(0.0:8.0), wall)
length(FG.Connectivity.neighbors_within(gw, 5, 3; ball = 3.5)),
length(FG.Connectivity.neighbors_within(gw, 5, 3; ball = 3.5,
                                        reach = FG.Connectivity.Connected(S.Moore(1))))

(28, 25)

Connectivity.Unrestricted is every cell within the radius, and the default. Connectivity.Connected is the connected component of that set containing the seed: it walks adjacency and prunes at the ball's edge, so it is a subset. The two agree whenever the ball is connected under the adjacency — always so on a maskless Cartesian StructuredGrid, where each index step moves monotonically in one coordinate — and Connected is strictly smaller wherever something separates two parts of the ball.

Neither is a cheaper version of the other. Take cells P (the seed), Q and R, adjacent only as P–Q–R, with d(P,Q) = 1.2r and d(P,R) = 0.8r. Walking outward from P and dropping anything farther than r stops at Q, so it never reaches R, which is inside the ball. Such a walk computes Connected; the ball is reached only by a spatial query.

Adjacency has to be named, because only a node set carries its own: Connected() means direction-1 adjacency on the index-space architectures and the stored neighbour lists on an UnstructuredGrid, and Connected(stencil) sets it explicitly for the former. On the latter a stencil has no index space to mean anything in, so it raises.

Querying neighbours

Three forms, in increasing order of how much they allocate:

out = Vector{Int}(undef, 8)
n = FG.Connectivity.neighbors!(out, grid, i, j; stencil = FG.Stencils.Moore(1))   # writes n indices
k = FG.Connectivity.nneighbors(grid, i, j)                            # just the count
it = FG.Grids.neighbors(grid, i, j)                            # lazy iterator, no array

FlowGeometries.Connectivity.StencilNeighbors{FlowGeometries.Grids.StructuredGrid{Float64, FlowGeometries.Geometry.SphericalGeometry{Float64}, 2, FlowGeometries.SphericalSampling.ClenshawCurtisSampling, Tuple{FlowGeometries.Grids.Periodic, FlowGeometries.Grids.Bounded}, Tuple{Vector{Float64}, Vector{Float64}}, FlowGeometries.Grids.SeparableMeasure{Float64, 2, Tuple{Vector{Float64}, Vector{Float64}}}, FlowGeometries.Grids.AllActive{2}}, 2, FlowGeometries.Stencils.Axial{1}}(StructuredGrid{Float64}(31×16), (3, 5), FlowGeometries.Stencils.Axial{1}(), true)

stencil is any Stencils shape — Axial(r), VonNeumann(r), Moore(r) (alias Vertex), Diagonal(r), Anisotropic(radii) or Custom(offsets) — at any radius and in any number of dimensions. A stencil is named by its type, so the neighbour iterator built from it resolves at compile time and allocates nothing per cell. active_only = true excludes masked-out cells from both ends of an edge.

Prefer neighbors! in hot loops: neighbors returns a lazy iterator, and neighbors! writes into a buffer you own.

CSR

For a whole graph at once:

conn = FG.Connectivity.build_connectivity(grid)                   # or (sampling, args...)
FG.Connectivity.nnodes(conn), FG.Connectivity.nedges(conn)
FG.Grids.neighbors(conn, i)                                # a view into the flat array

3-element view(::Vector{Int32}, 7:9) with eltype Int32:
  2
  4
 34

CSRConnectivity stores one flat neighbour array plus offsets, so a traversal reads contiguous memory and the whole graph is two allocations. Both buffers accept any Integer, and a builder picks the narrowest that holds both a node id and an offset: Int32 below two billion of either, halving the graph's memory and the bandwidth every traversal of it costs. Ti names a width explicitly.

Building from a sampling never constructs a grid:

FG.Connectivity.build_connectivity(FG.SphericalSampling.GaussLegendreSampling(), 1000)

FlowGeometries.Connectivity.CSRConnectivity{Vector{Int32}, Vector{Int32}}(Int32[1999, 2, 2000, 1, 3, 2001, 2, 4, 2002, 3  …  1996998, 1998997, 1998999, 1996999, 1998998, 1999000, 1997000, 1998999, 1997002, 1997001], Int32[1, 4, 7, 10, 13, 16, 19, 22, 25, 28  …  7991976, 7991979, 7991982, 7991985, 7991988, 7991991, 7991994, 7991997, 7992000, 7992003])

The neighbour graph of a tensor-product sampling is fixed by its axis lengths and longitude wrapping alone, so the axes are never evaluated. For Gauss–Legendre, evaluating them is an O(n²) root solve whose result the graph does not depend on.

Index topology directly

topo = FG.Connectivity.IndexTopology((nlon, nlat), (true, false), nothing)   # mask = nothing → all active
conn = FG.Connectivity.build_connectivity(topo; stencil = FG.Stencils.Axial(1))

FlowGeometries.Connectivity.CSRConnectivity{Vector{Int32}, Vector{Int32}}(Int32[31, 2, 32, 1, 3, 33, 2, 4, 34, 3  …  462, 493, 495, 463, 494, 496, 464, 495, 466, 465], Int32[1, 4, 7, 10, 13, 16, 19, 22, 25, 28  …  1896, 1899, 1902, 1905, 1908, 1911, 1914, 1917, 1920, 1923])

Useful when you have a shape and a wrapping rule but no grid, and enough to drive every stencil query in the package.

Dense and sparse adjacency

FG.Connectivity.adjacency_matrix(conn)                       # Matrix{Bool}
FG.Connectivity.adjacency_matrix!(A, conn)                   # into your buffer
# needs SparseArrays; see the Extensions page
# FG.Connectivity.sparse_adjacency_matrix(grid)

496×496 BitMatrix:
 0  1  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 1  0  1  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  1  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  1  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  0  1  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  0  1  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  1  0  1  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  1  0  1  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  1  0  1  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  1  0  1  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮        ⋱     ⋮              ⋮              ⋮
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  1  0  1  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  1  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  1  0  1  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  1  0  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  1  0  1  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  1  0  1  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  1  0  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  1  0  1
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  1  0

The sparse path builds CSC straight from the neighbour list by one counting pass and one placement pass — no coordinate triples, no sort, no permutation vector, and row indices come out ascending for free because the placement walks nodes in order.

For a symmetric adjacency with sorted rows the CSR and CSC arrays are the same arrays, so sparse_adjacency_matrix(grid) hands the connectivity's own buffers to the matrix, with no transpose and no second copy.

Whether that holds is read from the layout before the graph exists. Under an index stencil it is the stencil's own symmetry; a formula layout answers with Grids.has_symmetric_adjacency, which HEALPix, the cubed sphere, the icosahedral geodesic and Yin–Yang declare. A ring grid does not: adjacent rings of unequal width make its straddling relation directed, so it takes the transpose.

FG.Grids.has_symmetric_adjacency(FG.Grids.HEALPixGrid(4))   # licenses the shortcut, before the build
FG.Connectivity.is_symmetric_adjacency(conn)                # the same question of a built graph
FG.Connectivity.sort_neighbors!(conn)                       # order each node's block ascending

FlowGeometries.Connectivity.CSRConnectivity{Vector{Int32}, Vector{Int32}}(Int32[2, 31, 32, 1, 3, 33, 2, 4, 34, 3  …  494, 463, 493, 495, 464, 494, 496, 465, 466, 495], Int32[1, 4, 7, 10, 13, 16, 19, 22, 25, 28  …  1896, 1899, 1902, 1905, 1908, 1911, 1914, 1917, 1920, 1923])

Threading

Every builder takes an opt-in backend; see Performance.

using ComputationalBackends: ThreadedBackend
FG.Connectivity.build_connectivity(grid; backend = ThreadedBackend())

Results are bit-identical to serial; the test suite asserts it with ==.

Sampling-specific topology

FG.Connectivity.build_connectivity(FG.SphericalSampling.HEALPixSampling(64))              # RING face-table adjacency
FG.Connectivity.build_connectivity(FG.SphericalSampling.CubedSphereSampling(), n)         # with the gnomonic seam fold
FG.Connectivity.build_connectivity(FG.SphericalSampling.YinYangSampling(), nlon, nlat)
FG.Connectivity.healpix_neighbors!(out, nside, 0)             # one pixel; 0-based pixel ids

8

HEALPix RING neighbours follow the standard face-table algorithm (Górski et al. 2005; Reinecke 2003). The offsets are walked in ring order, so the emitted ids arrive already ascending for ~99% of pixels and the dedup pass has almost nothing to move.

Mask topology

The mask's own shape is a grid property, so it is answered here.

C = FG.Connectivity
geo = FG.Geometry.CartesianGeometry()
m = trues(9, 9)
m[4:6, 4:6] .= false          # an enclosed block
g = FG.Grids.StructuredGrid(geo, 0.0:1.0:8.0, 0.0:1.0:8.0, m)

C.count_holes(g), C.connected_components(g)[2]

(1, 1)

count(C.interior(g)), count(C.boundary_cells(g)), count(m)

(28, 44, 72)

interior is the active cells whose whole stencil is active and in range; boundary_cells is the rest of the active set, and the two partition it. count_holes counts the connected inactive regions fully enclosed by active cells — an estimate of the active region's first Betti number. A region that reaches a non-wrapping edge counts as outside, so wrapping a direction can turn an outside region into an enclosed one:

gp = FG.Grids.StructuredGrid(geo, 0.0:1.0:8.0, 0.0:1.0:8.0, m;
                             periodic = (true, true), period = (9.0, 9.0))
C.count_holes(gp)

1

All of it is dimension-generic:

m3 = trues(7, 7, 7); m3[3:5, 3:5, 3:5] .= false
C.count_holes(FG.Grids.StructuredGrid(geo, 0.0:1.0:6.0, 0.0:1.0:6.0, 0.0:1.0:6.0, m3))

1

Index helpers

FG.Connectivity.linear_index(grid, i, j)      # column-major, i fastest
FG.Connectivity.cartesian_index(grid, lin)    # the inverse

CartesianIndex(7, 1)