{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximate solvers for the Euler equations of gas dynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we discuss approximate solvers for the one-dimensional Euler equations:\n",
    "\n",
    "\\begin{align}\n",
    "    \\rho_t + (\\rho u)_x & = 0 \\\\\n",
    "    (\\rho u)_t + (\\rho u^2 + p)_x & = 0 \\\\\n",
    "    E_t + ((E+p)u)_x & = 0.\n",
    "\\end{align}\n",
    "\n",
    "As in [Euler_equations.ipynb](Euler_equations.ipynb), we focus on the case of an ideal gas, for which the total energy is given by\n",
    "\n",
    "\\begin{align} \\label{EA:EOS}\n",
    "    E = \\frac{p}{\\gamma-1} + \\frac{1}{2}\\rho u^2.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roe solver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first derive a Roe solver for the Euler equations, following the same approach as in [shallow_water_approximate_solvers.ipynb](shallow_water_approximate_solvers.ipynb).  Namely, we assume that $\\hat{A} = f'(\\hat{q})$ for some average state $\\hat{q}$, and impose the condition of conservation:\n",
    "\n",
    "\\begin{align} \\label{EA:cons}\n",
    "    f'(\\hat{q}) (q_r - q_l) & = f(q_r) - f(q_l).\n",
    "\\end{align}\n",
    "\n",
    "We will need the following quantities:\n",
    "\n",
    "\\begin{align}\n",
    "q & = \\begin{pmatrix} \\rho \\\\ \\rho u \\\\ E \\end{pmatrix}, & f(q) & = \\begin{pmatrix} \\rho u \\\\ \\rho u^2 + p \\\\ H u \\rho \\end{pmatrix}, \\\\\n",
    "f'(\\hat{q}) & = \\begin{pmatrix} \n",
    "                0 & 1 & 0 \\\\ \n",
    "                \\frac{\\gamma-3}{2}\\hat{u}^2 & (3-\\gamma)\\hat{u} & \\gamma-1 \\\\\n",
    "                \\frac{\\gamma-1}{2}\\hat{u}^3 - \\hat{u}\\hat{H} & \\hat{H} - (\\gamma-1)\\hat{u}^2 & \\gamma \\hat{u} \\end{pmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "Here $H = \\frac{E+p}{\\rho}$ is the enthalpy.  We have rewritten most expressions involving $E$ in terms of $H$ because it simplifies the derivation that follows.  We now solve \\eqref{EA:cons} to find $\\hat{u}$ and $\\hat{H}$.  It turns out that, for the case of a polytropic ideal gas, the average density $\\hat{\\rho}$ plays no role in the Roe solver.\n",
    "\n",
    "The first equation of \\eqref{EA:cons} is an identity, satisfied independently of our choice of $\\hat{q}$.  The second equation is (using \\eqref{EA:EOS})\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{\\gamma-3}{2}\\hat{u}^2 (\\rho_r - \\rho_l) + (3-\\gamma)\\hat{u}(\\rho_r u_r - \\rho_l u_l) + (\\gamma-1)\\left( \\frac{p_r-p_l}{\\gamma-1} + \\frac{1}{2}(\\rho_r u_r^2 - \\rho_l u_l^2) \\right) & = \\rho_r u_r^2 - \\rho_l u_l^2 + p_r - p_l,\n",
    "\\end{align}\n",
    "\n",
    "which simplifies to a quadratic equation for $\\hat{u}$:\n",
    "\n",
    "\\begin{align} \\label{EA:u_quadratic}\n",
    "    (\\rho_r - \\rho_l)\\hat{u}^2 - 2(\\rho_r u_r - \\rho_l u_l) \\hat{u} + (\\rho_r u_r^2 - \\rho_l u_l^2) & = 0,\n",
    "\\end{align}\n",
    "\n",
    "with roots\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{u}_\\pm & = \\frac{\\rho_r u_r - \\rho_l u_l \\mp \\sqrt{\\rho_r \\rho_l} (u_l - u_r)}{\\rho_r - \\rho_l} = \\frac{\\sqrt{\\rho_r} u_r \\pm \\sqrt{\\rho_l} u_l}{\\sqrt{\\rho_r}\\pm\\sqrt{\\rho_l}}\n",
    "\\end{align}\n",
    "\n",
    "Notice that this is identical to the Roe average of the velocity for the shallow water equations, if we replace the density $\\rho$ with depth $h$.  As before, we choose the root $u_+$ since it is well defined for all values of $\\rho_r, \\rho_l$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we find $\\hat{H}$ by solving the last equation of \\eqref{EA:cons}, which reads  \n",
    "\\begin{align}\n",
    "    \\left( \\frac{\\gamma-1}{2}\\hat{u}^3 - \\hat{u}\\hat{H} \\right)(\\rho_r - \\rho_l) + \\left( \\hat{H} - (\\gamma-1)\\hat{u}^2 \\right)(\\rho_r u_r - \\rho_l u_l) + \\gamma \\hat{u}(E_r - E_l) & = H_r u_r \\rho_r - H_l u_l \\rho_l.\n",
    "\\end{align}  \n",
    "We can simplify this using the equality $\\gamma E = \\rho H + \\frac{\\gamma-1}{2}\\rho u^2$ and solve for $\\hat{H}$ to find  \n",
    "\\begin{align}\n",
    "    \\hat{H} & = \\frac{\\rho_r H_r (u_r - \\hat{u}_+) - \\rho_l H_l (u_l - \\hat{u}_+)}{\\rho_r u_r - \\rho_l u_l - \\hat{u}_\\pm(\\rho_r -\\rho_l)} \\\\\n",
    "    & = \\frac{\\rho_r H_r (u_r - \\hat{u}_+) - \\rho_l H_l (u_l - \\hat{u}_+)}{\\pm\\sqrt{\\rho_r \\rho_l}(u_r-u_l)} \\\\\n",
    "    & = \\frac{\\rho_r H_r - \\rho_l H_l \\mp\\sqrt{\\rho_r \\rho_l}(H_r - H_l)}{\\rho_r - \\rho_l} \\\\\n",
    "    & = \\frac{\\sqrt{\\rho_r}H_r \\pm \\sqrt{\\rho_l} H_l}{\\sqrt{\\rho_r}\\pm\\sqrt{\\rho_l}}.\n",
    "\\end{align}  \n",
    "Once more, we take the plus sign in the final expression for $\\hat{H}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To implement the Roe solver, we also need the eigenvalues and eigenvectors of the averaged flux Jacobian $f'(\\hat{q})$.  These are just the eigenvalues of the true Jacobian, evaluated at the averaged state:\n",
    "\\begin{align}\n",
    "    \\lambda_1 & = \\hat{u} - \\hat{c}, & \\lambda_2 & = \\hat{u} & \\lambda_3 & = \\hat{u} + \\hat{c},\n",
    "\\end{align}  \n",
    "\\begin{align}\n",
    "r_1 & = \\begin{bmatrix} 1 \\\\ \\hat{u}-\\hat{c} \\\\ \\hat{H}-\\hat{u}\\hat{c}\\end{bmatrix} &\n",
    "r_2 & = \\begin{bmatrix} 1 \\\\ \\hat{u} \\\\ \\frac{1}{2}\\hat{u}^2 \\end{bmatrix} &\n",
    "r_3 & = \\begin{bmatrix} 1 \\\\ \\hat{u}+\\hat{c} \\\\ \\hat{H}+\\hat{u}\\hat{c}\\end{bmatrix}.\n",
    "\\end{align}\n",
    "Here $\\hat{c} = \\sqrt{(\\gamma-1)(\\hat{H}-\\hat{u}^2/2)}$.\n",
    "\n",
    "Solving the system of equations\n",
    "\\begin{align}\n",
    "q_r - q_l & = \\sum_{p=1}^3 \\wave^p = \\sum_{p=1}^3 \\alpha_p r_p\n",
    "\\end{align}  \n",
    "for the wave strengths gives\n",
    "\\begin{align}\n",
    "    \\alpha_2 & = \\delta_1 + (\\gamma-1)\\frac{\\hat{u}\\delta_2 - \\delta_3}{\\hat{c}^2} \\\\\n",
    "    \\alpha_3 & = \\frac{\\delta_2 + (\\hat{c}-\\hat{u})\\delta_1 - \\hat{c}\\alpha_2}{2\\hat{c}} \\\\\n",
    "    \\alpha_1 & = \\delta_1 - \\alpha_2 - \\alpha_3,\n",
    "\\end{align}  \n",
    "where $\\delta = q_r - q_l$.  We now have everything we need to implement the Roe solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['font.size'] = 8\n",
    "figsize =(8,4)\n",
    "mpl.rcParams['figure.figsize'] = figsize\n",
    "\n",
    "import numpy as np\n",
    "from exact_solvers import Euler\n",
    "from clawpack import riemann\n",
    "from utils import riemann_tools\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import namedtuple\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets\n",
    "import matplotlib\n",
    "Primitive_State = namedtuple('State', Euler.primitive_variables)\n",
    "Conserved_State = namedtuple('State', Euler.conserved_variables)\n",
    "gamma = 1.4\n",
    "problem_data = {}\n",
    "problem_data['gamma'] = gamma\n",
    "problem_data['gamma1'] = gamma - 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roe_averages(q_l, q_r, gamma=1.4):\n",
    "    rho_sqrt_l = np.sqrt(q_l[0])\n",
    "    rho_sqrt_r = np.sqrt(q_r[0])\n",
    "    p_l = (gamma-1.)*(q_l[2]-0.5*(q_l[1]**2)/q_l[0])\n",
    "    p_r = (gamma-1.)*(q_r[2]-0.5*(q_r[1]**2)/q_r[0])\n",
    "    denom = rho_sqrt_l + rho_sqrt_r\n",
    "    u_hat = (q_l[1]/rho_sqrt_l + q_r[1]/rho_sqrt_r)/denom\n",
    "    H_hat = ((q_l[2]+p_l)/rho_sqrt_l + (q_r[2]+p_r)/rho_sqrt_r)/denom\n",
    "    c_hat = np.sqrt((gamma-1)*(H_hat-0.5*u_hat**2))\n",
    "    \n",
    "    return u_hat, c_hat, H_hat\n",
    "    \n",
    "    \n",
    "def Euler_roe(q_l, q_r, gamma=1.4):\n",
    "    \"\"\"\n",
    "    Approximate Roe solver for the Euler equations.\n",
    "    \"\"\"\n",
    "    \n",
    "    rho_l = q_l[0]\n",
    "    rhou_l = q_l[1]\n",
    "    u_l = rhou_l/rho_l\n",
    "    rho_r = q_r[0]\n",
    "    rhou_r = q_r[1]\n",
    "    u_r = rhou_r/rho_r\n",
    "    \n",
    "    u_hat, c_hat, H_hat = roe_averages(q_l, q_r, gamma)\n",
    "    \n",
    "    delta = q_r - q_l\n",
    "    \n",
    "    s1 = u_hat - c_hat\n",
    "    s2 = u_hat\n",
    "    s3 = u_hat + c_hat\n",
    "    \n",
    "    alpha2 = (gamma-1.)/c_hat**2 *((H_hat-u_hat**2)*delta[0]+u_hat*delta[1]-delta[2])\n",
    "    alpha3 = (delta[1] + (c_hat - u_hat)*delta[0] - c_hat*alpha2) / (2.*c_hat)\n",
    "    alpha1 = delta[0] - alpha2 - alpha3\n",
    "    \n",
    "    r1 = np.array([1., u_hat-c_hat, H_hat - u_hat*c_hat])\n",
    "    r2 = np.array([1., u_hat, 0.5*u_hat**2])\n",
    "    q_l_star = q_l + alpha1*r1\n",
    "    q_r_star = q_l_star + alpha2*r2\n",
    "    \n",
    "    states = np.column_stack([q_l,q_l_star,q_r_star,q_r])\n",
    "    speeds = [s1, s2, s3]\n",
    "    wave_types = ['contact','contact', 'contact']\n",
    "    \n",
    "    def reval(xi):\n",
    "        rho = (xi<s1)*states[0,0] + (s1<=xi)*(xi<s2)*states[0,1] + (s2<=xi)*(xi<s3)*states[0,2] + (s3<=xi)*states[0,3]\n",
    "        mom = (xi<s1)*states[1,0] + (s1<=xi)*(xi<s2)*states[1,1] + (s2<=xi)*(xi<s3)*states[1,2] + (s3<=xi)*states[1,3]\n",
    "        E = (xi<s1)*states[2,0] + (s1<=xi)*(xi<s2)*states[2,1] + (s2<=xi)*(xi<s3)*states[2,2] + (s3<=xi)*states[2,3]\n",
    "        return rho, mom, E\n",
    "    \n",
    "    return states, speeds, reval, wave_types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examples\n",
    "\n",
    "Let's compare the Roe approximation to the exact solution.  As a first example, we use the Sod shock tube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left_state  = Primitive_State(Density = 3.,\n",
    "                              Velocity = 0.,\n",
    "                              Pressure = 3.)\n",
    "right_state = Primitive_State(Density = 1.,\n",
    "                              Velocity = 0.,\n",
    "                              Pressure = 1.)\n",
    "\n",
    "q_l = np.array(Euler.primitive_to_conservative(*left_state))\n",
    "q_r = np.array(Euler.primitive_to_conservative(*right_state))\n",
    "\n",
    "states, speeds, reval, wave_types = Euler_roe(q_l, q_r)\n",
    "ex_states, ex_speeds, ex_reval, ex_wt = Euler.exact_riemann_solution(q_l,q_r)\n",
    "\n",
    "plot_function = riemann_tools.make_plot_function([ex_states, states],\n",
    "                                                 [ex_speeds, speeds],\n",
    "                                                 [ex_reval, reval],\n",
    "                                                 [ex_wt, wave_types],\n",
    "                                                 variable_names=Euler.primitive_variables,\n",
    "                                                 derived_variables=Euler.cons_to_prim)\n",
    "\n",
    "interact(plot_function, t=widgets.FloatSlider(min=0,max=0.9,step=0.1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Euler.phase_plane_plot(q_l,q_r,approx_states=states,cons_inputs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that in the true solution the middle wave is a contact discontinuity and carries only a jump in the density. For that reason the three-dimensional phase space plot is generally shown projected onto the pressure-velocity plane as shown above: The two intermediate states in the true solution have the same pressure and velocity, and so are denoted by a single Middle state in the phase plane plot.  \n",
    "\n",
    "The Roe solver, on the other hand, generates a middle wave that carries a jump in all 3 variables  and there are two green dots appearing in the plot above for the two middle states (though the pressure jump is quite small in this example).  For a Riemann problem like this one with zero initial velocity on both sides, the Roe average velocity must also be zero, so the middle wave is stationary; this is of course not typically true in the exact solution, even when $u_l=u_r=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a second example.  Experiment with the initial states to explore how the Roe solution compares to the exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left_state  = Primitive_State(Density = 0.1,\n",
    "                              Velocity = 0.,\n",
    "                              Pressure = 0.1)\n",
    "right_state = Primitive_State(Density = 1.,\n",
    "                              Velocity = 1.,\n",
    "                              Pressure = 1.)\n",
    "\n",
    "q_l = np.array(Euler.primitive_to_conservative(*left_state))\n",
    "q_r = np.array(Euler.primitive_to_conservative(*right_state))\n",
    "\n",
    "states, speeds, reval, wave_types = Euler_roe(q_l, q_r)\n",
    "ex_states, ex_speeds, ex_reval, ex_wt = Euler.exact_riemann_solution(q_l,q_r)\n",
    "\n",
    "plot_function = riemann_tools.make_plot_function([ex_states, states],\n",
    "                                                 [ex_speeds, speeds],\n",
    "                                                 [ex_reval, reval],\n",
    "                                                 [ex_wt, wave_types],\n",
    "                                                 variable_names=Euler.primitive_variables,\n",
    "                                                 derived_variables=Euler.cons_to_prim)\n",
    "\n",
    "interact(plot_function, t=widgets.FloatSlider(min=0,max=0.9,step=0.1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Euler.phase_plane_plot(q_l,q_r,approx_states=states,cons_inputs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single-shock solution\n",
    "Next we demonstrate the exactness property of the Roe solver by applying it to a case where the left and right states are connected by a single shock wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 2.  # Mach number of the shock wave\n",
    "mu = 2*(M**2-1)/(M*(gamma+1.))\n",
    "right = Primitive_State(Density = 1.,\n",
    "                        Velocity = 0.,\n",
    "                        Pressure = 1.)\n",
    "c_r = np.sqrt(gamma*right.Pressure/right.Density)\n",
    "\n",
    "rho_l = right.Density * M/(M-mu)\n",
    "p_l = right.Pressure * ((2*M**2-1)*gamma+1)/(gamma+1)\n",
    "u_l = mu*c_r\n",
    "\n",
    "left = Primitive_State(Density = rho_l,\n",
    "                       Velocity = u_l,\n",
    "                       Pressure = p_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = np.array(Euler.primitive_to_conservative(*left))\n",
    "q_r = np.array(Euler.primitive_to_conservative(*right))\n",
    "\n",
    "states, speeds, reval, wave_types = Euler_roe(q_l, q_r)\n",
    "ex_states, ex_speeds, ex_reval, ex_wt = Euler.exact_riemann_solution(q_l,q_r)\n",
    "\n",
    "plot_function = riemann_tools.make_plot_function([ex_states, states],\n",
    "                                                 [ex_speeds, speeds],\n",
    "                                                 [ex_reval, reval],\n",
    "                                                 [ex_wt, wave_types],\n",
    "                                                 variable_names=Euler.primitive_variables,\n",
    "                                                 derived_variables=Euler.cons_to_prim)\n",
    "\n",
    "interact(plot_function, t=widgets.FloatSlider(min=0,max=0.9,step=0.1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Euler.phase_plane_plot(q_l,q_r,approx_states=states,cons_inputs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is evident that the solution consists of a single right-going shock.  The exact solution cannot be seen because it coincides exactly with the Roe solution.  The path of the shock in the first plot also cannot be seen since it is plotted under the path of the rightmost Roe solution wave.  The two solutions differ only in the wave speeds predicted for the other two waves, but since these waves have zero strength this makes no difference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transonic rarefactions and an entropy fix\n",
    "To be written."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HLLE Solver\n",
    "\n",
    "Recall that an HLL solver uses only two waves with a constant state between them.  The Euler equations are our first example for which the number of waves in the true solution is larger than the number of waves in the approximate solution.  As one might expect, this leads to noticeable inaccuracy in solutions produced by the solver.\n",
    "\n",
    "The left-going wave speed is chosen to be the minimum of the Roe speed for the 1-wave and the characterstic speed $\\lambda^1$ in the left state $q_\\ell$.  The right-going wave speed is chosen to be the maximum of the Roe speed for the 3-wave and the characteristic speed $\\lambda^3$ in the right state $q_r$.  Effectively, this means that\n",
    "\\begin{align}\n",
    "    s_1 & = \\min(u_l - c_l, \\hat{u}-\\hat{c}) \\\\\n",
    "    s_2 & = \\max(u_r + c_r, \\hat{u}+\\hat{c})\n",
    "\\end{align}\n",
    "\n",
    "Recall that once we have chosen these two wave speeds, conservation dictates the value of the intermediate state:  \n",
    "\\begin{align}  \\label{SWA:hll_middle_state}\n",
    "q_m = \\frac{f(q_r) - f(q_l) - s_2 q_r + s_1 q_l}{s_1 - s_2}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Euler_hll(q_l, q_r, gamma=1.4):\n",
    "    \"\"\"HLL approximate solver for the Euler equations.\"\"\"\n",
    "    \n",
    "    rho_l = q_l[0]\n",
    "    rhou_l = q_l[1]\n",
    "    u_l = rhou_l/rho_l\n",
    "    rho_r = q_r[0]\n",
    "    rhou_r = q_r[1]\n",
    "    u_r = rhou_r/rho_r\n",
    "    E_r = q_r[2]\n",
    "    E_l = q_l[2]\n",
    "    \n",
    "    u_hat, c_hat, H_hat = roe_averages(q_l, q_r, gamma)\n",
    "    p_r = (gamma-1.) * (E_r - rho_r*u_r**2/2.)\n",
    "    p_l = (gamma-1.) * (E_l - rho_l*u_l**2/2.)\n",
    "    H_r = (E_r+p_r) / rho_r\n",
    "    H_l = (E_l+p_l) / rho_l\n",
    "    c_r = np.sqrt((gamma-1.)*(H_r-u_r**2/2.))\n",
    "    c_l = np.sqrt((gamma-1.)*(H_l-u_l**2/2.))\n",
    "    \n",
    "    s1 = min(u_l-c_l,u_hat-c_hat)\n",
    "    s2 = max(u_r+c_r,u_hat+c_hat)\n",
    "    \n",
    "    rho_m = (rhou_r - rhou_l - s2*rho_r + s1*rho_l)/(s1-s2)\n",
    "    rhou_m = (rho_r*u_r**2 - rho_l*u_l**2 + p_r - p_l - s2*rhou_r + s1*rhou_l)/(s1-s2)\n",
    "    E_m = ( u_r*(E_r+p_r) - u_l*(E_l+p_l) -s2*E_r + s1*E_l)/(s1-s2)\n",
    "    q_m = np.array([rho_m, rhou_m, E_m])\n",
    "    \n",
    "    states = np.column_stack([q_l,q_m,q_r])\n",
    "    speeds = [s1, s2]\n",
    "    wave_types = ['contact','contact']\n",
    "    \n",
    "    def reval(xi):\n",
    "        rho  = (xi<s1)*rho_l + (s1<=xi)*(xi<=s2)*rho_m + (s2<xi)*rho_r\n",
    "        mom  = (xi<s1)*rhou_l + (s1<=xi)*(xi<=s2)*rhou_m + (s2<xi)*rhou_r\n",
    "        E  = (xi<s1)*E_l + (s1<=xi)*(xi<=s2)*E_m + (s2<xi)*E_r\n",
    "        return rho, mom, E\n",
    "\n",
    "    return states, speeds, reval, wave_types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left_state  = Primitive_State(Density = 3.,\n",
    "                              Velocity = 0.,\n",
    "                              Pressure = 3.)\n",
    "right_state = Primitive_State(Density = 1.,\n",
    "                              Velocity = 0.,\n",
    "                              Pressure = 1.)\n",
    "\n",
    "q_l = np.array(Euler.primitive_to_conservative(*left_state))\n",
    "q_r = np.array(Euler.primitive_to_conservative(*right_state))\n",
    "\n",
    "states, speeds, reval, wave_types = Euler_hll(q_l, q_r)\n",
    "ex_states, ex_speeds, ex_reval, ex_wt = Euler.exact_riemann_solution(q_l,q_r)\n",
    "\n",
    "plot_function = riemann_tools.make_plot_function([ex_states, states],\n",
    "                                                 [ex_speeds, speeds],\n",
    "                                                 [ex_reval, reval],\n",
    "                                                 [ex_wt, wave_types],\n",
    "                                                 variable_names=Euler.primitive_variables,\n",
    "                                                 derived_variables=Euler.cons_to_prim)\n",
    "\n",
    "interact(plot_function, t=widgets.FloatSlider(min=0,max=0.9,step=0.1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Euler.phase_plane_plot(q_l,q_r,approx_states=states,cons_inputs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preservation of positivity\n",
    "Just as we saw in the case of the shallow water equations, the Roe solver (or any linearized solver) for the Euler equations fails to preserve positivity of the pressure and/or density in some situations.  Here is one example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left_state  = Primitive_State(Density = 1.,\n",
    "                              Velocity = -5.,\n",
    "                              Pressure = 1.)\n",
    "right_state = Primitive_State(Density = 1.,\n",
    "                              Velocity = 1.,\n",
    "                              Pressure = 1.)\n",
    "\n",
    "\n",
    "q_l = np.array(Euler.primitive_to_conservative(*left_state))\n",
    "q_r = np.array(Euler.primitive_to_conservative(*right_state))\n",
    "\n",
    "states, speeds, reval, wave_types = Euler_roe(q_l, q_r)\n",
    "ex_states, ex_speeds, ex_reval, ex_wt = Euler.exact_riemann_solution(q_l,q_r)\n",
    "\n",
    "plot_function = riemann_tools.make_plot_function([ex_states, states],\n",
    "                                                 [ex_speeds, speeds],\n",
    "                                                 [ex_reval, reval],\n",
    "                                                 [ex_wt, wave_types],\n",
    "                                                 variable_names=Euler.primitive_variables,\n",
    "                                                 derived_variables=Euler.cons_to_prim)\n",
    "\n",
    "interact(plot_function, t=widgets.FloatSlider(min=0,max=0.9,step=0.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, in this example each Roe solver wave moves much more slowly than the leading edge of the corresponding true rarefaction.  In order to maintain conservation, this implies that the middle Roe state must have lower density than the true middle state.  This leads to a negative density.  Note that the velocity and pressure take huge values in the intermediate state.\n",
    "\n",
    "The HLL solver, on the other hand, guarantees positivity of the density and pressure.  Since the HLL wave speed in the case of a rarefaction is always the speed of the leading edge of the true rarefaction, and since the HLL solution is conservative, the density in a rarefaction will always be at least as great as that of the true solution.  This can be seen clearly in the example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "left_state  = Primitive_State(Density = 1.,\n",
    "                              Velocity = -10.,\n",
    "                              Pressure = 1.)\n",
    "right_state = Primitive_State(Density = 1.,\n",
    "                              Velocity = 1.,\n",
    "                              Pressure = 1.)\n",
    "\n",
    "\n",
    "q_l = np.array(Euler.primitive_to_conservative(*left_state))\n",
    "q_r = np.array(Euler.primitive_to_conservative(*right_state))\n",
    "\n",
    "states, speeds, reval, wave_types = Euler_hll(q_l, q_r)\n",
    "ex_states, ex_speeds, ex_reval, ex_wt = Euler.exact_riemann_solution(q_l,q_r)\n",
    "\n",
    "plot_function = riemann_tools.make_plot_function([ex_states, states],\n",
    "                                                 [ex_speeds, speeds],\n",
    "                                                 [ex_reval, reval],\n",
    "                                                 [ex_wt, wave_types],\n",
    "                                                 ['Exact','HLL'],\n",
    "                                                 layout='vertical',\n",
    "                                                 variable_names=Euler.primitive_variables,\n",
    "                                                 derived_variables=Euler.cons_to_prim)\n",
    "\n",
    "interact(plot_function, t=widgets.FloatSlider(min=0,max=0.9,step=0.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*I've used a vertical layout in the example above.  Which should we use in general?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Euler.phase_plane_plot(q_l,q_r,approx_states=states,cons_inputs=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}