Computational Seismology

Finite Differences - Grid-Staggering Elastic 1D

This notebook is part of the supplementary material to Computational Seismology: A Practical Introduction, Oxford University Press, 2016.

- Ashim Rijal (@ashimrijal)
- Heiner Igel (@heinerigel)

This exercise covers the following aspects:

- Solving velocity-stress formulation of 1D wave equation with finite difference method
- Understanding the grid-staggering in connection with finite difference solution to the elastic wave equation

** Please refer to the sections 4.6.2 and 4.6.3 from the book.**

The 1D wave equation (velocity-stress formulation) as a coupled system of two first-order partial differential equations

$$ \rho \partial_t v = \partial_x \sigma + f $$$$ \partial_t \sigma = \mu \partial_x v $$where,

$ \sigma $ is the stress,

$ \rho $ is the density,

$ v $ is the velocity,

$ \mu $ is the shear modulus, and

$ f $ is the source.

Grid- staggering is the concept in connection with finite-difference solutions to the elastic wave equation. The discrete velocity and stress are defined on a regular spaced grid in space and time. Then, partial derivatives are replaced with centered finite-difference approximations to first derivative. However, these are not defined at the grid points of a function but in-between the grid points. In grid staggering the following computational scheme is used

$$ \frac{v_i^{j+ \tfrac{1}{2}} - v_i^{j- \tfrac{1}{2}} }{dt} \ = \ \frac{1}{\rho_i}\frac{\sigma_{i + \tfrac{1}{2}}^j - \sigma_{i - \tfrac{1}{2}}^j }{dx} + \frac{f_i^j}{\rho_i} \ $$$$ \frac{\sigma_{i+\tfrac{1}{2}}^{j+1} - \sigma_{i+\tfrac{1}{2}}^j }{dt} \ = \ \mu_{i+\tfrac{1}{2}} \frac{v_{i + 1}^{j +\tfrac{1}{2}} - v_i^{j + \tfrac{1}{2}} }{dx} $$The extrapolation scheme becomes

$$ v_i^{j+ \tfrac{1}{2}} \ = \ \frac{dt}{\rho_i} \frac{\sigma_{i + \tfrac{1}{2}}^j - \sigma_{i - \tfrac{1}{2}}^j }{dx} \ + \ v_i^{j- \tfrac{1}{2}} + \frac{dt}{\rho_i} \ f_i^j \ $$$$ \sigma_{i+\tfrac{1}{2}}^{j+1} \ = \ dt \ \mu_{i+\tfrac{1}{2}} \frac{v_{i + 1}^{j +\tfrac{1}{2}} - v_i^{j + \tfrac{1}{2}} }{dx} \ + \ \sigma_{i+\tfrac{1}{2}}^j \ \ $$** Note that in the codes below we do not deal with the index fractions.**

First understand the codes below and run the simulation.

Then, improve the result using (4-point operator) for 1st derivative.

**Message: Once you become familiar with all the codes below you can go to the Cell tab on the toolbar and click Run All.**

In [ ]:

```
# Configuration step (Please run it before the simulation code!)
import numpy as np
import matplotlib
# Show Plot in The Notebook
matplotlib.use("nbagg")
import matplotlib.pyplot as plt
matplotlib.rcParams['figure.facecolor'] = 'w' # remove grey background
```

In [ ]:

```
# Initialization of parameters
# Simple finite difference solver
# Elastic wave equation
# 1-D regular staggered grid
# Basic parameters
nt = 1300 # number of time steps
nx = 1000 # number of grid points in x
c0 = 4500 # velocity (m/sec) (shear wave)
eps = 0.8 # stability limit
isnap = 2 # snapshot frequency
isx = round(nx/2) # source location
f0 = 1/15 # frequency (Hz)
xmax = 1000000. # maximum range (m)
rho0 = 2500. # density (kg/m**3)
mu0 = rho0*c0**2. # shear modulus (Pa)
nop = 4 # number of operator either 2 or 4
dx = xmax / (nx-1) # calculate space increment (m)
x = (np.arange(nx)*dx) # initialize space coordinates
dt = eps * dx/c0 # calculate time step from stability criterion(s)
# Source time function
t = (np.arange(0,nt) * dt) # initialize time axis
T0 = 1. / f0 # period
a = 4. / T0 # half-width (so called sigma)
t0 = T0 / dt
tmp = np.zeros(nt)
for it in range(nt):
t = (it - t0) * dt
tmp[it] = -2 * a * t * np.exp(-(a * t) ** 2) # derivative of Gaussian (so called sigma)
src = np.zeros(nt) # source
src[0:len(tmp)] = tmp
lam = c0 * T0 # wavelength
```

In [ ]:

```
# Extrapolation scheme and the plots
# Initialization of plot
# Initialization of fields
v = np.zeros(nx) # velocity
vnew = v
dv = v
s = np.zeros(nx) # stress
snew = s
ds = s
mu = np.zeros(nx) # shear modulus
rho = mu
rho = rho + rho0
mu = mu + mu0
# Print setup parameters
print("rho =", rho0, ", f_dom =", f0, ", stability limit =", eps, ", n_lambda", (lam/dx))
# Initialize the plot
title = "FD Elastic 1D staggered grid"
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
line1 = ax1.plot(x, v, color = "red", lw = 1.5)
line2 = ax2.plot(x, s, color = "blue", lw = 1.5)
ax1.set_ylabel('velocity (m/s)')
ax2.set_xlabel('x (m)')
ax2.set_ylabel('stress (Pa)')
plt.ion()
plt.show()
# Begin extrapolation and update the plot
for it in range (nt):
# Stress derivative
for i in range (2, nx-2):
ds[i] = (0.0416666 * s[i-1] - 1.125 * s[i] + 1.125 * s[i+1] - 0.0416666 * s[i+2]) / (dx)
# Velocity extrapolation
v = v + dt * ds / rho
# Add source term at isx
v[isx] = v[isx] + dt * src[it] / (dt * rho[isx])
# Velocity derivative
for i in range (2, nx-2):
dv[i] = (0.0416666 * v[i-2] - 1.125 * v[i-1] + 1.125 * v[i] - 0.0416666 * v[i+1]) / (dx)
# Stress extrapolation
s = s + dt * mu * dv
# Updating the plots
if not it % isnap:
for l in line1:
l.remove()
del l
for l in line2:
l.remove()
del l
line1 = ax1.plot(x, v, color = "red", lw = 1.5)
line2 = ax2.plot(x, s, color = "blue", lw = 1.5)
ax1.set_title(title + ", time step: %i" % (it))
plt.gcf().canvas.draw()
plt.ioff()
plt.show()
```

In [ ]:

```
```