%matplotlib inline

%config InlineBackend.figure_format = 'svg'
import numpy as np
from ipywidgets import widgets, fixed
from ipywidgets import interact
from exact_solvers import nonconvex
from exact_solvers import nonconvex_demos

"""
Explanation: Nonconvex scalar Conservation Laws
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 and Traffic_flow, 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.
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.
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.
Osher's Solution
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
$$
Q(\xi) = \begin{cases} 
    \text{argmin}{q\ell \leq q \leq q_r} [f(q) - \xi q]& \text{if} ~q_\ell\leq q_r,\
    \text{argmax}{q_r \leq q \leq q\ell} [f(q) - \xi q]& \text{if} ~q_r\leq q_\ell.\
\end{cases}
$$
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.
For more discussion, see also Section 16.1 of <cite data-cite="fvmhp"><a href="riemann.html#fvmhp">(LeVeque 2002)</a></cite>.
End of explanation
"""


nonconvex_demos.demo1()

"""
Explanation: If you wish to examine the Python code for this chapter, please see:

exact_solvers/nonconvex.py ...
   on github,
exact_solvers/nonconvex_demos.py ...
   on github.

Note that in 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, an additional function make_plot_function returns a plotting function for use in interactive widgets below.  
Traffic flow
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).
End of explanation
"""


f = lambda q: q*(1-q)
q_left = 0.1
q_right = 0.6
nonconvex_demos.plot_flux(f, q_left, q_right)

"""
Explanation: 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, the shock location is determined by the Rankine-Hugniot condition, or, equivalently, by the equal area condition discussed in Burgers.
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.
End of explanation
"""


a = 0.5
f_buckley_leverett = lambda q: q**2 / (q**2 + a*(1-q)**2)
q_left = 1.
q_right = 0.

nonconvex_demos.plot_flux(f_buckley_leverett, q_left, q_right)

"""
Explanation: Buckley-Leverett Equation
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
$$ 
f(q) = \frac{q^2}{q^2 + a(1-q)^2}
$$
where $a$ is some constant, $q=1$ corresponds to pure water and $q=0$ to pure oil, in a saturated porous medium. 
Here we make plots like the ones above, but for the Buckley-Leverett flux.
End of explanation
"""


q_left = 1.
q_right = 0.
plot_function = nonconvex_demos.make_plot_function(f_buckley_leverett, 
                 q_left, q_right, xi_left=-2, xi_right=2)

interact(plot_function, 
         t=widgets.FloatSlider(value=0.8,min=0,max=.9),
         fig=fixed(0));

"""
Explanation: 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.
The correct Riemann solution
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.
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.
End of explanation
"""


q_left = 0.
q_right = 1.
plot_function = nonconvex_demos.make_plot_function(f_buckley_leverett, 
                 q_left, q_right, xi_left=-2, xi_right=2)

interact(plot_function, 
         t=widgets.FloatSlider(value=0.8,min=0,max=.9),
         fig=fixed(0));

"""
Explanation: 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).
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.
Swapping left and right states
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.
End of explanation
"""


f_BL_leftward = lambda q: -f_buckley_leverett(q)
q_left = 0.
q_right = 1.
plot_function = nonconvex_demos.make_plot_function(f_BL_leftward, 
                 q_left, q_right, xi_left=-2, xi_right=2)

interact(plot_function, 
         t=widgets.FloatSlider(value=0.8,min=0,max=.9),
         fig=fixed(0));

"""
Explanation: 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)$.
Leftward flow
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.
End of explanation
"""


f_sin = lambda q: np.sin(q)
q_left = np.pi/4.
q_right = 15*np.pi/4.

nonconvex_demos.plot_flux(f_sin, q_left, q_right)

"""
Explanation: Sinusoidal flux
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.
End of explanation
"""


plot_function = \
    nonconvex_demos.make_plot_function(f_sin, q_left, q_right, -1.5, 1.5)

interact(plot_function, 
         t=widgets.FloatSlider(value=0.8,min=0.,max=.9),
         fig=fixed(0));

"""
Explanation: 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.
End of explanation
"""


f_sin = lambda q: np.sin(q)
q_left = 15*np.pi/4.
q_right = np.pi/4.

plot_function = \
    nonconvex_demos.make_plot_function(f_sin, q_left, q_right,
                                       -1.5, 1.5)

interact(plot_function, 
         t=widgets.FloatSlider(value=0.8,min=0.,max=.9),
         fig=fixed(0));

"""
Explanation: 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.
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:
End of explanation
"""


f = lambda q: 0.25*(1. - q)*np.sin(1.5*q)

plot_function = nonconvex_demos.make_plot_function_qsliders(f)

interact(plot_function, 
         t=widgets.FloatSlider(value=0.8,min=0.,max=.9),
         q_left=widgets.FloatSlider(value=-3.5,min=-4,max=4),
         q_right=widgets.FloatSlider(value=3.5,min=-4,max=4),
         fig=fixed(0));

"""
Explanation: Yet another example
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.
End of explanation
"""


