Chile2010b -- Example of gauges

This Jupyter notebook can be found in collection of Clawpack apps as the file $CLAW/apps/notebooks/geoclaw/chile2010b/chile2010b.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 example starts with a setrun.py file similar to what resulted after working through the $CLAW/apps/notebooks/geoclaw/chile2010a example, with a couple changes:

  • The Level 3 grids are now refined by a factor of 6 relative to Level 2.
  • A gauge has been added at the location of DART buoy 32412. DART (Deep Ocean Assessment and Reporting of Tsunamis) buoys transmit data from pressure gauges on the sea floor that can measure the water pressure accurately enough to detect long wave disturbances in the water surface level. Data from this site has been de-tided and also saved in this directory, in the file 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.

Version

Animation revised 2020-04-09 to run with v5.7.0

Notebook setup

You can skip to Exercise 1

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
from clawpack.clawutil import nbtools
from clawpack.visclaw import animation_tools
from IPython.display import HTML, Image
In [3]:
def show_anim(anim):
    html_version = HTML(anim.to_jshtml())
    #html_version = HTML(anim.to_html5_video())
    return html_version
In [4]:
nbtools.make_exe(new=True,verbose=False)
In [5]:
from setrun_original import setrun
In [6]:
rundata = setrun()
rundata.write()

Exercise 1

Run the code via make .plots and observe the results. The animation titled "Full Domain" should look like this:

In [7]:
outdir,plotdir = nbtools.make_output_and_plots(verbose=False)
anim = animation_tools.animate_from_plotdir(plotdir, figno=0)
show_anim(anim)
Out[7]:

Gauge plot

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 [8]:
Image('_plots/gauge32412fig300.png', width=400)
Out[8]:

Edge waves on the coast of Peru

The setplot.py file in this directory plots two figures each time frame. One shows the whole domain, the animation above. The second plot sets the axis limits so that it zooms in on the coast of Peru. The animation produces should look like what is seen below.

In [9]:
anim = animation_tools.animate_from_plotdir(plotdir, figno=1)
show_anim(anim)
Out[9]:

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.

Exercise 2.

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

Changing the output times

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 [10]:
# Change output times:
rundata.clawdata.output_style = 2
rundata.clawdata.output_times = linspace(3.5,8,19) * 3600.
In [11]:
# Change gauge location:
from clawpack.amrclaw.data import GaugeData
rundata.replace_data("gaugedata", GaugeData())  # clear old gauge data
rundata.gaugedata.gauges = []   # empty list
rundata.gaugedata.gauges.append([123, -77.3, -12.3, 3.5*3600, 1.e10])
In [12]:
rundata.write()
In [13]:
rundata.write()
outdir,plotdir = nbtools.make_output_and_plots(verbose=False)
In [14]:
anim = animation_tools.animate_from_plotdir(plotdir, figno=1)
show_anim(anim)
Out[14]:
In [15]:
Image('_plots/gauge0123fig300.png', width=400)
Out[15]:

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.

Viewing gauges in Python

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 [16]:
from setplot import setplot
plotdata = setplot()
plotdata.outdir = '_output'
g123 = plotdata.getgauge(123)

*** Could not load DART data file
Read in gauge 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 [17]:
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')
Out[17]:
Text(0.5, 1.0, 'Sea surface elevation at gauge 123')
In [18]:
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')
Out[18]:
Text(0, 0.5, 'meters per second')
In [ ]: