Automatic/algorithmic differentiation (AD)
SummationByPartsOperators.jl is written using generic Julia code whenever possible. This means that standard AD tools just work. For example, computing the Jacobian using ForwardDiff.jl is possible using the following code.
julia> using SummationByPartsOperators, ForwardDiffjulia> D = periodic_derivative_operator(derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 8)Periodic first-derivative operator of order 2 on a grid in [0.0, 1.0] using 8 nodes,stencils with 1 nodes to the left, 1 nodes to the right, and coefficients of Fornberg (1998) Calculation of Weights in Finite Difference Formulas. SIAM Rev. 40.3, pp. 685-691.julia> u = rand(size(D, 2));julia> J = ForwardDiff.jacobian(u -> D * u, u)8×8 Matrix{Float64}: 0.0 4.0 0.0 0.0 0.0 0.0 0.0 -4.0 -4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0 4.0 4.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0julia> J ≈ Matrix(D)trueThis works of course also for non-periodic SBP operators, e.g.,
julia> using SummationByPartsOperators, ForwardDiffjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 8)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] using 8 nodesand coefficients of Mattsson, Nordström (2004) Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics 199, pp. 503-540.julia> u = rand(size(D, 2));julia> J = ForwardDiff.jacobian(u -> D * u, u)8×8 Matrix{Float64}: -7.0 7.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.5 0.0 3.5 0.0 0.0 0.0 0.0 0.0 0.0 -3.5 0.0 3.5 0.0 0.0 0.0 0.0 0.0 0.0 -3.5 0.0 3.5 0.0 0.0 0.0 0.0 0.0 0.0 -3.5 0.0 3.5 0.0 0.0 0.0 0.0 0.0 0.0 -3.5 0.0 3.5 0.0 0.0 0.0 0.0 0.0 0.0 -3.5 0.0 3.5 0.0 0.0 0.0 0.0 0.0 0.0 -7.0 7.0julia> J ≈ Matrix(D)trueHowever, this does not work for Fourier derivative operators - and all other operators involving an FFT - since FFTW.jl cannot handle dual numbers and a simple reinterpret trick does also not help since FFTW.jl requires unit strides.
Jacobian-vector products
There is a neat trick that you can use if you are only interested in Jacobian-vector products. ForwardDiff.jl does not offer such a functionality at the time of writing, but Simon Byrne suggested the following implementation.
using SummationByPartsOperators, ForwardDiff, StructArraysfunction StructDual(x::AbstractVector{T}, w::AbstractVector{T}) where {T} @assert length(x) == length(w) # This was the original suggestion. However, it is currently not stable # under broadcasting. Thus, we use a slightly different version. # partials = StructArray{ForwardDiff.Partials{1, T}}( # (StructArray{Tuple{T}}( # (w,) # ),) # ) partials = reinterpret(reshape, ForwardDiff.Partials{1, T}, w) duals = StructArray{ForwardDiff.Dual{Nothing, T, 1}}((x, partials)) return dualsendfunction ForwardDiff.value(dx::StructArray{D}) where {D <: ForwardDiff.Dual} return dx.valueendfunction ForwardDiff.partials(dx::StructArray{<: ForwardDiff.Dual{Tag, T, 1}}, i) where {Tag, T} # This was the original suggestion. We need to update it (see above). # return getproperty(dx.partials.values, i) @assert i == 1 return reinterpret(reshape, T, dx.partials)endYou can use it as follows to compute the Jacobian-vector product
\[J_f(u) \cdot v\]
for the function $f$ given by
\[f(u) = D u.\]
D = fourier_derivative_operator(xmin = 0.0, xmax = 1.0, N = 8)u = randn(size(D, 2)); # the point `u` where we want to compute the derivativev = randn(size(D, 2)); # the direction `v` in which we want to compute the derivativeu_v = StructDual(u, v); # combined StructArray containing the value and directionf_df = D * u_v # compute the function value and its derivative@assert ForwardDiff.value(f_df) ≈ D * u@assert ForwardDiff.partials(f_df, 1) ≈ D * v # the Jacobian of `f(u) = D * u` is `D`You can of course also use this with nonlinear functions, e.g.,
f(u, D) = u .* (D * (u.^2))f_df = f(u_v, D)8-element StructArray(::Vector{Float64}, ::Vector{ForwardDiff.Partials{1, Float64}}) with eltype ForwardDiff.Dual{Nothing, Float64, 1}:
Dual{Nothing}(-44.50228965668829,36.069592079888885)
Dual{Nothing}(29.27320803803988,-21.185986662574216)
Dual{Nothing}(-1.1621101793030821,13.55104518741151)
Dual{Nothing}(-15.309927153749577,-13.789423204510712)
Dual{Nothing}(-28.124762639852136,1.1614449360505876)
Dual{Nothing}(11.657982157918777,22.243685846286176)
Dual{Nothing}(11.37583081093099,-29.208213204644206)
Dual{Nothing}(-14.08364373134745,-51.24794170650198)The Jacobian of this function is
using LinearAlgebraJ = Diagonal(D * u.^2) + 2 .* u .* Matrix(D) * Diagonal(u)@assert ForwardDiff.value(f_df) ≈ f(u, D)@assert ForwardDiff.partials(f_df, 1) ≈ J * vReproducibility
These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.status(["SummationByPartsOperators", "ForwardDiff", "StructArrays"], mode=PKGMODE_MANIFEST)Julia Version 1.10.11
Commit a2b11907d7b (2026-03-09 14:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × Intel(R) Xeon(R) 6973P-C
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, icelake-client)
Threads: 2 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl/docs/Manifest.toml`
[f6369f11] ForwardDiff v1.4.3
[09ab397b] StructArrays v0.7.3
[9f78cca6] SummationByPartsOperators v0.5.97-DEV `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`