Nonconvex_scalar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonconvex scalar Conservation Laws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scalar nonlinear conservation law $q_t + f(q)_x = 0$ is said to be \"genuinely nonlinear\" over some range of states between $q_\\ell$ and $q_r$ if $f''(q) \\neq 0$ for all $q$ between these values.  This is true in particular if $f(q)$ is any quadratic function; for example, the flux function for Burgers' equation and the LWR traffic flow model are genuinely nonlinear for all $q$.  The reason this is important is that it means that the characteristic speed $f'(q)$ is either monotonically increasing or monotonically decreasing over the interval between $q_\\ell$ and $q_r$.  As discussed already in chapters [Burgers](Burgers.ipynb) and [Traffic_flow](Traffic_flow.ipynb), if $f'(q)$ is increasing as $q$ varies from $q_\\ell$ to $q_r$ then the initial discontinuity of a Riemann problem spreads out into a smooth single-valued rarefaction wave.  On the other hand if $f'(q)$ is decreasing then the discontinuity propagates as a single shock wave with speed given by the Rankine-Hugoniot condition.  Hence the Riemann solution for a genuinely nonlinear scalar equation always consists of either a single rarefaction wave or a single shock.\n",
    "\n",
    "The solution can be much more complicated if $f''(q)$ vanishes somewhere between $q_\\ell$ and $q_r$, since this means that the characteristic speed may not vary monotonically.  As we will illustrate below, for a scalar conservation law of this type the solution to the Riemann problem might consist of multiple shocks and rarefaction waves.  These waves all originate from $x=0$ and the Riemann solution is still a similarity solution $q(x,t) = Q(x/t)$ for some single-valued function $Q(\\xi)$, but determining the correct set of waves is more challenging.\n",
    "\n",
    "In this chapter, we focus on *scalar nonconvex* problems and present results computed using an elegant form of the exact solution for a general scalar conservation laws due to Osher.  At the end of the chapter we comment briefly on implications for systems of equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Osher's Solution\n",
    "\n",
    "In this chapter we use Osher's general solution to the scalar nonlinear Riemann problem (valid also for nonconvex fluxes), using the formula from <cite data-cite=\"osher1984\"><a href=\"riemann.html#osher1984\">(Osher 1984)</a></cite>.  The Riemann solution is always a similarity solution $q(x,t) = Q(x/t)$ for all $t>0$ (constant on any ray $x=\\alpha t$ eminating from the origin, for any constant $\\alpha$).  The function $Q(\\xi)$ is given by\n",
    "\n",
    "$$\n",
    "Q(\\xi) = \\begin{cases} \n",
    "    \\text{argmin}_{q_\\ell \\leq q \\leq q_r} [f(q) - \\xi q]& \\text{if} ~q_\\ell\\leq q_r,\\\\\n",
    "    \\text{argmax}_{q_r \\leq q \\leq q_\\ell} [f(q) - \\xi q]& \\text{if} ~q_r\\leq q_\\ell.\\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Recall that $\\text{argmin}_{q_\\ell \\leq q \\leq q_r} G(q)$ returns the value of $q$ for which $G(q)$ is minimized over the indicated interval, while argmax returns the value of $q$ where $G(q)$ is maximized.\n",
    "For more discussion, see also Section 16.1 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque 2002)</a></cite>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = 'svg'\n",
    "import numpy as np\n",
    "from ipywidgets import widgets, fixed\n",
    "from ipywidgets import interact\n",
    "from exact_solvers import nonconvex\n",
    "from exact_solvers import nonconvex_demos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish to examine the Python code for this chapter, please see:\n",
    "\n",
    " - [exact_solvers/nonconvex.py](exact_solvers/nonconvex.py) ...\n",
    "   [on github,](https://github.com/clawpack/riemann_book/blob/FA16/exact_solvers/nonconvex.py)\n",
    " - [exact_solvers/nonconvex_demos.py](exact_solvers/nonconvex_demos.py) ...\n",
    "   [on github.](https://github.com/clawpack/riemann_book/blob/FA16/exact_solvers/nonconvex_demos.py)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in [exact_solvers/nonconvex.py](exact_solvers/nonconvex.py) we use the function `osher_solution` to define a function `nonconvex_solutions` that evaluates this solution at a set of `xi = x/t` values.  It also computes the possibly multi-valued solution that would be obtained by tracing characteristics, for plotting purposes.  In [exact_solvers/nonconvex_demos.py](exact_solvers/nonconvex_demos.py), an additional function `make_plot_function` returns a plotting function for use in interactive widgets below.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Traffic flow\n",
    "\n",
    "First we recall that the Riemann solution for a problem with convex flux consists of a single shock or rarefaction wave.  For example, consider the flux function $f(q) = q(1-q)$ from traffic flow (with $q$ now representing the density $\\rho$ that was used in [Traffic_flow](Traffic_flow.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nonconvex_demos.demo1()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot on the left above shows a case where the solution is a rarefaction wave that can be computed by tracing characteristics.  On the right we see the case for which tracing characteristics would give an multivalued solution (as a dashed line) whereas the correct Riemann solution consists of a shock wave (solid line).  As discussed in [Traffic_flow](Traffic_flow.ipynb), the shock location is determined by the Rankine-Hugniot condition, or, equivalently, by the equal area condition discussed in [Burgers](Burgers.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison with later examples, we also plot the quadratic flux function $f(q)$ and the linear characteristic speed $f'(q)$ for this range of $q$ values.  Plotting $q$ vs. the characteristic speed shows how we can interpret each value of $q$ in the jump discontinuity (represented by the dashed vertical line in the plot on the right below) as propagating to the left or right at its characteristic speed.  Since $f'(q)$ is linear in $q$, the rarefaction wave shown above is piecewise linear. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda q: q*(1-q)\n",
    "q_left = 0.1\n",
    "q_right = 0.6\n",
    "nonconvex_demos.plot_flux(f, q_left, q_right)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Buckley-Leverett Equation\n",
    "\n",
    "The Buckley-Leverett equation for two-phase flow is described in Section 16.1.1 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque 2002)</a></cite>.  It has the non-convex flux function\n",
    "\n",
    "$$ \n",
    "f(q) = \\frac{q^2}{q^2 + a(1-q)^2}\n",
    "$$\n",
    "where $a$ is some constant, $q=1$ corresponds to pure water and $q=0$ to pure oil, in a saturated porous medium. \n",
    "\n",
    "Here we make plots like the ones above, but for the Buckley-Leverett flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.5\n",
    "f_buckley_leverett = lambda q: q**2 / (q**2 + a*(1-q)**2)\n",
    "q_left = 1.\n",
    "q_right = 0.\n",
    "\n",
    "nonconvex_demos.plot_flux(f_buckley_leverett, q_left, q_right)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again the third plot above shows $q$ on the vertical axis and $f'(q)$ on the horizontal axis (it's the middle figure turned sideways).  You can think of this as showing the characteristic velocity for each point on a jump discontinuity from $q=0$ to $q=1$ (indicated by the dashed line), and hence a triple valued solution of the Riemann problem at $t=1$ when each $q$ value has propagated this far."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The correct Riemann solution\n",
    "\n",
    "Below we show this triple-valued solution together with the correct solution to the Riemann problem, with a shock wave inserted at the appropriate point (as computed using the Osher solution defined above).  Note that for this non-convex flux function the Riemann solution consists partly of a rarefaction wave together with a shock wave.\n",
    "\n",
    "In the plot on the right, we also show the flux function $f(q)$ as a red curve and the upper boundary of the convex hull of the set of points below the graph for $q_r \\leq q \\leq q_\\ell$.  Note that the convex hull boundary follows the flux function for the set of $q$ values corresponding to the rarefaction wave and then jumps from $q\\approx 0.6$ to $q=0$, corresponding to the shock wave.  See Section 16.1 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque 2002)</a></cite> for more discussion of this construction of the Riemann solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_left = 1.\n",
    "q_right = 0.\n",
    "plot_function = nonconvex_demos.make_plot_function(f_buckley_leverett, \n",
    "                 q_left, q_right, xi_left=-2, xi_right=2)\n",
    "\n",
    "interact(plot_function, \n",
    "         t=widgets.FloatSlider(value=0.8,min=0,max=.9),\n",
    "         fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note from the plot on the left above that the triple-valued solution suggested by tracing characteristics (the dashed line) has been partially replaced by a shock wave.  By conservation, the areas of the two regions cut off by the shock must cancel out.  Moreover, the shock speed coincides with the characteristic speed at the edge of the rarefaction wave that ends at the shock.  In terms of the flux function shown by the dashed curve in the right-most figure above, we see that the shock wave connects $q_r=0$ to the point where the slope of the linear segment of the solid line (which is the shock speed, by the Rankine-Hugoniot condition) agrees with the slope of the flux function (which is the characteristic speed at this edge of the rarefaction wave).\n",
    "\n",
    "Note that the correct Riemann solution lies along the upper *convex hull* of the flux function (i.e., the convex hull of the region below the curve $f(q)$).  We will see below that this is true more generally, even for flux functions with many local maxima or minima.  If $q_\\ell > q_r$ then the upper convex hull of the flux function between these points illustrates the Riemann solution.  On the other hand, if $q_\\ell < q_r$ then the lower convex hull gives the Riemann solution, as illustrated in the next example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Swapping left and right states\n",
    "Note what happens if we switch `q_left` and `q_right` in this problem.  This now corresponds to oil on the left pushing into water on the right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "q_left = 0.\n",
    "q_right = 1.\n",
    "plot_function = nonconvex_demos.make_plot_function(f_buckley_leverett, \n",
    "                 q_left, q_right, xi_left=-2, xi_right=2)\n",
    "\n",
    "interact(plot_function, \n",
    "         t=widgets.FloatSlider(value=0.8,min=0,max=.9),\n",
    "         fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again the shock replaces the triple-valued part of the solution in such a way that mass is conserved (equal areas are cut off), and the shock speed again agrees with the limiting characteristic speed at edge of the adjacent rarefaction wave.  Since $q_\\ell < q_r$ in this case, the correct Riemann solution corresponds to the lower convex hull of the flux function, i.e. the convex hull of the set of points lying above $f(q)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Leftward flow\n",
    "The Buckley-Leverett equation as written above simulates flow from left to right.   So if we wanted to model water on the right pushing leftward into oil on the left, we must negate the flux function, as is done in the next cell. Note that this gives the same solution as our original Riemann problem, but flipped in `x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_BL_leftward = lambda q: -f_buckley_leverett(q)\n",
    "q_left = 0.\n",
    "q_right = 1.\n",
    "plot_function = nonconvex_demos.make_plot_function(f_BL_leftward, \n",
    "                 q_left, q_right, xi_left=-2, xi_right=2)\n",
    "\n",
    "interact(plot_function, \n",
    "         t=widgets.FloatSlider(value=0.8,min=0,max=.9),\n",
    "         fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sinusoidal flux\n",
    "\n",
    "As another test, we consider the flux function $f(q) = \\sin(q)$, which is also used in Example 16.1 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque 2002)</a></cite>.  First we plot the flux function $f(q)$ and also the characteristic speed $df/dq$. We also plot $q$ as a function of $df/dq$.  This again helps us visualize whether the characteristic speed is positive or negative for each value of $q$ along the jump discontinuity in the Riemann data (again indicated by the vertical dashed line), and shows that trying to solve the Riemann problem by tracing characteristics alone would lead to a multi-valued solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_sin = lambda q: np.sin(q)\n",
    "q_left = np.pi/4.\n",
    "q_right = 15*np.pi/4.\n",
    "\n",
    "nonconvex_demos.plot_flux(f_sin, q_left, q_right)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a recreation of Figure 16.4 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque 2002)</a></cite>, illustrating where shocks must be inserted to make the Riemann solution single valued.  In the live notebook you can see how the solution evolves with time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_function = \\\n",
    "    nonconvex_demos.make_plot_function(f_sin, q_left, q_right, -1.5, 1.5)\n",
    "\n",
    "interact(plot_function, \n",
    "         t=widgets.FloatSlider(value=0.8,min=0.,max=.9),\n",
    "         fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, note that the shocks in the Riemann solution correspond to linear segments of the lower convex hull the flux function $f(q)$.  This is because we chose $q_\\ell < q_r$ in this example.\n",
    "\n",
    "If we switch the states so that $q_\\ell > q_r$, then the Riemann solution corresponds to the upper convex hull of the flux function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_sin = lambda q: np.sin(q)\n",
    "q_left = 15*np.pi/4.\n",
    "q_right = np.pi/4.\n",
    "\n",
    "plot_function = \\\n",
    "    nonconvex_demos.make_plot_function(f_sin, q_left, q_right,\n",
    "                                       -1.5, 1.5)\n",
    "\n",
    "interact(plot_function, \n",
    "         t=widgets.FloatSlider(value=0.8,min=0.,max=.9),\n",
    "         fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Yet another example\n",
    "\n",
    "Here's another example that you can run in a live notebook. Observe the collection of shock and rarefaction waves that result from this.  *In the notebook you can also adjust $q_\\ell$ and $q_r$.*  Notice how the structure of the solution changes as you vary them.  What do you think the Riemann solution looks like if you switch the left and right states?  Check if you are right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda q: 0.25*(1. - q)*np.sin(1.5*q)\n",
    "\n",
    "plot_function = nonconvex_demos.make_plot_function_qsliders(f)\n",
    "\n",
    "interact(plot_function, \n",
    "         t=widgets.FloatSlider(value=0.8,min=0.,max=.9),\n",
    "         q_left=widgets.FloatSlider(value=-3.5,min=-4,max=4),\n",
    "         q_right=widgets.FloatSlider(value=3.5,min=-4,max=4),\n",
    "         fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Experiment with this and other flux functions in this notebook!  But be warned that the plotting functions, as written, may break down if you try an example with too many oscillations in the flux."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implications for systems of equations\n",
    "\n",
    "For hyperoblic systems of $m$ equations, there are $m$ different characteristic fields. For classical nonlinear systems such as the shallow water equations or the Euler equations of gas dynamics, the Riemann solution generally consists of $m$ waves, each of which is either a discontinuity or a rarefaction wave.  The fact that there is only one wave in each family results from the fact that, for these systems, each characteristic speed either varies monotonically through the wave (the field is said to be *genuinely nonlinear*) or else is constant across the wave (in which case the field is called *linearly degenerate*).  In the former case the wave is either a single shock or rarefaction wave, and in the latter case the wave is a \"contact discontinuity\", as discussed further in [Shallow_tracer](Shallow_tracer.ipynb) and [Euler](Euler.ipynb).  As in the case of scalar advection, the characteristic speed is constant across a contact discontinuity ($f'(q) \\equiv 0$ between the left and right states). The definition of genuine nonlinearity and linear degeneracy for systems of equations is given in [Shallow_tracer](Shallow_tracer.ipynb).\n",
    "\n",
    "The nonconvex equations studied in this chapter illustrate that even for a *scalar* problem the Riemann solution becomes much more complicated if the problem fails to be genuinely nonlinear (e.g. Burgers' equation or the LWR model for traffic flow) or linearly degenerate (the advection equation).  *Systems* of equations that lack the corresponding properties are more difficult still and are beyond the scope of this book, although they do arise in some important applications, such as magnetohydrodynamics (MHD)."
   ]
  }
 ],
 "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.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}