Acoustics in one dimension

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 notebook can be found in $CLAW/apps/notebooks/classic/acoustics_1d_example1/acoustics_1d_example1.ipynb

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.

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/git/clawpack

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

from clawpack.clawutil import nbtools

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:

compile_output.txt

Make documentation files:

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:

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:

run_output_1.txt

Executing shell command:   make plots OUTDIR=_output_1 PLOTDIR=_plots_1
Done...  Check this file to see output:

plot_output_1.txt

View plots created at this link:

_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:

import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim

– +

Once Loop Reflect

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 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.

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 = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim

– +

Once Loop Reflect

Reflecting wall boundary conditions

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

rundata.clawdata.bc_lower = ['wall']
rundata.clawdata.bc_upper = ['wall']

rundata.write()

outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim

– +

Once Loop Reflect

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.

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

rundata.probdata.rho = 4.
rundata.write()
print_material()

outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim

The density rho = 4 and bulk modulus 4 give
  speed of sound c = 1
  impedance Z = 4

– +

Once Loop Reflect