Euler_approximate.ipynb

{
 "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](Euler.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": [
    "To examine the Python code for this chapter, and for the exact Riemann solution, see:\n",
    "\n",
    " - [exact_solvers/euler.py](exact_solvers/euler.py) ...\n",
    "   [on github.](https://github.com/clawpack/riemann_book/blob/FA16/exact_solvers/euler.py)"
   ]
  },
  {
   "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](Shallow_water_approximate.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_\\ell) & = f(q_r) - f(q_\\ell).\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 (\\ref{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 (\\ref{EA:cons}) is an identity, satisfied independently of our choice of $\\hat{q}$.  The second equation is (using (\\ref{EA:EOS}))\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{\\gamma-3}{2}\\hat{u}^2 (\\rho_r - \\rho_\\ell) + (3-\\gamma)\\hat{u}(\\rho_r u_r - \\rho_\\ell u_\\ell) \\\\ + (\\gamma-1)\\left( \\frac{p_r-p_\\ell}{\\gamma-1} + \\frac{1}{2}(\\rho_r u_r^2 - \\rho_\\ell u_\\ell^2) \\right) & = \\rho_r u_r^2 - \\rho_\\ell u_\\ell^2 + p_r - p_\\ell,\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_\\ell)\\hat{u}^2 - 2(\\rho_r u_r - \\rho_\\ell u_\\ell) \\hat{u} + (\\rho_r u_r^2 - \\rho_\\ell u_\\ell^2) & = 0,\n",
    "\\end{align}\n",
    "\n",
    "with roots\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{u}_\\pm & = \\frac{\\rho_r u_r - \\rho_\\ell u_\\ell \\mp \\sqrt{\\rho_r \\rho_\\ell} (u_\\ell - u_r)}{\\rho_r - \\rho_\\ell} = \\frac{\\sqrt{\\rho_r} u_r \\pm \\sqrt{\\rho_\\ell} u_\\ell}{\\sqrt{\\rho_r}\\pm\\sqrt{\\rho_\\ell}}\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_\\ell$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we find $\\hat{H}$ by solving the last equation of (\\ref{EA:cons}), which reads  \n",
    "\\begin{align}\n",
    "    \\left( \\frac{\\gamma-1}{2}\\hat{u}^3 - \\hat{u}\\hat{H} \\right)(\\rho_r - \\rho_\\ell) \\\\ + \\left( \\hat{H} - (\\gamma-1)\\hat{u}^2 \\right)(\\rho_r u_r - \\rho_\\ell u_\\ell) + \\gamma \\hat{u}(E_r - E_\\ell) & = H_r u_r \\rho_r - H_\\ell u_\\ell \\rho_\\ell.\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}_{\\pm} & = \\frac{\\rho_r H_r (u_r - \\hat{u}_+) - \\rho_\\ell H_\\ell (u_\\ell - \\hat{u}_+)}{\\rho_r u_r - \\rho_\\ell u_\\ell - \\hat{u}_\\pm(\\rho_r -\\rho_\\ell)} \\\\\n",
    "    & = \\frac{\\rho_r H_r (u_r - \\hat{u}_+) - \\rho_\\ell H_\\ell (u_\\ell - \\hat{u}_+)}{\\pm\\sqrt{\\rho_r \\rho_\\ell}(u_r-u_\\ell)} \\\\\n",
    "    & = \\frac{\\rho_r H_r - \\rho_\\ell H_\\ell \\mp\\sqrt{\\rho_r \\rho_\\ell}(H_r - H_\\ell)}{\\rho_r - \\rho_\\ell} \\\\\n",
    "    & = \\frac{\\sqrt{\\rho_r}H_r \\pm \\sqrt{\\rho_\\ell} H_\\ell}{\\sqrt{\\rho_r}\\pm\\sqrt{\\rho_\\ell}}.\n",
    "\\end{align}  \n",
    "Once more, we take the plus sign in the final expression for $\\hat{H}$, giving the Roe averages\n",
    "$$\n",
    "\\hat{u} = \\frac{\\sqrt{\\rho_r} u_r + \\sqrt{\\rho_\\ell} u_\\ell}{\\sqrt{\\rho_r} + \\sqrt{\\rho_\\ell}},\n",
    "\\qquad \\hat{H} = \\frac{\\sqrt{\\rho_r}H_r + \\sqrt{\\rho_\\ell} H_\\ell}{\\sqrt{\\rho_r} + \\sqrt{\\rho_\\ell}}.\n",
    "$$"
   ]
  },
  {
   "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_\\ell & = \\sum_{p=1}^3 {\\mathcal W}_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_\\ell$.  We now have everything we need to implement the Roe solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = 'svg'\n",
    "import numpy as np\n",
    "from exact_solvers import euler\n",
    "from utils import riemann_tools as rt\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets\n",
    "State = euler.Primitive_State"
   ]
  },
  {
   "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",
    "    dq = 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)*dq[0]+u_hat*dq[1]-dq[2])\n",
    "    alpha3 = (dq[1] + (c_hat - u_hat)*dq[0] - c_hat*alpha2) / (2.*c_hat)\n",
    "    alpha1 = dq[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] + \\\n",
    "              (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] + \\\n",
    "              (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] + \\\n",
    "              (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": [
    "An implementation of this solver for use in Clawpack can be found [here](https://github.com/clawpack/riemann/blob/FA16/src/rp1_euler_with_efix.f90).  Recall that an exact Riemann solver for the Euler equations appears in [exact_solvers/euler.py](exact_solvers/euler.py)."
   ]
  },
  {
   "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": [
    "def compare_solutions(left_state, right_state, solvers=['Exact','HLLE']):\n",
    "    q_l = np.array(euler.primitive_to_conservative(*left_state))\n",
    "    q_r = np.array(euler.primitive_to_conservative(*right))\n",
    "\n",
    "    outputs = []\n",
    "    states = {}\n",
    "\n",
    "    for solver in solvers:\n",
    "        if solver.lower() == 'exact':\n",
    "            outputs.append(euler.exact_riemann_solution(q_l,q_r))\n",
    "        if solver.lower() == 'hlle':\n",
    "            outputs.append(Euler_hlle(q_l, q_r))\n",
    "            states['hlle'] = outputs[-1][0]\n",
    "        if solver.lower() == 'roe':\n",
    "            outputs.append(Euler_roe(q_l, q_r))\n",
    "            states['roe'] = outputs[-1][0]\n",
    "\n",
    "    plot_function = \\\n",
    "        rt.make_plot_function([val[0] for val in outputs],\n",
    "                              [val[1] for val in outputs],\n",
    "                              [val[2] for val in outputs],\n",
    "                              [val[3] for val in outputs],\n",
    "                              solvers, layout='vertical',\n",
    "                              variable_names=euler.primitive_variables,\n",
    "                              derived_variables=euler.cons_to_prim,\n",
    "                              vertical_spacing=0.15,\n",
    "                              show_time_legend=True)\n",
    "    \n",
    "    interact(plot_function,\n",
    "             t=widgets.FloatSlider(min=0,max=0.9,step=0.1,value=0.4));\n",
    "    \n",
    "    return states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left  = State(Density = 3.,\n",
    "              Velocity = 0.,\n",
    "              Pressure = 3.)\n",
    "right = State(Density = 1.,\n",
    "              Velocity = 0.,\n",
    "              Pressure = 1.)\n",
    "\n",
    "states = compare_solutions(left, right, solvers=['Exact','Roe'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler.phase_plane_plot(left, right, approx_states=states['roe'])"
   ]
  },
  {
   "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_\\ell=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(Density = 0.1,\n",
    "              Velocity = 0.,\n",
    "              Pressure = 0.1)\n",
    "right = State(Density = 1.,\n",
    "              Velocity = 1.,\n",
    "              Pressure = 1.)\n",
    "\n",
    "states = compare_solutions(left, right, solvers=['Exact','Roe'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler.phase_plane_plot(left, right, approx_states=states['roe'])"
   ]
  },
  {
   "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",
    "gamma = 1.4\n",
    "mu = 2*(M**2-1)/(M*(gamma+1.))\n",
    "right = 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 = State(Density = rho_l,\n",
    "             Velocity = u_l,\n",
    "            Pressure = p_l)\n",
    "\n",
    "states = compare_solutions(left, right, solvers=['Exact','Roe'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler.phase_plane_plot(left, right, approx_states=states['roe'])"
   ]
  },
  {
   "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",
    "Here is an example of a Riemann problem whose solution includes a transonic 2-rarefaction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left  = State(Density = 0.1,\n",
    "              Velocity = -2.,\n",
    "              Pressure = 0.1)\n",
    "right = State(Density = 1.,\n",
    "              Velocity = -1.,\n",
    "              Pressure = 1.)\n",
    "\n",
    "states = compare_solutions(left, right, solvers=['Exact','Roe'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that in the exact solution, the right edge of the rarefaction travels to the right.  In the Roe solution, all waves travel to the left.  As in the case of the shallow water equations, here too this behavior can lead to unphysical solutions when this approximate solver is used in a numerical discretization.  In order to correct this, we can split the single wave into two when a transonic rarefaction is present, in a way similar to what is done in the shallow water equations.  We do not go into details here."
   ]
  },
  {
   "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",
    "Again following Einfeldt, 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_\\ell - c_\\ell, \\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_\\ell) - s_2 q_r + s_1 q_\\ell}{s_1 - s_2}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Euler_hlle(q_l, q_r, gamma=1.4):\n",
    "    \"\"\"HLLE 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 \\\n",
    "              + 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(Density = 3.,\n",
    "              Velocity = 0.,\n",
    "              Pressure = 3.)\n",
    "right = State(Density = 1.,\n",
    "              Velocity = 0.,\n",
    "              Pressure = 1.)\n",
    "    \n",
    "states = compare_solutions(left, right, solvers=['Exact','HLLE'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler.phase_plane_plot(left, right, approx_states=states['hlle'])"
   ]
  },
  {
   "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(Density = 1.,\n",
    "              Velocity = -5.,\n",
    "              Pressure = 1.)\n",
    "right = State(Density = 1.,\n",
    "              Velocity = 1.,\n",
    "              Pressure = 1.)\n",
    "\n",
    "states = compare_solutions(left, right, solvers=['Exact', 'Roe'])"
   ]
  },
  {
   "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 HLLE solver, on the other hand, guarantees positivity of the density and pressure.  Since the HLLE wave speed in the case of a rarefaction is always the speed of the leading edge of the true rarefaction, and since the HLLE 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": {},
   "outputs": [],
   "source": [
    "left  = State(Density = 1.,\n",
    "                    Velocity = -10.,\n",
    "                    Pressure = 1.)\n",
    "right = State(Density = 1.,\n",
    "                    Velocity = 1.,\n",
    "                    Pressure = 1.)\n",
    "\n",
    "states = compare_solutions(left, right, solvers=['Exact', 'HLLE']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler.phase_plane_plot(left,right,approx_states=states['hlle'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again recall that we are only considering a single Riemann solution in this chapter.  In [FV_compare](FV_compare.ipynb) we observe the effect of using these approximate solvers in a full discretization."
   ]
  }
 ],
 "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.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}