Enet et al benchmark
From Gerris
This benchmark for simulations of landslide-generated tsunamis using Gerris is based on these articles:
Enet, F., Grilli, S.T., Watts, P. - Laboratory experiments for tsunamis generated by underwater landslides: Comparison with numerical modeling
- Proc. of the 13th Offshore and Polar Engrg. Conf., ISOPE03, Honolulu, Hawaii 3:372--379, 2003
- Bibtex
Enet, F., Grilli, S.T. - Experimental study of tsunami generation by three-dimensional rigid underwater landslides
- Journal of Waterway, Port, Coastal, and Ocean Engineering 133:442, 2007
- Bibtex
Contents |
Generating the landslide body geometry
The moving boundary implementation GfsSolidMoving cannot deal (yet) with surfaces represented analytically (unlike GfsSolid). We need to generate a GTS surface instead. This can be done using the following script:
# generate a delaunay triangulation of the basic profile
awk -f enet.awk | delaunay > enet.gts
# extract a contour level at z = 0
gts2oogl -G -n -I 0 < enet.gts > iso
# at the points of this contour to the original dataset and retriangulate
( awk -f enet.awk ; cat iso ) | sed 's/4000/4367/' | delaunay -d > enet.gts
The script is a bit complicated because we need to circumvent another limitation in GfsSolidMoving: we want the solid to intersect (slightly) with the fixed slope (so that there is no gap), but in that case GfsSolidMoving will not move the solid vertices which are (just) outside of the fluid; this could give large deformations at the base of the solid as the other vertices (inside the fluid) start moving. The trick to avoid this is to add vertices which are very close to the boundary (but inside the fluid).
The awk script 'enet.awk' implements the analytical formula (1) of Enet and Grilli, 2007 as:
function sech(x) {
return 2./(exp(x) + exp(-x));
}
function acosh(x) {
return log(x + sqrt(x + 1.)*sqrt(x - 1.));
}
BEGIN {
# See Enet and Grilli, 2007, eqn (1)
w = 0.680
b = 0.395
T = 0.082
epsilon = 0.717
C = acosh(1./epsilon)
Kw = 2.*C/w
Kb = 2.*C/b
print 4000,0,0
for (x = -0.4; x <= 0.4; x += 0.01) {
for (y = -0.25; y <= 0.25; y += 0.01) {
z = T/(1. - epsilon)*(sech(Kw*x)*sech(Kb*y) - epsilon);
print x, y, z;
}
}
}
Help on 'delaunay' and 'gts2oogl' commands is displayed when using the '-h' option.
Numerical "wave gauge"
The experimental setup uses wave gauges to measure the free-surface height. In our 3D simulation we do not have direct access to the "interface height" (it may not even be well defined e.g. for breaking waves). We need to create a "numerical wave gauge". We use GfsOutputLocation to output model variables along a "vertical" transect along the "wire" of the wave gauge. We then filter the values to keep only those in cells containing the free surface. We then output the "height" of the center of the corresponding interface fragment. This is done using the following 'distance.awk' script:
BEGIN {
theta = 15.*3.14159265359/180.
tant = sin(theta)/cos(theta)
norm = sqrt (1. + tant*tant)
FS=" |:"
}
{
if ($1 == "#") {
# get column indices of the relevant fields (X,Y)
for (i = 1; i <= NF; i++) {
if ($i == "T")
iT = $(i-1);
else if ($i == "X")
iX = $(i-1);
else if ($i == "Y")
iY = $(i-1);
}
}
else {
T = $iT
if (T > 0. && T < 1.) {
# the cell contains an interface
x = $iX
y = $iY
print $1,(y - tant*x)/norm;
fflush (stdout);
}
}
}
The coordinates of the probe 'wires' themselves are generated using the following awk script. Note that since the x-axis of the numerical model is aligned with the slope (not the free surface) everything need to be rotated into this coordinate system (compared to the experimental coordinates).
BEGIN {
angle = 15.*3.14159265359/180.
for (y = -0.1; y <= 0.1; y+= 0.001) {
x = 0.544478770284516
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe1x"
x = 1.469
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0.35 > "probe2x"
x = 1.929
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe3x"
x = 1.929
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0.5 > "probe4x"
}
}
Note that the heights are not interpolated at the location of the wire itself i.e. the height measurement is not very accurate which explains the "noise" in the numerical results below (this is much larger than the real numerical noise due to e.g. VOF interface reconstruction).
Results
The Gerris parameter file below is used to generate the following results. This is a relatively straightforward file. The only trick is the use of a "reduced gravity" (or "piezometric pressure") to ensure accurate hydrostatic balance (the two SourceTension lines). More info on each object can be found in the syntax reference.
The simulation can be run like this:
% gerris3D -m enet.gfs
this will also generate a movie.mpg MPEG movie on-the-fly (provided you have ffmpeg installed on your system). If you want to visualise things interactively as it runs. Just add:
OutputSimulation { istep = 10 } stdout
and restart using
% gerris3D -m enet.gfs | gfsview3D movie.gfv
I ran the simulation for two resolutions, 8 and 9 (the 'LEVEL' macro at the top of the file). The simulation is setup to correspond to the shallower case of Enet and Grilli, 2007 i.e. d = 61 mm. A snapshot of the movie looks like:
which shows the wave about to break. This is already different from the experiment since according to Enet et al, 61 mm was the shallowest depth for which they didn't observe breaking. One can compare the wave gauge data. To do this I digitised the curves in Figure 13.a and 13.b. Plotting these with gnuplot e.g.
plot './probe1' every 10 w l, '../data/fig13a-a.dat', './finer/probe1' every 10 w l
gives
The finer resolution simulation shows smaller-amplitude oscillations in the "wake" of the landslide. This is probably due to a better resolution of the dissipative small-scale turbulent structures there (which are clearly seen in the experimental photos). Not quite enough resolution though to fully dissipate the initial wave as in the experimental timeseries (almost flat for t > 0.7 sec).
The most important discrepancy is the significantly larger initial trough in the model compared to the experiment (independently of spatial resolution). This is also consistent with the numerical results showing wave breaking whereas the experiment doesn't. Interestingly the initial "wave acceleration" is also lower in the numerics (the maximum trough is reached later) which is the opposite of what we would expect (larger acceleration causes larger troughs).
Aside from numerical errors (and trivial mistakes!), one hypothesis could be that imposing the motion of the solid (using the fitted experimental measurements) misses some subtle (or not so subtle) dynamic couplings between wave and solid motion. Given that the wave history is largely determined by the initial acceleration and that the acceleration is fitted over a longer period, it is not impossible that important "acceleration transients" are fitted away in this approach.
This probe shows a similar result. The model dissipation seems also to be too low compared to the experiment (in addition to the wave being to high). There is no explicit dissipation in the model at the moment, so it's actually reassuring (it's easy to add dissipation but difficult to remove it!).
Aside from checking that the simulation is setup properly (geometry, depth etc...) it would be interesting to see how the model compares for other initial submergence depths.
Other Depths
Gerris simulations were run for the other initial submersion depths. These are plotted against digitised results from Enet et al.
Note that the position of gauge 1 was variable as it was placed above the peak of the sliding block. The other gauges were at fixed positions.
A summary of the difference between the initial drawdown from the experiment and the model for the various depths and various probes are given below
Gauge 1
Depth (mm) | model time (s) | model \eta_0 (mm) | data time (s) | data \eta_0 (mm) | Ratio mod/dat | 61.00 | 0.52 | 17.2 | 0.48 | 13.0 | 1.32 | 80.00 | 0.50 | 12.8 | 0.47 | 9.2 | 1.39 | 100.00 | 0.50 | 9.6 | 0.48 | 7.8 | 1.23 | 120.00 | 0.51 | 7.7 | 0.49 | 5.1 | 1.51 | 140.00 | 0.53 | 6.0 | 0.49 | 4.4 | 1.36 | 149.00 | 0.52 | 5.8 | 0.47 | 4.2 | 1.38 | 189.00 | 0.52 | 4.0 | 0.49 | 3.1 | 1.29 |
Gauge 2
Depth (mm) | model time (s) | model \eta_0 (mm) | data time (s) | data \eta_0 (mm) | Ratio mod/dat | 61.00 | 1.59 | 16.5 | 1.59 | 12.2 | 1.36 | 120.00 | 1.34 | 11.7 | 1.38 | 8.6 | 1.35 | 189.00 | 1.07 | 6.0 | 1.12 | 4.4 | 1.36 |
Gauge 3
Depth (mm) | model time (s) | model \eta_0 (mm) | data time (s) | data \eta_0 (mm) | Ratio mod/dat | 61.00 | 1.91 | 20.0 | 1.91 | 15.1 | 1.32 | 189.00 | 1.45 | 11.7 | 1.48 | 9.2 | 1.27 |
Gauge 4
Depth (mm) | model time (s) | model \eta_0 (mm) | data time (s) | data \eta_0 (mm) | Ratio mod/dat | 61.00 | 1.96 | 10.6 | 1.92 | 8.5 | 1.24 | 120.00 | 1.74 | 8.7 | 1.74 | 6.2 | 1.42 | 189.00 | 1.49 | 5.7 | 1.47 | 4.5 | 1.26 |
Model time and data time are the times when the first trough is reached for the model and the data respectively. Model eta_0 and data eta_0 are the magnitude of that drawdown. The ratio is the ratio of the two values of eta_0 - model/data.
The model appears to overpredict the drawdown by around 30% (with a fair amount of variance)
There is some experimental variability in that the given values for the gauge 1 (x_0) agree with neither the calculated values of x_0 including a 4mm offset nor excluding it but generally lie somewhere between. The differences are small however.
Depth (mm) | 61 | 80 | 100 | 120 | 140 | 149 | 189 | x_0 (experiment) | 551 | 617 | 696 | 763 | 846 | 877 | 1017 | x_0 (no offset) | 545 | 615 | 690 | 765 | 839 | 873 | 1022 | x_0 (4mm offset) | 560 | 631 | 706 | 780 | 855 | 888 | 1038 |
There is also the issue that the slider is actually 4mm above the bottom of the tank. 9 levels is not fine enough to resolve this difference, 10 levels might just resolve it adequately. When I tried modelling it with 9 levels the solution blew up, I have not yet tried 10 levels. As an interim I also did a run where I increased the solid 4 mm higher and ran this to compare. These runs tended to have slightly larger draw down than the straight runs.
Gerris parameter file 'enet.gfs'
Define LEVEL 9
Define RHO MAX(0.01, (0.01*T + (1. - T)))
Define ANGLE (15.*M_PI/180.)
Define DEPTH 0.061
1 0 GfsSimulationMoving GfsBox GfsGEdge { x = 0.4 y = 0.5 } {
Refine 6
PhysicalParams { L = 3.7 }
VariableTracerVOF T
InitFraction T (y - tan(ANGLE)*x)
RefineSurface LEVEL (y - tan(ANGLE)*x)
# Solid (sphere (0., 0., 0., 0.5)) {
SolidMoving enet.gts {
rx = -90
ry = 90
# xp0(d)=(d*cos(theta)+T)/tan(theta)+sin(theta)*d
tx = 0.541714067835181
ty = 1e-6
flip = 1
} { level = LEVEL }
RefineSolid LEVEL
SurfaceBc U Dirichlet {
double ut = 1.70; /* m/s */
double a0 = 1.12; /* m/s^2 */
double t0 = ut/a0;
return ut*tanh(t/t0);
}
PhysicalParams { alpha = 1./RHO }
VariablePosition Y T y
VariablePosition X T x
SourceTension T -9.81*cos(ANGLE)*(0.01 - 1.) Y
SourceTension T 9.81*sin(ANGLE)*(0.01 - 1.) X
# SourceDiffusion U 1e-6
# SourceDiffusion V 1e-6
Time {
end = 2
# dtmax = 1e-2
}
AdaptVorticity { istep = 1 } {
cmax = 0.05
maxlevel = (T == 0. ? LEVEL : 0)
}
AdaptFunction { istep = 1 } {
cmax = 0
maxlevel = LEVEL
} (T > 0. && T < 1.)
OutputLocation { istep = 1 } { awk -f distance.awk > probe1 } probe1x { interpolate = 0 }
OutputLocation { istep = 1 } { awk -f distance.awk > probe2 } probe2x { interpolate = 0 }
OutputLocation { istep = 1 } { awk -f distance.awk > probe3 } probe3x { interpolate = 0 }
OutputLocation { istep = 1 } { awk -f distance.awk > probe4 } probe4x { interpolate = 0 }
OutputTime { istep = 1 } stderr
OutputBalance { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
OutputSimulation { step = 0.5 } sim-%g.gfs
GModule gfsview
OutputView { step = 2e-2 } { ppm2mpeg > movie.mpg } { width = 800 height = 600 } movie.gfv
}
GfsBox {
bottom = Boundary
# left = Boundary
# right = Boundary
}
GfsView parameter file 'movie.gfv'
Generated by GfsView
# GfsView 3D
View {
tx = -0.236536 ty = 0.137138
sx = 1 sy = 1 sz = 1
q0 = -0.275997 q1 = 0.228701 q2 = 0.327587 q3 = 0.87419
fov = 5.2916
r = 0.3 g = 0.4 b = 0.6
res = 1
lc = 0.001
reactivity = 0.1
}
VOF {
r = 1 g = 1 b = 1
shading = Constant
maxlevel = -1
font_size = 1
raster_font = 1
} {
n.x = 0 n.y = 0 n.z = 1
pos = 0
} (T > 0 && T < 1 ? (Y-tan(15.*M_PI/180.)*X)/sqrt(1.+tan(15.*M_PI/180.)*tan(15.*M_PI/180.)):0) {
amin = 0 min = -0.038413
amax = 0 max = 0.0524416
cmap = Jet
} T {
reversed = 1
use_scalar = 1
draw_edges = 0
interpolate = 0
}
Notes
There is a typo in Enet, 2003. Equations (1a) and (1b) are inconsistent. A consistent formulation (in gnuplot syntax), could be:
w = 0.7 b = 0.4 T = 0.08 r = 0.6 sech(x)=2./(exp(x) + exp(-x)) asech(x) = log(sqrt(1./x - 1.)*sqrt(1./x + 1.) + 1./x) Kw = 2./w*asech(1. - r) Kb = 2./b*asech(1. - r) max(a,b) = (a > b ? a : b) z(x,y) = max(0., T/r*(sech(Kw*x)*sech(Kb*y) - (1. - r))) splot [-0.4:0.4][-0.4:0.4]z(x,y) w l
Note the (1. - r) rather than (1 - r)/r in the article for Kw and Kb.
A different but consistent formula is given in Enet and Grilli, 2007.
