Acoustic Wave Modelling using Devito

Consider the acoustic wave equation given in 3D:

\[m(x,y,z)\frac{\partial^2 u}{\partial t^2} + \eta(x,y,z)\frac{\partial u}{\partial t}-\frac{\partial^2 u}{\partial x^2}-\frac{\partial^2 u}{\partial y^2}-\frac{\partial^2 u}{\partial z^2}= q\]

where Damp(\(\eta\)) is the dampening coefficient for absorbing boundary condition, \(m=\frac{1}{v(x,y,z)^2}\), \(v\) is the velocity , \(u\) is pressure field and \(q\) is the pressure source term.

First of all, we will set up seismic datas,

In this tutorial, Model, an instance of IGrid(), stores the origin, the position in meters of the position of index (0,0,0), as (0,0,0) in 3D, spacing, the size of the discrete grid, distance in meters between two consecutive grid points in each direction, dimensions, the size of the model in number of grid points and true velocity. The time stepping rate, dt, is derived from model.get_critical_dt().

Data, an instance of IShot() stores amplitudes of the source at each time step, source coordinates and receiver coordinates.:

from containers import IGrid, IShot

model = IGrid()
dimensions = (50, 50, 50)
model.shape = dimensions
origin = (0., 0., 0.)
spacing = (20., 20., 20.)

# True velocity
true_vp = np.ones(dimensions) + 2.0
true_vp[:, :, int(dimensions[0] / 2):int(dimensions[0])] = 4.5

model.create_model(origin, spacing, true_vp)

# Define seismic data.
data = IShot()

f0 = .010
dt = model.get_critical_dt()
t0 = 0.0
tn = 250.0
nt = int(1+(tn-t0)/dt)
h = model.get_spacing()


# Set up the source as Ricker wavelet for f0
def source(t, f0):
    r = (np.pi * f0 * (t - 1./f0))

    return (1-2.*r**2)*np.exp(-r**2)


time_series = source(np.linspace(t0, tn, nt), f0)
location = (origin[0] + dimensions[0] * spacing[0] * 0.5, 500,
                origin[1] + 2 * spacing[1])
data.set_source(time_series, dt, location)

receiver_coords = np.zeros((101, 3))
receiver_coords[:, 0] = np.linspace(50, 950, num=101)
receiver_coords[:, 1] = 500
receiver_coords[:, 2] = location[2]
data.set_receiver_pos(receiver_coords)
data.set_shape(nt, 101)

Then we set up the dampening coefficient.

from devito.interfaces import DenseData
self.damp = DenseData(name="damp", shape=self.model.get_shape_comp(),
                      dtype=self.dtype)
# Initialize damp by calling the function that can precompute damping
damp_boundary(self.damp.data)

DenseData is a devito data object used to store and manage spatially varying data.

damp_boundary() function initialises the damp on each grid point. The dampening values in data are stored in self.damp.data.

Then we will set up the source

srccoord = np.array(self.data.source_coords, dtype=self.dtype)[np.newaxis, :]
self.src = SourceLike(name="src", npoint=1, nt=data.traces.shape[1],
                      dt=self.dt, h=self.model.get_spacing(),
                      coordinates=srccoord, ndim=len(self.damp.shape),
                      dtype=self.dtype, nbpml=nbpml)
self.src.data[:] = data.get_source()[:, np.newaxis]

SourceLike is an object inheriting from PointData, a Devito data object for sparse point data as a Function symbol.

We initialize the source to at coordinates : srccoord and set its data to be data.get_source(). Receivers are initialized in a similar way

rec = SourceLike(name="rec", npoint=nrec, nt=nt, dt=dt, h=model.get_spacing(),
                     coordinates=data.receiver_coords, ndim=len(damp.shape),
                     dtype=damp.dtype, nbpml=model.nbpml)

Then, We set up the initial condition and allocate the grid for m. Value of m on each grid point is stored as a numpy array in m.data[:].

m = DenseData(name="m", shape=model.get_shape_comp(), dtype=damp.dtype)
m.data[:] = model.padm()

after that, we will initialize u

u = TimeData(name="u", shape=model.get_shape_comp(), time_dim=nt,
                 time_order=time_order, space_order=spc_order, save=True,
                 dtype=damp.dtype)

TimeData is a Devito data object used to store and manage time-varying as well as space-varying data

We initialize our grid to be of size model.get_shape_comp() which is a 3-D tuple. time_dim represents the size of the time dimension that dictates the leading dimension of the data buffer. time_order and space_order represent the discretization order for time and space respectively.

The next step is to generate the stencil to be solved by a devito.operator.Operator

The stencil is created according to Devito conventions. It uses a Sympy expression to represent the acoustic wave equation. Devito makes it easy to represent the equation in a finite-difference form by providing properties dt2, laplace, dt. We then generate the stencil by solving eqn for u.forward = u(t+dt,x,y,z), a symbol for the time-forward state of the function.

from devito import Operator
from sympy import Eq, solve, symbols
eqn = m*u.dt2-u.laplace+damp*u.dt
stencil = solve(eqn, u.forward)[0]

We plug the stencil in an Operator, as shown, as well as the the values of the spacing between cells h, the temporal spaces s, the number of timesteps nt, the spc_border which are the ghost points at the edge of the domain to touse the same stencil everywhere and so on. The forward argument represents propagation time direction. It’s set to be True if forwarding in time and False if backwarding in time.

s, h = symbols('s h')
subs = {s: model.get_critical_dt(), h: model.get_spacing()}
super(ForwardOperator, self).__init__(nt, m.shape,
                                            stencils=Eq(u.forward, stencil),
                                            subs=subs,
                                            spc_border=spc_order/2,
                                            time_order=time_order,
                                            forward=True,
                                            dtype=m.dtype,
                                            **kwargs)

After that, we will insert source and receiver terms into the input parameters to generate the C++ file that contains required inputs. For the output, we will add receivers so that it will output \(u\) on each receiver coordinate on all time steps. src.add(m, u) and red.read(u) will generate C iteration codes over points and they will be added into stencils in C++ file.

self.input_params += [src, src.coordinates, rec, rec.coordinates]
self.output_params += [rec]
self.propagator.time_loop_stencils_a = src.add(m, u) + rec.read(u)
self.propagator.add_devito_param(src)
self.propagator.add_devito_param(src.coordinates)
self.propagator.add_devito_param(rec)
self.propagator.add_devito_param(rec.coordinates)

To execute the generated Operator, we simply call apply(). As mentioned, it will output \(u\) on each receiver coordinates and \(u\) on each grid points for all time steps. The results can be found in u.data.

For a complete example of this code, check file acoustic_example.py in the examples folder.