# GfsRiver

(Difference between revisions)
 Revision as of 23:43, 14 May 2009Popinet (Talk | contribs)← Previous diff Revision as of 04:02, 28 August 2009Popinet (Talk | contribs) (Added description)Next diff → Line 1: Line 1: + GfsRiver solves the [[w:Shallow_water_equations|Saint-Venant equations]] (also known as the non-linear shallow-water equations) in conservative form. Wetting and drying, hydraulic jumps etc... are handled appropriately by the solver. + + The default variables are: + + ;U: horizontal component of the flux (dimension L^2/T) + ;V: vertical component of the flux + ;P: fluid depth + ;Zb: bed elevation + ;H: Zb + P i.e. the elevation of the free surface + + The value of the acceleration of gravity is set using [[GfsPhysicalParams]] (default to 1). +

## Revision as of 04:02, 28 August 2009

GfsRiver solves the Saint-Venant equations (also known as the non-linear shallow-water equations) in conservative form. Wetting and drying, hydraulic jumps etc... are handled appropriately by the solver.

The default variables are:

`U`
horizontal component of the flux (dimension L^2/T)
`V`
vertical component of the flux
`P`
fluid depth
`Zb`
bed elevation
`H`
`Zb + P` i.e. the elevation of the free surface

The value of the acceleration of gravity is set using GfsPhysicalParams (default to 1).

### Examples

• Dam break on complex topography
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
PhysicalParams { L = 8. }

# Define a 1D domain using cell masking
RefineSurface 9 (y - 8.*(0.5 - 1./512.))
InitMask {} (y < 8.*(0.5 - 1./512.))

# Set the topography Zb and the initial water surface elevation P
Init {} {
Zb = x*x/8.+cos(M_PI*x)/2.
P = {
double p = x > 0. ? 0.35 : x < -2. ? 1.9 : -1.;
return MAX (0., p - Zb);
}
}
PhysicalParams { g = 1. }

# Use a first-order scheme rather than the default second-order
# minmod limiter. This is just to add some numerical damping.
}

Time { end = 40 }
OutputProgress { istep = 10 } stderr
OutputScalarSum { istep = 10 } ke { v = (P > 0. ? U*U/P : 0.) }
OutputScalarSum { istep = 10 } vol { v = P }
OutputScalarNorm { istep = 10 } u { v = (P > 0. ? U/P : 0.) }

# use gfsplot/gnuplot for online visualisation and generation of figures
OutputSimulation { step = 0.3 } {
gfsplot "
set title sprintf('t = %4.1f', (t))
plot [-4.:4.]'-' u (x):(Zb):(H) w filledcu lc 3, \
'-' u (x):(Zb) w l lw 4 lc 1 lt 1
set term pngcairo size 800,600
set output sprintf('sim-%04.1f.png', (t))
replot
set term wxt noraise
"
} { format = text }

# Combine all the gif images into a gif animation using gifsicle
EventScript { start = end } {
echo "\rcreating animation...            " > /dev/stderr
sleep 10 # give a chance to gnuplot to catch up
for f in sim-*.png; do
convert \$f -trim +repage -bordercolor white \
-border 10 -resize 640x282! `basename \$f .png`.gif && rm -f \$f
done
gifsicle --colors 256 --optimize --delay 12 --loopcount=0 sim-*.gif > dam.gif && \
rm -f sim-*.gif
}
}
```

• Small amplitude solitary wave interacting with a parabolic hump
• ```2 1 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
Refine 8
Init {} {
# Parabolic hump
Zb = 0.8*exp(-5.*(x - 0.9)*(x - 0.9) - 50.*(y - 0.5)*(y - 0.5))
# Initial free surface and perturbation
P = (0.05 < x && x < 0.15 ? 1.01 : 1) - Zb
}
PhysicalParams { g = 1 }
AdvectionParams { cfl = 0.5 }
cmax = 1e-4
cfactor = 2
maxlevel = 8
minlevel = 6
} (P + Zb)
Time { end = 1.8 }
OutputTime { istep = 10 } stderr
OutputSimulation { istep = 10 } stdout
EventScript { istep = 10 } { echo "Save stdout { width = 640 height = 480 }" }
OutputSimulation { start = 0.6 step = 0.3 } sim-%g.gfs
EventScript { start = end } {
for i in 0.6 0.9 1.2 1.5 1.8; do
echo "Save stdout { format = EPS line_width = 0.2 }" | \
gfsview-batch2D sim-\$i.gfs isolines.gfv > iso-\$i.eps
echo "Save stdout { format = EPS line_width = 0.2 }" | \
gfsview-batch2D sim-\$i.gfs cells.gfv > cells-\$i.eps
done
echo "Save stdout { width = 1280 height = 960 }" | \
gfsview-batch2D sim-0.9.gfs hump.gfv | convert ppm:- hump.eps
}
}
```

• Shock reflection by a circular cylinder
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
Time { end = 0.3 }

PhysicalParams { L = 5 g = 9.81 }

# We use a sphere knowing that in 2D the resulting object will be
# a cross-section of the sphere at z = 0 i.e. a cylinder of radius
# 0.5
Solid (sphere(-0.5,0.,0.,0.5))

Init {} {
# Initial shock
P = (x <= -1. ? 3.505271526 : 1.)
U = (x <= -1. ? 22.049341608 : 0.)
}

Refine 9
cmax = 0.1
cfactor = 2
maxlevel = 9
} P

OutputTime { istep = 10 } stderr
OutputSimulation { istep = 10 } stdout

GModule gfsview
OutputView { step = 0.0025 } {
ppm2mpeg -b 3600K > depth.mpg
} { width = 400 height = 400 } depth.gfv
OutputView { start = end } {
convert ppm:- depth.eps
} { width = 400 height = 400 } depth.gfv
OutputView { start = end } {
convert ppm:- mesh.eps
} { width = 400 height = 400 } mesh.gfv
}
```

• Tsunami runup onto a complex three-dimensional beach
• ```2 1 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
# the domain is 3.402 m X 6.804 m
# units for time are seconds
PhysicalParams { L = 3.402 g = 9.81 }
Refine 6

# maintain the Zb1 variable using the GTS surface
VariableFunction Zb1 bathy.gts

# the initial water level is at z = 0, so the depth P is...
Init {} {
P = MAX (0., -Zb1)
}

# use a Sweby limiter rather than the default minmod which is too
# dissipative

# adapt down to 9 levels based on the slope of the (wet)
# free-surface and with a tolerance of 1 mm
cmax = 1e-3
cfactor = 2
maxlevel = 9
minlevel = 6
} (P < DRY ? 0. : P + Zb)

# at each timestep
Init { istep = 1 } {
# Add a "shelf" to simulate the wall on the right-hand-side boundary
Zb = (x > 5.448 ? 0.13535 : Zb1)
# read in the experimental timeseries in variable 'input'
input = input.cgd
# implicit quadratic bottom friction with coefficient 1e-3
U = (P > DRY ? U/(1. + dt*1e-3*Velocity/P) : 0.)
V = (P > DRY ? V/(1. + dt*1e-3*Velocity/P) : 0.)
P = (P > DRY ? P : 0.)
}

Time { end = 22.5 }
OutputTime { istep = 10 } stderr
OutputTiming { start = end } stderr
OutputSimulation { step = 1 } stdout
OutputSimulation { step = 1 } sim-%g.gfs
EventScript { step = 1 } { gzip -f sim-*.gfs }

# output data at probe locations
OutputLocation { istep = 1 } input 1e-3 1.7 0
OutputLocation { istep = 1 } p5 4.521 1.196 0
OutputLocation { istep = 1 } p7 4.521 1.696 0
OutputLocation { istep = 1 } p9 4.521 2.196 0

# kinetic energy
OutputScalarSum { istep = 1 } ke { v = (P > 0. ? Velocity2*P : 0.) }
# free-surface elevation
OutputScalarStats { istep = 1 } p { v = (Zb > 0. ? P : P + Zb) }
OutputScalarNorm { istep = 1 } u { v = Velocity }
OutputTime { istep = 1 } balance
OutputBalance { istep = 1 } balance

# generate movies
GModule gfsview
OutputView { start = 9 step = 0.0416 } { ppm2mpeg -s 640x480 > monai.mpg } {
width = 1280 height = 960
} 3D.gfv
OutputView { start = 14.63 end = 19.5 step = 0.033333333 } {
} {
width = 800 height = 1200
} {
dry = DRY
}
```

• The 2004 Indian Ocean tsunami
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
PhysicalParams {
# length of the domain (m)
L = LENGTH
# gravity is 9.81 m/s^2
g = 9.81
# from now on, units have been chosen to be metres and seconds
}
# use the longitude/latitude metric. x and y coordinates are now
# longitude and latitude in degrees
# shift the origin to where we want it
MapTransform { tx = LONGITUDE ty = LATITUDE }
# 10 hours
Time { end = 36000 }
# use the terrain module for topography
GModule terrain
# use the okada module for initial deformations
Refine 5
# maintain Zb as the terrain elevation defined by the samples in
# ETOPO1
# the 'etopo1' database needs to have been created beforehand and
# be accessible within GFS_TERRAIN_PATH, see the 'terrain' module
# documentation above
VariableTerrain Zb {
basename = etopo1
} {
# preserve lake-at-rest balance (but not volume) when
# reconstructing bathymetry
reconstruct = 1
}

# the default minmod limiter is too dissipative for tsunamis

# deformation from Grilli et al, 2007, Table 1
# 5 segments triggered at different times

# Segment 1
Init {} { D = 0 }
x = 94.57 y = 3.83
depth = 11.4857e3
strike = 323 dip = 12 rake = 90
length = 220e3 width = 130e3
U = 18
}
# Initial water level is at z = D
Init { start = 0 } { P = MAX (0., D - Zb) }

# Segment 2
EventList { start = 212 step = 6 end = 272 } {
Init {} { D = 0 }
x = 93.90 y = 5.22
depth = 11.4857e3
strike = 348 dip = 12 rake = 90
length = 150e3 width = 130e3
U = 23
}
}
# make sure the deformation is well resolved
AdaptGradient { start = 212 istep = 1 end = 272 } {
cmax = 0.05 cfactor = 2
minlevel = 5 maxlevel = LEVEL
} D
# add it to the current elevation field (only if wet)
Init { start = 272 } { P = (P < DRY ? P : MAX (0., P + D)) }

# Segment 3
EventList { start = 528 step = 6 end = 588 } {
Init {} { D = 0 }
x = 93.21 y = 7.41
depth = 12.525e3
strike = 338 dip = 12 rake = 90
length = 390e3 width = 120e3
U = 12
}
}
# make sure the deformation is well resolved
AdaptGradient { start = 528 istep = 1 end = 588 } {
cmax = 0.05 cfactor = 2
minlevel = 5 maxlevel = LEVEL
} D
# add it to the current elevation field (only if wet)
Init { start = 588 } { P = (P < DRY ? P : MAX (0., P + D)) }

# Segment 4
EventList { start = 853 step = 6 end = 913 } {
Init {} { D = 0 }
x = 92.60 y = 9.70
depth = 15.12419e3
strike = 356 dip = 12 rake = 90
length = 150e3 width = 95e3
U = 12
}
}
# make sure the deformation is well resolved
AdaptGradient { start = 853 istep = 1 end = 913 } {
cmax = 0.05 cfactor = 2
minlevel = 5 maxlevel = LEVEL
} D
# add it to the current elevation field (only if wet)
Init { start = 913 } { P = (P < DRY ? P : MAX (0., P + D)) }

# Segment 5
EventList { start = 1213 step = 6 end = 1273 } {
Init {} { D = 0 }
x = 92.87 y = 11.70
depth = 15.12419e3
strike = 10 dip = 12 rake = 90
length = 350e3 width = 95e3
U = 12
}
}
# make sure the deformation is well resolved
AdaptGradient { start = 1213 istep = 1 end = 1273 } {
cmax = 0.05 cfactor = 2
minlevel = 5 maxlevel = LEVEL
} D
# add it to the current elevation field (only if wet)
Init { start = 1273 } { P = (P < DRY ? P : MAX (0., P + D)) }

# we don't want the arrival time to be interpolated from coarse
# to fine, so we make it a "boolean" variable
VariableBoolean Atime
Init {} {
P = MAX (0., D - Zb)
Pmax = 0.
Atime = -1.
}

Init { istep = 1 } {
# elevation of the wet surface
Hwet = (P > DRY ? H : NODATA)
# maximum amplitude
Pmax = (P > DRY && H > Pmax ? H : Pmax)
# arrival time for an amplitude of +- 5 cm
Atime = (Atime >= 0. ? Atime : P > DRY && fabs (H) > 5e-2 ? t : -1.)
# implicit scheme for quadratic bottom friction with coefficient 1e-3
U = P > DRY ? U/(1. + dt*Velocity*1e-3/P) : 0
V = P > DRY ? V/(1. + dt*Velocity*1e-3/P) : 0
}

# with a tolerance of 5 cm
cmax = 0.05 cfactor = 2
minlevel = 5 maxlevel = LEVEL
} (P < DRY ? 0. : P + Zb)

# also keep enough resolution to resolve the maximum elevation
# field with a 5 cm discretisation error max
AdaptError { istep = 1 } {
cmax = 0.05
minlevel = 5 maxlevel = LEVEL
} Pmax

Global {
#define RESOLUTION (LENGTH/pow (2, MAXLEVEL))
#define close_enough(x,y) (distance(x,y,0.) < RESOLUTION)
}
# include a list of locations for which to output timeseries
Include output.gfs

# Jason profiles
OutputLocation { start = 6900 } jason.txt jason.xy

OutputTime { istep = 1 } stderr
OutputBalance { istep = 1 } stderr
OutputTiming { istep = 100 } stderr

# save a snapshot every 100 iterations
EventScript { istep = 100 } { rm -f snapshot-*.gfs }
OutputSimulation { istep = 100 } snapshot-%ld.gfs

# output solution on standard output every 20 timesteps
# for on-the-fly visualisation with GfsView
OutputSimulation { istep = 20 } stdout

# ouput solution every 900 seconds (15 minutes)
OutputSimulation { step = 900 } sim-%g.gfs
EventScript { step = 900} { gzip -f sim-*.gfs }

OutputScalarStats { istep = 1 } p { v = P }
OutputScalarNorm { istep = 1 } u { v = Velocity }
OutputScalarNorm { istep = 1 } U { v = U }
OutputScalarNorm { istep = 1 } V { v = V }
OutputScalarNorm { istep = 1 } hwet { v = Hwet }

# animation of the wave elevation
GModule gfsview
OutputView { step = 60 } { ppm2mpeg > h.mpg } {
width = 800 height = 700
} h.gfv
OutputView { start = 7200 } { convert ppm:- eps2:h.eps } {
width = 800 height = 700
} h.gfv
# animation of the level of refinement
OutputView { step = 60 } { ppm2mpeg > level.mpg } {
width = 800 height = 700
} level.gfv
OutputView { start = 7200 } { convert ppm:- eps2:level.eps } {
width = 800 height = 700
} level.gfv
# animation of the velocity
OutputPPM { step = 60 } { ppm2mpeg > velocity.mpg } {
maxlevel = 10
v = Velocity
min = 0 max = 0.25
}
# animation of the topography
OutputPPM { step = 60 } { ppm2mpeg > zb.mpg } {
maxlevel = 10
v = Zb
}
# graphics
OutputView { start = end } { convert ppm:- eps2:hmax.eps } {
width = 1600 height = 1600
} hmax.gfv
OutputView { start = end } { convert ppm:- eps2:hmax-detail.eps } {
width = 1600 height = 1600
} hmax-detail.gfv
EventScript { start = end } {
# timeseries at specific locations
gnuplot <<EOF
set term postscript eps color lw 2 size 5,2 18

set size ratio 0.4

set output 'hani.eps'
plot [3:5]'hanires.txt' u 1:(\\$2/100.) w lp t 'observed' ps 2, \
'hani.txt' u (\\$3/3600.):7 w l t 'modelled'

set output 'male.eps'
plot [3:5]'maleres.txt' u 1:(\\$2/100.) w lp t 'observed' ps 2, \
'male.txt' u (\\$3/3600.):7 w l t 'modelled'

set output 'gana.eps'
plot [3:5]'ganares.txt' u 1:(\\$2/100.) w lp t 'observed' ps 2, \
'gana.txt' u (\\$3/3600.):7 w l t 'modelled'

set output 'dieg.eps'
plot [3:5]'diegres.txt' u 1:(\\$2/100.) w lp t 'observed' ps 2, \
'dieg.txt' u (\\$3/3600.):7 w l t 'modelled'

set output 'colo.eps'
plot [2.5:4.5]'colores.txt' u 1:(\\$2/100.) w lp t 'observed' ps 2, \
'colo.txt' u (\\$3/3600.):7 w l t 'modelled'

set output 'jason.eps'
set xlabel 'Latitude (deg)'
set ylabel 'Elevation (m)'
plot [-6:20][-1:1]'jasonres.txt' u 2:4 w l t 'observed', \
'jason.txt' u 3:(\\$3 > 19.2 ? 0 : \\$8) w l t 'modelled'
EOF
}
} {
dry = DRY
}
```

• Poiseuille flow with multilayer Saint-Venant
• ```1 1 GfsRiver GfsBox GfsGEdge {} {
Layers NL
PhysicalParams { L = 1. }
Init {} { P = 0.5 }
Source U A*P/NL
EventStop { istep = 1 } U 1e-8 DU
#    OutputTime { istep = 1 } stderr
Time { iend = 20 }
OutputSimulation { start = end } {
awk '{
nl = NL;
if (\$1 == "#") {
for (i = 2; i <= NF; i++) {
split(\$i,a,":")
if (a[2] == "U0")
start = a[1];
}
}
else {
dz = \$4/nl;
emax = 0.;
for (i = 0; i < nl; i++) {
z = dz*(0.5+i)
u = A/(2.*NU)*(1./4 - (0.5 - z)*(0.5 - z))
eu = u - \$(start+i)/dz
if (eu < 0.)
eu = -eu;
if (eu > emax)
emax = eu;
}
print nl,emax
}
}'
} { format = text }
} {
# vertical viscosity
nu = NU
}
```

• Multi-layer Saint-Venant solver
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
Layers NL
PhysicalParams { L = RATIO g = 100./RE }
Refine LEVEL
InitMask {} (y < RATIO*(0.5 - 1./pow(2,LEVEL)))

Init {} { P = 1 }
Time { end = 1000 }
EventStop { step = 1 } U0 1e-9
#    OutputTime { step = 1 } stderr
#    OutputSimulation { start = end } end-RATIO-RE.txt { format = text }
#    OutputSimulation { start = end } end-RATIO-RE.gfs
OutputSimulation { start = end } {
awk '{
nl = NL;
if (\$1 == "#") {
for (i = 2; i <= NF; i++) {
split(\$i,a,":")
if (a[2] == "U0")
start = a[1];
}
}
else if (\$1 == RATIO/2**(LEVEL + 1)) {
dz = \$4/nl;
for (i = 0; i < nl; i++) {
z = dz*(0.5+i)
print \$1,\$2,z,\$(start+i)/dz
}
}
}' > uprof-RATIO-RE-NL
} { format = text }
} {
# vertical viscosity
nu = 1./RE
# Neumann condition at the surface
dut = 1.
}
```

• Wind-driven stratified lake
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
Layers NL
Refine LEVEL
InitMask {} (y < RATIO*(0.5 - 1./pow(2,LEVEL)))
VariableTracer DRHO
PhysicalParams {
L = RATIO g = 1./(ALPHA*RE)
alpha = 1./(1. + DRHO)
}
Init {} {
P = 1
DRHO = (z > 0.75 ? 0 : ALPHA/ALPHAT)*P/NL
}
Time { end = 100 }
OutputSimulation { start = end } end-NL.txt { format = text }
OutputSimulation { start = end } {
awk '{
nl = NL;
if (\$1 == "#") {
for (i = 2; i <= NF; i++) {
split(\$i,a,":")
if (a[2] == "U0")
start = a[1];
}
}
else if (\$1 == RATIO/2**(LEVEL + 1)) {
dz = \$4/nl;
for (i = 0; i < nl; i++) {
z = dz*(0.5+i)
print \$1,\$2,z,\$(start+i)/dz
}
}
}' > uprof-NL
} { format = text }
} {
# vertical viscosity
nu = 1./RE
# Neumann condition at the surface
dut = 1.
}
```

• ```1 0 GfsRiver GfsBox GfsGEdge {} {
Time { iend = 1580 dtmax = 1000 }
Refine 6
Global {
#define L0 1000e3
#define H0 1000
#define G 0.01
#define R0 100e3
#define ETA0 599.5
#define F0 1.0285e-4
}
PhysicalParams { L = L0 g = G }
Init {} {
# e-folding radius = 100 km
# umax = 0.5 m/s = sqrt(200)*exp(-1/2)
P = (H0 + ETA0*exp (-(x*x + y*y)/(R0*R0)))
U = 2.*G*ETA0*y/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))*P
V = - 2.*G*ETA0*x/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))*P
}
SourceCoriolis F0

OutputErrorNorm { istep = 1 } {
awk '{print \$3/86400. " " \$9; fflush (stdout); }' > e
} { v = P } {
s = (H0 + ETA0*exp (-(x*x + y*y)/(R0*R0)))
v = E
}
GModule gfsview
OutputView {  istart = 100 iend = 500 istep = 100 } error-%ld.eps { format = EPS } geo.gfv
OutputView {  istart = 1500 } error-%ld.eps { format = EPS } geo.gfv
EventScript { start = end } {
cat <<EOF | gnuplot
set term postscript eps lw 3 color solid 20
set output 'geo_error.eps'
set xlabel 'Time (days)'
set ylabel 'Maximum error on surface height (m)'
plot 'e.ref' t 'ref' w l, 'e' t '' w l
EOF
}
}
```

• Oscillations in a parabolic container
• ```1 0 GfsRiver GfsBox GfsGEdge {} {

# Analytical solution, see Sampson, Easton, Singh, 2006
Global {
static gdouble Psi (double x, double t) {
double p = sqrt (8.*G*h0)/a;
double s = sqrt (p*p - tau*tau)/2.;
return a*a*B*B*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) +
(tau*tau/4. - s*s)*cos (2.*s*t)) - B*B*exp(-tau*t)/(4.*G) -
exp (-tau*t/2.)/G*(B*s*cos (s*t) + tau*B/2.*sin (s*t))*x;
}
}

PhysicalParams { L = 10000 }
RefineSurface LEVEL (y - 10000.*(0.5 - 1./pow(2,LEVEL)))
InitMask {} (y < 10000.*(0.5 - 1./pow(2,LEVEL)))
Init {} {
Zb = h0*(x/a)*(x/a)
P = MAX (0., h0 + Psi (x, 0.) - Zb)
}
Init { istep = 1 } {
Pt = MAX (0., h0 + Psi (x, t) - Zb)
}
# Better convergence rates are obtained at lower CFL
# AdvectionParams { cfl = 0.1 }
PhysicalParams { g = G }
SourceCoriolis 0 tau
Time { end = 6000 }
OutputSimulation { start = 1500 } sim-LEVEL-%g.txt { format = text }
OutputSimulation { start = end } end-LEVEL.txt { format = text }
OutputScalarSum { istep = 10 } ke-LEVEL { v = (P > 0. ? U*U/P : 0.) }
OutputScalarSum { step = 50 } vol-LEVEL { v = P }
OutputScalarSum { step = 50 } U-LEVEL { v = U }
OutputErrorNorm { istep = 1 } error-LEVEL { v = P } {
s = Pt
v = DP
}
} {
# this is necessary to obtain good convergence rates at high
# resolutions
dry = 1e-4
}
```

• Parabolic container with embedded solid
• ```1 0 GfsRiver GfsBox GfsGEdge {} {

# Analytical solution, see Sampson, Easton, Singh, 2006
Global {
static gdouble Psi (double x, double y, double t) {
double p = sqrt (8.*G*h0)/a;
double s = sqrt (p*p - tau*tau)/2.;
return a*a*B*B*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) +
(tau*tau/4. - s*s)*cos (2.*s*t)) - B*B*exp(-tau*t)/(4.*G) -
exp (-tau*t/2.)/G*(B*s*cos (s*t) + tau*B/2.*sin (s*t))*R(x,y);
}
}

PhysicalParams { L = 10000 }
Refine LEVEL
Solid ({
double line1 = SLOPE*x - y + 28000./pow(2.,LEVEL);
double line2 = - SLOPE*x + y + 28000./pow(2.,LEVEL);
return union (line1, line2);
})
Init {} {
Zb = h0*pow(R(x,y), 2.)/(a*a)
P = MAX (0., h0 + Psi (x, y, 0.) - Zb)
}
Init { istep = 1 } {
Pt = MAX (0., h0 + Psi (x, y, t) - Zb)
}
PhysicalParams { g = G }
SourceCoriolis 0 tau
Time { end = 6000 }
OutputSimulation { start = 1500 } {
awk '
function atan(x)
{
return atan2(x,1.);
}
function pow(x,y)
{
return x**y;
}
{
if (\$1 == "#")
print \$0;
else {
printf ("%g", R(\$1,\$2));
for (i = 2; i <= NF; i++)
printf (" %s", \$i);
printf ("\n");
}
}' > sim-LEVEL-1500.txt
} { format = text }
OutputSimulation { start = end } end-LEVEL.txt { format = text }
OutputScalarSum { istep = 10 } ke-LEVEL { v = (P > 0. ? U*U/P : 0.) }
OutputScalarSum { step = 50 } vol-LEVEL { v = P }
OutputScalarSum { step = 50 } U-LEVEL { v = U*cos (ANGLE) + V*sin(ANGLE) }
OutputErrorNorm { istep = 1 } error-LEVEL { v = P } {
s = Pt
v = DP
}
} {
# this is necessary to obtain good convergence rates at high
# resolutions
dry = 1e-4
}
```

• Transcritical flow over a bump
• ```1 0 GfsRiver GfsBox GfsGEdge { x = 0.5 } {
PhysicalParams { L = 25 g = 9.81 }
Refine LEVEL
InitMask {} (y < 25.*(0.5 - 1./pow(2,LEVEL)))
Init {} {
Zb = MAX(0., 0.2 - 0.05*(x - 10.)*(x - 10.))
P = 0.33 - Zb
}
# Sweby seems to work better than minmod for this test case
}
EventStop { step = 1 } H 1e-5
Time { end = 1000 }
#    OutputTime { istep = 1 } stderr
OutputSimulation { start = end } {
gfsplot "
set term postscript eps color solid lw 2 16
set output 'h-LEVEL.eps'
set xlabel 'x'
set ylabel 'z'
plot [5:15]'swashes' u 1:4 t 'topography' w l, \
'swashes' u 1:6 t 'free surface (analytical)' w l, \
'-' u (x):(H) w p t 'free surface (numerical)'
"
} { format = text }
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' } { v = H } {
s = href.cgd
v = EH
}
OutputSimulation { start = end } end-LEVEL.txt { format = text }
}
```

• Transcritical flow with multiple layers
• ```1 0 GfsRiver GfsBox GfsGEdge { x = 0.5 } {
Global {
#define LENGTH 21.
#define WIDTH 5.75
#define HEIGHT 0.2
#define Q 1.
#define HE 0.6
#define G 9.81
}
Layers NL
PhysicalParams { L = LENGTH g = G }
Refine LEVEL
InitMask {} (y < LENGTH*(0.5 - 1./pow(2,LEVEL)))
Init {} {
Zb = MAX(0., HEIGHT*(1. - 1./(WIDTH/2.)/(WIDTH/2.)*(x - 10.)*(x - 10.)))
P = 0.6 - Zb
}
# Sweby seems to oscillate in this case
}
# Stop when stationary solution is reached
EventStop { step = 1 } H 1e-3
Time { end = 1000 }
# uncomment this for on-the-fly animation
# OutputSimulation { step = 0.1 } {
# 	gfsplot "
#       set term wxt noraise
#       set xlabel 'x'
#       set ylabel 'z'
#       set xrange [0:21]
#       set xtics 0,2,20
#       set key bottom right
#       plot '-' u (x):(Zb) t 'topography' w l, \
#            '-' u (x):(H) w p t 'free surface'
#     "
# } { format = text }
OutputLocation { start = end } { awk -v nl=NL -f uprof.awk > uprof-10-LEVEL-NL } 10. 10.49 0
OutputLocation { start = end } { awk -v nl=NL -f uprof.awk > uprof-15-LEVEL-NL } 15. 10.49 0
OutputLocation { start = end } { awk -v nl=NL -f uprof.awk > uprof-20-LEVEL-NL } 20. 10.49 0
OutputSimulation { start = end } {
awk '{ if (\$1 != "#") print \$1,\$8; }' > prof-LEVEL-NL
} { format = text }
OutputSimulation { start = end } end-LEVEL-NL { format = text }
} {
# vertical viscosity
nu = 0.01
# Gauckler-Manning-Strickler bottom friction, G = 9.81 m/s^2,
# Strickler coefficient S = 25 m^(1/3)/s.
k = (G/(25.*25.*pow(P,1./3.)))*Velocity
}
```

• Tsunami runup onto a plane beach
• ```1 0 GfsRiver GfsBox GfsGEdge { x = 0.3333333 } {
PhysicalParams { L = 60000 }

# Set a solid boundary close to the top boundary to limit the
# domain width to less than one cell (i.e. a 1D domain)
RefineSolid (MINLEVEL + (LEVEL - MINLEVEL)*(1. - x/50000.))
Solid (y - 60000.*(0.5 - 0.99/pow(2,LEVEL)))

# Set the topography Zb and the initial water surface elevation P
Init {} {
Zb = -x/10.
P = init.cgd
P = MAX (0., P - Zb)
}
PhysicalParams { g = 9.81 }

}

# Force the flow to stay 1D
Init { istep = 10 } { V = 0 }

Time { end = 220 }

# OutputTime { istep = 10 } stderr
OutputSimulation { start = 0 } sim-LEVEL-%g.txt { format = text }
OutputSimulation { start = 160 } sim-LEVEL-%g.txt { format = text }
OutputSimulation { start = 175 } sim-LEVEL-%g.txt { format = text }
OutputSimulation { start = 220 } sim-LEVEL-%g.txt { format = text }
} {
dry = DRY
time_order = 2
}
```

• Lake-at-rest balance in an inclined domain with cut cells
• ```1 0 GfsRiver GfsBox GfsGEdge { x = 0.5 } {
PhysicalParams { L = 10 g = 9.81 }

Refine 5
RefineSolid 7
Solid (sphere(5.,0.,0.,2.0))

Init {} {
Zb = x/10.
P = 12.1/(100. - 4.*M_PI)
}

Time { end = 200 }
#    OutputSimulation { istep = 10 } stdout
SourceCoriolis 0 1.0e-01
OutputScalarNorm { start = end } u { v = U }
OutputErrorNorm { start = end } ep { v = P } {
s = MAX(0, 0.5 - x/10.)
unbiased = 1
relative = 1
}
GModule gfsview
OutputView { start = end } still.eps { format = EPS } still.gfv
EventScript { start = end } {
status=0
if awk '{if (\$9 > 1e-5) { print "u: " \$9 > "/dev/stderr"; exit (1); }}' < u ; then :
else
status=\$GFS_STOP;
fi
if awk '{if (\$9 > 4e-3) { print "ep: " \$9 > "/dev/stderr"; exit (1); }}' < ep ; then :
else
status=\$GFS_STOP;
fi
exit \$status
}
}
GfsBox {
left = Boundary
right = Boundary
top = Boundary
bottom = Boundary
}
```

• Lake-at-rest balance in an inclined domain with bipolar metric
• ```8 8 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
PhysicalParams { g = 9.81 }

# Bipolar coordinates. see http://en.wikipedia.org/wiki/Bipolar_coordinates
Metric M {
x = 10.*sinh(rx/2. + 1.)/(cosh(rx/2. + 1.) - cos (M_PI*ry/4.))
y = 10.*sin(M_PI*ry/4.)/(cosh(rx/2. + 1.) - cos (M_PI*ry/4.))
}
Refine 4

Init {} {
Zb = x/10.
P = 12.1/(100. - 4.*M_PI)
}

Time { end = 200 }
#    OutputSimulation { istep = 10 } stdout
SourceCoriolis 0 1.0e-01
OutputScalarNorm { start = end } u { v = U }
OutputErrorNorm { start = end } ep { v = P } {
s = MAX(0, 1.26 - x/10.)
unbiased = 1
relative = 1
}
GModule gfsview
OutputScalarSum { start = 0 } vol { v = Zb }
OutputView { start = end } p.gnu { format = Gnuplot } p.gfv
OutputView { start = end } mesh.gnu { format = Gnuplot } mesh.gfv
EventScript { start = end } {
status=0
if gnuplot <<EOF; then :
set term postscript eps lw 2 solid color
set output 'still.eps'
unset key
unset logscale
unset border
unset xtics
unset ytics
unset xlabel
unset ylabel
set size ratio -1
plot 'mesh.gnu' w l lc 0, 'p.gnu' w l lc 1
EOF
else
status=\$GFS_STOP;
fi
if awk '{if (\$9 > 8e-7) { print "u: " \$9 > "/dev/stderr"; exit (1); }}' < u ; then :
else
status=\$GFS_STOP;
fi
if awk '{if (\$9 > 3e-4) { print "ep: " \$9 > "/dev/stderr"; exit (1); }}' < ep ; then :
else
status=\$GFS_STOP;
fi
exit \$status
}
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
GfsBox {
left = Boundary
right = Boundary
}
1 2 top
2 3 top
3 4 top
4 5 top
5 6 top
6 7 top
7 8 top
8 1 top
```

• Terrain reconstruction
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
PhysicalParams { L = 8 }
GModule terrain
RefineTerrain LEVEL H {
basename = terrain,terrain-high
} TRUE
VariableTerrain T {
basename = terrain,terrain-high
}
Time { end = 0 }
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' >> error-t } { v = T } {
s = {
double r = sqrt (x*x + y*y);
return r*r/8.+cos(M_PI*r)/2.;
}
w = (fabs (x) < 3.8 && fabs (y) < 3.8)
relative = 1
v = Et
}
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' >> error-h } { v = H0 } {
s = {
double r = sqrt (x*x + y*y);
return r*r/8.+cos(M_PI*r)/2.;
}
w = (fabs (x) < 3.8 && fabs (y) < 3.8)
relative = 1
v = Eh
}
OutputSimulation { start = end } {
awk '{
if (\$1 != "#") {
r = sqrt (\$1*\$1 + \$2*\$2);
print r,\$4,\$5,\$6,\$7;
}
}' | sort -n -k 1,2 > end-LEVEL.txt
} { format = text variables = T,H0,Et,Eh }
OutputSimulation { start = end } end-LEVEL.gfs
}
```

• Circular dam break on a sphere
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
PhysicalParams { L = LENGTH }
MetricLonLat M 1.
Refine 8
InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
Init {} { P = 0.2 + 1.8*P/LENGTH }
Time { end = 0.9 }
OutputTime { istep = 10 } stderr
GModule gfsview
OutputView { start = 0.3 step = 0.3 } isolines-%g.eps {
format = EPS line_width = 0.5
} isolines.gfv
OutputSimulation { step = 0.3 } sim-%g.txt { variables = U,V,P format = text }
OutputScalarSum { istep = 1 } sp { v = P }
#    OutputSimulation { istep = 10 } stdout
EventScript { start = end } {
status=0
for i in 0.3 0.6 0.9; do
if awk 'BEGIN { n1 = 0; pi = 3.14159265359 }{
if (\$1 != "#") {
c = cos(\$1*pi/180.)*cos(\$2*pi/180.);
d = atan2(sqrt(1. - c*c),c)*180./pi
i = int(d*2.)
x[i] += d
y[i] += \$6
n[i]++
x1[n1] = d;
y1[n1++] = \$6;
}
}END {
for (i = 0; i <= 180; i++)
if (n[i] > 0)
print x[i]/n[i], y[i]/n[i], n[i];
sum = 0.;
for (i = 0; i < n1; i++) {
j = int(x1[i]*2.)
d = y1[i] - y[j]/n[j];
sum += d*d;
}
scatter = sqrt(sum/n1);
if (scatter > 0.013) {
print scatter > "/dev/stderr";
exit (1);
}
}' < sim-\$i.txt > prof-\$i.txt ; then :
else
status=\$GFS_STOP;
fi

cat <<EOF | gnuplot
rdist(x,y)=acos(cos(x*pi/180.)*cos(y*pi/180.))*180./pi
cdist(x,y)=sqrt(x*x+y*y)
set xlabel 'Angular distance (degree)'
set ylabel 'Surface height'
set xtics 0,22.5,90
set ytics 0,0.25,0.75
set term postscript eps color lw 2 solid 20
set output 'p-\$i.eps'
plot [0:90][0:0.75]'sim-\$i.txt' u (rdist(\\$1,\\$2)):6 ps 0.25 pt 6 t '', 'prof-\$i.txt' w l lw 2 t ''
EOF

done
exit \$status
}
}
```

• Circular dam break on a rotating sphere
• ```1 0 GfsRiver GfsBox GfsGEdge {} {
PhysicalParams { L = LENGTH }
MetricLonLat M 1.
Refine 8
InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
Init {} { P = 0.2 + 1.8*P/LENGTH }
Time { end = 1.4 }
SourceCoriolis 10.*sin(y*M_PI/180.)
OutputTime { istep = 10 } stderr
GModule gfsview
OutputView { start = 0.4 step = 0.4 } isolines-%g.eps {
format = EPS line_width = 0.5
} isolines.gfv
OutputSimulation { step = 0.4 } sim-%g.txt { variables = U,V,P format = text }
#    OutputSimulation { istep = 10 } stdout
}
```

• Circular dam break on a ``cubed sphere''
• ```6 12 GfsRiver GfsBox GfsGEdge {} {
PhysicalParams { L = 2.*M_PI/4. }
MetricCubed M 8
Refine 8
InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
Init {} { P = 0.2 + 1.8*P/LENGTH }
Time { end = 1.5 }
AdaptGradient { istep = 1 } { cmax = 1e-2 maxlevel = 8 } P
#    OutputTime { istep = 10 } stderr
GModule gfsview
OutputView { start = 0.3 step = 0.3 } isolines-%g.gnu {
format = Gnuplot
} isolines.gfv
OutputSimulation { step = 0.3 } sim-%g.txt { variables = U,V,P format = text }
OutputScalarSum { istep = 1 } sp { v = P }
#    OutputSimulation { istep = 10 } stdout
EventScript { start = end } {
status=0
for i in 0.3 0.6 0.9 1.2 1.5; do
if awk 'BEGIN { n1 = 0; pi = 3.14159265359 }{
if (\$1 != "#") {
c = cos(\$1*pi/180.)*cos(\$2*pi/180.);
d = atan2(sqrt(1. - c*c),c)*180./pi
i = int(d*2.)
x[i] += d
y[i] += \$6
n[i]++
x1[n1] = d;
y1[n1++] = \$6;
}
}END {
for (i = 0; i <= 180; i++)
if (n[i] > 0)
print x[i]/n[i], y[i]/n[i], n[i];
sum = 0.;
for (i = 0; i < n1; i++) {
j = int(x1[i]*2.)
d = y1[i] - y[j]/n[j];
sum += d*d;
}
scatter = sqrt(sum/n1);
if (scatter > 0.015) {
print scatter > "/dev/stderr";
exit (1);
}
}' < sim-\$i.txt > prof-\$i.txt ; then :
else
status=\$GFS_STOP;
fi

gnuplot <<EOF
set term postscript eps color lw 1 18
set output 'isolines-\$i.eps'
set size ratio -1
set xtics -90,30,90
set ytics -90,30,90
unset key
plot [-90:90][-90:90]'isolines-\$i.gnu' w l

set term postscript eps color lw 2 solid 20
set output 'p-\$i.eps'
rdist(x,y)=acos(cos(x*pi/180.)*cos(y*pi/180.))*180./pi
cdist(x,y)=sqrt(x*x+y*y)
set xlabel 'Angular distance (degree)'
set ylabel 'Surface height'
set xtics 0,22.5,90
set ytics 0,0.25,0.75
plot [0:90][0:0.75]'sim-\$i.txt' u (rdist(\\$1,\\$2)):6 ps 0.25 pt 6 t '', 'prof-\$i.txt' w l lw 2 t ''
EOF

done
exit \$status
}
}
```

• Rossby--Haurwitz wave with Saint-Venant
• ```6 12 GfsRiver GfsBox GfsGEdge {} {
Global {
#define AR 6.37122e6
#define N 4.
#define Umax 50.
#define M (Umax/(N*AR))
#define K M

#define Omega 7.292e-5
#define DTR (M_PI/180.)
#define H0 8e3
#define G 9.81

// Williamson 1992, eq. (137)
static double u0 (double lambda, double phi) {
double cosphi = cos (phi), sinphi = sin (phi);
return AR*M*cosphi + AR*K*pow (cosphi, N - 1)*
cos (N*lambda)*(N*sinphi*sinphi - cosphi*cosphi);
}

// Williamson 1992, eq. (138)
static double v0 (double lambda, double phi) {
double cosphi = cos (phi), sinphi = sin (phi);
return - AR*K*N*pow (cosphi, N - 1)*sinphi*sin (N*lambda);
}

// Williamson 1992, eq. (139)
static double vorticity0 (double lambda, double phi) {
return 2.*M*sin(phi) - K*sin(phi)*pow (cos (phi), N)*(N*N + 3.*N + 2.)*cos (N*lambda);
}

// Williamson 1992, eqs. (140-143)
static double p0 (double lambda, double phi, double t) {
double nu = (N*(3. + N)*M - 2.*Omega)/((1. + N)*(2. + N));
lambda -= nu*t;
double cosphi = cos (phi);
double Aphi = M/2.*(2.*Omega + M)*cosphi*cosphi + K*K/4.*pow (cosphi, 2.*N)*
((N + 1.)*cosphi*cosphi + 2.*N*N - N - 2. - 2.*N*N*pow(cosphi, -2.));
double Bphi = 2.*(Omega + M)*K/((N + 1.)*(N + 2.))*pow(cosphi, N)*
(N*N + 2.*N + 2. - (N + 1.)*(N + 1.)*cosphi*cosphi);
double Cphi = K*K/4.*pow(cosphi,2.*N)*((N + 1.)*cosphi*cosphi - (N + 2.));
return AR*AR*(Aphi + Bphi*cos (N*lambda) + Cphi*cos (2.*N*lambda));
}
}
PhysicalParams { L = 2.*M_PI*AR/4. g = G }
MetricCubed M LEVEL
SourceCoriolis 2.*Omega*sin(y*DTR)

Init {} {
P = H0 + p0(x*DTR,y*DTR,t)/G
(U,V) = (u0(x*DTR,y*DTR)*P, v0(x*DTR,y*DTR)*P)
}

Refine LEVEL

}

# ~24 days
Time { end = 2073534 dtmax = 1e3 }
#    OutputTime { istep = 10 } stderr

OutputScalarNorm { istep = 10 } v-LEVEL { v = V }
OutputScalarSum { istep = 10 } ec-LEVEL { v = Velocity2 }
#    OutputScalarSum { istep = 10 } zeta-LEVEL { v = Vorticity }
OutputScalarSum { istep = 10 } p-LEVEL { v = P }
OutputErrorNorm { istep = 10 } eh-LEVEL { v = P } {
s = p0(x*DTR,y*DTR,t)/G
v = EH unbiased = 1 relative = 1
}
OutputSimulation { start = end } end-LEVEL.gfs
#    OutputSimulation { istep = 10 } stdout
#    OutputPPM { istep = 10 } eh.ppm { v = EH }
#    OutputPPM { istep = 10 } p.ppm { v = P }

GModule gfsview
OutputView { start = end } ehp-LEVEL.gnu { format = Gnuplot } ehp.gfv
OutputView { start = end } ehm-LEVEL.gnu { format = Gnuplot } ehm.gfv
OutputView { start = end } h-LEVEL.gnu { format = Gnuplot } h.gfv
OutputView { start = end } href-LEVEL.gnu { format = Gnuplot } href.gfv
}
```