Solver

Every solver is an instance of a Propagator subclass (recipes/solver.py). This page documents the public API: how to instantiate a solver, run propagation, and access results.

1. Instantiation

from recipes import recipes_registry

Solver = recipes_registry['iso-acoustic']   # or any other key
solver = Solver(model, options)

options is a plain dict; missing keys fall back to the defaults in the table below.

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 Ricker source frequency (Hz)
time_slots int time_order + 1 Wavefield rolling-buffer depth
t_sub int 4 Snapshot subsampling factor (every n-th step is saved)
save str 'lazy' Storage for saved wavefields: 'lazy' (host streaming), 'disk', or None
compression str|None None Wavefield compression ('cvxcompress', 'bitcomp')
coefficients str 'taylor' FD coefficient scheme ('taylor' or 'lsqr_drp')
interpolation str 'linear' Src/rec interpolation ('linear' or 'sinc')
interp_order int 2 Interpolation stencil order
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 a Gaussian initial condition
abox bool False Adaptive bounding-box subdomain (restricts computation to active wavefront)
Wdtype dtype np.float32 Wavefield data type (np.float16, np.float32, np.float64)
save_mode str 'all' Elastic adjoint save: 'all' (velocity + stress) or 'v' (velocity only, stress recomputed)
opt tuple auto Devito Operator optimisation options, e.g. ('advanced', {'openmp': True})

2. Forward propagation

# Returns a Devito PerformanceSummary (operator timings, GFlops, …)
summary = solver.forward()

# Save wavefields for a subsequent gradient computation
summary = solver.forward(save=True)

forward accepts the same keyword arguments as op.apply, which are forwarded verbatim after _prepare_kw adds the standard fields (dt, physical parameters, etc.). Useful overrides:

kwarg Description
nt Override the number of time steps for this run
dt Override the time step (default: model.critical_dt)
autotune Devito auto-tune mode for this call ('off', 'basic', 'aggressive')
src_field Restrict source injection to a named field
measurement False to skip receiver interpolation
opt Per-call operator optimisation dict (keyed by operator name)

3. Gradient computation

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

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

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

dm = solver.perturbation(param='vp')   # devito Function
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')

perturbation returns the same Function object on repeated calls, so gradients from multiple shots accumulate in-place. Zero it before a new accumulation:

solver.perturbation(param='vp').data[:] = 0

4. Accessing wavefields and acquisition objects

Current wavefields

solver.state is a dict of rolling-buffer TimeFunction objects. Keys are the field names defined by the solver:

Solver Keys
acoustic (2nd order) u
acoustic (1st order) p, v
TTI (Zhang, Fletcher, SA) u, v
elastic isotropic v, tau
elastic TTI v, tau 
u = solver.state['u']        # TimeFunction
p = solver.pressure          # symbolic pressure expression
p_data = solver.pressure_data  # numpy array of the last time slice

Saved wavefields

solver.saved_state is a dict of snapshot TimeFunction objects. Available only after forward(save=True). Keys are '<name>_save':

u_save = solver.saved_state['u_save']   # shape (nt//t_sub, nx, nz)

When compression is set the data is stored compressed; the compression ratio is available as:

cr = float(np.mean(solver.saved_state['u_save'].compression_ratio))

Source and receivers

src = solver.src   # devito RickerSource — 1 point at domain centre
rec = solver.rec   # devito Receiver   — surface array

# After forward():
traces = solver.rec.data          # (nt, nrec)

# Reposition the source before building the operator:
solver.src.coordinates.data[0] = [500., 20.]   # x, z in metres
Note

src and rec are cached_property objects. If you reposition them you must do so before the first call to forward, fwd_op, or any method that caches the operator. Clear the cache by deleting the cached attributes (del solver.src) if you need to rebuild with different coordinates.

5. Operator-level access

The Devito Operator objects are built lazily and memoised.

op_fwd = solver.fwd_op(save=True)   # forward Operator
op_jac = solver.jacobian_adj_op(params='vp')   # gradient Operator

Pass a per-operator opt dict to apply different optimisation levels to each pass:

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

6. Available solvers

Registry key Class Physics
iso-acoustic AcousticIsotropic 2nd-order (or 1st-order) isotropic acoustic; optional SLS viscosity
tti-zhang / vti-zhang ZhangTTI Zhang (2011) pseudo-acoustic TTI/VTI; optional SLS viscosity
tti-fletcher / vti-fletcher FletcherTTI Fletcher–Du–Fowler (2009) TTI/VTI
tti-selfadjoint / vti-selfadjoint SelfAdjointTTI Self-adjoint TTI/VTI (Bube et al.)
iso-elastic ElasticIsotropic Isotropic elastic velocity–stress (1st order)
tti-elastic / vti-elastic ElasticTTI Anisotropic (TTI/VTI) elastic

vti-* keys are aliases for their tti-* counterparts; the VTI regime is selected automatically when theta (and phi) are zero or absent in the model.

Back to top