GfsRiver

From Gerris

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

    # Set a solid boundary close to the top boundary to limit the
    # domain width to one cell (i.e. a 1D domain)
    Solid (y/8. + 1./pow(2,9) - 1e-6 - 0.5)

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

    # Save a text-formatted simulation
    OutputSimulation { step = 0.5 } sim-%g.txt { format = text }

    # Use gnuplot to create gif images
    EventScript { step = 0.5 } {
	time=`echo $GfsTime | awk '{printf("%4.1f\n", $1);}'`
	cat <<EOF | gnuplot
load 'dam.plot'
set title "t = $time"
set term postscript eps color 14
set output "sim.eps"
plot [-4.:4.]'sim-$GfsTime.txt' u 1:7:8 w filledcu lc 3, 'sim-0.txt' u 1:7 w l lw 4 lc 1 lt 1
EOF
	time=`echo $GfsTime | awk '{printf("%04.1f\n", $1);}'`
	convert -density 300 sim.eps -trim +repage -bordercolor white -border 10 -resize 640x282! sim-$time.gif
	rm -f sim.eps
    }

    # Combine all the gif images into a gif animation using gifsicle
    EventScript { start = end } {
	gifsicle --colors 256 --optimize --delay 25 --loopcount=0 sim-*.gif > dam.gif
	rm -f sim-*.gif sim-*.txt
    }
}

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

• 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 }
    RefineSolid LEVEL
    Solid (y/10000. + 1./pow(2,LEVEL) - 1e-6 - 0.5)
    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)
    }
    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 { istep = 10 } vol-LEVEL { v = P }
    OutputScalarSum { istep = 10 } U-LEVEL { v = U }
    OutputErrorNorm { istep = 1 } error-LEVEL { v = P } {
	s = Pt
	v = DP
    }
}

• 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
    OutputSimulation { step = 0.3 } sim-%g.gfs
    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.012) {
                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

	    echo "Save stdout { format = EPS line_width = 0.5 }" | \
		gfsview-batch2D sim-$i.gfs isolines.gfv > isolines-$i.eps
	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
    OutputSimulation { step = 0.4 } sim-%g.gfs
    OutputSimulation { step = 0.4 } sim-%g.txt { variables = U,V,P format = text }
#    OutputSimulation { istep = 10 } stdout
    EventScript { start = end } {
	for i in 0.4 0.8 1.2; do
	    echo "Save stdout { format = EPS line_width = 0.5 }" | \
		gfsview-batch2D sim-$i.gfs isolines.gfv > isolines-$i.eps
	done
    }
}