Traffic_variable_speed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LWR Traffic flow with varying speed limit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we again consider the LWR traffic model, but now we introduce a speed limit $v$:\n",
    "\n",
    "\\begin{align} \\label{LWR_vc:conslaw}\n",
    "\\rho_t + (v(x) \\rho(1-\\rho))_x = 0.\n",
    "\\end{align}  \n",
    "and we consider the Riemann problem with different speed limits to the left and right:\n",
    "\n",
    "$$\n",
    "v(x) = \\begin{cases} v_l & x<0 \\\\ v_r & x>0. \\end{cases}\n",
    "$$\n",
    "\n",
    "In [the earlier LWR chapter](Traffic_flow.ipynb), the speed limit was set to unity everywhere.  The variable-speed-limit case has been considered, for instance, in <cite data-cite=\"mochon1987analysis\"><a href=\"riemann.html#mochon1987analysis\">(Mochon, 1987)</a></cite> and in Chapter 16 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque, 2002)</a></cite>.  From a physical point of view, we might also imagine that $u_\\text{max}$ varies because of differing road conditions -- for instance, if part of the road is wet or foggy.\n",
    "\n",
    "\n",
    "Correspondingly, we have a flux function that is different on either side of $x=0$.  Let $f(\\rho,v) = v \\rho (1-\\rho)$ and define\n",
    "\n",
    "\\begin{align}\n",
    "    f_l & = f(\\rho_l,v_l) \\\\\n",
    "    f_r & = f(\\rho_r,v_r).\n",
    "\\end{align}\n",
    "\n",
    "The presence of a discontinuous flux function in the Riemann problem can introduce new difficulties, some of which will be illustrated here.  For a more detailed discussion of other systems with discontinuous flux, see <cite data-cite=\"Burger2008\"><a href=\"riemann.html#fvmhp\">(Burger et. al., 2008)</a></cite> and other papers in that special issue.\n",
    "\n",
    "Continuing with the problem at hand, if, say, $v_l=1$ and $v_r=2$ then the two flux functions look like this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['font.size'] = 8\n",
    "figsize =(8,4)\n",
    "mpl.rcParams['figure.figsize'] = figsize\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets, FloatSlider\n",
    "from utils import riemann_tools\n",
    "from exact_solvers import traffic_variable_speed\n",
    "from clawpack import pyclaw\n",
    "from clawpack import riemann"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho = np.linspace(0,1)\n",
    "f_l = rho*(1-rho); f_r = 2*rho*(1-rho)\n",
    "plt.plot(rho,f_l,'-c',rho,f_r,'-g'); plt.xlim(0,1); plt.ylim(0,0.6);\n",
    "plt.xlabel(r'$\\rho$'); plt.ylabel(r'$f(\\rho)$')\n",
    "plt.legend(['$f_l$','$f_r$']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Admissible regions\n",
    "We can view the solution of the Riemann problem as the task of finding a physically admissible path connecting the left and right states in the $\\rho-f(\\rho)$ plane.  The structure of the Riemann solution can be fairly complicated, but we can simplify things by recognizing a few \"rules\" that must govern the connecting path.  Two quick observations will rule out half of the available intermediate states:\n",
    "\n",
    "1. Since the characteristic speed is positive for $\\rho<1/2$, **we cannot reach any intermediate state with $\\rho<1/2$ on the $f_l$ curve** (since states on the $f_l$ curve must appear to the left of $x=0$).\n",
    "2. Similarly, since the characteristic speed is negative for $\\rho>1/2$, **we cannot reach any intermediate state with $\\rho>1/2$ on the $f_r$ curve**.\n",
    "\n",
    "To remind us of these facts, we'll plot the prohibited parts of the curves with dashed lines.  It is important to note that the initial states $\\rho_l$ and $\\rho_r$ **can** inhabit these regions.  However, if an initial state is on the inadmissible portion of the curve, then we cannot connect a rarefaction to it, since nearby states are inadmissible.  We must jump to it through a shock or (as discussed below) a stationary jump."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_l = rho*(1-rho); f_r = 2*rho*(1-rho)\n",
    "plt.plot(rho[26:],f_l[26:],'-c',rho[:26],f_r[:26],'-g');\n",
    "plt.xlim(0,1); plt.ylim(0,0.6);\n",
    "plt.plot(rho[:26],f_l[:26],'--k',rho[26:],f_r[26:],'--k',alpha=0.5);\n",
    "plt.xlim(0,1); plt.ylim(0,0.6);\n",
    "plt.xlabel(r'$\\rho$'); plt.ylabel(r'$f(\\rho)$')\n",
    "plt.legend(['$f_l$','$f_r$']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shocks and rarefactions\n",
    "We can also determine in advance what kind of waves can be connected to $\\rho_l$ or $\\rho_r$ within the admissible regions:\n",
    "\n",
    "1.  Any intermediate state on the $f_l$ curve with density greater than $\\rho_l$ must be connected to $\\rho_l$ by a shock.  Conversely, any intermediate state with density less than $\\rho_l$ must be connected to $\\rho_l$ by a rarefaction.\n",
    "2.  Similarly, any intermediate state on the $f_r$ curve with density greater than $\\rho_r$ must be connected to $\\rho_r$ by a rarefaction, and any intermediate state with density less than $\\rho_r$ must be connected to $\\rho_r$ by a shock.\n",
    "\n",
    "These conditions follow directly from the Lax entropy condition.  To remind us of these conditions, we'll plot in red the portions of the admissible curves that can be connected to the left or right state by a shock, and in blue those that can be connected by a rarefaction.  Here's an example of what that looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_l = 0.6; rho_r = 0.6\n",
    "v_l = 1.0; v_r = 0.8\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                        v_l,v_r,connect=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The stationary wave\n",
    "\n",
    "So far we've only discussed how to connect the left and right states to other states on the same flux curve, but to solve the Riemann problem we must connect points on the two different flux curves to each other.\n",
    "\n",
    "As discussed at length in [the earlier traffic flow chapter](Traffic_flow.ipynb), the flux of vehicles must be continuous everywhere -- and in particular, at $x=0$:\n",
    "\\begin{align} \\label{LWR_vc:fluxcont}\n",
    "f(\\rho(0^-,t),v_l)) = f(\\rho(0^+,t),v_r)) \\equiv f^*.\n",
    "\\end{align}\n",
    "\n",
    "If the speed limit changes abruptly at $x=0$, then (in order to have a continuous flux) the density must be discontinuous there.  Since there are two possible values of $\\rho$ corresponding to any given flux, we must carefully determine the intermediate states.  We can think of the discontinuity at $x=0$ as an additional wave with speed zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It may seem strange that the Riemann solution includes two waves when we are dealing with a scalar PDE.  One way to understand this is by adding a trivial evolution equation for the speed limit (which is now artificially considered a function of $t$):\n",
    "\n",
    "\\begin{align} \\label{LWR_vc:2by2}\n",
    "    \\rho_t + (v(x,t) \\rho(1-\\rho))_x & = 0 \\\\\n",
    "    v_t & = 0.\n",
    "\\end{align}\n",
    "\n",
    "The flux jacobian for this system is a $2\\times 2$ matrix with one eigenvalue equal to the characteristic speed $f'(\\rho)$ and the other equal to zero.  Thus we expect its solution to include a wave with speed zero.  We will not investigate the hyperbolic structure of this $2\\times 2$ system further here; instead the flux continuity condition \\eqref{LWR_vc:fluxcont} together with an admissibility condition (discussed below) will be enough to enable us to construct solutions to the Riemann problem based on the scalar equation form \\eqref{LWR_vc:conslaw}."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the flux is continuous at $x=0$, the stationary jump there corresponds to a horizontal line in the $\\rho-f(\\rho)$ plane plotted above.  This jump connects a point on one flux curve to a point on the other.  Two things are immediately clear:\n",
    "\n",
    "1. The value of $f^*$ must not lie fully above either of the flux curves; since the maximum flux is $v/4$, this means that $f^* \\le \\min(v_l,v_r)/4$.  If this condition were violated, it would be impossible to connect the two flux curves with a horizontal line at flux value $f^*$.\n",
    "2. For a given $f^*$ satisfying condition 1, there are always **two** values of $\\rho$ on each curve with flux $f^*$.  This means that if (say) we have determined the value $\\rho(x=0^-,t)$, there are two candidate values of $\\rho(x=0^+,t)$ -- one with positive characteristic velocity, and one with negative characteristic velocity.\n",
    "\n",
    "How can we determine which state to jump to?  We will apply an admissibility condition similar to the Lax entropy condition that allowed us to distinguish whether a wave should be a shock or rarefaction.  Recall that:\n",
    "\n",
    "1. For **shock waves**, characteristics near that of the shock must be impinging from both sides: $c_l>s>c_r$.\n",
    "2. For **rarefaction waves**, at any characteristic within the rarefaction, nearby characteristics must be **diverging** from both sides: $c_l<s<c_r$.\n",
    "3. For **contact waves**, nearby characteristics must be parallel to the contact wave on both sides: $c_l=s=c_r$.\n",
    "\n",
    "The stationary wave corresponds to the zero eigenvalue of the flux Jacobian of \\eqref{LWR_vc:2by2}; evidently\n",
    "all characteristics in this family have speed zero and are thus parallel to the stationary jump.  We cannot use those characteristics to distinguish the admissibility of the jump.  Instead, we impose the following admissibility criteria, related to the original characteristic speed of \\eqref{LWR_vc:conslaw}:\n",
    "\n",
    "1. The characteristic velocity does not change sign across the wave at $x=0$ (in a weak sense). That is, $f_l'(\\rho(x=0^-,t)) f_r'(\\rho(x=0^+,t)) \\ge 0$.\n",
    "2. If the characteristic velocity vanishes on one side of the stationary wave, the value on the other side must lie on the admissible portion of the flux curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the sonic point $\\rho=1/2$ plays such a crucial role, we will say a section of road is *congested* if $\\rho$ is greater than $1/2$ there.\n",
    "\n",
    "Finally, let us formalize one thing that we observed in the previous notebook:\n",
    "\n",
    ">In the LWR system, shocks always carry an increase in density (from left to right), while rarefactions always carry a decrease in density.\n",
    "\n",
    "Why is this?  The Lax entropy condition for a shock moving at speed $s$ tells us\n",
    "\n",
    "$$f'(\\rho_l) > s > f'(\\rho_r).$$\n",
    "\n",
    "In partcular, $f'(\\rho_l) > f'(\\rho_r)$, which for LWR means $1-2\\rho_l > 1-2\\rho_r$, or $\\rho_r > \\rho_l$.  This guarantees that characteristics are converging in the neighborhood of the shock.  At a rarefaction, characteristics must be diverging, so we must have the opposite inequality: $\\rho_r < \\rho_l$.\n",
    "\n",
    "Now let's look at some specific cases and see how to solve the Riemann problem by connecting the left and right states.  We'll start with the states already plotted above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case I: speed limit decrease with $f_l > f_r$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_l = 0.6; rho_r = 0.6\n",
    "v_l = 1.0; v_r = 0.8\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                        v_l,v_r,connect=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the rules determined above, we know that we can connect $\\rho_r$ to an intermediate state only by the stationary jump or by a shock.  If it is connected by a shock, the intermediate state must lie on the red part of the $f_r$ curve.  But then the stationary jump would need to either cross the midline (not allowed) or connect to the inadmissible left portion of the $f_l$ curve (also not allowed).  Hence $\\rho_r$ must be connected by the stationary wave to the admissible (red) portion of the $f_l$ curve.  It is then clear that $\\rho_l$ must be connected to the intermediate state by a shock.\n",
    "\n",
    "Here is what the solution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "def c(rho, xi, v):\n",
    "    return v*(1.-2*rho)\n",
    "    \n",
    "def make_plot_function(rho_l,rho_r,v_l,v_r,connect=True):\n",
    "    states, speeds, reval, wave_types = \\\n",
    "        traffic_variable_speed.exact_riemann_solution(rho_l,rho_r,\n",
    "                                                      v_l,v_r)\n",
    "    def plot_function(t):\n",
    "        ax = riemann_tools.plot_riemann(states,speeds,\n",
    "                          reval,wave_types,t=t,t_pointer=0,\n",
    "                          extra_axes=True,variable_names=['Density']);\n",
    "        riemann_tools.plot_characteristics(reval,c,(v_l,v_r),ax[0],\n",
    "                          extra_lines=[[[0,0],[0,1]]])\n",
    "        traffic_variable_speed.phase_plane_plot(rho_l,rho_r,v_l,v_r,\n",
    "                          axes=ax[2],connect=connect)\n",
    "        plt.show()\n",
    "    return plot_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "def plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r,connect=True):\n",
    "    plot_function = make_plot_function(rho_l,rho_r,v_l,v_r,connect)\n",
    "    interact(plot_function, \n",
    "             t=widgets.FloatSlider(value=0.5,min=0,max=.9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\rho-f$ plane is shown in the rightmost plot above.  The Riemann solution is indicated with a thick line; as usual we indicate shocks in red and rarefactions in blue; the stationary jump at $x=0$ is plotted in black.  Notice that the jump at the interface corresponds to a horizontal line, since the flux is the same on both sides of the interface.  Notice also that the jump does not cross the dashed line at $\\rho=1/2$.  The solution shown satisfies all of the admissibility criteria laid out above; indeed, it is the only solution that does so."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 2: speed limit decrease with $f_l < f_r$\n",
    "In the problem above, it was impossible to connect $\\rho_l$ directly to the $f_r$ curve since $f_l$ was greater than $v_r/4$.  Here is a problem with the same speed limits but a smaller value of $f_l$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_l = 0.9; rho_r = 0.6\n",
    "v_l = 1.0; v_r = 0.8\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                        v_l,v_r,connect=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we connect these states?  By the same reasoning used in the previous problem, $\\rho_r$ must be connected to the $f_l$ curve by the stationary wave.  Then it must be that $\\rho_l$ is connected to the intermediate state by a rarefaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, observe how the waves involved satisfy all of the required conditions.  No other set of waves does so."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 3: Transonic rarefaction\n",
    "So far we have seen solutions with a single shock or rarefaction; they are similar to solutions of the Riemann problem for the LWR model with fixed speed limit (see [this earlier example](Traffic_flow.ipynb#Example:-Traffic-jam) and [this earlier example](Traffic_flow.ipynb#Example:-green-light), except that they also include a stationary jump.  More interesting things can happen if we have a rarefaction in which the characteristic speed passes through zero (i.e., a *transonic rarefaction*).  In order to connect a rarefaction to both states in the variable-speed-limit setting, the densities $\\rho_l$ and $\\rho_r$ bounding the rarefaction wave must both lie on the admissible portions of the curve.  Thus we must have $\\rho_r \\le 1/2 \\le \\rho_l$; this is the same condition that was required for a transonic rarefaction in the basic LWR model.\n",
    "\n",
    "However, in the variable-speed-limit case there are additional conditions that must be satisfied.  Namely, the left and right fluxes must not be so great that we cannot balance them with the appropriate flux on the other side; i.e. we require\n",
    "\n",
    "\\begin{align}\n",
    "f_l & \\le v_r/4 \\label{vcLWR:sufficientA} \\\\\n",
    "f_r & \\le v_l/4.\\label{vcLWR:sufficientB}\n",
    "\\end{align}  \n",
    "This guarantees that neither the left nor right state lies above the entire (right or left) flux curve in the phase plane."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transonic rarefaction with $v_l > v_r$\n",
    "Here's an example of a Riemann problem that leads to a transonic rarefaction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_l = 0.9; rho_r = 0.2\n",
    "v_l = 1.0; v_r = 0.8\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                        v_l,v_r,connect=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exactly how will we connect these states?  In order to obtain an admissible stationary jump, there is only one possibility: we must connect $\\rho_r$ to the sonic point via a rarefaction, and then connect to an admissible point on the right side of the $f_l$ curve via the stationary jump.  Finally, the resulting state is connected to $\\rho_l$ by another rarefaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we have a stationary jump in $\\rho$ in order to maintain flux continuity at the interface, it is impossible that the characteristic velocity approach zero on *both* sides of the stationary jump.  Thus the transonic rarefaction is broken into two rarefactions, with an intermediate constant state lying either just to the left or just to the right of $x=0$.  Where this state lies depends on whether the speed limit increases or decreases at the interface; in the example above this \"shoulder\" appears on the right.  Next is an example in which it appears on the left."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transonic rarefaction with $v_r > v_l$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_l = 0.9; rho_r = 0.2\n",
    "v_l = 1.0; v_r = 1.2\n",
    "plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 4: Transonic rarefaction-shock\n",
    "In the last two examples above we required the conditions \\eqref{vcLWR:sufficientA}-\\eqref{vcLWR:sufficientB}, ensuring that neither initial state lies fully above the other flux curve.  What if one of these conditions is violated?\n",
    "\n",
    "Here is a case where $\\rho_l$ lies in the admissible region while $\\rho_r$ lies above the value $\\max(f_l)=v_l/4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_l = 0.7; rho_r = 0.6\n",
    "v_l = 1.0; v_r = 2.\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                        v_l,v_r,connect=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now $\\rho_r$ cannot be connected to anything via a rarefaction, and cannot be connected to a state on the $f_l$ curve by a stationary jump (i.e., a horizontal line).  The only possibility then is to connect $\\rho_r$ to the admissible (left) part of the $f_r$ curve, evidently by a shock.  But then how can we connect the intermediate state back to $\\rho_l$, which lies to the left of the midline?  This is only possible by connecting from the $f_r$ curve to the $f_l$ curve at the sonic point $\\rho=1/2$.  Since this connection corresponds to the stationary jump and must be a horizontal line, this condition determines the intermediate state on the $f_r$ curve.  Then $\\rho_l$ must be connected to the sonic point by a rarefaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Physically, this corresponds to the situation in which the downstream flux is so high that the influx of cars from upstream cannot possibly keep up.  The density just downstream of the jump must be lower than that just upstream.  Thus we see a left-going rarefaction as the upstream congestion decreases since downstream traffic flux is greater; a stationary jump at the interface, where cars instantaneously accelerate due to the higher speed limit; and a right-going shock carrying a drop in density (from right to left) as the road clears out since there is less traffic arriving from upstream."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are solving a scalar problem, and yet we see three waves!  Even writing it artificially as a system of two equations (as we did above) doesn't account for this.  What's going on?\n",
    "\n",
    "In fact, this case is not so different from that of the transonic rarefaction; observe that in the solution above we have half of a transonic rarefaction; i.e. a rarefaction in which one edge coincides with the $x=0$ characteristic.  However, the stationary jump connects this rarefaction to a state with characteristic velocity greater than that of $\\rho_r$.  Thus it's impossible to impose a rarefaction on the right side, and the other half of the transonic rarefaction must instead be a shock.  We can think of this rarefaction-shock pair as a single composite wave that might be called a *transonic rarefaction-shock*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the phase plane we see that the rarefaction raises the density as far as it can go on the lower flux curve.  In other words, the flux from the left at the interface is as high as it can possibly be.  Since the characteristic speed just to the left of $x=0$ approaches zero, the condition on the stationary wave that we introduced in the introduction allows a jump to either side of the $f_r$ curve; however, only the state to the left is admissible.  Observe that if the density increased at $x=0$ (by jumping to the right side), then the resulting intermediate state would have to be connected to $\\rho_r$ by a (*right-going*) rarefaction in which all characteristics are *left-going*, which makes no sense."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mathematically, the transonic rarefactions and transonic rarefaction-shock that we have just examined are examples of *resonant waves*.  In a strictly hyperbolic system, the characteristic speeds can always be ordered strictly so that\n",
    "\n",
    "$$\n",
    "\\lambda_1 < \\lambda_2 < \\cdots < \\lambda_m.\n",
    "$$\n",
    "\n",
    "However, if two characteristic speeds coincide, then the system is non-strictly hyperbolic.  Considering the $2 \\times 2$ form of the current model \\eqref{LWR_vc:2by2}, we see that strict hyperbolicity is lost when the characteristic speed $v \\cdot (1-2\\rho)$ vanishes.  The resulting resonant interaction between the characteristic families leads to different behavior than that of a strictly hyperbolic system; for instance, the appearance of additional waves in the Riemann problem here and (in some cases, like the last solution above) an increase in the total variation of the solution.  See for example <cite data-cite=\"isaacson1992nonlinear\"><a href=\"riemann.html#isaacson1992nonlinear\">(Isaacson & Temple, 1992)</a></cite> for a discussion of such problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A fictitious solution\n",
    "In the previous example, you might wonder why there is an intermediate state with $f'(\\rho)>0$ when both $f'(\\rho_r)<0$ and $f'(\\rho_l)<0$.  In order to connect these two states without such an intermediate state, we can imagine imposing only a jump at the interface and a right-going shock or rarefaction.  Here is what we get by imposing a shock:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states = np.array([[rho_l,0.88,rho_r]])\n",
    "speeds = np.array([0,0.5])\n",
    "def reval(xi):\n",
    "    rho = np.ones((1,len(xi)))\n",
    "    return rho\n",
    "wave_types = ['contact','shock']\n",
    "fig, ax = plt.subplots(1,2,figsize=figsize)\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,v_l,v_r,states,\n",
    "                       speeds,reval,wave_types,axes=ax[1],show=False)\n",
    "riemann_tools.plot_waves(states,speeds,reval,wave_types,ax=ax[0],\n",
    "                         t_pointer=False,t=0)\n",
    "riemann_tools.plot_characteristics(reval,c,(v_l,v_r),ax[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the shock is unphysical.  The characteristics to the left of the shock are emerging from it, rather than impinging on it; this violates the entropy condition.  It's also easy to see why we can't replace that shock with a rarefaction: the rarefaction would be right-going, but all the neighboring characteristics are left-going!  Finally, this solution involves connecting to an intermediate inadmissible state, which (according to our reasoning in [the beginning of the chapter](Traffic_variable_speed.ipynb#Admissible-regions)) is incorrect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 5: Transonic rarefaction-shock with speed limit decrease\n",
    "Finally, let us consider a situation similar to Case 4, but with the situation reversed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_l = 0.7; rho_r = 0.1\n",
    "v_l = 1.0; v_r = 0.4\n",
    "traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                        v_l,v_r,connect=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Through similar reasoning, it is evident that these states must be connected through a right-going rarefaction and a left-going shock."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_riemann_traffic_vc(rho_l,rho_r,v_l,v_r);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Physically in this case the flux upstream of the jump is so great that the right section of road (with lower speed limit) cannot let it all through fast enough.  Since there is then a high density of cars at $x=0$, we might expect a rarefaction as they spread out into the uncongested section of road.  That is indeed what happens.\n",
    "\n",
    "This may seem counterintuitive at first, but is a familiar experience.  For instance, when a highway narrows to fewer lanes, it often happens that a traffic jam forms in the region just before where it narrows.  Just as one reaches the point where it narrows, the jam clears and one is able to accelerate because the next section of road is less congested (even though it is narrower)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other solutions\n",
    "We've now considered all the essential types of solutions.  *In the live notebook, you can continue to experiment with this system using the interactive setup below.  See if you can determine the correct path in the $\\rho-f(\\rho)$ plane before revealing the result (using the checkbox).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "def plot_all(rho_l,rho_r,v_l,v_r,connect=False):\n",
    "    states, speeds, reval, wave_types = \\\n",
    "        traffic_variable_speed.exact_riemann_solution(rho_l,rho_r,v_l,v_r)\n",
    "    ax = riemann_tools.plot_riemann(states,speeds,reval,wave_types,\n",
    "                       t=0.5,extra_axes=True,variable_names=[r'$\\rho$']);\n",
    "    riemann_tools.plot_characteristics(reval,c,(v_l,v_r),ax[0],\n",
    "                                       extra_lines=[[[0,0],[0,1]]])\n",
    "    ax[1].set_ylim(0,1)\n",
    "    traffic_variable_speed.phase_plane_plot(rho_l,rho_r,v_l,v_r,\n",
    "                                            axes=ax[2],connect=connect)\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(plot_all,\n",
    "         rho_l=FloatSlider(min=0.,max=1.,step=0.01,value=0.4,\n",
    "                           description=r'$\\rho_l$'),\n",
    "         rho_r=FloatSlider(min=0.,max=1.,step=0.01,value=0.7,\n",
    "                           description=r'$\\rho_r$'),\n",
    "         v_l=FloatSlider(min=0.1,max=2.,value=1.),\n",
    "         v_r=FloatSlider(min=0.1,max=2.,value=0.4),\n",
    "         connect=widgets.Checkbox(description='Connect states')\n",
    "         );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximate solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to define a numerical solver, we need to compute fluctuations (for the first-order solver) as well as waves and speeds (for the second-order corrections).\n",
    "\n",
    "### Fluctuations\n",
    "It is straightforward to compute the fluctuations by first determining the flux at the interface (i.e., the Godunov flux) $f^*$.  We can then determine the fluctuations via\n",
    "\n",
    "\\begin{align}\n",
    "    {\\mathcal A}^-\\Delta Q & = f^* - f_l \\\\\n",
    "    {\\mathcal A}^+\\Delta Q & = f_r - f^*.\n",
    "\\end{align}\n",
    "\n",
    "In fact, $f^*$ is always equal to one of $\\{f_l, f_r, v_l/4, v_r/4\\}$.  We can use the criteria discussed above to determine which of these possibilities is correct.\n",
    "\n",
    "### Second-order corrections\n",
    "The fluctuations above are the exact fluctuations from the true solution of the Riemann problem.  We could take a similar approach and include the exact wave(s) appearing in the Riemann solution in order to form second-order corrections.  This may include one or two waves (we need not include the stationary jump at the interface since it doesn't modify the value of the solution in either cell).  However, a simpler approach discussed in <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque, 2002)</a></cite> seems to work well.  We use a single wave with strength $\\rho_r-\\rho_l$ and speed equal to the Rankine-Hugoniot average speed $(c_r + c_l)/2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "def test_solver(rho_l, rho_r, v_l, v_r, riemann_solver):\n",
    "    solver = pyclaw.ClawSolver1D(riemann_solver)\n",
    "\n",
    "    solver.bc_lower[0] = pyclaw.BC.extrap\n",
    "    solver.bc_upper[0] = pyclaw.BC.extrap\n",
    "    solver.aux_bc_lower[0] = pyclaw.BC.extrap\n",
    "    solver.aux_bc_upper[0] = pyclaw.BC.extrap\n",
    "\n",
    "    x = pyclaw.Dimension(-1.0,1.0,50,name='x')\n",
    "    domain = pyclaw.Domain(x)\n",
    "    num_aux = 1\n",
    "    state = pyclaw.State(domain,solver.num_eqn,num_aux)\n",
    "\n",
    "    grid = state.grid\n",
    "    xc=grid.p_centers[0]\n",
    "\n",
    "    state.q[0,:] = rho_l*(xc<0) + rho_r*(xc>=0.)\n",
    "    state.aux[0,:] = v_l*(xc<0) + v_r*(xc>=0.)  # Speed limit\n",
    "\n",
    "    claw = pyclaw.Controller()\n",
    "    claw.tfinal = 1.0\n",
    "    claw.solution = pyclaw.Solution(state,domain)\n",
    "    claw.solver = solver\n",
    "    claw.keep_copy = True\n",
    "    claw.verbosity = 0\n",
    "\n",
    "    claw.run()\n",
    "    return xc, claw.frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "rho_l = 0.7; rho_r = 0.5\n",
    "v_l = 1.0; v_r = 1.5\n",
    "t = 0.5\n",
    "\n",
    "states, speeds, reval, wave_types = \\\n",
    "    traffic_variable_speed.exact_riemann_solution(rho_l,rho_r,\n",
    "                                                  v_l,v_r)\n",
    "x, frames = test_solver(rho_l,rho_r,v_l,v_r,riemann.traffic_vc_1D)\n",
    "\n",
    "def plot_frame(t):\n",
    "    ax = riemann_tools.plot_riemann(states, speeds, reval,\n",
    "                       wave_types, t, layout='horizontal',\n",
    "                       extra_axes=True);\n",
    "    riemann_tools.plot_characteristics(reval,c,(v_l,v_r),ax[0],\n",
    "                                       extra_lines=[[[0,0],[0,1]]])\n",
    "    rho = frames[int(t*10)].q[0,:]\n",
    "    ax[1].plot(x,rho,'-sg')\n",
    "    ax[1].set_xlim(-0.7, 0.7)\n",
    "    ax[1].legend(['Exact','Approx'],loc='upper right',fontsize=10);\n",
    "    traffic_variable_speed.phase_plane_plot(rho_l,rho_r,\n",
    "                                            v_l,v_r,axes=ax[2])\n",
    "    plt.show()\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(plot_frame,\n",
    "         t=widgets.FloatSlider(min=0,max=1,step=0.1,value=0.5));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $f$-wave solvers\n",
    "The use of the full jump $\\rho_r-\\rho_l$ in the wave discussed above may seem strange since part of this jump remains stationary at the interface (and part of it may be sent in the opposite direction!)  For systems where the flux depends explicitly on $x$, it is often simplest to use the $f$-wave approach; for a scalar problem this corresponds to using the jump in the flux, $f_r-f_l$ in place of the product $s(\\rho_r-\\rho_l)$.\n",
    "\n",
    "This technique can be used in the present system, but it has a small drawback.  At the center of a transonic rarefaction, the cell average is $1/2$ and the $f$-wave $f_r-f_l$ will be essentially zero.  This causes the limiter to activate and reduce the method locally to first order, with the consequence that the solution is inaccurate near the sonic point.  The first approach discussed above (where the wave strength is $\\rho_r-\\rho_l$) is free from this deficiency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "rho_l = 0.7; rho_r = 0.2\n",
    "v_l = 1.0; v_r = 1.0\n",
    "t = 0.5\n",
    "\n",
    "states, speeds, reval, wave_types = \\\n",
    "    traffic_variable_speed.exact_riemann_solution(rho_l,rho_r,\n",
    "                                                  v_l,v_r)\n",
    "x, frames  = test_solver(rho_l,rho_r,v_l,v_r,\n",
    "                         riemann.traffic_vc_1D)\n",
    "x ,fframes = test_solver(rho_l,rho_r,v_l,v_r,\n",
    "                         riemann.traffic_vc_fwave_1D)\n",
    "\n",
    "def plot_frame(t):\n",
    "    rho = frames[int(t*10)].q[0,:]\n",
    "    rho2 = fframes[int(t*10)].q[0,:]\n",
    "    ax = riemann_tools.plot_riemann(states, speeds, \n",
    "                       reval, wave_types, t, layout='horizontal', \n",
    "                       extra_axes=True);\n",
    "    riemann_tools.plot_characteristics(reval,c,(v_l,v_r),ax[0])\n",
    "    ax[1].plot(x,rho,'-sg')\n",
    "    ax[1].plot(x,rho2,'-or')\n",
    "    ax[1].legend(['Exact','Approx','$f$-wave'],loc='upper right',\n",
    "                 fontsize=10);\n",
    "    ax[1].set_xlim(-0.2,0.2)\n",
    "    traffic_variable_speed.phase_plane_plot(rho_l,rho_r,v_l,v_r,\n",
    "                                            axes=ax[2])\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(plot_frame,\n",
    "         t=widgets.FloatSlider(min=0,max=1,step=0.1,value=0.5));"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}