Model

The WaveModel class (recipes/model.py) extends Devito’s SeismicModel with production-oriented features: mixed-precision physical parameters, lazy absorbing-boundary initialisation, water/earth layer subdomains, and immersed free-surface support via signed distance functions.

1. Built-in construction

layered_model — quick stratified model

from recipes import layered_model

model = layered_model(
    'iso-acoustic',          # physics string — controls which parameters are built
    shape=(101, 101),
    spacing=(10.0, 10.0),   # meters
    nbl=40,
    space_order=8,
    nlayers=4,
)

The physics string determines which secondary parameters are generated:

String contains Extra parameters created
tti or vti epsilon, delta, theta (and phi in 3-D)
elastic vs → converts to lam and mu
visco invqp (uniform 1/Q)

Extra keyword arguments forwarded to WaveModel:

Argument Default Description
density False Derive buoyancy b = 1/ρ from Gardner’s relation
water_layer False Add an isotropic acoustic water layer subdomain at the top
smooth False Gaussian-smooth Thomsen parameters after construction
Pdtype np.float32 Storage precision for physical parameters

WaveModel — direct construction

import numpy as np
from recipes.model import WaveModel

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

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

Constructor arguments

Argument Required Description
vp always P-wave velocity array (km/s)
shape always Grid dimensions without boundary padding
spacing always Grid spacing (m)
space_order always Finite-difference stencil order
nbl always Absorbing boundary layer width (cells)
origin always Physical coordinates of the first grid point
fs no Conventional flat free surface at the top
dtype no Grid floating-point type (default np.float32)
Pdtype no Physical parameter storage precision (default np.float32)
b no Buoyancy 1/ρ; scalar 1 = constant density
epsilon, delta, theta, phi no Thomsen anisotropy parameters
vs no S-wave velocity → lam and mu computed automatically
invqp no Inverse Q factor for viscoacoustic attenuation
sdf no Signed distance function for immersed free-surface topography
interp_order no Source/receiver interpolation order (default 2)
water_layer no Build water/earth subdomain partition

2. Using your own model

The solvers only require the model to expose a fixed set of attributes and methods. You can bring your own class as long as it satisfies this interface — you do not have to subclass WaveModel.

2.1 Mandatory attributes (all solvers)

These are accessed by Propagator unconditionally:

Attribute Type Description
grid devito.Grid Devito grid shared by all Functions
spacing tuple of float Physical grid spacing (m) per dimension
shape tuple of int Grid shape without boundary padding
domain_size tuple of float Physical extent: (shape[i]-1) * spacing[i]
space_order int FD stencil order
nbl int Absorbing layer width
critical_dt float Maximum stable time step (CFL condition)
spacing_map dict Maps symbolic spacings and dt to their numerical values; passed as subs to every Operator
damp Function or scalar Absorbing boundary field ('damp' convention: 0 inside, positive in layer; 'mask' convention: 1 inside, decreasing in layer). Lazily initialised — call model.initialize_damp(abc_type=...) before the first operator application
physical_params(**kwargs) method Returns a dict mapping each physical parameter Function to its data, suitable for passing to op.apply. Optional kwargs override individual values
fs bool Free-surface flag
has_sdf bool Whether an immersed boundary (SDF) is present

2.2 Physics-specific parameters

Each solver reads specific attributes from the model. The table below lists which attribute each physics requires. All of them must be either a devito.Function (array-valued) or a devito.Constant (scalar).

Attribute Solvers that use it Description
m acoustic, TTI Squared slowness 1/vp²
vp acoustic, TTI, elastic P-wave velocity (km/s)
b all Buoyancy 1/ρ; constant 1 for homogeneous density
lam elastic, elastic-TTI Lamé parameter λ
mu elastic, elastic-TTI Lamé parameter μ
vs elastic, elastic-TTI S-wave velocity (km/s)
mu0 elastic, elastic-TTI Scalar value of μ at the first grid point (used to detect fluid layers)
epsilon TTI, VTI Thomsen ε
delta TTI, VTI Thomsen δ
theta TTI (3-D and 2-D) Tilt angle (radians)
phi TTI 3-D only Azimuth angle (radians)
eps Fletcher TTI Alias for epsilon
sig Fletcher TTI Symmetry parameter σ (typically 0.5)
invqp viscoacoustic Inverse Q factor 1/Qp
is_visco acoustic, TTI Boolean — True if invqp is non-zero

2.3 Optional subdomain attributes

The following Grid subdomains are looked up by name at runtime. Provide them only when the corresponding feature is needed.

Subdomain name Accessed via Purpose
'coredomain' model.grid.subdomains Interior below an immersed free surface (SDF)
'boundarydomain' model.grid.subdomains Region near the immersed surface where modified stencils apply
'waterlayer' model.grid.subdomains Isotropic acoustic top region
'earthlayer' model.grid.subdomains Anisotropic / elastic region below the water

2.4 Immersed free surface

When has_sdf is True the solver also calls:

model._immersed_fs(boundary_subdomain, bcs, *stencil_groups)

This method returns modified stencil groups with boundary-corrected FD coefficients near the surface. Implement it only if you need immersed free-surface support; for flat free surfaces use model.fs = True and the mirror / mirror_source utilities instead.

2.5 Minimal custom model example

import numpy as np
from devito import Grid, Function, Constant
from devito.builtins import initialize_function


class MyModel:
    """
    Minimal acoustic model compatible with AcousticIsotropic.
    """

    def __init__(self, shape, spacing, vp, nbl=40, space_order=8):
        self.shape    = shape
        self.spacing  = spacing
        self.nbl      = nbl
        self.space_order = space_order
        self.fs       = False
        self.has_sdf  = False

        pad = [(nbl, nbl)] * (len(shape) - 1) + [(nbl, nbl)]
        padded = tuple(n + l + r for n, (l, r) in zip(shape, pad))
        self.grid = Grid(shape=padded,
                         extent=tuple((n-1)*d for n, d in zip(padded, spacing)),
                         dtype=np.float32)

        # P-wave velocity and derived quantities
        so = (space_order, space_order//2, space_order//2)
        self.vp  = Function(name='vp',  grid=self.grid, space_order=so)
        self.m   = Function(name='m',   grid=self.grid, space_order=so)
        self.b   = Constant(name='b',   value=1.0)
        self.damp = Function(name='damp', grid=self.grid, space_order=so)

        initialize_function(self.vp, vp, pad)
        self.m.data[:] = 1.0 / self.vp.data**2

    @property
    def domain_size(self):
        return tuple((n - 1) * d for n, d in zip(self.shape, self.spacing))

    @property
    def critical_dt(self):
        # Simple CFL estimate
        return 0.9 * min(self.spacing) / (np.sqrt(self.grid.dim) * self._vp_max)

    @property
    def spacing_map(self):
        mapper = dict(self.grid.spacing_map)
        mapper[self.grid.time_dim.spacing] = self.critical_dt
        return mapper

    def initialize_damp(self, abc_type='damp'):
        # Populate self.damp with your absorbing profile
        ...

    def physical_params(self, **kwargs):
        params = {self.vp: self.vp.data,
                  self.m:  self.m.data,
                  self.b:  float(self.b.data)}
        params.update(kwargs)
        return params

3. Absorbing boundaries

WaveModel initialises the damping field lazily: damp is allocated in __init__ but filled only the first time an operator is applied, via initialize_damp(abc_type). Two conventions are supported:

Convention _abc value on solver damp profile
Additive 'damp' 0 inside the domain, grows toward the edges
Multiplicative mask 'mask' 1 inside the domain, decreases toward the edges

The Propagator.run_propagation wrapper checks which convention the solver declares (_abc) and re-initialises damp automatically if the stored data does not match.

To skip absorbing boundaries entirely (e.g. for a pure kernel benchmark), pass nbl=0 or set model.damp = 0 (additive) / model.damp = 1 (mask).

4. Physical parameter precision

WaveModel stores parameters at the precision specified by Pdtype (default np.float32). When Pdtype=np.float16, DevitoPRO quantization is applied instead of raw float16 to maintain accuracy over the narrow value range typical of seismic parameters:

model = WaveModel(..., Pdtype=np.float16)

Wavefields use a separate precision controlled by the solver option Wdtype.

Back to top