Operators

AbstractMaskPolicy

What apply_stencil! does at the edge of the active region: BlankMasked, ShiftWithinRun or ReduceInRun.

A type, like the image and reach conventions in Connectivity: which cells carry a number and which carry masked is visible in the call and fixed at compile time.

BlankMasked()

Write masked at a cell that is inactive or whose stencil reads an inactive cell. The default, and the only policy that never invents a value: where the stencil cannot be formed from active data, there is no derivative.

Its cost is a dead band. Every active cell within nodes - 1 of a masked cell is blanked, so a five-point derivative loses two cells either side of every coastline.

GradientPlan

The geometry of a least-squares gradient, separated from any field: for each cell, the coefficient on each neighbour's difference from it, once per coordinate direction. Built by gradient_plan, applied by gradient!.

The two index buffers are typed independently and hold any Integer, so a plan built off a CSR keeps that CSR's own width — see Connectivity._index_type.

D is the number of directions the fit resolves, which is the grid's coordinate count: two on a surface — a (λ, φ) or (x, y) mesh — and three in a volume, where the neighbourhood spans a solid angle.

apply_stencil! covers the separable case. A CurvilinearGrid has no separable axis to build a scalar stencil along, and a node set's neighbours come from connectivity, so this is the gradient for both.

The construction, with tangent-plane displacements Δrₖ from Geometry.project_to_tangent_plane, differences Δfₖ = fₖ - f₀ and weights wₖ = 1/|Δrₖ|²: minimising Σₖ wₖ (∇f·Δrₖ - Δfₖ)² gives A ∇f = b with A = Σₖ wₖ Δrₖ⊗Δrₖ and b = Σₖ wₖ Δrₖ Δfₖ. A depends only on the geometry, so it is inverted once here and the per-neighbour coefficients A⁻¹ wₖ Δrₖ are what the plan stores; each apply is then one dot product per cell and allocates nothing.

What this buys over inverting an index-space Jacobian:

• exact for a linear field on any stencil, however skewed — the least-squares combination cancels the leading truncation term, which a 2×2 Jacobian inverse does not;
• second order on locally symmetric stencils, degrading toward first on strongly skewed cells;
• it reduces to the centred difference where the stencil is separable and orthogonal (A diagonal), so it agrees with apply_stencil! where both apply.

Where A is rank deficient — a direction with no data, at a boundary or beside a mask — the pseudo-inverse zeroes that component, under the rule apply_stencil! states for a mask: a value the active data does not determine is not produced.

ReduceInRun()

ShiftWithinRun, and where the run cannot hold nodes, use the largest window it can, down to order + 1 samples. masked below that, where no derivative of that order exists.

This trades accuracy order for coverage — a five-point scheme becomes three-point in a strait three cells wide — so it carries its own name. Ask for it when a value everywhere matters more than a uniform order.

Under this policy nodes is a ceiling, at the end of the axis as well as the end of a run: an axis with fewer than nodes samples uses as many as it has, and one with fewer than order + 1 is masked throughout. A single-latitude strip, a two-level column and a one-cell-wide channel are ordinary grids, and this asks for "second order where the axis allows it" with no clamping of nodes at the call site. The other two policies raise there, neither of them degrading.

ShiftWithinRun()

Shift the stencil to fit inside the run of active samples containing the cell, keeping the full node count — the same thing the stencil already does at the end of a bounded axis, with the end of the active run as the boundary. masked only where the run is shorter than nodes.

The accuracy order is therefore the same everywhere a value is written, as it is under Discretization.fd_weights's inward shift at a bounded end. On a run of at least nodes active samples the weights are identical to the unmasked ones, so the interior of an active region is bit-for-bit unchanged.

Discretization.axis_stencils(grid, dim; order=1, nodes=order+1) -> (indices, weights)

Discretization.axis_stencils for direction dim of grid, taking that direction's axis and wrap period from the grid.

The table depends on the grid alone, so a caller differencing many fields along the same direction builds it once and hands it to the (out, field, grid, indices, weights, dim) form. The (out, field, grid, dim) form above rebuilds it on every call.

Discretization.stencil_plan(grid, dim; order=1, nodes=order+1) -> AbstractStencilPlan

Discretization.stencil_plan for direction dim of grid, taking that direction's axis and wrap period from the grid.

The form to hold in a loop over fields. Which plan it is follows the axis's own spacing, so a uniform direction gets the register-resident weights and a stretched one the table.

_bary_weights(grid, geo, p0, nearest, active_only) -> (a, b, c, (w₁, w₂, w₃))

The cell of the grid's Grids.CellMesh containing p0, and the barycentric weights of p0 in it. a == 0 where there is none — outside the mesh, or against an inactive node.

Trying the cells incident on nearest is exhaustive: p0 lies in the Voronoi cell of its nearest node, and in a Delaunay triangulation that region is covered by the triangles incident on that node.

The weights are computed in the tangent plane at p0, where the three vertices sit at their Geometry.local_displacements and p0 itself is the origin.

_curl_location(sg, i) -> NTuple{N,AbstractLocation}

Where (∇×u)_i lives: Discretization.Center in direction i and Discretization.Face in every other — the vorticity point of an Arakawa C cell for that component. In two directions the scalar curl is the i = 3 case, giving (Face, Face), the corner.

_face_span(nc, periodic) -> Int

How many samples a direction of nc cells carries at Discretization.Face: nc where it wraps, its last face being its first, and nc + 1 where it does not.

_fusable(out, field, geo, mask) -> Bool

Whether the metric factor can be applied inside the sweep, which needs it to be constant across each contiguous direction-1 run the sweep writes.

For a geometry that declares direction 1 metric-invariant the factor depends on directions 2:N, so it is constant on every such run, whichever direction is differenced — the sweep scales run by run. A Cartesian metric is the identity and has nothing to apply, and the run addressing needs the linear layout.

_fuse_scale!(out, invh, R, start, nruns, runlen, rbase, masked) -> nothing

Scale a span the sweep has just written, run by run, while it is still in cache.

The span holds nruns contiguous runs of runlen cells along direction 1, and the metric factor is constant on each of them, direction 1 being metric-invariant. invh[1:R] holds one factor per spatial run in column-major order over directions 2:N, and a batch element repeats that cycle, so the run number is taken modulo R. rbase is the span's first run, counted the same way.

invh comes from _metric_scratch and is held across calls, so it may be longer than R: the count is the argument, never length(invh).

_mask_may_hide(grid, policy, active_only) -> Bool

Whether the neighbourhood test has anything to find: under BlankMasked an active-only k_nearest may have skipped an inactive cell nearer than the farthest it kept, and the test is a second range query over the same index. On an Grids.AllActive grid there is no inactive cell to skip, so the answer is known without the query.

_metric_row_factors!(buf, grid, Val(dim), Val(N)) -> Int

Write the inverse scale factor of every direction-1 run into buf, in column-major order over directions 2:N — the order the sweep writes them in, for any differenced direction. Returns how many.

Direction 1 is metric-invariant here (see _fusable), so one factor covers a whole run and any axis-1 coordinate serves as the point's first component.

A degenerate factor — one at or below the geometry's Discretization.metric_floor — is stored as zero, and _scale_span! blanks a run whose factor is zero.

_metric_scratch(T, n) -> Vector{T}

A task-held buffer of at least n factors, as Connectivity.ball_scratch is for a ball query.

The fused sweeps are serial and call no user code, so one buffer per task is enough; it grows to the largest grid the task has swept and is refilled per call.

_nodecount(nodes) -> Int

The node count, from either a Val or a plain Int.

One body then serves both. With a Val the trip count is a literal, so the loop unrolls and the weights reach registers; with an Int, the form a node count above the specialized set takes, it is an ordinary loop.

_plan_sweep_host!(out, field, plan, mask, masked, Val(dim), Val(N), Val(M), invh, R) -> out

The host sweep for a plan. Same nest as _stencil_sweep_host! — Cartesian range walked directly, nest split at dim, node count in the type — with a uniform plan's interior row a tuple in registers and the shifted end rows read from its own O(K²) table.

invh fuses the metric: each span is scaled as it is written, while it is in cache, saving the second full pass over out a separate _scale_by_metric! costs. Pass nothing for the plain sweep.

_scale_span!(out, start, len, inv_h, masked) -> nothing

Multiply the contiguous span out[start:start+len-1] by inv_h, or write masked across it when inv_h is zero — which is how a degenerate scale factor reaches here, derivative! having compared it against the geometry's own Discretization.metric_floor.

_stencil_sweep_host!(out, field, indices, weights, mask, masked, Val(dim), nodes, Val(N),
                     Val(M), invh, R) -> out

The host sweep. The index-parallel form serves a device launch; this one adds three things a per-work-item body cannot express:

• iterate the Cartesian range directly, with no index recovered per cell from a linear one;
• split the nest at dim, hoisting the stencil row — which depends only on the index along dim — out of the contiguous inner loop wherever the differenced direction is a slower-varying one;
• carry the node count in the type, so the innermost loop has a known trip count, unrolls, and holds the weights in registers.

The arithmetic and its order are identical to _stencil_cell!, so the two paths agree bit for bit.

invh fuses the metric factor into the sweep — see _fuse_scale!. Pass nothing for the plain sweep.

_sym_eigvals3(A, tol) -> (λ₁, λ₂, λ₃)

Eigenvalues of a symmetric 3×3, largest first, by the closed form for the characteristic cubic: shift by the mean eigenvalue, scale, and read the three roots off cos at thirds of a turn.

Eigenvalues only. A closed-form eigenvector of a 3×3 loses accuracy when two eigenvalues are close, and _sympinv3 needs none: a pseudo-inverse of a symmetric matrix is a polynomial in that matrix, and its coefficients are functions of the eigenvalues alone.

_sympinv(A::NTuple{D,NTuple{D,T}}, tol) -> NTuple{D,NTuple{D,T}}

Pseudo-inverse of a symmetric positive-semidefinite A given as its rows, dropping any eigendirection whose eigenvalue is at or below tol. D = 2 and D = 3, which are the tangent plane and the volume.

A direction is dropped, never regularised — see _sympinv2, which this calls at D = 2.

_sympinv2(a, b, c, tol) -> (p11, p12, p22)

Pseudo-inverse of the symmetric 2×2 [a b; b c], dropping any eigendirection whose eigenvalue is below tol. Closed form: a 2×2 symmetric eigenproblem has one.

A direction is dropped, never regularised. A stencil that carries no information along some tangent direction — every neighbour on one line, which happens at a boundary and beside a mask — leaves A singular there, and the data does not determine that component of the gradient. Inverting a nudged matrix answers with a number governed by the nudge.

_sympinv3(A, tol) -> NTuple{3,NTuple{3,T}}

Pseudo-inverse of a symmetric positive-semidefinite 3×3.

At every rank A⁺ is a polynomial in A with no constant term, so this is closed form and needs no eigenvector. Such a polynomial annihilates the null space and acts as 1/λ on the range, which is the Moore–Penrose conditions:

• rank 3: adj(A)/det(A), the ordinary inverse;
• rank 2: αA² + βA with α = -(λ₁+λ₂)/(λ₁λ₂)² and β = (λ₁² + λ₁λ₂ + λ₂²)/(λ₁λ₂)²;
• rank 1: A/λ₁², since A = λ₁vvᵀ there;
• rank 0: zero.

apply_stencil!(out, field, grid, indices, weights, dim; order=1, active_only=true,
               masked=zero, policy=BlankMasked(), backend=nothing) -> out

Apply a stencil table built by Discretization.axis_stencils. The mask, the wrap period and the axis all come from grid.

The axis coming too, any mask policy works here. The bare (indices, weights) form has no axis to rebuild a window from at a mask edge, so it accepts only BlankMasked.

This is the form to use in a loop over fields: the table is the same for all of them, and it is the one part of the work that depends on the grid alone.

apply_stencil!(out, field, grid, dim; order=1, nodes=order+1, active_only=true, masked=zero) -> out

apply_stencil! with the axis, wrap period and mask taken from grid, so a periodic direction wraps and an inactive cell is honoured without restating any of it.

Only a rectilinear direction has a 1-D axis to difference along, so this is a StructuredGrid method.

apply_stencil!(out, field, x, indices, weights, dim; order=1, period=nothing, mask=nothing,
               masked=zero, policy=BlankMasked(), backend=nothing) -> out

Apply a table built by Discretization.axis_stencils and keep the axis, so any mask policy works.

The table depends on the axis alone, so a caller differencing many fields along one direction builds it once. Degrading at a mask edge needs the axis to rebuild a window from, so the bare (indices, weights) form accepts only BlankMasked; this form takes both and serves every policy.

The split is the one the degrade path makes internally: the precomputed row is used wherever the window is intact, which is every cell away from a mask, and the axis is touched only where a window is rebuilt.

Building the table is O(n) against an O(n²) apply, so holding it across fields matters most on small grids. It also removes an O(n · nodes) allocation from every call, at any size.

apply_stencil!(out, field, x, dim; order=1, nodes=order+1, period=nothing,
               mask=nothing, masked=zero) -> out
apply_stencil!(out, field, indices, weights, dim; mask=nothing, masked=zero) -> out

Apply a weight set along direction dim of field, writing out[I] = Σ_q weights[I[dim], q] · field[…, indices[I[dim], q], …].

Every convention this needs is fixed by the package: the result sits at the same location as the input, so there is no staggering decision, and the stencil shifts inward at a bounded end and wraps on a periodic one, which is Discretization.fd_weights's stated boundary behaviour and removes the need for a halo. Operations that do take those choices live elsewhere — a staggered difference, or a multi-direction operator like a divergence or a curl, each of which needs a result location and a boundary-condition policy.

Pass the axis and an order to have the weights built for you, or precomputed indices/weights from Discretization.axis_stencils to reuse them across many fields.

With a mask, a cell is written as masked when it is inactive or when its stencil reads an inactive cell — the derivative there is not determined by the active data, so it is not invented. out and field may not alias.

apply_stencil!(out, field, plan, dim; mask=nothing, masked=zero, backend=nothing) -> out

Apply a held Discretization.stencil_plan along direction dim of field.

The primary form: the weights are built once, so a caller differencing many fields along one direction pays for them once, and nothing is allocated per call. The (out, field, x, dim; order, nodes) forms build a plan and call this, which is where their allocation comes from.

mask and masked behave as they do for the table forms: a cell whose stencil reads an inactive cell is written masked.

curl!(out, u1, u2, sg; masked = NaN) -> out
curl!(out, us, sg; masked = NaN) -> out
curl!(outs, us, sg; masked = NaN) -> outs

The curl of a C-staggered vector field, us[d] at the d-face location.

In two directions it is one scalar, written to out at the corner where both directions are at Discretization.Face:

(1/J)·[ ∂(h_2 u_2)/∂x_1 − ∂(h_1 u_1)/∂x_2 ]

In three it is a vector, component i written to outs[i] at that component's own vorticity point — see _curl_location — with (i, j, k) running cyclically:

(∇×u)_i = (1/(h_j·h_k))·[ ∂(h_k u_k)/∂x_j − ∂(h_j u_j)/∂x_k ]

Each term is one difference across one cell. Like the divergence this is the circulation around the cell over its area, so the circulation telescopes: the curl summed against the corner measures over a closed or wrapping domain vanishes to round-off. Being the same differences the gradient makes, curl(gradient(f)) is zero to round-off as well.

A point missing either difference — the outer faces of a bounded direction — is masked, as is one whose four bounding centres are not all active.

curl(u1, u2, sg; kwargs...) -> Array
curl(us, sg; kwargs...) -> Array | NTuple{3,Array}

curl! into fresh arrays: one at the corner in two directions, and one per component at its own vorticity point in three.

derivative!(out, field, grid, dim; order=1, nodes=order+1, policy=BlankMasked(),
            masked=zero, active_only=true, backend=nothing) -> out

The derivative with respect to distance along direction dim: apply_stencil! divided by the metric factor,

∂f/∂sᵈ = (1/hᵈ) · ∂f/∂ξᵈ,    hᵈ = Geometry.scale_factors(geo, p)[d]

On a Cartesian metric every hᵈ is 1 and this is apply_stencil! exactly, at no cost. Anywhere else it is the derivative a physical law is written in, singular points included.

Where the metric degenerates the derivative does not exist, and masked is written. Longitude at a pole is the case: h_λ = R cos φ → 0, so 1/h_λ diverges. The test scales with the geometry's own size and with the precision, |h| ≤ L·√eps(T). A fixed threshold serves neither: 1e-12 is below eps(Float32), and at Float32 on Earth's radius cos(Float32(π/2)) ≈ -4.4e-8 puts h_λ at the pole around 0.28 m.

No scale factor in this package depends on longitude, so hᵈ is constant along the first axis whichever direction is differenced. The scaling is applied once per remaining index and swept along that contiguous axis. (It is not generally constant along the differenced direction — on a spheroid h_φ = M(φ) varies with φ — so it is not hoisted that way.)

A divergence or a curl is still the caller's to assemble, needing a result location and a boundary policy this does not choose. Note the flux form when doing so: on a sphere

∇·u = (1/(R cos φ)) [ ∂u_λ/∂λ + ∂(u_φ cos φ)/∂φ ]

so the second term differentiates u_φ cos φ, not u_φ; taking two physical derivatives and adding them is a different, wrong expression.

derivative!(out, field, grid, indices, weights, dim; order=1, active_only=true, masked=zero,
            policy=BlankMasked(), backend=nothing) -> out

derivative! from a table the caller holds — the same reuse as the apply_stencil! form above, for the entry point a geometry-aware caller actually uses.

The metric fuses into the sweep here too, on the terms _fusable states.

derivative!(out, field, grid, plan, dim; active_only=true, masked=zero) -> out

derivative! from a held Discretization.stencil_plan.

The form to use in a loop: the weights and each row's metric factor depend on the grid alone, so both are built once here. Where the factor is constant across the span the sweep writes — see _fusable — it is applied to each row as that row is written, while it is still in cache, saving a second pass over out.

divergence!(out, us, sg; masked = NaN) -> out

The divergence of a C-staggered vector field — us[d] at the d-face location — written to out at the cell centres.

(1/J)·Σ_d [ (J/h_d)·u_d ] differenced across the cell: the finite-volume divergence, the net flux through the cell's own faces over its own volume. It is discretely conservative — summing it against the cell measures telescopes, two neighbours' shared face cancelling to the bit, so a closed or wrapping domain integrates to zero to round-off.

A centre whose bounding faces are not all active reads masked.

divergence(us, sg; kwargs...) -> Array

divergence! into a fresh array at the cell centres.

gradient!(outs::Tuple, field, plan) -> outs
gradient!(g1, g2, field, plan) -> (g1, g2)
gradient!(g1, g2, g3, field, plan) -> (g1, g2, g3)

Apply a GradientPlan: the components of ∇field at every cell, written one per output array. One dot product per cell over its neighbours, allocating nothing.

There must be exactly ncomponents of them. field and the outputs are indexed linearly, so an N-D array of the grid's shape works as is. The components are named by plan.names — (:λ, :φ) on a sphere, (:x, :y) on a plane, (:x, :y, :z) in a volume — and are per unit distance, the displacements being metric already.

gradient!(outs, f, sg; masked = NaN) -> outs

The gradient of a centre field f on a Grids.StaggeredGrid, component d written to outs[d] at the d-face location — the Arakawa C gradient.

Each component is one difference across one cell, evaluated where it lives: (f[i] − f[i−1]) / (h_d · (x[i] − x[i−1])), with h_d the scale factor at the face. No averaging enters, so the C arrangement's gradient is second order on a uniform mesh and free of the two-grid null space a collocated difference has.

The outer two faces of a bounded direction have a cell on one side only, so no difference exists there and they are written masked. That is where a boundary condition goes: pass masked = 0 for a zero-gradient edge, or write the edge yourself. A wrapping direction has no such face, its faces all lying between two cells.

A point whose cells are not all active is masked too, on the same rule.

gradient(f, sg; kwargs...) -> NTuple{N,Array}

gradient! into fresh arrays, one per direction at that direction's face location.

gradient_plan(grid; stencil=Stencils.Axial(1), active_only=true, conn=nothing) -> GradientPlan

Build the least-squares gradient of grid — the geometry of it, with no field involved. See GradientPlan for what it is and why it is that; apply it with gradient!.

This is the counterpart of apply_stencil! for the two architectures that have no separable axis to difference along: a CurvilinearGrid, whose neighbours come from its index topology, and an UnstructuredGrid, whose come from its stored adjacency. conn overrides the neighbour set; otherwise one is built, from stencil where the architecture takes one.

One direction per coordinate: a (λ, φ) or (x, y) surface resolves two, in its tangent plane, and a three-coordinate volume resolves three, on the local frame — see Geometry.local_displacement. The fixed-size symmetric solve behind it covers those two cases.

A masked cell gets no coefficients at all and reads zero gradient, and an inactive neighbour is not offered to the fit, on the same rule as everywhere else: not determined by the active data, so not invented.

interpolate(field, grid, p; policy=BlankMasked(), masked=NaN, …) -> value

The value of field at the coordinate p, which is how observational data arrives: a station, a float or a ship track carries a coordinate.

Discretization.interpolation_weights gives the weights along one axis; this composes them, on every layout.

• StructuredGrid — multilinear, the tensor product of the per-axis weights. A periodic direction interpolates across its seam, so a coordinate past the last sample wraps to the first.
• CurvilinearGrid, UnstructuredGrid — a weighted least-squares plane fitted to the k nearest cells in the tangent plane at p, which is exact for a linear field and reproduces a cell's own value at its centre. Falls back to the weighted mean where the fit is rank deficient, that being the part of it the data still determines.

p may be written any way a point is accepted elsewhere.

The mask policies say what an inactive contributor means, as they do for a stencil: BlankMasked — the default — returns masked if any contributor is inactive, and ReduceInRun renormalizes over the active ones. ShiftWithinRun has no meaning here, there being no window to shift, and says so.

A field carrying trailing batch axes beyond the grid's own — many tracers, or an ensemble, sharing one geometry — is evaluated for every element in one call: see interpolate! for the form that writes into a caller's buffer, which this one wraps.

interpolate!(out, field, grid, p; …) -> out

interpolate for a batched field, writing one value per batch element into out.

The bracketing cell — or, off a rectilinear grid, the neighbour set and the least-squares fit — depends on the point and the geometry alone, so it is solved once and applied to every element. One call is therefore less work than interpolate per slice.

interpolate!(out, field, grid, p; k=8, …) -> out

interpolate off a rectilinear grid for a field carrying trailing batch axes, writing one value per batch element.

The k nearest cells and the mask verdict are a property of the point and the geometry — and the k-d tree query is the expensive part — so they are solved once here and the tangent-plane fit is then applied to every element.

interpolate!(out, field, grid, p; active_only=true, masked=NaN, policy=BlankMasked()) -> out

interpolate at one coordinate for a field carrying trailing batch axes: out receives one value per batch element, in the order those axes are laid out.

The bracketing cell and its 2^N corner weights depend on the point and the grid alone, so they are solved once here and applied to every element.

ncomponents(plan::GradientPlan{D}) -> D

How many directions the plan resolves, and so how many outputs gradient! writes.