GitHub

Guide

This guide builds the StencilCore types by hand. In normal use you would obtain a stencil from StencilCalculus (by differentiation) and hand it to StencilAssembly (to assemble a matrix); here we look under the hood.

Lattice offsets

Offsets are type-level: a StaticShift is a normalized collection of StaticPair{D,O} (axis D, offset O). They support an algebra and print like lattice vectors, with ô the zero shift and ê₁ … ê₉ the unit shifts:

using StencilCore

-2ê₁              # StaticShift{Tuple{StaticPair{1,-2}}}
3ê₁ + ê₂          # StaticShift{Tuple{StaticPair{1,3}, StaticPair{2,1}}}
ê₁ - ê₁           # ô  (the zero shift)

Same-axis pairs are summed, zero offsets dropped, and the result is sorted by axis — so the representation is canonical.

The scalar algebra

A scalar — a timestep, a Reynolds number, a literal 2.0 — has no spatial extent. StencilCore models it directly with AbstractScalar{T} and a small CAS that mirrors the term-side one. Five concrete types:

  • Var{S,T} — a named substitution leaf. S is a Symbol (the name); T is the materialised element type. Declared with @var name T.
  • Constant{T} — a literal leaf carrying its val::T. T is any concrete type — Number, SVector, anything. This is what numeric (and non-Number) literals canonicalise to at the operator boundary.
  • Scalar{F,A,T} — an interior tree node fn(args…).
  • Null{T} / Unity{T} — type-level structural 0 / 1. Unity requires one(T) defined (Number, square SMatrix); its outer ctor routes an SVector through to its square Jacobian space.

The operator overloads build Scalar trees, with numeric literals canonicalising to Constant:

@var τ Float64
α = Constant(2.0)

τ * α + α              # Scalar(+, (Scalar(*, (τ, α)), α))
typeof(τ * α + α)      # Scalar{typeof(+), …, Float64}
eltype(τ * α + α)      # Float64

simplify post-walks a scalar tree with the same rule-rewriter shape as the term side. Identities are purely structural — Null + x → x, x * Unity → x (with an eltype-preservation gate), x * Null → Null — never inspecting .val. Folding combines values; it does not match on them.

simplify(Null{Float64}() + τ)              # τ
simplify(τ * Unity{Float64}())             # τ
simplify(Constant(2.0) + Constant(3.0))    # Constant(5.0)
simplify(τ * Constant(1.0))                # stays a Scalar — numerical `1`
                                           # is not a structural identity

differentiate on a scalar tree returns an AbstractScalar whose eltype is the Jacobian of eltype(s) w.r.t. eltype(v). For Number variables the Jacobian is a Number; for SVector{N} variables it is the square SMatrix{N, N}. Mixed shape-classes are rejected at the top level.

@var τ Float64
differentiate(sin(τ) * τ, τ)               # cos(τ)*τ + sin(τ)*1 simplified
differentiate(Constant(2.0), τ)            # Null{Float64}() — no dependence

# Vector-valued: ∂(2x)/∂x lands in the matching SMatrix Jacobian space.
@var x SVector{2, Float64}
differentiate(2x, x)

materialize reduces a scalar tree to a single value, substituting Var leaves from a NamedTuple:

materialize(τ * Constant(3.0), (τ = 4.0,))  # 12.0
materialize(Unity{SMatrix{2,2,Float64,4}}())  # the 2×2 identity matrix

These scalars become term coefficients once StencilCalculus wraps them in a Fill — see its docs for the bridge.

Generating methods from scalar expressions

@addmethod installs a typed Julia method whose body is derived from a scalar expression. Each Var argument declares both the argument name (from its S type parameter) and its declared type (from T):

using StencilCore

@var x Float64
@var y Int

foo = x ^ y
@addmethod myfun foo (y, x)
# equivalent to:
#   function myfun(y::Int, x::Float64)
#       x ^ y
#   end

myfun(2, 3.0)    # 9.0

Because @addmethod adds an ordinary Julia method, subsequent calls add more dispatch methods to the same function — hence the name:

using StaticArrays

@var A SMatrix{2, 2, Float64, 4}
@var u SVector{2, Float64}

bar = A * u
@addmethod myfun bar (A, u)
# adds:  function myfun(A::SMatrix{2,2,Float64,4}, u::SVector{2,Float64}); A * u; end

# Both methods coexist:
myfun(2, 3.0)                  # 9.0  (first method)
myfun(SMatrix{2,2}(1,0,0,2), SVector(1.0, 3.0))  # [1.0, 6.0]  (second method)

The argument list is a tuple (@addmethod fname expr (v1, v2, ...)), and the order determines the method signature. Every listed argument must be a Var; @addmethod raises an ArgumentError at call time if:

  • any Var in the expression is not listed as an argument (free variable),
  • an element in the argument list is not a Var,
  • the same name appears more than once, or
  • the expression calls fname itself — the generated method would recurse infinitely (this is why, e.g., getindex cannot be defined by indexing the same argument; the body c[i] lowers to getindex(c, i)).

Widening the emitted signature

The expression is built and validated at the Vars' concrete types, which keeps simplify and differentiate fully type-inferred. To emit a method that accepts a wider family of inputs, pass a NamedTuple of type overrides instead of the Var list — only the generated signature widens, the body is unchanged:

@var x Float64
@var y Int

foo = x ^ y
@addmethod myfun foo (y = Integer, x = AbstractFloat)
# adds:  function myfun(y::Integer, x::AbstractFloat); x ^ y; end

myfun(big(3), 2.0f0)   # works for any Integer / AbstractFloat

Each override must be a supertype of the Var's declared type (widening only). Note that StaticArrays element parameters are invariant (SVector{2,Float64} <: SVector{2,Number} is false), so shape widening uses the UnionAll form:

@var V SVector{2, Float64}
vv = V + V
@addmethod vfun vv (V = SVector{2, <:Number},)
vfun(SVector(1, 2))   # [2, 4]

Linear and star stencils

A LinearStencil has contiguous offsets along one mesh axis. Its coefficient is one SVector per cell, holding every per-offset coefficient on that column in ascending-offset order:

using StencilCore, StaticArrays
n = 5
# offsets 0,1 along axis 1; coefficients (-1, 1) at every column.
fwd = LinearStencil{1}(SUnitRange(0, 1), fill(SVector(-1.0, 1.0), n))

A StarStencil is the N-D star with symmetric reach −L … +L per axis, stored interlaced: a single SVector{2NL+1} per cell, in reverse-lex offset order, with the diagonal as the explicit middle slot — so a free diagonal term (Helmholtz, an unsteady term) has a home.

n1, n2 = 5, 4
# 2-D five-point Laplacian (L = 1 ⇒ M = 5):
# (axis2,−1), (axis1,−1), diagonal, (axis1,+1), (axis2,+1).
lap = StarStencil{1}(fill(SVector(-1.0, -1.0, 4.0, -1.0, -1.0), n1, n2))

The general stencil and narrowing

Stencil is the lingua franca: an explicit reverse-lex list of StaticShift offsets plus a matching SVector coefficient. It is not assembled directly — as_linear / as_star narrow it to an assemblable type by matching the offset pattern and reusing the coefficient verbatim:

st = Stencil(ColumnAccess, (-2ê₁, -ê₁, ô), fill(SVector(1.0, -4.0, 3.0), n))
ln = as_linear(st)        # LinearStencil{1, -2, 3, …}
ln.term === st.term       # true — verbatim copy

A Stencil whose offsets form the canonical star pattern narrows with as_star instead; mismatched patterns raise an ArgumentError.