{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Finite Differences - Grid-Staggering Elastic 1D
\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", "\n", "##### Authors:\n", "* Ashim Rijal ([@ashimrijal](https://github.com/ashimrijal))\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This exercise covers the following aspects:\n", "\n", "* Solving velocity-stress formulation of 1D wave equation with finite difference method\n", "* Understanding the grid-staggering in connection with finite difference solution to the elastic wave equation\n", "---"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Basic Equations\n", "** Please refer to the sections 4.6.2 and 4.6.3 from the book.**\n", "\n", "The 1D wave equation (velocity-stress formulation) as a coupled system of two first-order partial differential equations\n", "\n", "$$\n", "\\rho \\partial_t v = \\partial_x \\sigma + f \n", "$$\n", "$$\n", "\\partial_t \\sigma = \\mu \\partial_x v \n", "$$\n", "\n", "where,\n", "\n", "$\\sigma$ is the stress,\n", "\n", "$\\rho$ is the density,\n", "\n", "$v$ is the velocity,\n", "\n", "$\\mu$ is the shear modulus, and\n", "\n", "$f$ is the source.\n", "\n", "Grid- staggering is the concept in connection with finite-difference solutions to the elastic wave equation. \n", "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.\n", "In grid staggering the following computational scheme is used\n", "\n", "$$\n", "\\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} \\\n", "$$\n", "\n", "$$\n", "\\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}\n", "$$\n", "\n", "The extrapolation scheme becomes\n", "\n", "$$\n", "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 \\\n", "$$\n", "\n", "$$\n", "\\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 \\ \\\n", "$$\n", "\n", "\n", "** Note that in the codes below we do not deal with the index fractions.**"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Exercise\n", "First understand the codes below and run the simulation. \n", "\n", "Then, improve the result using (4-point operator) for 1st derivative.\n", "\n", "**Message: Once you become familiar with all the codes below you can go to the Cell tab on the toolbar and click Run All.**\n"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": , "lines_to_next_cell": 2}, "outputs": [], "source": ["# Configuration step (Please run it before the simulation code!)\n", "\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", "matplotlib.rcParams['figure.facecolor'] = 'w' # remove grey background"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": [], "lines_to_next_cell": 2}, "outputs": [], "source": ["# Initialization of parameters\n", "\n", "# Simple finite difference solver\n", "# Elastic wave equation\n", "# 1-D regular staggered grid\n", "\n", "# Basic parameters\n", "nt = 1300 # number of time steps\n", "nx = 1000 # number of grid points in x\n", "c0 = 4500 # velocity (m/sec) (shear wave)\n", "eps = 0.8 # stability limit\n", "isnap = 2 # snapshot frequency\n", "isx = round(nx/2) # source location\n", "f0 = 1/15 # frequency (Hz)\n", "xmax = 1000000. # maximum range (m)\n", "rho0 = 2500. # density (kg/m**3)\n", "mu0 = rho0*c0**2. # shear modulus (Pa)\n", "nop = 4 # number of operator either 2 or 4\n", "\n", "dx = xmax / (nx-1) # calculate space increment (m)\n", "x = (np.arange(nx)*dx) # initialize space coordinates \n", "dt = eps * dx/c0 # calculate time step from stability criterion(s)\n", "\n", "# Source time function\n", "t = (np.arange(0,nt) * dt) # initialize time axis\n", "T0 = 1. / f0 # period\n", "a = 4. / T0 # half-width (so called sigma)\n", "t0 = T0 / dt\n", "tmp = np.zeros(nt)\n", "for it in range(nt):\n", " t = (it - t0) * dt\n", " tmp[it] = -2 * a * t * np.exp(-(a * t) ** 2) # derivative of Gaussian (so called sigma)\n", "src = np.zeros(nt) # source\n", "src[0:len(tmp)] = tmp\n", "lam = c0 * T0 # wavelength"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": [], "lines_to_next_cell": 0}, "outputs": [], "source": ["# Extrapolation scheme and the plots\n", "\n", "# Initialization of plot\n", "\n", "# Initialization of fields\n", "v = np.zeros(nx) # velocity\n", "vnew = v\n", "dv = v\n", "\n", "s = np.zeros(nx) # stress\n", "snew = s\n", "ds = s\n", "\n", "mu = np.zeros(nx) # shear modulus\n", "rho = mu \n", "rho = rho + rho0\n", "mu = mu + mu0\n", "\n", "# Print setup parameters\n", "print(\"rho =\", rho0, \", f_dom =\", f0, \", stability limit =\", eps, \", n_lambda\", (lam/dx))\n", "\n", "# Initialize the plot\n", "title = \"FD Elastic 1D staggered grid\"\n", "fig = plt.figure(figsize=(12,8))\n", "ax1 = fig.add_subplot(2, 1, 1)\n", "ax2 = fig.add_subplot(2, 1, 2)\n", "line1 = ax1.plot(x, v, color = \"red\", lw = 1.5)\n", "line2 = ax2.plot(x, s, color = \"blue\", lw = 1.5)\n", "ax1.set_ylabel('velocity (m/s)')\n", "ax2.set_xlabel('x (m)')\n", "ax2.set_ylabel('stress (Pa)')\n", "plt.ion()\n", "plt.show()\n", "\n", "# Begin extrapolation and update the plot\n", "for it in range (nt):\n", " \n", " # Stress derivative\n", " for i in range (2, nx-2):\n", " ds[i] = (0.0416666 * s[i-1] - 1.125 * s[i] + 1.125 * s[i+1] - 0.0416666 * s[i+2]) / (dx)\n", " \n", " # Velocity extrapolation\n", " v = v + dt * ds / rho\n", " \n", " # Add source term at isx\n", " v[isx] = v[isx] + dt * src[it] / (dt * rho[isx])\n", " \n", " # Velocity derivative\n", " for i in range (2, nx-2):\n", "\n", " dv[i] = (0.0416666 * v[i-2] - 1.125 * v[i-1] + 1.125 * v[i] - 0.0416666 * v[i+1]) / (dx)\n", " \n", " # Stress extrapolation\n", " s = s + dt * mu * dv\n", " \n", " # Updating the plots\n", " if not it % isnap: \n", " for l in line1:\n", " l.remove()\n", " del l \n", " for l in line2:\n", " l.remove()\n", " del l \n", " line1 = ax1.plot(x, v, color = \"red\", lw = 1.5)\n", " line2 = ax2.plot(x, s, color = \"blue\", lw = 1.5)\n", " \n", " ax1.set_title(title + \", time step: %i\" % (it)) \n", " plt.gcf().canvas.draw()\n", "\n", "plt.ioff()\n", "plt.show()\n"]}, {"cell_type": "code", "execution_count": null, "metadata": {"lines_to_next_cell": 2}, "outputs": [], "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}