Building PDE Operators
The @operator macro lets you write PDE operators in mathematical notation. It translates symbolic expressions into composable operator objects, producing a single weight matrix that applies the full PDE operator in one matrix-vector multiply.
This is the recommended way to build operators the built-ins don't cover. If your operator can't be expressed here, see Custom Operators for the closure escape hatch.
using RadialBasisFunctions
using StaticArrays
x = rand(SVector{2,Float64}, 100)
f(p) = sin(p[1]) * cos(p[2])
u = f.(x)Recognized symbols
| Symbol | Meaning |
|---|---|
∇², Δ | Laplacian |
∂(dim) | First partial derivative in dimension dim |
∂²(dim) | Second partial derivative in dimension dim |
∇ ⋅ (κ * ∇) | Diffusion operator (scalar or vector κ) |
c ⋅ ∇ | Advection operator (vector c) |
f, I | Identity operator |
| Everything else | Scalar coefficient |
Standard arithmetic (+, -, *) and unary negation work as expected. Scalars can be literals, variables, or expressions like k^2 or c[1].
The macro produces rank-0 operators — see Understanding Rank (N) for what that means, and Custom Operators if you need rank 1.
Helmholtz Operator
The Helmholtz equation
k² = 4.0
op = @operator ∇² + k² * f
helm_op = op(x)
# Verify against separate built-in operators
expected = laplacian(x)(u) .+ k² .* u
maximum(abs, helm_op(u) .- expected)2.0976553827267708e-12Diffusion — Textbook Notation
The diffusion operator @operator macro recognizes the textbook form directly:
κ = [2.0, 0.5]
op = @operator ∇ ⋅ (κ * ∇)
diff_op = op(x)
# Verify against separate built-in operators
expected = κ[1] .* partial(x, 2, 1)(u) .+ κ[2] .* partial(x, 2, 2)(u)
maximum(abs, diff_op(u) .- expected)1.2363443602225743e-12Scalar
op = @operator ∇ ⋅ (3.0 * ∇)
diff_iso = op(x) # equivalent to 3∇²f
expected = 3.0 .* laplacian(x)(u)
maximum(abs, diff_iso(u) .- expected)4.035882739117369e-12Anisotropic Diffusion — Explicit Partials
The same anisotropic diffusion can also be written with explicit per-dimension coefficients:
κ_x = 2.0
κ_y = 0.5
op = @operator κ_x * ∂²(1) + κ_y * ∂²(2)
aniso_op = op(x)
# Verify against separate built-in operators
expected = κ_x .* partial(x, 2, 1)(u) .+ κ_y .* partial(x, 2, 2)(u)
maximum(abs, aniso_op(u) .- expected)1.2363443602225743e-12When
Advection-Diffusion
The steady advection-diffusion equation
ν = 0.01
c = SVector(1.0, 0.5)
op = @operator ν * ∇² - c ⋅ ∇
advdiff_op = op(x)
# Verify against separate built-in operators
expected = ν .* laplacian(x)(u) .- c[1] .* partial(x, 1, 1)(u) .- c[2] .* partial(x, 1, 2)(u)
maximum(abs, advdiff_op(u) .- expected)4.6629367034256575e-14Sharing Stencils
When multiple operators act on the same point set, precompute the neighbor list once and pass it to avoid redundant nearest-neighbor searches:
adjl = find_neighbors(x, 30)
helm_op = (@operator ∇² + k² * f)(x; adjl=adjl)
aniso_op = (@operator κ_x * ∂²(1) + κ_y * ∂²(2))(x; adjl=adjl)RadialBasisOperator
├─Operator: 2.0 × ∂ⁿf/∂xᵢ (n = 2, i = 1) + 0.5 × ∂ⁿf/∂xᵢ (n = 2, i = 2)
├─Data type: StaticArraysCore.SVector{2, Float64}
├─Number of points: 100
├─Stencil size: 30
└─Basis: Polyharmonic spline (r³) with degree 2 polynomial augmentationBoundary Conditions on PDE Operators
Operators built with @operator accept the hermite keyword just like the built-ins — see Boundary Conditions for how to enforce Dirichlet, Neumann, or Robin conditions.