Back to Details

devito

Compare original and translation side by side

🇺🇸

Original

English
🇨🇳

Translation

Chinese

Devito - Symbolic PDE Solver

Devito - 符号型偏微分方程求解器

Quick Reference

快速参考

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

Create grid

Create grid

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

Velocity model

Velocity model

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

Wavefield

Wavefield

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

Wave equation: d2p/dt2 = v^2 * laplacian(p)

Wave equation: d2p/dt2 = v^2 * laplacian(p)

stencil = Eq(p.forward, solve(p.dt2 - v**2 * p.laplace, p.forward))
stencil = Eq(p.forward, solve(p.dt2 - v**2 * p.laplace, p.forward))

Compile and run

Compile and run

op = Operator([stencil]) op(time_M=100, dt=0.5)
undefined
op = Operator([stencil]) op(time_M=100, dt=0.5)
undefined

Key Classes

核心类

ClassPurpose
Grid
Computational domain definition
Function
Spatial field on grid
TimeFunction
Time-dependent field
SparseTimeFunction
Point sources/receivers
Operator
Compiled computation kernel
类用途
Grid
计算域定义
Function
网格上的空间场
TimeFunction
随时间变化的场
SparseTimeFunction
点源/接收器
Operator
编译后的计算内核

Essential Operations

基础操作

Grid and Fields

网格与场

python
from devito import Grid, Function, TimeFunction
python
from devito import Grid, Function, TimeFunction

2D/3D Grid

2D/3D Grid

grid = Grid(shape=(nx, nz), extent=(x_size, z_size))
grid = Grid(shape=(nx, nz), extent=(x_size, z_size))

Velocity model (spatial field)

Velocity model (spatial field)

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

Wavefield (time-dependent)

Wavefield (time-dependent)

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

Source and Receivers

源与接收器

python
from examples.seismic import RickerSource, Receiver, TimeAxis

time_range = TimeAxis(start=0., stop=1000., step=dt)
python
from examples.seismic import RickerSource, Receiver, TimeAxis

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

Source

Source

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

Receivers

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.
undefined
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.
undefined

Build and Run

构建与运行

python
undefined
python
undefined

Wave equation

Wave equation

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

Compile and execute

Compile and execute

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

Results

Results

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

Symbolic Derivatives

符号导数

SyntaxDescription
p.dt
, 
p.dt2
First/second time derivative
p.dx
, 
p.dy
, 
p.dz
Spatial derivatives
p.laplace
Laplacian (auto-adapts to dims)
p.forward
p at t+dt (time stepping)
p.backward
p at t-dt (adjoint)
语法描述
p.dt
, 
p.dt2
一阶/二阶时间导数
p.dx
, 
p.dy
, 
p.dz
空间导数
p.laplace
拉普拉斯算子（自动适配维度）
p.forward
t+dt时刻的p（时间步进）
p.backward
t-dt时刻的p（伴随运算）

Stability and Accuracy

稳定性与精度

CFL Condition: 
dt < dx / (v_max * sqrt(ndim))
DimsMax dt
2Ddx / (v_max * 1.414)
3Ddx / (v_max * 1.732)
Space OrderStencil PointsError
23O(h^2)
45O(h^4)
89O(h^8)
Higher order = more accurate but slower. Use 4-8 for production.
CFL条件：
dt < dx / (v_max * sqrt(ndim))
维度最大dt
2Ddx / (v_max * 1.414)
3Ddx / (v_max * 1.732)
空间阶数模板点数误差
23O(h^2)
45O(h^4)
89O(h^8)
阶数越高，精度越高但速度越慢。生产环境建议使用4-8阶。

When to Use vs Alternatives

适用场景与替代方案

ScenarioRecommendation
Seismic wave propagation (acoustic/elastic)Devito - symbolic PDE, auto-optimized code
Full Waveform Inversion (FWI) or RTMDevito - adjoint support, GPU-ready
Legacy seismic processing pipelinesMadagascar - established, large script library
Simple 1D/2D wave demosCustom NumPy - no dependencies, easier to debug
General-purpose PDE solving (non-wave)FEniCS - FEM-based, broader PDE support
Production seismic imaging at scaleDevito - generates optimized C code, MPI support
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.
场景推荐方案
地震波传播（声波/弹性波）Devito - 符号型PDE，自动优化代码
全波形反演（FWI）或逆时偏移（RTM）Devito - 支持伴随运算，适配GPU
传统地震处理流水线Madagascar - 成熟稳定，拥有大量脚本库
简单1D/2D波演示自定义NumPy实现 - 无依赖，调试更简单
通用PDE求解（非波类）FEniCS - 基于有限元法，支持更广泛的PDE
大规模生产级地震成像Devito - 生成优化的C代码，支持MPI
选择Devito的场景：你需要高性能的有限差分波传播模拟，且希望通过符号化方式定义方程。它能从Python层面的数学表达式自动生成优化的C/OpenMP/GPU代码，非常适合FWI、RTM以及研究原型开发。
避免使用Devito的场景：你需要有限元法（请使用FEniCS），或者简单的教学示例（NumPy即可满足需求）。

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
  • 定义与速度模型匹配的
    Grid
    （包含形状和物理范围）
  • 创建速度
    Function
    并填充模型数值
  • 为波场创建
    TimeFunction
    （time_order=2，space_order≥4）
  • 验证CFL条件：
    dt < dx / (v_max * sqrt(ndim))
  • 构建波方程模板：
    Eq(p.forward, solve(p.dt2 - v**2 * p.laplace, p.forward))
  • 创建源（
    RickerSource
    ）和接收器，设置坐标
  • 添加源注入和接收器插值项
  • 用模板+源+接收器项编译
    Operator
  • 运行算子：
    op(time_M=nt-1, dt=dt)
  • 从
    rec.data
    提取炮集记录并绘图

References

参考资料

  • Operators and Stencils - Detailed operator construction
  • Performance Optimization - GPU execution and tuning
  • 算子与模板 - 算子构建详细说明
  • 性能优化 - GPU执行与调优

Scripts

脚本

  • scripts/acoustic_wave.py - Basic acoustic wave modeling
  • scripts/acoustic_wave.py - 基础声波模拟脚本