Computational Seismology
Finite Differences - Advection Diffusion Reaction Equation

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

##### Authors:¶

This exercise covers the following aspects:

• Solving the advection-diffusion-reaction equation with different schemes of finite difference method
• Investigation of different numerical solutions and finding the stable one

## Basic Equations¶

Please refer to the Exercise 4.20 on Page 105 from the book.

The advection-diffusion-reaction equation is given by

$$\partial_t C(x,t) = k \partial_x^2 C(x,t) + v \partial_x C(x,t) - R(x)C(x,t) + P(x)$$

Where,

$C(x,t)$ is the concentration,

$k$ is the diffusivity,

$R(x)$ is the reactivity,

$P(x)$ is the source, and

$v$ is the advection velocity.

We replace the derivatives as follows

$$\partial_x C(x,t) \ \approx \ \frac{C(x + \mathrm{d}x, t) - C(x -\mathrm{d}x,t)} {2 \mathrm{d}x}$$$$\partial_x^2 C(x,t) \ \approx \ \frac{C(x + \mathrm{d}x, t) - 2C(x,t) + C(x - \mathrm{d}x, t)}{\mathrm{d}x^2}$$$$\partial_t C(x,t) \ \approx \ \frac{C(x, t + \mathrm{d}t) - C(x,t)}{ \mathrm{d}t}$$

The R.H.S. of the advection-diffusion-reaction becomes

$$R.H.S. = dc = k \bigg[ \frac{C(x + \mathrm{d}x) -2 C(x) + C(x - \mathrm{d}x)}{\mathrm{d} x^2} \bigg] + v \bigg[ \frac{C(x + \mathrm{d}x) - C(x -\mathrm{d}x)}{2dx} \bigg] - RC + P$$

Now, solving for $C(x, t + \mathrm{d}t)$ the extrapolation scheme is

$$C(x, t + \mathrm{d}t) = dc \mathrm{d}t + C(t) = k \mathrm{d}t \bigg[ \frac{C(x + \mathrm{d}x) -2 C(x) + C(x - \mathrm{d}x)}{\mathrm{d} x^2} \bigg] + v \mathrm{d}t \bigg[ \frac{C(x + \mathrm{d}x) - C(x -\mathrm{d}x)}{2dx} \bigg] - RC \mathrm{d}t + P \mathrm{d}t + C(t)$$

Note that we used centered finite-difference (in velocity term) scheme.

### Exercise¶

First understand the following codes, run the simulation and interpret the solution (plot).

Then, implement the forward and the backward (in velocity term) finite-difference schemes. Investigate different numerical solutions. Which scheme is stable?

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 [1]:
# 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
from mpl_toolkits.mplot3d import Axes3D               # for 3D plot

In [2]:
# Initialization of parameters

k = 3*np.exp(4)                                       # eddie mixing (diffusivity) in km^2/year
v = 1500.                                             # velocity km/year
p0 = 0.01                                             # production rate in mol/year
r = 2*np.exp(-3)                                      # removal rate (reactivity) in mol/year

# Choose among "centered", "forward" or "backward" finite-difference scheme.
fd_type = "forward"
#fd_type = "centered"
#fd_type = "backward"

# Space discretization
nt = 300                                              # number of time steps
nx = 100                                              # space dimension

phi = 2*np.pi/(nx-1)*np.arange(0,nx)
xc = np.cos(phi)                                      # for visualization of eddie like circular plot
yc = np.sin(phi)                                      # for visualization of eddie like circular plot

dt = 0.01                                             # years
dx = 1000                                             # km
isnap = 5                                             # snapshot frequency for visualization

# Initialization of variables
c = np.zeros(nx)                                      # concentration
dc =np.zeros(nx)                                      # R.H.S. of equation
p = np.zeros(nx)                                      # source
p[int((nx-1)/2)] = p0                                 # injecting source (production rate)

In [3]:
# Extrapolation scheme and the plots

# Initialization of plot

fig = plt.figure()
ax = fig.gca(projection='3d')
lines = ax.plot(xc,yc,c, 'r-', label = 'concentration')
plt.legend(loc='best')
plt.ion()
plt.show()
title = "Concentration in"

# Begin extrapolation and update the plot
# Circular shifting is done by using numpy function roll in the codes below

for i in range(1,nt):
if fd_type == "centered":

dc = k * (np.roll(c,1) - 2 * c + np.roll(c,-1)) / (dx ** 2) + v * (np.roll(c,1) -
np.roll(c,-1)) / (2 * dx) + p - r * c

if fd_type == "backward":

dc = k * (np.roll(c,1) - 2 * c + np.roll(c,-1)) / (dx ** 2) + v * (np.roll(c,0) -
np.roll(c,-1)) / dx + p - r * c

if fd_type == "forward":

dc = k * (np.roll(c,1) - 2 * c + np.roll(c,-1)) / (dx ** 2) + v * (np.roll(c,1) -
np.roll(c,0)) / dx + p - r * c

# Time extrapolation scheme (Euler)
cnew = c + dt * dc

# The new presence is the current future!
c = cnew

c[0]=c[nx-1]                                      # Spatial boundary condition

# Updating the plots
if i % isnap == 0:
for l in lines:
l.remove()
del l
lines = ax.plot(xc,yc,c, 'r-')
plt.title(title + " time: %.2g years" % (i * dt), loc = 'left')
plt.gcf().canvas.draw()

plt.ioff()
plt.show()

In [ ]: