{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Finite Element 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", "\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", "This notebook presents a finite element code for the 1D elastic wave equation. Additionally, a solution using finite difference scheme is given for comparison.\n", "\n", "The problem of solving the wave 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", "using the finite element method is done after a series of steps performed on the above equation.\n", "\n", "1) We first obtain a weak form of the wave equation by integrating over the entire physical domain $D$ and at the same time multiplying by some basis $\\varphi_{i}$. \n", "\n", "2) Integration by parts and implementation of the stress-free boundary condition is performed.\n", "\n", "3) We approximate our unknown displacement field $u(x, t)$ by a sum over space-dependent basis functions $\\varphi_i$ weighted by time-dependent coefficients $u_i(t)$.\n", "\n", "\\begin{equation}\n", "u(x,t) \\ \\approx \\ \\overline{u}(x,t) \\ = \\ \\sum_{i=1}^{n} u_i(t) \\ \\varphi_i(x)\n", "\\end{equation}\n", "\n", "4) Utilize the same basis functions used to expand $u(x, t)$ as test functions in the weak form, this is the Galerkin principle.\n", "\n", "5) We can turn the continuous weak form into a system of linear equations by considering the approximated displacement field.\n", "\n", "\\begin{equation}\n", "\\mathbf{M}^T\\partial_t^2 \\mathbf{u} + \\mathbf{K}^T\\mathbf{u} = \\mathbf{f}\n", "\\end{equation}\n", "\n", "6) For the second time-derivative, we use a standard finite-difference approximation. Finally, we arrive at the explicit time extrapolation scheme.\n", "\n", "\\begin{equation}\n", "\\mathbf{u}(t + dt) = dt^2 (\\mathbf{M}^T)^{-1}[\\mathbf{f} - \\mathbf{K}^T\\mathbf{u}] + 2\\mathbf{u} - \\mathbf{u}(t-dt).\n", "\\end{equation}\n", "\n", "where $\\mathbf{M}$ is known as the mass matrix, and $\\mathbf{K}$ the stiffness matrix.\n", "\n", "7) As interpolating functions, we choose interpolants such that $\\varphi_{i}(x_{i}) = 1$ and zero elsewhere. Then, we transform the space coordinate into a local system. According to $\\xi = x \u2212 x_{i}$ and $h_{i} = x_{i+1} \u2212 x_{i}$, we have:\n", "\n", "

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

\n", "\n", "\\begin{equation}\n", " \\varphi_{i}(\\xi) =\n", " \\begin{cases}\n", " \\frac{\\xi}{h_{i-1}} + 1 & \\quad \\text{if} \\quad -h_{i-1} \\le \\xi \\le 0\\\\\n", " 1 + \\frac{\\xi}{h_{i}} & \\quad \\text{if} \\quad 0 \\le \\xi \\le h_{i}\\\\\n", " 0 & \\quad elsewhere\\\\\n", " \\end{cases}\n", "\\end{equation}\n", "\n", "with the corresponding derivatives\n", "\n", "\\begin{equation}\n", " \\partial_{\\xi}\\varphi_{i}(\\xi) =\n", " \\begin{cases}\n", " \\frac{1}{h_{i-1}} & \\quad \\text{if} \\quad -h_{i-1} \\le \\xi \\le 0\\\\\n", " -\\frac{1}{h_{i}} & \\quad \\text{if} \\quad 0 \\le \\xi \\le h_{i}\\\\\n", " 0 & \\quad elsewhere\\\\\n", " \\end{cases}\n", "\\end{equation}\n", "\n", "The figure on the left-hand side illustrates the shape of $\\varphi_{i}(\\xi)$ and $\\partial_{\\xi}\\varphi_{i}(\\xi)$ with varying $h$.\n", "\n", "Code implementation starts with the initialization of a particular setup of our problem. Then, we define the source that introduces perturbations following by initialization of the mass and stiffness matrices. Finally, time extrapolation is done."]}, {"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", "# Basic parameters\n", "nt = 2000 # Number of time steps\n", "vs = 3000 # Wave velocity [m/s] \n", "ro0 = 2500 # Density [kg/m^3]\n", "nx = 1000 # Number of grid points \n", "isx = 500 # Source location [m] \n", "xmax = 10000. # Maximum length\n", "eps = 0.5 # Stability limit\n", "iplot = 20 # Snapshot frequency\n", "\n", "dx = xmax/(nx-1) # calculate space increment\n", "x = np.arange(0, nx)*dx # initialize space coordinates\n", "x = np.transpose(x)\n", "\n", "h = np.diff(x) # Element sizes [m]\n", "\n", "# parameters\n", "ro = x*0 + ro0\n", "mu = x*0 + ro*vs**2\n", "\n", "# time step from stabiity criterion\n", "dt = 0.5*eps*dx/np.max(np.sqrt(mu/ro))\n", "# initialize time axis\n", "t = np.arange(1, nt+1)*dt \n", "\n", "# ---------------------------------------------------------------\n", "# Initialize fields\n", "# ---------------------------------------------------------------\n", "u = np.zeros(nx)\n", "uold = np.zeros(nx)\n", "unew = np.zeros(nx)\n", "\n", "p = np.zeros(nx)\n", "pold = np.zeros(nx)\n", "pnew = np.zeros(nx)"]}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 2. Source time function\n", "\n", "In 1D the propagating signal is an integral of the source time function. As we look for a Gaussian waveform, we initialize the source time function $f(t)$ using the first derivative of a Gaussian function.\n", "\n", "\\begin{equation}\n", "f(t) = -\\dfrac{2}{\\sigma^2}(t - t_0)e^{-\\dfrac{(t - t_0)^2}{\\sigma^2}}\n", "\\end{equation}\n", " \n", "#### Exercise 1 \n", "Initialize a source time function called 'src'. Use $\\sigma = 20 dt$ as Gaussian width, and time shift $t_0 = 3\\sigma$. Then, visualize the source in a given plot."]}, {"cell_type": "code", "execution_count": null, "metadata": {"tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# INITIALIZE THE SOURCE TIME FUCTION HERE!\n", "#################################################################\n", "\n", "# Source vector\n", "f = np.zeros(nx); f[isx:isx+1] = f[isx:isx+1] + 1.\n", "\n", "#################################################################\n", "# PLOT THE SOURCE TIME FUNCTION HERE!\n", "#################################################################"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": [], "solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 3. The Mass Matrix\n", "Having implemented the desired source, now we initialize the mass and stiffness matrices. In general, the mass matrix is given\n", "\n", "\\begin{equation}\n", "M_{ij} = \\int_{D} \\rho \\varphi_i \\varphi_j dx = \\int_{D_{\\xi}} \\rho \\varphi_i \\varphi_j d\\xi\n", "\\end{equation}\n", "\n", "next, the defined basis are introduced and some algebraic treatment is done to arrive at the explicit form of the mass matrix\n", "\n", "#### Exercise 2 \n", "Implement the mass matrix \n", "\n", "\\begin{equation}\n", "M_{ij} = \\frac{\\rho h}{6}\n", " \\begin{pmatrix}\n", " \\ddots & & & & 0\\\\\n", " 1 & 4 & 1 & & \\\\\n", " & 1 & 4 & 1 & \\\\\n", " & & 1 & 4 & 1\\\\\n", " 0 & & & & \\ddots\n", " \\end{pmatrix} \n", "\\end{equation}\n", "\n", "Compute the inverse mass matrix and display your result to visually inspect how it looks like"]}, {"cell_type": "code", "execution_count": null, "metadata": {"lines_to_next_cell": 2, "tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# IMPLEMENT THE MASS MATRIX HERE!\n", "#################################################################\n", "\n", "\n", "#################################################################\n", "# COMPUTE THE INVERSE MASS MATRIX HERE!\n", "#################################################################\n", "\n", "\n", "#################################################################\n", "# DISPLAY THE INVERSE MASS MATRIX HERE!\n", "#################################################################"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {"solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 4. The Stiffness matrix\n", "On the other hand, the general form of the stiffness matrix is\n", "\n", "\\begin{equation}\n", "K_{ij} = \\int_{D} \\mu \\partial_x\\varphi_i \\partial_x\\varphi_j dx = \\int_{D_{\\xi}} \\mu \\partial_\\xi\\varphi_i \\partial_\\xi\\varphi_j d\\xi\n", "\\end{equation} \n", "\n", "at this point, the defined basis are introduced. Again, with the help of some algebraic treatment, we arrive at the explicit form of the stiffness matrix\n", "\n", "#### Exercise 3 \n", "Implement the stiffness matrix \n", "\n", "\\begin{equation}\n", "K_{ij} = \\frac{\\mu}{h}\n", " \\begin{pmatrix}\n", " \\ddots & & & & 0\\\\\n", " -1 & 2 & -1 & & \\\\\n", " &-1 & 2 & -1 & \\\\\n", " & & -1 & 2 & -1\\\\\n", " 0 & & & & \\ddots\n", " \\end{pmatrix} \n", "\\end{equation}\n", "\n", "Display the stiffness matrix to visually inspect how it looks like"]}, {"cell_type": "code", "execution_count": null, "metadata": {"tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# IMPLEMENT THE STIFFNESS MATRIX HERE!\n", "#################################################################\n", "\n", "\n", "#################################################################\n", "# DISPLAY THE STIFFNESS MATRIX HERE!\n", "#################################################################"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": [], "solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### 5. Finite differences matrices\n", "We implement a finite difference scheme in order to compare with the finite elements solution. \n", "\n", "#### Exercise 4 \n", "Implement the finite differences matrices $M$ and $D$. Where $M$ is a diagonal mass matrix containing the inverse densities, and differentiation matrix \n", "\n", "\\begin{equation}\n", "D_{ij} = \\frac{\\mu}{dt^2}\n", " \\begin{pmatrix}\n", " -2 & 1 & & & \\\\\n", " 1 & -2 & 1 & & \\\\\n", " & & \\ddots & & \\\\\n", " & & 1 & -2 & 1\\\\\n", " & & & 1 & -2\n", " \\end{pmatrix} \n", "\\end{equation}\n", "\n", "Display both matrices to visually inspect how they look like"]}, {"cell_type": "code", "execution_count": null, "metadata": {"lines_to_next_cell": 2, "tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# INITIALIZE FINITE DIFFERENCES HERE!\n", "#################################################################\n", "\n", "\n", "#################################################################\n", "# DISPLAY THE DIFFERENCES MATRICES HERE!\n", "#################################################################\n"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {"solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### 6. Finite element solution \n", "\n", "Finally we implement the finite element solution using the computed mass $M$ and stiffness $K$ matrices together with a finite differences extrapolation scheme\n", "\n", "\\begin{equation}\n", "\\mathbf{u}(t + dt) = dt^2 (\\mathbf{M}^T)^{-1}[\\mathbf{f} - \\mathbf{K}^T\\mathbf{u}] + 2\\mathbf{u} - \\mathbf{u}(t-dt).\n", "\\end{equation}"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": [], "lines_to_next_cell": 0}, "outputs": [], "source": ["# Initialize animated plot\n", "# ---------------------------------------------------------------\n", "plt.figure(figsize=(12,4))\n", "\n", "line1 = plt.plot(x, u, 'k', lw=1.5, label='FEM')\n", "line2 = plt.plot(x, p, 'r', lw=1.5, label='FDM')\n", "plt.title('Finite elements 1D Animation', fontsize=16)\n", "plt.ylabel('Amplitude', fontsize=12)\n", "plt.xlabel('x (m)', fontsize=12)\n", "\n", "plt.ion() # set interective mode\n", "plt.show()\n", "\n", "# ---------------------------------------------------------------\n", "# Time extrapolation\n", "# ---------------------------------------------------------------\n", "for it in range(nt):\n", " # --------------------------------------\n", " # Finite Element Method\n", " unew = (dt**2) * Minv @ (f*src[it] - K @ u) + 2*u - uold \n", " uold, u = u, unew\n", " \n", " # --------------------------------------\n", " # Finite Difference Method\n", " pnew = (dt**2) * Mf @ (D @ p + f/dx*src[it]) + 2*p - pold\n", " pold, p = p, pnew\n", " \n", " # -------------------------------------- \n", " # Animation plot. Display both solutions\n", " if not it % iplot:\n", " for l in line1:\n", " l.remove()\n", " del l\n", " for l in line2:\n", " l.remove()\n", " del l\n", " line1 = plt.plot(x, u, 'k', lw=1.5, label='FEM')\n", " line2 = plt.plot(x, p, 'r', lw=1.5, label='FDM')\n", " plt.legend()\n", " plt.gcf().canvas.draw()"]}, {"cell_type": "code", "execution_count": null, "metadata": {"lines_to_next_cell": 2}, "outputs": [], "source": []}], "metadata": {"jupytext": {"encoding": "# -*- coding: utf-8 -*-"}, "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}