```
import numpy as np
import sympy as sp
from devito import Grid, TimeFunction
# Create our grid (computational domain)
= 10
Lx = Lx
Ly = 11
Nx = Nx
Ny = Lx/(Nx-1)
dx = dx
dy = Grid(shape=(Nx,Ny), extent=(Lx,Ly))
grid
# Define u(x,y,t) on this grid
= TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
u
# Define symbol for laplacian replacement
= sp.symbols('H') H
```

# 07 - Custom finite difference coefficients in Devito

## Introduction

When taking the numerical derivative of a function in Devito, the default behaviour is for ‘standard’ finite difference weights (obtained via a Taylor series expansion about the point of differentiation) to be applied. Consider the following example for some field \(u(\mathbf{x},t)\), where \(\mathbf{x}=(x,y)\). Let us define a computational domain/grid and differentiate our field with respect to \(x\).

Now, lets look at the output of \(\partial u/\partial x\):

`print(u.dx.evaluate)`

`-u(t, x, y)/h_x + u(t, x + h_x, y)/h_x`

By default the ‘standard’ Taylor series expansion result, where `h_x`

represents the \(x\)-direction grid spacing, is returned. However, there may be instances when a user wishes to use ‘non-standard’ weights when, for example, implementing a dispersion-relation-preserving (DRP) scheme. See e.g.

[1] Christopher K.W. Tam, Jay C. Webb (1993). ”Dispersion-Relation-Preserving Finite Difference Schemes for Computational Acoustics.” **J. Comput. Phys.**, 107(2), 262–281. https://doi.org/10.1006/jcph.1993.1142

for further details. The use of such modified weights is facilitated in Devito via the ‘symbolic’ finite difference coefficents functionality. Let us start by re-defining the function \(u(\mathbf{x},t)\) in the following manner:

`= TimeFunction(name='u', grid=grid, time_order=2, space_order=2, coefficients='symbolic') u `

Note the addition of the `coefficients='symbolic'`

keyword. Now, when printing \(\partial u/\partial x\) we obtain:

`print(u.dx.evaluate)`

`u(t, x, y)*W(x, 1, u(t, x, y), x) + u(t, x + h_x, y)*W(x + h_x, 1, u(t, x, y), x)`

Owing to the addition of the `coefficients='symbolic'`

keyword the weights have been replaced by sympy functions. Now, take for example the weight `W(x - h_x, 1, u(t, x, y), x)`

, the notation is as follows: * The first `x - h_x`

refers to the spatial location of the weight w.r.t. the evaluation point `x`

. * The `1`

refers to the order of the derivative. * `u(t, x, y)`

refers to the function with which the weight is associated. * Finally, the `x`

refers to the dimension along which the derivative is being taken.

Symbolic coefficients can then be manipulated using the Devito ‘Coefficient’ and ‘Substitutions’ objects. First, let us consider an example where we wish to replace the coefficients with a set of constants throughout the entire computational domain.

```
from devito import Coefficient, Substitutions # Import the Devito Coefficient and Substitutions objects
# Grab the grid spatial dimensions: Note x[0] will correspond to the x-direction and x[1] to y-direction
= grid.dimensions
x # Form a Coefficient object and then a replacement rules object (to pass to a Devito equation):
= Coefficient(1, u, x[0], np.array([-0.6, 0.1, 0.6]))
u_x_coeffs = Substitutions(u_x_coeffs) coeffs
```

Devito Coefficient ojects take arguments in the following order: 1. Derivative order (in the above example this is the first derivative) 2. Function to which the coefficients ‘belong’ (in the above example this is the time function `u`

) 3. Dimension on which coefficients will be applied (in the above example this is the x-direction) 4. Coefficient data. Since, in the above example, the coefficients have been applied as a 1-d numpy array replacement will occur at the equation level. (Note that other options are in development and will be the subject of future notebooks).

Now, lets form a Devito equation, pass it the Substitutions object, and take a look at the output:

```
from devito import Eq
= Eq(u.dt+u.dx, coefficients=coeffs)
eq print(eq.evaluate)
```

`Eq(-u(t, x, y)/dt + u(t + dt, x, y)/dt + 0.1*u(t, x, y) + 0.6*u(t, x + h_x, y), 0)`

We see that in the above equation the standard weights for the first derivative of `u`

in the \(x\)-direction have now been replaced with our user defined weights. Note that since no replacement rules were defined for the time derivative (`u.dt`

) standard weights have replaced the symbolic weights.

Now, let us consider a more complete example.

## Example: Finite difference modeling for a large velocity-contrast acousitc wave model

It is advised to read through the ‘Introduction to seismic modelling’ notebook located in devito/examples/seismic/tutorials/01_modelling.ipynb before proceeding with this example since much introductory material will be ommited here. The example now considered is based on an example introduced in

[2] Yang Liu (2013). ”Globally optimal finite-difference schemes based on least squares.” **GEOPHYSICS**, 78(4), 113–132. https://doi.org/10.1190/geo2012-0480.1.

See figure 18 of [2] for further details. Note that here we will simply use Devito to ‘reproduce’ the simulations leading to two results presented in the aforementioned figure. No analysis of the results will be carried out. The domain under consideration has a sptaial extent of \(2km \times 2km\) and, letting \(x\) be the horizontal coordinate and \(z\) the depth, a velocity profile such that \(v_1(x,z)=1500ms^{-1}\) for \(z\leq1200m\) and \(v_2(x,z)=4000ms^{-1}\) for \(z>1200m\).

```
#NBVAL_IGNORE_OUTPUT
from examples.seismic import Model, plot_velocity
%matplotlib inline
# Define a physical size
= 2000
Lx = Lx
Lz = 10
h = int(Lx/h)+1
Nx = Nx
Nz
= (Nx, Nz) # Number of grid point
shape = (h, h) # Grid spacing in m. The domain size is now 2km by 2km
spacing = (0., 0.)
origin
# Our scheme will be 10th order (or less) in space.
= 10
order
# Define a velocity profile. The velocity is in km/s
= np.empty(shape, dtype=np.float32)
v 121] = 1.5
v[:, :121:] = 4.0
v[:,
# With the velocity and model size defined, we can create the seismic model that
# encapsulates these properties. We also define the size of the absorbing layer as 10 grid points
= 10
nbl = Model(vp=v, origin=origin, shape=shape, spacing=spacing,
model =20, nbl=nbl, bcs="damp")
space_order
plot_velocity(model)
```

`Operator `initdamp` ran in 0.01 s`

The seismic wave source term will be modelled as a Ricker Wavelet with a peak-frequency of \(25\)Hz located at \((1000m,800m)\). Before applying the DRP scheme, we begin by generating a ‘reference’ solution using a spatially high-order standard finite difference scheme and time step well below the model’s critical time-step. The scheme will be 2nd order in time.

```
from examples.seismic import TimeAxis
= 0. # Simulation starts a t=0
t0 = 500. # Simulation lasts 0.5 seconds (500 ms)
tn = 1.0 # Time step of 0.2ms
dt
= TimeAxis(start=t0, stop=tn, step=dt)
time_range
#NBVAL_IGNORE_OUTPUT
from examples.seismic import RickerSource
= 0.025 # Source peak frequency is 25Hz (0.025 kHz)
f0 = RickerSource(name='src', grid=model.grid, f0=f0,
src =1, time_range=time_range)
npoint
# First, position source centrally in all dimensions, then set depth
0, :] = np.array(model.domain_size) * .5
src.coordinates.data[0, -1] = 800. # Depth is 800m
src.coordinates.data[
# We can plot the time signature to see the wavelet
src.show()
```

Now let us define our wavefield and PDE:

```
# Define the wavefield with the size of the model and the time dimension
= TimeFunction(name="u", grid=model.grid, time_order=2, space_order=order)
u
# We can now write the PDE
= model.m * u.dt2 - H + model.damp * u.dt
pde
# This discrete PDE can be solved in a time-marching way updating u(t+dt) from the previous time step
# Devito as a shortcut for u(t+dt) which is u.forward. We can then rewrite the PDE as
# a time marching updating equation known as a stencil using customized SymPy functions
from devito import solve
= Eq(u.forward, solve(pde, u.forward).subs({H: u.laplace})) stencil
```

```
# Finally we define the source injection and receiver read function to generate the corresponding code
= src.inject(field=u.forward, expr=src * dt**2 / model.m) src_term
```

Now, lets create the operator and execute the time marching scheme:

```
from devito import Operator
= Operator([stencil] + src_term, subs=model.spacing_map) op
```

```
#NBVAL_IGNORE_OUTPUT
=time_range.num-1, dt=dt) op(time
```

`Operator `Kernel` ran in 0.03 s`

```
PerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.02013399999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.004232999999999992, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
```

And plot the result:

```
# NBVAL_IGNORE_OUTPUT
#import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
= 2000
Lx = 2000
Lz
= nbl*h
abs_lay
= h
dx = dx
dz = np.mgrid[-abs_lay: Lx+abs_lay+1e-10: dx, -abs_lay: Lz+abs_lay+1e-10: dz]
X, Z
= 100
levels = 5
clip
= plt.figure(figsize=(14, 7))
fig = fig.add_subplot(111)
ax1 = ax1.imshow(u.data[0,:,:].T, vmin=-clip, vmax=clip, cmap=cm.seismic, extent=[0, Lx, 0, Lz])
cont
fig.colorbar(cont)'$x$')
ax1.set_xlabel('$z$')
ax1.set_ylabel('$u(x,z,500)$')
ax1.set_title(
plt.gca().invert_yaxis() plt.show()
```

We will now reimplement the above model applying the DRP scheme presented in [2].

First, since we wish to apply different custom FD coefficients in the upper on lower layers we need to define these two ‘subdomains’ using the `Devito SubDomain`

functionality:

```
from devito import SubDomain
# Define our 'upper' and 'lower' SubDomains:
class Upper(SubDomain):
= 'upper'
name def define(self, dimensions):
= dimensions
x, z # We want our upper layer to span the entire x-dimension and all
# but the bottom 80 (+boundary layer) cells in the z-direction, which is achieved via
# the following notation:
return {x: x, z: ('left', 80+nbl)}
class Lower(SubDomain):
= 'lower'
name def define(self, dimensions):
= dimensions
x, z # We want our lower layer to span the entire x-dimension and all
# but the top 121 (+boundary layer) cells in the z-direction.
return {x: x, z: ('right', 121+nbl)}
# Create these subdomains:
= Upper()
ur = Lower() lr
```

We now create our model incoporating these subdomains:

```
#NBVAL_IGNORE_OUTPUT
# Create our model passing it our 'upper' and 'lower' subdomains:
= Model(vp=v, origin=origin, shape=shape, spacing=spacing,
model =order, nbl=nbl, subdomains=(ur,lr), bcs="damp") space_order
```

`Operator `initdamp` ran in 0.01 s`

And re-define model related objects. Note that now our wave-field will be defined with `coefficients='symbolic'`

.

```
= 0. # Simulation starts a t=0
t0 = 500. # Simulation last 1 second (500 ms)
tn = 1.0 # Time step of 1.0ms
dt
= TimeAxis(start=t0, stop=tn, step=dt)
time_range
= 0.025 # Source peak frequency is 25Hz (0.025 kHz)
f0 = RickerSource(name='src', grid=model.grid, f0=f0,
src =1, time_range=time_range)
npoint
0, :] = np.array(model.domain_size) * .5
src.coordinates.data[0, -1] = 800. # Depth is 800m
src.coordinates.data[
# New wave-field
= TimeFunction(name="u_DRP", grid=model.grid, time_order=2, space_order=order, coefficients='symbolic') u_DRP
```

We now create a stencil for each of our ‘Upper’ and ‘Lower’ subdomains defining different custom FD weights within each of these subdomains.

```
# The underlying pde is the same in both subdomains
= model.m * u_DRP.dt2 - H + model.damp * u_DRP.dt
pde_DRP
# Define our custom FD coefficients:
= model.grid.dimensions
x, z # Upper layer
= np.array([ 2.00462e-03, -1.63274e-02, 7.72781e-02,
weights_u -3.15476e-01, 1.77768e+00, -3.05033e+00,
1.77768e+00, -3.15476e-01, 7.72781e-02,
-1.63274e-02, 2.00462e-03])
# Lower layer
= np.array([ 0. , 0. , 0.0274017,
weights_l -0.223818, 1.64875 , -2.90467,
1.64875 , -0.223818, 0.0274017,
0. , 0. ])
# Create the Devito Coefficient objects:
= Coefficient(2, u_DRP, x, weights_u/x.spacing**2)
ux_u_coeffs = Coefficient(2, u_DRP, z, weights_u/z.spacing**2)
uz_u_coeffs = Coefficient(2, u_DRP, x, weights_l/x.spacing**2)
ux_l_coeffs = Coefficient(2, u_DRP, z, weights_l/z.spacing**2)
uz_l_coeffs # And the replacement rules:
= Substitutions(ux_u_coeffs,uz_u_coeffs)
coeffs_u = Substitutions(ux_l_coeffs,uz_l_coeffs)
coeffs_l # Create a stencil for each subdomain:
= Eq(u_DRP.forward, solve(pde_DRP, u_DRP.forward).subs({H: u_DRP.laplace}),
stencil_u = model.grid.subdomains['upper'], coefficients=coeffs_u)
subdomain = Eq(u_DRP.forward, solve(pde_DRP, u_DRP.forward).subs({H: u_DRP.laplace}),
stencil_l = model.grid.subdomains['lower'], coefficients=coeffs_l)
subdomain
# Source term:
= src.inject(field=u_DRP.forward, expr=src * dt**2 / model.m)
src_term
# Create the operator, incoporating both upper and lower stencils:
= Operator([stencil_u, stencil_l] + src_term, subs=model.spacing_map) op
```

And now execute the operator:

```
#NBVAL_IGNORE_OUTPUT
=time_range.num-1, dt=dt) op(time
```

`Operator `Kernel` ran in 0.03 s`

```
PerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=0.01799599999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.0044569999999999966, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
```

And plot the new results:

```
# NBVAL_IGNORE_OUTPUT
= plt.figure(figsize=(14, 7))
fig = fig.add_subplot(111)
ax1 = ax1.imshow(u_DRP.data[0,:,:].T, vmin=-clip, vmax=clip, cmap=cm.seismic, extent=[0, Lx, 0, Lz])
cont
fig.colorbar(cont)0, Lx, 0, Lz])
ax1.axis(['$x$')
ax1.set_xlabel('$z$')
ax1.set_ylabel('$u_{DRP}(x,z,500)$')
ax1.set_title(
plt.gca().invert_yaxis() plt.show()
```

Finally, for comparison, lets plot the difference between the standard 20th order and optimized 10th order models:

```
# NBVAL_IGNORE_OUTPUT
= plt.figure(figsize=(14, 7))
fig = fig.add_subplot(111)
ax1 = ax1.imshow(u_DRP.data[0,:,:].T-u.data[0,:,:].T, vmin=-clip, vmax=clip, cmap=cm.seismic, extent=[0, Lx, 0, Lz])
cont
fig.colorbar(cont)'$x$')
ax1.set_xlabel('$z$')
ax1.set_ylabel(
plt.gca().invert_yaxis() plt.show()
```

```
#NBVAL_IGNORE_OUTPUT
# Wavefield norm checks
assert np.isclose(np.linalg.norm(u.data[-1]), 82.170, atol=0, rtol=1e-4)
assert np.isclose(np.linalg.norm(u_DRP.data[-1]), 83.624, atol=0, rtol=1e-4)
```