write_a_hydrocode.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# How To Write A Hydro Code\n",
    "\n",
    "Michael Zingale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are _many_ methods for solving the equations of hydrodynamics.  We will make some choices right from the start:\n",
    "\n",
    "  * We will consider **finite-volume methods**.  These are popular in astrophysics because they are based on the integral form of the conservative equations and properly conserve mass, momentum, and energy.\n",
    "  \n",
    "  * We will consider an **Eulerian** grid: the grid is fixed and the fluid moves through it.\n",
    "  \n",
    "  * We will be **explicit in time**: the new solution depends only on the previous state.\n",
    "  \n",
    "  * We will look at a simple 2nd order **method-of-lines** integration.  We do this for simplicity here, and will point out where things are commonly done differently.  This scheme has a much simpler spatial reconstruction than methods that do characteristic tracing and relies on an integrator (like a Runge-Kutta method) to advance in time.\n",
    "  \n",
    "  * We will work in 1-d.\n",
    "  \n",
    "  * We won't cover in detail how to write a Riemann solver (that's a math exercise as much as anything else and beyond the scope of this notebook).\n",
    "  \n",
    "  * We'll assume a gamma-law equation of state&mdash;this is often not the case in astrophysics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "    \n",
    "Much more in-depth details and derivations are given in my hydro notes available online: https://github.com/Open-Astrophysics-Bookshelf/numerical_exercises\n",
    "\n",
    "For a greater variety of methods, in 2-d, see the pyro code: https://github.com/python-hydro/pyro2 (ref: [Harpole et al. JOSS](http://joss.theoj.org/papers/10.21105/joss.01265))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "We'll focus on the Euler equations.  In 1-d, these are:\n",
    "\n",
    "\\begin{align*}\n",
    "  \\frac{\\partial \\rho}{\\partial t} + \\frac{\\partial (\\rho u)}{\\partial x} & = 0 \\\\\n",
    "  \\frac{\\partial (\\rho u)}{\\partial t} + \\frac{\\partial (\\rho u^2 + p)}{\\partial x} &= 0 \\\\\n",
    "  \\frac{\\partial (\\rho E)}{\\partial t} + \\frac{\\partial (u(\\rho E + p))}{\\partial x} &= 0 \\\\\n",
    "  \\end{align*}\n",
    "\n",
    "This is a set of (hyperbolic) partial differential equations.  To close the system, we need an equation of state, relating the specific internal energy, $e$, to the pressure:\n",
    "\\begin{align*}\n",
    "e &= E - \\frac{1}{2}u^2 \\\\\n",
    "p &= \\rho e (\\gamma - 1)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To solve these, we need to discretize the equations in both space and time.  We'll use grid-based methods (in addition to the finite-volume method we'll consider, this can include finite-difference and finite-element methods).  \n",
    "\n",
    "Our system of equations can be expressed in conservative form:\n",
    "$$ \\frac{\\partial U}{\\partial t} + \\frac{\\partial F(U)}{\\partial x} = 0$$\n",
    "where $U = (\\rho, \\rho u, \\rho E)^\\intercal$ and\n",
    "$$\n",
    "F(U) = \\left ( \\begin{array}{c} \\rho u \\\\ \\rho u^2 + p \\\\ u (\\rho E + p) \\end{array} \\right )$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In a finite-volume method, we store the state of the fluid in discrete volumes in space, and we can refer to this discretized state with an index.  To see this, we integrate the conservative law system in space over a volume $[x_{i-1/2},x_{i+1/2}]$:\n",
    "$$\\frac{\\partial \\langle U\\rangle_i}{\\partial t} = - \\frac{F_{i+1/2} - F_{i-1/2}}{\\Delta x}$$\n",
    "\n",
    "This is the form of the equations we will solve.  Here, $\\langle U\\rangle_i$ represents the average state of the fluid in a volume:\n",
    "$$\\langle U\\rangle_i = \\frac{1}{\\Delta x} \\int_{x_{i-1/2}}^{x_{i+1/2}} U(x) dx$$\n",
    "Visually, we usually think of this grid as:\n",
    "\n",
    "![FV grid](fv_grid.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The state on the grid represents an instance in time.  We evolve the state by computing the fluxes through the volumes.  These fluxes tell us how much the state changes in each volume over some small timestep, $\\Delta t$.  \n",
    "\n",
    "Our code will have the following structure:\n",
    "\n",
    "  * Create our numerical grid\n",
    "  \n",
    "  * Set the initial conditions\n",
    "  \n",
    "  * Main timestep evolution loop\n",
    "  \n",
    "    * Compute the timestep\n",
    "    \n",
    "    * Loop to advance one step (count depends on the number of stages in the integrator)\n",
    "  \n",
    "        * Reconstruct the state to interfaces\n",
    "    \n",
    "        * Solve Riemann problem to find the fluxes through the interface\n",
    "    \n",
    "        * Do a conservative update of the state to the stage\n",
    "    \n",
    "    * Output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Grid\n",
    "\n",
    "We'll manage our 1-d grid via a class `FVGrid`.  We will divide the domain into a number of zones (or volumes) that will store the state.  To implement boundary conditions, we traditionally use ghost cells--extra cells added to each end of the domain.  We'll consider a grid that looks like this:\n",
    "\n",
    "![grid w ghostcells](simplegrid_gc.png)\n",
    "\n",
    "We'll use the names `lo` and `hi` to refer to the first and last zone in our domain.  The domain boundaries are the bold lines shown above, and beyond that, on each end, we have ghost cells.\n",
    "\n",
    "The main information we need to setup the grid are the number of zones in the interior and the number of ghost cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To make life easier, we'll have a simple class with indices that we use to index the fluid state arrays.  We can pass this around and be sure that we are always accessing the correct fluid state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "class FluidVars:\n",
    "    \"\"\"A simple container that holds the integer indicies we will use to\n",
    "    refer to the different fluid components\"\"\"\n",
    "    def __init__(self, gamma=1.4, C=0.8):\n",
    "        self.nvar = 3\n",
    "    \n",
    "        # conserved variables\n",
    "        self.urho = 0\n",
    "        self.umx = 1\n",
    "        self.uener = 2\n",
    "    \n",
    "        # primitive variables\n",
    "        self.qrho = 0\n",
    "        self.qu = 1\n",
    "        self.qp = 2\n",
    "    \n",
    "        # EOS gamma\n",
    "        self.gamma = gamma\n",
    "        \n",
    "        # CFL number\n",
    "        self.C = C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This is the main class for managing the finite-volume grid.  In addition to holding coordinate information and knowing the bounds of the domain, it also can fill the ghost cells and give you a scratch array that lives on the same grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "class FVGrid:\n",
    "    \"\"\"The main finite-volume grid class for holding our fluid state.\"\"\"\n",
    "    \n",
    "    def __init__(self, nx, ng, xmin=0.0, xmax=1.0):\n",
    "\n",
    "        self.xmin = xmin\n",
    "        self.xmax = xmax\n",
    "        self.ng = ng\n",
    "        self.nx = nx\n",
    "\n",
    "        self.lo = ng\n",
    "        self.hi = ng+nx-1\n",
    "\n",
    "        # physical coords -- cell-centered\n",
    "        self.dx = (xmax - xmin)/(nx)\n",
    "        self.x = xmin + (np.arange(nx+2*ng)-ng+0.5)*self.dx\n",
    "\n",
    "    def scratch_array(self, nc=1):\n",
    "        \"\"\" return a scratch array dimensioned for our grid \"\"\"\n",
    "        return np.squeeze(np.zeros((self.nx+2*self.ng, nc), dtype=np.float64))\n",
    "\n",
    "    def fill_BCs(self, atmp):\n",
    "        \"\"\" fill all ghost cells with zero-gradient boundary conditions \"\"\"\n",
    "        if atmp.ndim == 2:\n",
    "            for n in range(atmp.shape[-1]):\n",
    "                atmp[0:self.lo, n] = atmp[self.lo, n]\n",
    "                atmp[self.hi+1:, n] = atmp[self.hi, n]            \n",
    "        else:\n",
    "            atmp[0:self.lo] = atmp[self.lo]\n",
    "            atmp[self.hi+1:] = atmp[self.hi]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Reconstruction\n",
    "\n",
    "We need to use the cell-averages to figure out what the fluid state is on the interfaces.  We'll _reconstruct_ the cell-averages as piecewise lines that give us the same average in the zone.  We then follow these lines to the interfaces to define the left and right state at each interface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Usually we work in terms of the **primitive variables**, $q = (\\rho, u, p)$.  So we first write a routine to do the algebraic transformation from conservative to primitive variables:\n",
    "\\begin{align}\n",
    "\\rho &= \\rho \\\\\n",
    "u &= \\frac{(\\rho u)}{\\rho} \\\\\n",
    "p &= \\left ( (\\rho E) - \\frac{1}{2} \\frac{(\\rho u)^2}{\\rho}\\right )(\\gamma - 1)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def cons_to_prim(grid, U):\n",
    "    \"\"\"take a conservative state U and return the corresponding primitive\n",
    "    variable state as a new array.\"\"\"\n",
    "    v = FluidVars()\n",
    "    q = grid.scratch_array(nc=v.nvar)\n",
    "\n",
    "    q[:, v.qrho] = U[:, v.urho]\n",
    "    q[:, v.qu] = U[:, v.umx]/U[:, v.urho]\n",
    "    rhoe = U[:, v.uener] - 0.5*q[:, v.qrho]*q[:, v.qu]**2\n",
    "    q[:, v.qp] = rhoe*(v.gamma - 1.0)\n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Next we need a routine to create the interface states.  Here's well construct a slope for each zone, $\\Delta q$ based on the average state in the neighboring zones.  This gives us a line representing the value of the fluid state as a function of position in each zone:\n",
    "$$q_i(x) = \\langle q\\rangle_i + \\frac{\\Delta q_i}{\\Delta x} (x - x_i)$$\n",
    "\n",
    "Note that there is a unique $q_i(x)$ for each zone&mdash;this is usually called _piecewise linear reconstruction_.  By design, the average of $q_i(x)$ over the zone is the cell-average, so it is conservative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We use this equation for a line to find the fluid state right at the interface.  For zone $i$, the line $q_i(x)$ gives you the right state on the left interface, $q_{i-1/2,R}$, and the left state on the right interface, $q_{i+1/2,L}$.  Visually this looks like:\n",
    "\n",
    "![finding interface states](riemann-mol.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "There's one additional wrinkle&mdash;2nd order codes tend to produce oscillations near discontinuities, so we usually need to _limit_ the slopes, $\\Delta q_i$, so we don't introduce new minima or maxima in the evolution.  We'll use the minmod limiter:\n",
    "\\begin{equation}                                                                                                         \n",
    "\\left . \\frac{\\partial a}{\\partial x} \\right |_i = \\mathtt{minmod} \\left (                                               \n",
    "  \\frac{a_i - a_{i-1}}{\\Delta x}, \\frac{a_{i+1} - a_i}{\\Delta x} \\right )                                                \n",
    "\\end{equation}\n",
    "with \n",
    "\\begin{equation}                                                                                                         \n",
    "\\mathtt{minmod}(a,b) = \\left \\{                                                                                          \n",
    "    \\begin{array}{ll}                                                                                                    \n",
    "    a & \\mathit{if~} |a| < |b| \\mathrm{~and~} a\\cdot b > 0 \\\\                                                            \n",
    "    b & \\mathit{if~} |b| < |a| \\mathrm{~and~} a\\cdot b > 0 \\\\                                                            \n",
    "    0 & \\mathit{otherwise}                                                                                               \n",
    "    \\end{array}                                                                                                          \n",
    "  \\right .                                                                                                               \n",
    "\\end{equation}  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def states(grid, U):\n",
    "    v = FluidVars()\n",
    "    q = cons_to_prim(grid, U)\n",
    "\n",
    "    # construct the slopes\n",
    "    dq = grid.scratch_array(nc=v.nvar)\n",
    "\n",
    "    for n in range(v.nvar):        \n",
    "        dl = grid.scratch_array()\n",
    "        dr = grid.scratch_array()\n",
    "\n",
    "        dl[grid.lo-1:grid.hi+2] = q[grid.lo:grid.hi+3,n] - q[grid.lo-1:grid.hi+2,n]\n",
    "        dr[grid.lo-1:grid.hi+2] = q[grid.lo-1:grid.hi+2,n] - q[grid.lo-2:grid.hi+1,n]\n",
    "\n",
    "        # these where's do a minmod()\n",
    "        d1 = np.where(np.fabs(dl) < np.fabs(dr), dl, dr)\n",
    "        dq[:, n] = np.where(dl*dr > 0.0, d1, 0.0)\n",
    "\n",
    "    # now make the states\n",
    "    q_l = grid.scratch_array(nc=v.nvar)\n",
    "    q_l[grid.lo:grid.hi+2, :] = q[grid.lo-1:grid.hi+1, :] + 0.5*dq[grid.lo-1:grid.hi+1, :]\n",
    "\n",
    "    q_r = grid.scratch_array(nc=v.nvar)\n",
    "    q_r[grid.lo:grid.hi+2, :] = q[grid.lo:grid.hi+2, :] - 0.5*dq[grid.lo:grid.hi+2, :]\n",
    "    \n",
    "    return q_l, q_r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Riemann problem and conservative update\n",
    "\n",
    "After doing our reconstruction, we are left with a left and right state on an interface.  To find the unique fluid state on the interface, we solve a _Riemann problem_, \n",
    "$$q_{i+1/2} = \\mathcal{R}(q_{i+1/2,L},q_{i+1/2,R})$$\n",
    "\n",
    "We could spend an entire day talking about how to solve the Riemann problem.  Well just summarize things here.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "At each interface, we have a left and right state.  Information about the jump across this interface will be carried away from the interface by the 3 hydrodynamic waves ($u$ and $u\\pm c$).\n",
    "\n",
    "![Riemann solution structure](riemann-waves.png)\n",
    "The solution to the Riemann problem that we need is the state on the interface--with that we can evaluate the flux through the interface.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To solve the Riemann problem, we need to know how much each variable changes across each of the three waves.  To complicate matters, the left and right waves can be either shocks or rarefactions.  The middle wave ($u$) is always a contact discontinuity (and of our primitive variables, only $\\rho$ jumps across it).\n",
    "\n",
    "For a gamma-law gas, we can write down analytic expressions for the change in the primitive variables across both a rarefaction and shock.  We can then solve these to find the state inbetween the left and right waves (the star state) and then compute the wave speeds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Finally, we can find the solution on the interface by determining which region we are in.\n",
    "![Riemann state](riemann_state.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We'll use an exact Riemann solver to find the solution on the interface.  There a lot of algebra involved in finding the expressions for the jumps across the waves and the wave speeds, which we'll skip (by see my notes).  Instead we'll just use this solver to give us the state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One we have the interface state, we can compute the fluxes using this state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cons_flux(state, v):\n",
    "    \"\"\" given an interface state, return the conservative flux\"\"\"\n",
    "    flux = np.zeros((v.nvar), dtype=np.float64)\n",
    "\n",
    "    flux[v.urho] = state.rho * state.u\n",
    "    flux[v.umx] = flux[v.urho] * state.u + state.p\n",
    "    flux[v.uener] = (0.5 * state.rho * state.u**2 +\n",
    "                     state.p/(v.gamma - 1.0) + state.p) * state.u\n",
    "    return flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on module riemann_exact:\n",
      "\n",
      "NAME\n",
      "    riemann_exact\n",
      "\n",
      "DESCRIPTION\n",
      "    An exact Riemann solver for the Euler equations with a gamma-law\n",
      "    gas.  The left and right states are stored as State objects.  We then\n",
      "    create a RiemannProblem object with the left and right state:\n",
      "    \n",
      "    > rp = RiemannProblem(left_state, right_state)\n",
      "    \n",
      "    Next we solve for the star state:\n",
      "    \n",
      "    > rp.find_star_state()\n",
      "    \n",
      "    Finally, we sample the solution to find the interface state, which\n",
      "    is returned as a State object:\n",
      "    \n",
      "    > q_int = rp.sample_solution()\n",
      "\n",
      "CLASSES\n",
      "    builtins.object\n",
      "        RiemannProblem\n",
      "        State\n",
      "    \n",
      "    class RiemannProblem(builtins.object)\n",
      "     |  RiemannProblem(left_state, right_state, gamma=1.4)\n",
      "     |  \n",
      "     |  a class to define a Riemann problem.  It takes a left\n",
      "     |  and right state.  Note: we assume a constant gamma\n",
      "     |  \n",
      "     |  Methods defined here:\n",
      "     |  \n",
      "     |  __init__(self, left_state, right_state, gamma=1.4)\n",
      "     |      Initialize self.  See help(type(self)) for accurate signature.\n",
      "     |  \n",
      "     |  find_star_state(self, p_min=0.001, p_max=1000.0)\n",
      "     |      root find the Hugoniot curve to find ustar, pstar\n",
      "     |  \n",
      "     |  rarefaction_solution(self, sgn, state)\n",
      "     |      return the interface solution considering a rarefaction wave\n",
      "     |  \n",
      "     |  sample_solution(self)\n",
      "     |      given the star state (ustar, pstar), find the state on the interface\n",
      "     |  \n",
      "     |  shock_solution(self, sgn, state)\n",
      "     |      return the interface solution considering a shock\n",
      "     |  \n",
      "     |  u_hugoniot(self, p, side)\n",
      "     |      define the Hugoniot curve, u(p).\n",
      "     |  \n",
      "     |  ----------------------------------------------------------------------\n",
      "     |  Data descriptors defined here:\n",
      "     |  \n",
      "     |  __dict__\n",
      "     |      dictionary for instance variables (if defined)\n",
      "     |  \n",
      "     |  __weakref__\n",
      "     |      list of weak references to the object (if defined)\n",
      "    \n",
      "    class State(builtins.object)\n",
      "     |  State(p=1.0, u=0.0, rho=1.0)\n",
      "     |  \n",
      "     |  a simple object to hold a primitive variable state\n",
      "     |  \n",
      "     |  Methods defined here:\n",
      "     |  \n",
      "     |  __init__(self, p=1.0, u=0.0, rho=1.0)\n",
      "     |      Initialize self.  See help(type(self)) for accurate signature.\n",
      "     |  \n",
      "     |  __str__(self)\n",
      "     |      Return str(self).\n",
      "     |  \n",
      "     |  ----------------------------------------------------------------------\n",
      "     |  Data descriptors defined here:\n",
      "     |  \n",
      "     |  __dict__\n",
      "     |      dictionary for instance variables (if defined)\n",
      "     |  \n",
      "     |  __weakref__\n",
      "     |      list of weak references to the object (if defined)\n",
      "\n",
      "FUNCTIONS\n",
      "    cons_flux(state, v)\n",
      "        given an interface state, return the conservative flux\n",
      "\n",
      "FILE\n",
      "    /home/zingale/classes/how_to_write_a_hydro_code/riemann_exact.py\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import riemann_exact as re\n",
    "help(re)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For a method-of-lines approach, we want to just compute the righthand side, $A = -\\partial F/\\partial x$.  Then we will turn our PDE into an ODE for time:\n",
    "$$\\frac{\\partial \\langle U\\rangle_i}{\\partial t} = -A_i = - \\frac{F_{i+1/2} - F_{i-1/2}}{\\Delta x}$$\n",
    "\n",
    "We can then use any ODE integration method, like Runge-Kutta to solve the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This routine will take the conserved state, $U$, construct the left and right states at all interfaces, solve the Riemann problem to get the unique state on the boundary, and then compute the advective term and return it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def make_flux_divergence(grid, U):\n",
    "    \n",
    "    v = FluidVars()\n",
    "    \n",
    "    # get the states\n",
    "    q_l, q_r = states(grid, U)\n",
    "\n",
    "    # now solve the Riemann problem\n",
    "    flux = grid.scratch_array(nc=v.nvar)\n",
    "    for i in range(grid.lo, grid.hi+2):\n",
    "        sl = re.State(rho=q_l[i,v.qrho], u=q_l[i,v.qu], p=q_l[i,v.qp])\n",
    "        sr = re.State(rho=q_r[i,v.qrho], u=q_r[i,v.qu], p=q_r[i,v.qp])\n",
    "        rp = re.RiemannProblem(sl, sr, gamma=v.gamma)\n",
    "        rp.find_star_state()\n",
    "        q_int = rp.sample_solution()\n",
    "        flux[i, :] = cons_flux(q_int, v)\n",
    "\n",
    "    A = grid.scratch_array(nc=v.nvar)\n",
    "    for n in range(v.nvar):\n",
    "        A[grid.lo:grid.hi+1, n] = (flux[grid.lo:grid.hi+1, n] -\n",
    "                                   flux[grid.lo+1:grid.hi+2, n])/grid.dx\n",
    "\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Timestep\n",
    "\n",
    "Explicit hydro codes have a restriction on the size of the timestep.  We cannot allow information to move more than one zone per step.  For the hydro equations, the speeds at which information travels are $u$ and $u \\pm c$, so we use the largest speed here to compute the timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def timestep(grid, U):\n",
    "\n",
    "    v = FluidVars()\n",
    "    \n",
    "    # compute the sound speed\n",
    "    q = cons_to_prim(grid, U)\n",
    "    c = grid.scratch_array()\n",
    "    c[grid.lo:grid.hi+1] = np.sqrt(v.gamma *\n",
    "                                   q[grid.lo:grid.hi+1,v.qp] /\n",
    "                                   q[grid.lo:grid.hi+1,v.qrho])\n",
    "\n",
    "    dt = v.C * grid.dx / (np.abs(q[grid.lo:grid.hi+1, v.qu]) +\n",
    "                          c[grid.lo:grid.hi+1]).max()\n",
    "    return dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Main driver\n",
    "\n",
    "This is the main driver.  For simplicity, I've hardcoded the initial conditions here for the standard Sod problem.  Usually those would be a separate routine.\n",
    "\n",
    "This does 2nd-order RK (or Euler's method) for the integration, and requires that we compute the advection terms twice to advance the solution by $\\Delta t$.  The update looks like:\n",
    "\\begin{align*}\n",
    "U^\\star &= U^n + \\frac{\\Delta t}{2} A(U^n) \\\\\n",
    "U^{n+1} &= U^n + \\Delta t A(U^\\star)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mol_solve(nx, tmax=1.0, init_cond=None):\n",
    "    \"\"\"Perform 2nd order MOL integration of the Euler equations.\n",
    "    You need to pass in a function foo(grid) that returns the \n",
    "    initial conserved fluid state.\"\"\"\n",
    "\n",
    "    grid = FVGrid(nx, 2)\n",
    "    v = FluidVars()\n",
    "    \n",
    "    U = init_cond(grid)\n",
    "    \n",
    "    t = 0.0\n",
    "    \n",
    "    while t < tmax:\n",
    "        dt = timestep(grid, U)\n",
    "        if t + dt > tmax:\n",
    "            dt = tmax - t\n",
    "\n",
    "        grid.fill_BCs(U)\n",
    "        k1 = make_flux_divergence(grid, U)\n",
    "\n",
    "        U_tmp = grid.scratch_array(nc=v.nvar)\n",
    "        for n in range(v.nvar):\n",
    "            U_tmp[:, n] = U[:, n] + 0.5 * dt * k1[:, n]\n",
    "\n",
    "        grid.fill_BCs(U_tmp)\n",
    "        k2 = make_flux_divergence(grid, U_tmp)\n",
    "\n",
    "        for n in range(v.nvar):\n",
    "            U[:, n] += dt * k2[:, n]\n",
    "\n",
    "        t += dt\n",
    "\n",
    "    return grid, U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example: Sod's problem\n",
    "\n",
    "The Sod problem is a standard test problem, consisting of a left and right state separated by an initial discontinuity.  As time evolves, a rightward moving shock and contact and leftward moving rarefaction form.\n",
    "\n",
    "One reason this problem is so popular is that you can find the exact solution (it's just the Riemann problem) and compare the performance of your code to the exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def sod(grid):\n",
    "    \n",
    "    v = FluidVars()\n",
    "    U = grid.scratch_array(nc=v.nvar)\n",
    "    \n",
    "    # setup initial conditions -- this is Sod's problem\n",
    "    rho_l = 1.0\n",
    "    u_l = 0.0\n",
    "    p_l = 1.0\n",
    "    rho_r = 0.125\n",
    "    u_r = 0.0\n",
    "    p_r = 0.1\n",
    "\n",
    "    idx_l = grid.x < 0.5\n",
    "    idx_r = grid.x >= 0.5\n",
    "\n",
    "    U[idx_l, v.urho] = rho_l\n",
    "    U[idx_l, v.umx] =  rho_l * u_l\n",
    "    U[idx_l, v.uener] = p_l/(v.gamma - 1.0) + 0.5 * rho_l * u_l**2\n",
    "\n",
    "    U[idx_r, v.urho] = rho_r\n",
    "    U[idx_r, v.umx] =  rho_r * u_r\n",
    "    U[idx_r, v.uener] = p_r/(v.gamma - 1.0) + 0.5 * rho_r * u_r**2\n",
    "    \n",
    "    return U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "g, U = mol_solve(128, tmax=0.2, init_cond=sod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.rcParams['figure.figsize'] = [8, 6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sod = np.genfromtxt(\"sod-exact.out\", skip_header=2, names=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f27746dec50>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD+CAYAAAAwAx7XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8lOWd9/HPb2ZyAhLOFeQUsJxCt0WNUsBaTzxFCNKqdbVKS1pX3a12wW5duwpdo7XauvVUty2PNti11bpon8qhVamth6DWoFglUUSICIZDOAaSkMNczx+TGe6EkAxhksnMfN+v17zuzNzXPXPdCflx5XedzDmHiIgkH1+8KyAiIl1DAV5EJEkpwIuIJCkFeBGRJKUALyKSpBTgRUSSlAK8iEiSUoAXEUlSCvAiIkkqEM8PHzRokMvNzY1nFUREEs7atWurnHODOyoX1wCfm5tLaWlpPKsgIpJwzOyjaMopRSMikqQU4EVEkpQCvIhIklKAFxFJUgrwIiJJSgFeRCRJKcCLiCSphA7w3u0GnXORB0AwGGxx9J5vXaa9c4l2fVvfGxFJTXGd6HQi7n1+AwfqGlhckMd9qz/gQG0Db27Zi5lx9thBrC7fQZ+MAAcPN3L+hE/x8sYqcOCgRZn2ziXa9ZNH9iMnM42FM8Zx2/L19M1KZ8EFYzEznHOYWRx/YiLS3Syalp6Z/QooAHY65z5zjDIzgfsBP/Cwc+6ujt43Pz/fdWYmq3OOohVlFJdUMH/aKAyjeE1F5PyAXmnsqWkgI+DjcGOQ/lkB9tY2tniPcJn2ziXi9ROHZpPu9/H21v18Y+pIDCMnK43qw43kZKZFAr6IJC4zW+ucy++wXJQB/mzgIPDrtgK8mfmBDcAMYCvwBnCFc66svfftbICHlkG+tXn+5+jLoU69bzLpneHn0OEmBmens6u6nlNH9sUwMtL8TD9lEA6HYdA/Fz771XhXV0SiFG2AjypF45x7ycxy2ylyJrDRObep+cOfAOYCRwV4M7sGuAZg5MiR0Xx8m8yMxQV5bQb4b/r/yGjfjk6/d9IIAmlAXfOx0nNuC7Roxw8aCydP7sbKiUhXi1UOfhjwsef5VmBKWwWdc0uAJRBqwXf2A8Mt+LZcUH9PZ9826U0cks0fvj2NH/7xPR5d8xHXTBnMTe9dhr18D/zjY/GunojEUEKOommdgy+cltvifN9emTThJxBIowk/OVkZNOFv8QiXae9cIl4/bmi/o8p6H+9ur+GURav51ZqtfH36KXxv7hReGnAJlC+HneXx+YGKSJeIVQt+GzDC83x482tdwszIyUyjcHpuZBRN4bTclB5F44KO+qCjvLKavKHZlFVWM7B3GrsPNRzz+3jrrAncvrKc32/+PH/r9QRpL92DXfpIV/3YRKSbRdXJCtCcg19xjE7WAKFO1vMJBfY3gK8559a3954n0skKtBj6570PMyMYDOLz+SLH1vfpLdPeuUS6/r7VGzhQ10hOZoD9NQ1gsHTNR0we3pfPjejLyne2U3Ww/qjvY+H0XBanP4579SGWnraMb150/jG/5yISf9F2skaVojGzx4FXgfFmttXMvtX8+iozO9k51whcDzwLlANPdhTcY8E73M/MIg8An8/X4ug937pMe+cS6fqFM8azuCCPhTPG84OLJtE3K53Cabk8/S/T8JmPqoP15A3N5vpzxrT4Pt46awI/PTSDoIPx257SJCmRJBHtKJorjvH6LM/Xq4BVMaqXdFI4+IcC/rjIXzk5WaGU1q2zJnDHqvdaXHPKLX8CYNaQs5hW/SzW1IDzp2m8vEiCS8hOVoleOEgvnDGORbMncseq9yguqaBwei4f/nBmi7LjZn4bq6nCvb+KohVl3Pv8hnhUWURiRAE+hfh8vkjndDjYe839UwYuZxgfPvsQxSUVHKhrULpGJIEpwKeYcEv+9pXlLVryeUOzeXd7DffvmcKY/W+wID+DxQV5StOIJDAF+BTkbckvLsjD7/ez4oazAHiy8RwA/rX/q5FFykQkMSXsapJyYrwdsM45bl8ZmuT0CYN4MfhZTn9tKX2++O/c/scPIitUikhiUQs+hYWDe3hWcN7QbABeyZlNTkMVRfc+oFy8SAJTgE9x3lnBK244i8LpuTy6ewK7XF+mH1gZSeMoFy+SeJSikRbpmvAKnf/b9EWu9S/ngrP7acMQkQSlFrwALdM1AE80nYvfHH/53X0Eg0GNixdJQArwAnBULn6LO4l1gc8ybtvTXPTAi8rFiyQgBXgB2s7FL6k5h+FWxeCdJcrFiyQg5eAlonUufmxJPjtdP67yr+a8gluUixdJMGrBSwveXHwjAR5vOpdzfet44KnVysWLJBgFeGmhdS7+8cbzcGakrfs1BQ++oly8SAJRgJcWWufiL5x+OqubTuMf/X9hY+Ue5eJFEohy8HKU1rn4ea/O4Ev+Umb6/sbigrkK7iIJQi14aZM3F18SnMTm4EnMCzxP0YqySHpGaRqRnk0BXtrkzcXPnz6GLWMu5wzfBl5d8xJFK8rU4SqSABTgpU3eXPyi2RN5re+F1Lk0bsh+keyMQGQ9eXW4ivRcFs9fzvz8fFdaWhq3z5eOeZcUXvezrzG26gU+f/hnHKSXOlxF4sTM1jrn8jsqpxa8tMu7iffkr3yXPlbHxf6XARTcRXo4BXiJinOOojczWRccwzf8z2EEIx2uStGI9EwaJikdinS4rvmI7MEXc2P1Pdw+aQe3lvhCKRyMnCzt+iTS06gFLx0Kd7jOnzaKQ6fMYafrx5Rdy5g/bRTrPt5P8Rp1tor0RGrBS1TCk58A/rrjUs6tfJiXXn2VTe5kdbaK9FBqwUvUzAwz45yv3US98/N1/3OAOltFeioFeDkuzjmK/lrF8uBULvW/RDY16mwV6aGUopGoeWe3XjLkUi7Z9wo/+fTfua6klzpbRXogteAlat7O1pwxZ1IaHMeZO5dROHW4OltFeiC14OW4eDtbl+2dR/7mRWx5/Q+sC56uzlaRHkYteDlu4c7WS6+8jko3gEL/nwB1tor0NArw0inOOYr++AH/0ziDs/zrmWBbtJSwSA+jAC/HzdvZ6vILafBl8t3s5yguqdBSwiI9iAK8HDfvUsLfmzuFNwcW8MXDL/KFk+q1lLBIDxJVJ6uZzQTuB/zAw865u9oosxC4GnDAO0Chc64uhnWVHsS7rd+ZV9yKe+Bppu1+irtfSAdQh6tID9BhC97M/MBDwIVAHnCFmeW1KjMM+A6Q75z7DKH/CC6PfXWlJ4ksJTxgNJY3h6/5/0xvagF1uIr0BNGkaM4ENjrnNjnn6oEngLltlAsAWWYWAHoBn8SumtKTOed4JFhAX6vhMv9fATS7VaQHiCZFMwz42PN8KzDFW8A5t83M7gG2ALXAc86559p6MzO7BrgGYOTIkZ2ps/QgkQ7Xdb2Z1ieP76SvxiZdza9KKjS7VSTOYtLJamb9CbXqRwMnA73N7Kq2yjrnljjn8p1z+YMHD47Fx0sceWe3/n3kPPrXVzJx30taSlikB4imBb8NGOF5Prz5Na8LgM3OuV0AZvY0MA14LBaVlJ4tMrs1OIHdP36IsR8W8736MYCps1UkjqJpwb8BjDWz0WaWTqjz9JlWZbYAnzezXhb6TT4fKI9tVaUnMzPMH2DA+QuZ7NvEGfY+oM5WkXjqMMA75xqB64FnCQXtJ51z6wHMbJWZneycex1YBrxJaIikD1jSZbWWHsk5x52Vp7LbZXNdYDmgzlaReIpqHLxzbhWwqo3XZ3m+/gHwg9hVTRJJpLP11e2MHjCXr9U8xr9PauBudbaKxI1mskpMeDtbP/70VRx0mZy14zF1torEkZYLlpjxLiX86o6LmbL9cb69fRZb3EnqbBWJA7XgJabCSwlPvXIRjfi51r8CONLZqha8SPdRgJeYc85R9OJenmo6m0v9LzKYvVplUiQOFOAlprxLCa8ecDkBmvj3vn+muKSCggdf0SqTIt1IAV5iyruU8MMLvkrZgPOZWbeKHA5SVlmtXLxIN1KAl5hbOGMciwvy8Pl8fOayH9DH6vi6/3lAE59EupMCvHSJcIdqUamfF5omUxj4E1nUaVs/kW6kAC9dwpuL35z3zwy0am7s/4q29RPpRgrw0iW8ufjCf7yMTdn5fLnmKU4dkqZt/US6iSY6SZfxbus3+pIibOksTt31DA+8cCGgbf1Euppa8NKlItv65U7H5X6B6wLLyaAe0OQnka6mAC/dwjnHo+mX8ynbxxX+FwCUixfpYkrRSJeLdLj+vT+n95nEdW45bw2eS3FJBa9v2h0ZHx9O54hIbKgFL13O2+E66Yo7GWJ7+dyu0J4xmvwk0nXUgpduEelwBdzIqfzzR8t5ouk86klTcBfpImrBS7cxMxzwPxmXM9T2cJn/r4B2fRLpKmrBS7eJ5OLfGUR+nzxuSltO5qR5PKxdn0S6hFrw0m2O7PqUyxuj/5mchirO2PV77fok0kXUgpdudWTXp0ls+umj5G9dysIPT6WGTHW2isSYWvDS7cK7Po2+7C4GWjWF/j8BmvgkEmsK8BIXzjmK1vXi+abTuDawghwOauKTSIwpwEu38640+ft+88mxGv6j3/Pa9UkkxhTgpdt5Jz79bOE83ul/AXNqn2Eg+zXxSSSGFOAlLlrs+nTlXWTQwL8EQrNbvcFdrXiRzlOAl7iJ7Pr0aj1PNZ3NVf7VDGU3ty1fH5n4pHy8SOdpmKTEjTcX/6+n34C/rISFgWXctGYgAIZRvKZCC5GJdJICvMSNNxe/oCAPnruWr776EI80XcjSNaEyyseLdJ5SNBJX4Vy8mWFf+C5k5nBz4PHIeQV3kc5TgJe4i3SoZvXn+YHzONf/NlN96wEtRCZyIpSikR4hnI//7Yf5vNxrMP/V7yn+7/iZFGshMpFOUwteeoRwPv6KaWMpGXEdJ9e8x6S9L2ghMpETYPH8hcnPz3elpaVx+3zpeZxz4ILs+MkU6g/t44L6e6gnTZ2tIh5mttY5l99RObXgpUcxM8zn56RL7makbxdX+lcDWohMpDMU4KXHcc5RVD6El5r+ge8Efk9fLUQm0ilRBXgzm2lm75vZRjO7+Rhl+pnZMjN7z8zKzWxqbKsqqcA7+emJ/teQwyGK+i7XQmQindBhgDczP/AQcCGQB1xhZnltFL0f+JNzbgLwOaA8lhWV1NBiIbIFV/HW4C8zu24lp9g2LUQmcpyiacGfCWx0zm1yztUDTwBzvQXMrC9wNvAIgHOu3jm3L9aVldTgXYjs9Pk/oYYMbgn8BtBCZCLHI5oAPwz42PN8a/NrXqOBXUCxmb1lZg+bWe+23szMrjGzUjMr3bVrV6cqLckvshDZX3byQOPFnOdfxxd9b2shMpHjEKtO1gBwGvBz59ypwCGgzVy9c26Jcy7fOZc/ePDgGH28JBtvLt6d8U/szhjOrYHHeGzNh9y2fD1Fy8uUjxfpQDQzWbcBIzzPhze/5rUV2Oqce735+TKOEeBFouHNxd9akAcTfszAJ77Glf4/s3RN6J9tOB8vIm2LpgX/BjDWzEabWTpwOfCMt4BzbjvwsZmNb37pfKAspjWVlNNiIbLxs3BjzmFhYBn9qAaIBHelakTa1mGAd841AtcDzxIaGfOkc249gJmtMrOTm4veAPzGzP4OTAbu7JoqSyqJdKgCP8/4Jn2o5cbAMgClakQ6ENViY865VcCqNl6f5fl6HdDh1FmR4xXJx78V4LQRlzBv1zKebPrikTXjp2nopEhbNJNVejxvPn7KN38CvQdze9pSjCAAi+doGQORtijAS0II5+PJ7Mf/G3wtp/o28lX/iwBaxkDkGLQevCSUohVlFL83gbF98rjZPU7FoHMpLqng9U27IzNdtX+rSIha8JIwjqRqRpP3rV/Sz2qZs/tXAFrGQKQNasFLQlk4Y1ykhe6m/BNXvvZLftd0Du+6MS2WFFaQF1ELXhJQOIjfffhidpPDHWnF+AgqFy/SilrwknAiwyZfq8I34Gpuqvkv/q3/K/y4xKdcvIiHWvCScLzDJv/tu7ewMedM5tU8yhB2Kxcv4qEALwkpsqSw388p85cQoInb0h4FtKSwSJgCvCSsyJLCJTXc13gJX/KX8iXfG1pSWKSZcvCSsLxLCn9z6nVsf6+U29xSLlgzCQDDKF5ToXy8pCwFeElY3lz8ooI8OP2X8PD5fC/wO36wphDQksKS2pSikYTWYknh4afDlGuZ51/NaRZKy2hJYUllCvCS8Lwdqj+qu5hKBvDjtCVkUK8lhSWlKcBLUgjn45e8XsULY2/h075PWBB4iqVrPgrl4bWksKQgBXhJCt58/FVXfhN36jyu8a/gc7YR0JLCkpoU4CVpRJYUBu5289hBf36S9ksyqNcyBpKSFOAl6RStKOMXr1Xxi+wbGOfbRlG/lRSXVFDw4CvKxUtKUYCXpOJN1fznjQt4a+BsLq1dxmftw6OWMVCQl2Rn8fxHnp+f70pLS+P2+ZK8IksK1+5lx12nst/15qL6Oyi/cy4+ny/SKZuTmcbCGePiXV2R42Jma51zHe6BrRa8JKXIMgarP+HmhqsZ79vKdwP/S8GDr0Ry8UrXSLLTTFZJSt5lDAqnfwVHFVevXcpfdkxmzH9UA5rlKslPLXhJSt5c/OKCPOxLP8QGjOGetF+QwyFAs1wl+SnAS9LyLmPg0nrx8OCbOYm9FKUVA2iWqyQ9pWgkqUVy8SvKKH67N58eXciXKx/hz02nsXRNqEx4lqtIslELXpKeN11zztV344afwR1pv2Iou4HQLFdQqkaSjwK8pITILFdfgAdz/o0ATfxX2s/xEVSqRpKWUjSSUopWlFH8ZhOjxt3I3C0/4rrgM/z3mi8DStVI8lELXlKGN1Vz0fybcJMu5sbAsiNrx8/JO2qbv2AwGLk+vA1guIUfPhc+es+3LtPeueO93lsfkfZoJquknPC/+bt//xpXrpuH35q4ov5W+g0bT0NTkLLKauZPG0VOZhqry3dwwcSTMIMDtY28uWUvZsbZYwexunwHfTICHDzcyPkTPsXLG6vAgYMWZdo7d7zXTx7ZLzL79rbl6+mblc6CC8ZGOpO1HHJqiHYmqwK8pBzvJKibJx/m6x8s4GCjcVX999ngRjCoTzqz/2EIf6vYS3llNROG9MEwyrdXR95jQK809tQ0kBHwcbgxSP+sAHtrG1t8TrhMe+c6c/3Eodmk+328vXU/35g6Ep/5yM4MUH24UUsvpAgFeJF23Pv8Bg7UNYRy7rveY8dDF5JNDdvdgBbl0gM+6huDx3iX+AvXL3zsm5XGwD4ZHLMdn5YJl/0aBozpzmpKjEUb4NXJKilp4YxxkVRN0etBVtcv5nr/H+hldS0LNsShcsej4chx9KDejB6aw/rKA6T5fYw/KTuU7gmXrT8IHzwH295UgE8RCvCS0o6sVzOFr86eT8GDr1BWWd3xhT3QdyZ+mpvKd0SWRV40eyK3ryw/krY5UAk/nQB1++NdVekmGkUjKcs7qiYcDMsqq5k4pA+fG5bTouzEIdktng/olQZARiD0K9Q/6+i2UrhMe+c6c/3EodlHlQV44IWNlFVWkzc0O3I/Lcb2Z/ULFazb1+b1knyiasGb2UzgfsAPPOycu+sY5fxAKbDNOVcQs1qKdJFwqiYc7OdPGwXA0jUfMXl4Xw43BSMdmuHOzXiNonFBR33QUd4cxMsqqxnUJ52qg/Ut7qmsspox//FHoNWKmYFM8KerBZ9COuxkbQ7aG4AZwFbgDeAK51xZG2VvBPKBnGgCvDpZpadxznHf6g8iHbDh4F+0oozsjEBkhIqZEQwG8fl8kWPr3yVvmfbOHc/1963ewIG6RnIyAxyobeTW2ROY87OSY6aVNv9oFsCRzU3WzYbxF8JFD8Tk+yXxEctO1jOBjc65Tc1v/AQwF2gR4M1sODAb+CFw43HXWKQHMLMWrfrwuPLwqpRePp+vxbGtMejRnDue6xfOGB+pWzAYjKSVCqfncuusCRQ8+Arl2w9Grr9t+XoMo3hNBYXTc3FZ/TClaFJGNAF+GPCx5/lWYEob5e4DbgLaThA2M7NrgGsARo4cGV0tRbpZ62DbkyYQhevi8/mO6kMo336QiUP6kBHws27rfpau+Qg4sgyDPdJXKZoUEpNRNGZWAOx0zq01s3PaK+ucWwIsgVCKJhafL5KqWvchhIO9mTH6+6si5RbPaV4XP7MfdmhXHGss3SmaUTTbgBGe58ObX/OaDlxkZhXAE8B5ZvZYTGooIu06kr45ssFJ0YqWXWRFK8oIBoO8uxv27a2KRzUlDqIJ8G8AY81stJmlA5cDz3gLOOe+75wb7pzLbT7/gnPuqpjXVkQ6FB7bn9c8nDJvaDbFJRUUPPgKb+2C9IYDWqgsRXQY4J1zjcD1wLNAOfCkc249gJmtMrOTu7aKIhItb6pmxQ1nUTg9NzLCpqyymtzhJ5MVPIgpwKcErUUjkoTCeXnnXItc/OY5H2LPL4KbP4bMnHbeQXqyaIdJaiarSBLy7kXrtXxDTegLDZVMCQrwIknIuyRy4fRcNv9oFvOnjWLlB7Wh87X7lIdPAQrwIknIm4tfXJDHfas/wDCmTWpeRbJunzYZTwEK8CJJKrLROHCgroHiNRXU+UMja558+V1tMp4CFOBFklh4uYXFBXkUTs/lf94OzWItfX9zpHXfk2bpSmwpwIukgHCQP+B6A5DDIQX3FKAAL5ICwp2u1WQRdEaO1VC0okzpmSSnAC+S5LwjauZPH4Nl5TB1qJ/ikorIEgaSnLRln0iS846oyckMsC/YmzOGGIUjc8nOCLTc1k+SilrwIilg4YxxLJo9kQN1jWyrS+eDLdtYNHsi1YcbNZomiakFL5IifD4fiwvy2LxhIHv2Vh21rZ86XJOPWvAiKcTMGD18GDkcirym4J68FOBFUohzjrd2OvrakQCv0TTJSwFeJEWER9O8sSPIAF8Nm+68kMLpuRpNk8SUgxdJEeHRNGNHDiPtk3ruWL6ORXMmA2g0TZJSC14khSycMY5zPjcWgD+8Ws7tK8s1miaJqQUvkmJ8vQYA8I3T+3FPSQXFJRWARtMkI7XgRVJNZl8Avj1lUIuXw8FdLfjkoQAvkmoy+wHw25f+3uLlohVlNDU1tVgnXsE+sSnAi6QY17wX6+tlm8kbGlofPm9oNsUlFZzxwz9TXFLB/tp6gsEgty1fz73Pb8A51+IBREbdhI/tlYn2eu9RTpxy8CIpxrL6AzDzlExmFp7F7SvLI3n4PTUNDOqTDsBX/nsNb2/dzzemjuQrD5UA4AiNxjl77CBWl++gT0aAg4cbOX/Cp3h5YxW4o8u0d857/cH6Jvqk+zlY3xQZzRPePFw6RwFeJNU05+Av/HQm1rx8QTjAA1QdrGfpmo8iz/+2eQ/l2w+2eIstuw+xp6aBjICPw41BKvfVsre2sc0y7Z0LX799fx17ahoY0CuNPTUNzJ82imAwGBm6ueCCsQr0naAAL5JqAhkQyMLefQpX9QHvbNvHvWmHjl1+N5DW6rUGz2tpQGM7Zdo7F76+4cgxI8vH8M1ZvPsgfLbqEKcM6s07ZZAW8DFh5FBsRhE0p5mkfQrwIqlo0pdxW15j74YScmob+WKvAP17pVO5v476xjjPaG0i9J8KMCXNR+YhHwfqGumf4bBdO3m6djIXXzY/njVMGArwIqnoK7/AgEef38CBugYWzZ5IUXMufuKQPmQE/Kzbuj/etYT60GH+tFEMqtvC9WVXQM0e5eajpFE0Iils4YxxLC7Iw+fzkZOZxvxpo5gyZiDrtu6PjLAZ2Lt1fgUG9Aq9lhEIhZD+WUe3FcNl2jsXvj78PHxsbemaj3jkzdB/OF+ZkKXgHiW14EVSXDhYhket3Lf6g8juT1NGD2RRwUQu/vmaDkfBnMgomvPGD+aF93dRVll9zL8g9tMHh2G1e7v9e5SoLJ5jTvPz811paWncPl9E2hZOgXiPXmZGMBjE5/NFju2Vieb6e59/n/21DUCoxT5/6iher9hDeWV15Lry3teSedrl2Ox7UjpNY2ZrnXP5HZVTC15EjhIOnK2PXj6fr8WxvTLRXL9wxvjIXxDzp43CMMorqymclsubW/aybV8t2+t7cXDjZiZ5hlBq9ctjUw5eRHoMM2PhjHH8YM4kcrJCG4UvKpjIaSP7U3WwnsNp/cgOVkcmZ2n1y/apBS8iPU440IfTMIvn5IHBJ3/LYvDeHRSXVGj1yyioBS8iPZY3RbS4II+99KG/hWbVKrh3TAFeRHq88HaD+1w2/QgFeO0l2zEFeBHp0cLBvbikgtEjh9PH6vjW1GGRvWQV5I9NOXgR6dHCe8kWTs/l3JMmwCdw67knEfSlkZOZpjRNO6IO8GY2E7gf8AMPO+fuanV+BPBr4CRC8xmWOOfuj2FdRSRFRTpc178PgNXuVQ4+ClGlaMzMDzwEXAjkAVeYWV6rYo3Ad51zecDngW+3UUZEpFPMDJr3k6V2TyS4K0VzbNHm4M8ENjrnNjnn6oEngLneAs65Sufcm81fVwPlwLDWb2Rm15hZqZmV7tq168RqLyIp5bG3Q7NaXU1ouclwfj68xaC0FG2AHwZ87Hm+lTaCd5iZ5QKnAq+3PuecW+Kcy3fO5Q8ePDj6mopISnPOsSvYG4AVr61v0fmqCU9ti3knq5n1AZ4CFjjnDsT6/UUkNZkZC+ZMgXeg7MOPuOH7qwA04akd0bbgtwEjPM+HN7/WgpmlEQruv3HOPX3i1RMROcLSe+MCmfSzIwuQKbgfW7QB/g1grJmNNrN04HLgGW8BC32HHwHKnXM/jW01RURCaZpqy6Y/R/aI1Vj4Y4sqwDvnGoHrgWcJdZ4+6ZxbD2Bmq8zsZGA6MA84z8zWNT9mdVG9RSTFhHPu2w5n8dmBTWz+0SwKp+dqwlM7os7BO+dWAavaeD0cxD8B9HeSiHSJ8ISnzL6DyM1pjKxPA2jC0zFoJquIJIyFM8bh9o7AdpYDRxYhU3Bvm9aiEZGEYlkDoGbPkecK7sekAC8iiaXXAKjdC8q5d0gBXkQSS9YAcE24un0tXlYn69EU4EUksTSvR/Pgyr9FgrqWLGibAryIJBSX1R+AP699j6IVZQT61y7qAAAFwUlEQVSDwRZLFgSDwTjXsOfQKBoRSSjWayAAl03qzS0lFRSXVAChJQtunTWB21eWk50R4Mb/M55gMIjP54scW6dxzCyqc7G83swix66mAC8iiSUrlKK54jO9ueXdIy9nZwSY87MSyiqryRuaDThWl++kT0aAg4cbOX/Cp3h5YxW40IYVZsbZYwexunxHu+dief3B+iayMwJUH24kJzONhTPGdem3SgFeRBJLcw7+udJyYErk5Qde2AjAgF5plFVWs31/HXtqGsgI+DjcGKRyXy17axtbvNWW3YfYU9PQ7rlYXR+uT97QbMoqqymcntvlLXkFeBFJKC4jB4ePKR8/wtq+TzGgdzoVVYeOFGgE0j1Hmo9Nnuety7Z37gSufz6Qz52NVwKwp6YBIBLcu2OClgK8iCQU8wd4ZdS/0Ld6A/9wcg7vfHKATe5gxxfGwRb3qTZf767ZtwrwIpJwvlD4w9DomZXlFH9SEelgDefge7qiFWXdEuQ1TFJEEpLP5yMnM43C6bksmj2RO1a9R1llNQN6pQFEjhmBUJjrn3V0ezZcpr1zsbo+/DxvaHa3rYCpFryIJKyFM8ZFOirDwT48SiU7IxDzUTCxGkUzZczAblkB0+I5vTc/P9+VlpbG7fNFJLm0HmeerOPgzWytcy6/o3JqwYtI0ggHzfDR5/O1OLYVVKM5F+vru2sFTOXgRUSSlAK8iEiSUoAXEUlSCvAiIklKAV5EJEkpwIuIJCkFeBGRJBXXiU5mtgv4qAveehBQ1QXv29Okyn1C6txrqtwnpM69dsV9jnLODe6oUFwDfFcxs9JoZnklulS5T0ide02V+4TUudd43qdSNCIiSUoBXkQkSSVrgF8S7wp0k1S5T0ide02V+4TUude43WdS5uBFRCR5W/AiIilPAV5EJEkpwIuIJKmEDfBmNtPM3jezjWZ2c2fLJIKO7sPMRpjZX8yszMzWm9m/xqOeJyran5eZ+c3sLTNb0Z31i5Uo/+32M7NlZvaemZWb2dTurmcsRHmvC5v/3b5rZo+bWWZ31/NEmdmvzGynmb3bTpnuj0fOuYR7AH7gQ2AMkA68DeQdb5lEeER5r0OB05q/zgY2JNq9Hs/PC7gR+C2wIt717qr7BB4Frm7+Oh3oF++6d8W9AsOAzUBW8/Mngfnxrnsn7vVs4DTg3RP5ucf6kagt+DOBjc65Tc65euAJYG4nyiSCDu/DOVfpnHuz+etqoJzQL04iiernZWbDgdnAw91cv1jp8D7NrC+hgPEIgHOu3jm3r9treuKi/R0MAFlmFgB6AZ90Yx1jwjn3ErCnnSJxiUeJGuCHAR97nm/l6IAWTZlEcFz3YWa5wKnA611aq9iL9j7vA24Cgt1RqS4QzX2OBnYBxc2pqIfNrHd3VTCGOrxX59w24B5gC1AJ7HfOPddtNew+cYlHiRrgpQ1m1gd4CljgnDsQ7/rEmpkVADudc2vjXZcuFiD05/7PnXOnAoeAhO1Dao+Z9SfUkh0NnAz0NrOr4lur5JGoAX4bMMLzfHjza8dbJhFEdR9mlkYouP/GOfd0N9UtlqK5z+nARWZWQehP3PPM7LHuqV7MRHOfW4GtzrnwX2HLCAX8RBPNvV4AbHbO7XLONQBPA9O6qX7dKT7xKN6dE53s0AgAmwj9rx/usJh0vGUS4RHlvRrwa+C+eNe3K++zVflzSMxO1qjuE3gZGN/89X8CP4l33bviXoEpwHpCuXcj1Ll8Q7zr3sn7zeXYnaxxiUcJ2YJ3zjUC1wPPEupQfNI5tx7AzFaZ2cntlUkk0dwroZbtPEIt2nXNj1lxq3QnRHmfCe847vMG4Ddm9ndgMnBnPOp7IqL8PX2d0F8obwLvEMoqJNwaNWb2OPAqMN7MtprZt5pfj2s80lo0IiJJKiFb8CIi0jEFeBGRJKUALyKSpBTgRUSSlAK8iEiSUoAXEUlSCvAiIklKAV5EJEn9fwc8ZQ7EgsWyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "v = FluidVars()\n",
    "plt.scatter(g.x, U[:,v.urho], marker=\"x\", color=\"C0\")\n",
    "plt.plot(sod[\"x\"], sod[\"rho\"], color=\"C1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "1. Run the problem without limiting the slopes to see how it compares\n",
    "\n",
    "2. Try a higher-order Runge-Kutta time integration methods to see how the problem changes\n",
    "\n",
    "3. Implement periodic boundary conditions and create a new set of initial conditions that just puts a low amplitude Gaussian pulse&mdash;this will create an acoustic wave that propagates through the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}