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 (km/s), as an array or read on demand
`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.

5. Models read on demand

Any physical parameter can be given as a file read on demand — a memory mapped array, an HDF5 dataset, a SEG-Y reader — instead of an array. With MPI each rank then reads only the window covering the slab of the grid it owns, so the model is never read in full nor held in memory. See MPI.