{
"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_l$ 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_l$ and $q_r$. As discussed already in chapters [Burgers.ipynb](Burgers.ipynb) and [Traffic_flow.ipynb](Traffic_flow.ipynb), if $f'(q)$ is increasing as $q$ varies from $q_l$ 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_l$ 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 right set of waves is more challenging. Below we present results computed using an elegant form of the solution to general scalar conservation laws due to Osher."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linear degeneracy\n",
"\n",
"The opposite extreme of genuine nonlinearity is \"linear degeneracy\". The scalar advection equation has $f(q) = aq$ with constant characteristic speed $f'(q)\\equiv a$ and so $f''(q) \\equiv 0$ for all values of $q$. Since it is a linear equation, it naturally fails to be genuinely nonlinear over any range of $q$ values.\n",
"\n",
"Recall that in the Riemann solution for the scalar advection equation (discussed in [Advection.ipynb](Advection.ipynb)) the characteristics are parallel to the propagating discontinuity in the $x$-$t$ plane. They are neither impinging on the discontinuity (as in a shock wave, where the characteristic speed decreases from $q_l$ to $q_r$) nor are they spreading out (as they would in a rarefaction wave with $f'(q)$ increasing). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implications to 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 (if the field is genuinely nonlinear) or else is constant across the wave (if the field is 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.ipynb](Shallow_tracer.ipynb) and [Euler.ipynb](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.ipynb](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 or linearly degenerate. *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 applcations, such as magnetohydrodynamics (MHD)."
]
},
{
"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 non-convex fluxes), using the formula from (Osher 1984). 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_l \\leq q \\leq q_r} [f(q) - \\xi q]& \\text{if} ~q_l\\leq q_r,\\\\\n",
" \\text{argmax}_{q_r \\leq q \\leq q_l} [f(q) - \\xi q]& \\text{if} ~q_r\\leq q_l.\\\\\n",
"\\end{cases}\n",
"$$\n",
"\n",
"Recall that $\\text{argmin}_{q_l \\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",
"\n",
"For more discussion, see also Section 16.1 of (LeVeque 2002)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"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",
" - [exact_solvers/nonconvex_demos.py](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 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.ipynb](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)."
]
},
{
"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 (LeVeque 2002). 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",
"Consider the Riemann problem for water intruding into oil, with $q_l=1$ and $q_r=0$."
]
},
{
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot the flux and its derivative"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide"
]
},
"outputs": [],
"source": [
"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_l$. 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 (LeVeque 2002) 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 along 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 solid line (the shock speed) agrees with the slope of the flux function (the characteristic speed at this edge of the rarefaction wave). The correct Riemann solution lies along the upper convex hull of the flux function."
]
},
{
"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. Since $q_l < q_r$ in this case, the correct Riemann solution corresponds to the lower convex hull of the flux function."
]
},
{
"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": [
"### Leftward flow\n",
"Note that 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, the flux function $f(q) = \\sin(q)$ is used in Example 16.1 of (LeVeque 2002).\n",
"\n",
"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 (LeVeque 2002), 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 boundary of the convex hull of the set of points that lie above the flux function $f(q)$. This is because we chose $q_l < q_r$ in this example.\n",
"\n",
"If we switch the states so that $q_l > q_r$, then the Riemann solution corresponds to the upper boundary of the convex hull of the set of points that lie below 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. Note the collection of shock and rarefaction waves that result from this. In the notebook you can adjust $q_l$ 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.\n",
"\n",
"Note this example doesn't work in the static html version due to sliders for $q_l$ and $q_r$."
]
},
{
"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 please note that the plotting functions, as written, may break down if you try an example with too many oscillations in the flux."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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": 1
}