PDE Operators Cookbook

Recipes for assembling PDE-specific differential operators via @operator, producing a single weight matrix that applies the full PDE operator in one matrix-vector multiply.

using RadialBasisFunctions
using StaticArrays

x = [SVector{2}(rand(2)) for _ in 1:100]
f(p) = sin(p[1]) * cos(p[2])
u = f.(x)

Helmholtz Operator

The Helmholtz equation    appears in acoustics, electromagnetics, and quantum mechanics. The operator combines a Laplacian with a scaled identity.

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)

1.234568003383174e-12

Diffusion — Textbook Notation

The diffusion operator   appears in heat conduction, mass transfer, and many other physical models. The @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.807887173299605e-12

Scalar produces an isotropic operator (scaled Laplacian):

op = @operator ∇ ⋅ (3.0 * ∇)
diff_iso = op(x)  # equivalent to 3∇²f
expected = 3.0 .* laplacian(x)(u)
maximum(abs, diff_iso(u) .- expected)

3.268496584496461e-12

Anisotropic 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.807887173299605e-12

When  , this reduces to a scaled Laplacian.

Advection-Diffusion

The steady advection-diffusion equation     balances viscous diffusion against transport by a velocity field. It appears in fluid dynamics, pollutant transport, and thermal convection.

ν = 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)

1.4099832412739488e-14

Sharing 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: Custom Operator
├─Data type: StaticArraysCore.SVector{2, Float64}
├─Number of points: 100
├─Stencil size: 30
└─Basis: Polyharmonic spline (r³) with degree 2 polynomial augmentation