{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Finite Differences - Seismometer 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", "* Ashim Rijal ([@ashimrijal](https://github.com/ashimrijal))\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This exercise covers the following aspects:\n", "\n", "* Solving seismometer equation with finite difference method\n", "* Getting familiar with seismometer response function\n", "---"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Basic Equations\n", "\n", "** Please refer to the Exercise 4.19 from the book.**\n", "\n", "The seismometer equation is given by\n", "$$\n", "\\ddot{x} + 2\\epsilon \\dot{x} + \\omega^2_0 x = - \\ddot{u}\n", "$$\n", "\n", "Where,\n", "\n", "$\\epsilon$ is the damping parameter,\n", "\n", "$\\omega_0$ is the eigenfrequency,\n", "\n", "$\\ddot{u}(t)$ is the ground motion by which the seismometer is excited, and\n", "\n", "$x(t)$ is the motion of the seismometer mass.\n", "\n", "We replace the time derivative by centered finite-differentiation\n", "$$\n", "\\dot{x} \\ \\approx \\ \\frac{x (t + \\mathrm{d}t) - x ( t- \\mathrm{d}t)} {2\\mathrm{d}t}\n", "$$\n", "\n", "$$\n", "\\ddot{x} \\ \\approx \\ \\frac{x( t+ \\mathrm{d}t)-2x(t) + x( t- \\mathrm{d}t)} {\\mathrm{d}t^2}\n", "$$\n", "\n", "Now, solving for $x(t+\\mathrm{d}t)$ the extrapolation scheme is\n", "\n", "$$\n", "x(t+\\mathrm{d}t) = \\frac{ - \\ddot{u}(t) \\mathrm{d}t^2 + (2-\\omega^2_0 \\mathrm{d}t^2) x(t) + (\\epsilon \\mathrm{d}t-1) x(t-\\mathrm{d}t)} {(1+\\epsilon \\mathrm{d}t)}\n", "$$\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Exercise\n", "\n", "** Part 1**\n", "\n", "While running the following cells frequency of forcing (the frequency of ground motion) and the damping parameter will be asked to enter. First try using undamped seismometer (i.e. h = 0) for some forcing frequency (eg. 0.1 Hz, 1 Hz, 2Hz, 3Hz, 4Hz, 5Hz, etc.) and interpret the results.\n", "\n", "** Part 2**\n", "\n", "Now try frequency of forcing of your choice (eg. 1 HZ) and try to search for suitable damping parameter (h).\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.**"]}, {"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", "f0 = 1. # eigenfrequency of seismometer (hertz)\n", "w = 2. * np.pi * f0 # in radians per second\n", "dt = .01 # time increment for numerical scheme\n", "isnap = 2 # snapshot frequency for visualization\n", "\n", "# Central frequency of forcing or ground motion (will be asked)\n", "fu0 = 1.0\n", "# Uncomment for interactivity.\n", "# fu0 = float(input('Give frequency of forcing (e.g. f=1 Hz) : '))\n", "\n", "# Damping parameter of seismometer (will be asked)\n", "h = 0.5\n", "# Uncomment for interactivity.\n", "# h = float(input('Give damping parameter (e.g. h=0.5) : '))\n", "\n", "# Initialize ground motion\n", "\n", "# Initialize parameters for ground motion\n", "p = 1. / fu0 # period\n", "nts = int(2. * p / dt) # time steps\n", "uii = np.zeros(nts) # ground motion\n", "t0 = p / dt # time (S)\n", "a = 4. / p # half-width (so called sigma)\n", "\n", "# we use derivative of a Gaussian as our ground motion\n", "for it in range(nts):\n", " t = (it - t0) * dt\n", " uii[it] = -2 * a * t * np.exp(-(a * t) ** 2)\n", "\n", "nt = int(round(5. * 1. / fu0 / dt)) # total number of time steps\n", "src = np.zeros(nt) # initialize an source array of zeros\n", "src[0:len(uii)] = uii # make derivative of gaussian as source\n", "\n", "# End initialization of ground motion\n", "\n", "\n", "# Initial conditions\n", "eps = h * w # damping factor\n", "x = np.zeros(nt)\n", "xnow = 0\n", "xold = 0\n", "x_vector = np.zeros(nt)"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": }, "outputs": [], "source": ["# Extrapolation scheme and the plots\n", "\n", "# Initialization of plots\n", "\n", "# lines1 will plot the seismometer response and lines2 for source function\n", "\n", "lines1 = plt.plot(np.dot((np.arange(1, nt+1)), dt), x_vector[0:nt] /\n", " np.max(np.abs(x[0:nt])), color = \"red\", lw = 1.5, label=\"Seismometer response\")\n", "lines2 = plt.plot(np.zeros(nt), color = 'blue',lw = 1.5, label=\"Ground Acceleration\")\n", "\n", "plt.title(\"At rest\")\n", "plt.axis([0, nt*dt, -1, 1])\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Displacement\")\n", "plt.legend(loc=\"upper right\")\n", "plt.ion()\n", "plt.show()\n", "\n", "# Begin extrapolation and update the plot\n", "\n", "# Extrapolation scheme (Centered finite-difference)\n", "\n", "for i in range(nt):\n", " if i == 0:\n", " xold = xnow\n", "\n", " xnew = (-src[i] * dt ** 2 + (2 - w ** 2 * dt ** 2) * xnow + (eps * dt - 1) * xold) / (1 + eps * dt)\n", "\n", " xold = xnow # for next loop present will be past\n", " xnow = xnew # for next step future will be present\n", " x[i] = xnow\n", " x_vector[i] = x[i]\n", "\n", "# Updating the plots\n", " if not i % isnap:\n", " for l in lines1:\n", " l.remove()\n", " del l\n", " for k in lines2:\n", " k.remove()\n", " del k\n", "\n", " lines1 = plt.plot(np.dot((np.arange (1, i+1)), dt), x_vector[0:i] /\n", " np.max(np.abs (x[0:nt])),color = \"red\",lw = 1.5, label=\"Seismometer response\")\n", " lines2 = plt.plot(np.dot((np.arange (1, i+1)), dt), src[0:i] /\n", " np.max(src[0:nt]), color = 'blue',lw = 1.5, label=\"Ground Acceleration\")\n", " plt.title(\"F0 = 1Hz, SRC = %.2f Hz, h = %.2f \" % (fu0, h))\n", " plt.gcf().canvas.draw()\n", "\n", "plt.ioff()\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}