Pressureless_flow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pressureless flow\n",
    "\n",
    "In this chapter we consider flow in the absence of pressure.\n",
    "\n",
    "## Shallow water flow in low gravity\n",
    "Recall the shallow water equations:\n",
    "\n",
    "\\begin{align}\n",
    "    h_t + (hu)_x & = 0 \\\\\n",
    "    (hu)_t + \\left(hu^2 + \\frac{1}{2}gh^2\\right)_x & = 0.\n",
    "\\end{align}\n",
    "\n",
    "These are very similar to the *isothermal flow* equations of gas dynamics:\n",
    "\n",
    "\\begin{align}\n",
    "    \\rho_t + (\\rho u)_x & = 0 \\\\\n",
    "    (\\rho u)_t + \\left(\\rho u^2 + a^2 \\rho \\right)_x & = 0\n",
    "\\end{align}\n",
    "\n",
    "Indeed, if we identify the shallow water depth $h$ with the isothermal gas density $\\rho$, then these systems differ only in the second term of the momentum flux.  In both systems this term represents pressure; in the first system the hydrostatic pressure $\\frac{1}{2}gh^2$ creates a net force from regions of greater depth to lower depth, while in the second system the pressure $a^2 \\rho$ increases linearly with density.\n",
    "\n",
    "In this chapter we investigate what happens when the pressure tends to zero; this corresponds to the limit $g\\to 0$ for shallow water and the limit $a \\to 0$ for isothermal flow.\n",
    "\n",
    "First, let's see what happens to the integral curves and Hugoniot loci as $g \\to 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from exact_solvers import shallow_water\n",
    "from collections import namedtuple\n",
    "from utils import riemann_tools\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets, Checkbox, IntSlider, FloatSlider\n",
    "State = namedtuple('State', shallow_water.conserved_variables)\n",
    "Primitive_State = namedtuple('PrimState', shallow_water.primitive_variables)\n",
    "plt.style.use('seaborn-talk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def connect_states(h_l=1.,u_l=-1.,h_r=1.,u_r=1.,logg=0.,plot_unphysical=True):\n",
    "    g = 10.**logg\n",
    "    q_l = np.array([h_l,h_l*u_l])\n",
    "    q_r = np.array([h_r,h_r*u_r])\n",
    "    fig, ax = plt.subplots(1,2,figsize=(14,6))\n",
    "    shallow_water.phase_plane_curves(q_l[0], q_l[1], 'qleft', g,\n",
    "                                     wave_family=1, ax=ax[0],\n",
    "                                     plot_unphysical=plot_unphysical)\n",
    "    shallow_water.phase_plane_curves(q_r[0], q_r[1], 'qright', g,\n",
    "                                     wave_family=2, ax=ax[0],\n",
    "                                     plot_unphysical=plot_unphysical)\n",
    "    shallow_water.phase_plane_curves(q_l[0], q_l[1], 'qleft', g,\n",
    "                                     wave_family=1, y_axis='hu',ax=ax[1],\n",
    "                                     plot_unphysical=plot_unphysical)\n",
    "    shallow_water.phase_plane_curves(q_r[0], q_r[1], 'qright', g,\n",
    "                                     wave_family=2, y_axis='hu',ax=ax[1],\n",
    "                                     plot_unphysical=plot_unphysical)\n",
    "    ax[0].set_title('h-u plane'); ax[1].set_title('h-hu plane'); \n",
    "    ax[0].set_xlim(0,3); ax[1].set_xlim(0,3);\n",
    "    ax[0].set_ylim(-10,10); ax[1].set_ylim(-10,10); \n",
    "    plt.tight_layout(); plt.show()\n",
    "    \n",
    "interact(connect_states,\n",
    "         h_l=widgets.FloatSlider(min=0.001,max=2,value=1),\n",
    "         u_l=widgets.FloatSlider(min=-5,max=5,value=-1),\n",
    "         h_r=widgets.FloatSlider(min=0.001,max=2,value=1),\n",
    "         u_r=widgets.FloatSlider(min=-5,max=5,value=1),\n",
    "         logg=widgets.FloatSlider(value=0,min=-5,max=2,\n",
    "                                  description='$\\log_{10}(g)$'));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def c1(q, xi, g=1.):\n",
    "    \"Characteristic speed for shallow water 1-waves.\"\n",
    "    h = q[0]\n",
    "    if h > 0:\n",
    "        u = q[1]/q[0]\n",
    "        return u - np.sqrt(g*h)\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "def c2(q, xi, g=1.):\n",
    "    \"Characteristic speed for shallow water 2-waves.\"\n",
    "    h = q[0]\n",
    "    if h > 0:\n",
    "        u = q[1]/q[0]\n",
    "        return u + np.sqrt(g*h)\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_plot_function(q_l,q_r,force_waves=None,extra_lines=None):\n",
    "    def plot_function(t,logg,plot_1_chars=False,plot_2_chars=False):\n",
    "        varnames = shallow_water.primitive_variables\n",
    "        prim = shallow_water.cons_to_prim\n",
    "        g = 10**logg\n",
    "        states, speeds, reval, wave_types = \\\n",
    "            shallow_water.exact_riemann_solution(q_l,q_r,g,\n",
    "                                                 force_waves=force_waves)\n",
    "        ax = riemann_tools.plot_riemann(states,speeds,reval,wave_types,t=t,\n",
    "                                        t_pointer=0,extra_axes=True,\n",
    "                                        variable_names=varnames,fill=[0],\n",
    "                                        derived_variables=prim);\n",
    "        if plot_1_chars:\n",
    "            riemann_tools.plot_characteristics(reval,c1,(g,g),axes=ax[0],\n",
    "                                               extra_lines=extra_lines)\n",
    "        if plot_2_chars:\n",
    "            riemann_tools.plot_characteristics(reval,c2,(g,g),axes=ax[0],\n",
    "                                               extra_lines=extra_lines)\n",
    "        shallow_water.phase_plane_plot(q_l,q_r,g,ax=ax[3],\n",
    "                                       force_waves=force_waves,y_axis='u')\n",
    "        plt.show()\n",
    "    return plot_function\n",
    "\n",
    "def plot_riemann_SW(q_l,q_r,force_waves=None,extra_lines=None):\n",
    "    plot_function = make_plot_function(q_l,q_r,force_waves,extra_lines)\n",
    "    interact(plot_function, t=FloatSlider(value=0.1,min=0,max=.9),\n",
    "             plot_1_chars=Checkbox(description='1-characteristics',value=False),\n",
    "             plot_2_chars=Checkbox(description='2-characteristics'),\n",
    "             logg=IntSlider(value=0,min=-10,max=2,description='$\\log_{10}(g)$'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the following behaviors as $g \\to 0$:\n",
    "\n",
    "- The integral curves and hugoniot loci become parallel.\n",
    "- All four curves become horizontal lines in the $h-u$ plane."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens to the hyperbolic structure of this system as $g \\to 0$?  Recall that the flux jacobian is\n",
    "\n",
    "\\begin{align}\n",
    "f'(q) & = \\begin{pmatrix} 0 & 1 \\\\ -(q_2/q_1)^2 + g q_1 & 2 q_2/q_1 \\end{pmatrix} \n",
    "        = \\begin{pmatrix} 0 & 1 \\\\ -u^2 + g h & 2 u \\end{pmatrix},\n",
    "\\end{align}\n",
    "\n",
    "with eigenvalues\n",
    "\n",
    "\\begin{align}\n",
    "    \\lambda_1 & = u - \\sqrt{gh} & \\lambda_2 & = u + \\sqrt{gh},\n",
    "\\end{align}\n",
    "\n",
    "and corresponding eigenvectors\n",
    "\n",
    "\\begin{align}\n",
    "    r_1 & = \\begin{bmatrix} 1 \\\\ u-\\sqrt{gh} \\end{bmatrix} &\n",
    "    r_2 & = \\begin{bmatrix} 1 \\\\ u+\\sqrt{gh} \\end{bmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "As $g \\to 0$, the eigenvalues both approach a single value of $\\lambda_0 = u$ and the eigenvectors both approach\n",
    "\n",
    "\\begin{align}\n",
    "    r_0 & = \\begin{bmatrix} 1 \\\\ u \\end{bmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "Thus for $g=0$ there is a single eigenvalue $\\lambda_0=u$ with algebraic multiplicity 2 but geometric multiplicity 1.  The flux Jacobian is defective and the system is no longer hyperbolic.  The integral curves are everywhere parallel to the eigenvector $r_0$, but this eigenvector points in a direction of constant $u$ and depends only on $u$, so the integral curves in the $h-hu$ plane are straight lines of constant $u$, all passing through the origin.  \n",
    "\n",
    "What about the Hugoniot loci?  Recall from the shallow water chapter that a shock with a jump in the depth of $\\alpha = h-h_*$ must have a jump in the momentum of\n",
    "\n",
    "\\begin{align}\n",
    "    h u & = h_* u_* + \\alpha \\left[u_* \\pm \\sqrt{gh_* \\left(1+\\frac{\\alpha}{h_*}\\right)\\left(1+\\frac{\\alpha}{2h_*}\\right)}\\right]\n",
    "\\end{align}\n",
    "\n",
    "Setting $g=0$ we obtain $hu = h_* u_* + (h-h_*)u_*$, so that $u=u_*$.  Thus the Hugoniot loci are the same as the integral curves and are lines of constant $u$.\n",
    "\n",
    "In the Riemann problem, if the left and right states have the same initial velocity, then we can still connect them as $g \\to 0$; in the limit, everything just moves along at that constant velocity.  Here's an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l  = State(Depth = 3.,\n",
    "             Momentum = 1.)\n",
    "q_r = State(Depth = 1.,\n",
    "            Momentum = 1./3)\n",
    "\n",
    "plot_riemann_SW(q_l,q_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose we have a Riemann problem whose initial states have different velocities.  The only way to connect them is through a pair of rarefactions with a dry state in between.  Recall from our earlier treatment of the shallow water equations that a middle dry state occurred if $u_l + 2 \\sqrt{gh_l} < u_r-2\\sqrt{gh_r}$; for $g=0$ this reduces simply to the condition $u_l < u_r$.\n",
    "Here's an example with $u_l < u_r$; use the second slider to see what happens as $g \\to 0$.\n",
    "\n",
    "Also, turn on the plots of the two characteristic families and notice how they become parallel as $g \\to 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l  = State(Depth = 3.,\n",
    "             Momentum = -1.)\n",
    "q_r = State(Depth = 1.,\n",
    "            Momentum = 2.)\n",
    "\n",
    "plot_riemann_SW(q_l,q_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens if $u_l \\ge u_r$, so that the double-rarefaction solution is unphysical?  The answer can be found by considering what happens when $g$ is small but nonzero.  Here's an example to play with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l  = State(Depth = 3.,\n",
    "             Momentum = 1.)\n",
    "q_r = State(Depth = 1.,\n",
    "            Momentum = 0.1)\n",
    "\n",
    "plot_riemann_SW(q_l,q_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep a close eye on the scale of the h-axis in the phase plane plot, and notice again that as $g \\to 0$ the hugoniot loci become nearly parallel (horizontal in the $h-u$ plane).  This means that their intersection occurs at a very large value of $h$, so the middle state depth goes to $\\infty$ as $g\\to 0$.  Meanwhile, the speed of both shocks approaches a single value, so the region occupied by the middle state gets ever narrower.  In the limit, the solution for the depth involves a delta function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make this even clearer, here's an example with uniform initial depth and both flows directed inward, very similar to the first example we considered in [the initial chapter on shallow water](Shallow_water.ipynb#All-shock-Riemann-solutions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l  = State(Depth = 1.,\n",
    "             Momentum = 2.)\n",
    "q_r = State(Depth = 1.,\n",
    "            Momentum = -1.)\n",
    "\n",
    "plot_riemann_SW(q_l,q_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's going on here, physically, when $g$ is very small? At the interface, water is colliding.  With the gravitational pressure term present, this leads to the formation of outgoing shock waves that redistribute water away from the interface.  But without that term, there is no force to equilibrate the depth and each parcel of water just flows with its initial velocity, until eventually it collides with water of a different velocity.  Water from both sides accumulates at the interface, leading to what is known as a *delta shock*.\n",
    "\n",
    "If you drag the time slider for different fixed values of $g$ in the last example, you'll see that since $g$ is still nonzero there are two shock waves and the near-delta-function region expands with time, but the rate of this expansion is smaller for smaller values of $g$.  For very small values, the region between the shocks cannot be resolved on the scale of the plots above.  Indeed, the width of the line plotted there significantly exaggerates the width of this region when $g$ is very small.\n",
    "\n",
    "There are two import details we haven't yet explained: first, the delta shock is moving; second, it is growing in time.  Movement of the delta shock is necessary in order to conserve momentum; since the water from the left is arriving with a different speed than that from the right, if the delta shock were stationary then momentum would not be conserved.  Growth of the delta shock is necessary in order to conserve mass.\n",
    "\n",
    "Let $h_\\delta$ denote the mass of the delta shock.  How must $h_\\delta$ grow in order to conserve mass?  Let $\\hat{u}$ denote the velocity of the delta shock, with $u_l \\ge \\hat{u} \\ge u_r$.  Then over a unit time interval, a quantity $h_l(u_l-\\hat{u})$ of water arrives from the left, while $h_r(\\hat{u}-u_r)$ arrives from the right.  Thus we must have\n",
    "\n",
    "$$\n",
    "    h_\\delta(t) = t \\left( h_l u_l - h_r u_r + \\hat{u}(h_r-h_l) \\right).\n",
    "$$\n",
    "\n",
    "At what speed must the water in the delta shock move in order to conserve momentum?  The rate of momentum flowing into the delta shock from the left is $h_l u_l (u_l-\\hat{u})$ and the rate of momentum flowing into the delta shock from the right is $h_r u_r (\\hat{u}-u_r)$, so we must have\n",
    "\n",
    "$$\n",
    "    h_\\delta(t) \\hat{u} = t \\left( h_l u_l^2 - h_r u_r^2 + \\hat{u}(h_r u_r - h_l u_l) \\right).\n",
    "$$\n",
    "\n",
    "Combining these two conservation conditions, we find the speed of the delta shock:\n",
    "\n",
    "$$\n",
    "    \\hat{u} = \\frac{\\sqrt{h_l}u_l + \\sqrt{h_r}{u_r}}{\\sqrt{h_l}+\\sqrt{h_r}}\n",
    "$$\n",
    "\n",
    "and its rate of growth:\n",
    "\n",
    "$$\n",
    "    h_\\delta(t) = t \\sqrt{h_l h_r}(u_l - u_r).\n",
    "$$\n",
    "\n",
    "Note that for any non-zero value of $g$, the length of the interval occupied by the middle state grows in time while its amplitude is constant.  In the limit $g\\to 0$ the interval becomes infinitesimal so the amplitude must change to conserve mass.  Notice also that this change in amplitude implies that the Riemann solution is no longer a similarity solution when $g=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## High-speed flows\n",
    "Above we have considered the limit as gravity tends to zero in the shallow water equations (or equivalently as the sound speed goes to zero in the isothermal flow equations).  As with any mathematical analysis involving a \"small\" parameter, one should ask \"small relative to what\"?\n",
    "\n",
    "By non-dimensionalizing the shallow water equations, one obtains the system\n",
    "\n",
    "\\begin{align}\n",
    "    h_t + (hu)_x & = 0 \\\\\n",
    "    (hu)_t + \\left(hu^2 + \\frac{1}{2F^2}h^2\\right)_x & = 0.\n",
    "\\end{align}\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "    F = \\frac{U}{\\sqrt{gH}}\n",
    "$$\n",
    "\n",
    "is the ratio of a typical velocity $U$ to the speed of gravity waves with typical depth $H$.  From this form we can see that in fact we have studied the *large Froude number limit* of the shallow water equations, and we should expect to see something similar to a delta shock form whenever the Froude number is large.  We can approach this limit in ways other than taking $g \\to 0$.  For instance, if the velocity is very large and $\\sqrt{gh}=O(1)$, then we obtain (approximately) delta-shock solutions like those seen above.  Here is an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l  = State(Depth = 1.,\n",
    "             Momentum = 100.)\n",
    "q_r = State(Depth = 1.,\n",
    "            Momentum = 0.)\n",
    "\n",
    "plot_riemann_SW(q_l,q_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, a non-dimensionalization of the isothermal equations shows that the pressureless flow equations correspond to the large Mach number limit (i.e. when the typical flow velocity is much greater than the sound speed)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "Delta shocks arising in pressureless flows have been studied by many authors; see Section 16.3 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque, 2002)</a></cite> and references therein, as well as <cite data-cite=\"leveque2004dynamics\"><a href=\"riemann.html#leveque2004dynamics\">(LeVeque, 2004)</a></cite>.  For an examination of the large Froude number limit in shallow water flow, see <cite data-cite=\"edwards2008non\"><a href=\"riemann.html#fvmhp\">(Edwards et. al., 2008)</a></cite>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "(1) In the example below, the initial velocity is zero for both states.  Observe that multiplying the value of $g$ by a factor $\\alpha$ is equivalent to rescaling all the velocities by $\\sqrt{\\alpha}$.  Prove that this scaling must hold for any Riemann solution of this system when $u_l = u_r = 0$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l  = State(Depth = 3.,\n",
    "             Momentum = 0.)\n",
    "q_r = State(Depth = 1.,\n",
    "            Momentum = 0.)\n",
    "\n",
    "plot_riemann_SW(q_l,q_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(2) Here we have considered the large Froude number limit.  What happens in the low Froude number limit?  What form do the integral curves and Hugoniot loci take?  What happens to the hyperbolic structure of the system?  To the solution of the Riemann problem?  You may wish to modify the examples above to allow for much larger values of $g$ in order to check your answers."
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}