devito

name: devito description: | Symbolic PDE solver with automatic code generation for finite-difference computations. Use when Claude needs to: (1) Perform seismic wave propagation modeling, (2) Implement acoustic or elastic wave equations, (3) Run forward modeling for shot gathers, (4) Set up Full Waveform Inversion (FWI) workflows, (5) Implement Reverse Time Migration (RTM), (6) Create absorbing boundary conditions, (7) Generate optimized stencil code for CPUs/GPUs, (8) Solve custom PDEs with finite differences. version: 1.0.0 author: Geoscience Skills license: MIT tags: [PDE Solver, Wave Propagation, Seismic Modelling, FWI, RTM, Finite Difference] dependencies: [devito>=4.8.0, numpy] complements: [pylops, simpeg, segyio] workflow_role: modelling

Devito - Symbolic PDE Solver

Quick Reference

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

# Create grid
grid = Grid(shape=(101, 101), extent=(1000., 1000.))

# Velocity model
v = Function(name='v', grid=grid, space_order=4)
v.data[:] = 1500.

# Wavefield
p = TimeFunction(name='p', grid=grid, time_order=2, space_order=4)

# Wave equation: d2p/dt2 = v^2 * laplacian(p)
stencil = Eq(p.forward, solve(p.dt2 - v**2 * p.laplace, p.forward))

# Compile and run
op = Operator([stencil])
op(time_M=100, dt=0.5)

Key Classes

Essential Operations

Grid and Fields

from devito import Grid, Function, TimeFunction

# 2D/3D Grid
grid = Grid(shape=(nx, nz), extent=(x_size, z_size))

# Velocity model (spatial field)
v = Function(name='v', grid=grid, space_order=4)
v.data[:] = 1500.

# Wavefield (time-dependent)
p = TimeFunction(name='p', grid=grid, time_order=2, space_order=4)

Source and Receivers

from examples.seismic import RickerSource, Receiver, TimeAxis

time_range = TimeAxis(start=0., stop=1000., step=dt)

# Source
src = RickerSource(name='src', grid=grid, f0=10., npoint=1, time_range=time_range)
src.coordinates.data[0, :] = [500., 20.]

# Receivers
rec = Receiver(name='rec', grid=grid, npoint=101, time_range=time_range)
rec.coordinates.data[:, 0] = np.linspace(0., 1000., 101)
rec.coordinates.data[:, 1] = 20.

Build and Run

# Wave equation
stencil = Eq(p.forward, solve(p.dt2 - v**2 * p.laplace, p.forward))
src_term = src.inject(field=p.forward, expr=src * dt**2 * v**2)
rec_term = rec.interpolate(expr=p)

# Compile and execute
op = Operator([stencil] + src_term + rec_term)
op(time_M=nt-1, dt=dt)

# Results
shot_record = rec.data        # (nt, nrec)
snapshot = p.data[0]          # Current wavefield

Symbolic Derivatives

Stability and Accuracy

CFL Condition: dt < dx / (v_max * sqrt(ndim))

Higher order = more accurate but slower. Use 4-8 for production.

When to Use vs Alternatives

Choose Devito when: You need high-performance finite-difference wave propagation with symbolic equation specification. It auto-generates optimized C/OpenMP/GPU code from Python-level math, making it ideal for FWI, RTM, and research prototyping.

Avoid Devito when: You need finite-element methods (use FEniCS), or simple pedagogical examples where NumPy suffices.

Common Workflows

Acoustic wave forward modelling with sources and receivers

• Define Grid with shape and physical extent matching the velocity model
• Create velocity Function and populate with model values
• Create TimeFunction for the wavefield (time_order=2, space_order=4+)
• Verify CFL condition: dt < dx / (v_max * sqrt(ndim))
• Build wave equation stencil: Eq(p.forward, solve(p.dt2 - v**2 * p.laplace, p.forward))
• Create source (RickerSource) and receivers, set coordinates
• Add source injection and receiver interpolation terms
• Compile Operator with stencil + source + receiver terms
• Run operator: op(time_M=nt-1, dt=dt)
• Extract shot record from rec.data and plot

References

• Operators and Stencils - Detailed operator construction
• Performance Optimization - GPU execution and tuning

Scripts

• scripts/acoustic_wave.py - Basic acoustic wave modeling