Euler_equations_TammannEOS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nearly-incompressible flow: the Tammann equation of state\n",
    "\n",
    "In this chapter, we consider the one-dimensional Riemann problem for the Euler equations with the Tammann equation of state (EOS). The Tammann EOS, also known as stiffened gas EOS, is used to model almost-incompressible fluids. It can be physically understood as modeling an ideal gas under very high ambient pressures. The Riemann problem we want to solve consists of the one-dimensional Euler equations,\n",
    "\n",
    "\\begin{align*}\n",
    "    \\left[\\begin{array}{c} \\rho \\\\   \\rho u \\\\    E   \\end{array} \\right]_t +\n",
    "    \\left[\\begin{array}{c} \\rho u \\\\   \\rho u^2 + p \\\\     u(E+p)   \\end{array} \\right]_x\n",
    "=   0,\n",
    "\\end{align*}\n",
    "where $\\rho$ is density, $u$ the velocity, $E$ the internal energy, and $p$ the pressure change relative to an ambient pressure $p_{\\infty}$. As usual, the subcripts $x,t$ denote partial derivatives. The initial conditions are given by the left and right constant states $q_l=[\\rho_l, \\rho_l u_l, E_l]$ and $q_r=[\\rho_r, \\rho_r u_r, E_r]$. The Tamann EOS is given by \n",
    "\n",
    "\\begin{align*}\n",
    "p (\\rho,e)=\\rho e(\\gamma -1) - \\gamma p_{\\infty},\n",
    "\\end{align*}\n",
    "where $e$ is the specific internal energy. For a general equation of state $p (\\rho,e)$, the speed of sound is given by the isentropic derivative of the pressure with respect to density (keeping entropy constant),\n",
    "\n",
    "\\begin{align*}\n",
    "  c = \\sqrt{\\frac{\\partial p(\\rho,e)}{\\partial \\rho} \\bigg|_s} = \\sqrt{\\frac{\\partial p(\\rho,e)}{\\partial \\rho} + \\frac{p(\\rho,e)}{\\rho^2}\\frac{\\partial p(\\rho,e)}{\\partial e}}.\n",
    "\\end{align*}\n",
    "In the case of the Tammann EOS, we obtain\n",
    "\n",
    "\\begin{align*}\n",
    " c= \\sqrt{\\gamma\\frac{p+p_{\\infty }}{\\rho}}.\n",
    "\\end{align*}\n",
    "Note that this is almost the same as the [sound speed in an ideal gas](Euler_equations.ipynb#Characteristic-velocities) except for the additional $p_{\\infty}$ term. This means that the fluid behaves similarly to an ideal gas that is already at pressure $p_{\\infty }$. When $p_{\\infty } \\gg 1$, the Tammann EOS emulates an ideal gas at a very high pressure, which can be used to model almost-incompressible fluids like water. They are almost-incompressible because a very large change is pressure is required to produce even a small change in density.  Note the sound speed (as well as the density and any other physical quantity) is well defined for negative values of $p$, with the only limitation that $p+p_{\\infty}>0$. We will revisit this issue in [the last section of this chapter](#Strong-rarefaction-waves).\n",
    "\n",
    "Using the equation of state, the state of the system can also be written in the primitive variables $[\\rho, u, p]$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact solution\n",
    "\n",
    "The exact solution to the Riemann problem with the Tammann EOS when the fluid parameters $\\gamma$ and $p_\\infty$ do not vary can be found in <cite data-cite=\"Ivings1998\"><a href=\"riemann.html#Ivings1998\">(Ivings & Toro 1998)</a></cite>. For the more general case when the left and right states may have different values of $\\gamma$ and $p_\\infty$, the solution can be found in the supplemental material of the work <cite data-cite=\"delrazo2017\"><a href=\"riemann.html#delrazo2017\">(del Razo, 2017)</a></cite>. In this chapter we will focus on the general case. The presentation follows that of <cite data-cite=\"delrazo2017\"><a href=\"riemann.html#delrazo2017\">(del Razo, 2017)</a></cite>.\n",
    "\n",
    "Just as in [the earlier chapter on ideal gas flow](Euler_equations.ipynb), the solution of the Riemann problem considered here consists of three waves.  The 1-wave and 3-wave are genuinely nonlinear (rarefactions or shocks), while the 2-wave is a linearly degenerate contact discontinuity. The solution thus includes four constant regions with solutions $q_l, q_{*l}, q_{*r}, q_r$, separated by the three waves. Just as in the case of the ideal gas, the velocity and pressure do not change across the contact wave, so we denote them simply by $u_*$ and $p_*$. Note that the interface between the fluids that were initially to the left and right of $x=0$ moves along the contact discontinuity, so if the fluid parameters $\\gamma$ and $p_\\infty$ vary between the left and right states then the this jump in parameters also moves with the contact discontinuity.  Consequently, the 1-wave propagates through left-state fluid and the 3-waves propagates through right-state fluid.  (This will be true even in the supersonic case, where all the wave speeds have the same sign.) This allows us to determine the Hugoniot loci and integral curves for these two nonlinear waves using the fact that fluid parameters $\\gamma$ and $p_\\infty$ are always constant across each of these waves.\n",
    "\n",
    "In order to determine whether the 1-wave and 3-wave are rarefactions or shocks, we will need to create a function of the \n",
    "middle state pressure $p_*$ that ensures that the velocity $u_*$ across the contact discontinuity is consistent. Since we know that the left state velocity $u_l$ should be connected by a rarefaction or shock to $u_*$, we write $u_*=u_l + [u]_1$, where $[u]_1$ is the jump of the velocity across the 1-wave. In a similar manner, we also \n",
    "know the 3-wave should be a shock or rarefaction, so we write $u_*= u_r - [u]_3$.  Therefore, it is useful to define\n",
    "\n",
    "\\begin{align}\n",
    "  \\phi_l(p) &= u_l - \\mathcal{F}_l(p),   \\label{CDvel-1}  \\\\\n",
    "  \\phi_r(p) &= u_r + \\mathcal{F}_r(p).   \\label{CDvel-2}\n",
    "\\end{align}\n",
    "\n",
    "where the value of $\\mathcal{F}_{l,r}(p) = -[u]_{1,3}$ depends on whether the wave in question is a shock or a rarefaction (signs were \n",
    "chosen for notational consistency). As we expect these two equations to yield the same contact discontinuity \n",
    "velocity $u_* = \\phi_l(p_*) = \\phi_r(p_*)$, we have \n",
    "\n",
    "\\begin{align}\n",
    " \\Phi(p_*)= \\phi_r(p_*) - \\phi_l(p_*) = 0.\n",
    " \\label{Phi} \n",
    "\\end{align}\n",
    "\n",
    "This nonlinear equation will yield the pressure $p_*$ that provides consistency between the type of waves (rarefactions or shocks), their speeds and the contact discontinuity velocity $u_*$. As we mentioned before, the shape of $\\phi_k(p_*)$ will depend on whether the states are connected by a shock wave or rarefaction. Once \n",
    "the $p_*$ has been found, the contact discontinuity velocity can be found from the expressions just derived. However, it is not yet clear how to calculate the density or the speeds of the 1-wave \n",
    "and 3-wave. It will become obvious how to obtain these quantities once we write the explicit equations for our system further below.\n",
    "\n",
    "Before coding the exact solution, we should first note that having a rarefaction or shock in the 1-wave \n",
    "and 3-wave will depend on the pressure $p_*$, in [the same way it did for an ideal gas](Euler_equations.ipynb#Entropy-condition).  If the pressure is higher on the side toward which the wave is propagating (i.e., in front of the wave), it must be a rarefaction. If the pressure is lower in front of the wave, it must be a shock. In the Euler equations, this yields four possible cases for the value $\\Phi(p_*)$, just as in the solution with an ideal gas EOS:\n",
    "\n",
    " * 1-rarefaction, 3-rarefaction: $p_*< p_l$ and $p_*<p_r$  \n",
    " $\\Phi(p_*)= \\phi_r^R(p_*) - \\phi_l^R(p_*)$,\n",
    " * 1-shock,       3-rarefaction $p_l \\le p_* \\le p_r$  \n",
    " $\\Phi(p_*)= \\phi_r^R(p_*) - \\phi_l^S(p_*)$,\n",
    " * 1-rarefaction,       3-shock $p_r \\le p_* \\le p_l$  \n",
    " $\\Phi(p_*)= \\phi_r^S(p_*) - \\phi_l^R(p_*)$,\n",
    " * 1-shock, 3-shock: $p_* > p_l$ and $p_*>p_r$  \n",
    " $\\Phi(p_*)= \\phi_r^S(p_*) - \\phi_l^S(p_*)$,\n",
    " \n",
    "where the index $S,R$ indicates whether the function $\\phi$ was obtained by using the Rankine-Hugoniot equations to connect states by shocks or the Riemann invariants to connect them by rarefactions respectively.\n",
    "\n",
    "We now derive the functions $\\phi_{k}^{\\mu}$ for all four cases with $k=l,r$ (left and right) and $\\mu=R,S$ (rarefaction and shock).  We also show how to obtain the density and the missing wave speeds. We write $s_n$ for the speed of the wave belonging to family $n$ ($n=1,2,3$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rankine-Hugoniot conditions for shock waves\n",
    "Let us first consider the case in which both the 1-wave and the 3-wave are shocks. The velocity of each\n",
    "shock will be given by the Rankine-Hugoniot conditions. These conditions will also allow us to calculate $\\phi_k^S(p_*)$ with $k=l,r$.  The following analysis is written to be applicable to both the 1-wave \n",
    "velocity $s_1$ and the 3-wave velocity $s_3$, by considering $s_n$ and $q_k$, with $n=1, k=l$ for the 1-wave and $s=3, k=r$ for the 3-wave.\n",
    "\n",
    "The Rankine-Hugoniot conditions for \n",
    "the Euler equations can be written as (see e.g. <cite data-cite=\"Ivings1998\"><a href=\"riemann.html#Ivings1998\">(Ivings & Toro 1998)</a></cite>),\n",
    "\n",
    "\\begin{align}\n",
    " \\rho_k \\omega_k = \\rho_{*k} \\omega_*, \\label{RH-mass} \\\\\n",
    " \\rho_k \\omega_k^2 + p_k = \\rho_{*k} \\omega_*^2 + p_{*k}, \\label{RH-mom}  \\\\\n",
    " \\frac{1}{2} \\omega_k^2 + h_k = \\frac{1}{2} \\omega_*^2 + h_{*k}, \\label{RH-ene}\n",
    "\\end{align}\n",
    "where $\\omega_k = u_k - s_n$ , $\\omega_* = u_* - s_n$, and the specific enthalpy is given \n",
    "by $h = e+(p+p_\\infty)/\\rho$ with $e$ the specific internal energy that relates to the internal energy \n",
    "of our original variables by $E = \\rho e + \\rho u^2/2$. We will use these relations to find $\\phi_k^S(p_*)$ and the wave speeds $s_n$. \n",
    "\n",
    "#### Finding $\\phi_l^S(p_*)$ and $\\phi_r^S(p_*)$ and $s_1$ and $s_3$: \n",
    "We start by defining the mass fluxes $F_k$ for $k=l,r$ as\n",
    "\n",
    "\\begin{align}\n",
    " F_l = \\rho_l \\omega_l = \\rho_{*l} \\omega_* \\label{Ql} \\\\\n",
    " F_r = -\\rho_r \\omega_r = -\\rho_{*r} \\omega_*. \\label{Qr} \n",
    "\\end{align}\n",
    "\n",
    "As $\\omega_k = u_k - s_n$, from these two equations we can obtain the wave speeds in terms of $F_l$ and $F_r$,\n",
    "\n",
    "\\begin{align}\n",
    " s_1 = u_l - \\frac{F_l}{\\rho_l}, \\  \\  \\\n",
    " s_3 = u_r + \\frac{F_r}{\\rho_r}. \\label{RH-speeds} \n",
    "\\end{align}\n",
    "\n",
    "We still need to find $F_l$ and $F_r$, so we substitute equations (\\ref{Ql}) and (\\ref{Qr}) \n",
    "into (\\ref{RH-mom}) to obtain\n",
    "\n",
    "\\begin{align}\n",
    "  F_l &= \\frac{\\tilde{p}_{*l} -\\tilde{p}_l}{\\omega_l - \\omega_*} = \\frac{p_* -p_l}{u_l - u_*} \\label{Ql-1} \\\\\n",
    "  F_r &= -\\frac{\\tilde{p}_{*r} -\\tilde{p}_r}{\\omega_r - \\omega_*} = -\\frac{p_* -p_r}{u_r - u_*}, \\label{Qr-1} \n",
    "\\end{align}\n",
    "where $\\tilde{p}_\\kappa = p_\\kappa + p_{\\infty \\kappa}$ is defined to simplify future notation \n",
    "with $\\kappa = l, *l, *r$ and $r$. Note that $\\tilde{p}_{*l} \\neq \\tilde{p}_{*r}$ and \n",
    "that $\\tilde{p}_{*k} -\\tilde{p}_k=p_* -p_k$ since $p_*=p_{*l}=p_{*r}$ and $p_{\\infty k} = p_{\\infty *k}$  ($k=l,r)$.\n",
    "Solving for $u_*$ we obtain the equations\n",
    "\n",
    "\\begin{align}\n",
    " u_* = u_l - \\frac{\\tilde{p}_{*l} - \\tilde{p}_l}{F_l} = \\phi_l^S(p_*) \\\\\n",
    " u_* = u_r + \\frac{\\tilde{p}_{*r} - \\tilde{p}_r}{F_r} = \\phi_r^S(p_*) \n",
    " \\label{RH-phis}\n",
    "\\end{align}\n",
    "\n",
    "Comparing to equations (\\ref{CDvel-1}, \\ref{CDvel-2}), we notice $\\mathcal{F}_k(p_*) = \\frac{p_* - p_k}{F_k}$. We also notice we have almost\n",
    "obtained the $\\phi$ functions we are looking for, though we still need to find $F_l$ and $F_r$ in terms of known \n",
    "variables.\n",
    "\n",
    "#### Finding $F_l$ and $F_r$: \n",
    "\n",
    "From equations (\\ref{Ql}) and  (\\ref{Ql-1}), we know\n",
    "that\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\tilde{p}_{*l} -\\tilde{p}_l}{\\omega_l - \\omega_*} = \\rho_l \\omega_l.\n",
    "\\end{align*}\n",
    "\n",
    "Solving for $\\omega_*$, substituting the solution into (\\ref{Ql}) and subtituting the $w_l$ for $F_l/\\rho_l$, \n",
    "we obtain a new equation that we can solve for $F_l$ that yields\n",
    "\n",
    "\\begin{align}\n",
    " F_k = \\sqrt{\\rho_k \\rho_{*k} \\frac{\\tilde{p}_{*k} -\\tilde{p}_k}{\\rho_{*k} - \\rho_k}}. \\label{Qk} \n",
    "\\end{align}\n",
    "with $k=l,r$, since we repeated the same process for $F_r$ and obtained exactly the same equation. \n",
    "However, we still don't know $\\rho_{*k}$, for this we will need our third Rankine-Hugoniot condition (\\ref{RH-ene}).\n",
    "\n",
    "#### Finding $\\rho_{*l}$ and $\\rho_{*r}$: \n",
    "\n",
    "From equation (\\ref{RH-ene}), we can obtain\n",
    "\n",
    "\\begin{align}\n",
    " h_{*k} - h_k &= \\frac{1}{2}\\left(w_k^2 - w_*^2\\right) \\nonumber\\\\\n",
    "              &= \\frac{1}{2}\\left(\\pm \\frac{F_k^2}{\\rho_k^2} \\mp \\frac{F_k^2}{\\rho_{*k}^2}\\right) \\nonumber \\\\\n",
    "              &= \\frac{1}{2}\\left(\\frac{1}{\\rho_k} + \\frac{1}{\\rho_{*k}}\\right)(\\tilde{p}_{*k}-\\tilde{p}_k), \\label{delta-h} \n",
    "\\end{align}\n",
    "\n",
    "where the sign above is used for $k=l$ and the one below for $k=r$, and we used equations \n",
    "(\\ref{Ql}) and (\\ref{Qr}) for the second line and (\\ref{Qk}) for the third line. We can now substitute \n",
    "the specific enthalpy $h=\\gamma \\tilde{p}/(\\rho(\\gamma -1))$\n",
    "in equation (\\ref{delta-h}) to obtain\n",
    "\n",
    "\\begin{align*}\n",
    " \\frac{\\gamma_{*k}}{\\gamma_{*k} - 1}\\frac{\\tilde{p}_{*k}}{\\rho_{*k}} - \\frac{\\gamma_{k}}{\\gamma_{k} - 1}\\frac{\\tilde{p}_k}{\\rho_{k}}\n",
    " =\\frac{1}{2}\\left(\\frac{1}{\\rho_k} + \\frac{1}{\\rho_{*k}}\\right)(\\tilde{p}_{*k}-\\tilde{p}_k).\n",
    "\\end{align*}\n",
    "\n",
    "As the interface is the contact discontinuity, the jump in the parameters is only across the contact discontinuity, so \n",
    "$\\gamma_{*k}=\\gamma_k$. Now we can solve for the unknown density:\n",
    "\n",
    "\\begin{align}\n",
    " \\rho_{*k} = \\rho_k \\left(\\frac{\\frac{\\tilde{p_*}}{\\tilde{p_k}} + \n",
    " \\frac{\\gamma_k-1}{\\gamma_k+1}}{\\frac{\\tilde{p_*}}{\\tilde{p_k}}\\frac{\\gamma_k-1}{\\gamma_k+1}+1}\\right). \n",
    " \\label{rho-k}\n",
    "\\end{align}\n",
    "\n",
    "Substituting this result into equation (\\ref{Qk}), we obtain $F_k$ in terms of $p_*$ and known variables:\n",
    "\n",
    "\\begin{align}\n",
    "  F_k = \\sqrt{\\rho_k \\frac{\\tilde{p}_{*k} + \\tilde{p}_{k}\\frac{\\gamma_k - 1}{\\gamma_k +1}}{\\frac{2}{\\gamma_k+1}}}.\n",
    "\\label{Qk-f}\n",
    "\\end{align}\n",
    "\n",
    "Subsitituting (\\ref{Qk-f}) into (\\ref{RH-phis}), we can calculate the nonlinear functions $\\phi^S_{l,r}(p_*)$\n",
    "in terms of known variables, which yields\n",
    "\n",
    "\\begin{align}\n",
    "\\phi_l^S(p_*) = u_l - (p - p_l)\\sqrt{\\frac{\\alpha_l}{p_* + p_{\\infty l} + \\beta_l}} \\\\\n",
    "\\phi_r^S(p_*) = u_r + (p - p_r)\\sqrt{\\frac{\\alpha_r}{p_* + p_{\\infty r} + \\beta_r}}\n",
    "\\end{align}\n",
    "with $\\alpha_k = 2.0/((\\gamma_k + 1.0)\\rho_l)$ and $\\beta_k = (p_l + p_{\\infty l})(\\gamma_l - 1.0) / (\\gamma_l + 1.0)$. The functions $\\phi_k^S$ define the Hugoniot Loci, and they allow us to construct equation (\\ref{Phi}) and solve it using a Newton method or other root finder in order to obtain the value of $p_*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Riemann invariants for rarefaction waves\n",
    "\n",
    "Riemann invariants are variables that remain constant through simple waves such as centered rarefactions. The Riemann invariants across \n",
    "the 2-wave (contact discontinuity) are the pressure $p_*$ and the normal velocity $u_*$. The Riemann invariants for the 1-wave and 3-wave are \n",
    "the entropy and the quantities\n",
    "\n",
    "\\begin{align}\n",
    " u_l+\\frac{2c_l}{\\gamma_l-1} = u_* + \\frac{2c_{*l}}{\\gamma_l-1}, \\label{RI1} \\\\\n",
    " u_r-\\frac{2c_r}{\\gamma_r-1} = u_* - \\frac{2c_{*r}}{\\gamma_r-1}. \\label{RI2} \n",
    "\\end{align}\n",
    "The speed of sound $c_k$ with the Tammann EOS is given by\n",
    "\n",
    "\\begin{align}\n",
    " c_k= \\sqrt{\\gamma_k\\frac{p_k+p_{\\infty k}}{\\rho_k}}.  \\label{sndsp} \n",
    "\\end{align}\n",
    "\n",
    "As the entropy is invariant along the rarefaction wave, we can use the Tammann EOS isentropic relation to obtain the density in the middle states,\n",
    "\n",
    "\\begin{align}\n",
    " \\rho_{*k} = \\rho_k \\left(\\frac{\\tilde{p}_{*k}}{\\tilde{p}_k}\\right)^{1/\\gamma}.  \\label{isentrop}\n",
    "\\end{align}\n",
    "\n",
    "Solving (\\ref{RI1}) and (\\ref{RI2}) for $u_*$ and using equations (\\ref{sndsp}) and (\\ref{isentrop}), we immediately obtain\n",
    "\n",
    "\\begin{align*}\n",
    " u_* = u_l + \\frac{2c_l}{\\gamma_l-1} \\left[\n",
    " 1-\\left(\\frac{\\tilde{p}_{*l}}{\\tilde{p}_l}\\right)^{\\frac{\\gamma_l-1}{2\\gamma_l}}\n",
    " \\right] = \\phi_l^R(p_*), \\\\\n",
    " u_* = u_r - \\frac{2c_r}{\\gamma_r-1} \\left[\n",
    " 1-\\left(\\frac{\\tilde{p}_{*r}}{\\tilde{p}_r}\\right)^{\\frac{\\gamma_r-1}{2\\gamma_r}}\n",
    " \\right] = \\phi_r^R(p_*).\n",
    "\\end{align*}\n",
    "\n",
    "These functions define the integral curves, and they indicate which pairs of states can be connected by a rarefaction. Note that since $\\tilde{p}_\\kappa = p_\\kappa + p_{\\infty \\kappa}$, when we compare to equations (\\ref{CDvel-1}) and (\\ref{CDvel-2}), we obtain \n",
    "the functions $\\phi_{l,r}^R$. The rarefaction head velocities are given by $u_l-c_l$ \n",
    "and $u_r+c_r$; the tail velocities are $u_*-c_{*l}$ and $u_*+c_{*r}$.\n",
    "\n",
    "In order to compute the complete structure of the rarefaction wave <cite data-cite=\"Ivings1998\"><a href=\"riemann.html#Ivings1998\">(Ivings & Toro 1998, </a></cite> <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">LeVeque 2002)</a></cite>, we can use the \n",
    "Riemann invariants from Eqs. \\eqref{RI1} and \\eqref{RI2}, along with Eq. \\eqref{sndsp} and the isentropic \n",
    "relation from Eq. \\eqref{isentrop}. The solution for the 1-rarefaction wave along the rays $x/t=\\xi= u_{rar1} - c_{rar1}$ is then\n",
    "\n",
    "\\begin{align*}\n",
    " u_{rar1}(\\xi) &= \\frac{u_l(\\gamma_l - 1) + 2(\\xi + c_l)}{\\gamma_l + 1}, \\\\\n",
    " \\rho_{rar1}(\\xi) &= \\rho_l \\left[\\frac{u_{rar1}(\\xi) - \\xi}{c_l}\\right]^{\\frac{2}{\\gamma_l -1}}, \\\\\n",
    " p_{rar1}(\\xi) &= \\tilde{p}_l \\left[\\frac{u_{rar1}(\\xi) - \\xi}{c_l}\\right]^{\\frac{2\\gamma_l}{\\gamma_l -1}} - p_{\\infty l},\n",
    "\\end{align*}\n",
    "\n",
    "and for a 3-rarefaction wave along the rays $x/t=\\xi= u_{rar3} + c_{rar3}$ is,\n",
    "\n",
    "\\begin{align*}\n",
    " u_{rar3}(\\xi) &= \\frac{u_r(\\gamma_r - 1) + 2(\\xi - c_r)}{\\gamma_r + 1}, \\\\\n",
    " \\rho_{rar3}(\\xi) &= \\rho_r \\left[\\frac{\\xi - u_{rar3}(\\xi)}{c_r}\\right]^{\\frac{2}{\\gamma_r -1}}, \\\\\n",
    " p_{rar3}(\\xi) &= \\tilde{p}_r \\left[\\frac{\\xi - u_{rar3}(\\xi)}{c_r}\\right]^{\\frac{2\\gamma_r}{\\gamma_r -1}} - p_{\\infty r}.\n",
    "\\end{align*}\n",
    "\n",
    "Now that we know the integral curves given by $\\phi^{R}_{l,r}$ (for the rarefactions), we can construct the function $\\Phi(p_*)$ from (\\ref{Phi}) for any of the 4 possible scenarios. \n",
    "\n",
    "\\begin{align*}\n",
    "\\Phi(p_*) & = \\begin{cases}\n",
    " \\phi_r^R(p_*) - \\phi_l^R(p_*) & \\text{if } p_* < p_l \\text{ and } p_*<p_r \\\\\n",
    " \\phi_r^R(p_*) - \\phi_l^S(p_*) & \\text{if } p_l \\le p_* \\le p_r \\\\\n",
    " \\phi_r^S(p_*) - \\phi_l^R(p_*) & \\text{if } p_r \\le p_* \\le p_l \\\\\n",
    " \\phi_r^S(p_*) - \\phi_l^S(p_*) & \\text{if } p_* > p_l \\text{ and } p_* > p_r\n",
    "              \\end{cases}\n",
    "\\end{align*}\n",
    "\n",
    "The value of $p_*$ will be \n",
    "found by numerically finding the root of $\\Phi(p_*)=0$. Once $p_*$ is found, $u_*$, $\\rho_{*l}$, $\\rho_{*r}$, $s_1$ \n",
    "and $s_3$ can be found using the relations above. As we \n",
    "know the three wave speeds $s_1$, $s_2$ and $s_3$ and the values of the primitive variables $[\\rho,u,p]$ in all the 4 states, we have solved the Riemann problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution in the pressure-velocity plane\n",
    "\n",
    "As we did in [the chapter on ideal gas flow](Euler_equations.ipynb), here we visualize the solution in the plane whose axes are pressure and velocity. We have a left and a right state such that each one is connected to some middle state by a shock or a rarefaction, depending on the entropy conditions. The connection is established by curves in the phase plane: Hugoniot loci in the case of shocks and integral curves in the case of rarefactions. The intersection between these curves yields the middle state, which is the most important element of the solution to the Riemann problem. Let's look at an interactive plot illustrating the solution in the phase plane. The initial states and the middle state are indicated as black points; the Hugoniot loci are plotted in red and the integral curves in blue. The physical solution satisfying the entropy conditions is plotted as a continous line, whereas the unphysical solution is shown as a dashed line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets\n",
    "import matplotlib.pyplot as plt\n",
    "from exact_solvers import euler_tammann, interactive_pplanes\n",
    "from utils import riemann_tools\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial states [density, velocity, pressure]\n",
    "ql = [600.0, 10.0, 50000.0]\n",
    "qr = [50.0, -10.0, 25000.0]\n",
    "# Tammann EOS parameters [gamma, pinf]\n",
    "paramsl = [1.4, 0.0]\n",
    "paramsr = [7.0, 100.0]\n",
    "\n",
    "interactive_pplanes.euler_tammann_interactive_phase_plane(ql,qr,paramsl,paramsr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full solution\n",
    "In order to get the full solution, we need some additional quantities. Equations (\\ref{RH-phis}) will yield the contact discontinuity speed $s_2=u_*$ in terms of $p_*$. We can also calculate $F_l$ and $F_r$ from (\\ref{Qk-f}), and we can substitute in (\\ref{RH-speeds}) to obtain the corresponding wave speeds. With this, we can provide the full solution of the Riemann problem. In the next section, we will show the full solution for different examples.   we define a function `plot_riemann_solution` to plot the exact solution of the Riemann problem and parameters that we will use for the different examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "def plot_riemann_solution(q_l, q_r, auxl, auxr):\n",
    "    ex_states, ex_speeds, reval, wave_types, varsout = \\\n",
    "        euler_tammann.exact_riemann_solution(q_l ,q_r, auxl, auxr, \n",
    "                                             varin='primitive',\n",
    "                                             varout='primitive')\n",
    "    varnames = euler_tammann.primitive_variables\n",
    "    plot_function = \\\n",
    "            riemann_tools.make_plot_function(ex_states, ex_speeds, \n",
    "                                             reval, wave_types, \n",
    "                                             layout='vertical', \n",
    "                                             variable_names=varnames,\n",
    "                                             aux=(np.array(auxl), np.array(auxr)),\n",
    "                                             plot_chars=[euler_tammann.lambda1,\n",
    "                                                         euler_tammann.lambda2,\n",
    "                                                         euler_tammann.lambda3],\n",
    "                                            contact_index=1)\n",
    "        \n",
    "    dropdown = widgets.Dropdown(options=[None,1,2,3],\n",
    "                                description='Show characteristics')\n",
    "    return interact(plot_function, \n",
    "                    t=widgets.FloatSlider(value=0.0,min=0,max=1.0), \n",
    "                    which_char=dropdown);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples of Riemann Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1: Sod shock tube in water\n",
    "We first consider the classic Sod shock tube problem with higher density and pressure on the left state than on the right state and both sides initially at rest. The main difference with the classic problem is that in this case the shock tube is filled with water and not with air. The propagation of waves in water can be accurately modeled using $\\gamma = 7.15$, $p_{\\infty} = 3\\times 10^8$ and $\\rho\\approx 1000 kg/m^3$; see <cite data-cite=\"del2016computational\"><a href=\"riemann.html#del2016computational\">(del Razo 2016)</a></cite>. A value of $p_{atm}=101325 Pa$ is used for Atmospheric pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gamma_water = 7.15\n",
    "gamma_air = 1.4\n",
    "pinf_water = 3.e8\n",
    "pinf_air = 0.\n",
    "patm = 101325.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We take a very large jump in pressure so that it's clear the left-going wave is a rarefaction. With a more modest jump the rarefaction wave barely spreads and appears more like a discontinuity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = [1500.0, 0.0, 3000*patm]\n",
    "q_r = [1000.0, 0.0, patm]\n",
    "gamma = [gamma_water, gamma_water]\n",
    "pinf = [pinf_water, pinf_water]\n",
    "\n",
    "plot_riemann_solution(q_l, q_r, [gamma[0],pinf[0]], [gamma[1],pinf[1]]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the small density difference accross left-going wave and the very compressed rarefaction fan\n",
    "are consequence of the low compressibility of water."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2: Transmission and reflection at an air-water interface\n",
    "We can use the Euler equations with the Tammann EOS to model wave propagation in almost-incompressible media, such as water. Since we have the solution for different EOS parameters in the left and right states, we can model interaction of different fluids across the interface. For instance, we can use this model to study the Riemann problem across an air-water interface. Air is modeled by the ideal gas EOS, which corresponds to $\\gamma=1.4$ and $p_{\\infty} = 0$, and water can be modeled using the same parameters as the example above. \n",
    "\n",
    "Consider the problem of a shock traveling from left to right through the air and hitting the interface at time 0. Then this wave should be partly transmitted as a shock wave moving to the right in the water, and partly reflected as a lef-going wave in the air.  If we suppose the shock wave hits the interface at time $t=0$, then we can calculate the transmitted and reflected waves by solving the Riemann problem starting at $t=0$ for a case with a positive pressure and velocity in the left state and ambient water conditions in the right state.\n",
    "\n",
    "As we want to model an air-water interface, the left state has only air while the right state has only water with their corresponding densities. Although in a realistic two or three dimensional setting, the water would spill into the left state, we can suppose that the water is being held in place by a container with a thin wall. In this setting, the current model has been shown to be a valid approximation <cite data-cite=\"delrazo2017\"><a href=\"riemann.html#delrazo2017\">(del Razo 2017)</a></cite>. The solution has the following form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = [1.0, 350.0, 300*patm] # Initial condition for shock in air\n",
    "q_r = [1000.0, 0.0, patm] # Initial condition in still water\n",
    "gamma = [gamma_air, gamma_water]\n",
    "pinf = [pinf_air, pinf_water]\n",
    "\n",
    "plot_riemann_solution(q_l, q_r, [gamma[0],pinf[0]], [gamma[1],pinf[1]]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note we exaggerated again the initial jump in pressure, so the jump in the density and velocity of the right-going wave could be appreciated in the plots; otherwise, it seems there is no jump in the right-going wave. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3: Water-air interface\n",
    "Here we consider a setup similar to that of Problem 2, but we invert the materials. The left region has water, while the right region has air."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = [1000.0, 350.0, 300*patm] # Initial condition for shock in air\n",
    "q_r = [1.0, 0.0, patm] # Initial condition in still water\n",
    "gamma = [gamma_water, gamma_air]\n",
    "pinf = [pinf_water, pinf_air]\n",
    "\n",
    "plot_riemann_solution(q_l, q_r, [gamma[0],pinf[0]], [gamma[1],pinf[1]]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we exaggerated the jump in pressure to be able to observe the jumps of the solution. Nonetheless, the density jump in the right-going wave is still not perceptible. A more exaggerated initial condition can make this jump more evident."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Strong rarefaction waves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we mentioned when we introduced the Tammann EOS, the relevant physical quantities remain well-defined for \n",
    "negative values of $p$, as long as $p+p_{\\infty}>0$. For instance, the densities \\eqref{rho-k}, \\eqref{isentrop} and the sound speed \\eqref{sndsp} are all well defined in this case. This is consistent with the interpretation of the Tammann EOS as a shift of the zero pressure point of an ideal gas.\n",
    "\n",
    "As an example of this scenario, we consider a strong symmetric expansion wave in water. Analogous to the Euler equations using an ideal gas EOS, we consider equal densities and pressures, and equal and opposite velocities, with the initial states moving away from each other. The result is two rarefaction waves (the contact has zero strength). Notice that the value in the pressure plot is negative; it reaches a value of about $-2.8\\times 10^8$, so $p+p_{\\infty}$ is still positive.  Of course, it is questionable whether the model remains accurate when $p+p_{\\infty}$ is nearly zero, as it is in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity = 350\n",
    "q_l = [1000.0, -velocity, 2*patm]\n",
    "q_r = [1000.0, velocity, 2*patm]\n",
    "auxl = [gamma_water, pinf_water]\n",
    "auxr = [gamma_water, pinf_water]\n",
    "ex_states, ex_speeds, reval, wave_types, varsout = \\\n",
    "        euler_tammann.exact_riemann_solution(q_l ,q_r, auxl, auxr, \n",
    "                                             varin='primitive',\n",
    "                                             varout='primitive')\n",
    "    \n",
    "print('Middle-state pressure p = %10.3e' % ex_states[2,1])\n",
    "print('        p_inf for water = %10.3e,' % pinf_water)\n",
    "print('                p+p_inf = %10.3e' % (ex_states[2,1]+pinf_water))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_riemann_solution(q_l, q_r, auxl, auxr);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}