GfsRiver

From Gerris

Jump to: navigation, search

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
x-component of the (depth-averaged) flux (dimension L^2/T)
V
y-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 = [ GfsFunction ]
  nu = [ GfsFunction ]
  dut = [ GfsFunction ]
}
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 = ... - 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

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.

See also

References

S. Popinet - Quadtree-adaptive tsunami modelling

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

S. Popinet - Adaptive modelling of long-distance wave propagation and fine-scale flooding during the Tohoku tsunami

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.
    AdvectionParams {
       # gradient = gfs_center_minmod_gradient
	gradient = none
    }

    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 "
          load 'dam.plot'
          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 }
    AdaptGradient { istep = 1 } { 
	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
    AdaptGradient { istep = 1 } {
        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
    AdvectionParams { gradient = gfs_center_sweby_gradient }

    # adapt down to 9 levels based on the slope of the (wet)
    # free-surface and with a tolerance of 1 mm
    AdaptGradient { istep = 1 } {
        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 } { 
	ppm2mpeg -s 400x600 > overhead.mpg 
    } { 
	width = 800 height = 1200
    } overhead.gfv
} {
    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
    MetricLonLat M RADIUS
    # 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
    GModule okada
    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
    AdvectionParams { gradient = gfs_center_sweby_gradient }

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

    # Segment 1
    Init {} { D = 0 }
    InitOkada D {
	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 }
	InitOkada D {
	    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 }
	InitOkada D {
	    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 }
	InitOkada D {
	    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 }
	InitOkada D {
	    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
    }

    # adapt on gradient of wet surface elevation
    # with a tolerance of 5 cm
    AdaptGradient { istep = 1 } {
	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.
}

  • Geostrophic adjustment with Saint-Venant

    • 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
  }
  AdvectionParams { gradient = gfs_center_gradient }
  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
    }
    AdvectionParams { 
	# Sweby seems to work better than minmod for this test case
	gradient = gfs_center_sweby_gradient 
    }
    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
    }
    AdvectionParams { 
	# Sweby seems to oscillate in this case
#	gradient = gfs_center_sweby_gradient 
    }
    # 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 }

    AdvectionParams {
	gradient = gfs_center_sweby_gradient
    }

    # 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

    AdvectionParams { 
	gradient = gfs_center_gradient 
#	gradient = none
    }

    # ~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
}

Retrieved from "http://gfs.sourceforge.net/wiki/index.php/GfsRiver"
Views
Navigation
communication