Scalar equation with exponential flux

In this notebook we solve the scalar equation $q_t + f(q)_x = 0$ with flux function $f(q) = e^q$ and initial data

$ q(x,0) = \left\{ \begin{array}{ll} 2 & \text{if}~ 0\leq x \leq 1,\\ 0 & \text{otherwise.} \end{array}\right. $

Sample notebook adapted from Homework 3 solutions.

Note: Getting this notebook to work might require updating Clawpack

    cd $CLAW/clawutil
    git checkout master
    git pull

    cd $CLAW/visclaw
    git checkout master
    git pull



In [ ]:

    
%pylab inline

The solution initially consists of a shock wave at $x=1$ moving at speed $s_0 = (f(2)-f(0))/(2-0) = (e^2-1)/2$ and a centered rarefaction at $x=0$. In the rarefaction wave the solution takes the form $q_{rare}(x,t) = \log(x/t)$. The leading edge of the rarefaction wave is moving at speed $f'(2) = e^2$ and catches up to the shock at time $t_s = 2/(e^2 + 1)$. For $t \leq t_s$, the solution is

$q(x,t) = \left\{ \begin{array}{ll} \log(x/t) & \text{if}~ t \leq x \leq e^2 t,\\ 2 & \text{if}~ e^2 t\leq x \leq 1+s_0t,\\ 0 & \text{otherwise.} \end{array}\right. $

For $t > t_s$, the rarefaction is interacting with the shock and the shock speed varies. If $x_s(t)$ is the shock position, then it satisfies and ODE determined by the fact that the shock speed is given by the Rankine-Hugoniot condition at each instant in time:

$x_s' = \frac{f(q_{rare}(x_s(t),t)) - f(0)}{q_{rare}(x_s(t),t) - 0} = \frac{x_s(t)/t - 1}{\log(x_s(t)/t)}$

with the intial condition $x_s(t_s) = 1 + s_0 t_s$.

This ODE can be solved in terms of the Lambert W-function $W(z)$ defined to be the solution to $We^W = z$:

$x_s(t) = \frac{2-t}{W((2-t)/et)}.$

This can also be determined by conservation since $x_s(t)$ must be the value for which $\int_t^{x_s(t)} \log(x/t) dx = 2$. Evaluating this integral gives $x_s \log(x_s) - x_s - x_s\log(t) + t = 2$ and some manipulations allow this nonlinear equation to be solved in terms of the W-function as above (or it could be solved numerically, e.g. with the Python function scipy.optimize.fsolve).

Assuming a way to evaluate $x_s(t)$, the solution for $t > t_s$ can be written as

$q(x,t) = \left\{ \begin{array}{ll} \log(x/t) & \text{if}~ t \leq x \leq x_s(t),\\ 0 & \text{otherwise.} \end{array}\right. $

Evaluating from this form requires solving a nonlinear equation to determine $x_s(t)$ for each $t > t_s$, or having an implementation of the Lambert W-function. In Python, one can use scipy.special.lambertw. This is implemented in the function below.



In [ ]:

    
def q(x,t):
    from scipy.special import lambertw
    q = zeros(x.shape)
    q = where(x>t, log(x/t), q)  # rarefaction
    ts = 2./(exp(2.)+1)  # time rarefaction catches shock
    if t<ts:
        q = where(x<exp(2)*t, q, 2.) # flat top
        q = where(x < 1+t*(exp(2)-1)/2., q, 0.)  # cut off at shock
    else:
        xs = (2-t) / (lambertw((2-t)/(exp(1.)*t)))
        q = where(x < xs, q, 0.) # cut off at shock
    return q

Plot the solution for several times between 0 and 1:



In [ ]:

    
xx = linspace(-1,5,500)
for t in linspace(0,1,11):
    qq = q(xx,t)
    plot(xx,qq)
ylim(-0.5,2.5)

Another way to compute the true solution on a set of x points

We can avoid having to use the W-function by simply evaluating the rarefaction wave function everywhere to the right of $x=t$ and then reset the value of $q$ to zero at points where the integral would evaluate to a value greater than 2, e.g. when $x\log(x/t) - x + t > 2$:



In [ ]:

    
def q2(x,t):
    ts = 2./(exp(2.)+1)  # time rarefaction catches shock
    q = zeros(x.shape)
    q = where(x>t, log(x/t), q)  # rarefaction for x>t, ignoring shock
    if t<ts:
        q = where(x<exp(2)*t, q, 2.) # flat top
        q = where(x < 1+t*(exp(2)-1)/2., q, 0.)  # cut off at shock
    else:
        # cut off at shock, ignoring log values at invalid arguments:
        q = where(logical_and(x>0, x*log(x/t) - x + t > 2), 0., q)
    return q



In [ ]:

    
xx = linspace(-1,5,500)
for t in linspace(0,1,11):
    qq = q2(xx,t)
    plot(xx,qq)
ylim(-0.5,2.5)

This latter approach is used in the setplot.py routine used in the Clawpack code below.

Clawpack solution

Check that the CLAW environment variable is set. (It must be set in the Unix shell before starting the notebook server).



In [ ]:

    
try:
    import clawpack
    location = clawpack.__file__.replace('clawpack/__init__.pyc','')
    print "Using Clawpack from ",location
except:
    print "*** Problem importing Clawpack -- check if environment variable set"

Import some modules needed below...



In [ ]:

    
import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
from clawpack.clawutil import nbtools
from IPython.display import FileLink

Compile the Fortran code:



In [ ]:

    
nbtools.make_exe(new=True)

Make documentation files:



In [ ]:

    
nbtools.make_htmls()

Run code, and plot results...



In [ ]:

    
outdir,plotdir = nbtools.make_output_and_plots(verbose=True)

Add the animation to the notebook:



In [ ]:

    
anim = J.make_anim(plotdir, figno=1, figsize=(8,6))
anim



In [ ]: