{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Finite Differences Method - Acoustic Waves in 2D
\n", "
\n", "
\n", "
"]}, {"cell_type": "markdown", "metadata": {}, "source": ["\n", "

\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", "\n", "##### Authors:\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))\n", "* Florian W\u00f6lfl ([@flo-woelfl](https://github.com/flo-woelft))\n", "* Lion Krischer ([@krischer](https://github.com/krischer))\n", "\n", "---"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This exercise covers the following aspects:\n", "\n", "* presenting you with an implementation of the 2D acoustic wave equation \n", "* allowing you to explore the benefits of using high-order finite-difference operators\n", "* understanding the concepts of stability (Courant criterion)\n", "* exploration of numerical dispersion and numerical grid anisotropy\n", "* changing the earth model and exploring some effects of structural heterogeneities (e.g., fault zones)\n", "\n", "---"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Basic Equations\n", "\n", "The acoustic wave equation in 2D is \n", "$$\n", "\\ddot{p}(x,z,t) \\ = \\ c(x,z)^2 (\\partial_x^2 p(x,z,t) + \\partial_z^2 p(x,z,t)) \\ + s(x,z,t)\n", "$$\n", "\n", "and we replace the time-dependent (upper index time, lower indices space) part by\n", "\n", "$$\n", " \\frac{p_{j,k}^{n+1} - 2 p_{j,k}^n + p_{j,k}^{n-1}}{\\mathrm{d}t^2} \\ = \\ c_j^2 ( \\partial_x^2 p + \\partial_z^2 p) \\ + s_{j,k}^n\n", "$$\n", "\n", "solving for $p_{j,k}^{n+1}$. \n", "The extrapolation scheme is\n", "$$\n", "p_{j,k}^{n+1} \\ = \\ c_j^2 \\mathrm{d}t^2 \\left[ \\partial_x^2 p + \\partial_z^2 p \\right]\n", "+ 2p_{j,k}^n - p_{j,k}^{n-1} + \\mathrm{d}t^2 s_{j,k}^n\n", "$$\n", "The space derivatives are determined by \n", "\n", "$$\n", "\\partial_x^2 p \\ = \\ \\frac{p_{j+1,k}^{n} - 2 p_{j,k}^n + p_{j-1,k}^{n}}{\\mathrm{d}x^2}\n", "$$\n", "$$\n", "\\partial_z^2 p \\ = \\ \\frac{p_{j,k+1}^{n} - 2 p_{j,k}^n + p_{j,k-1}^{n}}{\\mathrm{d}z^2} \n", "$$\n", "\n", "---"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercises\n", "\n", "### 1. Getting started\n", "Before you start it is good practice to immediately make a copy of the original notebook (e.g., X_orig.ipynb). \n", "Run the simulation code. Relate the time extrapolation loop with the numerical algorithm we developed in the course. Understand the input parameters for the simulation and the plots that are generated. Modify source and receiver locations and observe the effects on the seismograms. \n", "\n", "### 2. Stability\n", "Introduce a new parameter (e.g., eps) and calculate the Courant criterion. Determine numerically the stability limit of the code as accurately as possible by increasing the time step. Print the max value of the pressure field at each time step and observe the evolution of it in the case of stable and unstable simulations. \n", "(Hint: The Courant criterion is defined as $eps = (velocity * dt) / dx$ . With this information you can calculate the maximum possible, stable time step. )\n", "\n", "### 3. High-order operators\n", "Extend the code by adding the option to use a 5-point difference operator (see problem 1 of exercise sheet). Compare simulations with the 3-point and 5-point operator. Is the stability limit still the same? Make it an option to change between 3-pt and 5-pt operator. Estimate the number of points per wavelength and investigate the accuracy of the simulation by looking for signs of numerical dispersion in the resulting seismograms. The 5-pt weights are: $[-1/12, 4/3, -5/2, 4/3, -1/12]/dx^2$. \n", "\n", "### 4. Numerical anisotropy\n", "Increase the frequency of the wavefield by varying f0. Investigate the angular dependence of the wavefield. Why does the wavefield look anisotropic? Which direction is the most accurate and why? What happens if you set the source time function to a spike (zero everywhere except one element with value 1). \n", "\n", "### 5. Heterogeneous models\n", "Now let us explore the power of the finite-difference method by varying the internal structure of the model. Here we can only modify the velocity c that can vary at each grid point (any restrictions?). Here are some suggestions. Investigate the influence of the structure by analysing the snapshots and the seismograms. \n", "\n", "* Add a low(high) velocity layer near the surface. Put the source at zs=2.\n", "* Add a vertical low velocity zone (fault zone) of a certain width (e.g. 10 grid points), and discuss the resulting wavefield\n", "* Simulate topography by setting the pressure to 0 above the surface. Use a Gaussian hill shape or a random topography.\n", "* etc. \n", "\n", "### 6. Source-receiver reciprocity\n", " Initialize a strongly heterogeneous 2D velocity model of your choice and simulate waves propagating from an internal source point ($x_s, z_s$) to an internal receiver ($x_r, z_r$). Show that by reversing source and receiver you obtain the same seismogram.\n", "\n", "### 7. Time reversal\n", "Time reversal. Define in an arbitrary 2D velocity model a source at the centre of the domain an a receiver circle at an appropriate distance around the source. Simulate a wavefield, record it at the receiver ring and store the results. Reverse the synthetic seismograms and inject the as sources at the receiver points. What happens? Do you know examples where this principle is used? \n", "\n", "---"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": }, "outputs": [], "source": ["# This is a configuration step for the exercise. Please run it before the simulation code! \n", "%matplotlib notebook\n", "import numpy as np\n", "import matplotlib.pyplot as plt"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Below is the 2D acoustic simulation code:"]}, {"cell_type": "code", "execution_count": null, "metadata": {"tags": ["exercise"]}, "outputs": [], "source": ["# Simple finite difference solver \n", "# Acoustic wave equation p_tt = c^2 p_xx + src\n", "# 2-D regular grid\n", "\n", "nx = 200 # grid points in x\n", "nz = 200 # grid points in z\n", "nt = 750 # number of time steps\n", "dx = 10.0 # grid increment in x\n", "dt = 0.001 # Time step\n", "c0 = 3000.0 # velocity (can be an array)\n", "isx = nx // 2 # source index x\n", "isz = nz // 2 # source index z\n", "ist = 100 # shifting of source time function\n", "f0 = 100.0 # dominant frequency of source (Hz)\n", "isnap = 10 # snapshot frequency\n", "T = 1.0 / f0 # dominant period\n", "nop = 3 # length of operator\n", "\n", "# Model type, available are \"homogeneous\", \"fault_zone\",\n", "# \"surface_low_velocity_zone\", \"random\", \"topography\",\n", "# \"slab\"\n", "model_type = \"homogeneous\"\n", "\n", "# Receiver locations\n", "irx = np.array([60, 80, 100, 120, 140])\n", "irz = np.array([5, 5, 5, 5, 5])\n", "seis = np.zeros((len(irx), nt))\n", "\n", "# Initialize pressure at different time steps and the second\n", "# derivatives in each direction\n", "p = np.zeros((nz, nx))\n", "pold = np.zeros((nz, nx))\n", "pnew = np.zeros((nz, nx))\n", "pxx = np.zeros((nz, nx))\n", "pzz = np.zeros((nz, nx))\n", "\n", "# Initialize velocity model\n", "c = np.zeros((nz, nx))\n", "\n", "if model_type == \"homogeneous\":\n", " c += c0\n", "elif model_type == \"fault_zone\":\n", " c += c0\n", " c[:, nx // 2 - 5: nx // 2 + 5] *= 0.8 \n", "elif model_type == \"surface_low_velocity_zone\":\n", " c += c0\n", " c[1:10,:] *= 0.8\n", "elif model_type == \"random\":\n", " pert = 0.4\n", " r = 2.0 * (np.random.rand(nz, nx) - 0.5) * pert\n", " c += c0 * (1 + r) \n", "elif model_type == \"topography\":\n", " c += c0\n", " c[0 : 10, 10 : 50] = 0 \n", " c[0 : 10, 105 : 115] = 0 \n", " c[0 : 30, 145 : 170] = 0\n", " c[10 : 40, 20 : 40] = 0\n", " c[0 : 15, 50 : 105] *= 0.8 \n", "elif model_type == \"slab\":\n", " c += c0\n", " c[110 : 125, 0 : 125] = 1.4 * c0\n", " for i in range(110, 180):\n", " c[i , i-5 : i + 15 ] = 1.4 * c0\n", "else:\n", " raise NotImplementedError\n", " \n", "cmax = c.max()\n", "\n", "# Source time function Gaussian, nt + 1 as we loose the last one by diff\n", "src = np.empty(nt + 1)\n", "for it in range(nt):\n", " src[it] = np.exp(-1.0 / T ** 2 * ((it - ist) * dt) ** 2)\n", "# Take the first derivative\n", "src = np.diff(src) / dt\n", "src[nt - 1] = 0\n", "\n", "v = max([np.abs(src.min()), np.abs(src.max())])\n", "# Initialize animated plot\n", "image = plt.imshow(pnew, interpolation='nearest', animated=True,\n", " vmin=-v, vmax=+v, cmap=plt.cm.RdBu)\n", "\n", "# Plot the receivers\n", "for x, z in zip(irx, irz):\n", " plt.text(x, z, '+')\n", "\n", "plt.text(isx, isz, 'o')\n", "plt.colorbar()\n", "plt.xlabel('ix')\n", "plt.ylabel('iz')\n", "\n", "\n", "plt.ion()\n", "plt.show()\n", "\n", "\n", "# required for seismograms\n", "ir = np.arange(len(irx))\n", "\n", "\n", "#################################################\n", "# CALCULATE AND PRINT COURANT CRITERION HERE\n", "#################################################\n", "\n", "\n", "# Time extrapolation\n", "for it in range(nt):\n", " if nop==3:\n", " # calculate partial derivatives, be careful around the boundaries\n", " for i in range(1, nx - 1):\n", " pzz[:, i] = (p[:, i + 1] - 2 * p[:, i] + p[:, i - 1]) / dx ** 2\n", " for j in range(1, nz - 1):\n", " pxx[j, :] = (p[j - 1, :] - 2 * p[j, :] + p[j + 1, :]) / dx ** 2\n", "\n", " if nop==5:\n", " #################################################\n", " # IMPLEMENT 5 POINT CENTERED DIFFERENCES HERE!\n", " #################################################\n", " pass\n", "\n", " # Time extrapolation\n", " pnew = 2 * p - pold + dt ** 2 * c ** 2 * (pxx + pzz)\n", " # Add source term at isx, isz\n", " pnew[isz, isx] = pnew[isz, isx] + src[it]\n", "\n", " # Plot every isnap-th iteration\n", " if it % isnap == 0: # you can change the speed of the plot by increasing the plotting interval\n", " \n", " plt.title(\"Max P: %.2f\" % p.max())\n", " image.set_data(pnew)\n", " plt.gcf().canvas.draw()\n", "\n", " pold, p = p, pnew\n", "\n", " # Save seismograms\n", " seis[ir, it] = p[irz[ir], irx[ir]]"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": [], "tags": ["solution"]}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["The cell below allows you to plot source time function, seismic velocites, and the resulting seismograms in windows inside the notebook. Remember to rerun after you simulated again!"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": , "lines_to_next_cell": 0}, "outputs": [], "source": ["# Plot the source time function and the seismograms \n", "\n", "plt.ioff()\n", "plt.figure(figsize=(12, 12))\n", "\n", "plt.subplot(221)\n", "time = np.arange(nt) * dt\n", "plt.plot(time, src)\n", "plt.title('Source time function')\n", "plt.xlabel('Time (s) ')\n", "plt.ylabel('Source amplitude ')\n", "\n", "plt.subplot(222)\n", "ymax = seis.ravel().max() \n", "for ir in range(len(seis)):\n", " plt.plot(time, seis[ir, :] + ymax * ir)\n", " plt.xlabel('Time (s)')\n", " plt.ylabel('Amplitude')\n", "\n", "plt.subplot(223)\n", "ymax = seis.ravel().max()\n", "for ir in range(len(seis)):\n", " plt.plot(time, seis[ir, :] + ymax * ir)\n", " plt.xlabel('Time (s)')\n", " plt.ylabel('Amplitude')\n", "\n", "plt.subplot(224)\n", "# The velocity model is influenced by the Earth model above\n", "plt.title('Velocity Model')\n", "plt.imshow(c)\n", "plt.xlabel('ix')\n", "plt.ylabel('iz')\n", "plt.colorbar()\n", "\n", "plt.show()"]}, {"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}