Discretization

The geometric inputs a numerical method needs from a grid: where a point falls, the weights that interpolate to it, where a cell's faces are, the metric factors of the coordinate system, and the finite-difference weights of any stencil.

All of it is a function of coordinates alone. Applying weights to a field is not here — that needs a result location, a boundary-condition policy and a halo convention, which are the caller's to choose, not this package's to impose.

Centres and faces

using FlowGeometries: FlowGeometries as FG
D = FG.Discretization
A = FG.Axes

x = [0.0, 1.0, 3.0, 6.0]
D.faces(x)

5-element Vector{Float64}:
 -0.5
  0.5
  2.0
  4.5
  7.5

N centres give N+1 faces, at the midpoints, with the outermost two extrapolated a half-cell beyond the end centres. Note that centres do not determine faces — the system is underdetermined by one — so that midpoint rule is a stated convention, the same one the curvilinear corner reconstruction uses.

A uniform axis keeps its uniformity, so nothing is lost by asking for its faces:

u = A.UniformAxis(0.0, 0.5, 6)
A.isuniform(D.faces(u)), A.spacing(D.faces(u)), length(D.faces(u))

(true, 0.5, 7)

D.nodes(u, D.Center()) === u, length(D.nodes(u, D.Face()))

(true, 7)

Gaps and cell widths

local_spacing is the gap either side of one sample, and cell_width the width of one cell. Both are a scalar subtraction of two stored numbers — no array is built — so they are the forms to call per grid point, where faces would materialize the whole axis to answer about one cell.

xs = cumsum([0.0, 1.0, 0.3, 2.5, 0.7, 4.0])
D.local_spacing(xs, 3), D.cell_width(xs, 3)

((0.30000000000000004, 2.5), 1.4)

cell_width is exactly the distance between the cell's two faces — that is what it means:

f = D.faces(xs)
[D.cell_width(xs, i) for i in eachindex(xs)] ≈ abs.(diff(f))

true

The gaps are signed and the width is not. A derivative needs the sign, since it distinguishes an axis stored increasing from one stored decreasing; a width is a length and cannot be negative. Reverse the axis and cell i moves to n+1-i, where the gaps come back negated and swapped — the same two neighbours, now on the other side — while the width is unchanged:

xr = reverse(xs)                              # cell 3 of xs is cell 4 of xr
D.local_spacing(xr, 4), D.cell_width(xr, 4)

((-2.5, -0.30000000000000004), 1.4)

At a bounded end the outward gap is 0 and the caller falls back to a one-sided stencil. Given a period it wraps instead, which is what makes a seam no different from the interior:

λ = collect(range(0, 2π * (1 - 1/8); length = 8))
D.local_spacing(λ, 8), D.local_spacing(λ, 8, 2π)

((0.7853981633974483, 0.0), (0.7853981633974483, 0.7853981633974483))

cell_widths is the whole axis at once. A uniform axis returns an Axes.ConstantVector — one number and a length — so nothing is materialized for it either:

D.cell_widths(u), D.cell_widths(xs)

(ConstantVector{Float64}(0.5, 6), [1.0, 0.65, 1.4, 1.6, 2.35, 4.0])

On a grid, Grids.local_spacing, Grids.cell_width and Grids.cell_widths take the direction's axis and its wrap period from the grid, so a periodic seam is right without being asked for — see the Grids page.

Point location

D.locate(u, 0.2), D.locate(u, 0.4), D.locate(u, -5.0)

(1, 2, 0)

locate returns the cell containing a coordinate, or 0 outside it. Cells are the intervals between faces, so cell 1 of u spans [-0.25, 0.25). It is O(1) on a uniform axis — a direct payoff of the spacing living in the type — and O(log n) by bisection on a stretched one. Both storage orders work.

D.locate(xs, 2.0), D.nearest_index(xs, 2.0)

(3, 3)

nearest_index always returns a valid index, clamping rather than reporting "outside"; exact ties go to the lower index. It brackets and compares rather than scanning, so it is O(log n) on a stretched axis and O(1) on a uniform one, and both paths compare the same two samples — a uniform axis and its collect hold the same numbers, so they must not disagree about which is nearest.

Interpolation weights

i, w = D.interpolation_weights(xs, 2.0)
i, w, sum(w)

(3, (0.72, 0.27999999999999997), 1.0)

Weights only — applying them is w[1]*f[i] + w[2]*f[i+1] at the call site. For higher order on an arbitrarily spaced axis, lagrange_weights is exact for polynomials up to nodes-1:

idx, wl = D.lagrange_weights(xs, 2.0, 4)
idx, sum(wl)

(2:5, 1.0)

Finite-difference weights

fd_weights is the recursion of Fornberg (1988), Math. Comp. 51, 699–706. One recursion covers every case: any derivative order, any node count (hence any order of accuracy), any evaluation point, and arbitrarily spaced nodes.

D.fd_weights([-1.0, 0.0, 1.0], 0.0, 1)      # centred first difference

3-element Vector{Float64}:
 -0.5
  0.0
  0.5

D.fd_weights([-1.0, 0.0, 1.0], 0.0, 2)      # centred second difference

3-element Vector{Float64}:
  1.0
 -2.0
  1.0

D.fd_weights([-2.0, -1.0, 0.0, 1.0, 2.0], 0.0, 1)   # fourth-order first derivative

5-element Vector{Float64}:
  0.08333333333333333
 -0.6666666666666666
  0.0
  0.6666666666666666
 -0.08333333333333333

With m nodes the weights are exact for polynomials of degree m-1, so the accuracy order is m - order. Nothing about that requires equal spacing:

D.fd_weights([-1.7, -0.4, 0.9], 0.0, 1)     # same order, unequal nodes

3-element Vector{Float64}:
 -0.14792899408284027
 -0.4733727810650888
  0.621301775147929

The axis form centres the stencil on a sample and shifts it inward at a boundary, so the node count — and therefore the accuracy order — is the same everywhere, rather than degrading at the two ends:

xg = collect(range(0.0, 1.0; length = 11))
idx1, w1 = D.fd_weights(xg, 1, 1, 5)        # at the first sample
idxm, wm = D.fd_weights(xg, 6, 1, 5)        # in the interior
idx1, idxm

(1:5, 4:8)

FG.Geometry.nonuniform_first_derivative is the three-node, first-derivative case of the same thing. Fed the gaps above it is a complete derivative at a point, which is the shape an operator assembled at the call site takes — a divergence, a curl, a staggered difference are each this plus the metric factors and whatever boundary policy the caller chose:

g(x) = 3x^2 - 2x + 5                        # exact for a quadratic, on any spacing
maximum(abs(FG.Geometry.nonuniform_first_derivative(
                g(xs[i-1]), g(xs[i]), g(xs[i+1]), D.local_spacing(xs, i)...) - (6xs[i] - 2))
        for i in 2:length(xs)-1)

2.6645352591003757e-15

The gaps are signed for exactly this reason: the same expression on a descending axis differentiates with respect to the coordinate, not the index, without the caller correcting anything:

maximum(abs(FG.Geometry.nonuniform_first_derivative(
                g(xr[i-1]), g(xr[i]), g(xr[i+1]), D.local_spacing(xr, i)...) - (6xr[i] - 2))
        for i in 2:length(xr)-1)

3.552713678800501e-15

Applying a weight set

apply_stencil! is the one function here that touches a field, and only along a single direction with the result left where the input was. That case needs no convention the package has not already fixed: nothing to stagger, fd_weights' inward shift at a bounded end, wrapping on a periodic one — so no halo either.

x = collect(range(0.0, 2.0; length = 11))
f = @. 3x^2 - 2x + 5
out = similar(f)
D.apply_stencil!(out, f, x, 1; order = 1, nodes = 3)
maximum(abs, out .- (6x .- 2))          # exact for a quadratic, ends included

7.105427357601002e-15

A stretched axis is equally exact, because the weights are built per sample rather than one set reused:

xs = [0.0, 0.11, 0.37, 0.9, 1.05, 1.6, 1.62, 2.0]
outs = similar(xs)
D.apply_stencil!(outs, (@. 3xs^2 - 2xs + 5), xs, 1; order = 1, nodes = 3)
maximum(abs, outs .- (6xs .- 2))

1.3677947663381929e-13

Given a period the stencil stays centred and wraps, carrying the wrapped samples' coordinates across the seam so the spacing there is the true one — the seam is then no worse than the interior:

λ = collect(range(0, 2π; length = 65)[1:64])
o = similar(λ)
D.apply_stencil!(o, sin.(λ), λ, 1; order = 1, nodes = 5, period = 2π)
maximum(abs, o .- cos.(λ)), abs(o[1] - cos(λ[1]))

(3.0930005776586e-6, 3.093000576104288e-6)

The grid form takes the axis, the wrap period and the mask from the grid, so a periodic direction wraps without being told. Where a mask bites, a value whose stencil would read an inactive cell is written as masked rather than invented:

geo = FG.Geometry.CartesianGeometry()
X = collect(range(0.0, 1.0; length = 9)); Y = collect(range(0.0, 2.0; length = 7))
F = [xi^2 + 3yi for xi in X, yi in Y]
mk = trues(9, 7); mk[5, 3] = false
gm = FG.Grids.StructuredGrid(geo, X, Y, mk)
Om = similar(F)
D.apply_stencil!(Om, F, gm, 1; order = 1, nodes = 3, masked = NaN)
Om[3:7, 3]                               # the masked cell and the two that read it

5-element Vector{Float64}:
   0.5
 NaN
 NaN
 NaN
   1.5

What a mask edge does to the stencil

Blanking is the default, and it is not free: a cell is blanked when any sample its window reads is inactive, so a single masked cell takes out up to nodes - 1 cells either side of it. A five-point derivative on a short axis can be annihilated outright by one hole.

The window already shifts inward at the end of an axis, precisely so the node count — and the accuracy order — survives a boundary. The end of an active run is the same situation, and ShiftWithinRun treats it that way:

xr = collect(0.0:1.0:6.0)
mr = trues(7, 1); mr[4, 1] = false
gr = FG.Grids.StructuredGrid(geo, xr, [0.0], mr)
fr = reshape(collect(0.0:6.0), 7, 1)         # f = x, so df/dx is exactly 1

blanked  = zeros(7, 1); shifted = zeros(7, 1)
D.apply_stencil!(blanked, fr, gr, 1; order = 1, nodes = 5, masked = NaN)
D.apply_stencil!(shifted, fr, gr, 1; order = 1, nodes = 5, masked = NaN,
                 policy = D.ShiftWithinRun())
vec(blanked), vec(shifted)

([NaN, NaN, NaN, NaN, NaN, NaN, NaN], [NaN, NaN, NaN, NaN, NaN, NaN, NaN])

Every cell is blanked in the first, and every active cell is exact in the second. The policies are:

ReduceInRun is the only one that will silently lower the accuracy order, which is why it is a separate policy rather than a fallback inside ShiftWithinRun: a narrow strait genuinely wants an answer at reduced order, and everywhere else wants to be told the run is too short.

A cell whose window already lies inside its run reuses the precomputed row, so a run's interior is bit-for-bit what BlankMasked gives, and only the cells within nodes - 1 of an edge cost anything extra. A run that wraps a periodic seam is one run, not two.

Build the weights once with axis_stencils to reuse them across many fields; applying a precomputed set allocates nothing.

idx, w = D.axis_stencils(X, 1, 3)
size(idx), size(w)

((9, 3), (9, 3))

The physical derivative

apply_stencil! differentiates with respect to a coordinate. On any curved geometry that is not the derivative a physical law is written in — that one is per unit distance, ∂f/∂sᵈ = (1/hᵈ)·∂f/∂ξᵈ. derivative! is the two together:

R   = 6.371e6
sph = FG.Geometry.SphericalGeometry(R)
lon = collect(range(0, 2π * (1 - 1/48); length = 48))
lat = collect(range(-π/2, π/2; length = 25))          # both poles are rows of this grid
gs  = FG.Grids.StructuredGrid(sph, lon, lat)
fs  = [sin(φ) for _ in lon, φ in lat]                 # ∂/∂north should be cos(φ)/R
dn  = zeros(48, 25)
D.derivative!(dn, fs, gs, 2; order = 1, nodes = 5, masked = NaN)
dn[1, 12], cos(lat[12]) / R

(1.5561688561297248e-7, 1.5561840548953231e-7)

Where the metric degenerates the derivative does not exist, and masked is written rather than a number invented. Longitude at a pole is that case — h_λ = R cos φ → 0:

de = zeros(48, 25)
D.derivative!(de, [sin(λ) * cos(φ) for λ in lon, φ in lat], gs, 1;
              order = 1, nodes = 5, masked = NaN)
all(isnan, de[:, 1]), all(isnan, de[:, 25]), any(isnan, de[:, 2:24])

(true, true, false)

The threshold is metric_floor, L·√eps(T) — relative to the geometry's size and to the element type. An absolute constant cannot serve both: 1e-12 is below eps(Float32), and in Float32 cos(Float32(π/2)) ≈ -4.4e-8, so h_λ at the pole is around 0.28 metres and a fixed small threshold never fires.

A divergence or a curl remains the caller's to assemble, needing a result location and a boundary-condition policy this does not choose. When you do, note the flux form — on a sphere

\[\nabla\cdot\mathbf{u} = \frac{1}{R\cos\varphi}\left[\frac{\partial u_\lambda}{\partial\lambda} + \frac{\partial (u_\varphi\cos\varphi)}{\partial\varphi}\right]\]

so the second term differentiates u_φ·cos φ, not u_φ. Adding two physical derivatives is a different expression, and a wrong one.

Off a rectilinear grid: the least-squares gradient

apply_stencil! needs a separable axis to difference along. A CurvilinearGrid has none, and a node set's neighbours come from connectivity rather than an index offset, so neither has a stencil. Connectivity.gradient_plan builds a least-squares gradient for them instead, and gradient! applies it:

cart = FG.Geometry.CartesianGeometry{Float64}()
nn   = 14
xg   = [t + 0.35u for t in range(0, 10; length = nn), u in range(0, 6; length = nn)]   # sheared
yg   = [u - 0.2t  for t in range(0, 10; length = nn), u in range(0, 6; length = nn)]
cgrid = FG.Grids.CurvilinearGrid(cart, xg, yg, trues(nn, nn); measure = fill(1.0, nn, nn))
plan  = FG.Connectivity.gradient_plan(cgrid)

fld = 2.0 .* xg .- 3.0 .* yg .+ 7        # ∇ = (2, -3) everywhere
g1, g2 = zeros(nn, nn), zeros(nn, nn)
D.gradient!(g1, g2, fld, plan)
maximum(abs, g1 .- 2), maximum(abs, g2 .+ 3)

(5.551115123125783e-15, 5.329070518200751e-15)

Exact for a linear field on any stencil, however skewed — the least-squares combination cancels the leading truncation term, which inverting a 2×2 index-space Jacobian does not. It is second order on locally symmetric stencils and degrades toward first on strongly skewed cells, and where the stencil is separable and orthogonal it is the centred difference, so it agrees with apply_stencil! where both apply.

A depends only on the geometry, so the plan holds the per-neighbour coefficients and each apply is one dot product per cell, allocating nothing. Where A is rank deficient — every neighbour on one line, at a boundary or beside a mask — that component is zeroed rather than invented, the same rule apply_stencil! states for a mask.

Evaluating a field at a coordinate

Observational data has a coordinate, not a cell index. interpolate answers for one:

xa = collect(range(0, 2; length = 11)); ya = collect(range(-1, 3; length = 9))
ga = FG.Grids.StructuredGrid(cart, xa, ya)
fa = [2xi - 3yj + 0.5xi * yj + 7 for xi in xa, yj in ya]
D.interpolate(fa, ga, (0.7, 1.1)), 2*0.7 - 3*1.1 + 0.5*0.7*1.1 + 7

(5.484999999999999, 5.484999999999999)

Multilinear on a rectilinear grid; a weighted least-squares plane on a curvilinear grid or a node set, which is exact for a linear field and reproduces a cell's own value at its centre. A periodic direction interpolates across its seam rather than clamping at the last sample — composing the per-axis weights by hand gets that wrong by half a cell everywhere on the seam.

The mask policies mean here what they mean for a stencil: BlankMasked returns masked when a contributor is inactive, ReduceInRun renormalizes over the active ones, and ShiftWithinRun is refused, there being no window to shift.

Anything that does need a convention the package has not chosen — a staggered difference, or a multi-direction operator like a divergence or a curl, which also need a result location and a boundary-condition policy — is assembled at the call site from these weights and the metric factors below.

Metric factors

G = FG.Geometry
sph = G.SphericalGeometry(2.0)
G.scale_factors(sph, (0.0, π/3)), G.scale_factors(sph, (0.0, π/3, 5.0))

((1.0000000000000002, 2.0), (2.5000000000000004, 5.0, 1.0))

scale_factors gives the physical length of a unit coordinate step in each direction — (R cosφ, R) on a sphere's surface, (r cosφ, r, 1) with a radius direction, and 1 throughout for a Cartesian metric. That is what turns a coordinate derivative into a physical one, ∂/∂sᵈ = (1/hᵈ)·∂/∂ξᵈ, so a divergence or a curl is assembled from these plus fd_weights with its conventions stated at the call site.

G.jacobian(sph, (0.0, 0.5)), 4cos(0.5)

(3.510330247561491, 3.510330247561491)

jacobian is their product: the volume element per unit coordinate volume, and the quantity the grid's own cell measure is built from.

An ellipsoidal metric supplies its own curvature radii, and nothing else in the stack changes:

wgs = G.SpheroidGeometry()
G.scale_factors(wgs, (0.0, 0.0)), G.meridional_radius(wgs, π/2) > wgs.a

((6.378137e6, 6.3354393272928195e6), true)