1  Introduction to Devito

This chapter introduces Devito, a domain-specific language (DSL) for solving partial differential equations using finite differences. We begin with the motivation for symbolic PDE specification, then work through a complete example using the 1D wave equation.

1.1 What is Devito?

Devito is a Python-based domain-specific language (DSL) for expressing and solving partial differential equations using finite difference methods. Rather than writing low-level loops that update arrays at each time step, you write the mathematical equations symbolically and let Devito generate optimized code automatically.

1.1.1 The Traditional Approach

Consider solving the 1D diffusion equation: \[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]

A traditional NumPy implementation might look like:

import numpy as np

# Parameters
Nx, Nt = 100, 1000
dx, dt = 0.01, 0.0001
alpha = 1.0
F = alpha * dt / dx**2  # Fourier number

# Initialize
u = np.zeros(Nx + 1)
u_new = np.zeros(Nx + 1)
u[Nx//2] = 1.0  # Initial impulse

# Time stepping loop
for n in range(Nt):
    for i in range(1, Nx):
        u_new[i] = u[i] + F * (u[i+1] - 2*u[i] + u[i-1])
    u, u_new = u_new, u  # Swap arrays

This approach has several limitations:

  1. Error-prone: Manual index arithmetic is easy to get wrong
  2. Hard to optimize: Achieving good performance requires expertise in vectorization, parallelization, and cache optimization
  3. Dimension-specific: The code must be rewritten for 2D or 3D problems
  4. Not portable: Optimizations for one architecture don’t transfer to others

1.1.2 The Devito Approach

With Devito, the same problem becomes:

from devito import Grid, TimeFunction, Eq, Operator, solve, Constant

# Problem parameters
Nx = 100
L = 1.0
alpha = 1.0  # diffusion coefficient
F = 0.5      # Fourier number (for stability, F <= 0.5)

# Compute dt from stability condition: F = alpha * dt / dx^2
dx = L / Nx
dt = F * dx**2 / alpha

# Create computational grid
grid = Grid(shape=(Nx + 1,), extent=(L,))

# Define the unknown field
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)

# Set initial condition
u.data[0, Nx // 2] = 1.0

# Define the PDE symbolically and solve for u.forward
a = Constant(name='a')
pde = u.dt - a * u.dx2
update = Eq(u.forward, solve(pde, u.forward))

# Create and run the operator
op = Operator([update])
op(time=1000, dt=dt, a=alpha)

This approach offers significant advantages:

  1. Mathematical clarity: The PDE u.dt - a * u.dx2 = 0 is written symbolically, and Devito derives the update formula automatically using solve()
  2. Automatic optimization: Devito generates C code with loop tiling, SIMD vectorization, and OpenMP parallelization
  3. Dimension-agnostic: The same code structure works for 1D, 2D, or 3D
  4. Portable performance: Generated code adapts to the target architecture

1.1.3 How Devito Works

Devito’s workflow consists of three stages:

Python DSL → Symbolic Processing → C Code Generation → Compilation → Execution

  1. Symbolic representation: Your Python code creates SymPy expressions that represent the PDE and its discretization

  2. Code generation: Devito analyzes the expressions and generates optimized C code with appropriate loop structures

  3. Just-in-time compilation: The C code is compiled (and cached) the first time the operator runs

  4. Execution: Subsequent runs use the cached compiled code for maximum performance

1.1.4 When to Use Devito

Devito excels at:

  • Explicit time-stepping schemes: Forward Euler, leapfrog, Runge-Kutta
  • Structured grids: Regular Cartesian meshes in 1D, 2D, or 3D
  • Stencil computations: Any PDE discretized with finite differences
  • Large-scale problems: Where performance optimization matters

Common applications include:

  • Wave propagation (acoustic, elastic, electromagnetic)
  • Heat conduction and diffusion
  • Computational fluid dynamics
  • Seismic imaging (reverse time migration, full waveform inversion)

1.1.5 Installation

Devito can be installed via pip:

pip install devito

For this book, we recommend installing the optional dependencies as well:

pip install devito[extras]

This includes visualization tools and additional solvers that we’ll use in later chapters.

1.1.6 What You’ll Learn

In this chapter, you will:

  1. Solve your first PDE (the 1D wave equation) using Devito
  2. Understand the core abstractions: Grid, Function, TimeFunction, Eq, and Operator
  3. Implement boundary conditions in Devito
  4. Verify your numerical solutions using convergence testing

1.2 Your First PDE: The 1D Wave Equation

We begin our exploration of Devito with the one-dimensional wave equation, a fundamental PDE that describes vibrations in strings, sound waves in tubes, and many other physical phenomena.

1.2.1 The Mathematical Model

The 1D wave equation is: \[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \tag{1.1}\]

where:

  • \(u(x, t)\) is the displacement at position \(x\) and time \(t\)
  • \(c\) is the wave speed (a constant)

We solve this on a domain \(x \in [0, L]\) for \(t \in [0, T]\) with:

  • Initial conditions: \(u(x, 0) = I(x)\) and \(\frac{\partial u}{\partial t}(x, 0) = 0\)
  • Boundary conditions: \(u(0, t) = u(L, t) = 0\) (fixed ends)

1.2.2 Finite Difference Discretization

Using central differences in both space and time, we approximate: \[ \frac{\partial^2 u}{\partial t^2} \approx \frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\Delta t^2} \] \[ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2} \]

Substituting into 1.1 and solving for \(u_i^{n+1}\): \[ u_i^{n+1} = 2u_i^n - u_i^{n-1} + C^2(u_{i+1}^n - 2u_i^n + u_{i-1}^n) \tag{1.2}\]

where \(C = c\Delta t/\Delta x\) is the Courant number. The scheme is stable for \(C \le 1\).

1.2.3 The Devito Implementation

Let’s implement this step by step:

from devito import Grid, TimeFunction, Eq, Operator
import numpy as np

# Problem parameters
L = 1.0       # Domain length
c = 1.0       # Wave speed
T = 1.0       # Final time
Nx = 100      # Number of grid points
C = 0.5       # Courant number (for stability)

# Derived parameters
dx = L / Nx
dt = C * dx / c
Nt = int(T / dt)

# Create the computational grid
grid = Grid(shape=(Nx + 1,), extent=(L,))

# Create a time-varying field
# time_order=2 because we have second derivative in time
# space_order=2 for standard second-order accuracy
u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)

# Set initial condition: a Gaussian pulse
x = grid.dimensions[0]
x_coord = 0.5 * L  # Center of domain
sigma = 0.1        # Width of pulse
u.data[0, :] = np.exp(-((np.linspace(0, L, Nx+1) - x_coord)**2) / (2*sigma**2))
u.data[1, :] = u.data[0, :]  # Zero initial velocity

# Define the update equation
# u.forward is u at time n+1, u is at time n, u.backward is at time n-1
# u.dx2 is the second spatial derivative
eq = Eq(u.forward, 2*u - u.backward + (c*dt)**2 * u.dx2)

# Create the operator
op = Operator([eq])

# Run the simulation
op(time=Nt, dt=dt)

# The solution is now in u.data
print(f"Simulation complete: {Nt} time steps")
print(f"Max amplitude at t={T}: {np.max(np.abs(u.data[0, :])):.6f}")

1.2.4 Understanding the Code

Let’s examine each component:

Grid creation:

grid = Grid(shape=(Nx + 1,), extent=(L,))

This creates a 1D grid with Nx + 1 points spanning a domain of length L. The grid spacing is automatically computed as dx = L / Nx.

TimeFunction:

u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
  • name='u': The symbolic name for this field
  • time_order=2: We need values at three time levels (\(n-1\), \(n\), \(n+1\)) for the second time derivative
  • space_order=2: Use second-order accurate spatial stencils

Initial conditions:

u.data[0, :] = ...  # u at t=0
u.data[1, :] = ...  # u at t=dt (for zero initial velocity, same as t=0)

The data attribute provides direct access to the underlying NumPy arrays. Index 0 and 1 represent the two most recent time levels.

Update equation:

eq = Eq(u.forward, 2*u - u.backward + (c*dt)**2 * u.dx2)
  • u.forward: The solution at the next time step (\(u^{n+1}\))
  • u: The solution at the current time step (\(u^n\))
  • u.backward: The solution at the previous time step (\(u^{n-1}\))
  • u.dx2: The second spatial derivative, computed using finite differences

Operator and execution:

op = Operator([eq])
op(time=Nt, dt=dt)

The Operator compiles the equations into optimized C code. Calling it runs the time-stepping loop for Nt steps with time increment dt.

1.2.5 Visualizing the Solution

import matplotlib.pyplot as plt

# Get spatial coordinates
x_vals = np.linspace(0, L, Nx + 1)

# Plot the solution at the final time
plt.figure(figsize=(10, 4))
plt.plot(x_vals, u.data[0, :], 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('u')
plt.title(f'Wave equation solution at t = {T}')
plt.grid(True)
plt.show()

1.2.6 The CFL Condition

The Courant-Friedrichs-Lewy (CFL) condition states that for stability: \[ C = \frac{c \Delta t}{\Delta x} \le 1 \]

Physically, this means information cannot travel more than one grid cell per time step. If \(C > 1\), the numerical solution will grow without bound.

Exercise: Try running the code with C = 1.5 and observe what happens.

1.2.7 What Devito Does Behind the Scenes

When you create the Operator, Devito:

  1. Analyzes the symbolic equations
  2. Determines the stencil pattern and data dependencies
  3. Generates optimized C code with:
    • Proper loop ordering for cache efficiency
    • SIMD vectorization where possible
    • OpenMP parallelization for multi-core execution
  4. Compiles the code and caches the result

You can inspect the generated code:

print(op.ccode)

This reveals the low-level implementation that Devito creates automatically.

1.3 Core Devito Abstractions

Devito provides a small set of powerful abstractions for expressing PDEs. Understanding these building blocks is essential for writing effective Devito code.

1.3.1 Grid: The Computational Domain

The Grid defines the discrete domain on which we solve our PDE:

from devito import Grid

# 1D grid: 101 points over [0, 1]
grid_1d = Grid(shape=(101,), extent=(1.0,))

# 2D grid: 101x101 points over [0, 1] x [0, 1]
grid_2d = Grid(shape=(101, 101), extent=(1.0, 1.0))

# 3D grid: 51x51x51 points over [0, 2] x [0, 2] x [0, 2]
grid_3d = Grid(shape=(51, 51, 51), extent=(2.0, 2.0, 2.0))

Key properties:

  • shape: Number of grid points in each dimension
  • extent: Physical size of the domain
  • dimensions: Symbolic dimension objects (x, y, z)
  • spacing: Grid spacing in each dimension (computed automatically)
grid = Grid(shape=(101, 101), extent=(1.0, 1.0))
x, y = grid.dimensions      # Symbolic dimensions
dx, dy = grid.spacing       # Symbolic spacing (h_x, h_y)
print(f"Grid spacing: dx={float(dx)}, dy={float(dy)}")

1.3.2 Function: Static Fields

A Function represents a field that does not change during time-stepping. Use it for material properties, source terms, or any spatially-varying coefficient:

from devito import Function

grid = Grid(shape=(101,), extent=(1.0,))

# Wave velocity field
c = Function(name='c', grid=grid)
c.data[:] = 1500.0  # Constant velocity (m/s)

# Spatially varying velocity
import numpy as np
x_vals = np.linspace(0, 1, 101)
c.data[:] = 1500 + 500 * x_vals  # Linear velocity gradient

The space_order parameter controls the stencil width for derivatives:

# Higher-order derivatives need wider stencils
c = Function(name='c', grid=grid, space_order=4)

1.3.3 TimeFunction: Time-Varying Fields

A TimeFunction represents the solution field that evolves in time:

from devito import TimeFunction

grid = Grid(shape=(101,), extent=(1.0,))

# For first-order time derivatives (diffusion equation)
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)

# For second-order time derivatives (wave equation)
u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)

Key parameters:

  • time_order: Number of time levels needed (1 for first derivative, 2 for second)
  • space_order: Accuracy order for spatial derivatives

Time indexing shortcuts:

Syntax Meaning Mathematical notation
u Current time level \(u^n\)
u.forward Next time level \(u^{n+1}\)
u.backward Previous time level \(u^{n-1}\)
u.dt First time derivative \(\partial u/\partial t\)
u.dt2 Second time derivative \(\partial^2 u/\partial t^2\)

1.3.4 Derivative Notation

Devito provides intuitive notation for spatial derivatives:

Syntax Meaning Stencil
u.dx \(\partial u/\partial x\) Centered difference
u.dy \(\partial u/\partial y\) Centered difference
u.dx2 \(\partial^2 u/\partial x^2\) Second derivative
u.dy2 \(\partial^2 u/\partial y^2\) Second derivative
u.laplace \(\nabla^2 u\) Laplacian (dimension-agnostic)

The laplace operator is particularly useful because it works in any number of dimensions:

# These are equivalent for 2D:
laplacian_explicit = u.dx2 + u.dy2
laplacian_auto = u.laplace

# In 3D, u.laplace automatically becomes u.dx2 + u.dy2 + u.dz2

1.3.5 Eq: Defining Equations

The Eq class creates symbolic equations:

from devito import Eq

# Explicit update: u^{n+1} = expression
update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.laplace)

# Using the solve() helper for implicit forms
from devito import solve

pde = u.dt2 - c**2 * u.laplace  # The PDE residual
update = Eq(u.forward, solve(pde, u.forward))

The solve() function is useful when the update formula is complex. It symbolically solves for the target variable.

1.3.6 Operator: Compilation and Execution

The Operator takes a list of equations and generates executable code:

from devito import Operator

# Single equation
op = Operator([update])

# Multiple equations (e.g., with boundary conditions)
op = Operator([update, bc_left, bc_right])

# Run for Nt time steps
op(time=Nt, dt=dt)

The operator compiles the equations into optimized C code on first execution. Subsequent calls reuse the cached compiled code.

1.3.7 Complete Example: 2D Diffusion

Let’s put these abstractions together for a 2D diffusion problem:

from devito import Grid, TimeFunction, Eq, Operator
import numpy as np

# Create a 2D grid
grid = Grid(shape=(101, 101), extent=(1.0, 1.0))

# Time-varying field (first-order in time for diffusion)
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)

# Parameters
alpha = 0.1  # Diffusion coefficient
dx = 1.0 / 100
F = 0.25     # Fourier number (for stability)
dt = F * dx**2 / alpha

# Initial condition: hot spot in the center
u.data[0, 45:55, 45:55] = 1.0

# The diffusion equation: u_t = alpha * (u_xx + u_yy)
# Using .laplace for dimension-agnostic code
eq = Eq(u.forward, u + alpha * dt * u.laplace)

# Create and run
op = Operator([eq])
op(time=500, dt=dt)

# Visualize
import matplotlib.pyplot as plt
plt.imshow(u.data[0, :, :], origin='lower', cmap='hot')
plt.colorbar(label='Temperature')
plt.title('2D Diffusion')
plt.show()

1.3.8 Summary of Core Abstractions

Abstraction Purpose Key Parameters
Grid Define computational domain shape, extent
Function Static fields (coefficients) name, grid, space_order
TimeFunction Time-varying fields name, grid, time_order, space_order
Eq Define equations LHS, RHS
Operator Compile and execute List of equations

These five abstractions form the foundation of all Devito programs. In the following sections, we’ll see how to handle boundary conditions and verify our numerical solutions.

1.4 Boundary Conditions in Devito

Properly implementing boundary conditions is crucial for accurate PDE solutions. Devito provides several approaches, each suited to different situations.

1.4.1 Dirichlet Boundary Conditions

Dirichlet conditions specify the solution value at the boundary: \[ u(0, t) = g_0(t), \quad u(L, t) = g_L(t) \]

Method 1: Explicit equations

The most direct approach adds equations that set boundary values:

from devito import Grid, TimeFunction, Eq, Operator

grid = Grid(shape=(101,), extent=(1.0,))
u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)

# Get the time dimension for indexing
t = grid.stepping_dim

# Interior update (wave equation)
update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.dx2)

# Boundary conditions: u = 0 at both ends
bc_left = Eq(u[t+1, 0], 0)
bc_right = Eq(u[t+1, 100], 0)

# Include all equations in the operator
op = Operator([update, bc_left, bc_right])

Method 2: Using subdomain

For interior-only updates, use subdomain=grid.interior:

# Update only interior points (automatically excludes boundaries)
update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.dx2,
            subdomain=grid.interior)

# Set boundaries explicitly
bc_left = Eq(u[t+1, 0], 0)
bc_right = Eq(u[t+1, 100], 0)

op = Operator([update, bc_left, bc_right])

The subdomain=grid.interior approach is often cleaner because it explicitly separates the physics (interior PDE) from the boundary treatment.

1.4.2 Neumann Boundary Conditions

Neumann conditions specify the derivative at the boundary: \[ \frac{\partial u}{\partial x}(0, t) = h_0(t), \quad \frac{\partial u}{\partial x}(L, t) = h_L(t) \]

For a zero-flux condition (\(\partial u/\partial x = 0\)), we use the ghost point method. The central difference at the boundary requires a point outside the domain: \[ \frac{\partial u}{\partial x}\bigg|_{i=0} \approx \frac{u_1 - u_{-1}}{2\Delta x} = 0 \]

This gives \(u_{-1} = u_1\), which we substitute into the interior equation:

grid = Grid(shape=(101,), extent=(1.0,))
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
x = grid.dimensions[0]
t = grid.stepping_dim

# Interior update (diffusion equation)
update = Eq(u.forward, u + alpha * dt * u.dx2, subdomain=grid.interior)

# Neumann BC at left (du/dx = 0): use one-sided update
# u_new[0] = u[0] + alpha*dt * 2*(u[1] - u[0])/dx^2
dx = grid.spacing[0]
bc_left = Eq(u[t+1, 0], u[t, 0] + alpha * dt * 2 * (u[t, 1] - u[t, 0]) / dx**2)

# Neumann BC at right (du/dx = 0)
bc_right = Eq(u[t+1, 100], u[t, 100] + alpha * dt * 2 * (u[t, 99] - u[t, 100]) / dx**2)

op = Operator([update, bc_left, bc_right])

1.4.3 Mixed Boundary Conditions

Often we have different conditions on different boundaries:

# Dirichlet on left, Neumann on right
bc_left = Eq(u[t+1, 0], 0)  # u(0,t) = 0
bc_right = Eq(u[t+1, 100], u[t+1, 99])  # du/dx(L,t) = 0 (copy from interior)

op = Operator([update, bc_left, bc_right])

1.4.4 2D Boundary Conditions

For 2D problems, boundary conditions apply to all four edges:

grid = Grid(shape=(101, 101), extent=(1.0, 1.0))
u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)

x, y = grid.dimensions
t = grid.stepping_dim
Nx, Ny = 100, 100

# Interior update
update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.laplace,
            subdomain=grid.interior)

# Dirichlet BCs on all four edges
bc_left = Eq(u[t+1, 0, y], 0)
bc_right = Eq(u[t+1, Nx, y], 0)
bc_bottom = Eq(u[t+1, x, 0], 0)
bc_top = Eq(u[t+1, x, Ny], 0)

op = Operator([update, bc_left, bc_right, bc_bottom, bc_top])

1.4.5 Time-Dependent Boundary Conditions

For boundaries that vary in time, use the time index:

from devito import Constant

# Time-varying amplitude
A = Constant(name='A')

# Sinusoidal forcing at left boundary
# u(0, t) = A * sin(omega * t)
import sympy as sp
omega = 2 * sp.pi  # Angular frequency

# The time value at step n
t_val = t * dt  # Symbolic time value

bc_left = Eq(u[t+1, 0], A * sp.sin(omega * t_val))

# Set the amplitude before running
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt, A=1.0)  # Pass A as keyword argument

1.4.6 Absorbing Boundary Conditions

For wave equations, we often want waves to exit the domain without reflection. A simple first-order absorbing condition is: \[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0 \quad \text{at } x = L \]

This can be discretized as:

# Absorbing BC at right boundary (waves traveling right)
dx = grid.spacing[0]
bc_right_absorbing = Eq(
    u[t+1, Nx],
    u[t, Nx] - c * dt / dx * (u[t, Nx] - u[t, Nx-1])
)

More sophisticated absorbing conditions use damping layers (sponges) near the boundaries. This is covered in detail in Section 2.13.10.

1.4.7 Periodic Boundary Conditions

For periodic domains, the solution wraps around: \[ u(0, t) = u(L, t) \]

Devito doesn’t directly support periodic BCs, but they can be implemented by copying values:

# Periodic BCs: u[0] = u[Nx-1], u[Nx] = u[1]
bc_periodic_left = Eq(u[t+1, 0], u[t+1, Nx-1])
bc_periodic_right = Eq(u[t+1, Nx], u[t+1, 1])

Note: The order of equations matters. Update the interior first, then copy for periodicity.

1.4.8 Best Practices

  1. Use subdomain=grid.interior for interior updates to clearly separate physics from boundary treatment

  2. Check boundary equation order: Boundary equations should typically come after interior updates in the operator

  3. Verify boundary values: After running, check that boundaries have the expected values

  4. Test with known solutions: Use problems with analytical solutions to verify boundary condition implementation

1.4.9 Example: Complete Wave Equation Solver

Here’s a complete example combining interior updates with boundary conditions:

from devito import Grid, TimeFunction, Eq, Operator
import numpy as np

# Setup
L, c, T = 1.0, 1.0, 2.0
Nx = 100
C = 0.9  # Courant number
dx = L / Nx
dt = C * dx / c
Nt = int(T / dt)

# Grid and field
grid = Grid(shape=(Nx + 1,), extent=(L,))
u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
t = grid.stepping_dim

# Initial condition: plucked string
x_vals = np.linspace(0, L, Nx + 1)
u.data[0, :] = np.sin(np.pi * x_vals)
u.data[1, :] = u.data[0, :]  # Zero initial velocity

# Equations
update = Eq(u.forward, 2*u - u.backward + (c*dt)**2 * u.dx2,
            subdomain=grid.interior)
bc_left = Eq(u[t+1, 0], 0)
bc_right = Eq(u[t+1, Nx], 0)

# Solve
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt)

# Verify: solution should return to initial shape at t = 2L/c
print(f"Initial max: {np.max(u.data[1, :]):.6f}")
print(f"Final max: {np.max(u.data[0, :]):.6f}")

For a string with fixed ends and initial shape \(\sin(\pi x)\), the solution oscillates with period \(2L/c\). After one period, it should return to the initial configuration.

1.5 Verification and Convergence Testing

How do we know our numerical solution is correct? Verification is the process of confirming that our code correctly solves the mathematical equations we intended. This section introduces key verification techniques.

1.5.1 The Importance of Verification

Numerical codes can produce plausible-looking but incorrect results due to:

  • Programming errors (typos, off-by-one errors)
  • Incorrect boundary condition implementation
  • Stability violations
  • Insufficient resolution

Systematic verification catches these problems before they corrupt scientific results.

1.5.2 Convergence Rate Testing

The most powerful verification technique is convergence rate testing. For a scheme with truncation error \(O(\Delta x^p)\), the error should decrease as: \[ E(\Delta x) \approx C \Delta x^p \]

By measuring errors at different resolutions, we can estimate \(p\): \[ p \approx \frac{\log(E_1/E_2)}{\log(\Delta x_1/\Delta x_2)} \]

If the measured rate matches the theoretical order, we have strong evidence the implementation is correct.

1.5.3 Implementing a Convergence Test

import numpy as np
from devito import Grid, TimeFunction, Eq, Operator

def solve_wave_equation(Nx, L=1.0, T=0.5, c=1.0, C=0.5):
    """Solve 1D wave equation and return error vs exact solution."""

    dx = L / Nx
    dt = C * dx / c
    Nt = int(T / dt)

    grid = Grid(shape=(Nx + 1,), extent=(L,))
    u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
    t_dim = grid.stepping_dim

    # Initial condition: sin(pi*x)
    x_vals = np.linspace(0, L, Nx + 1)
    u.data[0, :] = np.sin(np.pi * x_vals)
    u.data[1, :] = np.sin(np.pi * x_vals) * np.cos(np.pi * c * dt)

    # Wave equation
    update = Eq(u.forward, 2*u - u.backward + (c*dt)**2 * u.dx2,
                subdomain=grid.interior)
    bc_left = Eq(u[t_dim+1, 0], 0)
    bc_right = Eq(u[t_dim+1, Nx], 0)

    op = Operator([update, bc_left, bc_right])
    op(time=Nt, dt=dt)

    # Exact solution: u(x,t) = sin(pi*x)*cos(pi*c*t)
    t_final = Nt * dt
    u_exact = np.sin(np.pi * x_vals) * np.cos(np.pi * c * t_final)

    # Return max error
    error = np.max(np.abs(u.data[0, :] - u_exact))
    return error, dx


def convergence_test(grid_sizes):
    """Run convergence test and compute rates."""

    errors = []
    dx_values = []

    for Nx in grid_sizes:
        error, dx = solve_wave_equation(Nx)
        errors.append(error)
        dx_values.append(dx)
        print(f"Nx = {Nx:4d}, dx = {dx:.6f}, error = {error:.6e}")

    # Compute convergence rates
    rates = []
    for i in range(len(errors) - 1):
        rate = np.log(errors[i] / errors[i+1]) / np.log(dx_values[i] / dx_values[i+1])
        rates.append(rate)

    print("\nConvergence rates:")
    for i, rate in enumerate(rates):
        print(f"  {grid_sizes[i]} -> {grid_sizes[i+1]}: rate = {rate:.2f}")

    return errors, dx_values, rates


# Run the test
grid_sizes = [20, 40, 80, 160, 320]
errors, dx_values, rates = convergence_test(grid_sizes)

# Check: rates should be close to 2 for second-order scheme
expected_rate = 2.0
assert all(abs(r - expected_rate) < 0.2 for r in rates), \
    f"Convergence rates {rates} differ from expected {expected_rate}"

1.5.4 Method of Manufactured Solutions (MMS)

For problems without analytical solutions, we use the Method of Manufactured Solutions:

  1. Choose a solution \(u_{\text{mms}}(x, t)\) (any smooth function)
  2. Compute the source term by substituting into the PDE
  3. Solve the modified PDE with the computed source
  4. Compare the numerical solution to \(u_{\text{mms}}\)

Example: Diffusion equation

Let’s verify a diffusion solver using MMS:

import sympy as sp

# Symbolic variables
x_sym, t_sym = sp.symbols('x t')
alpha_sym = sp.Symbol('alpha')

# Manufactured solution (arbitrary smooth function)
u_mms = sp.sin(sp.pi * x_sym) * sp.exp(-t_sym)

# Compute required source term: f = u_t - alpha * u_xx
u_t = sp.diff(u_mms, t_sym)
u_xx = sp.diff(u_mms, x_sym, 2)
f_mms = u_t - alpha_sym * u_xx

print("Manufactured solution:")
print(f"  u_mms = {u_mms}")
print(f"Required source term:")
print(f"  f = {sp.simplify(f_mms)}")

Now implement the solver with this source term:

from devito import Grid, TimeFunction, Function, Eq, Operator
import numpy as np

def solve_diffusion_mms(Nx, alpha=1.0, T=0.5, F=0.4):
    """Solve diffusion with MMS source term."""

    L = 1.0
    dx = L / Nx
    dt = F * dx**2 / alpha
    Nt = int(T / dt)

    grid = Grid(shape=(Nx + 1,), extent=(L,))
    u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
    t_dim = grid.stepping_dim

    # Spatial coordinates for evaluation
    x_vals = np.linspace(0, L, Nx + 1)

    # MMS: u = sin(pi*x) * exp(-t)
    # Source: f = sin(pi*x) * exp(-t) * (alpha*pi^2 - 1)
    def u_exact(x, t):
        return np.sin(np.pi * x) * np.exp(-t)

    def f_source(x, t):
        return np.sin(np.pi * x) * np.exp(-t) * (alpha * np.pi**2 - 1)

    # Initial condition from MMS
    u.data[0, :] = u_exact(x_vals, 0)

    # We need to add source term at each time step
    # For simplicity, use time-lagged source
    f = Function(name='f', grid=grid)

    # Update equation with source
    update = Eq(u.forward, u + alpha * dt * u.dx2 + dt * f,
                subdomain=grid.interior)
    bc_left = Eq(u[t_dim+1, 0], 0)  # u_mms(0,t) = 0
    bc_right = Eq(u[t_dim+1, Nx], 0)  # u_mms(1,t) = 0

    op = Operator([update, bc_left, bc_right])

    # Time stepping with source update
    for n in range(Nt):
        t_current = n * dt
        f.data[:] = f_source(x_vals, t_current)
        op(time=1, dt=dt)

    # Compare to exact solution
    t_final = Nt * dt
    u_exact_final = u_exact(x_vals, t_final)
    error = np.max(np.abs(u.data[0, :] - u_exact_final))

    return error, dx


# Convergence test with MMS
print("MMS Convergence Test for Diffusion Equation:")
grid_sizes = [20, 40, 80, 160]
errors = []
dx_vals = []

for Nx in grid_sizes:
    error, dx = solve_diffusion_mms(Nx)
    errors.append(error)
    dx_vals.append(dx)
    print(f"Nx = {Nx:4d}, error = {error:.6e}")

# Compute rates
for i in range(len(errors) - 1):
    rate = np.log(errors[i] / errors[i+1]) / np.log(2)
    print(f"Rate {grid_sizes[i]}->{grid_sizes[i+1]}: {rate:.2f}")

1.5.5 Quick Verification Checks

Before running full convergence tests, use these quick checks:

1. Conservation properties

For problems that should conserve mass or energy:

# Check mass conservation for diffusion with Neumann BCs
mass_initial = np.sum(u.data[1, :]) * dx
mass_final = np.sum(u.data[0, :]) * dx
print(f"Mass change: {abs(mass_final - mass_initial):.2e}")

2. Symmetry

For symmetric initial conditions and domains:

# Check symmetry is preserved
u_left = u.data[0, :Nx//2]
u_right = u.data[0, Nx//2+1:][::-1]  # Reversed
symmetry_error = np.max(np.abs(u_left - u_right))
print(f"Symmetry error: {symmetry_error:.2e}")

3. Steady state

For problems with known steady states:

# Run to steady state and check
u_steady_numerical = u.data[0, :]
u_steady_exact = ...  # Known analytical steady state
error = np.max(np.abs(u_steady_numerical - u_steady_exact))

1.5.6 Debugging Tips

When convergence tests fail:

  1. Check boundary conditions: Are they correctly implemented? Plot the solution near boundaries.

  2. Check stability: Is the CFL/Fourier number within limits? Try smaller time steps.

  3. Check initial conditions: Are they set correctly? Verify u.data[0, :] and u.data[1, :].

  4. Inspect generated code: Use print(op.ccode) to see what Devito actually computes.

  5. Test components separately: Verify spatial derivatives work on known functions before testing full PDE.

1.5.7 Summary

Verification is essential for trustworthy numerical results:

Technique When to Use What It Checks
Convergence testing Always Correct order of accuracy
MMS No analytical solution Correct PDE implementation
Conservation Physics requires it No spurious sources/sinks
Symmetry Symmetric problems Consistent treatment

A well-verified code gives confidence that results represent the physics, not numerical artifacts.