01_Introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](figs/COMPGEOP_TITLE.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>Module 1: Introduction to Geophysical Computing</center></h1>\n",
    "\n",
    "\n",
    "For the modern professional geophysicists, the use of computers and algorithms is ubiquitous in all aspects of the daily routine.  Whether using a \"canned\" commercial geophysical software application to or developing new code to address a particular task, employing/implementing some sort of computational algorithm is almost always at the heart of any geophysical task.\n",
    "\n",
    "## The challenge of canned software\n",
    "\n",
    "One of the key challenges geophysicsts face when using commercial geophysical software is that it is often difficult to determine what is going on \"under the hood\".  In many instances it may seem that geophysical results are found in the following way:\n",
    "\n",
    "$$ \n",
    "\\boxed{\\begin{array}{c}\\mathrm{Field}\\\\\\mathrm{Geophysical}\\\\\\mathrm{Data}\\end{array}}\n",
    "\\Longrightarrow \n",
    "\\boxed{\\begin{array}{c}\\mathrm{Magic\\,Black}\\\\\\mathrm{Box\\,of\\, Geophysical}\\\\\\mathrm{Data\\,Processing}\\end{array}}\n",
    "\\Longrightarrow \n",
    "\\boxed{\\begin{array}{c}\\mathrm{Top-notch}\\\\\\mathrm{Geophysical}\\\\\\mathrm{Result}\\end{array}}\n",
    "$$\n",
    "\n",
    "Often there are numerous parameters that go into the blackbox, various combinations of which can greatly affect the output results! It can be extremely frustrating for users - especially when the documentation and usage examples cannot be found, and there is no open-source code to help figure it out.\n",
    "\n",
    "## Python and Scientific Libraries\n",
    "\n",
    "One of the major goals of this course is to strip away some of the **magic black box** nature of how students are commonly introduced to geophysical computing, and to provide a better conceptual idea of what is going on when software is applied to tackle geophysical problems.  That said, it is not our goal to \"reinvent the wheel\" on absolutely every algorithm because this would take far too long and one would never finish finding solutions to the problem at hand.  \n",
    "\n",
    "The approach we are taking in this course tries to find a middle way.  In particular, we will focus on the **Python** language and leverage its open-source, diverse and well validated computational toolkits including **Numpy** and **Scipy** as well as the strong plotting library **Matplotlib**.  \n",
    "\n",
    "However, to fully appreciate the theory behind some of these tools, it is important to code up some algorithms by oneself to help deepen the understanding of just what the algorithm is doing. Moreover, there will be scenarios where one has to develop a full computer program for a specialized task for which no community-based solution exists.  In these cases, Python may not be the optimal language in which to develop and thus there will be the need to include algorithms directly embedded into the code. Having a deeper understanding of the geophysical algoritms will be greatly helpful in these scenarios. Overall, it is our hope that this will assist students in designing better geophysical computing solutions that will facilitate achiement of the geophysical data processing goals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Course Aims\n",
    "\n",
    "The course is intended to provide senior-level undergraduate and first-year graduate students with material that aimed at improving the algorithmic, programming and computing skills, while simultaneously enhancing skills in handling geophysical data sets.  In some ways, these skills are those that are not really taught in any one particular geophysics course; however, they are often the \"glue\" that students need to help stick together disparate components learned in specific courses into a workable solution.\n",
    "\n",
    "If we were to broadly classify the material presented in this course, it could be broken down into the following three skill sets:\n",
    "\n",
    "* Handling Geophysical Data: \n",
    "    * **Data Interpolation**\n",
    "    * **Regression (Model Fitting)**\n",
    "    \n",
    "    \n",
    "* Applied Numerical Algorithms\n",
    "    * **Numerical Quadrature (Integration)**\n",
    "    * **Applied Linear Algebra for Geophysical Problems**\n",
    "    * **Solutions to Ordinary Differential Equations (ODEs)** \n",
    "    * **Differentiation and Discretization**\n",
    "    \n",
    "    \n",
    "* Numerical Solutions of Geophysical Partial Differential Equations (PDEs): \n",
    "    * **Elliptical PDEs (Laplace and Poisson Equations)**\n",
    "    * **Parabolic PDEs (2D Heat Flow)**\n",
    "    * **Hyperbolic PDEs (2D Acoustic wave equation)**\n",
    "    \n",
    "We briefly summarize each of these sections below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 2: Interpolation\n",
    "\n",
    "Geophysicists handle data that are acquired at regular and irregular (spatial and/or temporal intervals); however, most of the time the locations where data are acquired are insufficient in number or not in all of the desired locations required to make the corresponding geophysical or geological interpretation.   \n",
    "\n",
    "Let's look a the following example.  For simplicitly, let's say that we have 30+ elevation data points that were acquired over a 1 km$^2$ area.  An important question that will no doubt be asked is: *Based on this data, can we determine the elevation profile throughout the entire area?*\n",
    "\n",
    "<img src=\"figs/2D_Interpolation_intro.png\" width=\"500\">\n",
    "\n",
    "The good news is that yes we can create high-density maps using the acquired data points; however, the challenge is that there is no one unique map that we can create from the data that would be arguably correct.  To do this, we we will explore some of the fundamentals of **interpolation** and see a number of different ways in which plausible maps can be created in 1D and 2D.  Many of these algorithms even can be exported to higher dimensions for those cases where 3D, 4D and even 5D interpolation is required; however, these fall beyond the scope of this module.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 3: Numerical Quadrature / Integration\n",
    "\n",
    "One of the fundamental building blocks of modern calculus is the **integral**. We have all no doubt sat though numerous classes where the instructors went over analytic evaluation of 1D, 2D and perhaps 3D integrals\n",
    "\n",
    "$$\n",
    "J_1 = \\int_a^b f(x) \\, dx, \\quad \\quad J_2 = \\int_c^d\\int_a^b g(x,y)\\,dx \\,dy, \\quad \\quad J_3 = \\int_e^f \\int_c^d\\int_a^b h(x,y,z)\\, dx\\, dy \\,dz\n",
    "$$\n",
    "\n",
    "using many different tricks. However, for these scenarios a key constraint is that the functions to be evaluated, here $f(x)$, $g(x,y)$ and $h(x,y,z)$ were always specified *analytically*.  This represents a strong constraint on your ability to evaluate integrals!\n",
    "\n",
    "Here's a couple of interesting examples:\n",
    "\n",
    "<table><tr>\n",
    "<td> <img src=\"figs/Integration_1D_Intro.png\" alt=\"Drawing\" style=\"width: 400px;\"/> </td>\n",
    "<td> <img src=\"figs/comet_Interp_2D_Intro.png\" alt=\"Drawing\" style=\"width: 400px;\"/> </td>\n",
    "</tr></table>\n",
    "\n",
    "The left-hand side figure represents a frequency histogram of values of some recorded data set.  We see that some scientists have decided that this histogram can be represented by the model curve shown in red. Now, let's interpret these at probability density functions (PDFs) and say that we wanted to integrate this data and/or model curve over the interval $J=[0,a]$ to find out what the probabilty of something happening on the interval $J$. Arguably, it would easier to integrate the red curve if it were given by some analytic function.  However, how accurate would this actually be given that the model doesn't exact match data!\n",
    "\n",
    "The right-hand figure shows the an image of the Rosetta comet. Let's say that you were involved in a space mining operation that needed some key information regarding the comet.  If the operation is run entirely by solar power, an important question would be how much of the 2D surface could be facing the sun and could be covered by solar panels? Another question would be if the comet contains 1% per volume of iron ore, how much iron ore is there available in the comet?  Evidently, these questions would require evaluating 2D and 3D integrals over some irregular surface that is not easily described by analytic functions $g(x,y)$ and $h(x,y,z)$!  \n",
    "\n",
    "To address scenarios like this one, we will be exploring some of the **numerical quadrature** tools that are available in the python Scipy library."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 4: Applied Linear Algebra for Geophysical Problems\n",
    "\n",
    "In your academic career you have probably taken a semester-long course in Linear Algebra.  It is equally likely that this course was pretty abstract and didn't really make a strong connection with geophysics.  However, applied linear algebra comes up all of the time in geophysical computing - especially in the context of linearized system of equations. \n",
    "\n",
    "In many situations we will investigate this semester, the geophysical problem can be represeted by a straightforward matrix equation:\n",
    "\n",
    "$${\\bf Gm}= {\\bf d}$$\n",
    "\n",
    "where ${\\bf m}$ is some geophysical **model parameter** (e.g., acoustic wave speed, thermal conductivity, electric charge), ${\\bf d}$ is some sort of geophysical **data** (e.g., acoustic pressure, heat distibution, electric potential), and ${\\bf G}$ is a numerical representation of the **physics** and often **experimental geometry** that is used to forward model ${\\bf d}$ given ${\\bf m}$.  The following figure illustrates this concept using graphical depictions of ${\\bf G}$, ${\\bf m}$ and ${\\bf d}$.\n",
    "\n",
    "<img src=\"figs/Gm_d_Inversion.jpeg\" width=\"750\">\n",
    "\n",
    "A good portion of this course looks at how one can solve these types of forward modeling equations in an efficient and straightforward manner using different methods of **applied linear algebra**.  We will also examine scenarios where solutions cannot actually be generated due to various numerical challenges. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 5: Regression and Model Fitting\n",
    "\n",
    "Another important skill set to develop is fitting models to a given data set. It's likely that we've all done this with using the linear regression tool in Excel that allows us to find the parameters $m$ and $d$ such that the linear relation \n",
    "\n",
    "$$y=mx+d$$ \n",
    "\n",
    "is optimal in the least-squares sense.  You may have played around with higher order polynomials up to and including the same order of fit as the number of data points!\n",
    "\n",
    "While it may be numerically possible to perfectly fit an underlying model (except perhaps when interpolating!), this is seldom the optimal thing to do. One important reason for this is that data $f(x)$ almost *noisy* to some degree such that it is better represented by \n",
    "\n",
    "$$y=f(x+\\Delta x_{N})+\\Delta f$$\n",
    "\n",
    "where $\\Delta x_N$ is the noise in the position measurement and $\\Delta f_N$ is the noise in the measurement (e.g., instrument fluctuaions).  \n",
    "\n",
    "Thus, we need to develop tools and experience in fitting data sets like those in the figure below.\n",
    "\n",
    "<table><tr>\n",
    "<td> <img src=\"figs/Polynomial_Fit_9th_Order_Regression.png\" alt=\"Drawing\" style=\"width: 400px;\"/> </td>\n",
    "<td> <img src=\"figs/2D_REGRESSION.gif\" alt=\"Drawing\" style=\"width: 400px;\"/> </td>\n",
    "</tr></table>\n",
    "\n",
    "The left-hand image shows a number of data points (dots) along with three different curves that are fit to higher orders.  The right-hand image shows data acquired over an 2D area with significant scatter along with a 2D curve fit to it.  Of course we may be interested in fitting models to much higher order data in order to be able to generate some more complex model that we could possibly hope to visualize in 1D or 2D images. This type of analysis is proving to be fundamental in **data science** and **machine learning** applications in the geosciences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 6: Ordinary Differential Equations (ODEs)\n",
    "\n",
    "The next module in the course starts to move in a slightly different direction: we being examining numerical approaches to solving differential equations.  We first start with ordinary differential equations (ODEs), which again should be one of the math courses you have taken during your first two years.  Your \"diffy-e-q\" course probably largely focused on increasingly complex analytic solutions to increasingly more involved ODEs.  You probably memorized the formulae and the solutions approaches for your exams, and then promptly forgot these solution approaches after walking out of the exam hall.  (Don't worry, your instructors did this too ...)\n",
    "\n",
    "However, as you no doubt expect, very powerful numerical solvers for ODEs that can be used in place of going back to your Differential Equations textbook.  Importamtly, in many situations these routines can solve ODEs that cannot be solved by analytic methods!  Take for example the problem of a driven damped pendulum. The equation of motion for the angle $\\theta=\\theta(t,\\Omega)$ that the pendulum makes with the vertical is given by\n",
    "\n",
    "$$\\frac{d^2\\theta}{dt^2} = -\\frac{1}{Q} \\frac{d\\theta}{dt} + \\sin\\theta + d \\cos\\Omega t$$\n",
    "\n",
    "where $t$ is time, $Q$ is the quality factor (i.e., how slow or fast the oscillations die off), $d$ is the forcing amplitude, and $\\Omega$ is the driving frequency of the forcing. The ODE is **nonlinear** owing to the $\\sin\\theta$ term and there one needs numerical approaches to solve this system!\n",
    "\n",
    "The following figure shows the solution computed using Scipy ODE solvers in only a few lines of code.  The lower panel shows a complex driving frequency $\\Omega$, while upper panel shows the computed solution of the pendulum.\n",
    "\n",
    "<img src=\"figs/odePend.png\" width=\"500\">\n",
    "\n",
    "ODE solvers can be quite useful in geophysical and geoscience problems - especially in scenarios where you have a whole number of coupled ODEs!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 7: Numerical Differentiation / Discretization\n",
    "\n",
    "The next module continues our exploration of numerical methods to used to solve equations - especially partial differential equations (PDEs) such as the acoustic wave equation \n",
    "\n",
    "$$\\frac{\\partial^2 u }{\\partial x^2} = \\frac{1}{c^2}\\frac{\\partial^2 u}{\\partial t^2}$$\n",
    "\n",
    "To solve these equations numerically in a computer, we need to solve PDEs on a discrete solution grid that necessarily requires the process of **discretizing the physical system** including the partial differential operators (e.g., $\\frac{\\partial^2 }{\\partial x^2}$) comprising the PDE.  This is commonly done using **finite-difference** approximations, which commonly lead to numerical **stencils** and very **efficient** solution algorithms.\n",
    "\n",
    "In fact, for those of you who have taken digital signal processing, you hopefully will already be familiar with the concept of discretization. To remind you, the left panel below shows an image of the moon. The center panel shows this picture after applying a **convolutional Laplacian filter** for the purpose of emphasizing the discontinuous structure. The right panel shows the **2D finite-difference approximation** of the continuous Laplacian operators $\\nabla^2$.\n",
    "\n",
    "<table><tr>\n",
    "<td> <img src=\"figs/Moon_discretization.jpg\" alt=\"Drawing\" style=\"width: 500px;\"/> </td>\n",
    "<td> <img src=\"figs/Lapfilt.png\" alt=\"Drawing\" style=\"width: 200px;\"/> </td>\n",
    "</tr></table>\n",
    "\n",
    "However, we will go beyond this to explore how different discretizations affect the **stabilty** and **accuracy** of numerical solution of **time-dependent** PDEs.  In this module we will be taking a deeper look at the numerical approximations and discretization with the goal of developing the skills required to solve PDEs of interest like those in Modules 8, 9 and 10."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 8: Elliptical PDEs (Solutions to Laplace and Poisson Equations)\n",
    "\n",
    "Having tackled issues pertaining to discretization of partial differential operators and thereby PDEs, we are now ready to move on to actually solving some PDEs of interest in geophysics: **elliptical PDEs**.  One of the most straightforward 2D PDEs to solve is Poisson's equation for a potential surface $U(x,y)$:\n",
    "\n",
    "$$\\nabla^2 U = \\frac{\\partial^2 U}{\\partial x^2} + \\frac{\\partial^2 U}{\\partial y^2} = f(x,y)$$\n",
    "\n",
    "which reduced to Laplace's equation in the instance where $f(x,y)=0$ throughout the solution domain.  You may recall that one must specify the **boundary conditions** prior to calcuating solutions to these types of equations. \n",
    "\n",
    "The example below shows the potential surface solution related to the distribution of either four point masses (gravitational potential) or the response to four positive point charges (electrical potential).\n",
    "\n",
    "<img src=\"figs/solution_mixt.png\" width=\"600\">\n",
    "\n",
    "In this module we will be exploring solutions to such equations using some of the approaches we investigated in the applied linear algebra section (Module 4)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 9: Parabolic PDEs (Solutions to Heat Flow / Diffusion Equations)\n",
    "\n",
    "The next set of PDEs we will study fall into the family of **parabolic** equations.  In particular, we are looking to generate numerical solutions of 2D heat flow and 2D diffusion equations, both of which can be written in the following form:\n",
    "\n",
    "$${\\frac {\\partial \\phi (x,y,t)}{\\partial t}}=\\nabla \\cdot {\\big [}D(x,y)\\ \\nabla \\phi (x,y,t){\\big ]}$$\n",
    "\n",
    "where $\\phi (x,y,t)$ is, e.g., the distribution of heat through the 2D solution domain, and $D(x,y)$ is the heterogeneous thermal conductivity field.  Computing solutions requires setting both the **initial condition** (i.e., at time $t=0s$) and the **boundary conditions** at the edges of the computational domain.\n",
    "\n",
    "The left panels shows example below shows the time evolution of the heat distribution $U(x,y,t)$ for a square computational domain where the four edges are held at $0^\\circ$C and all interior points initially start at $100^\\circ$C. The right panel shows a cross-section through the solution and is easier to see how the $U(x,y=0.5,t)$ evolves over the computed solution time.\n",
    "\n",
    "<img src=\"figs/2dheat.gif\" width=\"800\">\n",
    "\n",
    "In this module we will be developing numerical solutions to the 2D heat flow equation, validating them against analytical solutions to a PDE system, and then using your validated code to solve some interesting geophysical problems for which **no analytical solutions exist**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module 10: Hyperbolic PDEs (Solutions to Acoustic Wave Equation)\n",
    "\n",
    "The final set of PDEs we will study fall into the family of **hyperbolic** equations.  In particular, we are looking to generate numerical solutions of 2D acoustic wave equation, which is written as the following:\n",
    "\n",
    "$$\\left[ \\frac{ \\partial^2 }{\\partial x^2}+\\frac{ \\partial^2 }{\\partial x^2} -\\frac{1}{c^2(x,y)} {\\frac {\\partial^2}{\\partial t^2}}\\right] U(x,y,t) = f(x,y,t)$$\n",
    "\n",
    "where $U(x,y,t)$ is the displacement of the acoustic wave disturbance, and $c(x,y)$ is the heterogeneous acoustic wave speed, and $f(x,y,t)$ is the force source distribution (e.g., what is causing the acoustic wave disturbance).  Computing solutions requires setting the initial conditions and the boundary conditions at the edges of the computational domain.\n",
    "\n",
    "The example below shows a numerical solution to the 2D acoustic wave equation that is modeling the temporal and spatial evolution of a taut square drum head (i.e., clamped boundaries with zero displacement) and some initial amplitude and/or velocity distribution.\n",
    "\n",
    "<img src=\"figs/2dwave.gif\" width=\"400\">\n",
    "\n",
    "In this module we will be developing numerical solutions to the 2D acoustic, validating them against analytical solutions to a PDE system like in the example above, and then using our validated code to solve some interesting geophysical wave propagation problems for which there are no analytical solution.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}