Quickstart

Yonderdrake adds two things to Firedrake: markers you put inside an ordinary UFL residual, and a stepper that knows how to evolve them. Everything else is Firedrake.

A time-fractional ODE to start with

Solve \(D_C^{0.6}u+u=0\), \(u(0)=1\):

import firedrake as fd
from yonderdrake import (
    BirkSong,
    CaputoDerivative,
    Diethelm,
    FractionalTimeStepper,
    FullHistory,
)

mesh = fd.UnitIntervalMesh(1)
V = fd.FunctionSpace(mesh, "CG", 1)
u = fd.Function(V).assign(1.0)
v = fd.TestFunction(V)
t = fd.Constant(0.0)
dt = fd.Constant(0.01)

F = (fd.inner(CaputoDerivative(u, 0.6), v) + fd.inner(u, v)) * fd.dx
stepper = FractionalTimeStepper(F, BirkSong(48), t, dt, u)

for _ in range(100):
    stepper.advance()
    t.assign(t + dt)          # advance() updates u and memory, not t

print(f"u({float(t):.2f}) = {u.dat.data_ro[0]:.6f}")

The final value should be close to 0.41. The exact solution is the Mittag-Leffler value \(E_{0.6}(-1)\). Timestep and representation errors explain the remaining difference.

The quickstart computation over \(0\leq t\leq1\), compared with \(E_{0.6}(-t^{0.6})\). The lower panel shows the absolute error from the same 48-mode Birk-Song calculation with \(\Delta t=0.01\).

Three choices appear here:

Choice	Here	Alternatives
Which derivative	CaputoDerivative	RiemannLiouvilleDerivative
How memory is represented	BirkSong(48)	Diethelm, FullHistory
How memory is advanced	default Recurrence	AuxiliaryODE

Swapping a representation means changing one argument:

stepper = FractionalTimeStepper(F, Diethelm(48), t, dt, u)
stepper = FractionalTimeStepper(F, FullHistory(), t, dt, u)

Stick with BirkSong or Diethelm unless you are reproducing a specific paper. See use the defaults. For single-exponential memory rather than a fractional derivative, see Exponential memory and Caputo-Fabrizio.

Adding space

The marker must wrap the stepped field u directly. Fixed spatial operators go around it. Time-fractional diffusion is then just UFL:

F = (
    fd.inner(CaputoDerivative(u, alpha), v)
    + kappa * fd.inner(fd.grad(u), fd.grad(v))
    - fd.inner(source, v)
) * fd.dx

A localized field spreading across a Pac-Man-shaped domain under a Caputo time derivative and an ordinary spatial Laplacian. The field lives on a 2D domain, but only the time derivative is fractional in this example.

For a spatially fractional operator, use one of the three spatial operators as an ordinary term in the residual:

from yonderdrake import SpectralFractionalLaplacian

bc = fd.DirichletBC(V, 0.0, "on_boundary")
Lu = SpectralFractionalLaplacian(u, 0.4, bcs=bc)
F = (fd.inner(CaputoDerivative(u, alpha), v) + fd.inner(Lu, v)) * fd.dx

SpectralFractionalLaplacian, RieszFractionalLaplacian, and PeriodicFractionalLaplacian are different operators, not three approximations of one. Choose deliberately: which one you want.

Controlling numerical error

Fractional problems have more independent error sources than classical ones. Refine one at a time:

Symptom	Control	Note
Memory poorly resolved	mode count, rate_scale	independent of dt
Time discretization	dt	independent of mode count
Spatial discretization	mesh	refine last
Spectral quadrature	sinc_truncation_target	model estimate, not an FE bound
Riesz quadrature	quadrature_degree	separate from compression tolerance
Periodic Fourier	uniform grid size	keep represented modes below Nyquist
Algebraic	PETSc tolerances	keep below the discretization error

Separating error sources gives a separate refinement sequence for each operator.

Refining the quickstart

Rerun the first example with a finer memory representation and a smaller timestep:

u.assign(1.0)
t.assign(0.0)
dt.assign(0.001)
stepper = FractionalTimeStepper(F, BirkSong(256), t, dt, u)

for _ in range(1000):
    stepper.advance()
    t.assign(t + dt)

The original calculation uses 48 modes and \(\Delta t=0.01\). The refined calculation uses 256 modes and \(\Delta t=0.001\). Refining both independent controls reduces the error over the full interval.

The refinement costs more work. On our test MacBook Pro with an Apple M3 Pro, the original 100-step loop took about 0.05 seconds. The refined 1,000-step loop took about 1.10 seconds. These are medians of seven warm runs. This is the usual trade-off between accuracy and computational cost.

Common setup problems

Symptom	What to check
The marker is rejected	Wrap the stepped field u directly. Put fixed spatial operators around the marker.
The solution advances but displayed time does not	Update the caller-owned t only after each successful advance().
A periodic operator rejects the mesh	Use a fully periodic, uniform Q1 interval, quadrilateral rectangle, or hexahedral box with at least three cells per direction.
A Riesz operator rejects the space	Use scalar CG1 or CG2 on affine 2D triangles and provide complete homogeneous boundary conditions when \(s\geq1/2\).
Sinc construction emits a precision warning	The requested target is below useful float64 resolution. Choose a realistic truncation target and measure mesh error separately.
An imaging demo reaches its iteration cap	Inspect the reported optimizer reason and gradient norm, then raise --inverse-iterations or relax --inverse-tolerance deliberately.

Exact supported combinations and validation rules are collected in API and supported scope.

Next

• Guides and examples: worked examples, one feature at a time.

• Gallery: larger applications, with animations.

• Mathematics and methods: definitions, methods, and their papers.

• API and supported scope: signatures and exact limitations.