{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
The Chebyshev Pseudospectral Method - Elastic Waves in 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", "* David Vargas ([@dvargas](https://github.com/davofis))\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))\n", "---"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Basic Equations\n", "This notebook presents the numerical solution for the 1D elastic wave equation using the Chebyshev Pseudospectral Method. We depart from the equation \n", "\n", "\\begin{equation}\n", "\\rho(x) \\partial_t^2 u(x,t) = \\partial_x (\\mu(x) \\partial_x u(x,t)) + f(x,t),\n", "\\end{equation}\n", "\n", "and use a standard 3-point finite-difference operator to approximate the time derivatives. Then, the displacement field is extrapolated as\n", "\n", "\\begin{equation}\n", "\\rho_i\\frac{u_{i}^{j+1} - 2u_{i}^{j} + u_{i}^{j-1}}{dt^2}= \\partial_x (\\mu(x) \\partial_x u(x,t))_{i}^{j} + f_{i}^{j}\n", "\\end{equation}\n", "\n", "An alternative way of performing space derivatives of a function defined on the Chebyshev collocation points is to define a derivative matrix $D_{ij}$\n", "\n", "\\begin{equation}\n", "D_{ij} =\n", " \\begin{cases}\n", " -\\frac{2 N^2 + 1}{6} \\hspace{1.5cm} \\text{for i = j = N}\\\\\n", " -\\frac{1}{2} \\frac{x_i}{1-x_i^2} \\hspace{1.5cm} \\text{for i = j = 1,2,...,N-1}\\\\\n", " \\frac{c_i}{c_j} \\frac{(-1)^{i+j}}{x_i - x_j} \\hspace{1.5cm} \\text{for i $\\neq$ j = 0,1,...,N}\n", " \\end{cases}\n", "\\end{equation}\n", "\n", "where $N+1$ is the number of Chebyshev collocation points $\\ x_i = cos(i\\pi / N)$, $\\ i=0,...,N$ and the $c_i$ are given as\n", "\n", "$$c_i = 2 \\hspace{1.5cm} \\text{for i = 0 or N}$$\n", "$$c_i = 1 \\hspace{1.5cm} \\text{otherwise}$$\n", "\n", "This differentiation matrix allows us to write the derivative of the function $f_i = f(x_i)$ (possibly depending on time) simply as\n", "\n", "$$\\partial_x u_i = D_{ij} \\ u_j$$\n", "\n", "where the right-hand side is a matrix-vector product, and the Einstein summation convention applies."]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": , "lines_to_next_cell": 2}, "outputs": [], "source": ["# This is a configuration step for the exercise. Please run it before calculating the derivative!\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", "from ricker import ricker "]}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 1. Chebyshev derivative method\n", "\n", "#### Exercise\n", "Define a python function call \"get_cheby_matrix(nx)\" that initializes the Chebyshev derivative matrix $D_{ij}$, call this function and display the Chebyshev derivative matrix. "]}, {"cell_type": "code", "execution_count": null, "metadata": {"lines_to_next_cell": 2, "tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# IMPLEMENT THE CHEBYSHEV DERIVATIVE MATRIX METHOD HERE!\n", "################################################################# \n", "\n", "# Call the chebyshev differentiation matrix\n", "# ---------------------------------------------------------------\n", "#D_ij = \n", "\n", "# ---------------------------------------------------------------\n", "# Display Differentiation Matrix\n", "# ---------------------------------------------------------------"]}, {"cell_type": "code", "execution_count": null, "metadata": {"solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### 2. Initialization of setup"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# Basic parameters\n", "# ---------------------------------------------------------------\n", "#nt = 5000 # number of time steps\n", "tmax = 0.0006\n", "eps = 1.4 # stability limit\n", "isx = 100\n", "lw = 0.7\n", "ft = 10\n", "f0 = 100000 # dominant frequency\n", "iplot = 20 # Snapshot frequency\n", "\n", "# material parameters\n", "rho = 2500.\n", "c = 3000.\n", "mu = rho*c**2\n", "\n", "# space domain\n", "nx = 100 # number of grid points in x 199\n", "xs = np.floor(nx/2) # source location\n", "xr = np.floor(nx*0.8)\n", "x = np.zeros(nx+1) \n", "\n", "# initialization of pressure fields\n", "p = np.zeros(nx+1) \n", "pnew = np.zeros(nx+1)\n", "pold = np.zeros(nx+1)\n", "d2p = np.zeros(nx+1) \n", "\n", "for ix in range(0,nx+1):\n", " x[ix] = np.cos(ix * np.pi / nx) \n", "dxmin = min(abs(np.diff(x)))\n", "dxmax = max(abs(np.diff(x)))\n", "\n", "dt = eps*dxmin/c # calculate time step from stability criterion\n", "nt = int(round(tmax/dt))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### 3. Source Initialization "]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# source time function\n", "# ---------------------------------------------------------------\n", "t = np.arange(1, nt+1)*dt # initialize time axis\n", "T0 = 1./f0\n", "tmp = ricker(dt, T0)\n", "isrc = tmp\n", "tmp = np.diff(tmp)\n", "src = np.zeros(nt) \n", "src[0:np.size(tmp)] = tmp\n", "\n", "#spatial source function\n", "# ---------------------------------------------------------------\n", "sigma = 1.5*dxmax\n", "x0 = x[int(xs)]\n", "sg = np.exp(-1/sigma**2*(x-x0)**2)\n", "sg = sg/max(sg) "]}, {"cell_type": "markdown", "metadata": {}, "source": ["### 4. Time Extrapolation\n", "\n", "Now we time extrapolate using the previously defined get_cheby_matrix(nx) method to call the differentiation matrix. The discrete values of the numerical simulation are indicated by dots in the animation, they represent the Chebyshev collocation points. Observe how the wavefield near the domain center is less dense than towards the boundaries."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# Initialize animated plot\n", "# ---------------------------------------------------------------\n", "plt.figure(figsize=(10,6))\n", "line = plt.plot(x, p, 'k.', lw=2)\n", "plt.title('Chebyshev Method - 1D Elastic wave', size=16)\n", "plt.xlabel(' x(m)', size=14)\n", "plt.ylabel(' Amplitude ', size=14)\n", "\n", "plt.ion() # set interective mode\n", "plt.show()\n", "# ---------------------------------------------------------------\n", "# Time extrapolation\n", "# ---------------------------------------------------------------\n", "# Differentiation matrix\n", "D = get_cheby_matrix(nx)\n", "for it in range(nt):\n", " # Space derivatives\n", " dp = np.dot(D, p.T)\n", " dp = mu/rho * dp\n", " dp = D @ dp\n", " \n", " # Time extrapolation \n", " pnew = 2*p - pold + np.transpose(dp) * dt**2\n", " \n", " # Source injection\n", " pnew = pnew + sg*src[it]*dt**2/rho\n", " \n", " # Remapping\n", " pold, p = p, pnew\n", " p = 0; p[nx] = 0 # set boundaries pressure free \n", "\n", " # -------------------------------------- \n", " # Animation plot. Display solution\n", " if not it % iplot: \n", " for l in line:\n", " l.remove()\n", " del l \n", " \n", " # -------------------------------------- \n", " # Display lines\n", " line = plt.plot(x, p, 'k.', lw=1.5)\n", " plt.gcf().canvas.draw()"]}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}