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 (Luporini et al. 2020) 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 arraysThis approach has several limitations:
- Error-prone: Manual index arithmetic is easy to get wrong
- Hard to optimize: Achieving good performance requires expertise in vectorization, parallelization, and cache optimization
- Dimension-specific: The code must be rewritten for 2D or 3D problems
- Not portable: Optimizations for one architecture don’t transfer to others
1.1.2 The Devito Approach
With Devito, the same problem becomes:
This code is included from src/book_snippets/what_is_devito_diffusion.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Constant, Eq, Grid, Operator, TimeFunction, solve
# 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)
# Used by tests
RESULT = float(np.max(u.data[0, :]))
This approach offers significant advantages:
- Mathematical clarity: The PDE
u.dt - a * u.dx2 = 0is written symbolically, and Devito derives the update formula automatically usingsolve() - Automatic optimization: Devito generates C code with loop tiling, SIMD vectorization, and OpenMP parallelization
- Dimension-agnostic: The same code structure works for 1D, 2D, or 3D
- 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 → ExecutionSymbolic representation: Your Python code creates SymPy expressions that represent the PDE and its discretization
Code generation: Devito analyzes the expressions and generates optimized C code with appropriate loop structures
Just-in-time compilation: The C code is compiled (and cached) the first time the operator runs
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) (Louboutin et al. 2019)
1.1.5 Installation
Devito can be installed via pip:
pip install devitoFor 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:
- Solve your first PDE (the 1D wave equation) using Devito
- Understand the core abstractions:
Grid,Function,TimeFunction,Eq, andOperator - Implement boundary conditions in Devito
- 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:
This code is included from src/book_snippets/first_pde_wave1d.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# 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,))
t_dim = grid.stepping_dim
# Create a time-varying field (2nd order in time, 2nd order in space)
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
# Initial condition: Gaussian pulse
x_coords = np.linspace(0, L, Nx + 1)
x0 = 0.5 * L
sigma = 0.1
u0 = np.exp(-((x_coords - x0) ** 2) / (2 * sigma**2))
u.data[0, :] = u0
# First step for zero initial velocity (second-order accurate)
u_xx_0 = np.zeros_like(u0)
u_xx_0[1:-1] = (u0[2:] - 2 * u0[1:-1] + u0[:-2]) / dx**2
u1 = u0 + 0.5 * dt**2 * c**2 * u_xx_0
u1[0] = 0.0
u1[-1] = 0.0
u.data[1, :] = u1
# Update equation (interior) + Dirichlet boundaries
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.0)
bc_right = Eq(u[t_dim + 1, Nx], 0.0)
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt)
# Used by tests
RESULT = float(np.max(np.abs(u.data[0, :])))
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 fieldtime_order=2: We need values at three time levels (\(n-1\), \(n\), \(n+1\)) for the second time derivativespace_order=2: Use second-order accurate spatial stencils
Initial conditions:
u.data[0, :] = ... # u at t=0
u.data[1, :] = ... # u at t=dt (computed from u_t(x,0)=0)The data attribute provides direct access to the underlying NumPy arrays. Index 0 and 1 represent the two most recent time levels.
For the 2nd-order wave scheme, “zero initial velocity” does not mean u.data[1, :] = u.data[0, :] if you want to keep 2nd-order accuracy at the first step. Instead, the included (tested) snippet computes the first time level using the spatial second derivative at t=0:
u1 = u0 + 0.5 * dt**2 * c**2 * u_xx_0and enforces the fixed-end boundary conditions at that first step.
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:
- Analyzes the symbolic equations
- Determines the stencil pattern and data dependencies
- Generates optimized C code with:
- Proper loop ordering for cache efficiency
- SIMD vectorization where possible
- OpenMP parallelization for multi-core execution
- 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 dimensionextent: Physical size of the domaindimensions: 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 gradientThe 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.dz21.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. Boundary conditions fall into two broad categories:
| Category | Type | Description | Typical use |
|---|---|---|---|
| Physical | Dirichlet | Prescribed value \(u = g\) | Fixed ends, walls |
| Neumann | Prescribed flux \(\partial u/\partial n = h\) | Insulation, symmetry | |
| Mixed/Robin | \(\alpha u + \beta \partial u/\partial n = g\) | Convective heat transfer | |
| Periodic | \(u(0) = u(L)\) | Infinite domains, wrapping | |
| Computational | First-order ABC | \(u_t + c\, u_n = 0\) | Open boundaries (1D) |
| Damping layer | Sponge zone with \(\gamma u_t\) | Open boundaries (2D/3D) | |
| PML | Complex coordinate stretching | High-accuracy open boundaries |
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:
This code is included from src/book_snippets/boundary_dirichlet_wave.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# Setup
L, c, T = 1.0, 1.0, 0.2
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,))
t = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
# 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 (demo)
# 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.0)
bc_right = Eq(u[t + 1, Nx], 0.0)
# Solve
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt)
# Used by tests
RESULT = float(max(abs(u.data[0, 0]), abs(u.data[0, -1])))
Method 2: Using subdomain
The snippet above already uses subdomain=grid.interior to keep the interior PDE update separate from boundary treatment.
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:
This code is included from src/book_snippets/neumann_bc_diffusion_1d.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# Neumann boundary conditions: du/dx = 0 at both ends for diffusion.
L = 1.0
Nx = 100
alpha = 1.0
F = 0.4 # stable for Forward Euler diffusion in 1D when F <= 0.5
dx = L / Nx
dt = F * dx**2 / alpha
Nt = 25
grid = Grid(shape=(Nx + 1,), extent=(L,))
t = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=1, space_order=2)
x = np.linspace(0.0, L, Nx + 1)
u.data[0, :] = np.exp(-((x - 0.5) ** 2) / (2 * 0.05**2))
update = Eq(u.forward, u + alpha * dt * u.dx2, subdomain=grid.interior)
dx_sym = grid.spacing[0]
bc_left = Eq(
u[t + 1, 0],
u[t, 0] + alpha * dt * 2.0 * (u[t, 1] - u[t, 0]) / dx_sym**2,
)
bc_right = Eq(
u[t + 1, Nx],
u[t, Nx] + alpha * dt * 2.0 * (u[t, Nx - 1] - u[t, Nx]) / dx_sym**2,
)
op = Operator([update, bc_left, bc_right])
for _ in range(Nt):
op.apply(time_m=0, time_M=0)
u.data[0, :] = u.data[1, :]
grad_left = float(abs(u.data[0, 1] - u.data[0, 0]))
grad_right = float(abs(u.data[0, -1] - u.data[0, -2]))
RESULT = max(grad_left, grad_right)
1.4.3 Mixed Boundary Conditions
Often we have different conditions on different boundaries:
This code is included from src/book_snippets/mixed_bc_diffusion_1d.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# Mixed boundary conditions: Dirichlet on left, Neumann (copy) on right.
L = 1.0
Nx = 80
alpha = 1.0
F = 0.4
dx = L / Nx
dt = F * dx**2 / alpha
grid = Grid(shape=(Nx + 1,), extent=(L,))
t = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=1, space_order=2)
x = np.linspace(0.0, L, Nx + 1)
u.data[0, :] = np.exp(-((x - 0.25) ** 2) / (2 * 0.05**2))
update = Eq(u.forward, u + alpha * dt * u.dx2, subdomain=grid.interior)
bc_left = Eq(u[t + 1, 0], 0.0)
bc_right = Eq(u[t + 1, Nx], u[t + 1, Nx - 1]) # du/dx = 0 (copy trick)
op = Operator([update, bc_left, bc_right])
op.apply(time_m=0, time_M=0)
RESULT = {
"left_boundary": float(u.data[1, 0]),
"right_copy_error": float(abs(u.data[1, -1] - u.data[1, -2])),
}
1.4.4 2D Boundary Conditions
For 2D problems, boundary conditions apply to all four edges:
This code is included from src/book_snippets/bc_2d_dirichlet_wave.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# 2D Dirichlet boundary conditions on all edges (wave equation).
Lx = 1.0
Ly = 1.0
Nx = 51
Ny = 51
c = 1.0
C = 0.5
dx = Lx / (Nx - 1)
dt = C * dx / c
Nt = 5
grid = Grid(shape=(Ny, Nx), extent=(Ly, Lx))
t = grid.stepping_dim
x, y = grid.dimensions
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
xx = np.linspace(0.0, Lx, Nx)
yy = np.linspace(0.0, Ly, Ny)
X, Y = np.meshgrid(xx, yy)
u0 = np.exp(-((X - 0.5) ** 2 + (Y - 0.5) ** 2) / (2 * 0.08**2))
u.data[0, :, :] = u0
u.data[1, :, :] = u0 # demo first step
update = Eq(u.forward, 2 * u - u.backward + (c * dt) ** 2 * u.laplace, subdomain=grid.interior)
bc_left = Eq(u[t + 1, x, 0], 0.0)
bc_right = Eq(u[t + 1, x, Ny - 1], 0.0)
bc_bottom = Eq(u[t + 1, 0, y], 0.0)
bc_top = Eq(u[t + 1, Nx - 1, y], 0.0)
op = Operator([update, bc_left, bc_right, bc_bottom, bc_top])
op(time=Nt, dt=dt)
edges = [
u.data[0, :, 0],
u.data[0, :, -1],
u.data[0, 0, :],
u.data[0, -1, :],
]
RESULT = float(max(np.max(np.abs(e)) for e in edges))
1.4.5 Time-Dependent Boundary Conditions
For boundaries that vary in time, use the time index:
This code is included from src/book_snippets/time_dependent_bc_sine.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# Time-dependent Dirichlet boundary condition: u(0,t) = A*sin(omega*t).
# For time-varying BCs, we loop manually and update the boundary each step.
L = 1.0
Nx = 80
c = 1.0
C = 0.9
dx = L / Nx
dt = C * dx / c
Nt = 10
grid = Grid(shape=(Nx + 1,), extent=(L,))
t_dim = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
u.data[:] = 0.0
# Interior update (wave equation)
update = Eq(u.forward, 2 * u - u.backward + (c * dt) ** 2 * u.dx2, subdomain=grid.interior)
# Time-independent right BC
bc_right = Eq(u[t_dim + 1, Nx], 0.0)
# Create operator without the time-dependent left BC
op = Operator([update, bc_right])
# Amplitude and frequency
A = 1.0
omega = 2 * np.pi
# Time-stepping loop with manual BC update
for n in range(Nt):
t_val = n * dt
# Set time-dependent BC at left boundary
u.data[(n + 1) % 3, 0] = A * np.sin(omega * t_val)
op(time=1, dt=dt)
# Check that the left boundary has non-zero values (was driven by sine)
RESULT = float(np.max(np.abs(u.data[:, 0])))
1.4.6 Absorbing / Open Boundary Conditions
For wave equations on truncated domains, we need boundaries that absorb outgoing waves without spurious reflections. Several approaches exist, ranging from simple to sophisticated:
First-order ABC (Engquist and Majda 1977; Clayton and Engquist 1977). Based on the one-way wave equation: \[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0 \quad \text{at } x = L \]
This is exact for normally incident waves but reflects oblique waves in 2D/3D. It can be discretized as:
This code is included from src/book_snippets/absorbing_bc_right_wave.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# First-order absorbing boundary condition at the right boundary (1D wave).
L = 1.0
Nx = 200
c = 1.0
C = 0.9
dx = L / Nx
dt = C * dx / c
Nt = 10
grid = Grid(shape=(Nx + 1,), extent=(L,))
t = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
x = np.linspace(0.0, L, Nx + 1)
u.data[0, :] = np.exp(-((x - 0.8) ** 2) / (2 * 0.03**2))
u.data[1, :] = u.data[0, :]
update = Eq(u.forward, 2 * u - u.backward + (c * dt) ** 2 * u.dx2, subdomain=grid.interior)
bc_left = Eq(u[t + 1, 0], 0.0)
dx_sym = grid.spacing[0]
bc_right_absorbing = Eq(u[t + 1, Nx], u[t, Nx] - c * dt / dx_sym * (u[t, Nx] - u[t, Nx - 1]))
op = Operator([update, bc_left, bc_right_absorbing])
op(time=Nt, dt=dt)
RESULT = float(np.max(np.abs(u.data[0, :])))
Damping layers (sponge zones) (Cerjan et al. 1985; Sochacki et al. 1987). Add a dissipative term \(\gamma(\mathbf{x})\, u_t\) to the PDE in a region near the boundary. The damping coefficient ramps from zero in the interior to a maximum at the boundary. Effective with 20–40 cells of absorbing layer.
Perfectly Matched Layer (PML) (Berenger 1994). Uses complex coordinate stretching to create a layer with zero theoretical reflection at any angle and frequency. Requires auxiliary field variables but achieves excellent absorption with only 10–15 cells. The CPML variant (Roden and Gedney 2000) improves stability.
Higher-order ABCs (Higdon 1986, 1987). Generalize the first-order condition to absorb waves at multiple angles simultaneously. Combined with a thin damping layer, the hybrid HABC-Higdon method achieves up to 99% reflection reduction (Dolci et al. 2022).
For detailed implementations and comparisons of all these methods, see Section 2.50.
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:
This code is included from src/book_snippets/periodic_bc_advection_1d.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Constant, Eq, Grid, Operator, TimeFunction
# Periodic boundary conditions using copy equations (1D advection).
L = 1.0
Nx = 80
c = 1.0
C = 0.8
dx = L / Nx
dt = C * dx / c
grid = Grid(shape=(Nx + 1,), extent=(L,))
(x_dim,) = grid.dimensions
t = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=1, space_order=1)
x = np.linspace(0.0, L, Nx + 1)
u.data[0, :] = np.exp(-0.5 * ((x - 0.25) / 0.05) ** 2)
u.data[1, :] = u.data[0, :]
courant = Constant(name="C", value=C)
update = Eq(u.forward, u - courant * (u - u.subs(x_dim, x_dim - x_dim.spacing)))
bc_left = Eq(u[t + 1, 0], u[t, Nx])
bc_right = Eq(u[t + 1, Nx], u[t + 1, 0])
op = Operator([update, bc_left, bc_right])
op.apply(time_m=0, time_M=0, dt=dt)
RESULT = float(abs(u.data[1, 0] - u.data[1, -1]))
Note: The order of equations matters. Update the interior first, then copy for periodicity.
1.4.8 Best Practices
Use
subdomain=grid.interiorfor interior updates to clearly separate physics from boundary treatmentCheck boundary equation order: Boundary equations should typically come after interior updates in the operator
Verify boundary values: After running, check that boundaries have the expected values
Test with known solutions: Use problems with analytical solutions to verify boundary condition implementation
Verify ABC effectiveness: For absorbing boundaries, always compare against a reference solution on a larger domain to confirm that reflections are below an acceptable threshold (see Section 2.50.7)
1.4.9 Example: Complete Wave Equation Solver
Here’s a complete example combining interior updates with boundary conditions:
This code is included from src/book_snippets/boundary_dirichlet_wave.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
# Setup
L, c, T = 1.0, 1.0, 0.2
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,))
t = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
# 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 (demo)
# 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.0)
bc_right = Eq(u[t + 1, Nx], 0.0)
# Solve
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt)
# Used by tests
RESULT = float(max(abs(u.data[0, 0]), abs(u.data[0, -1])))
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 (Roache 2009) 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 Floating-Point Precision and Convergence Testing
Convergence rate testing requires measuring how the error \(E(\Delta x)\) decreases across several grid refinements. The total error has two components: \[ E_{\text{total}} = C \Delta x^p + E_{\text{roundoff}} \tag{1.3}\] As \(\Delta x \to 0\), the discretization term \(C \Delta x^p\) shrinks but eventually reaches the round-off error floor set by the floating-point precision. The machine epsilon for single precision (FP32) is approximately \(1.2 \times 10^{-7}\), while for double precision (FP64) it is approximately \(2.2 \times 10^{-16}\).
For a second-order scheme (\(p=2\)), the discretization error is \(O(\Delta x^2)\). Consider the number of useful refinement levels before round-off dominates:
| Grid points | \(\Delta x\) | \(O(\Delta x^2)\) | FP32 resolved? | FP64 resolved? |
|---|---|---|---|---|
| 10 | \(10^{-1}\) | \(10^{-2}\) | Yes | Yes |
| 100 | \(10^{-2}\) | \(10^{-4}\) | Yes | Yes |
| 1000 | \(10^{-3}\) | \(10^{-6}\) | Marginal | Yes |
| 10000 | \(10^{-4}\) | \(10^{-8}\) | No | Yes |
With FP32, only 2–3 useful refinement levels are available before round-off noise corrupts the convergence rate estimate. With FP64, 5–7 levels are available — enough to robustly establish the asymptotic convergence rate. Roy (Roy 2005) emphasises that the observed order of accuracy is only meaningful when the solution is in the asymptotic range, where higher-order error terms are negligible. FP32 often lacks the headroom to both reach the asymptotic range and have sufficient refinement levels before hitting the round-off floor.
All solver functions in this book default to dtype=np.float64 (double precision) because code verification requires measuring errors across many orders of magnitude. Users may pass dtype=np.float32 for production runs where throughput matters more than verification precision.
1.5.4 Implementing a Convergence Test
This code is included from src/book_snippets/verification_convergence_wave.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
def solve_wave_equation(Nx, L=1.0, T=0.5, c=1.0, C=0.5):
dx = L / Nx
dt = C * dx / c
Nt = int(T / dt)
grid = Grid(shape=(Nx + 1,), extent=(L,))
t_dim = grid.stepping_dim
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=2)
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)
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.0)
bc_right = Eq(u[t_dim + 1, Nx], 0.0)
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt)
t_final = Nt * dt
u_exact = np.sin(np.pi * x_vals) * np.cos(np.pi * c * t_final)
# For time_order=2, buffer has 3 slots; final solution is at Nt % 3
final_idx = Nt % 3
error = float(np.max(np.abs(u.data[final_idx, :] - u_exact)))
return error, dx
def convergence_test(grid_sizes):
errors = []
dx_values = []
for Nx in grid_sizes:
error, dx = solve_wave_equation(Nx)
errors.append(error)
dx_values.append(dx)
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(float(rate))
return rates
# Use grid sizes that avoid numerical resonance issues
RESULT = convergence_test([25, 50, 100, 200])
1.5.5 Method of Manufactured Solutions (MMS)
For problems without analytical solutions, we use the Method of Manufactured Solutions:
- Choose a solution \(u_{\text{mms}}(x, t)\) (any smooth function)
- Compute the source term by substituting into the PDE
- Solve the modified PDE with the computed source
- Compare the numerical solution to \(u_{\text{mms}}\)
Example: Computing MMS source terms
For the diffusion equation \(u_t = \alpha u_{xx}\), we can compute the required source term for any manufactured solution:
This code is included from src/book_snippets/verification_mms_symbolic.py and exercised by pytest (see tests/test_book_snippets.py).
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
# Verify the expressions
RESULT = {
"u_mms": str(u_mms),
"f_mms": str(sp.simplify(f_mms)),
}
Practical approach: eigenfunction solutions
For diffusion problems, a simpler approach uses exact eigenfunction solutions that require no source term. The solution \(u(x,t) = \sin(\pi x) e^{-\alpha \pi^2 t}\) satisfies both the PDE and homogeneous boundary conditions:
This code is included from src/book_snippets/verification_mms_diffusion.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Constant, Eq, Grid, Operator, TimeFunction, solve
def solve_diffusion_exact(Nx, alpha=1.0, T=0.1, F=0.4):
"""Solve diffusion equation and compare with exact eigenfunction solution.
Uses exact solution: u(x,t) = sin(pi*x) * exp(-alpha*pi^2*t)
which satisfies u_t = alpha * u_xx with u(0,t) = u(L,t) = 0.
"""
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
x_vals = np.linspace(0, L, Nx + 1)
# Exact solution: eigenfunction of diffusion operator
def u_exact(x, t):
return np.sin(np.pi * x) * np.exp(-alpha * np.pi**2 * t)
# Initial condition from exact solution
u.data[0, :] = u_exact(x_vals, 0)
# Diffusion equation: u_t = alpha * u_xx
a = Constant(name="a")
pde = u.dt - a * u.dx2
update = Eq(u.forward, solve(pde, u.forward), subdomain=grid.interior)
bc_left = Eq(u[t_dim + 1, 0], 0.0)
bc_right = Eq(u[t_dim + 1, Nx], 0.0)
op = Operator([update, bc_left, bc_right])
# Run all time steps at once
op(time=Nt, dt=dt, a=alpha)
# Compare to exact solution
t_final = Nt * dt
u_exact_final = u_exact(x_vals, t_final)
# Determine which buffer has the final solution
final_idx = Nt % 2
error = float(np.max(np.abs(u.data[final_idx, :] - u_exact_final)))
return error, dx
def convergence_test_mms(grid_sizes):
"""Run MMS convergence test for diffusion equation."""
errors = []
dx_vals = []
for Nx in grid_sizes:
error, dx = solve_diffusion_exact(Nx)
errors.append(error)
dx_vals.append(dx)
# Compute rates
rates = []
for i in range(len(errors) - 1):
rate = np.log(errors[i] / errors[i + 1]) / np.log(2)
rates.append(float(rate))
return rates
RESULT = convergence_test_mms([20, 40, 80, 160])
1.5.6 Quick Verification Checks
Before running full convergence tests, use these quick checks for conservation and symmetry properties:
This code is included from src/book_snippets/verification_quick_checks.py and exercised by pytest (see tests/test_book_snippets.py).
import numpy as np
from devito import Eq, Grid, Operator, TimeFunction
def check_mass_conservation(Nx=50, alpha=1.0, T=0.1, F=0.4):
"""Check mass conservation for diffusion with Neumann BCs (approximated)."""
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)
# Symmetric initial condition
x_vals = np.linspace(0, L, Nx + 1)
u.data[0, :] = np.exp(-((x_vals - 0.5) ** 2) / 0.01)
# Diffusion with zero-flux BCs (approximate via copying)
t_dim = grid.stepping_dim
update = Eq(u.forward, u + alpha * dt * u.dx2, subdomain=grid.interior)
bc_left = Eq(u[t_dim + 1, 0], u[t_dim + 1, 1])
bc_right = Eq(u[t_dim + 1, Nx], u[t_dim + 1, Nx - 1])
op = Operator([update, bc_left, bc_right])
mass_initial = float(np.sum(u.data[0, :]) * dx)
op(time=Nt, dt=dt)
mass_final = float(np.sum(u.data[0, :]) * dx)
return abs(mass_final - mass_initial)
def check_symmetry(Nx=50, alpha=1.0, T=0.1, F=0.4):
"""Check symmetry preservation for symmetric initial conditions."""
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)
# Symmetric initial condition centered at L/2
x_vals = np.linspace(0, L, Nx + 1)
u.data[0, :] = np.exp(-((x_vals - 0.5) ** 2) / 0.01)
t_dim = grid.stepping_dim
update = Eq(u.forward, u + alpha * dt * u.dx2, subdomain=grid.interior)
bc_left = Eq(u[t_dim + 1, 0], 0.0)
bc_right = Eq(u[t_dim + 1, Nx], 0.0)
op = Operator([update, bc_left, bc_right])
op(time=Nt, dt=dt)
# Check symmetry: left half vs reversed right half
u_left = u.data[0, : Nx // 2]
u_right = u.data[0, Nx // 2 + 1 :][::-1]
symmetry_error = float(np.max(np.abs(u_left - u_right)))
return symmetry_error
RESULT = {
"mass_change": check_mass_conservation(),
"symmetry_error": check_symmetry(),
}
1.5.7 Debugging Tips
When convergence tests fail:
Check boundary conditions: Are they correctly implemented? Plot the solution near boundaries.
Check stability: Is the CFL/Fourier number within limits? Try smaller time steps.
Check initial conditions: Are they set correctly? Verify
u.data[0, :]andu.data[1, :].Inspect generated code: Use
print(op.ccode)to see what Devito actually computes.Test components separately: Verify spatial derivatives work on known functions before testing full PDE.
1.5.8 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.