Advection in one dimension¶

This Jupyter notebook can be found in collection of Clawpack apps as the file $CLAW/apps/notebooks/classic/advection_1d/advection_1d.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.

This notebook demonstrates running the Clawpack Fortran code and plotting results from a Jupyter notebook, and illustration of the behavior of various methods and limiters on the same problem.

The 1-dimensional advection equation $q_t + uq_x = 0$ is solved in the interval $0 \leq x \leq 1$ with periodic boundary conditions.

The constant advection velocity $u$ is specified in setrun.py but can be changed below.

Version¶

Animation revised to run with v5.7.0

Have plots appear inline in notebook:

%pylab inline

Populating the interactive namespace from numpy and matplotlib

from __future__ import print_function  # so notebook works in Python2 or Python3

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/clawpack_src/clawpack_master

Module with functions used to execute system commands and capture output:¶

from clawpack.clawutil import nbtools
from clawpack.visclaw import animation_tools
from IPython.display import HTML

Inline animations:¶

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

Executing shell command:   make new
Done...  Check this file to see output:

Make documentation files:¶

nbtools.make_htmls()

See the README.html file for links to input files...

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.

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:

Or you may have to copy and paste:
file:///Users/rjl/clawpack_src/clawpack_master/apps/notebooks/classic/advection_1d/_plots_1/_PlotIndex.html

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);
show_anim(anim)

Using figno = 1

Adjust some parameters to explore the methods in Clawpack¶

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 order is currently set to ",rundata.clawdata.order)
print("The limiter is currently set to ",rundata.clawdata.limiter)

The order is currently set to  2
The limiter is currently set to  ['mc']

First order upwind method¶

Switch to the first order method and write out the data. Then rerun the code. Note that the results are much more diffusive.

rundata.clawdata.order = 1
rundata.write()

outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = animation_tools.animate_from_plotdir(plotdir);
show_anim(anim)

Using figno = 1

Lax-Wendroff method¶

Switch back to the second-order method and use no limiter so the high-resolution method reduces to Lax-Wendroff.

A list of 1 limiter is required since this is a scalar problem with only 1 wave in the Riemann solution, so we set the limiter to ['none'].

Note that the Lax-Wendroff method is dispersive and introduces spurious oscillations, particularly near the discontinuities in the solution.

rundata.clawdata.order = 2
rundata.clawdata.limiter = ['none']
rundata.write()

outdir, plotdir = nbtools.make_output_and_plots(label='LW',verbose=False)
anim = animation_tools.animate_from_plotdir(plotdir);
show_anim(anim)

Using figno = 1

Things to try:¶

Use print rundata.clawdata to see what parameters can be set. Also print rundata.probdata to see what problem-specific paramaters are available.
Try the Lax-Wendroff method with 400 grid cells rather than the 200 used by default. Do the oscillations decrease in amplitude? The parameter to change is rundata.clawdata.num_cells, which is a list of the number of cells in each direction (only a single entry for this 1-dimensional problem).
Set rundata.clawdata.tfinal = 20. so that the 20 frames displayed are at time $1, 2, \ldots, 20$ in order to investigate how the solution degrades over time. Note that at these times the true solution agrees with the initial data because of the periodic boundary conditions.
Try upwind and various limiters: none, minmod, mc, superbee.
The size of the timesteps used is controlled by the parameter rundata.clawdata.cfl_desired specifying the desired Courant number. By default it is 0.9. What happens if you change it to 0.1? Do any of the above methods perform better? What if you change it to 1.0?