Tsunami runup example
Validation of the AnuGA implementation of the shallow water wave equation. This script sets up Okushiri Island benchmark as published at the Third International Workshop on Long-Wave Runup Models.
The validation data is available from our anuga-clinic repository and the original data is available online where a detailed description of the problem is also available.
Setup Notebook for Visualisation and Animation
We are using the format of a jupyter notebook. As such we need to setup inline matplotlib plotting and animation.
[1]:
import os
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Allow inline jshtml animations
from matplotlib import rc
rc('animation', html='jshtml')
Import ANUGA
We assume that anuga has been installed. If so we can import anuga.
[2]:
import anuga
import pathlib
# Root of the anuga_core repository (works for editable installs and Codespaces)
ANUGA_CORE = str(pathlib.Path(anuga.__file__).parent.parent)
data_dir = os.path.join(ANUGA_CORE, 'examples', 'data', 'okushiri')
print(f"Data directory: {data_dir}")
Data directory: /home/steve/anuga_core/examples/data/okushiri
The Code
This code is taken from run_okishuri.py.
First we define the location of our data files, which specify:
the incoming tsunami wave,
the bathymetry file,
the measured stage height at 3 gauge locations.
[3]:
# Incoming boundary wave (m)
boundary_filename = anuga.join(data_dir, 'Benchmark_2_input.txt')
# Digital Elevation Model (x,y,z) (m)
bathymetry_filename = anuga.join(data_dir, 'Benchmark_2_Bathymetry.xya')
# Observed timeseries (cm)
gauge_filename = anuga.join(data_dir, 'output_ch5-7-9.txt')
Load Barthymetry Data
Using in the barthymetry data provided from the workshop. We need to reshape the data and form a raster (x,y,Z).
[4]:
xya = np.loadtxt(bathymetry_filename, skiprows=1, delimiter=',')
X = xya[:,0].reshape(393,244)
Y = xya[:,1].reshape(393,244)
Z = xya[:,2].reshape(393,244)
plt.contourf(X,Y,Z, 20, cmap=plt.get_cmap('gist_earth'));
plt.title('Barthymetry')
plt.colorbar();
# Create raster tuple
x = X[:,0]
y = Y[0,:]
Zr = np.flipud(Z.T)
raster = x,y,Zr
Load Incoming Wave
We will apply an incoming wave on the left boundary. So first we load the data from the boundary_filename file.
From the data we form an interpolation functon called wave_function, which will be used to specify the boundary condition on the left.
And we also plot the function. The units of the data in the file are metres, and the scale of the experimental setup is 1 in 400.
[5]:
bdry = np.loadtxt(boundary_filename, skiprows=1)
bdry_t = bdry[:,0]
bdry_v = bdry[:,1]
import scipy.interpolate
#wave_function = scipy.interpolate.interp1d(bdry_t, bdry_v, kind='zero', fill_value='extrapolate')
#wave_function = scipy.interpolate.interp1d(bdry_t, bdry_v, kind='zero', fill_value=0.0)
wave_function = scipy.interpolate.interp1d(bdry_t, bdry_v, kind='zero',
bounds_error=False, fill_value=0.0)
t = np.linspace(0.0,25.0, 100)
plt.plot(t,wave_function(t)*400);
plt.xlabel('Seconds')
plt.ylabel('Metres')
plt.title('Incoming Wave');
Setup Domain
We define the domain for our simulation. This object encapsulates the mesh for our problem, which is defined by setting up a bounding polygon and associated tagged boundary. We use the base_resolution variable to set the maximum area of the triangles of our mesh.
At the end we use matplotlib to visualise the mesh associated with the domain.
[6]:
base_resolution = 0.01
#base_resolution = 0.0025
#base_resolution = 0.0005
# Basic geometry and bounding polygon
xleft = 0
xright = 5.448
ybottom = 0
ytop = 3.402
point_sw = [xleft, ybottom]
point_se = [xright, ybottom]
point_nw = [xleft, ytop]
point_ne = [xright, ytop]
bounding_polygon = [point_se,
point_ne,
point_nw,
point_sw]
domain = anuga.create_domain_from_regions(bounding_polygon,
boundary_tags={'wall': [0, 1, 3],
'wave': [2]},
maximum_triangle_area=base_resolution,
use_cache=False,
verbose=False)
domain.set_name('okushiri') # Name of output sww file
domain.set_minimum_storable_height(0.001) # Don't store w-z < 0.001m
#domain.set_flow_algorithm('DE1')
domain.set_flow_algorithm('DE_ader2')
print ('Number of Elements ', domain.number_of_elements)
domain.set_plotter(min_depth=0.001)
domain.triplot(linewidth = 0.4);
Number of Elements 2884
Setup Quantities
We use the raster created earlier to set the quantity called elevation. We also set the stage and the Mannings friction.
We also visualise the elevation quantity.
[7]:
domain.set_quantity('elevation',raster=raster, location='centroids')
domain.set_quantity('stage', 0.0)
domain.set_quantity('friction', 0.0025)
fig, ax, im = domain.tripcolor(facecolors = domain.elev,
edgecolors='k',
cmap='gist_earth')
plt.colorbar(im);
Setup Boundary Conditions
Excuse the verbose boundary type name Transmissive_n_momentum_zero_t_momentum_set_stage_boundary, but we use that to set the incoming wave boundary condition.
On the other boundaries we will have just reflective boundaries.
[8]:
Bts = anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, wave_function)
Br = anuga.Reflective_boundary(domain)
domain.set_boundary({'wave': Bts, 'wall': Br})
Setup Interogation Variables
We will record the stage at the 3 gauge locations and at the Monai valley.
[9]:
yieldstep = 0.5
finaltime = 25.0
nt = int(finaltime/yieldstep)+1
# area for gulleys
x1 = 4.85
x2 = 5.25
y1 = 2.05
y2 = 1.85
# indices in gulley area
x = domain.centroid_coordinates[:,0]
y = domain.centroid_coordinates[:,1]
v = np.sqrt( (x-x1)**2 + (y-y1)**2 ) + np.sqrt( (x-x2)**2 + (y-y2)**2 ) < 0.5
# Gauge and bounday locations
gauge = [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9
bdyloc = [0.00001, 2.5]
g5_id = domain.get_triangle_containing_point(gauge[0])
g7_id = domain.get_triangle_containing_point(gauge[1])
g9_id = domain.get_triangle_containing_point(gauge[2])
bc_id = domain.get_triangle_containing_point(bdyloc)
# Arrays to store data
meanstage = np.nan*np.ones((nt,))
g5 = np.nan*np.ones((nt,))
g7 = np.nan*np.ones((nt,))
g9 = np.nan*np.ones((nt,))
bc = np.nan*np.ones((nt,))
Evolve
[10]:
stage_qoi = domain.stage[v]
k = 0
for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
domain.print_timestepping_statistics()
# stage
stage_qoi = domain.stage[v]
meanstage[k] = np.mean(stage_qoi)
g5[k] = domain.stage[g5_id]
g7[k] = domain.stage[g7_id]
g9[k] = domain.stage[g9_id]
bc[k] = domain.stage[bc_id]
k = k+1
#domain.save_depth_frame()
# Read in the png files stored during the evolve loop
#domain.make_depth_animation()
Time = 0.0000 (sec), steps=0 (0s), elapsed (0s), eta (??), mem=212MB
Time = 0.5000 (sec), delta t in [0.01553112, 0.01553127] (s), steps=33 (0s), elapsed (0s), eta (4s), mem=212MB
Time = 1.0000 (sec), delta t in [0.01552375, 0.01553124] (s), steps=33 (0s), elapsed (0s), eta (3s), mem=212MB
Time = 1.5000 (sec), delta t in [0.01552365, 0.01552589] (s), steps=33 (0s), elapsed (0s), eta (3s), mem=212MB
Time = 2.0000 (sec), delta t in [0.01552457, 0.01552603] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 2.5000 (sec), delta t in [0.01552105, 0.01552455] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 3.0000 (sec), delta t in [0.01552092, 0.01553030] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 3.5000 (sec), delta t in [0.01551921, 0.01553072] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 4.0000 (sec), delta t in [0.01550542, 0.01551912] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 4.5000 (sec), delta t in [0.01548079, 0.01550530] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 5.0000 (sec), delta t in [0.01544411, 0.01548050] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 5.5000 (sec), delta t in [0.01542346, 0.01544378] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 6.0000 (sec), delta t in [0.01541570, 0.01542334] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 6.5000 (sec), delta t in [0.01540800, 0.01541559] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 7.0000 (sec), delta t in [0.01540714, 0.01540935] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=212MB
Time = 7.5000 (sec), delta t in [0.01540943, 0.01541761] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=213MB
Time = 8.0000 (sec), delta t in [0.01541780, 0.01545141] (s), steps=33 (0s), elapsed (0s), eta (2s), mem=213MB
Time = 8.5000 (sec), delta t in [0.01545216, 0.01553199] (s), steps=33 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 9.0000 (sec), delta t in [0.01548369, 0.01558370] (s), steps=33 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 9.5000 (sec), delta t in [0.01514326, 0.01548245] (s), steps=33 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 10.0000 (sec), delta t in [0.01456064, 0.01513442] (s), steps=34 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 10.5000 (sec), delta t in [0.01398208, 0.01454857] (s), steps=36 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 11.0000 (sec), delta t in [0.01358486, 0.01398114] (s), steps=37 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 11.5000 (sec), delta t in [0.01328858, 0.01358214] (s), steps=38 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 12.0000 (sec), delta t in [0.01311119, 0.01328720] (s), steps=38 (0s), elapsed (1s), eta (1s), mem=213MB
Time = 12.5000 (sec), delta t in [0.01304766, 0.01310843] (s), steps=39 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 13.0000 (sec), delta t in [0.01304748, 0.01310301] (s), steps=39 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 13.5000 (sec), delta t in [0.01310388, 0.01328438] (s), steps=38 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 14.0000 (sec), delta t in [0.01329039, 0.01357066] (s), steps=38 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 14.5000 (sec), delta t in [0.01357279, 0.01390341] (s), steps=37 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 15.0000 (sec), delta t in [0.01390763, 0.01432602] (s), steps=36 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 15.5000 (sec), delta t in [0.01433227, 0.01485791] (s), steps=35 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 16.0000 (sec), delta t in [0.01486278, 0.01510000] (s), steps=34 (0s), elapsed (1s), eta (1s), mem=214MB
Time = 16.5000 (sec), delta t in [0.01475298, 0.01495813] (s), steps=34 (0s), elapsed (1s), eta (0s), mem=214MB
Time = 17.0000 (sec), delta t in [0.01463034, 0.01475001] (s), steps=35 (0s), elapsed (1s), eta (0s), mem=215MB
Time = 17.5000 (sec), delta t in [0.01456372, 0.01463022] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 18.0000 (sec), delta t in [0.01453605, 0.01456346] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 18.5000 (sec), delta t in [0.01453170, 0.01454518] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 19.0000 (sec), delta t in [0.01454555, 0.01456639] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 19.5000 (sec), delta t in [0.01456657, 0.01462151] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 20.0000 (sec), delta t in [0.01462224, 0.01472712] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 20.5000 (sec), delta t in [0.01447670, 0.01482411] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 21.0000 (sec), delta t in [0.01155465, 0.01485605] (s), steps=38 (0s), elapsed (2s), eta (0s), mem=215MB
Time = 21.5000 (sec), delta t in [0.01159757, 0.01268591] (s), steps=41 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 22.0000 (sec), delta t in [0.01269690, 0.01312095] (s), steps=39 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 22.5000 (sec), delta t in [0.01312396, 0.01344051] (s), steps=38 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 23.0000 (sec), delta t in [0.01333954, 0.01373016] (s), steps=38 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 23.5000 (sec), delta t in [0.01354149, 0.01377782] (s), steps=37 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 24.0000 (sec), delta t in [0.01373701, 0.01392646] (s), steps=37 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 24.5000 (sec), delta t in [0.01392762, 0.01426491] (s), steps=36 (0s), elapsed (2s), eta (0s), mem=216MB
Time = 25.0000 (sec), delta t in [0.01426842, 0.01443465] (s), steps=35 (0s), elapsed (2s), eta (0s), mem=216MB
Animation Using swwfile
Read in the sww file and then iterate through the time slices to produce an animation.
[11]:
swwplotter = anuga.SWW_plotter('okushiri.sww', min_depth = 0.001)
n = len(swwplotter.time)
for k in range(n):
swwplotter.save_stage_frame(frame=k, vmin=-0.04, vmax = 0.04)
pass
swwplotter.make_stage_animation()
[11]:
View Time Series
[12]:
old_figsize = plt.rcParams['figure.figsize']
plt.rcParams['figure.figsize'] = [12, 5]
gauge = np.loadtxt(gauge_filename, skiprows=1)
gauge_t = gauge[:,0]
gauge_5 = gauge[:,1]
gauge_7 = gauge[:,2]
gauge_9 = gauge[:,3]
nt = int(finaltime/yieldstep)+1
import scipy
gauge_5_f = scipy.interpolate.interp1d(gauge_t, gauge_5, kind='zero', fill_value='extrapolate')
gauge_7_f = scipy.interpolate.interp1d(gauge_t, gauge_7, kind='zero', fill_value='extrapolate')
gauge_9_f = scipy.interpolate.interp1d(gauge_t, gauge_9, kind='zero', fill_value='extrapolate')
t = np.linspace(0.0,finaltime, nt)
tt= np.linspace(0.0,finaltime, nt)
# Flume has a 400 scaling, but gauge data is in centemeters (hence the factor 4)
plt.subplot(1,5,1)
plt.plot(t,gauge_5_f(t)*4)
plt.plot(tt,g5*400)
plt.title('Gauge 5')
plt.subplot(1,5,2)
plt.plot(t,gauge_7_f(t)*4)
plt.plot(tt,g7*400)
plt.title('Gauge 7')
plt.subplot(1,5,3)
plt.plot(t,gauge_9_f(t)*4)
plt.plot(tt,g9*400)
plt.title('Gauge 9')
plt.subplot(1,5,4)
plt.plot(t,wave_function(t)*400)
plt.plot(tt,bc*400)
plt.title('Boundary')
plt.subplot(1,5,5)
plt.plot(tt,meanstage*400)
plt.title('Runup');
plt.rcParams['figure.figsize'] = old_figsize
[ ]: