{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Finite Volume Method - 1D Elastic Wave Equation
\n", "
\n", "
\n", "
"]}, {"cell_type": "markdown", "metadata": {}, "source": ["

\n", " \n", "\n", "\n", "

\n", "\n", "\n", "---\n", "\n", "This notebook is part of the supplementary material \n", "to [Computational Seismology: A Practical Introduction](https://global.oup.com/academic/product/computational-seismology-9780198717416?cc=de&lang=en&#), \n", "Oxford University Press, 2016.\n", "\n", "##### Authors:\n", "* David Vargas ([@dvargas](https://github.com/davofis))\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Basic Equations\n", "\n", "The source free elastic wave equation can be written in terms of a coupled first-order system \n", "\n", "\\begin{align}\n", "\\partial_t \\sigma - \\mu \\partial_x v & = 0 \\\\ \n", "\\partial_t v - \\frac{1}{\\rho} \\partial_x \\sigma & = 0\n", "\\end{align}\n", "\n", "with $\\rho$ the density and $\\mu$ the shear modulus. This equation in matrix-vector notation follows\n", "\n", "\\begin{equation}\n", "\\partial_t \\mathbf{Q} + \\mathbf{A} \\partial_x \\mathbf{Q} = 0\n", "\\end{equation}\n", "\n", "where $\\mathbf{Q} = (\\sigma, v)$ is the vector of unknowns and the matrix $\\mathbf{A}$ contains the parameters $\\rho$ and $\\mu$. The above matrix equation is analogous to the advection equation $\\partial_t q + a \\partial_x q = 0$. Although it is a coupled system, diagonalization of $\\mathbf{A} = \\mathbf{R}^{-1}\\mathbf{\\Lambda}\\mathbf{R}$ allows us to implement all elements developed for the solution of the advection equation in terms of fluxes. It turns out that the decoupled version is\n", "\n", "\\begin{equation}\n", "\\partial_t \\mathbf{W} + \\mathbf{\\Lambda} \\partial_x \\mathbf{W} = 0\n", "\\end{equation}\n", "\n", "where the eigenvector matrix $\\mathbf{R}$ and the diagonal matrix of eigenvalues $\\mathbf{\\Lambda}$ are given\n", "\n", "\\begin{equation}\n", "\\mathbf{W} = \\mathbf{R}^{-1}\\mathbf{Q}\n", "\\qquad\\text{,}\\qquad\n", "\\mathbf{\\Lambda}=\n", " \\begin{pmatrix}\n", " -c & 0 \\\\\n", " 0 & c \n", " \\end{pmatrix}\n", "\\qquad\\text{,}\\qquad\n", "\\mathbf{R} = \n", " \\begin{pmatrix}\n", " Z & -Z \\\\\n", " 1 & 1 \n", " \\end{pmatrix}\n", "\\qquad\\text{and}\\qquad\n", "\\mathbf{R}^{-1} = \\frac{1}{2Z}\n", " \\begin{pmatrix}\n", " 1 & Z \\\\\n", " -1 & Z \n", " \\end{pmatrix}\n", "\\end{equation}\n", "\n", "Here $Z = \\rho c$ with $c = \\sqrt{\\mu/\\rho}$ represents the seismic impedance. \n", "\n", "This notebook implements both upwind and Lax-Wendroff schemes for solving the free source version of the elastic wave equation in a homogeneous media. To keep the problem simple we use as spatial initial condition a Gauss function with half-width $\\sigma$\n", "\n", "\\begin{equation}\n", "Q(x,t=0) = e^{-1/\\sigma^2 (x - x_{o})^2}\n", "\\end{equation}"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": }, "outputs": [], "source": ["# Import all necessary libraries, this is a configuration step for the exercise.\n", "# Please run it before the simulation code!\n", "import numpy as np\n", "import matplotlib\n", "# Show Plot in The Notebook\n", "matplotlib.use(\"nbagg\")\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### 1. Initialization of setup"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": []}, "outputs": [], "source": ["# Initialization of setup\n", "# --------------------------------------------------------------------------\n", "nx = 800 # number of grid points \n", "c = 2500 # acoustic velocity in m/s\n", "ro = 2500 # density in kg/m^3\n", "Z = ro*c # impedance\n", "mu = ro*c**2 # shear modulus\n", "\n", "xmax = 10000 # Length in m \n", "eps = 0.5 # CFL\n", "tmax = 2.0 # simulation time in s\n", "isnap = 10 # plotting rate\n", "sig = 200 # argument in the inital condition\n", "x0 = 5000 # position of the initial condition\n", "\n", "imethod = 'upwind' # 'Lax-Wendroff', 'upwind'\n", "\n", "# Initialize Space\n", "x, dx = np.linspace(0,xmax,nx,retstep=True)\n", "\n", "# use wave based CFL criterion\n", "dt = eps*dx/c # calculate time step from stability criterion\n", "\n", "# Simulation time\n", "nt = int(np.floor(tmax/dt))\n", "\n", "# Initialize wave fields\n", "Q = np.zeros((2,nx))\n", "Qnew = np.zeros((2,nx))\n", "Qa = np.zeros((2,nx))"]}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 2. Initial condition\n", "\n", "Seismic disturbances are introduced through specification of a particular spatial initial condition, in this case we use a Gaussian function. \n", "\n", "#### Exercise 1 \n", "Implement the spatial initial condition given by:\n", "\n", "\\begin{equation}\n", "Q(x,t=0) = e^{-1/\\sigma^2 (x - x_{o})^2}\n", "\\end{equation}\n", "\n", "Then, visualize the initial condition in a given plot."]}, {"cell_type": "code", "execution_count": null, "metadata": {"lines_to_next_cell": 2, "tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# INITIALIZE THE SOURCE TIME FUCTION HERE!\n", "#################################################################\n", "\n", "\n", "#################################################################\n", "# PLOT THE SOURCE TIME FUNCTION HERE!\n", "#################################################################"]}, {"cell_type": "code", "execution_count": null, "metadata": {"solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 3. Solution for the homogeneous problem\n", "\n", "#### Upwind finite volume scheme\n", "We decompose the solution into right propagating $\\mathbf{\\Lambda}^{+}$ and left propagating eigenvalues $\\mathbf{\\Lambda}^{-}$ where\n", "\n", "\\begin{equation}\n", "\\mathbf{\\Lambda}^{+}=\n", " \\begin{pmatrix}\n", " -c & 0 \\\\\n", " 0 & 0 \n", " \\end{pmatrix}\n", "\\qquad\\text{,}\\qquad\n", "\\mathbf{\\Lambda}^{-}=\n", " \\begin{pmatrix}\n", " 0 & 0 \\\\\n", " 0 & c \n", " \\end{pmatrix}\n", "\\qquad\\text{and}\\qquad\n", "\\mathbf{A}^{\\pm} = \\mathbf{R}^{-1}\\mathbf{\\Lambda}^{\\pm}\\mathbf{R}\n", "\\end{equation}\n", "\n", "This strategy allows us to formulate an upwind finite volume scheme for any hyperbolic system as \n", "\n", "\\begin{equation}\n", "\\mathbf{Q}_{i}^{n+1} = \\mathbf{Q}_{i}^{n} - \\frac{dt}{dx}(\\mathbf{A}^{+}\\Delta\\mathbf{Q}_{l} - \\mathbf{A}^{-}\\Delta\\mathbf{Q}_{r})\n", "\\end{equation}\n", "\n", "with corresponding flux term given by\n", "\n", "\\begin{equation}\n", "\\mathbf{F}_{l} = \\mathbf{A}^{+}\\Delta\\mathbf{Q}_{l}\n", "\\qquad\\text{,}\\qquad\n", "\\mathbf{F}_{r} = \\mathbf{A}^{-}\\Delta\\mathbf{Q}_{r}\n", "\\end{equation}\n", "\n", "#### Lax-Wendroff finite volume scheme\n", "\n", "The upwind solution presents a strong diffusive behaviour. In this sense, the Lax-Wendroff perform better, with the advantage that it is not needed to decompose the eigenvalues into right and left propagations. Here the matrix $\\mathbf{A}$ can be used in its original form. The Lax-Wendroff follows\n", "\n", "\\begin{equation}\n", "\\mathbf{Q}_{i}^{n+1} = \\mathbf{Q}_{i}^{n} - \\frac{dt}{2dx}\\mathbf{A}(\\mathbf{Q}_{i+1}^{n} - \\mathbf{Q}_{i-1}^{n}) + \\frac{1}{2}\\Big(\\frac{dt}{dx}\\Big)^2\\mathbf{A}^2(\\mathbf{Q}_{i+1}^{n} - 2\\mathbf{Q}_{i}^{n} + \\mathbf{Q}_{i-1}^{n})\n", "\\end{equation} \n", "\n", "#### Exercise 2\n", "Initialize all relevant matrices, i.e $R$, $R^{-1}$, $\\mathbf{\\Lambda}^{+}$, $\\mathbf{\\Lambda}^{-}$, $\\mathbf{A}^{+}$, $\\mathbf{A}^{-}$, $\\mathbf{A}$. "]}, {"cell_type": "code", "execution_count": null, "metadata": {"tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# INITIALIZE ALL MATRICES HERE!\n", "#################################################################\n", "# R = \n", "# Rinv = \n", "# Lp = \n", "# Lm = \n", "# Ap = \n", "# Am = \n", "# A = "]}, {"cell_type": "code", "execution_count": null, "metadata": {"solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### 4. Finite Volumes solution"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": []}, "outputs": [], "source": ["# Initialize animated plot\n", "# ---------------------------------------------------------------\n", "fig = plt.figure(figsize=(10,6))\n", "ax1 = fig.add_subplot(2,1,1)\n", "ax2 = fig.add_subplot(2,1,2)\n", "line1 = ax1.plot(x, Q[0,:], 'k', x, Qa[0,:], 'r--')\n", "line2 = ax2.plot(x, Q[1,:], 'k', x, Qa[1,:], 'r--')\n", "ax1.set_ylabel('Stress')\n", "ax2.set_ylabel('Velocity')\n", "ax2.set_xlabel(' x ')\n", "plt.suptitle('Homogeneous F. volume - %s method'%imethod, size=16)\n", "\n", "plt.ion() # set interective mode\n", "plt.show()\n", "\n", "# ---------------------------------------------------------------\n", "# Time extrapolation\n", "# --------------------------------------------------------------- \n", "for i in range(nt): \n", " if imethod =='Lax-Wendroff': \n", " for j in range(1,nx-1):\n", " dQ1 = Q[:,j+1] - Q[:,j-1]\n", " dQ2 = Q[:,j-1] - 2*Q[:,j] + Q[:,j+1]\n", " Qnew[:,j] = Q[:,j] - 0.5*dt/dx*(A @ dQ1)\\\n", " + 1./2.*(dt/dx)**2 * (A @ A) @ dQ2 # Eq. 8.56\n", " \n", " # Absorbing boundary conditions\n", " Qnew[:,0] = Qnew[:,1]\n", " Qnew[:,nx-1] = Qnew[:,nx-2]\n", "\n", " elif imethod == 'upwind': \n", " for j in range(1,nx-1):\n", " dQl = Q[:,j] - Q[:,j-1]\n", " dQr = Q[:,j+1] - Q[:,j]\n", " Qnew[:,j] = Q[:,j] - dt/dx * (Ap @ dQl + Am @ dQr) # Eq. 8.54\n", " \n", " # Absorbing boundary conditions \n", " Qnew[:,0] = Qnew[:,1]\n", " Qnew[:,nx-1] = Qnew[:,nx-2]\n", " else:\n", " raise NotImplementedError\n", "\n", " Q, Qnew = Qnew, Q\n", " \n", " # -------------------------------------- \n", " # Animation plot. Display solution\n", " if not i % isnap: \n", " for l in line1:\n", " l.remove()\n", " del l \n", " for l in line2:\n", " l.remove()\n", " del l \n", " # -------------------------------------- \n", " # Analytical solution (stress i.c.)\n", " Qa[0,:] = 1./2.*(np.exp(-1./sig**2 * (x-x0 + c*i*dt)**2)\\\n", " + np.exp(-1./sig**2 * (x-x0-c*i*dt)**2))\n", "\n", " Qa[1,:] = 1/(2*Z)*(np.exp(-1./sig**2 * (x-x0+c*i*dt)**2)\\\n", " - np.exp(-1./sig**2 * (x-x0-c*i*dt)**2))\n", " \n", " # -------------------------------------- \n", " # Display lines\n", " line1 = ax1.plot(x, Q[0,:], 'k', x, Qa[0,:], 'r--', lw=1.5)\n", " line2 = ax2.plot(x, Q[1,:], 'k', x, Qa[1,:], 'r--', lw=1.5)\n", " plt.legend(iter(line2), ('F. Volume', 'Analytic'))\n", " plt.gcf().canvas.draw()\n"]}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}