Implementation of a Devito self adjoint variable density visco- acoustic isotropic modeling operator
-- Nonlinear Ops --

This operator is contributed by Chevron Energy Technology Company (2020)¶

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

Introduction¶

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 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:

1. Implementation of a Devito self adjoint variable density visco- acoustic isotropic modeling operator -- Nonlinear Ops¶

• Implement the nonlinear modeling operations.
• sa_01_iso_implementation1.ipynb

2. Implementation of a Devito self adjoint variable density visco- acoustic isotropic modeling operator -- Linearized Ops¶

• Implement the linearized (Jacobian) forward and adjoint modeling operations.
• sa_02_iso_implementation2.ipynb

3. Implementation of a Devito self adjoint variable density visco- acoustic isotropic modeling operator -- Correctness Testing¶

• Tests the correctness of the implemented operators.
• sa_03_iso_correctness.ipynb

There are similar series of notebooks implementing and testing operators for VTI and TTI anisotropy (README.md).

Below we introduce the 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.

Outline¶

1. Define symbols
2. Introduce the SA wave equation
3. Show generation of skew symmetric derivatives and prove correctness with unit test
4. Derive the time update equation used to implement the nonlinear forward modeling operator
5. Create the Devito grid and model fields
6. Define a function to implement the attenuation profile ($\omega\ /\ Q$)
7. Create the Devito operator
8. Run the Devito operator
9. Plot the resulting wavefields
10. References

Table of symbols¶

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$

A word about notation¶

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$.

Self adjoint variable density visco- acoustic wave equation¶

Our 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:

• the nonlinear forward
• the linearized forward (Jacobian forward)
• the linearized adjoint (Jacobian adjoint)

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 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 self adjoint (SA) wave equation.

Imports¶

We have grouped all imports used in this notebook here for consistency.

Unit test demonstrating skew symmetry for shifted derivatives¶

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:

$$ \begin{aligned} f_2 = \overrightarrow{\partial_x}\ f_1 \\[5pt] g_2 = \overleftarrow{\partial_x}\ g_1 \end{aligned} $$

And then verify the dot products are equivalent. Recall that skew symmetry introduces a minus sign in this equality:

$$ 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 $$

Show the finite difference operators and generated code¶

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.

Doing some algebra to solve for the time update¶

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.

1. Numerical approximation for $\partial_{tt}\ u$, solved for for $u(t+\Delta_t)$¶

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} $$

2. Numerical approximation for $\overleftarrow{\partial_{t}}\ u$¶

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} $$

3. Solve the wave equation for $\partial_{tt}$¶

$$ \begin{aligned} \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\\[10pt] \partial_{tt}\ u &= \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] - \frac{\omega_c}{Q} \overleftarrow{\partial_{t}}\ u \end{aligned} $$

4. Plug in $\overleftarrow{\partial_t} u$ and $\partial_{tt} u$ into the time update equation¶

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} $$

5. Simplify ...¶

$$ \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 \left(2 -\ \Delta_t\ \frac{\omega_c}{Q} \right) u(t) + \left(\Delta_t\ \frac{\omega_c}{Q} - 1 \right) u(t-\Delta_t) \end{aligned} $$

6. et voila ...¶

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.

Instantiate the Devito grid for a two dimensional problem¶

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.

Define velocity and buoyancy model parameters¶

We have the following constants and fields from our self adjoint wave equation that we define as time invariant using Functions:

Symbol	Description
$$m(x,z)$$	Acoustic velocity
$$b(x,z)=\frac{1}{\rho(x,z)}$$	Buoyancy (reciprocal density)

Define the simulation time range¶

In this notebook we run 2 seconds of simulation using the sample rate related to the CFL condition as implemented in examples/seismic/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.

Define the acquisition geometry: locations of sources and receivers¶

source:

• X coordinate: center of the model: dx*(nx//2)
• Z coordinate: center of the model: dz*(nz//2)
• We use a 10 Hz center frequency RickerSource wavelet as defined in examples/seismic/source.py

receivers:

• X coordinate: center of the model: dx*(nx//2)
• Z coordinate: vertical line from top to bottom of model
• We use a vertical line of Receivers as defined with a PointSource in examples/seismic/source.py

Plot velocity and density models¶

Next we plot the velocity and density models for illustration.

• The demarcation between interior and absorbing boundary is shown with a dotted white line
• The source is shown as a large red asterisk
• The extent of the receiver array is shown with a thick black line

Create and plot the $\frac{\omega_c}{Q}$ model used for dissipation only attenuation¶

We have two remaining constants and fields from our SA 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.

Define the pressure wavefield as a TimeFunction¶

We 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.

Define the source injection and receiver extraction¶

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.

Finally, the Devito operator¶

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.

$$ \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 \left(2 -\ \Delta_t\ \frac{\omega_c}{Q} \right) u(t) + \left(\Delta_t\ \frac{\omega_c}{Q} - 1 \right) u(t-\Delta_t) \end{aligned} $$

Impact of hardwiring the grid spacing on operation count¶

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.

Print the arguments to the Devito operator¶

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.

Print the generated c code for review¶

We use print(op) to output the generated c code for review.

Run the operator for the Q=25 and Q=100 models¶

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.

Plot the computed Q=25 and Q=100 wavefields¶

Plot the computed Q=25 and Q=100 receiver gathers¶

Show the output from Devito solving for the stencil¶

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.

Discussion¶

This concludes the implementation of the nonlinear forward operator. This series continues in the next notebook that describes the implementation of the Jacobian linearized forward and adjoint operators.

sa_02_iso_implementation2.ipynb

References¶

• 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