This notebook demonstrates running the Clawpack Fortran code and plotting results from a Jupyter notebook, and illustration of the effect of the entropy fix on the results.
The 1-dimensional Burgers' equation $q_t + f(q)_x = 0$ with $f(q) = \frac 1 2 q^2$ is solved in the interval $-2 \leq x \leq 2$ with periodic boundary conditions.
The entropy fix switch efix
is specified in setrun.py but can be changed below.
Have plots appear inline in notebook:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
from __future__ import print_function
Check that the CLAW environment variable is set. (It must be set in the Unix shell before starting the notebook server).
import os
try:
CLAW = os.environ['CLAW']
print("Using Clawpack from ", CLAW)
except:
print("*** Environment variable CLAW must be set to run code")
Using Clawpack from /Users/rjl/git/clawpack
from clawpack.clawutil import nbtools
nbtools.make_exe(new=True) # new=True ==> force recompilation of all code
Executing shell command: make new Done... Check this file to see output:
nbtools.make_htmls()
See the README.html file for links to input files...
First create data files needed for the Fortran code, using parameters specified in setrun.py:
nbtools.make_data(verbose=False)
Now run the code and produce plots. Specifying a label insures the resulting plot directory will persist after later runs are done below.
outdir,plotdir = nbtools.make_output_and_plots(label='1')
Executing shell command: make output OUTDIR=_output_1 Done... Check this file to see output:
Executing shell command: make plots OUTDIR=_output_1 PLOTDIR=_plots_1 Done... Check this file to see output:
View plots created at this link:
Clicking on the _PlotIndex link above, you can view an animation of the results.
After creating all png files in the _plots directory, these can also be combined in an animation that is displayed inline:
import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
The animation above was computed using the default parameter values specified in setrun.py
, which specified using the high-resolution method with the MC limiter.
See the README.html file for a link to setrun.py
.
We can adjust the parameters by reading in the default values, changing one or more, and then writing the data out for the Fortran code to use:
import setrun
rundata = setrun.setrun()
print("The efix switch is currently set to ",rundata.probdata.efix)
print("The order is currently set to ",rundata.clawdata.order)
print("The limiter is currently set to ",rundata.clawdata.limiter)
The efix switch is currently set to True The order is currently set to 2 The limiter is currently set to ['mc']
Switch to the first order method and write out the data. Then rerun the code. Note that the results are a bit more diffusive. Also note that even with the entropy fix, there is a mild discontinuity right at the sonic point. The amplitude of this would decrease with $\Delta x$, however, and this method does converge to the correct weak solution as the grid is refined.
rundata = setrun.setrun()
rundata.clawdata.order = 1
rundata.write()
outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
With the entropy fix turned off, however, there is a substantial stationary jump at $x=0$ from $q=-1$ to $q=1$ in the computed solution that does not go away if you refine the grid:
rundata = setrun.setrun()
rundata.clawdata.order = 1
rundata.probdata.efix = False
rundata.clawdata.num_cells[0] = 400
rundata.write()
outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
The high-resolution method that includes slopes and limiters smooths out the solution so that the correct weak solution is obtained in the limit. But it doesn't look as good as the original results shown above for the case when the entropy fix is used.
Here we also run out to a later time. Note that the true solution as defined in setplot.py
is only valid until the shock catches up to the rarefaction wave. If you run to a longer time, no true solution is plotted at later times.
rundata = setrun.setrun()
rundata.clawdata.tfinal = 2.
rundata.probdata.efix = False
rundata.write()
outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim