<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Riemann problem for Euler equations with the Tammann EOS\n",
    "We consider the one dimensional Riemann problem for the Euler equations with the Tammann 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",
    "\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",
    "\n",
    "\n",
    "where $\\rho$ is density, $u$ the velocity, $E$ the internal energy, and $p$ the pressure; the subcripts $x,t$ denote partial derivatives. The Tamann EOS is given by \n",
    "\n",
    "\n",
    "\\begin{align*}\n",
    "p=\\rho e(\\gamma -1) - \\gamma p_{\\infty},\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "where $e$ is the specific internal energy.  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]$. 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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Exact solution\n",
    "\n",
    "The exact solution to the Riemann problem with a Tammann EOS can be found in <cite data-cite=\"Ivings1998\"><a href=\"riemann.html#Ivings1998\">(Ivings & Toro 1998)<a></cite>. In the supplemental material of the work <cite data-cite=\"delrazo2017\"><a href=\"riemann.html#delrazo2017\">(del Razo, 2017)<a></cite>, the reader can find an extension of the analytic solution to the case where the parameters $\\gamma$ and $p_\\infty$ in the Tammann EOS have a discontinuity. In this chapter we will focus in the latter, since it is a more general case. We will reproduce the derivation of the exact solution while at the same time showing some sricpts and plots to illustrate the structure of the solution. \n",
    "\n",
    "As we know from the previous examples, the solution of the Euler equations will consist 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 system will have four different solution states, $q_l,q_{*l},q_{*r},q_r$ separated by the three waves. In order to figure out if 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. As we know the velocity on the left state $u_l$ should be connected by a rarefaction or shock to $u_*$, we can \n",
    "calculate $u_*=u_l + [u]_1$ with $[u]_1$ 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 can calculate $u_*= u_r - [u]_3$; therefore, it is useful to define\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "  \\phi_l(p_*) = u_* = u_l - \\mathcal{F}_l(p_*),   \\label{CDvel-1}  \\\\\n",
    "  \\phi_r(p_*) = u_* = u_r + \\mathcal{F}_r(p_*).   \\label{CDvel-2} \n",
    "\\end{align}\n",
    "\n",
    "\n",
    "where $\\mathcal{F}_{l,r}(p_*) = -[u]_{1,3}$ will change form depending if it's a shock or a rarefaction (signs were \n",
    "chosen for notation consistency). As we expect these two equations yield the same contact discontinuity \n",
    "velocity $u_*$, then \n",
    "\n",
    "\n",
    "\\begin{align}\n",
    " \\Phi(p_*)= \\phi_r(p_*) - \\phi_l(p_*) = 0.\n",
    " \\label{Phi} \n",
    "\\end{align}\n",
    "\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_*$. How can we know which one, can be answered by simple physical \n",
    "intuition. If the pressure is higher on the side toward which the wave is propagating, it will yield a rarefaction. If the pressure is lower, it will 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 if the $\\phi$ was obtained by using the Rankine-Hugoniot equations to connect \n",
    "states by shocks or the Riemann invariants to connect them by rarefactions respectively.\n",
    "\n",
    "We will 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 will also show how to obtain the density and the missing wave speeds. We denote the speed of the 1-wave, $s_1$, the 2-wave, $s_2$, and the 3-wave $s_3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Rankine-Hugoniot conditions for shock waves\n",
    "As we know the 1-wave and the 3-wave could each be a shock. In that case, the velocity of the wave, i.e. the \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 writing $s_n$, with $n=1,3$.\n",
    "\n",
    "The Rankine-Hugoniot conditions are in general given by $s_n\\left(q_k - q_{*k}\\right) = f(q_k) - f(q_{*k})$, \n",
    "where $q$ is the vector state variable, $f(q)$ the vector state flux and $s_n$ the shock velocity. For \n",
    "the Euler equations this can be easily rewritten as (<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",
    "\n",
    "where $k=l,r$, $\\omega_k = u_k - s_n$ , $\\omega_* = u_* - s_n$, with $n=1$ when $k=l$ and $n=3$ when $k=r$, 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 can start by defining \n",
    "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 equation (\\ref{Ql}) and (\\ref{Qr}) \n",
    "into (\\ref{RH-mom}) to immediately 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",
    "\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",
    "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",
    "\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",
    "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\n",
    "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",
    "Replacing 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 $\\phi^S_{l,r}$ nonlinear functions of $p_*$\n",
    "in terms of known variables, which yields\n",
    "\n",
    "\\begin{align}\n",
    "\\phi_l^S(p_*) = u_l - (p - pl)\\sqrt{\\frac{\\alpha_l}{p_* + p_{\\infty l} + \\beta_l}} \\\\\n",
    "\\phi_r^S(p_*) = u_r + (p - pr)\\sqrt{\\frac{\\alpha_r}{p_* + p_{\\infty r} + \\beta_r}}\n",
    "\\end{align}\n",
    "\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 $\\phi_k^S$ functions are called 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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "\n",
    "### Riemann invariants for rarefaction waves\n",
    "\n",
    "Riemann invariants are variables that remain constant through simple waves such as rarefactions. The Riemann invariants across \n",
    "the 2-wave 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",
    "\n",
    "correspondingly. 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, 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 are called the integral curves, and they tell you which values are acceptable to connect two states with a rarefaction curve. Note that as $\\tilde{p}_\\kappa = p_\\kappa + p_{\\infty \\kappa}$, when we compare to equations (\\ref{CDvel-1} and \\ref{CDvel-2}), we obtain \n",
    "the $\\phi_{l,r}^R$ functions. The rarefaction head velocities will be given by $u_l-c_l$ \n",
    "and $u_r+c_r$; the tail velocities will be $u_*-c_{*l}$ and $u_*+c_{*r}$. For numerical purposes, a simple \n",
    "approximate velocity is provided for $s_1$ and $s_3$ as the average between the head and tail velocity. \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. \\ref{RI1} and \\ref{RI2}, along with Eq. \\ref{sndsp} and the isentropic \n",
    "relation from Eq. \\ref{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^{s,r}_{l,r}$ (for the rarefactions), we can construct the function $\\Phi(p_*)$ from (\\ref{Phi}) for any of the 4 possible scenarios. The value of $p_*$ will be \n",
    "found by numerically finding the roots of $\\Phi(p_*)=0$. Note which case to employ to calculate $\\Phi(p_*)$ \n",
    "might change in each iteration of the root finder. 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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Solution in the pressure-velocity plane\n",
    "\n",
    "The easiest way to visualize the solution is in the phase 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. Lets 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. Note that, in this case, the phase plane is three dimensional, but it is easierly understood when projecting into a two dimensional phase space. Nonetheless, the variable $\\rho$ in the flattened dimension can still be modified as an additional paramater of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interact, widgets\n",
    "import matplotlib.pyplot as plt\n",
    "from exact_solvers import euler_tammann, interactive_pplanes\n",
    "from utils import riemann_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Full solution\n",
    "In order to get the full solution, we still need a couple more 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. It should be noted that this model is to be handled carefully because it could yield unphysical results (negative pressures) for some initial conditions. The cases where this happens are in the limit of validity for modeling almost incompressible media with the Tammann EOS. We will explore these limits in Problem 4. Here we define a function 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
   },
   "outputs": [],
   "source": [
    "def plot_riemann_solution(q_l, q_r, gamma, pinf):\n",
    "    ex_states, ex_speeds, reval, wave_types, varsout = \\\n",
    "        euler_tammann.exact_riemann_solution(q_l ,q_r, gamma, pinf, \n",
    "                                             varin='primitive', varout='primitive')\n",
    "\n",
    "    plot_function = riemann_tools.make_plot_function(ex_states, ex_speeds, \n",
    "                                                     reval, wave_types, \n",
    "                                                     layout='vertical', \n",
    "                                                     variable_names=euler_tammann.primitive_variables)\n",
    "\n",
    "    return interact(plot_function, t=widgets.FloatSlider(value=0.0,min=0,max=1.0));\n",
    "\n",
    "gamma_water = 7.15\n",
    "gamma_air = 1.4\n",
    "pinf_water = 3.e8\n",
    "pinf_air = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Examples of Riemann Solutions\n",
    "\n",
    "### Problem 1: Shock tube in water\n",
    "We first consider a shock tube filled with water, and a shock traveling from left two right. 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>. Atmospheric pressure is modeled at $p_{atm}=101325 Pa$ and speed in $m/s$. The density difference between the states is very small because water is almost incompressible. Note the small waves generated in the density plot and the very compressed rarefaction fan, which are consequence of the low compressibility in water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1010.0, 0.0, 3*101325.0]\n",
    "q_r = [1000.0, 0.0, 101325.0]\n",
    "gamma = [gamma_water, gamma_water]\n",
    "pinf = [pinf_water, pinf_water]\n",
    "\n",
    "plot_riemann_solution(q_l, q_r, gamma, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Problem 2: Air-water interface in shock tube\n",
    "We can use the Euler equations with the Tammann EOS to model wave propagation in almost-incompressible media, such as water. As we provided the solution with different EOS parameters on the left and right states, it can couple different materials across the interface. For instance, we can use this model to study the propagation of waves from air into water. 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 before. Consider a shock traveling from left to right; the left state has only air while the right state has only water. Although in a realistic two or three dimensional setting, the water would spill into the left state, we can consider 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=\"del2016computational0\"><a href=\"riemann.html#del2016computational0\">(del Razo 2016)<a></cite>. The solution has the following form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1.0, 350.0, 200000.0] # Initial condition for shock in air\n",
    "q_r = [1000.0, 0.0, 101325.0] # 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, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Problem 3: Water-air interface in shock tube\n",
    "The same setup as in Problem 2, but now the shock moves from water (on the left) into air (on the right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1000.0, 350.0, 200000.0] # Initial condition for shock in air\n",
    "q_r = [1.0, 0.0, 101325.0] # 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, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Problem 4: Strong symmetric expansion in water (unphysical)\n",
    "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). However, notice the value in the pressure plot is several orders of magnitude below zero. This is an unphysical result; however, it is mathematically correct. The main problem is that we have not yet established how our model can deal wih vacuum states. It is important to be aware that vacuum states can be observed frequently when using the Tammann EOS to model almost incompressible materials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1000.0, -350.0, 2*101325.0] # Initial condition for shock in air\n",
    "q_r = [1000.0, 350.0, 2*101325.0] # Initial condition in still water\n",
    "gamma = [gamma_water, gamma_water]\n",
    "pinf = [pinf_water, pinf_water]\n",
    "\n",
    "plot_riemann_solution(q_l, q_r, gamma, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true<commit_msg>Added section in negative pressure<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Riemann problem for Euler equations with the Tammann EOS\n",
    "We consider the one dimensional Riemann problem for the Euler equations with the Tammann 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",
    "\n",
    "where $\\rho$ is density, $u$ the velocity, $E$ the internal energy, and $p$ the pressure; 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",
    "\n",
    "where $e$ is the specific internal energy. For a general equation of state $p (\\rho,e)$, the sound of speed is given by the isentropic derivative (keeping entropy constant) of the pressure,\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",
    "\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",
    "\n",
    "Note this is almost the same as the sound speed in an ideal gas 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. 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.\n",
    "\n",
    "Using the equation of state, the state of the system can also be written in the primitive variables $[\\rho, u, p]$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Exact solution\n",
    "\n",
    "The exact solution to the Riemann problem with a Tammann EOS can be found in <cite data-cite=\"Ivings1998\"><a href=\"riemann.html#Ivings1998\">(Ivings & Toro 1998)<a></cite>. In the supplemental material of the work <cite data-cite=\"delrazo2017\"><a href=\"riemann.html#delrazo2017\">(del Razo, 2017)<a></cite>, the reader can find an extension of the analytic solution to the case where the parameters $\\gamma$ and $p_\\infty$ in the Tammann EOS have a discontinuity. In this chapter we will focus in the latter, since it is a more general case. We will reproduce the derivation of the exact solution while at the same time showing some sricpts and plots to illustrate the structure of the solution. \n",
    "\n",
    "As we know from the previous examples, the solution of the Euler equations will consist 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 system will have four different solution states, $q_l,q_{*l},q_{*r},q_r$ separated by the three waves. In order to figure out if 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. As we know the velocity on the left state $u_l$ should be connected by a rarefaction or shock to $u_*$, we can \n",
    "calculate $u_*=u_l + [u]_1$ with $[u]_1$ 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 can calculate $u_*= u_r - [u]_3$; therefore, it is useful to define\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "  \\phi_l(p_*) = u_* = u_l - \\mathcal{F}_l(p_*),   \\label{CDvel-1}  \\\\\n",
    "  \\phi_r(p_*) = u_* = u_r + \\mathcal{F}_r(p_*).   \\label{CDvel-2} \n",
    "\\end{align}\n",
    "\n",
    "\n",
    "where $\\mathcal{F}_{l,r}(p_*) = -[u]_{1,3}$ will change form depending if it's a shock or a rarefaction (signs were \n",
    "chosen for notation consistency). As we expect these two equations yield the same contact discontinuity \n",
    "velocity $u_*$, then \n",
    "\n",
    "\n",
    "\\begin{align}\n",
    " \\Phi(p_*)= \\phi_r(p_*) - \\phi_l(p_*) = 0.\n",
    " \\label{Phi} \n",
    "\\end{align}\n",
    "\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_*$. How can we know which one, can be answered by simple physical \n",
    "intuition. If the pressure is higher on the side toward which the wave is propagating, it will yield a rarefaction. If the pressure is lower, it will 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 if the $\\phi$ was obtained by using the Rankine-Hugoniot equations to connect \n",
    "states by shocks or the Riemann invariants to connect them by rarefactions respectively.\n",
    "\n",
    "We will 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 will also show how to obtain the density and the missing wave speeds. We denote the speed of the 1-wave, $s_1$, the 2-wave, $s_2$, and the 3-wave $s_3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Rankine-Hugoniot conditions for shock waves\n",
    "As we know the 1-wave and the 3-wave could each be a shock. In that case, the velocity of the wave, i.e. the \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 writing $s_n$, with $n=1,3$.\n",
    "\n",
    "The Rankine-Hugoniot conditions are in general given by $s_n\\left(q_k - q_{*k}\\right) = f(q_k) - f(q_{*k})$, \n",
    "where $q$ is the vector state variable, $f(q)$ the vector state flux and $s_n$ the shock velocity. For \n",
    "the Euler equations this can be easily rewritten as (<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",
    "\n",
    "where $k=l,r$, $\\omega_k = u_k - s_n$ , $\\omega_* = u_* - s_n$, with $n=1$ when $k=l$ and $n=3$ when $k=r$, 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 can start by defining \n",
    "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 equation (\\ref{Ql}) and (\\ref{Qr}) \n",
    "into (\\ref{RH-mom}) to immediately 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",
    "\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",
    "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",
    "\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",
    "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\n",
    "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",
    "Replacing 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 $\\phi^S_{l,r}$ nonlinear functions of $p_*$\n",
    "in terms of known variables, which yields\n",
    "\n",
    "\\begin{align}\n",
    "\\phi_l^S(p_*) = u_l - (p - pl)\\sqrt{\\frac{\\alpha_l}{p_* + p_{\\infty l} + \\beta_l}} \\\\\n",
    "\\phi_r^S(p_*) = u_r + (p - pr)\\sqrt{\\frac{\\alpha_r}{p_* + p_{\\infty r} + \\beta_r}}\n",
    "\\end{align}\n",
    "\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 $\\phi_k^S$ functions are called 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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "\n",
    "### Riemann invariants for rarefaction waves\n",
    "\n",
    "Riemann invariants are variables that remain constant through simple waves such as rarefactions. The Riemann invariants across \n",
    "the 2-wave 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",
    "\n",
    "correspondingly. 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, 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 are called the integral curves, and they tell you which values are acceptable to connect two states with a rarefaction curve. Note that as $\\tilde{p}_\\kappa = p_\\kappa + p_{\\infty \\kappa}$, when we compare to equations (\\ref{CDvel-1} and \\ref{CDvel-2}), we obtain \n",
    "the $\\phi_{l,r}^R$ functions. The rarefaction head velocities will be given by $u_l-c_l$ \n",
    "and $u_r+c_r$; the tail velocities will be $u_*-c_{*l}$ and $u_*+c_{*r}$. For numerical purposes, a simple \n",
    "approximate velocity is provided for $s_1$ and $s_3$ as the average between the head and tail velocity. \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. \\ref{RI1} and \\ref{RI2}, along with Eq. \\ref{sndsp} and the isentropic \n",
    "relation from Eq. \\ref{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^{s,r}_{l,r}$ (for the rarefactions), we can construct the function $\\Phi(p_*)$ from (\\ref{Phi}) for any of the 4 possible scenarios. The value of $p_*$ will be \n",
    "found by numerically finding the roots of $\\Phi(p_*)=0$. Note which case to employ to calculate $\\Phi(p_*)$ \n",
    "might change in each iteration of the root finder. 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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Solution in the pressure-velocity plane\n",
    "\n",
    "The easiest way to visualize the solution is in the phase 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. Lets 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. Note that, in this case, the phase plane is three dimensional, but it is easierly understood when projecting into a two dimensional phase space. Nonetheless, the variable $\\rho$ in the flattened dimension can still be modified as an additional paramater of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interact, widgets\n",
    "import matplotlib.pyplot as plt\n",
    "from exact_solvers import euler_tammann, interactive_pplanes\n",
    "from utils import riemann_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Full solution\n",
    "In order to get the full solution, we still need a couple more 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. It should be noted that this model is to be handled carefully because it could yield unphysical results (negative pressures) for some initial conditions. The cases where this happens are in the limit of validity for modeling almost incompressible media with the Tammann EOS. We will explore these limits in Problem 4. Here we define a function 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,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def plot_riemann_solution(q_l, q_r, gamma, pinf):\n",
    "    ex_states, ex_speeds, reval, wave_types, varsout = \\\n",
    "        euler_tammann.exact_riemann_solution(q_l ,q_r, gamma, pinf, \n",
    "                                             varin='primitive', varout='primitive')\n",
    "\n",
    "    plot_function = riemann_tools.make_plot_function(ex_states, ex_speeds, \n",
    "                                                     reval, wave_types, \n",
    "                                                     layout='vertical', \n",
    "                                                     variable_names=euler_tammann.primitive_variables)\n",
    "\n",
    "    return interact(plot_function, t=widgets.FloatSlider(value=0.0,min=0,max=1.0));\n",
    "\n",
    "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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Examples of Riemann Solutions\n",
    "\n",
    "### Problem 1: Shock tube in water\n",
    "We first consider a shock tube filled with water, and a shock traveling from left two right. 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>. Atmospheric pressure is modeled at $p_{atm}=101325 Pa$ and speed in $m/s$. The density difference between the states is very small because water is almost incompressible. Note the small waves generated in the density plot and the very compressed rarefaction fan, which are consequence of the low compressibility in water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1010.0, 0.0, 3*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, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Problem 2: Air-water interface in shock tube\n",
    "We can use the Euler equations with the Tammann EOS to model wave propagation in almost-incompressible media, such as water. As we provided the solution with different EOS parameters on the left and right states, it can couple different materials across the interface. For instance, we can use this model to study the propagation of waves from air into water. 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 before. Consider a shock traveling from left to right; the left state has only air while the right state has only water. Although in a realistic two or three dimensional setting, the water would spill into the left state, we can consider 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=\"del2016computational0\"><a href=\"riemann.html#del2016computational0\">(del Razo 2016)<a></cite>. The solution has the following form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1.0, 350.0, 2*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, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Problem 3: Water-air interface in shock tube\n",
    "The same setup as in Problem 2, but now the shock moves from water (on the left) into air (on the right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q_l = [1000.0, 350.0, 2*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, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Negative pressure and vacuum states\n",
    "\n",
    "As we mentioned when we introduced the Tammann EOS, the relvant physical quantities remain well defined for \n",
    "negative pressures, as long as $p+p_{\\infty}>0$. For instance, the densities \\ref{rho-k}, \\ref{isentrop} and the sound speed \\ref{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. Nonetheless, negative pressures remain unphysical; therefore, if $p<0$ while $p+p_{\\infty}>0$, one should determine on a case by case basis if the solution to the Riemann problem is physically relevant or not. In order to determine this, one should ask if the Tammann EOS is an accurate model for your specific problem and what are your variables of interest since pressure is not physically accurate when approaching the negative regime. \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). However, notice the value in the pressure plot is negative. This is exactly the case where $p<0$ and $p+p_{\\infty}>0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "velocity = 350\n",
    "q_l = [1000.0, -velocity, 2*patm] # Initial condition for shock in air\n",
    "q_r = [1000.0, velocity, 2*patm] # Initial condition in still water\n",
    "gamma = [gamma_water, gamma_water]\n",
    "pinf = [pinf_water, pinf_water]\n",
    "\n",
    "plot_riemann_solution(q_l, q_r, gamma, pinf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It may be possible to extend the solution to cases where the Riemann solver yields $p+p_{\\infty} \\le 0$. As in the Euler equations chapter, the vacuum states arise in regions with zero density. As vacuum states can only connected through rarefaction waves, the middle state densities can only be given by \\ref{isentrop}. The middle state density will only yield zero when $\\tilde{p_*} = p_*+p_{\\infty}=0$ or $p_*=-p_{\\infty}$. The convergence failure of the Riemann solution for this scenario can be observed if the velocity in the example above is increased one order of magnitude. How to deal with vacuum states when using a Tammann EOS with a jump in the parameters is beyond the scope of this book."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
{'fraction_non_alphanumeric': 0.1883184473843261, 'fraction_numerical': 0.01418642132831912, 'mean_word_length': 2.8588364861363047, 'pattern_counts': {'":': 223, '<': 65, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '103468', 'n_tokens_mistral': 21958, 'n_tokens_neox': 20845, 'n_words': 6832}