{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Lamb's problem
\n", "
\n", "
\n", "
"]}, {"cell_type": "markdown", "metadata": {}, "source": ["

\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))\n", "\n", "(Based on the original code by Lane Johnson)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Basic Equations\n", "\n", "The fundamental analytical solution to the three-dimensional Lamb\u2019s problem, the problem of determining the elastic disturbance resulting from a point force in a homogeneous half space, is implemented in this Ipython notebook. This solution provides fundamental information used as benchmark for comparison with entirely numerical solutions. A setup of the fundamental problem is illustrated below. The figure on the right hand side is published in  (Figure 1. System of coordinates)\n", "\n", "

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

\n", "\n", "Simulations of 3D elastic wave propagation need to be validated by the use of analytical solutions. In order to evaluate how healthy a numerical solution is, one may recreate conditions for which analytical solutions exist with the aim of reproducing and compare the different results.\n", "\n", "We which to find the displacement wavefield $\\mathbf{u}(\\mathbf{x},t)$ at some distance $\\mathbf{x}$ from a seismic source with $\\mathbf{F} = f_1\\mathbf{\\hat{x}_1} + f_2\\mathbf{\\hat{x}_2} + f_3\\mathbf{\\hat{x}_3}$.\n", "\n", "For a uniform elastic material and a Cartesian co-ordinate system the equation for the conservation of linear momentum can be written\n", "\n", "\\begin{align*}\n", "\\rho(x) \\frac{\\partial^2}{\\partial t^2} \\mathbf{u(\\mathbf{x},t)} = (\\lambda + \\mu)\\nabla(\\nabla\\mathbf{u(\\mathbf{x},t)}) + \\mu\\nabla^2 \\mathbf{u(\\mathbf{x},t)} + \\mathbf{f(\\mathbf{x},t)}\n", "\\end{align*}\n", "\n", "We will consider the case where the source function is localized in both time and space\n", "\n", "\\begin{align*}\n", "\\mathbf{f(\\mathbf{x},t)} = (f_1\\mathbf{\\hat{x}_1} + f_2\\mathbf{\\hat{x}_2} + f_3\\mathbf{\\hat{x}_3})\\delta(x_1 - x^{'}_{1})\\delta(x_2 - x^{'}_{2})\\delta(x_3 - x^{'}_{3})\\delta(t - t^{'})\n", "\\end{align*}\n", "\n", "For such a source we will refer to the displacement solution as a Green\u2019s function, and use the standard notation\n", "\n", "\\begin{align*}\n", "\\mathbf{u(\\mathbf{x},t)} = g_1(\\mathbf{x},t;\\mathbf{x^{'}},t^{'})\\mathbf{\\hat{x}_1} + g_2(\\mathbf{x},t;\\mathbf{x^{'}},t^{'})\\mathbf{\\hat{x}_2} + g_3(\\mathbf{x},t;\\mathbf{x^{'}},t^{'})\\mathbf{\\hat{x}_3}\n", "\\end{align*}\n", "\n", "The complete solution is found after applying the Laplace transform to the elastic wave equation, implementing the stress-free boundary condition, defining some transformations, and performing some algebraic manoeuvres. Then, the Green's function at the free surface is given:\n", "\n", "\\begin{align*}\n", "\\begin{split}\n", "\\mathbf{G}(x_1,x_2,0,t;0,0,x^{'}_{3},0) & = \\dfrac{1}{\\pi^2\\mu r} \\dfrac{\\partial}{\\partial t}\\int_{0}^{((t/r)^2 - \\alpha^{-2})^{1/2}}\\mathbf{H}(t-r/\\alpha)\\mathbb{R}[\\eta_\\alpha\\sigma^{-1}((t/r)^2 - \\alpha^{-2} - p^2)^{-1/2}\\mathbf{M}(q,p,0,t,x^{'}_{3})\\mathbf{F}] dp \\\\\n", " & + \\dfrac{1}{\\pi^2\\mu r} \\dfrac{\\partial}{\\partial t}\\int_{0}^{p_2}\\mathbf{H}(t-t_2)\\mathbb{R}[\\eta_\\beta\\sigma^{-1}((t/r)^2 - \\beta^{-2} - p^2)^{-1/2}\\mathbf{N}(q,p,0,t,x^{'}_{3})\\mathbf{F}] dp\n", "\\end{split}\n", "\\end{align*}\n", "\n", "Details on the involved terms are found in the original paper . The Green's $\\mathbf{G}$ function consist of three components of displacement evolving from the application of three components of force $\\mathbf{F}$. If we assume that each component of $\\mathbf{F}$ provokes three components of displacement, then $\\mathbf{G}$ is composed by nine independent components that correspond one to one to the matrices $\\mathbf{M}$ and $\\mathbf{N}$. Without losing generality it is shown that among them four are equal zero, and we end up only with five possible components. \n", "\n", "

\n", "  Eduardo Kausel - Lamb's problem at its simplest, 2012

\n", " \n", "

\n", "  Lane R. Johnson - Green\u2019s Function for Lamb\u2019s Problem, 1974

\n", "\n"]}, {"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.pyplot as plt\n", "import os\n", "from ricker import ricker\n", "\n", "# Show the plots in the Notebook.\n", "plt.switch_backend(\"nbagg\")"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": }, "outputs": [], "source": ["# Compile the source code (needs gfortran!)\n", "!rm -rf lamb.exe output.txt\n", "!gfortran canhfs.for -fcheck=all -o lamb.exe"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Calling the original FORTRAN code"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# Initialization of setup:\n", "# Figure 4 in Lane R. Johnson - Green\u2019s Function for Lamb\u2019s Problem, 1974\n", "# is reproduced when the following parameters are given\n", "# -----------------------------------------------------------------------------\n", "r = 10.0 # km\n", "vp = 8.0 # P-wave velocity km/s\n", "vs = 4.62 # s-wave velocity km/s\n", "rho = 3.3 # Density kg/m^3 \n", "nt = 512 # Number of time steps\n", "dt = 0.01 # Time step s\n", "h = 0.2 # Source position km (0.01 to reproduce Fig 2.16 of the book)\n", "ti = 0.0 # Initial time s\n", "\n", "var = [vp, vs, rho, nt, dt, h, r, ti]\n", "\n", "# -----------------------------------------------------------------------------\n", "# Execute fortran code\n", "# -----------------------------------------------------------------------------\n", "with open('input.txt', 'w') as f:\n", " for i in var:\n", " print(i, file=f, end=' ') # Write input for fortran code\n", "\n", "f.close()\n", "\n", "os.system(\"./lamb.exe\") # Code execution\n", "\n", "# -----------------------------------------------------------------------------\n", "# Load the solution\n", "# -----------------------------------------------------------------------------\n", "G = np.genfromtxt('output.txt')\n", "\n", "u_rx = G[:,0] # Radial displacement owing to horizontal load\n", "u_tx = G[:,1] # Tangential displacement due to horizontal load\n", "u_zx = G[:,2] # Vertical displacement owing to horizontal load\n", "\n", "u_rz = G[:,3] # Radial displacement owing to a vertical load\n", "u_zz = G[:,4] # Vertical displacement owing to vertical load\n", "\n", "t = np.linspace(dt, nt*dt, nt) # Time axis\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Visualization of the Green's function"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": }, "outputs": [], "source": ["# Plotting\n", "# -----------------------------------------------------------------------------\n", "seis = [u_rx, u_tx, u_zx, u_rz, u_zz] # Collection of seismograms\n", "labels = ['$u_{rx}(t) [cm]$','$u_{tx}(t)[cm]$','$u_{zx}(t)[cm]$','$u_{rz}(t)[cm]$','$u_{zz}(t)[cm]$']\n", "cols = ['b','r','k','g','c']\n", "\n", "# Initialize animated plot\n", "fig = plt.figure(figsize=(12,8), dpi=80)\n", "\n", "fig.suptitle(\"Green's Function for Lamb's problem\", fontsize=16)\n", "\n", "plt.ion() # set interective mode\n", "plt.show()\n", "\n", "for i in range(5): \n", " st = seis[i]\n", " ax = fig.add_subplot(2, 3, i+1)\n", " ax.plot(t, st, lw = 1.5, color=cols[i]) \n", " ax.set_xlabel('Time(s)')\n", " ax.text(0.8*nt*dt, 0.8*max(st), labels[i], fontsize=16)\n", " plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n", " \n", " ax.spines['left'].set_position('zero')\n", " ax.spines['right'].set_color('none')\n", " ax.spines['bottom'].set_position('zero')\n", " ax.spines['top'].set_color('none')\n", " ax.spines['left'].set_smart_bounds(True)\n", " ax.spines['bottom'].set_smart_bounds(True)\n", " ax.xaxis.set_ticks_position('bottom')\n", " ax.yaxis.set_ticks_position('left')\n", "\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {"solution2": "hidden", "solution2_first": true}, "source": ["### Convolution\n", "\n", "Let $S(t)$ be a general source time function, then the displacent seismogram is given in terms of the Green's function $G$ via\n", "\n", "\\begin{equation}\n", "u(\\mathbf{x},t) = G(\\mathbf{x},t; \\mathbf{x}',t') \\ast S(t)\n", "\\end{equation}\n", "\n", "#### Exercise\n", "Compute the convolution of the source time function 'ricker' with the Green's function of a Vertical displacement due to vertical loads. Plot the resulting displacement."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# call the source time function\n", "T = 1/5 # Period\n", "src = ricker(dt,T)\n", "\n", "# Normalize source time function\n", "src = src/max(src)\n", "\n", "# Initialize source time function\n", "f = np.zeros(nt)\n", "f[0:int(2 * T/dt)] = src"]}, {"cell_type": "code", "execution_count": null, "metadata": {"tags": ["exercise"]}, "outputs": [], "source": ["#################################################################\n", "# COMPUTE THE CONVOLUTION HERE!\n", "#################################################################\n", "\n", "\n", "#################################################################\n", "# PLOT THE SEISMOGRAMS HERE!\n", "#################################################################"]}, {"cell_type": "code", "execution_count": null, "metadata": {"solution2": "hidden", "tags": ["solution"]}, "outputs": [], "source": []}], "metadata": {"jupytext": {"encoding": "# -*- coding: utf-8 -*-"}, "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}