Ex12-Linear-Advection-new.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 12: Linear Advection\n",
    "\n",
    "After implementing different schemes for diffusion in the last exercise, we are now only a small step away from solving the entire transport equation. Recall that the entire equation is:\n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t} + \\nabla \\cdot \\left(\\vec{v} u\\right) - \\nabla \\cdot \\left( \\kappa \\nabla u \\right) + S = 0$$\n",
    "\n",
    "Today, we are going to implement steps to solve the advective term in a transient simulation. To get started, we will now only consider the 1-D case of (incompressible) linear advection as:\n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t} + v_x \\frac{\\partial u}{\\partial x} = 0$$\n",
    "\n",
    "The main difference to the previous examples of diffusion is that, now, we only have first derivatives to solve. Let's recall the standard finite difference schemes to discretise first-order derivatives at a location $i$:\n",
    "\n",
    "Forward difference:\n",
    "\n",
    "$$\\left(\\frac{\\partial u}{\\partial x}\\right)_i = \\frac{u_{i+1} - u_i}{\\Delta x} + \\mathcal{O}(\\Delta x)$$\n",
    "\n",
    "Backward difference:\n",
    "\n",
    "$$\\left(\\frac{\\partial u}{\\partial x}\\right)_i = \\frac{u_{i} - u_{i-1}}{\\Delta x} + \\mathcal{O}(\\Delta x)$$\n",
    "\n",
    "\n",
    "Central difference:\n",
    "\n",
    "$$\\left(\\frac{\\partial u}{\\partial x}\\right)_i = \\frac{u_{i+1} - u_{i-1}}{2 \\Delta x} + \\mathcal{O}(\\Delta x^2)$$\n",
    "\n",
    "As the central difference scheme is more accurate than the one-sided difference schemes, it is usually the first choice in FD discretisations of first derivatives. So, let's implement it here!\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#from IPython.core.display import HTML\n",
    "#css_file = 'nre2_style.css'\n",
    "#HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "plt.rcParams['figure.figsize'] = [10., 7.]\n",
    "plt.rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the case of diffusion, we now also have the choice to define the point in our time interval where we estimate the spatial gradient. We will first start with an explicit implementation, as this is the straight-forward way to solve the problem.\n",
    "\n",
    "We are implementing the standard \"Forward-Time-Central-Space\" (FTCS) scheme:\n",
    "\n",
    "$$\\frac{u_i^{k+1} - u_i^k}{\\Delta t} = v \\frac{u_{i+1} - u_{i-1}}{2 \\Delta x}$$\n",
    "\n",
    "If we sort this equation to have unknown terms (at time $k+1$) on the left, and known terms (at time $k$) on the right as before, we obtain a nice explicit equation:\n",
    "\n",
    "$$u_i^{k+1} = u_i^k - \\frac{1}{2} \\frac{v \\Delta t}{\\Delta x} (u_{i+1}^k - u_{i-1}^k)$$ \n",
    "\n",
    "Or, using the definition for the Courant number $Cr = \\frac{v \\Delta t}{\\Delta x}$:\n",
    "\n",
    "$$u_i^{k+1} = u_i^k - \\frac{1}{2} Cr (u_{i+1}^k - u_{i-1}^k)$$ \n",
    "\n",
    "With all your knowledge of Python, it will be very easy for you to implement this equation in a function to solve a linear advection problem.\n",
    "\n",
    "As initial condition, we will use gaussian peak as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define model dimension\n",
    "L = 100\n",
    "nx = 100\n",
    "dx = L/nx\n",
    "x = np.arange(dx/2., L+dx/2., dx)\n",
    "u_init = 1 * np.exp(-((np.linspace(0,1,nx)-0.3) / 0.1)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.1, 1.1)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VfWd9/H3N4QQRC6igKuGOwqC3AIEAgpHwCFRLFqq\n1dE6tS12zapO5+laj9qn48h0TVfb6UWlOig6dQamigUvoFIRlCNyCYS7SLjIRZA7BEGC3JLf80dO\nMKaEBDhn/87Z5/Na6yzOPtnZ+5O9wif7/M6+mHMOEREJrwzfAUREJLFU9CIiIaeiFxEJORW9iEjI\nqehFREJORS8iEnKZQa7MzHQsp4jIBXDO2YV+b+B79M45PZzj8ccf954hWR7aFtoW2hbnflwsDd2I\niIScil5EJORU9J5EIhHfEZKGtsVXtC2+om0RPxaP8Z96r8zMBbk+EZEwMDNcKn0YKyIiwVLRi4iE\nnIpeRCTkVPQiIiGnohcRCTkVvYhIyKnoRURCTkUvIhJydRa9mf2Xme01szXnmGeCmW0ys1Vm1ie+\nEUVE5GLUZ4/+RWBUbV80s0Kgs3PuauBHwLNxyiYiInFQZ9E75xYAh84xyxhgcmzeJUBzM2sTn3gi\nInKx4jFGfxWwo9r0zthrkgJOnz7N0aNHfccQkQQK9A5TAOPHjz/zPBKJ6Ap1HmzatIkpU6awaNEi\nli5dSmZmJtOmTWPEiBG+o4kIEI1GiUajcVteva5eaWbtgTedc73O8rVngXnOuVdi0+uBYc65vWeZ\nV1ev9OzQoUPk5uZy++23M3LkSAYNGsTatWu54447eOqpp7jrrrt8RxSRGi726pX13aO32ONsZgI/\nBl4xs0HA52crefHPOcf999/PmDFj+MMf/nDm9aFDhzJ37lxuvvlm9uzZwz//8z97TCki8VZn0ZvZ\nS0AEuNzMtgOPA1mAc85Ncs7NMrObzewToAy4P5GB5cJNmDCBnTt38pe//OVvvtazZ08WLlxIJBKh\na9euFBYWekgoIomgG4+kieLiYm655RaKioro1KlTrfNNmzaN//iP/2Dp0qWYXfA7RRGJI914ROpU\nXl7OPffcw8SJE89Z8gBjx47lxIkTvPXWWwGlE5FEU9GngbfffpuWLVsyduzYOufNyMjg3/7t3/jX\nf/1X9O5LJBxU9Gng6aef5sEHH6z3/LfddhtmxhtvvJHAVCISFI3Rh9z69esZNmwY27dvp1GjRvX+\nvjfffJOf//znrFq1iowM7Q+I+KQxejmnZ555hnHjxp1XyQOMHj2aRo0a8dprryUomYgERXv0IXbk\nyBE6dOjA6tWradu27Xl//9SpU/nv//5v3nnnnQSkE5H60h691GrKlCkMHz78gkoe4NZbb2Xx4sUc\nOHAgzslEJEgq+pByzp33h7A1NWnShMLCQl599dU4JhORoKnoQ2rBggVkZGQwbNiwi1rOd77zHaZO\nnRqnVCLig4o+pN544w3uvPPOiz67tbCwkFWrVrF79+44JRORoKnoQ8g5x4wZM/jmN7950cvKzs7m\n1ltvZdq0aXFIJiI+qOhDqKSkhFOnTtGnT3xu33vXXXfxyiuvxGVZIhI8FX0IVe3Nx+uiZCNHjmTD\nhg1s3749LssTkWCp6ENo5syZcRm2qZKVlcXtt99+1ssbi0jyU9GHzJ49e85c9iCeNHwjkrpU9CHz\n1ltvMWrUKLKysuK63KFDh7Jx40b2798f1+WKSOKp6ENmxowZjBkzJu7LbdiwIcOGDeO9996L+7JF\nJLFU9CFSVlbGBx98kLDbAN50003MmTMnIcsWkcRR0YfInDlzyMvLo0WLFglZflXR68J0IqlFRR8i\nb7/9NrfeemvClt+1a1ecc2zcuDFh6xCR+FPRh8i8efMYMWJEwpZvZhq+EUlBKvqQ2LFjB0eOHKFH\njx4JXY+KXiT1qOhDYt68eUQikbidDVubESNGEI1GOXXqVELXIyLxo6IPiWg0SiQSSfh6WrduTefO\nnVmyZEnC1yUi8aGiD4l58+Zx4403BrIuDd+IpBYVfQhs27aNL7/8km7dugWyPhW9SGpR0YdA1bBN\nosfnq1x//fV89NFHHD58OJD1icjFUdGHQFDj81Wys7MZOHAgH374YWDrFJELV6+iN7MCM1tvZhvN\n7JGzfL2Zmc00s1Vm9pGZfS/uSeWsnHOBjs9Xuf7661m4cGGg6xSRC1Nn0ZtZBvA0MAroAdxtZjUH\ng38MfOyc6wPcCPzezDLjHVb+1tatWzl58iTXXHNNoOsdMmQICxYsCHSdInJh6rNHnwdscs596pw7\nBUwFal4e0QFNY8+bAgedc6fjF1NqE41GufHGGwMbn68yaNAgVq5cyYkTJwJdr4icv/oU/VXAjmrT\nn8Veq+5poLuZ7QJWAz+JTzypS9WJUkFr2rQpXbt2Zfny5YGvW0TOT7yGV0YBK51zw82sMzDHzHo5\n547WnHH8+PFnnkciES8lFSbRaJTHHnvMy7qHDBnCwoULGTx4sJf1i4RVNBolGo3GbXlW1yVnzWwQ\nMN45VxCbfhRwzrnfVJvnLeBXzrmFsen3gEecc8tqLMvpErfxs2PHDvr378+ePXsCH7oB+Mtf/sKf\n//xnZsyYEfi6RdKJmeGcu+D/5PUZuikGuphZezPLAu4CZtaY51NgZCxQG+AaYMuFhpL6WbRoEfn5\n+V5KHr7ao9cfb5HkVmfRO+fKgQeBd4GPganOuRIz+5GZPRCb7d+BwWa2BpgDPOycK01UaKm0ePFi\n8vPzva3/qquuomnTpmzYsMFbBhGpW73G6J1z7wBda7z2XLXnu6kcp5cALV68mN/+9rdeM1x//fUs\nWLAgsMsviMj505mxKerLL79k7dq19O/f32sOHU8vkvxU9Clq+fLldO/enUsuucRrDp0hK5L8VPQp\nyvf4fJXu3btz8OBB9uzZ4zuKiNRCRZ+iFi1alBTHr2dkZJCfn6+9epEkpqJPQc65pNmjh68OsxSR\n5KSiT0Hbtm0jIyODdu3a+Y4CVF73RrcWFEleKvoUtHjxYgYPHuztRKma+vfvz6pVq3TDcJEkpaJP\nQVVnxCaLZs2a0bFjR9asWeM7ioichYo+BSXT+HyVgQMHavhGJEmp6FNMWVkZ69evJzc313eUr1HR\niyQvFX2KWbZsGb169SI7O9t3lK9R0YskLxV9ilmyZAkDBw70HeNv9OjRg507d3Lo0CHfUUSkBhV9\nilm6dCl5eXm+Y/yNzMxMcnNzKS4u9h1FRGpQ0aeYZC160PCNSLJS0aeQ3bt3c/ToUTp37uw7ylmp\n6EWSk4o+hRQXF5OXl5c0J0rVVFX0uuOUSHJR0aeQZB62AcjJySErK4utW7f6jiIi1ajoU0iyFz1o\n+EYkGanoU0RFRQXFxcUMGDDAd5RzUtGLJB8VfYr45JNPaN68OW3atPEd5ZxU9CLJR0WfIlJh2AYq\nr2S5Zs0aTp486TuKiMSo6FNEqhT9pZdeSseOHVm7dq3vKCISo6JPEalS9AB5eXk6Q1YkiajoU8DJ\nkyf56KOPku6KlbUZMGAAS5cu9R1DRGJU9ClgzZo1dO7cmUsvvdR3lHoZMGCA9uhFkoiKPgWk0rAN\nQK9evdi8eTNlZWW+o4gIKvqUsHTp0qQ/fr66rKwsrrvuOlasWOE7ioigok8JqbZHDxq+EUkm9Sp6\nMysws/VmttHMHqllnoiZrTSztWY2L74x09eRI0f49NNPue6663xHOS868kYkedRZ9GaWATwNjAJ6\nAHebWbca8zQHngFGO+euA+5IQNa0tHz5cnr37k3Dhg19Rzkv2qMXSR712aPPAzY55z51zp0CpgJj\naszz98CrzrmdAM65A/GNmb5S4fo2Z9O1a1f279/PwYMHfUcRSXv1KfqrgB3Vpj+LvVbdNUBLM5tn\nZsVm9t14BUx3VdegTzUZGRn069ePZcuW+Y4ikvbi9WFsJpALFAIFwGNm1iVOy05rqXbETXUavhFJ\nDpn1mGcn0K7adE7steo+Aw44544Dx81sPtAb+KTmwsaPH3/meSQSIRKJnF/iNLJv3z4OHz5Mly6p\n+TdzwIABTJ482XcMkZQTjUaJRqNxW57Vdds3M2sAbABGALuBpcDdzrmSavN0A/5I5d58I2AJ8B3n\n3Loay3K6zVz9vf322zz55JPMmTPHd5QLsn37dgYOHMiuXbuS9vaHIqnAzHDOXfB/ojqHbpxz5cCD\nwLvAx8BU51yJmf3IzB6IzbMemA2sAYqASTVLXs5fKg/bALRt25aKigp27qz5BlBEglSfoRucc+8A\nXWu89lyN6d8Bv4tfNCkuLmbcuHG+Y1wwMyMvL48lS5aQk5PjO45I2tKZsUnKOZeyR9xUpw9kRfxT\n0Sepbdu20bBhQ666quaRrKlFZ8iK+KeiT1KpeqJUTQMGDGDZsmVUVFT4jiKStlT0SSosRX/55ZfT\nqlUrNmzY4DuKSNpS0SepMIzPV8nLy9Mdp0Q8UtEnofLyclasWEH//v19R4kLFb2IXyr6JFRSUkKb\nNm1o2bKl7yhxoaIX8UtFn4SWLFnCwIEDfceIm759+7Ju3TqOHz/uO4pIWlLRJ6GwFX3jxo3p2rUr\nq1ev9h1FJC2p6JNQ2IoeNHwj4pOKPsmUlZWxadMmevfu7TtKXKnoRfxR0SeZ5cuX07NnTxo1auQ7\nSlyp6EX8UdEnmTAO2wBce+217N69m0OHDvmOIpJ2VPRJZunSpaEs+gYNGtC3b1/dWlDEAxV9klmy\nZElozoitqeqSxSISLBV9Etm9ezdlZWUpe+vAugwcOFDj9CIeqOiTSNXefFhvuzdo0CAWL16Mbicp\nEiwVfRIJ6/h8lZycHLKzs9m8ebPvKCJpRUWfRMI8Pl8lPz+foqIi3zFE0oqKPkmUl5ezbNmy0Bd9\n1fCNiARHRZ8kNmzYQKtWrbjiiit8R0mo/Px8Fb1IwFT0SaKoqCjU4/NV+vbty4YNGygrK/MdRSRt\nqOiTxKJFixg8eLDvGAmXnZ1Nz549deKUSIBU9Eli4cKFDBkyxHeMQAwaNEgfyIoESEWfBA4ePMjO\nnTu57rrrfEcJhMbpRYKlok8CRUVF5OXlkZmZ6TtKIKqKXidOiQRDRZ8E0mnYBqBt27Y0aNCAbdu2\n+Y4ikhZU9EkgXT6IrWJmGr4RCVC9it7MCsxsvZltNLNHzjHfADM7ZWbfil/EcDt16hTLli1j0KBB\nvqMESidOiQSnzqI3swzgaWAU0AO428y61TLfr4HZ8Q4ZZqtXr6ZTp040b97cd5RA6VIIIsGpzx59\nHrDJOfepc+4UMBUYc5b5HgKmA/vimC/0Fi5cmFbDNlX69evHunXrdOKUSADqU/RXATuqTX8We+0M\nM/sGcJtzbiIQzmvsJki6jc9Xady4MX369NFevUgA4vVh7JNA9bF7lX09LVq0KK2OuKnuhhtu4MMP\nP/QdQyT06nPg9k6gXbXpnNhr1fUHplrlHTOuAArN7JRzbmbNhY0fP/7M80gkQiQSOc/I4bF9+3ZO\nnjxJp06dfEfxYujQofzud7/zHUMk6USjUaLRaNyWZ3WdtGJmDYANwAhgN7AUuNs5V1LL/C8Cbzrn\nXjvL15xOkvnK1KlTeeWVV3j99dd9R/Hi888/p23bthw8eJCsrCzfcUSSlpnhnLvgkZI6h26cc+XA\ng8C7wMfAVOdciZn9yMweONu3XGiYdJOu4/NVWrRoQZcuXVi+fLnvKCKhVq9z7p1z7wBda7z2XC3z\nfj8OudLC/PnzmThxou8YXlWN0+fn5/uOIhJaOjPWk9LSUrZs2UL//v19R/Fq6NChzJ8/33cMkVBT\n0Xsyf/58Bg8eTMOGDX1H8eqGG25g4cKFVFRU+I4iEloqek+i0WhaH3FUpU2bNrRq1Yq1a9f6jiIS\nWip6T1T0X9HwjUhiqeg9KC0tZevWrfTr1893lKSgE6dEEktF74HG57+uao9e51iIJIaK3gMN23xd\nhw4daNCgAZs3b/YdRSSUVPQeqOi/zsyIRCLMmzfPdxSRUFLRB6xqfD43N9d3lKQycuRI5s6d6zuG\nSCip6AOm8fmzGzlyJO+9956OpxdJABV9wDRsc3Y5OTm0atWKlStX+o4iEjoq+oCp6Gt30003MWfO\nHN8xREJHRR+g/fv3a3z+HFT0Iomhog/Q7NmzGT58uMbnaxGJRFiyZAnHjh3zHUUkVFT0AXrnnXco\nLCz0HSNpNW3alD59+rBgwQLfUURCRUUfkIqKCmbPnk1BQYHvKElNwzci8aeiD8jy5ctp3bo17dq1\nq3vmNKaiF4k/FX1A/vrXv2pvvh7y8vLYtm0be/fu9R1FJDRU9AHR+Hz9ZGZmMmzYMN5//33fUURC\nQ0UfgNLSUtauXcsNN9zgO0pKuOmmm5g9e7bvGCKhoaIPwJw5cxg2bBiNGjXyHSUljB49mlmzZlFe\nXu47ikgoqOgDoPH589OhQwe+8Y1vsHjxYt9RREJBRZ9gFRUVGp+/AGPGjOGNN97wHUMkFFT0CbZ6\n9WqaNWtGp06dfEdJKWPGjGHGjBm665RIHKjoE2zmzJnccsstvmOknL59+3LixAlKSkp8RxFJeSr6\nBJs2bRrf/va3fcdIOWam4RuROFHRJ1BJSQmHDh0iPz/fd5SUVDV8IyIXR0WfQNOnT2fs2LFkZGgz\nX4hhw4axadMmdu3a5TuKSEpTAyXQ9OnTueOOO3zHSFkNGzakoKCAN99803cUkZRWr6I3swIzW29m\nG83skbN8/e/NbHXsscDMesY/amrZuHEj+/btY/Dgwb6jpLTbbrtN4/QiF6nOojezDOBpYBTQA7jb\nzLrVmG0LMNQ51xv4d+D5eAdNNVXDNg0aNPAdJaUVFBSwcOFCDh065DuKSMqqzx59HrDJOfepc+4U\nMBUYU30G51yRc+5wbLIIuCq+MVOPjraJj2bNmjFq1CimTZvmO4pIyqpP0V8F7Kg2/RnnLvIfAn+9\nmFCp7pNPPmHXrl26iFmc3HfffUyePNl3DJGUlRnPhZnZjcD9wPW1zTN+/PgzzyORCJFIJJ4RksL0\n6dP51re+pWGbOCkoKOAHP/gBmzdvpnPnzr7jiCRcNBolGo3GbXlW1ynmZjYIGO+cK4hNPwo459xv\naszXC3gVKHDOba5lWS7sp7Q75+jduzcTJkwI5R8xX37yk59w2WWXfW1HQSRdmBnOObvQ76/P0E0x\n0MXM2ptZFnAXMLNGiHZUlvx3ayv5dLF06VKOHTvG0KFDfUcJlfvuu48pU6bo2jciF6DOonfOlQMP\nAu8CHwNTnXMlZvYjM3sgNttjQEvgP81spZktTVjiJPf888/zwx/+UCdJxVlubi7Z2dksWrTIdxSR\nlFPn0E1cVxbyoZsjR47Qvn17SkpKuPLKK33HCZ3f/OY3bNmyheeee853FJFABTF0I/X08ssvM3z4\ncJV8gtxzzz1Mnz6d48eP+44iklJU9HH0/PPPM27cON8xQisnJ4e+ffvy+uuv+44iklJU9HGycuVK\n9u/fz0033eQ7Sqg99NBDPPnkk/pQVuQ8qOjjpOpDWB07n1ijR4+mtLRUH8qKnAd9GBsHX3zxBe3b\nt2fNmjXk5OT4jhN6zzzzDO+//z6vvvqq7ygigbjYD2NV9HHw29/+lhUrVvDyyy/7jpIWysrKaN++\nPUuWLNGZspIWVPSeHTt2jM6dO/Puu+/Ss2faX505MD/72c84duwYTz31lO8oIgmnovdswoQJRKNR\nXnvtNd9R0srOnTvp2bMnW7ZsoUWLFr7jiCSUit6jEydO0LlzZ2bMmEG/fv18x0k73/3ud+nZsycP\nP/yw7ygiCaWi92jSpEm88cYbzJo1y3eUtPTRRx8xcuRINm7cSPPmzX3HEUkYFb0np06d4pprruGl\nl14iPz/fd5y09f3vf582bdrwq1/9yncUkYRR0Xvywgsv8NJLL/H+++/7jpLWdu7cSa9evVixYgXt\n27f3HUckIVT0Hhw8eJAePXowa9YscnNzfcdJe48//jibN2/mf//3f31HEUkIFb0HDzzwANnZ2UyY\nMMF3FAGOHj3KNddcw4wZMxgwYIDvOCJxp6IP2OLFi/n2t7/NunXr9AFgEnnhhReYPHkyH3zwAWYX\n/P9BJCnpMsUBOn36NP/4j//I73//e5V8krn//vs5duwYzz77rO8oIklHe/Tn4YknnuDtt99mzpw5\n2mtMQhs2bGDIkCF8+OGHXHvttb7jiMSNhm4Csnz5cgoKCli0aBFXX3217zhSi0mTJjFx4kSKiopo\n1KiR7zgicaGhmwDs37+fsWPH8uyzz6rkk9y4cePo0KED//Iv/+I7ikjS0B59HU6fPk1BQQH9+/fn\n17/+te84Ug8HDhygT58+TJo0iZtvvtl3HJGLpj36BPv5z3+OmfHLX/7SdxSppyuuuIJp06bxve99\nTzcoEUFFf05/+MMfmD59Oi+//LLuHJVi8vPzmTJlCrfffjurV6/2HUfEKxX9WTjn+MUvfsGzzz5L\nNBrliiuu8B1JLsCoUaN4+umnKSwsZNOmTb7jiHiT6TtAsnHO8eijjzJr1izmz5/PlVde6TuSXIQ7\n7riDI0eOEIlEmDZtGoMHD/YdSSRw2qOv5vDhw9x777289957RKNRlXxI/OAHP+D555/n9ttv57nn\nniPVDggQuVgq+pgPP/yQPn360KxZMz744AMuv/xy35Ekjm6++WYWLFjAhAkTGDduHF988YXvSCKB\nSfuiP3ToEA8//DB33nknf/zjH5k4cSJNmjTxHUsS4Oqrr6aoqIjy8nK6devG5MmTqaio8B1LJOHS\ntugPHz7ML37xC66++moOHDjAqlWrGD16tO9YkmBNmzblxRdf5LXXXuOZZ55h8ODBzJ07V8M5Emr1\nKnozKzCz9Wa20cweqWWeCWa2ycxWmVmf+MaMD+ccRUVFPPTQQ3Tp0oXNmzdTVFTEn/70J9q0aeM7\nngRo4MCBLF68mB//+Mf89Kc/pVu3bjzxxBOUlpb6jiYSd3WeGWtmGcBGYASwCygG7nLOra82TyHw\noHPuFjMbCDzlnBt0lmUFfmbs559/zoIFC/jggw94/fXXadCgAffeey/33HMPnTp1CjSLJCfnHAsX\nLmTixInMnDmT/v37U1hYSGFhIT169CAjI23f+EqSSPhFzcxsEPC4c64wNv0o4Jxzv6k2z7PAPOfc\nK7HpEiDinNtbY1kJKXrnHAcPHuSzzz5j27ZtrFu3jnXr1rF27Vo2b95MXl4ekUiEwsJC+vXrpytP\nSq3KysqIRqPMmjWL2bNns2/fPnr37k1ubi7du3enY8eOdOzYkXbt2umiaRKYIIp+LDDKOfdAbPpe\nIM8590/V5nkT+JVzblFsei7wsHNuRY1lueLiYrZt20ZFRQXOOSoqKqioqKC8vJzy8nJOnz595nHy\n5ElOnDjBiRMnOH78OGVlZZSVlXH06FEOHTpEaWkppaWl7N27l8aNG5OTk0O7du3o3r37mUefPn3I\nysq60O0jaa60tJRVq1axfPlyNmzYwNatW9m6dSufffYZl1xyCa1bt6ZVq1Y0b96cZs2a0axZM5o0\naULjxo3Jzs4mOzubhg0bnnlkZmbSoEGDM4+MjAwyMjIws795AH/z77meS/g0bdqUSCRy0UUf+AlT\njz32GJ9++ikArVu35sorrzzzC9+gQQMyMzPPPLKysmjUqBGNGjXisssuIycnhyZNmnDppZdy2WWX\n0bJlS1q2bEnr1q11pIwkRMuWLRk+fDjDhw//2uvOOT7//HP27dvHvn37OHz4MF988QWHDx/m2LFj\nHD9+nC+//JIvvviCU6dOnXlU7cxU7dhU39lxzp15VK2j+r/nei7hcuDAAQ4ePEiLFi3icmG++g7d\njHfOFcSm6zN0sx4YFtTQjYhImAVx9cpioIuZtTezLOAuYGaNeWYC98UCDQI+r1nyIiLiR51DN865\ncjN7EHiXyj8M/+WcKzGzH1V+2U1yzs0ys5vN7BOgDLg/sbFFRKS+dOMREZEkpxuPiIjIOanoRURC\nTkUvIhJyKnoRkZBT0YuIhJyKXkQk5FT0IiIhp6IXEQk5Fb2ISMip6EVEQk5F70k0GvUdIWloW3xF\n2+Ir2hbxo6L3RL/EX9G2+Iq2xVe0LeJHRS8iEnIqehGRkAv8MsWBrUxEJEQSenNwERFJbRq6EREJ\nORW9iEjIBVb0ZlZgZuvNbKOZPRLUepOBmeWY2ftm9rGZfWRm/xR7/TIze9fMNpjZbDNr7jtrEMws\nw8xWmNnM2HS6bofmZjbNzEpivxsD03hb/B8zW2tma8zsz2aWlU7bwsz+y8z2mtmaaq/V+vOb2c/M\nbFPsd+fv6lp+IEVvZhnA08AooAdwt5l1C2LdSeI08FPnXA8gH/hx7Od/FJjrnOsKvA/8zGPGIP0E\nWFdtOl23w1PALOfctUBvYD1puC3M7BvAQ0Cuc64XkAncTXptixep7Mfqzvrzm1l34E7gWqAQ+E8z\nO+cHtUHt0ecBm5xznzrnTgFTgTEBrds759we59yq2POjQAmQQ+U2+J/YbP8D3OYnYXDMLAe4GXih\n2svpuB2aATc4514EcM6dds4dJg23RUwDoImZZQKNgZ2k0bZwzi0ADtV4ubaf/5vA1NjvzDZgE5Ud\nW6ugiv4qYEe16c9ir6UdM+sA9AGKgDbOub1Q+ccAaO0vWWCeAP4vUP1wr3TcDh2BA2b2YmwYa5KZ\nXUIabgvn3C7g98B2Kgv+sHNuLmm4LWpoXcvPX7NPd1JHn+rD2ACZ2aXAdOAnsT37mse2hvpYVzO7\nBdgbe3dzrreaod4OMZlALvCMcy4XKKPyrXpa/U4AmFkLKvde2wPfoHLP/h7ScFvU4YJ//qCKfifQ\nrtp0Tuy1tBF7SzodmOKcmxF7ea+ZtYl9/Upgn698ARkCfNPMtgAvA8PNbAqwJ822A1S+q93hnFsW\nm36VyuJPt98JgJHAFudcqXOuHHgdGEx6bovqavv5dwJtq81XZ58GVfTFQBcza29mWcBdwMyA1p0s\n/gSsc849Ve21mcD3Ys//AZhR85vCxDn3/5xz7Zxznaj8HXjfOfdd4E3SaDsAxN6S7zCza2IvjQA+\nJs1+J2K2VWmFAAAAx0lEQVS2A4PMLDv2oeIIKj+sT7dtYXz9nW5tP/9M4K7YkUkdgS7A0nMu2TkX\nyAMoADZQ+cHBo0GtNxkeVO7JlgOrgJXAitj2aAnMjW2Xd4EWvrMGuE2GATNjz9NyO1B5pE1x7Pfi\nNaB5Gm+Lx6k8SGENlR88NkynbQG8BOwCTlD5h+9+4LLafn4qj8D5JLbN/q6u5esSCCIiIacPY0VE\nQk5FLyIScip6EZGQU9GLiIScil5EJORU9CIiIaeiFxEJORW9iEjI/X9lGn8iPpQtigAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xad95978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(x, u_init, 'k-', lw=1)\n",
    "ax.set_ylim([-0.1,1.1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would like to know how this \"peak\" of $u$ is advancing over time in a fluid field with velocity $v = 2.0$. Let's say that we would like to determine the total movement for a time $T = 20$. \n",
    "\n",
    "In fact, you surely realise that you don't need a numerical method to solve this problem and that we can directly obtain the analytical solution.\n",
    "\n",
    "**Your task**: calculate the advected peak and create a plot with initial and advected peak:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# your code here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical solution\n",
    "\n",
    "So, let's now implement the FTCS solution for this problem. We now need to decide on an additional aspect: the number of time steps (or, respectively, the time step size).\n",
    "\n",
    "As a start, let's take 40 time steps: $n_t = 40$.\n",
    "\n",
    "**Your task**: Write a function to solve this problem numerically\n",
    "\n",
    "Suggestion: write a separate function to solve one time step with the FTCS method so that you can exchange it later (you might remember why from the lecture):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# your code here\n",
    "def ftcs(u, v, dt, dx):\n",
    "    \"\"\"Advance one time step with FTCS scheme\"\"\"\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-9-4a96c140d11f>, line 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-9-4a96c140d11f>\"\u001b[1;36m, line \u001b[1;32m3\u001b[0m\n\u001b[1;33m    nt =\u001b[0m\n\u001b[1;37m         ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# define additional parameter and perform iteration over time:\n",
    "\n",
    "nt = \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a plot of the result and compare it to the analytical solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I hope that the result looks different to what you expected from the first consideration of the problem, above!\n",
    "\n",
    "Of course, after following the lecture you will know what happened here: the FTCS scheme is, unfortunately, *unconditionally unstable*! So, whatever we do, we will not get a reasonable result!\n",
    "\n",
    "To solve this problem with a first-order finite difference scheme, we therefore have to revert to the one-sided differences and live with the fact that we only get first-order accuracy.\n",
    "\n",
    "When considering these schemes, we have two choices now: forward difference and backward difference. Which one is the logical one to use in this case? How do you make this choice - and why?\n",
    "\n",
    "**Your task**: write a function that always chooses the right one-sided difference scheme and then performs one iteration in time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# your code here\n",
    "def upwind(u, v, dt, dx):\n",
    "    \"\"\"Advance one time step with upwind scheme\"\"\"\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Perform simulations with the upwind scheme\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a plot of the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implicit solutions\n",
    "\n",
    "We talked about the Courant criterion and the famous Courant-Friedrichs-Levy (CFL) condition in the lecture. We encountered a similar condition before: in the case of the diffusion equation, we had to consider a criterion, as well (**Question: what was it?**). \n",
    "\n",
    "We could resolve a similar issue for the diffusion equation by using an implicit approach. Let's try the same thing here.\n",
    "\n",
    "We will implement a Crank-Nicolson scheme with equal weighting between beginning ($k$) and end ($k+1$) of the time step to determine the spatial gradient. Writing down the equation and putting $k+1$ terms on the left, and $k$ terms on the right, as before, we obtain:\n",
    "\n",
    "$$-\\frac{1}{4} Cr u_{i-1}^{k+1} + u_i^{k+1} + \\frac{1}{4} Cr u_{i+1}^{k+1} = \n",
    "\\frac{1}{4} Cr u_{i-1}^k + u_i^k - \\frac{1}{4} Cr u_{i+1}^k$$\n",
    "\n",
    "Where $Cr$ is the Courant number. This scheme is, in fact, unconditionally stable and second order accurate.\n",
    "\n",
    "**Your task**: use a C-N method to solve the advection problem.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further ideas\n",
    "\n",
    "We also talked about other ways to increase accuracy in explicit schemes in the lecture, specifically:\n",
    "\n",
    "- the \"Leapfrog\" scheme\n",
    "- the \"Lax-Wendroff\" scheme\n",
    "- Use a C-N scheme, but with a higher-order approximation in space\n",
    "\n",
    "Implement one of those schemes and compare the accuracy to the simple upwind scheme above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}