This example starts with a setrun.py
file similar to what resulted after working through the chile2010a
example, with a couple changes:
32412_notide.txt
.In this example we will first run the code as it is set up to view the gauge results.
Then we will replace this with a gauge off the coast of Peru to better observe edge waves that are trapped on the continental shelf in this region. We will also see a different way to specify output times for the animation, and how to work with gauge data in Python.
You can skip to Exercise 1
In [1]:
%pylab inline
In [2]:
from clawpack.clawutil import nbtools
import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
from IPython.display import Image
In [3]:
nbtools.make_exe(new=True,verbose=False)
In [4]:
from setrun_original import setrun
In [5]:
rundata = setrun()
rundata.write()
In [6]:
outdir,plotdir = nbtools.make_output_and_plots(verbose=False)
J.make_anim(plotdir, figno=0, figsize=(6,4))
From the _plots/_PlotIndex.html
page you can click on a figure that shows the gauge results. You should see something like the plot below.
Note that the file setplot.py
contains the code that causes this plot to show the de-tided observations at the DART location (the red curve) along with the GeoClaw solution. The first signal in the red curve is the seismic wave that is picked up by the pressure gauge at the DART location. The seismic wave travels much faster than the tsunami and arrives about 3 hours earlier.
The GeoClaw result captures the main waves but none of the details. This is partly because of the relatively coarse grid used for this example, but more because the earthquake source used in this example is not a very good representation of the actual earthquake. Better source models are available, see the documentation for some links to other data, and the this notebook for some information on how to work with this data in GeoClaw.
In [7]:
Image('_plots/gauge32412fig300.png', width=400)
In [8]:
J.make_anim(plotdir, figno=1, figsize=(6,4))
Note that there are large waves (dark red and blue) that propagate up and down along the coast. This is because the continental shelf along the coast here is fairly broad and traps waves due to reflections at the steep continental slope. This example is adapted from an Acta Numerica paper (by R.J. LeVeque, D.L. George, and M.J. Berger) that describes the GeoClaw algorithms in some detail along with some basic tsunami science. The paper and a more detailed animation of these edge waves can be found at http://www.clawpack.org/links/an11.
To explore this in more detail, let's remove the DART gauge and instead put a gauge along the coast.
In setrun.py
replace the line
rundata.gaugedata.gauges.append([32412, -86.392, -17.975, 0., 1.e10])
by
rundata.gaugedata.gauges.append([123, -77.3, -12.3, 3.5*3600, 1.e10])
(Or you could include both lines, to record time series at both gauges).
The format of the list we are appending to the list of gauges is [gaugeno, x, y, t1, t2]
. Here we use gaugeno = 123
which is an arbitray number (up to 5 digits) used to identify this gauge. x, y
are the longitude and latitude of the gauge location, and t1, t2
specify the time range over which we want to record at this gauge. In our case we start at 3.5 hours and go to 1e10
(effectively infinite, i.e. to the end of the computation).
We also want to go out to a later time to see the oscillations we expect at this gauge, and we want to output more frequently to see how they propagate in looking at the animation made from time frames. But waves don't reach this coast until about 3.5 hours into the simulation, so rather than using equally spaced times from t = 0
, let's use a different output_style
that allows setting a specific set of output times:
Change the line
clawdata.output_style = 1
to
clawdata.output_style = 2
Note that in this case the setrun.py
will execute the lines
elif clawdata.output_style == 2:
# Specify a list of output times.
clawdata.output_times = np.linspace(3.5,8,19) * 3600.
which will specify 17 output frames, every 15 minutes from 3.5 hours to 8 hours.
Below we make the same changes in the notebook and show the resulting plots.
In [9]:
rundata.gaugedata.gauges = [] # empty list
rundata.gaugedata.gauges.append([123, -77.3, -12.3, 3.5*3600, 1.e10])
rundata.clawdata.output_style = 2
rundata.clawdata.output_times = linspace(3.5,8,19) * 3600.
In [10]:
rundata.write()
outdir,plotdir = nbtools.make_output_and_plots(verbose=False)
J.make_anim(plotdir, figno=1, figsize=(6,4))
In [11]:
Image('_plots/gauge0123fig300.png', width=400)
Note that the way gauges work in GeoClaw is to output the elevation at this gauge (also momenta) at each time step, interpolating to each gauge location from the finest AMR grid available at this time in this region. In this example we started recording at this gauge at time 3.5 hours, when there was already a Level 3 grid in this region. If we had started recording earlier, this gauge might have been in a cell that was dry on a coarser grid (since it is so close to the coast) in which case the elevation of the land (possibly 10s or 100s of meters) would have been recorded at earlier times. If you see big jump discontinuities in gauge plots, this is probably what is happening. You can adjust the time over which plotting is done in setplot.py
.
The file setplot.py
can be modified to plot different things at the gauges. We don't do that here, and instead illustrate how to directly plot things in Python. The file plot_gauge123.py
in this directory also implements this. So if you type
python plot_gauge123.py
two png files should be produced with plots similar to those shown below.
First we read the plot data from the most recent GeoClaw run:
In [13]:
from setplot import setplot
plotdata = setplot()
plotdata.outdir = '_output'
g123 = plotdata.getgauge(123)
The object g123
has an attribute t
that is time in seconds, and q
that contains columns corresponding to water depth, x-momentum, y-momentum, and eta = surface elevation.
Here's how we can make some plots:
In [21]:
t = g123.t / 3600. # convert to hours
eta = g123.q[3,:] # eta = h + B (depth plus bathymetry)
plot(t,eta)
xlim(3,8)
xlabel('Hours since earthquake')
ylabel('meters')
title('Sea surface elevation at gauge 123')
In [23]:
h = g123.q[0,:] # water depth
u = g123.q[1,:] / h # u = east-west velocity
v = g123.q[2,:] / h # v = north-south velocity
plot(t,u,'g',label='u')
plot(t,v,'r',label='v')
legend()
title('Velocities at gauge 123')
xlabel('Hours since earthquake')
ylabel('meters per second')
In [ ]: