devito

Devito - 符号 PDE 求解器

快速参考

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)

主要类

类	用途
Grid	定义计算域
Function	网格上的空间场
TimeFunction	随时间变化的场
SparseTimeFunction	点源/接收器
Operator	编译后的计算内核

基本操作

网格和场

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)

源和接收器

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.

构建和运行

# 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

符号导数

语法	描述
p.dt, p.dt2	一阶/二阶时间导数
p.dx, p.dy, p.dz	空间导数
p.laplace	拉普拉斯算子（自动适配维度）
p.forward	t+dt 时刻的 p（时间步进）
p.backward	t-dt 时刻的 p（伴随）

稳定性和精度

CFL 条件：dt < dx / (v_max * sqrt(ndim))

维度	最大 dt
2D	dx / (v_max * 1.414)
3D	dx / (v_max * 1.732)

空间阶数	模板点数	误差
2	3	O(h^2)
4	5	O(h^4)
8	9	O(h^8)

阶数越高，精度越高，但速度越慢。生产环境中建议使用 4-8 阶。

适用场景与替代方案

场景	建议
地震波传播（声学/弹性）	Devito - 符号 PDE，自动优化代码
全波形反演（FWI）或逆时偏移（RTM）	Devito - 支持伴随计算，可用于 GPU
遗留地震处理流水线	Madagascar - 成熟稳定，脚本库丰富
简单的 1D/2D 波动演示	Custom NumPy - 无依赖，更易调试
通用 PDE 求解（非波动问题）	FEniCS - 基于 FEM，支持更广泛的 PDE
大规模生产级地震成像	Devito - 生成优化的 C 代码，支持 MPI

选择 Devito 的情况：需要高性能有限差分波传播，并使用符号方程进行描述。它可以根据 Python 层面的数学表达式自动生成优化的 C/OpenMP/GPU 代码，非常适合 FWI、RTM 和研究原型开发。

避免使用 Devito 的情况：需要有限元方法（使用 FEniCS），或只需 NumPy 就能实现的简单教学示例。

常见工作流

使用源和接收器进行声波正演建模

• 定义 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 提取炮集记录并绘图

参考资料

• 算子和模板 - 算子构造详解
• 性能优化 - GPU 执行和调优

脚本

• scripts/acoustic_wave.py - 基本声波建模