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	1 階／2 階時間微分
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 - 音響波モデリングの基本例