Demo of running AMRClaw Fortran code and plotting results¶
This Jupyter notebook can be found in collection of Clawpack apps as the file $CLAW/apps/notebooks/amrclaw/advection_2d_square/amrclaw_advection_2d_square.ipynb.
To run this notebook, install Clawpack, and clone the apps repository.
A static view of this and other notebooks can be found in the Clawpack Gallery of Jupyter notebooks.
The is adapted from $CLAW/amrclaw/examples/advection_2d_square to illustrate the use of IPython notebooks.
The setrun.py and setplot.py files are taken from that example but some parameters are modified in the notebook below.
The 2-dimensional advection equation $q_t + uq_x + vq_y = 0$ is solved in the unit square with periodic boundary conditions. The constant advection velocities $u$ and $v$ are specified in setrun.py but can be changed below.
Have plots appear inline in notebook:
%pylab inline
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")
Some tools for running and plotting in the notebook:¶
from clawpack.clawutil import nbtools
from clawpack.visclaw import animation_tools
from IPython.display import HTML
Using anim.to_jshtml() gives animations similar to what you see in the html files if you do make plots, but you may prefer the anim.to_html5_video() option. See the matplotlib.animation.Animation documentation for more information, also on how to save an animation as a separate file.
def show_anim(anim):
html_version = HTML(anim.to_jshtml())
#html_version = HTML(anim.to_html5_video())
return html_version
Compile the code:¶
nbtools.make_exe(new=True) # new=True ==> force recompilation of all code
Make documentation files:¶
nbtools.make_htmls()
Run the code and plot results using the setrun.py and setplot.py files in this directory:¶
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.
from importlib import reload
reload(nbtools)
outdir,plotdir = nbtools.make_output_and_plots(label='1',verbose=True)
Display the animation inline:¶
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:
anim = animation_tools.animate_from_plotdir(plotdir, figno=0);
show_anim(anim)
Illustrate how to adjust some parameters and rerun the code:¶
First read in the current rundata. (See the README.html file for a link to setrun.py).
import setrun
rundata = setrun.setrun()
Here we adjust the mesh to be $40 \times 40$ and set two gauges at $(0.6, 0.4)$ and $(0.6,0.8)$.
rundata.clawdata.num_cells = [40,40]
rundata.gaugedata.gauges = [[1, 0.6, 0.4, 0., 10.], [2, 0.6, 0.8, 0., 10.]]
The advection velocities can also be changed (in setrun.py, $(u,v) = (0.5,1)$.
rundata.probdata.u = -1.0
rundata.probdata.v = 1.0
Change the limiter if you want to experiment with other methods. A list of 1 limiter is required since this is a scalar problem with only 1 wave in the Riemann solution.
Here we set the limiter to ['none'] to see the oscillations that arise if Lax-Wendroff is used with no limiter.
rundata.clawdata.limiter = ['none']
Write the data out (for the Fortran code to read in) and run the code:¶
rundata.write()
outdir, plotdir = nbtools.make_output_and_plots(label='2', verbose=False)
Changing plot parameters and plotting inline¶
In addition to producing a _plots directory using the parameters in setplot.py, plot parameters can be changed and plots produced inline as the examples below illustrate...
import setplot
plotdata = setplot.setplot()
plotdata.outdir = outdir
Plot one frame of the solution:¶
plotdata.plotframe(2)
Adjust some plot parameters:¶
plotdata.showitems() # shows the currently defined figures, axes, items
Change the color limits for the pcolor plot to show the overshoots and undershoots with Lax-Wendroff. Also suppress showing the other plots...
plotitem = plotdata.getitem('ITEM1','AXES1','pcolor')
plotitem.pcolor_cmin = -0.2
plotitem.pcolor_cmax = 1.2
plotdata.getfigure('contour').show = False
plotdata.getfigure('cells').show = False
plotdata.plotframe(2)
Explore the gauge output:¶
The overshoots and undershoots with Lax-Wendroff are even more visible in the gauge output...
for gaugeno in [1,2]:
gauge = plotdata.getgauge(gaugeno)
q = gauge.q[0,:]
t = gauge.t
qmax = q.max()
tmax = t[q.argmax()]
print("Gauge %s at %s has a maximum value of q = %7.5f at t = %7.5f" \
% (gaugeno, gauge.location, qmax, tmax))
plot(t,q,label="Gauge %s" % gaugeno)
legend()
xlim(0,3)
ylim(-0.1,1.1)
xlabel('time')
ylabel('q at gauge')
title("Gauge output")
