Acoustics in one dimension

This Jupyter notebook can be found in collection of Clawpack apps as the file $CLAW/apps/notebooks/classic/acoustics_1d_example1/acoustics_1d_example1.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. This example is based on the one in $CLAW/classic/examples/acoustics_1d_example1.

The 1-dimensional acoustics equations $q_t + Aq_x = 0$ is solved in the interval $-1 \leq x \leq 1$ with extrapolation boundary conditions.

The density rho and bulk modulus K are specified in setrun.py but can be changed below.

Version

Animation revised to run with v5.7.0

Have plots appear inline in notebook:

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
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).

In [3]:
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:

In [4]:
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.

In [5]:
def show_anim(anim):
    html_version = HTML(anim.to_jshtml())
    #html_version = HTML(anim.to_html5_video())
    return html_version

Compile the code:

In [6]:
nbtools.make_exe(new=True)  # new=True ==> force recompilation of all code

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

Make documentation files:

In [7]:
nbtools.make_htmls()

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

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:

In [8]:
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.

In [9]:
outdir,plotdir = nbtools.make_output_and_plots(label='1')

Executing shell command:   make output OUTDIR=nb_output_1
Done...  Check this file to see output:
run_output_1.txt
Executing shell command:   make plots OUTDIR=nb_output_1 PLOTDIR=nb_plots_1
Done...  Check this file to see output:
plot_output_1.txt
View plots created at this link:
nb_plots_1/_PlotIndex.html
Or you may have to copy and paste:
file:///Users/rjl/clawpack_src/clawpack_master/apps/notebooks/classic/acoustics_1d_example1/nb_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:

In [10]:
anim = animation_tools.animate_from_plotdir(plotdir);
show_anim(anim)

Using figno = 1
Out[10]:

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:

In [11]:
import setrun
rundata = setrun.setrun()

print("The left boundary condition is currently set to ",rundata.clawdata.bc_lower)
print("The right boundary condition is currently set to ",rundata.clawdata.bc_upper)

The left boundary condition is currently set to  ['extrap']
The right boundary condition is currently set to  ['extrap']

Periodic boundary conditions

Change the boundary conditions and write out the data. Then rerun the code.

In [12]:
rundata.clawdata.bc_lower = ['periodic']
rundata.clawdata.bc_upper = ['periodic']

# also extend the time over which the solution is computed and plotted:
rundata.clawdata.num_output_times = 40
rundata.clawdata.tfinal = 2.0

rundata.write()

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

anim = animation_tools.animate_from_plotdir(plotdir);
show_anim(anim)

Using figno = 1
Out[12]:

Reflecting wall boundary conditions

Change the boundary conditions and write out the data. Then rerun the code.

In [13]:
rundata.clawdata.bc_lower = ['wall']
rundata.clawdata.bc_upper = ['wall']

rundata.write()

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

Using figno = 1
Out[13]:

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.

  • Change the density rho or the bulk modulus K and see what effect this has.

In [14]:
def print_material():
    rho = rundata.probdata.rho
    K = rundata.probdata.K
    Z = sqrt(K*rho)
    c = sqrt(K/rho)
    print("The density rho = %g and bulk modulus %g give" % (rho,K))
    print("  speed of sound c = %g" % c)
    print("  impedance Z = %g" % Z)
print_material()

The density rho = 1 and bulk modulus 4 give
  speed of sound c = 2
  impedance Z = 2
In [15]:
rundata.probdata.rho = 4.
rundata.write()
print_material()

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

The density rho = 4 and bulk modulus 4 give
  speed of sound c = 1
  impedance Z = 4
Using figno = 1
Out[15]:
In [ ]: