Liu et al benchmark

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 18:40, 27 March 2011
Delauxs (Talk | contribs)
(3D Results)
← Previous diff
Current revision
Delauxs (Talk | contribs)
(Numerical "wave gauge")
Line 71: Line 71:
== Numerical "wave gauge" == == Numerical "wave gauge" ==
[ [
-The setup of the numerical gauges is identical to that of [[Enet_et_al_benchmark]]. Please see [[Enet_et_al_benchmark]] for more details. Only the location of the gauges and the angle of the slope angle have to be modified. Consequently the 'distance.awk' which is used to calculate the surface elevation becomes:+The setup of the numerical gauges is slightly different to that of [[Enet_et_al_benchmark]]. The location of the gauges, angle and length have been modified. In an attempt to mimic the behavior of a capacitance gauge, the height is defined as the area of the gauge wire that is in contact with the sea-water.
<source lang=bash> <source lang=bash>
-BEGIN {+BEGIN {
- +
- tant = 0.4572/0.9144+
- norm = sqrt (1. + tant*tant)+
FS=" |:" FS=" |:"
 + sum=0
} }
{ {
 + if (NR > 2 && t0 != $1) {
 + print $1,0.2-sum;
 + fflush (stdout);
 + sum=0;
 + }
 +
if ($1 == "#") { if ($1 == "#") {
# get column indices of the relevant fields (X,Y) # get column indices of the relevant fields (X,Y)
Line 86: Line 90:
if ($i == "T") if ($i == "T")
iT = $(i-1); iT = $(i-1);
- else if ($i == "X") 
- iX = $(i-1); 
- else if ($i == "Y") 
- iY = $(i-1); 
} }
} }
else { else {
- T = $iT+ if ($1 == t0) {
- if (T > 0. && T < 1.) {+ T = $iT
- # the cell contains an interface+ sum = sum + 0.001*(T+T0)/2.;
- x = $iX+
- y = $iY+
- print $1,(y - tant*x)/norm;+
- fflush (stdout);+
} }
} }
 +
 + T0=$iT; X0=$2; Y0=$3; t0=$1
} }
</source> </source>
Line 110: Line 108:
BEGIN { BEGIN {
angle = atan2(0.4572,0.9144) angle = atan2(0.4572,0.9144)
- for (y = -0.1; y <= 0.1; y+= 0.001) {+ for (y = -0.2; y <= 0.2; y+= 0.001) {
x = 1.83 x = 1.83
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe1x" print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe1x"

Current revision

This benchmark for simulations of landslide-generated tsunamis using Gerris is based on these articles:

P. L.-F. Liu, T.R. Wu, F. Raichlen, C.E. Synolakis, J. C. Borrero - Runup and rundown generated by three-dimensional sliding masses

Journal of Fluid Mechanics 536:107-144, 2005
Bibtex

S. Abadie, D. Morichon, S. Grilli, S. Glockner - Numerical simulation of waves generated by landslides using a multiple-fluid Navier-Stokes model

Coastal Engineering 57:779-794, 2010
Bibtex

We focus here on one of the configuration that was retained as one of the benchmark test of the third international workshop on long-wave runup models (run 30). We look at the wave created by the free fall of a half-rectangle slide along a uniform slope.

The design of this page is inspired from that dedicated to the Enet_et_al_benchmark.

Contents

Generating the landslide body geometry

A GTS surface representation of the slide has to be generated to be fed to the moving boundary implementation GfsSolidMoving. This can be done using the following script:

# Generates the points of the triangulation
awk -f liu.awk > points.tmp
# Generates and rotates the triangulation
( wc points.tmp | awk '{print $1" 0 0"}'; cat points.tmp ) | delaunay -d | transform --rx=-90 > liu.gts

The awk script 'liu.awk' allows to generate the points of the triangulation. The notations a and b correspond to that of (Liu et al., 2005)

BEGIN {
a=0.4572
b=0.9144
alpha=atan2(a,b);
h=a*cos(alpha);
l1=b*cos(alpha);
l2=a*sin(alpha);
 
for (k=0; k <= 100; k++) {
for (i=-10; i <= 100; i++) { printf("%5.4f %5.4f %5.4f\n",l1*i/100, -0.5+k/100 , h*i/100) }
for (i=1; i <= 110; i++) { printf("%5.4f %5.4f %5.4f\n",l1+l2*i/100, -0.5+k/100, h*(1-i/100)) }
}
 
 
dx=1.0e-4
k=0.
j=-10.
for (i=-10; i <= 100; i++) { printf("%5.4f %5.4f %5.4f\n",l1*i/100, -0.5+k/100-dx,h*j/100) }
for (i=1; i <= 110; i++) { printf("%5.4f %5.4f %5.4f\n",l1+l2*i/100, -0.5+k/100-dx,h*j/100) }
 
k=100.
for (i=-10; i <= 100; i++) { printf("%5.4f %5.4f %5.4f\n",l1*i/100, -0.5+k/100+dx,h*j/100) }
for (i=1; i <= 110; i++) { printf("%5.4f %5.4f %5.4f\n",l1+l2*i/100, -0.5+k/100+dx,h*j/100) }
}

Numerical "wave gauge"

[ The setup of the numerical gauges is slightly different to that of Enet_et_al_benchmark. The location of the gauges, angle and length have been modified. In an attempt to mimic the behavior of a capacitance gauge, the height is defined as the area of the gauge wire that is in contact with the sea-water.

BEGIN {  
FS=" |:"
sum=0
}
{
if (NR > 2 && t0 != $1) {
print $1,0.2-sum;
fflush (stdout);
sum=0;
}
 
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 ($1 == t0) {
T = $iT
sum = sum + 0.001*(T+T0)/2.;
}
}
 
T0=$iT; X0=$2; Y0=$3; t0=$1
}

and the coordinates of the probe 'wires' can be generated using the awk script:

BEGIN {
angle = atan2(0.4572,0.9144)
for (y = -0.2; y <= 0.2; y+= 0.001) {
x = 1.83
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe1x"
x = 1.2446
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0.635 > "probe2x"
x = 0
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe3x"
x = 0
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0.305 > "probe4x"
x = 0
print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0.61 > "probe5x"
}
}

3D Results

So far, the full simulations I have managed to run give reasonable results. Those can be commented in the light of Stephane's previous comment on the Enet_et_al_benchmark. Wave gauge comparisons are possibly better for this test case, but we can observe the same trends in the results i.e. the over-prediction of the wave amplitude and the simulation is slightly late on the first trough. While this might not appear to be very satisfactory, this is as good as the results of (Abadie et al., 2010).

Wave generation is generally mainly governed by pressure forces and as a first approximation it was assumed that the fluids here are inviscid. In order to assess whether the validation on this benchmark test can be improved to more runs will be made: one with viscous fluids and a second one with a better resolution.

Evolution in time of the VOF interface for a simulation with adaptive refinement down to level 8
0.6 s
Enlarge
0.6 s
1.0 s
Enlarge
1.0 s
1.4 s
Enlarge
1.4 s
1.8 s
Enlarge
1.8 s
2.2 s
Enlarge
2.2 s
2.6 s
Enlarge
2.6 s


Comparison of the evolution of wave height with time with the measurements of (Liu et al., 2005) at two locations
Gauge 1
Enlarge
Gauge 1
Gauge 2
Enlarge
Gauge 2

Adding viscosity to the simulation seem to have essentially an influence on the post wave breaking part of the simulation and in improves the results. It is interesting to see that the "noise" associated with wave breaking almost vanishes except around time 1.6 s for the first gauge wich also corresponds to a secondary peak in (Liu et al., 2005)'s measurements.

Increasing the resolution changes little the first half of the simulation. The secondary peak is still captured by gauge 1 with the right timing and amplitude. The amplitude of the wave at time 2 is around 15% higher than for a refinement down to level 8 for gauge 1, but does not vary much for gauge 2. For times larger than 2 seconds the wave high is slightly different than for lower resolutions.

Comparison of the evolution of wave height with time with the measurements of (Liu et al., 2005) at two locations
Gauge 1
Enlarge
Gauge 1
Gauge 2
Enlarge
Gauge 2

2D/3D Comparisons

Preliminary simulations were run in two-dimension. Great agreement is not expected but this is a useful step to go through to assess the behavior of the code and tune the simulation. Even though the scale are slightly different between the plots of the 3D and 2D simulation presented hereafter, it is possible to see that 2D and 3D simulations remains qualitatively more or less similar during the first half of the simulation. Then wave-breaking occurs with accompanies with a dampening of the wave height in the 3D case, which is not the case in the 2D simulation.

Globally the wave height in 2D are much larger than in 3D or in the laboratory experiments which is not a surprise as momentum cannot travel in the transverse direction. The strong wave-breaking found in 2D makes the wave height at the different gauges hard to assess, but the arrival time of the first wave seem to be predicted well enough while the wave amplitude from the start are two to three time large than what is found in reality.


Comparison between 2D simulations (right) and the central cross-section of the 3D simulations for similar resolution. From top to bottom snapshots are for time = 0.6, 1.0, 1.4, 1.8, 2.2 and 2.6 s.
Enlarge

Gerris parameter file 'liu.gfs'

The time series of the position of the center of gravity of the solid during the experiment can be downloaded here.

The simulation can be run using: gerris3D -m liu.gfs

Define LEVEL 8
Define LMIN 4
Define RHO MAX(0.01, (0.01*T + (1. - T)))
Define ANGLE (atan2(0.4572,0.9144))
 
2 1 GfsSimulationMoving GfsBox GfsGEdge { y = 0.5 } {
 
# Sets the duration of the simulation to 3 time units
Time { end = 3 dtmax = 0.01}
 
# The CFL tolerance is set to 0.4
AdvectionParams { cfl = 0.4 }
 
Global { double interface (double x, double y) {
return ((y) - 0.4572/0.9144*x);
}
 
}
 
# Initial refinement
Refine LMIN
RefineSurface LEVEL (interface(x,y))
 
# Physical length of a box
PhysicalParams { L = 5. }
 
# Tank bottom
Solid ((y) - 0.4572/0.9144*(x- 1.0223302793129*2.44/0.4572))
 
# Fluid
VariableTracerVOF T
InitFraction T (interface (x,y))
 
# We lower the tolerance on the Poisson solver
ApproxProjectionParams { tolerance = 1e-5 }
ProjectionParams { tolerance = 1e-5 }
 
# Definition of functions that read the file containing the time evolution of the position
# of the solid object and derive the corresponding boundary conditions.
Global {
double position (double t) {
FILE * fin = fopen ("Benchmark_4_run30.txt","r");
double t0=0, p0=0, t1=0, p1=0;
 
g_assert(fscanf(fin,"%lf %lf\n",&t1,&p1));
 
while (t1 <= t) {
t0=t1; p0=p1;
g_assert(fscanf(fin,"%lf %lf\n",&t1,&p1));
 
}
 
fclose (fin);
 
if (t1 == t) {
return p1/100;
}
else {
return (p0*(t1-t)+p1*(t-t0))/(100*(t1-t0));
}
}
 
double instant_vel (double t) {
FILE * fin = fopen ("Benchmark_4_run30.txt","r");
double t0=0, p0=0, t1=0, p1=0;
 
g_assert(fscanf(fin,"%lf %lf\n",&t1,&p1));
 
while (t1 <= t) {
t0=t1; p0=p1;
g_assert(fscanf(fin,"%lf %lf\n",&t1,&p1));
 
}
 
fclose (fin);
 
return (p1-p0)/(t1-t0);
}
 
double vel_x (double t, double dt) {
if (dt == 0.)
return (instant_vel (t));
else
return ((position(t+dt)-position(t))/dt);
}
}
 
# Setup of the moving slide
SolidMoving liu.gts {
tx = 0.14142135623731
ty = 1e-6
sz = 0.6525
flip = 1
} { level = LEVEL }
 
#
SurfaceBc U Dirichlet (vel_x (t, dt))
 
# Setup of mesh refinement on the solid boundaries
RefineSolid ( x < 1.5 ? LEVEL : LMIN)
 
# Density
PhysicalParams { alpha = 1./RHO }
 
 
# Use the reduced gravity approach
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
 
 
AdaptFunction { istart = 1 istep = 1 } {
cmax = 0
maxlevel = (x < 4. ? LEVEL : LEVEL - 2)
minlevel = (T < 1. ? LMIN:0)
} (T > 0 && T < 1)
 
 
AdaptVorticity { istep = 1 } {
cmax = 1.0
maxlevel = LEVEL
minlevel = LMIN
}
 
# Wave gauges
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 }
OutputLocation { istep = 1 } { awk -f distance.awk > probe4 } probe5x { interpolate = 0 }
 
# Output control statistics on the simulation
OutputTime { istep = 1 } stderr
OutputBalance { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
 
# Output simulation files every 0.1 second
OutputSimulation { step = 0.1 end = 3} sim-%g.gfs
 
# Viscosity using a fully implicit scheme
SourceViscosity {} (1./(T/1.78e-5+(1-T)/1.015e-3)) { beta = 1}
 
# Generates a movie
GModule gfsview
OutputView { step = 2e-2 } { ppm2mpeg > movie.mpg } { width = 800 height = 600 } view3D.gfv
 
}
GfsBox {
left = BoundaryOutflow
top = BoundaryOutflow
}
GfsBox {
top = BoundaryOutflow
}
1 2 right

GfsView3D file 'view3D.gfv'

Generated by GfsView


# GfsView 3D
View {
  tx = -0.303493 ty = 0.0588803
  sx = 1 sy = 1 sz = 1
  q0 = -0.175687 q1 = 0.0609815 q2 = 0.294233 q3 = 0.937466
  fov = 4.37797
  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
} P {
  amin = 1
  amax = 1
  cmap = Jet
} T {
  reversed = 1
  use_scalar = 0
  draw_edges = 0
  interpolate = 0
}
Solid {
  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
} P {
  amin = 1
  amax = 1
  cmap = Jet
} {
  reversed = 0
  use_scalar = 0
}

Personal tools
communication