{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Acoustics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we consider our first *system* of hyperbolic conservation laws.  We study the acoustics equations that were introduced briefly in [Introduction.ipynb](Introduction.ipynb).  We first describe the physical context of this system and then investigate its characteristic structure and the solution to the Riemann problem.  This system is described in more detail in Chapter 3 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque 2002)</a></cite>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Physical setting\n",
    "The linear acoustic equations describe the propagation of small perturbations in a fluid, such as sound waves.  In [Advection.ipynb](Advection.ipynb) we derived the one-dimensional continuity equation, which describes mass conservation:  \n",
    "\\begin{align} \\label{Ac:continuity}\n",
    "    \\rho_t + (\\rho u)_x & = 0.\n",
    "\\end{align}  \n",
    "For more realistic fluid models, we need another equation that determines the velocity $u$.  This typically takes the form of a conservation law for momentum $\\rho u$.  Momentum, like density, is transported by fluid motion with corresponding flux $\\rho u$.  Additionally, any difference in pressure will also lead to a flux of momentum that is proportional to the pressure difference.  Thus the momentum equation takes the form  \n",
    "\\begin{align} \\label{Ac:mom_cons}\n",
    "(\\rho u)_t + (\\rho u^2 + P(\\rho))_x & = 0.\n",
    "\\end{align}  \n",
    "Here we have assumed that the pressure $P$ is a function of the density.  We will consider fully nonlinear fluid motions with more general equation of state in [Euler_equations.ipynb](Euler_equations.ipynb) and [Euler_equations_TammannEOS.ipynb](Euler_equations_TammannEOS.ipynb).  For now we investigate the behavior of small perturbations in the system above.\n",
    "\n",
    "Together, \\eqref{Ac:continuity} and \\eqref{Ac:mom_cons} form a hyperbolic system $q_t+f(q)_x=0$ with  \n",
    "\\begin{align*}\n",
    "q & = \\begin{bmatrix} \\rho \\\\ \\rho u \\end{bmatrix} & \n",
    "f(q) & = \\begin{bmatrix} \\rho u \\\\ \\rho u^2 + P(\\rho) \\end{bmatrix}\n",
    "\\end{align*}  \n",
    "Throughout this book we will make use of the quasilinear form of a given hyperbolic system:\n",
    "$$q_t + f'(q) q_x = 0.$$  \n",
    "Here $f'(q)$ denotes the Jacobian of the flux $f$ with respect to the conserved variables $q$.  In the present system, as is often the case, $f$ is most naturally written in terms of so-called primitive variables (in this case $\\rho$ and $u$) rather than in terms of the conserved variables $q$.  In order to find the flux Jacobian (and thus the quasilinear form), we first write $f$ in terms of $(q_1,q_2) = (\\rho, \\rho u)$:  \n",
    "\\begin{align}\n",
    "f(q) & = \\begin{bmatrix} q_2 \\\\ q_2^2/q_1 + P(q_1) \\end{bmatrix}.\n",
    "\\end{align}  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can differentiate to find the flux Jacobian:  \n",
    "\\begin{align*}\n",
    "f'(q) & = \\begin{bmatrix} \\partial f_1/\\partial q_1 & \\partial f_1/\\partial q_2 \\\\\n",
    "                          \\partial f_2/\\partial q_1 & \\partial f_2/\\partial q_2 \\end{bmatrix} \\\\\n",
    "      & = \\begin{bmatrix} 0 & 1 \\\\ -q_2^2/q_1^2 + P'(q_1) & 2q_2/q_1 \\end{bmatrix} \\\\\n",
    "      & = \\begin{bmatrix} 0 & 1 \\\\ P'(\\rho)-u^2 & 2u \\end{bmatrix}.\n",
    "\\end{align*}\n",
    "\n",
    "Thus small perturbations to an ambient fluid state $\\rho_0, u_0$ evolve according to the linearized equations $q_t + f'(q_0) q_x = 0$, or  \n",
    "\\begin{align*}\n",
    "\\rho_t + (\\rho u)_x & = 0 \\\\\n",
    "(\\rho u)_t + (P'(\\rho_0)-u_0^2)\\rho_x + 2u_0(\\rho u)_x & = 0.\n",
    "\\end{align*}  \n",
    "Continuing with the assumption that $\\rho-\\rho_0$ and $\\rho u - \\rho_0 u_0$ are proportional to a small parameter $\\epsilon$, and discarding terms of order $\\epsilon^2$ and higher, one obtains the linear hyperbolic system  \n",
    "\\begin{align*}\n",
    "p_t + u_0 p_x + P'(\\rho_0) u_x & = 0 \\\\\n",
    "u_t + \\frac{1}{\\rho_0} p_x + u_0 u_x & = 0,\n",
    "\\end{align*}\n",
    "where $p(x,t)$ is the pressure as a function of $x$ and $t$.\n",
    "*Do we want to do anything with this more general system?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the ambient fluid is at rest (i.e. $u_0=0$) and the pressure is directly proportional to the density, then this simplifies to\n",
    "\\begin{align} \\label{Ac:main}\n",
    " \\left[ \\begin{array}{c}\n",
    "p \\\\\n",
    "u \n",
    "\\end{array} \\right]_t\n",
    "+  \\underbrace{\\left[ \\begin{array}{cc}\n",
    "0 & K_0 \\\\\n",
    "1/\\rho_0 & 0  \\\\\n",
    "\\end{array} \\right]}_{\\mathbf{A}}\n",
    "\\left[ \\begin{array}{c}\n",
    "p \\\\\n",
    "u \\end{array} \\right]_x = 0,\n",
    "\\end{align}\n",
    "where $K_0=\\rho_0 P'(\\rho_0)$ is referred to as the bulk modulus of compressibility."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the rest of this chapter we work with \\eqref{Ac:main} and let $q=[p,u]^T$.  Then we can write \\eqref{Ac:main} as $q_t + A q_x = 0$.  For simplicity, we also drop the subscripts on $K, \\rho$.  Direct calculation reveals that the eigenvectors of $A$ are\n",
    "\\begin{align}\n",
    "\\lambda_{1,2} & = \\pm \\sqrt{{K}/{\\rho}} = \\pm c,\n",
    "\\end{align}\n",
    "while the right eigenvectors are\n",
    "\\begin{align*}\n",
    "r_1 = \\begin{bmatrix}\\begin{array}{c}-\\sqrt{K\\rho}\\\\1\\end{array}\\end{bmatrix}, \\qquad r_2 = \\begin{bmatrix}\\begin{array}{c}\\sqrt{K\\rho}\\\\1\\end{array}\\end{bmatrix}.\n",
    "\\end{align*}\n",
    "\n",
    "Defining $R = [r_1 r_2]$ and $\\Lambda = diag(\\lambda_1, \\lambda_2)$, we have $AR = R\\Lambda$, or $A = R \\Lambda R^{-1}$.  Substituting this into \\eqref{Ac:main} yields\n",
    "\\begin{align*}\n",
    "q_t + A q_x & = 0 \\\\\n",
    "q_t + R \\Lambda R^{-1} q_x & = 0 \\\\\n",
    "R^{-1}q_t + \\Lambda R^{-1} q_x & = 0 \\\\\n",
    "w_t + \\Lambda w_x & = 0,\n",
    "\\end{align*}\n",
    "where we have introduced the *characteristic variables* $w=R^{-1}q$.  The last system above is simply a pair of decoupled advection equations for $w_1$ and $w_2$, with velocities $\\lambda_1$ and $\\lambda_2$.  Thus we see that the eigenvalues of $A$ are the velocities at which information propagates in the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution by characteristics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The discussion above suggests a strategy for solving the Cauchy problem:\n",
    "\n",
    "1. Decompose the initial data $(p(x,0), u(x,0))$ into characteristic variables $(w_1^0(x),w_2^0(x,0))$ using the relation $w = R^{-1}q$.\n",
    "2. Evolve the characteristic variables: $w_p(x,t) = w_p^0(x-\\lambda_p t)$.\n",
    "3. Transform back to the physical variables: $q = Rw$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in this process amounts to expressing the vector $q$ in the basis given by $r_1, r_2$.  Solving the system $$Rw=q$$ yields $$w = \\alpha_1 r_1 + \\alpha_2 r_2,$$ where\n",
    "\\begin{align*}\n",
    "\\alpha_1 = \\frac{- p + Z u}{2Z}, \\ \\ \\ \\ \\ \\\n",
    "\\alpha_2 = \\frac{ p + Z u}{2Z}.\n",
    "\\end{align*}\n",
    "\n",
    "We visualize this below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "from exact_solvers import acoustics, interactive_pplanes\n",
    "from utils import riemann_tools\n",
    "colors = ['g','orange']\n",
    "import seaborn as sns\n",
    "sns.set_style('white',{'legend.frameon':'True'});\n",
    "#sns.set_context('paper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decompose_q(p,u,K,rho):\n",
    "    # Should also print the eigenvectors and the values w_1, w_2\n",
    "    Z = np.sqrt(K*rho)\n",
    "    fig, axes = plt.subplots(1,2,figsize=(10,6))\n",
    "    axes[0].arrow(0,0,-Z,1,head_width=0.05, head_length=0.1, color=colors[0])\n",
    "    axes[0].arrow(0,0,Z,1, head_width=0.05, head_length=0.1, color=colors[1])\n",
    "    l1 = axes[0].plot([],[],colors[0])\n",
    "    l2 = axes[0].plot([],[],'-',color=colors[1])\n",
    "    axes[0].set_xlim(-2,2)\n",
    "    axes[0].set_ylim(-2,2)\n",
    "    axes[0].set_aspect('equal')\n",
    "    axes[0].set_title('Eigenvectors')\n",
    "    axes[0].legend(['$r_1$','$r_2$'],loc=3)\n",
    "    axes[0].plot([0,0],[-2,2],'--k',alpha=0.2)\n",
    "    axes[0].plot([-2,2],[0,0],'--k',alpha=0.2)\n",
    "\n",
    "    \n",
    "    axes[1].plot([0,p],[0,u],'k',lw=3)    \n",
    "    alpha1 = (Z*u-p)/(2.*Z)\n",
    "    alpha2 = (Z*u+p)/(2.*Z)\n",
    "    axes[1].plot([0,-Z*alpha1],[0,1*alpha1], color=colors[0], lw=3)\n",
    "    axes[1].plot([-Z*alpha1,-Z*alpha1+Z*alpha2],[1*alpha1,alpha1+alpha2], color=colors[1], lw=3)\n",
    "    axes[1].set_xlim(-1.2,1.2)\n",
    "    axes[1].set_ylim(-1.2,1.2)\n",
    "    axes[1].set_aspect('equal')\n",
    "    axes[1].legend(['$q$',r'$\\alpha_1 r_1$',r'$\\alpha_2 r_2$'],loc='best')\n",
    "    axes[1].plot([0,0],[-2,2],'--k',alpha=0.2)\n",
    "    axes[1].plot([-2,2],[0,0],'--k',alpha=0.2)\n",
    "\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(decompose_q,p=widgets.FloatSlider(min=-1,max=1,value=1.),\n",
    "                     u=widgets.FloatSlider(min=-1,max=1,value=0.),\n",
    "                     rho=widgets.FloatSlider(min=0.1,max=2,value=1.,description=r'$\\rho$'),\n",
    "                     K=widgets.FloatSlider(min=0.1,max=2,value=1.));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the second and third steps, we exactly evolve the characteristic variables and then transform back to the original variables.  This is visualized below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def char_solution(K, rho, t):\n",
    "    fig, axes = plt.subplots(1,2,figsize=(12,6))\n",
    "    c = np.sqrt(K/rho)\n",
    "    x = np.linspace(-2*c-1,2*c+1,41)\n",
    "    tt = np.linspace(0,2,20)\n",
    "    for ix in x:\n",
    "        axes[0].plot(ix-c*tt,tt,'-k',lw=0.5,color=colors[0])\n",
    "        axes[0].plot(ix+c*tt,tt,'-k',lw=0.5,color=colors[1])\n",
    "    axes[0].set_xlim(-1,1)\n",
    "    axes[0].set_ylim(-0.2,2)\n",
    "    axes[0].set_title('Characteristics')\n",
    "    axes[0].set_xlabel('$x$')\n",
    "    axes[0].set_ylabel('$t$')\n",
    "    \n",
    "    xx = np.linspace(-2*c-1,2*c+1,1000)\n",
    "    w120 = lambda x: -0.1*np.exp(-50*x**2)\n",
    "    w220 = lambda x:  0.1*np.exp(-50*x**2)\n",
    "    spacing = 1\n",
    "    l1, = axes[0].plot(xx,w120(xx+c*spacing*t)+spacing*t,color=colors[0],lw=2,label='$w_{12}$')\n",
    "    l2, = axes[0].plot(xx,w220(xx-c*spacing*t)+spacing*t,color=colors[1],lw=2,label='$w_{22}$')\n",
    "    axes[0].legend(handles=[l1,l2], loc=4)\n",
    "    axes[1].plot(xx,w120(xx-c*spacing*t)+w220(xx+c*spacing*t)+spacing*t,'-k',lw=2)\n",
    "    axes[1].set_xlim(-1,1)\n",
    "    axes[1].set_ylim(-0.2,2)\n",
    "    axes[1].set_title('Velocity')\n",
    "    axes[1].set_xlabel('$x$')\n",
    "    axes[1].set_ylabel('$t$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(char_solution, \n",
    "         rho=widgets.FloatSlider(min=0.1,max=2,value=1.,description=r'$\\rho$'),\n",
    "         K=widgets.FloatSlider(min=0.1,max=2,value=1.),\n",
    "         t=widgets.FloatSlider(min=0.,max=2,value=0.));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how in the characteristic variables (plotted on the left), each part of the solution simply advects (translates)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Riemann problem\n",
    "\n",
    "Now that we know how to solve the Cauchy problem, solution of the Riemann problem is merely a special case.  We have the special initial data  \n",
    "\\begin{align*}\n",
    "q(x,0) = \\begin{cases}\n",
    "q_\\ell & \\text{if   } x \\le 0, \\\\\n",
    "q_r & \\text{if   } x > 0.\n",
    "\\end{cases}\n",
    "\\end{align*}  \n",
    "We can proceed as before, by decomposing into characteristic components, advecting, and then transforming back.  But since we know the solution will be constant almost everywhere, it's even simpler to just decompose the jump $\\Delta q = q_r - q_\\ell$ in terms of the characteristic variables, and advect the two resulting jumps $\\Delta w_1$ and $\\Delta w_2$:  \n",
    "\\begin{align*}\n",
    "q_r-q_\\ell = \\Delta q = \\alpha_1 r_1 + \\alpha_2 r_2,\n",
    "\\end{align*}  \n",
    "Since $R\\Delta w = \\Delta q$, we have  \n",
    "\\begin{align*}\n",
    "\\alpha_1 = \\frac{-\\Delta p + Z\\Delta u}{2Z}, \\ \\ \\ \\ \\ \\\n",
    "\\alpha_2 = \\frac{\\Delta p + Z\\Delta u}{2Z}.\n",
    "\\end{align*}  \n",
    "\n",
    "Thus the solution has the structure depicted below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Fig 4.1: Structure of the Riemann solution for acoustics.](./figures/acoustics_xt_plane.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three constant states are related by the jumps:   \n",
    "\\begin{align}\n",
    "q_m = q_\\ell + \\alpha_1 r_1 = q_r - \\alpha_2 r_2.\n",
    "\\label{eq:acussol}\n",
    "\\end{align}  \n",
    "Note that the form of the eigenvectors shows that in any propagating discontinuity, the jump in $p$ is $\\pm Z$ times the jump in $u$.  More generally, the eigenvectors of the coefficient matrix of a linear hyperbolic system reveal the relation between jumps in the different components of $q$ across a wave propagating with speed given by the corresponding eigenvalue.  For acoustics, the impedance is the physical parameter that determines this relation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simple solution\n",
    "Here we provide some very simple initial data, and we call the linear Riemann solver. This will output the three states $q_l$, $q_m$ and $q_r$, and the speeds of the two waves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = 1.               # density\n",
    "bulk = 4.              # bulk modulus\n",
    "c = np.sqrt(bulk/rho)  # sound speed\n",
    "Z = np.sqrt(bulk*rho)     # impedance\n",
    "\n",
    "print(\"Density   rho = %g,  Bulk modulus K = %g\" % (rho,bulk))\n",
    "print(\"Sound speed c = %g,  Impedance    Z = %g\" % (c,Z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ql = np.array([1,2])  # Left state\n",
    "qr = np.array([2,-2])  # Right state\n",
    "states, speeds, reval = acoustics.exact_riemann_solution(ql ,qr, [rho, bulk])\n",
    "print(\"The states ql, qm and qr are: \")\n",
    "print(states)\n",
    "print(\" \")\n",
    "print(\"The left and right wave speeds are:\")\n",
    "print(speeds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also show the structure of the solution in the $p-u$ phase plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppplot=interactive_pplanes.acoustics_phase_plane_plot()\n",
    "ppplot(ql[0],ql[1],qr[0],qr[1],rho,bulk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive solution in the phase plane\n",
    "As we already showed in the previous section, one way to understand the solution to a linear system like acoustics equations is by looking at the phase plane. The middle state $q_m$ generated between the two waves must be connected to $q_l$ by a multiple of the first eigenvector and to $q_r$ by a multiple of the second eigenvector,as it is evident from equation (\\ref{eq:acussol}). Therefore, the solution of the Riemann\n",
    "problem is nothing more but the intersection of the line generated by the first eigenvector going through $q_l$ with the line generated by the second eigenvector going through $q_r$. This is easier to understand visually, like in the interactive plot we show next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial states [pressure, velocity]\n",
    "ql = [10.0, -5.0]\n",
    "qr = [40.0, 5.0]\n",
    "\n",
    "# Acoustic eqs. parameters\n",
    "rho = 2.0\n",
    "bulk = 1.0\n",
    "\n",
    "interactive_pplanes.acoustics_interactive_phase_plane(ql,qr,rho,bulk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the eigenvectors are given in terms of the impedance $Z$, which depends on the density $\\rho$\n",
    "and the bulk modulus $K$. Therefore, when $\\rho$ and $K$ are modified the eigenvectors change and consequently the slope of the lines changes as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples\n",
    "We will begin by defining a function that calls the exact solver in [exact_solvers/acoustics.py](exact_solvers/acoustics.py) and plots the solution for different interesting examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_riemann_solution(ql, qr, aux):\n",
    "    ex_states, ex_speeds, reval = acoustics.exact_riemann_solution(ql ,qr, aux)\n",
    "\n",
    "    plot_function = riemann_tools.make_plot_function(ex_states, ex_speeds, reval, layout='vertical',\n",
    "                                                    variable_names=['pressure', 'velocity'],\n",
    "                                                    aux=(np.array(aux),np.array(aux)), \n",
    "                                                    plot_chars=[acoustics.lambda1,\n",
    "                                                                acoustics.lambda2])\n",
    "\n",
    "    return interact(plot_function, t=widgets.FloatSlider(value=0.0,min=0,max=1.0),\n",
    "                    which_char=widgets.Dropdown(options=[None,1,2],description='Show characteristics'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shock tube\n",
    "\n",
    "If the velocity is 0 in both initial states (the shock tube problem) then the resulting Riemann solution consists of pressure jumps of equal magnitude propagating in each direction, with equal and opposite jumps in velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ql = np.array([5,0])\n",
    "qr = np.array([1,0])\n",
    "rho = 1.0\n",
    "bulk = 4.0\n",
    "plot_riemann_solution(ql, qr, [rho, bulk]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution in the phase plane\n",
    "The solution in the phase plane is given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppplot(ql[0],ql[1],qr[0],qr[1],rho,bulk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reflection at a wall\n",
    "\n",
    "As another example, suppose the pressure is initially the same in the left and right states, while the velocities are non-zero with $u_r = -u_\\ell > 0$.  Particles are converging from both sides and if the initial states have this symmetry, then the result is a middle state $q_m$ in which the velocity is 0 (and the pressure is higher than on either side)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ql = np.array([3,2])  \n",
    "qr = np.array([3,-2])  \n",
    "rho = 1.0\n",
    "bulk = 20.0\n",
    "plot_riemann_solution(ql, qr, [rho, bulk]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution in the phase plane and the particle trajectories are the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppplot=interactive_pplanes.acoustics_phase_plane_plot()\n",
    "rho = 1\n",
    "bulk = 1\n",
    "ppplot(ql[0],ql[1],qr[0],qr[1],rho,bulk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you discard half the solution (for $x>0$ or for $x<0$) then what you see can be viewed as the solution to a problem with fluid streaming at constant velocity toward a solid wall.  The result is an acoustic wave that moves away from the wall, and the fluid behind the shock has been decelerated to velocity 0, i.e. it is stationary at the wall.\n",
    "\n",
    "This type of Riemann solution is critical when imposing solid wall boundary conditions in a numerical method. If ghost cells are introduced outside the domain and the state in the ghost cell set by reflecting the interior solution with the symmetry seen here (equal pressure, negated velocity), then the solution to the Riemann problem at the cell interfaces yields a solution that satisfies the desired boundary conditions. \n",
    "\n",
    "Note it is possible to extend the Riemann problem solution for acoustic equations to cases where there are different materials on the left and right side. This is useful to solve the acoustic wave propagation across interfaces or heterogeneous media, and it will be explored further in the section on acoustic equations for heterogeneous media.\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}