13_finite_volume.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "\n",
    "Material heavily borrowed from David I. Ketcheson's [HyperPython](https://github.com/ketch/HyperPython) under a similar license as the one above.  Mistakes introduced are the responsibility of Kyle Mandli so complain to him.\n",
    "\n",
    "If you are interested in the animation tools used there is a tutorial on their use at \n",
    "\n",
    "https://nbviewer.jupyter.org/url/www.clawpack.org/_static/notebooks/animation_tools_demo.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import print_function\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Note that this module is found in the base directory\n",
    "import utils.animation as animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Finite Volume Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Hyperbolic Conservation Laws\n",
    "\n",
    "Before we specifically discuss finite volume methods we will explore the class of PDEs that are *hyperbolic conservation laws*.  These equations generally describe the movement of some physical quantity that should be conserved.  Examples of these quantities include mass, momentum, or energy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Conservation of Mass\n",
    "\n",
    "Imagine a fluid flowing in a narrow tube.  We'll use $q$ to indicate the density of the fluid and $u$ to indicate its velocity.  Both of these are functions of space and time: $q = q(x,t)$; $u=u(x,t)$.  The total mass in the section of tube $[x_1,x_2]$ is\n",
    "\n",
    "\\begin{equation}\n",
    "\\int_{x_1}^{x_2} q(x,t) dx.\n",
    "\\end{equation}\n",
    "\n",
    "This total mass can change in time due to fluid flowing in or out of this section of the tube.  We call the rate of flow the *flux*, and represent it with the function $f(q)$.  Thus the net rate of flow of mass into (or out of) the interval $[x_1,x_2]$ at time $t$ is\n",
    "\n",
    "$$f(q(x_1,t)) - f(q(x_2,t)).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We just said that this rate of flow must equal the time rate of change of total mass; i.e.\n",
    "$$\n",
    "    \\frac{d}{dt} \\int_{x_1}^{x_2} q(x,t) dx = f(q(x_1,t)) - f(q(x_2,t)).\n",
    "$$\n",
    "\n",
    "Assuming the functions are continuous (enough) we can write \n",
    "$$\n",
    "    \\int_{x_1}^{x_2} \\frac{\\partial}{\\partial x} f(q) dx = f(q(x_2,t)) - f(q(x_1,t))\n",
    "$$ \n",
    "and therefore rewrite out conservation law as\n",
    "$$\n",
    "    \\frac{d}{dt} \\int_{x_1}^{x_2} q(x,t) dx = -\\int_{x_1}^{x_2} \\frac{\\partial}{\\partial x} f(q) dx.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Under certain smoothness assumptions on $q$, we can move the time derivative inside the integral.  We'll also put everything on the left side, to obtain\n",
    "$$\n",
    "    \\int_{x_1}^{x_2} \\left(\\frac{\\partial}{\\partial t}q(x,t) + \\frac{\\partial}{\\partial x} f(q)\\right) dx = 0.\n",
    "$$\n",
    "\n",
    "Since this integral is zero for *any* choice of $x_1,x_2$, it must be that the integrand (the expression in parentheses) is actually zero *everywhere*!  Therefore we can write the **differential conservation law**\n",
    "$$\n",
    "    q(x,t)_t + f(q(x,t))_x = 0.\n",
    "$$\n",
    "\n",
    "This equation expresses the fact that the total mass is conserved -- since locally the mass can change only due to a net inflow or outflow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Advection \n",
    "\n",
    "In order to solve the conservation law above, we need an expression for the flux, $f$.  The rate of flow is just mass times velocity: $f=u q$.  Thus we obtain the **continuity equation**\n",
    "$$\n",
    "    q_t + (uq)_x = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In general, we need another equation to determine the velocity $u(x,t)$.  Later we will look at the full equations of fluid dynamics, but for now let's consider the simplest case, in which all of the fluid flows at a single, constant velocity $u(x,t)=a$.  Then the continuity equation becomes the *advection equation*\n",
    "\n",
    "$$q_t + a q_x = 0.$$\n",
    "\n",
    "We already saw that the solution to the advection equation is $$q(x,t) = q_0(x-at).$$  Let's plot the solution of the advection equation on the interval $[0,1]$ for the initial condition $$q_0(x) = e^{-2(x-1/2)^2}.$$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Grid\n",
    "x = numpy.linspace(0, 1, 1000) \n",
    "\n",
    "# Advection speed\n",
    "a = 1.0                    \n",
    "\n",
    "# Initial condition\n",
    "def q_0(x):\n",
    "    return numpy.exp(-200.*(x-0.2)**2)\n",
    "\n",
    "figs = []\n",
    "for t in numpy.linspace(0,1):\n",
    "    figs.append(plt.figure(figsize=(8, 4)))\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, q_0(x - a * t), linewidth=2)\n",
    "    axes.set_xlim((0, 1))\n",
    "    axes.set_ylim((-0.1, 1.1))\n",
    "    axes.set_xlabel(\"$x$\")\n",
    "    axes.set_ylabel(\"$q(x,t)$\")\n",
    "    axes.set_title(\"Solution to the advection equation\")\n",
    "    plt.close(figs[-1])\n",
    "    \n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As you can see, the initial pulse just moves to the right at speed $a$ as time advances.  This isn't very interesting, but it captures the most important feature of hyperbolic equations: waves travel at finite speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Finite Volume Method for Advection\n",
    "\n",
    "We can easily solve the advection equation exactly.  But the advection equation is a prototype for more complicated conservation laws that we will only be able to solve approximately by using numerical methods.  In order to better understand these methods, we will discuss them first in the context of the advection equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For simplicity, we'll suppose that we wish to solve the advection equation on the interval $[0,1]$.  We introduce a set of equally spaced *grid cells* of width $\\Delta x$, and write $x_i$ to mean the center of cell $i$.  Thus the first cell is the interval $[0,\\Delta x]$ and $x_1=\\Delta x/2$.  We will also write $x_{i-1/2}$ or $x_{i+1/2}$ to denote the left or right boundary of cell $i$, respectively.\n",
    "\n",
    "We write $Q_i$ to denote the *average* value of the solution over cell $i$:\n",
    "\n",
    "$$Q_i = \\frac{1}{\\Delta x} \\int_{x_{i-1/2}}^{x_{i-1/2}} q \\ dx.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The simplest finite volume method is obtained by supposing that the solution is actually *equal* to $Q_i$ over all of cell $i$.\n",
    "\n",
    "![](./images/finite_volume.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Suppose $a>0$.  Then the flux into cell $i$ from the left is $a Q_{i-1}$ and the flux out of cell $i$ to the right is $a Q_i$.  Then our integral conservation law reads\n",
    "\n",
    "$$Q_i'(t) = -\\frac{a}{\\Delta x}\\left(Q_i - Q_{i-1}\\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Applying a forward difference in time we obtain the *upwind method*\n",
    "$$\n",
    "    Q^{n+1}_i = Q^n_i -\\frac{a \\Delta t}{\\Delta x}\\left(Q_i - Q_{i-1}\\right).\n",
    "$$\n",
    "This of course is also the same upwind method we derived when considering finite differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "a = 1.0\n",
    "\n",
    "m = 100\n",
    "delta_x = 1.0 / float(m)\n",
    "x = numpy.arange(-delta_x / 2.0, 1.0 + delta_x / 2.0, delta_x)\n",
    "\n",
    "t = 0.0\n",
    "T = 0.5\n",
    "CFL = 0.8\n",
    "delta_t = CFL * delta_x / a\n",
    "\n",
    "# Initial data\n",
    "def q_0(x):\n",
    "    return numpy.exp(-200.*(x-0.2)**2)\n",
    "Q = q_0(x)\n",
    "Q_new = Q.copy()\n",
    "\n",
    "figs = []\n",
    "while t < T:\n",
    "    # Incoming wave\n",
    "#     Q_new[0] = numpy.sin(20.0 * numpy.pi * t)\n",
    "    # Extrapolation BCs\n",
    "    Q_new[0] = Q_new[1]\n",
    "    \n",
    "    Q_new[1:] = Q[1:] - a * delta_t / delta_x * (Q[1:] - Q[:-1])\n",
    "    \n",
    "    Q = Q_new.copy()\n",
    "    t = t + delta_t\n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, q_0(x - a * t), 'k', linewidth=2, label=\"True Solution\")\n",
    "    axes.plot(x, Q, 'o-', linewidth=2, label=\"Upwind\")\n",
    "    axes.set_title('t = %s' % str(t))\n",
    "    axes.set_xlim((x[0], x[-1]))\n",
    "    axes.set_ylim((-1.1, 1.1))\n",
    "    axes.set_xlabel(\"$x$\")\n",
    "    axes.set_ylabel(\"$q(x,t)$\")\n",
    "    axes.set_title(\"Upwind Method\")\n",
    "    axes.legend(loc=4)\n",
    "    plt.close(figs[-1])\n",
    "\n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can see that we have picked a $\\Delta t$ based on the CFL condition we derived for finite difference methods.  The same applies here (we are even using the same method).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Notice how we set up a grid that contains an extra cell at each end, outside of the problem domain $[0,1]$.  These are called *ghost cells* and are often useful in handling the solution at the grid boundaries.\n",
    "\n",
    "![](./images/ghost-cell.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The technique we have used to set the ghost cell values above, by copying the last value inside the grid to the ghost cells, is known as *zero-order extrapolation*.  It is useful for allowing waves to pass out of the domain (so-called *non-reflecting* boundaries).  Note that we don't actually need the ghost cell at the right end for the upwind method, but for other methods we will.\n",
    "\n",
    "The upwind method is simple, but it is not very accurate.  Notice how the computed solution becomes wider and shorter over time.  This behavior is referred to as *dissipation*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Exercises\n",
    "\n",
    "Now do the following with the code above:\n",
    "\n",
    "1. Set $m=1000$ or more and notice that it takes some time to compute the solution.  Rewrite the inner loop (over $i$) as a single line with no loop, using numpy slicing.  For large values of $m$, the code with slicing is much faster.\n",
    "1. Notice that the last step of the simulation goes past time $T$.  Modify the code so that the last step is adjusted to exactly reach $T$.\n",
    "2. Change the code so that animation of the solution versus time is plotted.  You will want to accumulate frames of the solution in a list and then use the same kind of code we used above to animate the exact solution.\n",
    "3. Add some code to plot the exact solution.\n",
    "4. Change the left boundary condition so that there is a sinusoidal wave coming in from the left:\n",
    "$$u(0,t) = \\sin(20 \\pi t).$$\n",
    "What do you notice about the sinusoid as it moves into the domain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Alternative Derivations\n",
    "\n",
    "Another way to look view finite volume methods is more directly as a discretization of the PDE:\n",
    "$$\n",
    "    \\frac{d}{dt} \\int_{x_1}^{x_2} q(x,t) dx = f(q(x_1,t)) - f(q(x_2,t)).\n",
    "$$\n",
    "to\n",
    "$$\n",
    "    Q_i^{n+1} = Q^n_i - \\frac{\\Delta t}{\\Delta x} \\left [ F^n_{i+1/2} - F^n_{i-1/2} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Stencil](./images/flux_stencil.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Averaging view of advection](./images/ave_demo.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If the flux from the left of cell $i$ says that the value $Q^n_{i-1}$ moves in at a speed $\\lambda$ then it will cover $\\lambda \\Delta t$ of the total call $i$'s length leaving $\\Delta x - \\lambda \\Delta t$ the original value $Q^n_i$.  Averaging these two values to find $Q^{n+1}_{i}$ leads to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    Q^{n+1}_i &= \\frac{Q^n_i (\\Delta x - \\lambda \\Delta t) + Q^n_{i-1} \\lambda \\Delta t}{\\Delta x} \\\\\n",
    "    &=Q^n_i - \\frac{\\lambda \\Delta t}{\\Delta x} [Q^n_i - Q^n_{i-1}],\n",
    "\\end{aligned}$$\n",
    "\n",
    "exactly the up-wind method!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Traffic Flow\n",
    "\n",
    "As an illustrative example we will now consider a particular scalar conservation law, that of a simplified model for traffic flow.  This model is sometimes referred to as the Lighthill-Whitham-Richards (or LWR) traffic model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "From the continuity equation above we know:\n",
    "\n",
    "$$q_t + (uq)_x = 0.$$\n",
    "\n",
    "Now we will think of $q$ as the density of cars on a road, traveling with velocity $u$.  Note that we're not keeping track of the individual cars, but just of the average number of cars per unit length of road.  Thus $q=0$ represents an empty stretch of road, and we can choose the units so that $q=1$ represents bumper-to-bumper traffic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We'll also choose units so that the speed limit is $u_\\text{max}=1$, and assume that drivers never go faster than this (yeah, right!)  If we assume that drivers always travel at a single uniform velocity, we obtain once again the advection equation.  But we all know that's not accurate in practice -- cars go faster in light traffic and slower when there is congestion.  The simplest way to incorporate this effect is to make the velocity a linearly decreasing function of the density:\n",
    "\n",
    "$$u(q) = 1 - q.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Notice that $u$ goes to zero as $q$ approaches the maximum density of 1, while $u$ goes to the maximum value of 1 as traffic density goes to zero.  Obviously, both $q$ and $u$ should always stay in the interval $[0,1]$.\n",
    "\n",
    "![](./images/LWR-Velocity.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Combining the two equations above, our conservation law says\n",
    "\n",
    "$$q_t + (q (1-q))_x = 0.$$\n",
    "\n",
    "The function $q(1-q)$ is the flux, or the rate of flow of cars.  Notice how the flux is zero when there are no cars and also when the road is completely full.  The maximum flow of traffic actually occurs when the road is half full:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "q = numpy.linspace(0, 1)\n",
    "f = q * (1.0 - q)\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(q,f,linewidth=2)\n",
    "axes.set_xlabel('q')\n",
    "axes.set_ylabel('$f(q)$')\n",
    "axes.set_xlim((0, 1))\n",
    "axes.set_ylim((0.0, 0.3))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This equation is fundamentally different from the advection equation because the flux is **nonlinear**.  This fact will have dramatic consequences for both the behavior of solutions and our numerical methods.  But we can superficially make this equation look like the advection equation by using the chain rule to write\n",
    "\n",
    "$$f(q)_x = f'(q) q_x = (1-2q)q_x.$$\n",
    "\n",
    "Then we have\n",
    "\n",
    "$$q_t + (1-2q)q_x = 0.$$\n",
    "\n",
    "This is like the advection equation, but with a velocity $1-2q$ that depends on the density of cars.  The value $f'(q)=1-2q$ is referred to as the *characteristic speed*.  This characteristic speed is not the speed at which cars move (notice that it can even be negative!)  Rather, it is the speed at which *information* is transmitted along the road."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Shock waves (traffic jams)\n",
    "\n",
    "Let's use the *Lax-Friedrichs method* to solve the LWR traffic model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "m = 100\n",
    "delta_x = 1.0 / m\n",
    "x = numpy.arange(-delta_x / 2.0, 1.0 + delta_x / 2.0, delta_x)\n",
    "\n",
    "t = 0.0\n",
    "T = 0.5\n",
    "CFL = 0.9\n",
    "delta_t = CFL * delta_x\n",
    "\n",
    "Q = 0.9 * numpy.exp(-100.0 * (x - 0.5)**2)\n",
    "Q_new = numpy.empty(Q.shape)\n",
    "Q_list = [Q]\n",
    "\n",
    "figs = []\n",
    "while t < T:\n",
    "    \n",
    "    Q_new[1:-1] = 0.5 * (Q[:-2] + Q[2:]) - 0.5 * delta_t / delta_x \\\n",
    "                                               * (Q[2:] * (1.0 - Q[2:]) - Q[:-2] * (1.0 - Q[:-2]))\n",
    "    \n",
    "    Q = Q_new.copy()\n",
    "    Q[0] = Q[-2]\n",
    "    Q[-1] = Q[1]\n",
    "    t += delta_t\n",
    "    \n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, Q, 'o-', linewidth=2)\n",
    "    axes.set_title('Lax-Friedrichs Method - Traffic Flow at t = %s' % str(t))\n",
    "    axes.set_xlabel(\"$x$\")\n",
    "    axes.set_ylabel(\"$q(x,t)$\")\n",
    "    axes.set_xlim((x[0], x[-1]))\n",
    "    axes.set_ylim((-0.1, 1.1))\n",
    "    plt.close(figs[-1])\n",
    "    \n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Starting from an initially smooth distribution of cars, we see that the solution evolves into something discontinuous!  What's going on?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Naturally, the cars in the low-density region behind the peak initially drive faster than the cars in the most congested region, ahead of them.  This leads to a steepening of the rear part of the traffic flow, and eventually there is a traffic jam: cars coming from the back have an open road, until suddenly they reach a highly congested area and must slam on the brakes.  Sound familiar?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Exercises\n",
    "\n",
    "1. Recall our discussion of the CFL condition for advection.  To apply the CFL condition to a nonlinear problem, we use the characteristic velocity (in this case $1-2q$) in place of the velocity $a$.  For this problem, how large can we take the time step without violating the CFL condition? \n",
    "\n",
    "1. In the cell below, set things up to find out what happens when cars reach a red light.  Place the red light at $x=1$ and set $q=1$ in the interval $0.8<x<1$ near the light (the cars will be stopped, bumper to bumper there).  In the rest of the domain, set $q=0.2$, corresponding to light traffic.  In the ghost cell at the right boundary, set $q=1$; in the ghost cell at the left set $q=0.2$.  After you have computed the solution, check it against your inuition.  Does it make sense?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Speed of a shock wave: the Rankine-Hugoniot conditions\n",
    "\n",
    "In the second exercise above, you should have seen a shock wave (i.e., a discontinuity) that moves to the left as more and more cars pile up behind the red light.  How quickly does this discontinuity move to the left?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can figure it out by putting an imaginary line at the location of the shock.  Let $q_l$ be the density of cars just to the left of the line, and let $q_r$ be the density of cars just to the right.  Imagine for a moment that the line is stationary.  Then the rate of cars reaching the line from the left is $f(q_l)$ and the rate of cars departing from the line to the right is $f(q_r)$.  If the line really were stationary, we would need to have $f(q_l)-f(q_r)=0$ to avoid cars accumulating at the line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![](./images/shock_diagram_traffic_a.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "However, the shock is not stationary, so the line is moving.  Let $s$ be the speed of the shock.  Then as the line moves to the left, some cars that were to the left are now to the right of the line.  The rate of cars removed from the left is $s q_l$ and the rate of cars added on the right is $s q_r$.  So in order to avoid an accumulation of cars at the shock, these two effects need to be balanced:\n",
    "$$\n",
    "    f(q_l) - f(q_r) = s(q_l - q_r).\n",
    "$$\n",
    "\n",
    "This condition is known as the *Rankine-Hugoniot condition*, and it holds for any shock wave in the solution of any hyperbolic system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![](./images/shock_diagram_traffic_b.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another more analytical way to get at the Rankine-Hugoniot conditions is to consider a rectangle $[x_1, x_1 + \\Delta x] \\times [t_1, t_1 + \\Delta t]$ which is split by a shock with speed $s$.\n",
    "\n",
    "Applying the integral conservation law to this region assuming that $q$ is constant in each region we have\n",
    "$$\\begin{aligned}\n",
    "    \\int^{x_1 + \\Delta x}_{x_1} q(x, t_1 + \\Delta t) dx - \\int^{x_1 + \\Delta x}_{x_1} q(x, t_1) dx &= \\int^{t_1 + \\Delta t}_{t_1} q(x_1, t) dt - \\int^{t_1 + \\Delta t}_{t_1} q(x_1 + \\Delta x, t) dt \\\\\n",
    "    \\Delta x q_r - \\Delta x q_\\ell &= \\Delta t f(q_\\ell) - \\Delta t f(q_r) + \\mathcal{O}(\\Delta t^2) \\\\\n",
    "    \\frac{\\Delta x}{\\Delta t} (q_r - q_\\ell) &=  (f(q_\\ell) -  f(q_r)) + \\mathcal{O}(\\Delta t)\n",
    "\\end{aligned}$$\n",
    "then letting $\\Delta t \\rightarrow 0$ we have\n",
    "$$\n",
    "    s (q_r - q_\\ell) =  f(q_r) -  f(q_\\ell)\n",
    "$$\n",
    "noting that $-s \\Delta t = \\Delta x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Exercises\n",
    "\n",
    "1. Recall our discussion of the CFL condition for advection.  To apply the CFL condition to a nonlinear problem, we use the characteristic velocity (in this case $1-2q$) in place of the velocity $a$.  For this problem, how large can we take the time step without violating the CFL condition? \n",
    "\n",
    "1. In the cell below, set things up to find out what happens when cars reach a red light.  Place the red light at $x=1$ and set $q=1$ in the interval $0.8<x<1$ near the light (the cars will be stopped, bumper to bumper there).  In the rest of the domain, set $q=0.2$, corresponding to light traffic.  In the ghost cell at the right boundary, set $q=1$; in the ghost cell at the left set $q=0.2$.  After you have computed the solution, check it against your inuition.  Does it make sense?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Rarefaction Waves\n",
    "\n",
    "Why didn't we see a shock wave in the third exercise above?  Initially, the solution is discontinuous, but after time zero it becomes continuous.\n",
    "\n",
    "Physically, the reason is clear: after the light turns green, the cars in front accelerate and spread out; then the cars behind them accelerate, and so forth.  This kind of expansion wave is referred to as a *rarefaction wave* by analogy to fluid dynamics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "How can we determine whether an initial discontinuity will lead to a shock or a rarefaction?\n",
    "- Shocks appear in regions where characteristics overlap:\n",
    "\n",
    "![](./images/entropy_condition_shock.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Rarefactions appear in regions where characteristics are spreading out:\n",
    "\n",
    "![](./images/entropy_condition_rarefaction.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "More precisely, if the value to the left of a shock is $q_l$ and the value to the right is $q_r$, then it must be that $f'(q_l)>f'(q_r)$.  In fact the shock speed must lie between these characteristic speeds:\n",
    "\n",
    "$$f'(q_l) > s > f'(q_r).$$\n",
    "\n",
    "We say that the characteristics *impinge* on the shock.  This is known as the *entropy condition*, because in fluid dynamics such a shock obeys the 2nd law of thermodynamics.\n",
    "\n",
    "On the other hand, if $f'(q_l)< f'(q_r)$, then a rarefaction wave results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Riemann Problem\n",
    "\n",
    "The problems you set up in the red and green light exercises above are examples of what is known as a Riemann problem.  A Riemann problem consists of  a hyperbolic conservation law, together with piecewise-constant initial data with a single point of discontinuity:\n",
    "\n",
    "$$q(x,0) = \\begin{cases} q_l & x < x_0 \\\\ q_r & x > x_0 \\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In fact, our entire finite volume method is built on approximately solving Riemann problems.  Since the solution is approximated by a constant function in each cell, a Riemann problem arises at each interface between two cells.  The upwind method for advection is based on exactly solving the Riemann problem to find the flux (of course, the Riemann problem for advection is extremely simple to solve).  The Lax-Friedrichs flux is based on approximating the solution to the Riemann problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For nonlinear systems of hyperbolic PDEs, like the Euler equations that we will consider later, the solution of the Riemann problem can be quite complicated.  Computational methods for solving such systems are usually based on approximate Riemann solvers.  We won't go into the theory of Riemann problems or approximate Riemann solvers in this course."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## High-Resolution Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The methods we have used so far (the *upwind method* and the *Lax-Friedrichs method*) are both dissipative.  This means that over time they artificially smear out the solution -- especially shocks.  Furthermore, both of these methods are only *first order accurate*, meaning that if we reduce the values of $\\Delta t$ and $\\Delta x$ by a factor of two, the overall error decreases (only) by a factor of two.  We can do better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Reducing Diffusion (Dissipation)\n",
    "\n",
    "The first step in improving our accuracy is to use a more accurate representation of $q(x)$ at each step.  Instead of assuming that $q$ is piecewise-constant, we will now approximate it by a piecewise-linear function:\n",
    "\n",
    "$$q(x) = Q^n_i + \\sigma^n_i (x-x_i).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here $\\sigma_i$ represents the slope in cell $i$.  The most obvious choice to ensure that this results in a second-order accurate approximation is to take the centered approximation\n",
    "\n",
    "$$\\sigma^n_i = \\frac{Q^n_{i+1} - Q^n_{i-1}}{2\\Delta x}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We use this to obtain values at the cell interfaces:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "q^+_{i - 1/2} & = Q_i - \\sigma_i \\frac{\\Delta x}{2} \\\\\n",
    "q^-_{i + 1/2} & = Q_i + \\sigma_i \\frac{\\Delta x}{2}.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![](./images/linear_reconstruction.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We'll use these interface values to approximate the flux, based on the *Lax-Friedrichs flux* (related to the Lax-Friedrichs method):\n",
    "\n",
    "$$F_{i - 1/2} = \\frac{1}{2} \\left( f(q^-_{i - 1/2}) + f(q^+_{i - 1/2}) - \\frac{\\Delta t}{\\Delta x} (q^+_{i - 1/2} - q^-_{i - 1/2})\\right)$$\n",
    "\n",
    "This provides second-order accuracy in space.  We also need to make the method second-order accurate in time.  We can do so by using a second-order Runge--Kutta method.  Then the full method is\n",
    "\n",
    "$$\\begin{aligned}\n",
    "Q^*_i & = Q^n_i - \\frac{\\Delta t}{\\Delta x}\\left( F^n_{i - 1/2} - F^n_{i - 1/2} \\right) \\\\\n",
    "Q^{n+1}_i & = \\frac{1}{2} Q^n_i + \\frac{1}{2}\\left( Q^*_i - \\frac{\\Delta t}{\\Delta x}\\left( F^*_{i + 1/2} - F^*_{i - 1/2} \\right) \\right) \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def f(q):\n",
    "    return q * (1.0 - q)\n",
    "\n",
    "m = 100\n",
    "delta_x = 1.0 / float(m)\n",
    "x = numpy.arange(-3.0 * delta_x / 2.0, 1.0 + 5.0 * delta_x / 2.0, delta_x)\n",
    "\n",
    "t = 0. # Initial time\n",
    "T = 0.5 # Final time\n",
    "CFL = 0.9\n",
    "delta_t = CFL * delta_x  # Time step\n",
    "\n",
    "Q_HR = 0.9 * numpy.exp(-100.0 * (x - 0.5)**2)\n",
    "Q_new_HR = numpy.zeros(m + 4)\n",
    "Q_star = numpy.zeros(m + 4)\n",
    "sigma = numpy.zeros(m + 4)\n",
    "F = numpy.zeros(m + 4)\n",
    "\n",
    "Q_LF = 0.9 * numpy.exp(-100.0 * (x - 0.5)**2)\n",
    "Q_new_LF = numpy.zeros(m + 4)\n",
    "\n",
    "figs = []\n",
    "while t < T:\n",
    "    \n",
    "    # High-resolution version\n",
    "    sigma[1:-1]  = (Q_HR[2:] - Q_HR[:-2]) / (2.0 * delta_x)\n",
    "    q_plus  = Q_HR[1:-1] - sigma[1:-1] * delta_x / 2.0\n",
    "    q_minus = Q_HR[:-2] + sigma[:-2]  * delta_x / 2.0\n",
    "    F[1:-1] = 0.5 * (f(q_plus) + f(q_minus) - delta_x / delta_t * (q_plus - q_minus) )\n",
    "    \n",
    "    Q_star[2:-2] = Q_HR[2:-2] - delta_t / delta_x * (F[3:-1] - F[2:-2])\n",
    "    Q_star[0:2] = Q_star[2]\n",
    "    Q_star[-2:] = Q_star[-3]\n",
    "    \n",
    "    sigma[1:-1]  = (Q_star[2:] - Q_star[:-2]) / (2.0 * delta_x)\n",
    "    q_plus  = Q_star[1:-1] - sigma[1:-1] * delta_x / 2.0\n",
    "    q_minus = Q_star[:-2] + sigma[:-2]  * delta_x / 2.0\n",
    "    F[1:-1] = 0.5 * (f(q_plus) + f(q_minus) - delta_x / delta_t * (q_plus - q_minus))\n",
    "    \n",
    "    Q_new_HR[2:-2] = 0.5 * Q_HR[2:-2] + 0.5 * (Q_star[2:-2] - delta_t / delta_x * (F[3:-1] - F[2:-2]))\n",
    "        \n",
    "    Q_HR = Q_new_HR.copy()\n",
    "    Q_HR[0:2] = Q_HR[2]\n",
    "    Q_HR[-2:] = Q_HR[-3]\n",
    "    t += delta_t\n",
    "    \n",
    "    # LF Version\n",
    "    Q_new_LF[1:-1] = 0.5 * (Q_LF[:-2] + Q_LF[2:]) - 0.5 * delta_t / delta_x * (Q_LF[2:] *(1.0 - Q_LF[2:]) - Q_LF[:-2] * (1.0 - Q_LF[:-2]))\n",
    "    \n",
    "    Q_LF = Q_new_LF.copy()\n",
    "    Q_LF[0] = Q_LF[-2]\n",
    "    Q_LF[-1] = Q_LF[1]\n",
    "    t += delta_t\n",
    "    \n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, Q_HR, 'o-', linewidth=2, label=\"High Resolution\")\n",
    "    axes.plot(x, Q_LF, 'x', linewidth=2, label=\"Original\")\n",
    "    axes.set_title('t = %s' % str(t))\n",
    "    axes.set_xlim((x[0], x[-1]))\n",
    "    axes.set_ylim((-0.1, 1.1))\n",
    "    axes.legend()\n",
    "    plt.close(figs[-1])\n",
    "    \n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The shock wave is much sharper now, but we have a new problem. Do you see the little dip behind the shock?  If you look closely, you'll see that the solution is actually negative there!  Obviously, a negative density of cars makes no sense.  What's more, the solution shouldn't dip there -- and it shouldn't have that funny bump just in front of the shock either.  What to do?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Slope Limiting\n",
    "\n",
    "The spurious oscillations in our solution are not a particular feature of the method we've chosen.  In fact, *any* second-order (or higher) method that computes $q^\\pm_{i + 1/2}$ as a linear function of the cell averages will have oscillations (this is a famous result known as *Godunov's Theorem*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can get around this difficulty by choosing the slope $\\sigma$ as a *nonlinear* function of the cell averages.  In particular, to avoid oscillations we can choose the smaller of the two one-sided slopes.  Let $\\Delta Q_{i - \\Delta x} = Q_i - Q_{i-1}$.  Then we use the slope\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\sigma_i & = \\text{minmod}(\\Delta Q_{i-1/2},\\Delta Q_{i + 1/2})/\\Delta x \\\\\n",
    "& = \\begin{cases} \\min(\\Delta Q_{i-1/2}, \\Delta Q_{i + 1/2})/\\Delta x & \\text{ if } \\Delta Q_{i-1/2}, \\Delta Q_{i + 1/2} > 0 \\\\\n",
    "\\max(\\Delta Q_{i - 1/2}, \\Delta Q_{i + 1/2})/\\Delta x & \\text{ if } \\Delta Q_{i-1/2}, \\Delta Q_{i + 1/2} < 0 \\\\\n",
    "0 & \\text{ if } \\Delta Q_{i-1/2}\\cdot \\Delta Q_{i + 1/2} < 0.\n",
    "\\end{cases}\n",
    "\\end{aligned}$$\n",
    "\n",
    "This choice of slope is known as the minimum-modulus, or *minmod* slope."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Local Lax-Friedrichs flux\n",
    "\n",
    "The Lax-Friedrichs flux ensures that our solution is stable, but it does so by adding a lot of dissipation everywhere.  In fact, we could get away with using less dissipation over most of the domain.  A variant that does this is called the *local Lax-Friedrichs flux*.  It is little more accurate than the Lax-Friedrichs flux because it uses the local characteristic speeds to determine how much dissipation is needed at each interface.  It is\n",
    "\n",
    "$$F_{i - \\Delta x} = \\frac{1}{2} \\left( f(q^-_{i - 1/2}) + f(q^+_{i - 1/2}) - \\alpha_{i - 1/2}\\frac{\\Delta x}{\\Delta t} (q^+_{i - 1/2} - q^-_{i - 1/2})\\right)$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\alpha_{i - 1/2} = \\min(\\left|f'(q^-_{i - 1/2})\\right|,\\left|f'(q^+_{i - 1/2})\\right|).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Exercise\n",
    "\n",
    "In the cell below, copy and modify the second-order method above to use the minmod slope and local Lax-Friedrichs flux.  \n",
    "\n",
    "*Hint 1*: You will want to use the functions `numpy.minimum` and `numpy.maximum`, which compare two arrays elementwise (not `numpy.min`, which finds the minimum of a single array).\n",
    "\n",
    "*Hint 2*: to avoid using a loop for the slope computation, note that the minmod function can be written as\n",
    "$$\n",
    "\\text{minmod}(a,b) = \\frac{\\text{sign}(a)+\\text{sign}(b)}{2} \\min(|a|,|b|).\n",
    "$$\n",
    "The signum function is implemented as `numpy.sign()`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def minmod(a, b):\n",
    "    return (numpy.sign(a) + numpy.sign(b)) / 2.0 * numpy.minimum(numpy.abs(a), numpy.abs(b))\n",
    "\n",
    "def f(q):\n",
    "    return q * (1.0 - q)\n",
    "\n",
    "def f_prime(q):\n",
    "    return 1.0 - 2.0 * q\n",
    "\n",
    "m = 100\n",
    "delta_x = 1.0 / float(m)\n",
    "x = numpy.arange(-3.0 * delta_x / 2.0, 1.0 + 5.0 * delta_x / 2.0, delta_x)\n",
    "\n",
    "t = 0. # Initial time\n",
    "T = 0.5 # Final time\n",
    "# T = delta_t * 20\n",
    "CFL = 0.9\n",
    "delta_t = CFL * delta_x  # Time step\n",
    "\n",
    "Q = 0.9 * numpy.exp(-100.0 * (x - 0.5)**2)\n",
    "Q_new = numpy.zeros(m + 4)\n",
    "Q_star = numpy.zeros(m + 4)\n",
    "sigma = numpy.zeros(m + 4)\n",
    "F = numpy.zeros(m + 4)\n",
    "\n",
    "Q_LF = 0.9 * numpy.exp(-100.0 * (x - 0.5)**2)\n",
    "Q_new_LF = numpy.zeros(m + 4)\n",
    "\n",
    "figs = []\n",
    "while t < T:\n",
    "    \n",
    "    sigma[1:-1]  = minmod(Q[1:-1] - Q[:-2], Q[2:] - Q[1:-1]) / delta_x\n",
    "    Q_plus  = Q[1:-1] - sigma[1:-1] * delta_x / 2.0\n",
    "    Q_minus = Q[:-2] + sigma[:-2]  * delta_x / 2.0\n",
    "    alpha = numpy.minimum(numpy.abs(f_prime(Q_minus)), numpy.abs(f_prime(Q_plus)))\n",
    "    F[1:-1] = 0.5 * (f(Q_plus) + f(Q_minus) - alpha * delta_x / delta_t * (Q_plus - Q_minus) )\n",
    "    \n",
    "    Q_star[2:-2] = Q[2:-2] - delta_t / delta_x * (F[3:-1] - F[2:-2])\n",
    "    Q_star[0:2] = Q_star[2]\n",
    "    Q_star[-2:] = Q_star[-3]\n",
    "    \n",
    "    sigma[1:-1]  = (Q_star[2:] - Q_star[:-2]) / (2.0 * delta_x)\n",
    "    q_plus  = Q_star[1:-1] - sigma[1:-1] * delta_x / 2.0\n",
    "    q_minus = Q_star[:-2] + sigma[:-2]  * delta_x / 2.0\n",
    "    F[1:-1] = 0.5 * (f(Q_plus) + f(Q_minus) - delta_x / delta_t * (Q_plus - Q_minus))\n",
    "    \n",
    "    Q_new[2:-2] = 0.5 * Q[2:-2] + 0.5 * (Q_star[2:-2] - delta_t / delta_x * (F[3:-1] - F[2:-2]))\n",
    "        \n",
    "    Q = Q_new.copy()\n",
    "    Q[0:2] = Q[2]\n",
    "    Q[-2:] = Q[-3]\n",
    "    t += delta_t\n",
    "        \n",
    "    # LF Version\n",
    "    Q_new_LF[1:-1] = 0.5 * (Q_LF[:-2] + Q_LF[2:]) - 0.5 * delta_t / delta_x * (Q_LF[2:] *(1.0 - Q_LF[2:]) - Q_LF[:-2] * (1.0 - Q_LF[:-2]))\n",
    "    \n",
    "    Q_LF = Q_new_LF.copy()\n",
    "    Q_LF[0] = Q_LF[-2]\n",
    "    Q_LF[-1] = Q_LF[1]\n",
    "    t += delta_t\n",
    "    \n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, Q, 'o-', linewidth=2, label=\"High Resolution\")\n",
    "    axes.plot(x, Q_LF, 'x', linewidth=2, label=\"Original\")\n",
    "    axes.set_title('t = %s' % str(t))\n",
    "    axes.set_xlim((x[0], x[-1]))\n",
    "    axes.set_ylim((-0.1, 1.1))\n",
    "    axes.legend()\n",
    "    plt.close(figs[-1])\n",
    "    \n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As you can see, this method keeps the shock fairly sharp and avoids the creation of negative solution values.  This method falls into the class of [MUSCL]() schemes and is proven to avoid oscillations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### High-Order Methods\n",
    "\n",
    "Just like the method we implemented above, most methods that are more than first-order accurate consist of three components:\n",
    "1. **Reconstruction**: A method for computing interface values $q^\\pm_{i - \\Delta x}$ from cell averages $Q_i$.  This involves some kind of limiting in order to avoid oscillations.  Higher-order reconstruction is often done using [weighted essentially non-oscillatory (WENO)](http://www.dam.brown.edu/scicomp/media/report_files/BrownSC-2006-21.ps.gz) methods.\n",
    "2. **Numerical flux**:  An approximation of the flux, computed based on the interface values $q^\\pm_{i - \\Delta x}$.  The Lax-Friedrichs flux above is one of the simplest.  Much more accurate fluxes can be computed using an exact or approximate [Riemann solver](http://en.wikipedia.org/wiki/Riemann_solver).\n",
    "3. **Time integrator**: In order to get high-order accuracy in time, usually a [Runge-Kutta method](http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) is used.  [Strong stability preserving methods](http://www.davidketcheson.info/assets/papers/sspreview.pdf) are particularly popular.\n",
    "\n",
    "It's possible to devise methods of very high order by increasing the order of the polynomial reconstruction and of the time integrator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Reconstruction](./images/reconstruction_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Reconstruction-constants](./images/reconstruction_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Reconstruction-linear](./images/reconstruction_3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Reconstruction-quadratic](./images/reconstruction_4.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are many different ways to reconstruct and limit solutions, try out a few.  This requires a Clawpack installation, see the [Clawpack webpage](http://www.clawpack.org) for installation instructions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Demo limiters\n",
    "from clawpack import pyclaw\n",
    "from clawpack import riemann\n",
    "\n",
    "# Pick out Riemann solver\n",
    "riemann_solver = riemann.advection_1D_py.advection_1D\n",
    "\n",
    "# Create solver\n",
    "solver = pyclaw.ClawSolver1D(riemann_solver)\n",
    "solver.kernel_language = \"Python\"\n",
    "solver.bc_lower[0] = pyclaw.BC.periodic\n",
    "solver.bc_upper[0] = pyclaw.BC.periodic\n",
    "\n",
    "# Limiters avaiable\n",
    "solver.limiters = 4\n",
    "# To pick a Godunov upwind method set order to 1\n",
    "solver.order = 2\n",
    "limiter_name = [\"None\", \"minmod\", \"superbee\", \"Van Leer\", \"MC\", \"Beam-Warming\", \"Fromm\",\n",
    "                \"Albada 2\", \"Albada 3\", \"Van Leer with Klein sharpening\"]\n",
    "\n",
    "# Set up domain\n",
    "x = pyclaw.Dimension(0.0, 1.0, 100, name='x')\n",
    "domain = pyclaw.Domain(x)\n",
    "state = pyclaw.State(domain, solver.num_eqn)\n",
    "\n",
    "# Problem specific data\n",
    "state.problem_data['u'] = 1.0\n",
    "\n",
    "# Initial condition\n",
    "xc = state.grid.x.centers\n",
    "beta = 100; gamma = 0; x0 = 0.75\n",
    "state.q[0, :] = numpy.exp(-beta * (xc-x0)**2) * numpy.cos(gamma * (xc - x0))\n",
    "state.q[0, :] += 1.0 * (0.2 < xc) * (xc < 0.4)\n",
    "\n",
    "# Set up time stepping and other run parameters\n",
    "claw = pyclaw.Controller()\n",
    "claw.keep_copy = True\n",
    "claw.solution = pyclaw.Solution(state, domain)\n",
    "claw.solver = solver\n",
    "claw.num_output_times = 30\n",
    "claw.tfinal = 30.0\n",
    "\n",
    "claw.run()\n",
    "\n",
    "figs = []\n",
    "for frame in claw.frames:\n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(xc, frame.q[0, :], 'o-', linewidth=2)\n",
    "    axes.plot(xc, claw.frames[0].q[0, :], 'k--')\n",
    "    axes.set_xlim((0, 1))\n",
    "    axes.set_ylim((-0.2, 1.2))\n",
    "    axes.set_title(\"Limiter = %s\" % limiter_name[solver.limiters])\n",
    "    plt.close(figs[-1])\n",
    "    \n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Fluid Dynamics\n",
    "\n",
    "As a more complex example of a system of conservation laws we will now consider the Euler equations of compressible gas dynamics.  These consist of conservation laws for *mass, momentum*, and *energy*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Mass Conservation\n",
    "\n",
    "We will use $\\rho(x,t)$ to denote the fluid density and $u(x,t)$ for its velocity.  Then the equation for conservation of mass is just the **continuity equation** we discussed in [Lesson 1](Lesson_01_Advection.ipynb):\n",
    "\n",
    "$$\\rho_t + (\\rho u)_x = 0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Momentum Conservation\n",
    "\n",
    "The momentum is given by the product of density and velocity, i.e. $\\rho u$.  The momentum flux has two components.  First, the momentum is transported in the same way that the density is; this flux is given by the momentum times the density; i.e. $\\rho u^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To understand the second term in the momentum flux, we must realize that a fluid is made up of many tiny molecules.  The density and velocity we are modeling are average values over some small region of space.  The individual molecules in that region are not all moving with exactly velocity $u$; that's just their average.  Each molecule also has some additional random velocity component.  These random velocities are what accounts for the **pressure** of the fluid, which we'll denote by $p$.  These velocity components also lead to a net flux of momentum.  Thus the momentum conservation equation is\n",
    "\n",
    "$$(\\rho u)_t + (\\rho u^2 + p)_x = 0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Energy Conservation\n",
    "\n",
    "The energy has two components: internal energy $\\rho e$ and kinetic energy $\\rho u^2/2$:\n",
    "\n",
    "$$E = \\rho e + \\frac{1}{2}\\rho u^2.$$\n",
    "\n",
    "Like the momentum flux, the energy flux involves both bulk transport ($Eu$) and transport due to pressure ($pu$):\n",
    "\n",
    "$$E_t + (u(E+p)) = 0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Equation of State\n",
    "\n",
    "You may have noticed that we have 4 unknowns (density, momentum, energy, and pressure) but only 3 conservation laws.  We need one more relation to close the system.  That relation, known as the equation of state, expresses how the pressure is related to the other quantities.  We'll focus on the case of an ideal gas, for which\n",
    "\n",
    "$$p = \\rho e (\\gamma-1).$$\n",
    "\n",
    "Here $\\gamma$ is the ratio of specific heats, which for air is approximately 1.4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Full Equations\n",
    "\n",
    "We can write the three conservation laws as a single system $q_t + f(q)_x = 0$ by defining\n",
    "\\begin{align}\n",
    "q & = \\begin{pmatrix} \\rho \\\\ \\rho u \\\\ E\\end{pmatrix}, & \n",
    "f(q) & = \\begin{pmatrix} \\rho u \\\\ \\rho u^2 + p \\\\ u(E+p)\\end{pmatrix}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In three dimensions, the equations are similar.  We have two additional velocity components $v, w$, and their corresponding fluxes.  Additionally, we have to account for fluxes in the $y$ and $z$ directions.  We can write the full system as\n",
    "\n",
    "$$ q_t + f(q)_x + g(q)_y + h(q)_z = 0$$\n",
    "\n",
    "with\n",
    "\n",
    "\\begin{align}\n",
    "q & = \\begin{pmatrix} \\rho \\\\ \\rho u \\\\ \\rho v \\\\ \\rho w \\\\ E\\end{pmatrix}, &\n",
    "f(q) & = \\begin{pmatrix} \\rho u \\\\ \\rho u^2 + p \\\\ \\rho u v \\\\ \\rho u w \\\\ u(E+p)\\end{pmatrix} &\n",
    "g(q) & = \\begin{pmatrix} \\rho v \\\\ \\rho uv \\\\ \\rho v^2 + p \\\\ \\rho v w \\\\ v(E+p)\\end{pmatrix} &\n",
    "h(q) & = \\begin{pmatrix} \\rho w \\\\ \\rho uw \\\\ \\rho vw \\\\ \\rho w^2 + p \\\\ w(E+p)\\end{pmatrix}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving the Euler Equations\n",
    "\n",
    "These equations can be solved in a manner similar to what we used for advection and traffic flow.  As you might guess, computing the flux gets significantly more complicated since we now have 3 (or 5) equations and more complicated flux expressions.\n",
    "\n",
    "Implementing a solver for the Euler equations from scratch would be a lot of fun, but to save some time we'll use a  package called [PyClaw](http://clawpack.github.io/doc/pyclaw/), which is part of the [Clawpack](http://clawpack.github.io/) software (Clawpack stands for Conservation LAWs PACKage).  PyClaw allows us to quickly and easily set up and solve problems modeled by hyperbolic PDEs.\n",
    "\n",
    "Now let's get started.  First, import the parts of Clawpack that we'll use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from clawpack import pyclaw\n",
    "from clawpack import riemann"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Setting up a problem\n",
    "To solve a problem, we'll need to create the following:\n",
    "\n",
    "- A **domain** over which to solve the problem\n",
    "- A **solution**, where we will provide the initial data.  After running, the solution will contain -- you guessed it! -- the solution.\n",
    "- A **solver**, which is responsible for actually evolving the solution in time.  Here we'll need to specify the equations to be solved and the boundary conditions.\n",
    "- A **controller**, which handles the running, output, and can be used for plotting\n",
    "\n",
    "This might sound complicated at first, but stick with me."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's start by creating a controller and specifying the simulation end time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "claw = pyclaw.Controller()\n",
    "claw.tfinal = 0.1                    # Simulation end time\n",
    "\n",
    "claw.keep_copy = True          # Keep solution data in memory for plotting\n",
    "claw.output_format = None    # Don't write solution data to file\n",
    "claw.num_output_times = 50  # Write 50 output frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Riemann solvers\n",
    "\n",
    "The method used to compute the flux between each pair of cells is referred to as a *Riemann solver*.  By specifying a Riemann solver, we will specify the system of PDEs that we want to solve.  So far we have only used very simple approximate Riemann solvers.  Clawpack includes much more sophisticated Riemann solvers for many hyperbolic systems.\n",
    "\n",
    "Place your cursor at the end of the line in the box below and hit 'Tab' (for autocompletion).  You'll see a dropdown list of all the Riemann solvers currently available in PyClaw.  The ones with 'py' at the end of the name are written in pure Python; the others are written in Fortran and wrapped with f2py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dir(riemann)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We'll start with a simple 1D problem, using the Riemann solver `riemann.euler_with_efix_1D`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "riemann_solver = riemann.euler_with_efix_1D\n",
    "claw.solver = pyclaw.ClawSolver1D(riemann_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We also need to specify boundary conditions.  We'll use extrapolation boundary conditions, so that waves simply pass out of the domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "claw.solver.all_bcs = pyclaw.BC.extrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The problem domain\n",
    "Next we need to specify the domain and the grid.  We'll solve on the unit line $[0,1]$ using 100 grid cells.  Note that each argument to the Domain constructor must be a tuple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "domain = pyclaw.Domain( (0.,), (1.,), (100,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The initial solution\n",
    "Next we create a solution object that belongs to the controller and extends over the domain we specified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The initial solution\n",
    "Next we create a solution object that belongs to the controller and extends over the domain we specified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "claw.solution = pyclaw.Solution(3,domain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The initial data is specified in an array named `solution.q`.  The density is contained in `q[0,:]`, the momentum in `q[1,:]`, and the energy in `q[2,:]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x=domain.grid.x.centers # grid cell centers\n",
    "gam = 1.4 # ratio of specific heats\n",
    "\n",
    "rho_left = 1.0; rho_right = 0.125\n",
    "p_left = 1.0; p_right = 0.1\n",
    "\n",
    "claw.solution.q[0,:] = (x<0.5)*rho_left + (x>=0.5)*rho_right\n",
    "claw.solution.q[1,:] = 0.\n",
    "claw.solution.q[2,:] = ((x<0.5)*p_left + (x>=0.5)*p_right)/(gam-1.0)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, claw.solution.q[0,:],'-o')\n",
    "axes.set_xlim((0.0, 1.0))\n",
    "axes.set_ylim((-0.1, 1.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This problem is known as the *Sod shock-tube*.  It amounts to setting up a tube with a thin separator between a high-pressure, high-density region and a low-pressure, low-density region, then suddenly removing the separator.\n",
    "\n",
    "Next we need to specify the value of $\\gamma$, the ratio of specific heats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "problem_data = claw.solution.problem_data\n",
    "problem_data['gamma'] = 1.4\n",
    "problem_data['gamma1'] = 0.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Finally, let's run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "claw.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Plotting\n",
    "Now we'll plot the results, which are contained in a list called `claw.frames`.  It's simple to plot a single frame with matplotlib:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "figs = []\n",
    "for frame in claw.frames:\n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, frame.q[0, :], 'o-', linewidth=2)\n",
    "    axes.set_xlim((0, 1))\n",
    "    axes.set_ylim((-0.2, 1.2))\n",
    "    plt.close(figs[-1])\n",
    "\n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Waves\n",
    "\n",
    "In the solution, 3 waves are visible:\n",
    "1. A **shock wave** moving rapidly to the right as the low-density fluid is compressed.\n",
    "2. A **rarefaction** wave moving to the left as the high-density fluid expands.\n",
    "3. A **contact discontinuity** moving more slowly to the right.  This discontinuity in the density separates the region containing fluid that started in the high-pressure region and fluid that started in the low-pressure region.\n",
    "\n",
    "In fact, the solution of any Riemann problem consists of some combination of these three types of waves.  In the Euler equations, one of the waves is always a contact discontinuity, but each of the other two waves may be a shock or a rarefaction, depending on the left and right states."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Putting it all together\n",
    "\n",
    "For convenience, all of the code from the cells above to set up and run the shocktube problem is pasted together below.  Play around with the code.  You might:\n",
    "- Increase the number of grid points to see what the solution converges to.  Notice that the code still runs pretty fast even for larger grids.  This is because the bottom layer of code in PyClaw is compiled Fortran, not Python.\n",
    "- Change the initial left and right states, or set up a completely different initial condition.  See if you can generate a solution with two shock waves, or two rarefaction waves (some physical intuition is helpful here).\n",
    "- Change the ratio of specific heats\n",
    "- Make the boundaries periodic, so that there is a second shock wave moving left from $x=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from clawpack import pyclaw\n",
    "from clawpack import riemann\n",
    "\n",
    "claw = pyclaw.Controller()\n",
    "claw.tfinal = 0.2\n",
    "\n",
    "claw.keep_copy = True       # Keep solution data in memory for plotting\n",
    "claw.output_format = None   # Don't write solution data to file\n",
    "claw.num_output_times = 50  # Write 50 output frames\n",
    "\n",
    "riemann_solver = riemann.euler_with_efix_1D\n",
    "claw.solver = pyclaw.ClawSolver1D(riemann_solver)\n",
    "\n",
    "claw.solver.all_bcs = pyclaw.BC.extrap\n",
    "\n",
    "domain = pyclaw.Domain( (0.,), (1.,), (100,))\n",
    "x=domain.grid.x.centers # grid cell centers\n",
    "\n",
    "claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)\n",
    "\n",
    "gam = 1.4 # ratio of specific heats\n",
    "claw.solution.problem_data['gamma'] = gam\n",
    "claw.solution.problem_data['gamma1'] = gam-1.0\n",
    "\n",
    "rho_left = 1.0; rho_right = 0.125\n",
    "p_left = 1.0; p_right = 0.1\n",
    "\n",
    "claw.solution.q[0,:] = (x<0.5)*rho_left + (x>=0.5)*rho_right\n",
    "claw.solution.q[1,:] = 0.\n",
    "claw.solution.q[2,:] = ((x<0.5)*p_left + (x>=0.5)*p_right)/(gam-1.0)\n",
    "\n",
    "status = claw.run()\n",
    "\n",
    "figs = []\n",
    "for frame in claw.frames:\n",
    "    figs.append(plt.figure())\n",
    "    axes = figs[-1].add_subplot(1, 1, 1)\n",
    "    axes.plot(x, frame.q[0, :], 'o-', linewidth=2)\n",
    "    axes.set_xlim((0, 1))\n",
    "    axes.set_ylim((-0.2, 1.2))\n",
    "    plt.close(figs[-1])\n",
    "    \n",
    "# Animate the solution\n",
    "images = animation.make_images(figs)\n",
    "animation.JSAnimate_images(images)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.4"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}