GfsRiver

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 00:00, 14 June 2012
Popinet (Talk | contribs)
(Bibliography)
← Previous diff
Revision as of 04:35, 14 June 2012
Popinet (Talk | contribs)
(See also)
Next diff →
Line 41: Line 41:
*[[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
== References == == References ==

Revision as of 04:35, 14 June 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.

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

The syntax in parameter files is

[ GfsSimulation ] {
  time_order = 2
  dry = 1e-6
}
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).

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
    }
    

Personal tools
communication