{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One version of the [advection equation](https://en.wikipedia.org/wiki/Advection) is\n",
    "\n",
    "$$\\partial_t q(x, t) = c \\, \\partial_x q(x, t)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's simulate this numerically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discretization using finite differences "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discretizing this with finite differences (first order in time and second order centered in space) leads to:\n",
    "\n",
    "$$\n",
    "q_i^ {n+1} = q_i^ {n} + \\frac{c \\Delta t}{2 \\Delta x} (q_{i+1}^n - q_{i-1}^n)\n",
    "$$\n",
    "\n",
    "On the edges, I used a first order discretization (I know, that's kind of silly...)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's implement the initial, conditions and the main loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = np.linspace(0, 1, num=151)\n",
    "dx = (x[1] - x[0]) / x.size\n",
    "dt = 0.00001\n",
    "c = -1.\n",
    "Courant = abs(c) * dt / dx\n",
    "print(f\"Courant number: {Courant}\")\n",
    "\n",
    "def f(x, x0=0.5):\n",
    "    return np.exp(-30 * (x - x0)**2)\n",
    "\n",
    "q0 = f(x)\n",
    "\n",
    "def step(q_curr, c, dt, dx):\n",
    "    q_next = q_curr.copy()\n",
    "    q_next[1:-1] += c * dt / (2 * dx) * (q_curr[2:] - q_curr[:-2])\n",
    "    q_next[0] += c * dt / dx * (q_curr[1] - q_curr[0]) \n",
    "    q_next[-1] += c * dt / dx * (q_curr[-1] - q_curr[-2]) \n",
    "    return q_next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def step_loop(q_curr, c, dt, dx):\n",
    "    q_next = np.zeros_like(q_curr)\n",
    "    for i in range(1, x.size - 1):\n",
    "        q_next[i] = q_curr[i] + c * dt / (2 * dx) * (q_curr[i+1] - q_curr[i-1])\n",
    "    # boundaries\n",
    "    q_next[0] = q_curr[0] + c * dt / dx * (q_curr[1] - q_curr[0]) \n",
    "    q_next[-1] = q_curr[-1] + c * dt / dx * (q_curr[-1] - q_curr[-2]) \n",
    "    return q_next"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_curr = q0.copy()\n",
    "q_next = np.zeros_like(q_curr)\n",
    "\n",
    "snapshots = {}\n",
    "for i in range(700):\n",
    "    t = i * dt\n",
    "    if i % 10 == 0:\n",
    "        snapshots[t] = q_curr.copy()    \n",
    "    q_next = step_loop(q_curr, c, dt, dx)\n",
    "    q_curr, q_next = q_next, q_curr\n",
    "    \n",
    "len(snapshots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's do an animated plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "line, = ax.plot([], [], label='numerical')\n",
    "line2, = ax.plot(x, q0, label='initial')\n",
    "\n",
    "def setup():\n",
    "    ax.set_xlim((x.min(), x.max()))\n",
    "    ax.set_ylim((-0.2, 1.2))\n",
    "    ax.legend(loc='upper right')\n",
    "    \n",
    "def update(frame):\n",
    "    key = list(snapshots.keys())[frame]\n",
    "    ydata = snapshots[key]\n",
    "    line.set_data(x, ydata)\n",
    "    ax.set_title(f'time: {key:.2e}')\n",
    "    return line, line2\n",
    "\n",
    "anim = FuncAnimation(fig, update, frames=range(len(snapshots)), init_func=setup)\n",
    "plt.close(fig)\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears that the scheme is not stable. It seems to be stable at first but the solutions grow and explode!\n",
    "\n",
    "\n",
    "To highlight this feature a little better, let's plot the maximum and minimum value over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(list(snapshots.keys()), [vals.max() for vals in snapshots.values()], label='max')\n",
    "ax.plot(list(snapshots.keys()), [vals.min() for vals in snapshots.values()], label='min')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What appears in this graph is that an exponential explosion happens at longer propagation times.\n",
    "\n",
    "Why is that?\n",
    "\n",
    "A good explanation by Langtangen is provided here: http://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book012.html#simplest-scheme-forward-in-time-centered-in-space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the explanation go? We consider a plane wave solution of the equation. "
   ]
  }
 ],
 "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.7.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}