In this sequence of computer problems you will implement a simple example of seismic imaging by *reverse time migration*. This problem sequence will guide you step-by-step through implementing everything you need to do 1D seismic imaging.

These exercises are best solved using a single python script file, which I will call `cp1_sol.py`

(but you can name it whatever you wish).

I suggest you develop your solutions inside this file, and then run the script in an IPython console. This is accomplished by starting an IPython console, from an command prompt, with the command `ipython --pylab`

.

Then, from the IPython console, you can run your script with the command

`%run cp1_sol`

You will need a wave velocity model for these problems. These are provided in msri_tools.py and can be accessed by downloading the file to the directory you will work in. From an `ipython --pylab`

console you can access the model function by executing the command `from msri_tools import model`

. You can also put the line

`from msri_tools import model`

at the top of your python scripts.

In this exercise, the physical configuration in which we are working will have a significant amount of problem setup data (or meta data, if you will). It will be convenient to have a central place to store this. In Python we could use a `class`

as if we had a MATLAB or C style struct, however it is much easier to use Python's built-in dictionary type `dict`

. This configuration dictionary `config`

will be passed to all of your routines and should contain the physical and numerical/computational parameters of the problem, e.g., physical domain, number of points, time step, CFL factor, etc.

In [1]:

```
config = dict()
```

In this exercise, you will implement a solver for the 1D scalar acoustic wave equation with absorbing boundary conditions,

\[ \begin{align*} (\frac{1}{c(x)^2}\partial_{tt}-\partial_{xx})u(x,t) & = f(x,t), \\\\ (\frac{1}{c(0)}\partial_t-\partial_x)u(0,t) & = 0, \\\\ (\frac{1}{c(1)}\partial_t+\partial_x)u(1,t) & = 0, \\\\ u(x,t) & = 0 \quad\text{for}\quad t \le 0, \end{align*} \]

where the last two equations are the absorbing boundary conditions, \(x \in [0,1]\) and \(t \in [0,T]\).

In our notation we write that solving this PDE is equivalent to applying a nonlinear operator \(\mathcal{F}\) to a model parameter \(m\), and so for

\[ m(x) = \frac{1}{c(x)^2},\\\\ \mathcal{F}[m] = u. \]

We define our source functions as \(f(x,t) = w(t)\delta(x-x_0)\), where \(w\) is the time profile and the \(\delta\) indicates that we will use point sources. For seismic, the time profile is not known and is estimated as part of the inverse problem. However, it is common to model source signals with the negative second derivative of a Gaussian, also known as the Ricker Wavelet,

\[ w(t) = (1-2\pi^2\nu_0^2t^2)e^{-\pi^2\nu_0^2t^2}, \]

where \(\nu_0\) is known as the characteristic or peak frequency (in Hz), because it is the \(|\hat w|\) attains its maximum at that point. It is also important that this function is causal (\(w(t) = 0\) for \(\le 0\)), so we introduce a time shift \(t_0\),

\[ w(t) = (1-2\pi^2\nu_0^2(t-t_0)^2)e^{-\pi^2\nu_0^2(t-t_0)^2}. \]

**Problem 1:** Write a Python function `ricker(t, config)`

which implements the Ricker Wavelet, taking a time `t`

in seconds and a peak frequency `nu0`

in Hz and returns the value of the wavelet at time `t`

. You can guarantee causality by setting \(t_0= 6\sigma\) for \(\sigma = \tfrac{1}{\pi\nu_0\sqrt{2}}\), the standard deviation of the underlying Gaussian. You may also want to implement an optional threshold to prevent excessively small numbers.

Plot your function for \(t = 0, \dots, T=0.5\) at \(\nu_0 = 10\textrm{Hz}\) and label the plot.

In [2]:

```
# Import my solutions
from cp1_sol import ricker
config['nu0'] = 10 #Hz
ts = linspace(0, 0.5, 1000)
ws = ricker(ts, config)
plot(ts, ws, color='green',
label=r'$\nu_0 =\,{0}$Hz'.format(config['nu0']),
linewidth=2)
xlabel(r'$t$', fontsize=18)
ylabel(r'$w(t)$', fontsize=18)
title('Ricker Wavelet', fontsize=22)
legend();
```

**Problem 2:** Write a Python function `point_source(value, position, config)`

which takes a value `value`

, a source location `position`

, and uses the range of the spatial domain `x_limits`

, and the number of points `nx`

from the configuration to implement a numerical approximation to the \(\delta\). This function should return an array with `value`

at the correct index, correctly evaluating \(w(t)\delta(x-x_s)\) for `value = ricker(t)`

. Be careful with your implementation of the numerical delta, as \(\int \delta(x-x_s)\textrm{d}x = 1\).

In [3]:

```
# Domain parameters
config['x_limits'] = [0.0, 1.0]
config['nx'] = 201
config['dx'] = (config['x_limits'][1] - config['x_limits'][0]) / (config['nx']-1)
# Source parameter
config['x_s'] = 0.1
```

The scalar acoustic wave equation can be written as

\[ \begin{align*} (M\partial_{tt} + A\partial_t + K)u(x,t) = f(x,t), \end{align*} \]

where \(K = K_x + K_{xx}\) is called the stiffness matrix and contains the spatial derivatives, \(A\) is the attenuation matrix and relates to the first time derivatives in the boundary conditions, and \(M\) is the mass matrix and relates to the second time derivatives in the bulk.

**Problem 3:** Write a Python function which `construct_matrices(C, config)`

which constructs the matrices \(M\), \(A\), \(K\) for a given model velocity `C`

. Use a second order accurate finite difference for the second spatial derivative and and use an 'upwinded' forward or backward first order difference scheme for the first spatial derivatives.

It may be helpful to write down precisely what the differential equation looks like at the each interesting point (the two boundaries and some point in the middle of the domain) for your discretized wavefield \(u(j\Delta x, n \Delta t)\). Hint: this is at one time step, so \(n\) is fixed. Which \(j\) are relevant at each of the spatial points?

In [4]:

```
from msri_tools import model
from cp1_sol import construct_matrices
# Load the model
C, C0 = model(config)
# Build an example set of matrices
M, A, K = construct_matrices(C, config)
```

We can discretize the time derivatives using the usual second-order accurate finite difference approximation,

\[ \partial_{tt}u(x,t) \approx \frac{u(x,t-\Delta t) - 2u(x,t) + u(x,t+\Delta t)}{\Delta t^2}, \]

which will result in the explicit 'leap-frog' scheme. For explicit methods, the stability of the scheme is restricted by the Courant-Friedrichs-Lewy (CFL) condition,

\[ \Delta t \le \alpha\frac{\Delta x}{c_\text{max}}. \]

**Problem 4:** Write a Python function `leap_frog(C, sources, config)`

which takes a velocity `C`

, the spatial parameters, a time step `dt`

, a number of time steps `nt`

, and a list of source wavefields (one for each time step) and returns the time series of wavefields \(u\). Use \(\alpha = \dfrac{1}{6}\) and \(x_s = 0.1\).

In [5]:

```
from cp1_sol import point_source, leap_frog
# Set CFL safety constant
config['alpha'] = 1.0/6.0
# Define time step parameters
config['T'] = 3 #seconds
config['dt'] = config['alpha'] * config['dx'] / C.max()
config['nt'] = int(config['T']/config['dt'])
dt = config['dt']
sources = [point_source(ricker(i*dt, config), config['x_s'], config) \
for i in xrange(config['nt'])]
us = leap_frog(C, sources, config)
```

**Problem 5:** Write a Python function `plot_space_time(us, config)`

, using the matplotlib command `imshow`

, to plot and label the space-time diagram for a wavefield \(u(x,t)\). Hint: the matplotlib `xticks`

and `yticks`

functions will be useful. Try to use a gray-scale color map and consider an optional argument to set the title.

In [6]:

```
from cp1_sol import plot_space_time
plot_space_time(us, config, title=r'u(x,t)')
```

**Bonus Problem 1:** Derive how you might use your `leap_frog`

function and periodic boundary conditions to design a 4th order accurate, in both space and time, scheme for solving the wave equation.

**Bonus Problem 2:** The implementation of time stepping used in this exercise is not the most efficient approach for implementing time stepping, particularly in higher dimensions. Why? What might be a faster way to implement the time stepping?

**Bonus Problem 3:** The wave equation can be solved using an ODE integrator. Change the formulation of the wave equation so that this is poddible. Write a function that uses the built-in SciPy ODE integrator to do your time stepping.