Enet et al benchmark

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 00:59, 23 November 2010
EmilyMLane (Talk | contribs)
(Other Depths)
← Previous diff
Current revision
EmilyMLane (Talk | contribs)
(Other Depths)
Line 187: Line 187:
Depth (mm) | model time (s) | model \eta_0 (mm) | data time (s) | data \eta_0 (mm) | Ratio mod/dat | 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 | 12.9 | 1.33 |+ 61.00 | 0.52 | 17.2 | 0.48 | 13.0 | 1.32 |
80.00 | 0.50 | 12.8 | 0.47 | 9.2 | 1.39 | 80.00 | 0.50 | 12.8 | 0.47 | 9.2 | 1.39 |
- 100.00 | 0.50 | 9.6 | 0.48 | 7.9 | 1.21 |+ 100.00 | 0.50 | 9.6 | 0.48 | 7.8 | 1.23 |
- 120.00 | 0.51 | 7.7 | 0.49 | 4.8 | 1.59 |+ 120.00 | 0.51 | 7.7 | 0.49 | 5.1 | 1.51 |
- 140.00 | 0.53 | 6.0 | 0.49 | 4.0 | 1.52 |+ 140.00 | 0.53 | 6.0 | 0.49 | 4.4 | 1.36 |
- 149.00 | 0.52 | 5.8 | 0.47 | 3.9 | 1.49 |+ 149.00 | 0.52 | 5.8 | 0.47 | 4.2 | 1.38 |
- 189.00 | 0.52 | 4.0 | 0.49 | 3.0 | 1.35 |+ 189.00 | 0.52 | 4.0 | 0.49 | 3.1 | 1.29 |
Gauge 2 Gauge 2
Line 219: Line 219:
The model appears to overpredict the drawdown by around 30% (with a fair amount of variance) 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 gauges (x_0) agree with neither the calculated values of x_0+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
-according to the setup shown in the paper nor the values of x_i as shown in the paper but generally lie somewhere between.+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 | Depth (mm) | 61 | 80 | 100 | 120 | 140 | 149 | 189 |
x_0 (experiment) | 551 | 617 | 696 | 763 | 846 | 877 | 1017 | x_0 (experiment) | 551 | 617 | 696 | 763 | 846 | 877 | 1017 |
x_0 (no offset) | 545 | 615 | 690 | 765 | 839 | 873 | 1022 | x_0 (no offset) | 545 | 615 | 690 | 765 | 839 | 873 | 1022 |
- x_i (no offset) | 542 | 615 | 692 | 770 | 847 | 882 | 1036 | 
x_0 (4mm offset) | 560 | 631 | 706 | 780 | 855 | 888 | 1038 | x_0 (4mm offset) | 560 | 631 | 706 | 780 | 855 | 888 | 1038 |
- x_i (4mm offset) | 557 | 630 | 707 | 785 | 862 | 897 | 1051 | 
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. 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.

Current revision

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:

t = 1 s, 9 levels
Enlarge
t = 1 s, 9 levels

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

probe 1
Enlarge
probe 1

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.

probe 2
Enlarge
probe 2

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.

probe 1
Enlarge
probe 1

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.

probe 2
Enlarge
probe 2
probe 3
Enlarge
probe 3
probe 4
Enlarge
probe 4

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.

Personal tools
communication