Scripting API

This page covers the full Python API: building a model, creating a solver, running propagation, and accessing results.

1. Building a model

layered_model — quick setup

layered_model builds a stratified WaveModel and initialises all physical parameters required for the requested physics.

import numpy as np
from recipes import layered_model

# 2-D, isotropic acoustic
model = layered_model(
    'iso-acoustic',
    shape=(101, 101),
    spacing=(10.0, 10.0),   # meters
    nbl=40,
    space_order=8,
    nlayers=4,
)

The first argument is a physics string; it controls which secondary parameters are generated:

String contains Parameters added
tti or vti epsilon, delta, theta (and phi in 3-D)
elastic vs → internally converted to lam and mu
visco invqp (1/Q, uniform)

Additional keyword arguments:

Argument Default Description
density False Derive buoyancy b = 1/ρ from Gardner’s relation
water_layer False Add an isotropic water layer at the top
smooth False Gaussian-smooth Thomsen parameters after creation
Pdtype np.float32 Physical parameter storage precision

WaveModel — full control

For production use, construct WaveModel directly with your own velocity and parameter arrays.

from recipes.model import WaveModel

shape   = (201, 201)
spacing = (10., 10.)

vp  = np.ones(shape, dtype=np.float32) * 2.0   # km/s
eps = np.zeros(shape, dtype=np.float32)
dlt = np.zeros(shape, dtype=np.float32)
tht = np.zeros(shape, dtype=np.float32)

model = WaveModel(
    origin=(0., 0.),
    spacing=spacing,
    shape=shape,
    space_order=8,
    vp=vp,
    nbl=40,
    fs=False,
    dtype=np.float32,
    Pdtype=np.float32,
    epsilon=eps,
    delta=dlt,
    theta=tht,
)

WaveModel constructor arguments

Argument Description
vp P-wave velocity array (km/s), shape (nx, [ny,] nz)
shape Grid dimensions without boundary padding
spacing Grid spacing (m)
nbl Absorbing boundary layer width
fs Conventional free surface (True/False)
dtype Grid data type
Pdtype Physical parameter storage precision
b Buoyancy 1/ρ; scalar 1 means constant density
epsilon, delta, theta, phi Thomsen anisotropy parameters
vs S-wave velocity — automatically computes lam and mu
invqp Inverse Q factor for viscoacoustic attenuation
sdf Signed distance function for immersed-boundary topography
interp_order Point source/receiver interpolation order
water_layer Build water/earth subdomain partition

After construction, physical parameters are accessible as attributes:

model.vp          # devito Function
model.m           # 1/vp² (derived)
model.b           # buoyancy
model.epsilon     # (TTI only)
model.critical_dt # maximum stable time step
model.damp        # absorbing boundary field (initialised lazily)

2. Creating a solver

Pick a solver class from recipes_registry and pass the model and an options dictionary:

from recipes import recipes_registry

options = {
    'space_order': 8,
    'nt': 600,
    'f0': 0.025,          # Hz — dominant Ricker frequency
    'Wdtype': np.float32,
    'time_slots': 3,      # rolling wavefield buffer depth
    't_sub': 4,           # snapshot subsampling factor
}

Solver = recipes_registry['iso-acoustic']
solver = Solver(model, options)

Solver options reference

Key Type Default Description
space_order int 4 Spatial FD stencil order
nt int 100 Number of time steps
f0 float 0.010 Dominant source frequency (Hz)
t_sub int 4 Snapshot subsampling factor
time_slots int time_order+1 Rolling wavefield buffer depth
save str 'lazy' Storage strategy: 'lazy' (host streaming), 'disk', or None
compression str|None None Wavefield compression ('cvxcompress', 'bitcomp')
coefficients str 'taylor' FD coefficient type ('taylor' or 'lsqr_drp')
interpolation str 'linear' Source/receiver interpolation ('linear' or 'sinc')
transient bool False Use transient TimeFunction to reduce memory traffic
factorization int 1 Stiffness-tensor precompute mode for TTI elastic (0/1/2)
fd_kernel_only bool False Skip source and receivers; use Gaussian IC
abox bool False Adaptive bounding-box subdomain (active wavefield region)
Wdtype dtype np.float32 Wavefield data type
save_mode str 'all' Elastic: 'all' saves velocity and stress; 'v' saves velocity only
opt tuple auto Devito Operator optimisation options, e.g. ('advanced', {'openmp': True})

3. Running propagation

Forward pass

# Returns a Devito PerformanceSummary
summary = solver.forward()
print(summary)

To save wavefields for gradient computation:

summary = solver.forward(save=True)

Accessing results

# Last pressure snapshot as a numpy array
p = solver.pressure_data

# Recorded traces at receivers — shape (nt, nrec)
traces = solver.rec.data

Source and receiver objects

The default acquisition uses a single Ricker point source placed at the centre of the domain and a surface receiver array:

src = solver.src   # devito RickerSource
rec = solver.rec   # devito Receiver

# Reposition the source
src.coordinates.data[0] = [500., 20.]   # [x, z] in meters

# Inspect the receiver layout
print(rec.coordinates.data)             # (nrec, ndim) array

Both objects are built lazily and cached. If you need non-default acquisition geometries, override acquisition_eq — see Extending the Framework.

Gradient computation

After a forward pass with save=True, inject residuals into rec.data and call jacobian_adjoint:

solver.rec.data[:] = observed - solver.rec.data   # residuals

summary = solver.jacobian_adjoint(params='vp')

# Gradient is accumulated in a persistent Function
dm = solver.perturbation(param='vp')
print(dm.data.shape)

Multiple parameters in one adjoint pass:

solver.jacobian_adjoint(params=('vp', 'epsilon'))
dm_vp  = solver.perturbation(param='vp')
dm_eps = solver.perturbation(param='epsilon')

4. Wavefield compression

When compression is set, saved wavefields are stored in compressed form. The compression ratio is available after the forward pass:

import numpy as np

options = {
    'compression': 'cvxcompress',
    'save': 'lazy',
    'Wdtype': np.float16,
}
solver = recipes_registry['iso-acoustic'](model, options)
solver.forward(save=True)

cr = solver.saved_state['u_save'].compression_ratio
print(f"Compression ratio: {float(np.mean(cr)):.2f}")

5. Per-operator optimisation

Pass an opt dictionary keyed by operator name to apply different optimisation levels to the forward and gradient passes:

opt = {
    'ForwardAcousticIsotropic':    ('advanced', {'openmp': True}),
    'JacobianAdjAcousticIsotropic': ('advanced', {}),
}
solver.forward(save=True, opt=opt)
solver.jacobian_adjoint(params='vp', opt=opt)
Back to top