Time-fractional PDEs

A fractional time derivative composes with ordinary Firedrake spatial forms. The marker stays on the stepped field. The spatial operator stays outside it.

Time-fractional diffusion

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

time_fractional_diffusion.py uses the ordinary spatial Laplacian with the Diethelm representation. Refine modes at a fixed small dt, then dt at a fixed large mode count, then the mesh.

Caputo-Wismer-Kelvin damping

The fractionally damped acoustic wave equation is

\[ u_{tt}-c^2\Delta u-b\Delta(D_C^\alpha u)=0, \qquad 0<\alpha<1 . \]

For homogeneous Dirichlet data, integrate both Laplacians by parts and place the marker directly on the stepped field:

Du = CaputoDerivative(u, alpha)
F = (
    inner(u_tt, v)
    + c**2 * inner(grad(u), grad(v))
    + b * inner(grad(Du), grad(v))
) * dx
stepper = FractionalTimeStepper(F, Diethelm(32), t, dt, u, bcs=bc)

The last term is \(-b\Delta(D_C^\alpha u)\), a fractional time derivative of a classical Laplacian, not a fractional spatial Laplacian. It is well defined here because on a fixed domain, for a sufficiently regular field and a time-independent linear spatial operator,

\[ \Delta(D_C^\alpha u)=D_C^\alpha(\Delta u). \]

Do not wrap a strong CG Laplacian in the marker. Keep the marker on the evolving field. RiemannLiouvilleDerivative uses the same weak form and adds its exact trace term.

Several materials at once

With fixed indicator functions \(\chi_m\), one residual may carry several orders, each with its own representation:

F_layered = sum(
    b_m * chi_m * inner(grad(CaputoDerivative(u, alpha_m)), grad(v)) * dx
    for b_m, chi_m, alpha_m in materials
)

All terms share the field and the time increment. The BrainWeb and skullball demos in Gallery run this form under MPI. For reusable 2D and 3D wave and sensor-array helpers, see Caputo-Wismer imaging.