This operator is based on simplfications of the systems presented in:
Self-adjoint, energy-conserving second-order pseudoacoustic systems for VTI and TTI media for reverse time migration and full-waveform inversion (2016)
Kenneth Bube, John Washbourne, Raymond Ergas, and Tamas Nemeth
SEG Technical Program Expanded Abstracts
https://library.seg.org/doi/10.1190/segam2016-13878451.1
The goal of this tutorial set is to generate and prove correctness of modeling and inversion capability in Devito for variable density visco- acoustics using an energy conserving form of the wave equation. We describe how the linearization of the energy conserving skew self adjoint system with respect to modeling parameters allows using the same modeling system for all nonlinear and linearized forward and adjoint finite difference evolutions. There are three notebooks in this series:
forward and adjoint modeling operations.There are similar series of notebooks implementing and testing operators for VTI and TTI anisotropy (README.md).
Below we introduce the skew self adjoint form of the scalar isotropic variable density visco- acoustic wave equation with a simple form of dissipation only Q attenuation. This dissipation only (no dispersion) attenuation term $\left (\frac{\displaystyle \omega}{Q}\ \partial_t\ u \right)$ is an approximation of a Maxwell Body -- that is to say viscoelasticity approximated with a spring and dashpot in series. In practice this approach for attenuating outgoing waves is very similar to the Cerjan style damping in absorbing boundaries used elsewhere in Devito (References).
The derivation of the attenuation model is not in scope for this tutorial, but one important point is that the physics in the absorbing boundary region and the interior of the model are unified, allowing the same modeling equations to be used everywhere, with physical Q values in the interior tapering to non-physical small Q at the boundaries to attenuate outgoing waves.
| Symbol | Description | Dimensionality |
|---|---|---|
| $\omega_c = 2 \pi f_c$ | center angular frequency | constant |
| $m(x,y,z)$ | P wave velocity | function of space |
| $b(x,y,z)$ | buoyancy $(1 / \rho)$ | function of space |
| $Q(x,y,z)$ | Attenuation at frequency $\omega_c$ | function of space |
| $u(t,x,y,z)$ | Pressure wavefield | function of time and space |
| $q(t,x,y,z)$ | Source term | function of time, localized in space |
| $\overleftarrow{\partial_t}$ | shifted first derivative wrt $t$ | shifted 1/2 sample backward in time |
| $\partial_{tt}$ | centered second derivative wrt $t$ | centered in time |
| $\overrightarrow{\partial_x},\ \overrightarrow{\partial_y},\ \overrightarrow{\partial_z}$ | + shifted first derivative wrt $x,y,z$ | shifted 1/2 sample forward in space |
| $\overleftarrow{\partial_x},\ \overleftarrow{\partial_y},\ \overleftarrow{\partial_z}$ | - shifted first derivative wrt $x,y,z$ | shifted 1/2 sample backward in space |
| $\Delta_t, \Delta_x, \Delta_y, \Delta_z$ | sampling rates for $t, x, y , z$ | $t, x, y , z$ |
We use the arrow symbols over derivatives $\overrightarrow{\partial_x}$ as a shorthand notation to indicate that the derivative is taken at a shifted location. For example:
$\overrightarrow{\partial_x}\ u(t,x,y,z)$ indicates that the $x$ derivative of $u(t,x,y,z)$ is taken at $u(t,x+\frac{\Delta x}{2},y,z)$.
$\overleftarrow{\partial_z}\ u(t,x,y,z)$ indicates that the $z$ derivative of $u(t,x,y,z)$ is taken at $u(t,x,y,z-\frac{\Delta z}{2})$.
$\overleftarrow{\partial_t}\ u(t,x,y,z)$ indicates that the $t$ derivative of $u(t,x,y,z)$ is taken at $u(t-\frac{\Delta_t}{2},x,y,z)$.
We usually drop the $(t,x,y,z)$ notation from wavefield variables unless required for clarity of exposition, so that $u(t,x,y,z)$ becomes $u$.
Our skew self adjoint wave equation is written:
$$ \frac{b}{m^2} \left( \frac{\omega_c}{Q} \overleftarrow{\partial_t}\ u + \partial_{tt}\ u \right) = \overleftarrow{\partial_x}\left(b\ \overrightarrow{\partial_x}\ u \right) + \overleftarrow{\partial_y}\left(b\ \overrightarrow{\partial_y}\ u \right) + \overleftarrow{\partial_z}\left(b\ \overrightarrow{\partial_z}\ u \right) + q $$An advantage of this form is that the same system can be used to provide stable modes of propagation for all operations needed in quasi- Newton optimization:
This advantage is more important for anisotropic operators, where widely utilized non energy conserving formulations can provide unstable adjoints and thus unstable gradients for anisotropy parameters.
The skew self adjoint formulation is evident in the shifted spatial derivatives, with the derivative on the right side $\overrightarrow{\partial}$ shifting forward in space one-half cell, and the derivative on the left side $\overleftarrow{\partial}$ shifting backward in space one-half cell.
$\overrightarrow{\partial}$ and $\overleftarrow{\partial}$ are anti-symmetric (also known as skew symmetric), meaning that for two random vectors $x_1$ and $x_2$, correctly implemented numerical derivatives will have the following property:
$$ x_2 \cdot \left( \overrightarrow{\partial_x}\ x_1 \right) \approx -\ x_1 \cdot \left( \overleftarrow{\partial_x}\ x_2 \right) $$Below we will demonstrate this skew symmetry with a simple unit test on Devito generated derivatives.
In the following notebooks in this series, material parameters sandwiched between the derivatives -- including anisotropy parameters -- become much more interesting, but here buoyancy $b$ is the only material parameter between derivatives in our skew self adjoint (SSA) wave equation.
In [1]:
import numpy as np
from examples.seismic import RickerSource, Receiver, TimeAxis
from devito import (Grid, Function, TimeFunction, SpaceDimension, Constant,
Eq, Operator, solve, configuration)
from devito.finite_differences import Derivative
from devito.builtins import gaussian_smooth
from examples.seismic.skew_self_adjoint import compute_critical_dt, setup_w_over_q
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
from timeit import default_timer as timer
# These lines force images to be displayed in the notebook, and scale up fonts
%matplotlib inline
mpl.rc('font', size=14)
# Make white background for plots, not transparent
plt.rcParams['figure.facecolor'] = 'white'
# Set the default language to openmp
configuration['language'] = 'openmp'
# Set logging to debug, captures statistics on the performance of operators
configuration['log-level'] = 'DEBUG'
As noted above, we prove with a small 1D unit test and 8th order spatial operator that the Devito shifted first derivatives are skew symmetric. This anti-symmetry can be demonstrated for the forward and backward half cell shift first derivative operators $\overrightarrow{\partial}$ and $\overleftarrow{\partial}$ with two random vectors $x_1$ and $x_2$ by verifying the dot product test as written above.
We will use Devito to implement the following two equations with an Operator:
And then verify the dot products are equivalent:
$$ f_1 \cdot g_2 \approx g_1 \cdot f_2 $$We use the following test for relative error (note the flipped signs in numerator and denominator due to anti- symmetry):
$$ \frac{\displaystyle f_1 \cdot g_2 + g_1 \cdot f_2} {\displaystyle f_1 \cdot g_2 - g_1 \cdot f_2}\ <\ \epsilon $$
In [2]:
#NBVAL_INGNORE_OUTPUT
# Make 1D grid to test derivatives
n = 101
d = 1.0
shape = (n, )
spacing = (1 / (n-1), )
origin = (0., )
extent = (d * (n-1), )
dtype = np.float64
# Initialize Devito grid and Functions for input(f1,g1) and output(f2,g2)
# Note that space_order=8 allows us to use an 8th order finite difference
# operator by properly setting up grid accesses with halo cells
grid1d = Grid(shape=shape, extent=extent, origin=origin, dtype=dtype)
x = grid1d.dimensions[0]
f1 = Function(name='f1', grid=grid1d, space_order=8)
f2 = Function(name='f2', grid=grid1d, space_order=8)
g1 = Function(name='g1', grid=grid1d, space_order=8)
g2 = Function(name='g2', grid=grid1d, space_order=8)
# Fill f1 and g1 with random values in [-1,+1]
f1.data[:] = -1 + 2 * np.random.rand(n,)
g1.data[:] = -1 + 2 * np.random.rand(n,)
# Equation defining: [f2 = forward 1/2 cell shift derivative applied to f1]
equation_f2 = Eq(f2, f1.dx(x0=x+0.5*x.spacing))
# Equation defining: [g2 = backward 1/2 cell shift derivative applied to g1]
equation_g2 = Eq(g2, g1.dx(x0=x-0.5*x.spacing))
# Define an Operator to implement these equations and execute
op = Operator([equation_f2, equation_g2])
op()
# Compute the dot products and the relative error
f1g2 = np.dot(f1.data, g2.data)
g1f2 = np.dot(g1.data, f2.data)
diff = (f1g2+g1f2)/(f1g2-g1f2)
tol = 100 * np.finfo(np.float32).eps
print("f1g2, g1f2, diff, tol; %+.6e %+.6e %+.6e %+.6e" % (f1g2, g1f2, diff, tol))
# At last the unit test
# Assert these dot products are float epsilon close in relative error
assert diff < 100 * np.finfo(np.float32).eps
You can inspect the finite difference coefficients and locations for evaluation with the code shown below.
For your reference, the finite difference coefficients seen in the first two stanzas below are exactly the coefficients generated in Table 2 of Fornberg's paper Generation of Finite Difference Formulas on Arbitrarily Spaced Grids linked below (References).
Note that you don't need to inspect the generated code, but this does provide the option to use this highly optimized code in applications that do not need or require python. If you inspect the code you will notice hallmarks of highly optimized c code, including pragmas for vectorization, and decorations for pointer restriction and alignment.
In [3]:
#NBVAL_INGNORE_OUTPUT
# Show the FD coefficients generated by Devito
# for the forward 1/2 cell shifted first derivative operator
print("\n\nForward +1/2 cell shift;")
print("..................................")
print(f1.dx(x0=x+0.5*x.spacing).evaluate)
# Show the FD coefficients generated by Devito
# for the backward 1/2 cell shifted first derivative operator
print("\n\nBackward -1/2 cell shift;")
print("..................................")
print(f1.dx(x0=x-0.5*x.spacing).evaluate)
# Show code generated by Devito for applying the derivatives
print("\n\nGenerated c code;")
print("..................................")
print(op.ccode)
The next step in implementing our Devito modeling operator is to define the equation used to update the pressure wavefield as a function of time. What follows is a bit of algebra using the wave equation and finite difference approximations to time derivatives to express the pressure wavefield forward in time $u(t+\Delta_t)$ as a function of the current $u(t)$ and previous $u(t-\Delta_t)$ pressure wavefields.
The second order accurate centered approximation to the second time derivative involves three wavefields: $u(t-\Delta_t)$, $u(t)$, and $u(t+\Delta_t)$. In order to advance our finite difference solution in time, we solve for $u(t+\Delta_t)$.
$$ \begin{aligned} \partial_{tt}\ u &= \frac{\displaystyle u(t+\Delta_t) - 2\ u(t) + u(t-\Delta_t)}{\displaystyle \Delta_t^2} \\[5pt] u(t+\Delta_t)\ &= \Delta_t^2\ \partial_{tt}\ u + 2\ u(t) - u(t-\Delta_t) \end{aligned} $$The argument for using a backward approximation is a bit hand wavy, but goes like this: a centered or forward approximation for $\partial_{t}\ u$ would involve the term $u(t+\Delta_t)$, and hence $u(t+\Delta_t)$ would appear at two places in our time update equation below, essentially making the form implicit (although it would be easy to solve for $u(t+\Delta_t)$).
We are interested in explicit time stepping and the correct behavior of the attenuation term, and so prefer the backward approximation for $\overleftarrow{\partial_{t}}\ u$. Our experience is that the use of the backward difference is more stable than forward or centered.
The first order accurate backward approximation to the first time derivative involves two wavefields: $u(t-\Delta_t)$, and $u(t)$. We can use this expression as is.
$$ \overleftarrow{\partial_{t}}\ u = \frac{\displaystyle u(t) - u(t-\Delta_t)}{\displaystyle \Delta_t} $$Next we plug in the right hand sides for $\partial_{tt}\ u$ and $\overleftarrow{\partial_{t}}\ u$ into the the time update expression for $u(t+\Delta_t)$ from step 2.
$$ \begin{aligned} u(t+\Delta_t) &= \Delta_t^2 \frac{m^2}{b} \left[ \overleftarrow{\partial_x}\left(b\ \overrightarrow{\partial_x}\ u \right) + \overleftarrow{\partial_y}\left(b\ \overrightarrow{\partial_y}\ u \right) + \overleftarrow{\partial_z}\left(b\ \overrightarrow{\partial_z}\ u \right) + q \right] \\[10pt] & \quad -\ \Delta_t^2 \frac{\omega_c}{Q} \left( \frac{\displaystyle u(t) - u(t-\Delta_t)} {\displaystyle \Delta_t} \right) + 2\ u(t) - u(t-\Delta_t) \end{aligned} $$The last equation is how we update the pressure wavefield at each time step, and depends on $u(t)$ and $u(t-\Delta_t)$.
The main work of the finite difference explicit time stepping is evaluating the nested spatial derivative operators on the RHS of this equation. The particular advantage of Devito symbolic optimization is that Devito is able to solve for the complicated expressions that result from substituting the discrete forms of high order numerical finite difference approximations for these nested spatial derivatives.
We have now completed the maths required to implement the modeling operator. The remainder of this notebook deals with setting up and using the required Devito objects.
Define the dimensions and coordinates for the model. The computational domain of the model is surrounded by an absorbing boundary region where we implement boundary conditions to eliminate outgoing waves. We define the sizes for the interior of the model nx and nz, the width of the absorbing boundary region npad, and the sizes for the entire model padded with absorbing boundaries become nxpad = nx + 2*npad and nzpad = nz + 2*npad.
In [4]:
# Define dimensions for the interior of the model
nx,nz = 751,751
dx,dz = 10.0,10.0 # Grid spacing in m
shape = (nx, nz) # Number of grid points
spacing = (dx, dz) # Domain size is now 5 km by 5 km
origin = (0., 0.) # Origin of coordinate system, specified in m.
extent = tuple([s*(n-1) for s, n in zip(spacing, shape)])
# Define dimensions for the model padded with absorbing boundaries
npad = 50 # number of points in absorbing boundary region (all sides)
nxpad,nzpad = nx+2*npad, nz+2*npad
shape_pad = np.array(shape) + 2 * npad
origin_pad = tuple([o - s*npad for o, s in zip(origin, spacing)])
extent_pad = tuple([s*(n-1) for s, n in zip(spacing, shape_pad)])
# Define the dimensions
# Note if you do not specify dimensions, you get in order x,y,z
x = SpaceDimension(name='x', spacing=Constant(name='h_x',
value=extent_pad[0]/(shape_pad[0]-1)))
z = SpaceDimension(name='z', spacing=Constant(name='h_z',
value=extent_pad[1]/(shape_pad[1]-1)))
# Initialize the Devito grid
grid = Grid(extent=extent_pad, shape=shape_pad, origin=origin_pad,
dimensions=(x, z), dtype=dtype)
print("shape; ", shape)
print("origin; ", origin)
print("spacing; ", spacing)
print("extent; ", extent)
print("")
print("shape_pad; ", shape_pad)
print("origin_pad; ", origin_pad)
print("extent_pad; ", extent_pad)
print("")
print("grid.shape; ", grid.shape)
print("grid.extent; ", grid.extent)
print("grid.spacing_map;", grid.spacing_map)
In [5]:
# Create the velocity and buoyancy fields.
# - We use a wholespace velocity of 1500 m/s
# - We use a wholespace density of 1 g/cm^3
# - These are scalar fields so we use Function to define them
# - We specify space_order to establish the appropriate size halo on the edges
space_order = 8
# Wholespace velocity
m = Function(name='m', grid=grid, space_order=space_order)
m.data[:] = 1.5
# Constant density
b = Function(name='b', grid=grid, space_order=space_order)
b.data[:,:] = 1.0 / 1.0
In this notebook we run 2 seconds of simulation using the sample rate related to the CFL condition as implemented in examples/seismic/skew_self_adjoint/utils.py.
Important note smaller Q values in highly viscous media may require smaller temporal sampling rates than a non-viscous medium to achieve dispersion free propagation. This is a cost of the visco- acoustic modification we use here.
We also use the convenience TimeRange as defined in examples/seismic/source.py.
In [6]:
t0 = dtype(0.) # Simulation time start
tn = dtype(2000.) # Simulation time end (1 second = 1000 msec)
dt = compute_critical_dt(m)
time_range = TimeAxis(start=t0, stop=tn, step=dt)
print("Time min, max, dt, num; %10.6f %10.6f %10.6f %d" % (t0, tn, dt, int(tn//dt) + 1))
print("time_range; ", time_range)
source:
examples/seismic/source.pyreceivers:
PointSource in examples/seismic/source.py
In [7]:
# Source in the center of the model at 10 Hz center frequency
fpeak = 0.010
src = RickerSource(name='src', grid=grid, f0=fpeak, npoint=1, time_range=time_range)
src.coordinates.data[0,0] = dx * (nx//2)
src.coordinates.data[0,1] = dz * (nz//2)
# line of receivers along the right edge of the model
rec = Receiver(name='rec', grid=grid, npoint=nz, time_range=time_range)
rec.coordinates.data[:,0] = dx * (nx//2)
rec.coordinates.data[:,1] = np.linspace(0.0, dz*(nz-1), nz)
print("src_coordinate X; %+12.4f" % (src.coordinates.data[0,0]))
print("src_coordinate Z; %+12.4f" % (src.coordinates.data[0,1]))
print("rec_coordinates X min/max; %+12.4f %+12.4f" % \
(np.min(rec.coordinates.data[:,0]), np.max(rec.coordinates.data[:,0])))
print("rec_coordinates Z min/max; %+12.4f %+12.4f" % \
(np.min(rec.coordinates.data[:,1]), np.max(rec.coordinates.data[:,1])))
# We can plot the time signature to see the wavelet
src.show()
In [8]:
# note: flip sense of second dimension to make the plot positive downwards
plt_extent = [origin_pad[0], origin_pad[0] + extent_pad[0],
origin_pad[1] + extent_pad[1], origin_pad[1]]
vmin, vmax = 1.4, 1.7
dmin, dmax = 0.9, 1.1
plt.figure(figsize=(12,8))
plt.subplot(1, 2, 1)
plt.imshow(np.transpose(m.data), cmap=cm.jet,
vmin=vmin, vmax=vmax, extent=plt_extent)
plt.colorbar(orientation='horizontal', label='Velocity (m/msec)')
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'white', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot(rec.coordinates.data[:, 0], rec.coordinates.data[:, 1], \
'black', linestyle='-', label="Receiver")
plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \
'red', linestyle='None', marker='*', markersize=15, label="Source")
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("Velocity w/ absorbing boundary")
plt.subplot(1, 2, 2)
plt.imshow(np.transpose(1 / b.data), cmap=cm.jet,
vmin=dmin, vmax=dmax, extent=plt_extent)
plt.colorbar(orientation='horizontal', label='Density (m^3/kg)')
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'white', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot(rec.coordinates.data[:, 0], rec.coordinates.data[:, 1], \
'black', linestyle='-', label="Receiver")
plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \
'red', linestyle='None', marker='*', markersize=15, label="Source")
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("Density w/ absorbing boundary")
plt.tight_layout()
None
We have two remaining constants and fields from our SSA wave equation that we need to define:
| Symbol | Description |
|---|---|
| $$\omega_c = 2 \pi f_c$$ | Center angular frequency |
| $$\frac{1}{Q(x,z)}$$ | Inverse Q model used in the modeling system |
The absorbing boundary condition strategy we use is designed to eliminate any corners or edges in the attenuation profile. We do this by making Q a function of distance from the nearest boundary.
We have implemented the function setup_w_over_q for 2D and 3D fields in the file utils.py, and will use it below. In Devito these fields are type Function, a concrete implementation of AbstractFunction.
Feel free to inspect the source at utils.py, which uses Devito's symbolic math to write a nonlinear equation describing the absorbing boundary for dispatch to automatic code generation.
Note that we will generate two Q models, one with strong attenuation (a Q value of 25) and one with moderate attenuation (a Q value of 100) -- in order to demonstrate the impact of attenuation in the plots near the end of this notebook.
In [9]:
#NBVAL_INGNORE_OUTPUT
# Initialize the attenuation profile for Q=25 and Q=100 models
w = 2.0 * np.pi * fpeak
print("w,fpeak; ", w, fpeak)
qmin = 0.1
wOverQ_025 = Function(name='wOverQ_025', grid=grid, space_order=space_order)
wOverQ_100 = Function(name='wOverQ_100', grid=grid, space_order=space_order)
setup_w_over_q(wOverQ_025, w, qmin, 25.0, npad)
setup_w_over_q(wOverQ_100, w, qmin, 100.0, npad)
# Plot the log of the generated Q profile
q025 = np.log10(w / wOverQ_025.data)
q100 = np.log10(w / wOverQ_100.data)
lmin, lmax = np.log10(qmin), np.log10(100)
plt.figure(figsize=(12,8))
plt.subplot(1, 2, 1)
plt.imshow(np.transpose(q025.data), cmap=cm.jet,
vmin=lmin, vmax=lmax, extent=plt_extent)
plt.colorbar(orientation='horizontal', label='log10(Q)')
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'white', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'white', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot(rec.coordinates.data[:, 0], rec.coordinates.data[:, 1], \
'black', linestyle='-', label="Receiver")
plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \
'red', linestyle='None', marker='*', markersize=15, label="Source")
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("log10 of $Q=25$ model")
plt.subplot(1, 2, 2)
plt.imshow(np.transpose(q100.data), cmap=cm.jet,
vmin=lmin, vmax=lmax, extent=plt_extent)
plt.colorbar(orientation='horizontal', label='log10(Q)')
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'white', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'white', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot(rec.coordinates.data[:, 0], rec.coordinates.data[:, 1], \
'black', linestyle='-', label="Receiver")
plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \
'red', linestyle='None', marker='*', markersize=15, label="Source")
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("log10 of $Q=100$ model")
plt.tight_layout()
None
TimeFunctionWe specify the time_order as 2, which allocates 3 time steps in the pressure wavefield. As described elsewhere, Devito will use "cyclic indexing" to index into this multi-dimensional array, mneaning that via the modulo operator, the time indices $[0, 1, 2, 3, 4, 5, ...]$ are mapped into the modulo indices $[0, 1, 2, 0, 1, 2, ...]$
This FAQ entry explains in more detail.
In [10]:
# Define the TimeFunction
u = TimeFunction(name="u", grid=grid, time_order=2, space_order=space_order)
# Get the symbols for dimensions for t, x, z
# We need these below in order to write the source injection and the
t,x,z = u.dimensions
If you examine the equation for the time update we derived above you will see that the source $q$ is scaled by the term $(\Delta_t^2 m^2\ /\ b)$. You will see that scaling term in the source injection below. For $\Delta_t^2$ we use the time dimension spacing symbol t.spacing**2.
Note that source injection and receiver extraction are accomplished via linear interpolation, as implemented in SparseTimeFunction in sparse.py.
In [11]:
# Source injection, with appropriate scaling
src_term = src.inject(field=u.forward, expr=src * t.spacing**2 * m**2 / b)
# Receiver extraction
rec_term = rec.interpolate(expr=u.forward)
We next transcribe the time update expression we derived above into a Devito Eq. Then we add that expression with the source injection and receiver extraction and build an Operator that will generate the c code for performing the modeling.
We copy the time update expression from above for clarity. Note we omit $q$ because we will be explicitly injecting the source using src_term defined immediately above. However, for the linearized Born forward modeling operation the $q$ term is an appropriately scaled field, as shown in the next notebook in this series.
In [12]:
#NBVAL_INGNORE_OUTPUT
# Generate the time update equation and operator for Q=25 model
eq_time_update = (t.spacing**2 * m**2 / b) * \
((b * u.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) + \
(b * u.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2)) + \
(2 - t.spacing * wOverQ_025) * u + \
(t.spacing * wOverQ_025 - 1) * u.backward
stencil = Eq(u.forward, eq_time_update)
# Update the dimension spacing_map to include the time dimension
# These symbols will be replaced with the relevant scalars by the Operator
spacing_map = grid.spacing_map
spacing_map.update({t.spacing : dt})
print("spacing_map; ", spacing_map)
# op = Operator([stencil] + src_term + rec_term)
op = Operator([stencil] + src_term + rec_term, subs=spacing_map)
The argument subs=spacing_map passed to the operator substitutes values for the temporal and spatial dimensions into the expressions before code generation. This reduces the number of floating point operations executed by the kernel by pre-evaluating certain coefficients, and possibly absorbing the spacing scalars from the denominators of the numerical finite difference approximations into the finite difference coefficients.
If you run the two cases of passing/not passing the subs=spacing_map argument by commenting/un-commenting the last two lines of the cell immediately above, you can inspect the difference in computed flop count for the operator. This is reported by setting Devito logging configuration['log-level'] = 'DEBUG' and is reported during Devito symbolic optimization with the output line Flops reduction after symbolic optimization. Note also if you inspect the generated code for the two cases, you will see extra calling parameters are required for the case without the substitution. We have compiled the flop count for 2D and 3D operators into the table below.
| Dimensionality | Passing subs | Flops reduction | Delta |
|---|---|---|---|
| 2D | False | 588 --> 81 | |
| 2D | True | 300 --> 68 | 13.7% |
| 3D | False | 875 --> 116 | |
| 3D | True | 442 --> 95 | 18.1% |
Note the gain in performance is around 14% for this example in 2D, and around 18% in 3D.
We use op.arguments() to print the arguments to the operator. As noted above depending on the use of subs=spacing_map you will see different arguments here. In the case of no subs=spacing_map argument to the operator, you will see arguments for the dimensional spacing constants as parameters to the operator, including h_x, h_z, and dt.
In [13]:
# NBVAL_IGNORE_OUTPUT
op.arguments()
Out[13]:
In [14]:
# NBVAL_IGNORE_OUTPUT
print(op)
By setting Devito logging configuration['log-level'] = 'DEBUG' we have enabled output of statistics related to the performance of the operator, which you will see below when the operator runs.
We will run the Operator once with the Q model as defined wOverQ_025, and then run a second time passing the wOverQ_100 Q model. For the second run with the different Q model, we take advantage of the placeholder design patten in the Devito Operator.
For more information on this see the FAQ entry.
In [15]:
#NBVAL_INGNORE_OUTPUT
# Run the operator for the Q=25 model
print("m min/max; %+12.6e %+12.6e" % (np.min(m.data), np.max(m.data)))
print("b min/max; %+12.6e %+12.6e" % (np.min(b.data), np.max(b.data)))
print("wOverQ_025 min/max; %+12.6e %+12.6e" % (np.min(wOverQ_025.data), np.max(wOverQ_025.data)))
print("wOverQ_100 min/max; %+12.6e %+12.6e" % (np.min(wOverQ_100.data), np.max(wOverQ_100.data)))
print(time_range)
u.data[:] = 0
op(time=time_range.num-1)
# summary = op(time=time_range.num-1, h_x=dx, h_z=dz, dt=dt)
# Save the Q=25 results and run the Q=100 case
import copy
pQ25 = copy.copy(u)
recQ25 = copy.copy(rec)
u.data[:] = 0
op(time=time_range.num-1, wOverQ_025=wOverQ_100)
print("Q= 25 receiver data min/max; %+12.6e %+12.6e" %\
(np.min(recQ25.data[:]), np.max(recQ25.data[:])))
print("Q=100 receiver data min/max; %+12.6e %+12.6e" %\
(np.min(rec.data[:]), np.max(rec.data[:])))
In [16]:
#NBVAL_INGNORE_OUTPUT
# Plot the two wavefields, normalized to Q=100 (the larger amplitude)
amax_Q25 = 1.0 * np.max(np.abs(pQ25.data[1,:,:]))
amax_Q100 = 1.0 * np.max(np.abs(u.data[1,:,:]))
print("amax Q= 25; %12.6f" % (amax_Q25))
print("amax Q=100; %12.6f" % (amax_Q100))
plt.figure(figsize=(12,8))
plt.subplot(1, 2, 1)
plt.imshow(np.transpose(pQ25.data[1,:,:] / amax_Q100), cmap="seismic",
vmin=-1, vmax=+1, extent=plt_extent)
plt.colorbar(orientation='horizontal', label='Amplitude')
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'black', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \
'red', linestyle='None', marker='*', markersize=15, label="Source")
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("Data for $Q=25$ model")
plt.subplot(1, 2, 2)
plt.imshow(np.transpose(u.data[1,:,:] / amax_Q100), cmap="seismic",
vmin=-1, vmax=+1, extent=plt_extent)
plt.colorbar(orientation='horizontal', label='Amplitude')
plt.plot([origin[0], origin[0], extent[0], extent[0], origin[0]],
[origin[1], extent[1], extent[1], origin[1], origin[1]],
'black', linewidth=4, linestyle=':', label="Absorbing Boundary")
plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \
'red', linestyle='None', marker='*', markersize=15, label="Source")
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("Data for $Q=100$ model")
plt.tight_layout()
None
In [17]:
#NBVAL_INGNORE_OUTPUT
# Plot the two receiver gathers, normalized to Q=100 (the larger amplitude)
amax_Q25 = 0.1 * np.max(np.abs(recQ25.data[:]))
amax_Q100 = 0.1 * np.max(np.abs(rec.data[:]))
print("amax Q= 25; %12.6f" % (amax_Q25))
print("amax Q=100; %12.6f" % (amax_Q100))
plt.figure(figsize=(12,8))
plt.subplot(1, 2, 1)
plt.imshow(recQ25.data[:,:] / amax_Q100, cmap="seismic",
vmin=-1, vmax=+1, extent=plt_extent, aspect="auto")
plt.colorbar(orientation='horizontal', label='Amplitude')
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("Receiver gather for $Q=25$ model")
plt.subplot(1, 2, 2)
plt.imshow(rec.data[:,:] / amax_Q100, cmap="seismic",
vmin=-1, vmax=+1, extent=plt_extent, aspect="auto")
plt.colorbar(orientation='horizontal', label='Amplitude')
plt.xlabel("X Coordinate (m)")
plt.ylabel("Z Coordinate (m)")
plt.title("Receiver gather for $Q=100$ model")
plt.tight_layout()
None
Note this takes a long time ... about 50 seconds, but obviates the need to solve for the time update expression as we did above.
If you would like to see the time update equation as generated by Devito symbolic optimization, uncomment the lines for the solve below.
In [18]:
#NBVAL_INGNORE_OUTPUT
# Define the partial_differential equation
# Note the backward shifted time derivative is obtained via u.dt(x0=t-0.5*t.spacing)
pde = (b / m**2) * (wOverQ_100 * u.dt(x0=t-0.5*t.spacing) + u.dt2) -\
(b * u.dx(x0=x+0.5*x.spacing)).dx(x0=x-0.5*x.spacing) -\
(b * u.dz(x0=z+0.5*z.spacing)).dz(x0=z-0.5*z.spacing)
# Uncomment the next 5 lines to see the equation as generated by Devito
# t1 = timer()
# stencil = Eq(u.forward, solve(pde, u.forward))
# t2 = timer()
# print("solve ran in %.4f seconds." % (t2-t1))
# stencil
A nonreflecting boundary condition for discrete acoustic and elastic wave equations (1985)
Charles Cerjan, Dan Kosloff, Ronnie Kosloff, and Moshe Resheq
Geophysics, Vol. 50, No. 4
https://library.seg.org/doi/pdfplus/10.1190/segam2016-13878451.1
Generation of Finite Difference Formulas on Arbitrarily Spaced Grids (1988)
Bengt Fornberg
Mathematics of Computation, Vol. 51, No. 184
http://dx.doi.org/10.1090/S0025-5718-1988-0935077-0
https://web.njit.edu/~jiang/math712/fornberg.pdf
Self-adjoint, energy-conserving second-order pseudoacoustic systems for VTI and TTI media for reverse time migration and full-waveform inversion (2016)
Kenneth Bube, John Washbourne, Raymond Ergas, and Tamas Nemeth
SEG Technical Program Expanded Abstracts
https://library.seg.org/doi/10.1190/segam2016-13878451.1