GfsRiver

From Gerris

Revision as of 23:43, 14 May 2009 Popinet (Talk | contribs) ← Previous diff		Revision as of 04:02, 28 August 2009 Popinet (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:
		+
		+	;<code>U</code>: horizontal component of the flux (dimension L^2/T)
		+	;<code>V</code>: vertical component of the flux
		+	;<code>P</code>: fluid depth
		+	;<code>Zb</code>: bed elevation
		+	;<code>H</code>: <code>Zb + P</code> i.e. the elevation of the free surface
		+
		+	The value of the acceleration of gravity is set using [[GfsPhysicalParams]] (default to 1).
		+
	<examples/>		<examples/>

---

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