GfsRiver

(Difference between revisions)
 Revision as of 00:00, 14 June 2012Popinet (Talk | contribs) (→Bibliography)← Previous diff Revision as of 23:03, 31 July 2012Popinet (Talk | contribs) (info on multilayers)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. + 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. Multiple layers can be added using [[GfsLayers]]. The default variables are: The default variables are: - ;U: horizontal component of the flux (dimension L^2/T) + ;U: horizontal component of the (depth-averaged) flux (dimension L^2/T) - ;V: vertical component of the flux + ;V: vertical component of the (depth-averaged) flux ;P: fluid depth ;P: fluid depth ;Zb: bed elevation ;Zb: bed elevation Line 16: Line 16: time_order = 2 time_order = 2 dry = 1e-6 dry = 1e-6 + k = 1 + nu = 1e-2 + dut = 1 } } GfsBox {} GfsBox {} Line 24: Line 27: ; time_order: the order of the temporal discretisation, either 1 or 2. The default is 2. ; time_order: the order of the temporal discretisation, either 1 or 2. The default is 2. ; dry: the depth threshold below which a cell is considered to be "dry". The default is 1e-6 times the box length scale (L parameter in [[GfsPhysicalParams]]). ; dry: the depth threshold below which a cell is considered to be "dry". The default is 1e-6 times the box length scale (L parameter in [[GfsPhysicalParams]]). + ;k: Navier friction coefficient. The default is zero. The boundary condition on the bottom is ∂t U = ... - kU/P. + ;nu: for more than one layer (see [[GfsLayers]]), the vertical kinematic viscosity. Default is zero. + ;dut: for more than one layer, the vertical derivative of the horizontal velocity on the top layer. Default is zero. == Slope limiters == == Slope limiters == Line 37: Line 43: == See also == == See also == - *[[Karamea flood tutorial]] + *[[GfsLayers]] — Setting the number of layers for a multi-layer Saint-Venant model - *[[11th March Japanese tsunami]] + *[[GfsDischargeElevation]] — Elevation for a given discharge *[[GfsDischargeElevation]] — Elevation for a given discharge *[[GfsSourcePipe]] — Pipes and culverts models for GfsRiver *[[GfsSourcePipe]] — Pipes and culverts models for GfsRiver + *[[GfsSourceCulvert]] — Culvert model based on Boyd, 1987 + *[[Karamea flood tutorial]] + *[[11th March Japanese tsunami]] == References == == References ==

Revision as of 23:03, 31 July 2012

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. Multiple layers can be added using GfsLayers.

The default variables are:

`U`
horizontal component of the (depth-averaged) flux (dimension L^2/T)
`V`
vertical component of the (depth-averaged) 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).

The syntax in parameter files is

```[ GfsSimulation ] {
time_order = 2
dry = 1e-6
k = 1
nu = 1e-2
dut = 1
}
GfsBox {}
...
```

where the parameter block is optional. The optional parameters are:

`time_order`
the order of the temporal discretisation, either 1 or 2. The default is 2.
`dry`
the depth threshold below which a cell is considered to be "dry". The default is 1e-6 times the box length scale (`L` parameter in GfsPhysicalParams).
`k`
Navier friction coefficient. The default is zero. The boundary condition on the bottom is ∂t `U` = ... - `k``U/P`.
`nu`
for more than one layer (see GfsLayers), the vertical kinematic viscosity. Default is zero.
`dut`
for more than one layer, the vertical derivative of the horizontal velocity on the top layer. Default is zero.

Slope limiters

GfsRiver uses slope-limiting to deal with shocks while preserving second-order spatial accuracy away from shocks. By default GfsRiver uses a minmod limiter which is the most dissipative (and hence very stable) in the classical family of TVD limiters. This is appropriate for shallow flows with friction parameterisation (e.g. rivers) but can be overly dissipative for long-wave propagation in deep water (e.g. tsunamis). A better option for such flows is the Sweby limiter which offers a good compromise between stability and dissipation. This limiter can be selected using:

```AdvectionParams { gradient = gfs_center_sweby_gradient }
```

Note also that the overall scheme is stable only for CFL numbers smaller than 0.5 (independently of the choice of limiter).

The scheme can also be made first-order in space by setting `gradient` to `none`.

References

Ocean Dynamics 61(9):1261-1285, 2011
Bibtex

Hyunuk An, Soonyoung Yu - Well-balanced shallow water flow simulation on quadtree cut cell grids

Advances in Water Resources 39:60-70, 2012
Bibtex

Natural Hazards and Earth System Sciences , 2012
Bibtex

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
}
```