[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:57 UTC 2009
The following commit has been merged in the upstream branch:
commit 5b46073bc7927ab9dda739b9fdbc43032723be06
Author: Stephane Popinet <popinet at users.sf.net>
Date: Sat Jan 10 12:39:03 2009 +1100
New 'parabola' test case for GfsRiver
darcs-hash:20090110013903-d4795-6f3cf66d4efbe83a77d3317e6621d4c3e71a8bd3.gz
diff --git a/test/Makefile.am b/test/Makefile.am
index 1b91588..0b9d2a8 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -24,7 +24,8 @@ TESTDIRS = \
oscillation \
geo \
waves \
- nz
+ nz \
+ parabola
EXTRA_DIST = \
template.tex \
diff --git a/test/parabola/error.ref b/test/parabola/error.ref
new file mode 100644
index 0000000..5f68662
--- /dev/null
+++ b/test/parabola/error.ref
@@ -0,0 +1,5 @@
+5 0.00230918 0.00460721 0.04365
+6 0.00100378 0.00197881 0.0257
+7 0.0006263 0.00109614 0.01443
+8 0.000685491 0.00122904 0.02707
+9 0.000773271 0.00146527 0.0418
diff --git a/test/parabola/order.ref b/test/parabola/order.ref
new file mode 100644
index 0000000..ddc2c59
--- /dev/null
+++ b/test/parabola/order.ref
@@ -0,0 +1,4 @@
+6 1.3963 1.35575 0.829894
+7 0.496801 0.691639 0.348922
+8 0.00249267 0.0323125 -0.604711
+9 -0.188674 -0.275253 -0.454448
diff --git a/test/parabola/parabola.gfs b/test/parabola/parabola.gfs
new file mode 100644
index 0000000..17063b1
--- /dev/null
+++ b/test/parabola/parabola.gfs
@@ -0,0 +1,120 @@
+# Title: Oscillations in a parabolic container
+#
+# Description:
+#
+# Analytical solutions for the damped oscillations of a liquid in a
+# parabolic container have been derived by Sampson et al
+# \cite{sampson2006,liang2008}. Figure \ref{elevation} illustrates the
+# numerical and analytical solutions at $t = 1500$ seconds. Wetting
+# and drying occur at the two moving contact points and hydrostatic
+# balance is approached as time passes.
+#
+# Figure \ref{u0} gives the analytical and numerical solutions for the
+# horizontal component of velocity (which is spatially
+# constant). Figures \ref{error} and \ref{error-u} give the relative
+# errors in surface elevation and horizontal velocity respectively, as
+# functions of spatial resolution. Figures \ref{order} and
+# \ref{order-u} give the corresponding convergence orders.
+#
+# The errors are generally small, however convergence is not reached
+# with increasing resolution. This is due to errors in velocity at the
+# contact points.
+#
+# \begin{figure}[htbp]
+# \caption{\label{elevation}Solution at $t = 1500$ seconds. Six levels of refinement.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{elevation.eps}
+# \end{center}
+# \end{figure}
+#
+# \begin{figure}[htbp]
+# \caption{\label{u0}Time evolution of the (spatially constant)
+# horizontal velocity. Seven levels of refinement.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{u0.eps}
+# \end{center}
+# \end{figure}
+#
+# \begin{figure}[htbp]
+# \caption{\label{error}Evolution of the relative elevation error norms as functions of resolution.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{error.eps}
+# \end{center}
+# \end{figure}
+#
+# \begin{figure}[htbp]
+# \caption{\label{order}Corresponding convergence order.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{order.eps}
+# \end{center}
+# \end{figure}
+#
+# \begin{figure}[htbp]
+# \caption{\label{error-u}Evolution of the relative velocity error
+# L2-norm as a function of resolution.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{error-u.eps}
+# \end{center}
+# \end{figure}
+#
+# \begin{figure}[htbp]
+# \caption{\label{order-u}Corresponding convergence order.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{order-u.eps}
+# \end{center}
+# \end{figure}
+#
+# Author: St\'ephane Popinet
+# Command: sh parabola.sh
+# Version: 1.3.1
+# Required files: parabola.sh error.ref order.ref
+# Generated files: elevation.eps error.eps error-u.eps order.eps order-u.eps u0.eps
+#
+Define h0 10.
+Define a 3000.
+Define tau 1e-3
+Define B 5
+Define G 9.81
+
+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 }
+ Init { istep = 1 } {
+ U = (P <= 0. ? 0. : U/(1. + dt*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
+ }
+}
+GfsBox {
+ left = Boundary
+ right = Boundary
+}
diff --git a/test/parabola/parabola.sh b/test/parabola/parabola.sh
new file mode 100644
index 0000000..3150a77
--- /dev/null
+++ b/test/parabola/parabola.sh
@@ -0,0 +1,162 @@
+levels="5 6 7 8 9"
+
+if ! $donotrun; then
+ for level in $levels; do
+ if gerris2D -DLEVEL=$level parabola.gfs; then :
+ else
+ exit 1
+ fi
+ done
+fi
+
+for level in $levels; do
+ awk -v level=$level 'BEGIN{
+ s1 = 0.;
+ s2 = 0.;
+ smax = 0.;
+ n = 0;
+ h0 = 10.;
+ }{
+ n++;
+ s1 += $5;
+ s2 += $7*$7;
+ if ($9 > smax) smax = $9;
+ }END { print level, s1/n/h0, sqrt(s2/n)/h0, smax/h0; }' < error-$level
+done > error
+
+for level in $levels; do
+ if paste U-$level vol-$level | awk -v level=$level 'BEGIN {
+ sum = 0.;
+ n = 0;
+
+ h0 = 10.
+ a = 3000.
+ tau = 1e-3
+ B = 5.
+ G = 9.81
+ p = sqrt (8.*G*h0)/a
+ s = sqrt (p*p - tau*tau)/2.
+ } {
+ u0 = $5/$10;
+ t = $3;
+ ref = B*exp (-tau*t/2.)*sin (s*t);
+ sum += (u0 - ref)*(u0 - ref);
+ n += 1;
+ }
+ END {
+ print level, sqrt (sum/n);
+ }'; then
+ :
+ else
+ exit 1;
+ fi
+done > error-u
+
+if awk '
+BEGIN { n = 0 }
+{
+ l[n] = $1; n1[n] = $2; n2[n] = $3; ni[n++] = $4;
+}
+END {
+ for (i = 1; i < n; i++)
+ print l[i] " " log(n1[i-1]/n1[i])/log(2.) " " log(n2[i-1]/n2[i])/log(2.) " " log(ni[i-1]/ni[i])/log(2.);
+}' < error > order; then :
+else
+ exit 1
+fi
+
+if awk '
+BEGIN { n = 0 }
+{
+ l[n] = $1; n2[n++] = $2;
+}
+END {
+ for (i = 1; i < n; i++)
+ print l[i] " " log(n2[i-1]/n2[i])/log(2.);
+}' < error-u > order-u; then :
+else
+ exit 1
+fi
+
+if cat <<EOF | gnuplot ; then :
+ set term postscript eps color lw 3 solid 20
+
+ set output 'error.eps'
+ set xlabel 'Level'
+ set ylabel 'Relative error norms'
+ set key
+ set logscale y
+ set xtics 0,1
+ set grid
+ plot 'error.ref' u 1:2 t '1 (ref)' w lp, \
+ 'error.ref' u 1:3 t '2 (ref)' w lp, \
+ 'error.ref' u 1:4 t 'max (ref)' w lp, \
+ 'error' u 1:2 t '1' w lp, \
+ 'error' u 1:3 t '2' w lp, \
+ 'error' u 1:4 t 'max' w lp
+
+ set output 'error-u.eps'
+ set ylabel 'u0 relative error L2 norm'
+ plot 'error-u' u 1:(\$2/4.) t '' w lp
+
+ set output 'order.eps'
+ set xlabel 'Level'
+ set ylabel 'Order'
+ unset logscale
+ set ytics 0,1
+ plot [][-1:2] 'order.ref' u 1:2 t '1 (ref)' w lp, \
+ 'order.ref' u 1:3 t '2 (ref)' w lp, \
+ 'order.ref' u 1:4 t 'max (ref)' w lp, \
+ 'order' u 1:2 t '1' w lp, \
+ 'order' u 1:3 t '2' w lp, \
+ 'order' u 1:4 t 'max' w lp
+
+ set output 'order-u.eps'
+ plot [][-1:2]'order-u' u 1:2 t '' w lp
+
+ set output 'u0.eps'
+
+ h0 = 10.
+ a = 3000.
+ tau = 1e-3
+ B = 5.
+ G = 9.81
+ p = sqrt (8.*G*h0)/a
+ s = sqrt (p*p - tau*tau)/2.
+ u0(t) = B*exp (-tau*t/2.)*sin (s*t)
+
+ set xtics auto
+ set ytics auto
+ unset grid
+ set ylabel 'u0'
+ set xlabel 'Time'
+ plot u0(x) t 'Analytical', '< paste U-7 vol-7' u 3:(\$5/\$10) every 2 w p t 'Gerris'
+
+ set output 'elevation.eps'
+ set xlabel 'x (m)'
+ set ylabel 'z (m)'
+ t = 1500
+ psi(x) = 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 + h0
+ bed(x) = h0*(x/a)**2
+ set key top center
+ plot [-5000:5000] \
+ 'sim-6-1500.txt' u 1:8:(\$4+\$8) w filledcu lc 3 t 'Numerical', \
+ psi(x) > bed(x) ? psi(x) : bed(x) lc 2 t 'Analytical', \
+ bed(x) lw 3 lc 1 lt 1 t 'Bed profile'
+EOF
+else
+ exit 1
+fi
+
+if cat <<EOF | python ; then :
+from check import *
+from sys import *
+if (Curve('error',1,4) - Curve('error.ref',1,4)).max() > 1e-4:
+ print (Curve('error',1,4) - Curve('error.ref',1,4)).max()
+ exit(1)
+EOF
+else
+ exit 1
+fi
diff --git a/test/template.tex b/test/template.tex
index 9db0d2c..65f91f1 100644
--- a/test/template.tex
+++ b/test/template.tex
@@ -99,6 +99,10 @@ branch only.
\input{waves/adaptive/adaptive.tex}
\input{nz/nz.tex}
+\section{Saint-Venant}
+
+\input{parabola/parabola.tex}
+
\bibliographystyle{plain}
\bibliography{tests}
diff --git a/test/tests.bib b/test/tests.bib
index 5298c44..c29b3d0 100644
--- a/test/tests.bib
+++ b/test/tests.bib
@@ -109,6 +109,14 @@
pages = {1931-1951}
}
+ at Article{liang2008,
+ title = {Adaptive quadtree simulation of shallow flows with wet--dry fronts over complex topography},
+ author = {Q. Liang and A.G.L. Borthwick},
+ journal= {Computers and Fluids},
+ year = {2008},
+ publisher = {Elsevier}
+}
+
@Article{lynch87,
author = {D. R. Lynch and F. E. Werner},
title = {3-{D} hydrodynamics on finite elements. {P}art I: linearized harmonic model},
@@ -207,6 +215,14 @@
url = {http://marine.rutgers.edu/po/index.php?model=test-problems&title=circle}
}
+ at Article{sampson2006,
+ title = {Moving boundary shallow water flow above parabolic bottom topography},
+ author = {J. Sampson and A. Easton and M. Singh},
+ journal = {ANZIAM J},
+ volume = {47},
+ year = {2006}
+}
+
@Article{torres00,
author = {D. J. Torres and J. U. Brackbill},
title = {The Point-Set method: front-tracking without connectivity},
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list