[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:54:36 UTC 2009
The following commit has been merged in the upstream branch:
commit 2a88393a8d850e118c100c482db227d701e4708d
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Sep 6 12:05:28 2007 +1000
New inviscid droplet oscillation test case
darcs-hash:20070906020528-d4795-2056bdd047981e8955e971b822b0ffc5d1556a6c.gz
diff --git a/test/Makefile.am b/test/Makefile.am
index 6b087b5..49319fe 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -19,7 +19,8 @@ TESTDIRS = \
dumbell \
nz \
shear \
- kinetic
+ kinetic \
+ oscillation
EXTRA_DIST = \
template.tex \
diff --git a/test/oscillation/fit.ref b/test/oscillation/fit.ref
new file mode 100644
index 0000000..9a5f962
--- /dev/null
+++ b/test/oscillation/fit.ref
@@ -0,0 +1,4 @@
+5 0.000273799399566086 5.49284377610436 160.012068757225
+6 0.000276489458821021 3.15635419121691 152.638147630339
+7 0.000289109113132265 2.37334656363034 154.294252244876
+8 0.000291520326836409 1.64059036642349 154.473757484488
diff --git a/test/oscillation/oscillation.gfs b/test/oscillation/oscillation.gfs
new file mode 100644
index 0000000..45f954e
--- /dev/null
+++ b/test/oscillation/oscillation.gfs
@@ -0,0 +1,115 @@
+# Title: Shape oscillation of an inviscid droplet
+#
+# Description:
+#
+# A two-dimensional elliptical droplet (density ratio 1/1000) is
+# released in a large domain. Under the effect of surface-tension the
+# shape of the droplet oscillates around its (circular) equilibrium
+# shape. The fluids inside and outside the droplet are inviscid so
+# ideally no damping of the oscillations should occur. As illustrated
+# on Figure \ref{kinetic} some damping occurs in the simulation due to
+# numerical dissipation.
+#
+# This simulation is also a stringent test case of the accuracy of the
+# surface tension representation as no explicit viscosity can damp
+# eventual spurious currents. Onset of spurious-current-driven
+# instability can be seen on the coarsest mesh (in Figure \ref{kinetic}
+# the kinetic energy behaves strangely for the 32x32 resolution).
+#
+# \begin{figure}[htbp]
+# \caption{\label{kinetic}Evolution of the kinetic energy as a function
+# of time for the spatial resolutions indicated in the legend. The
+# black lines are fitted decreasing exponential functions.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{k.eps}
+# \end{center}
+# \end{figure}
+#
+# The initial shape of the droplet is given by:
+# $$r(\theta) = r_0 + \alpha\cos(n\theta)$$
+# The oscillation frequency is then \cite{torres00}:
+# $$\omega_n^2={(n^3-n)\sigma\over (\rho_d+\rho_e)r_0^3}$$
+#
+# A comparison between the theoretical and numerical values of the
+# frequency is given in Figure \ref{frequency}.
+#
+# \begin{figure}[htbp]
+# \caption{\label{frequency}Relative error in the oscillation
+# frequency as a function of resolution.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{frequency.eps}
+# \end{center}
+# \end{figure}
+#
+# The amount of numerical damping can be estimated by computing an
+# equivalent viscosity. With viscosity, kinetic energy is expected to
+# decrease as:
+# $$\exp(-C\nu/D^2t)$$
+# where $C$ is a constant, $\nu$ the viscosity and $D$ the droplet
+# diameter. Using curve fitting the damping coefficient $b=C\nu/D^2$
+# can be estimated (black curves on Figure \ref{kinetic}). An
+# equivalent Laplace number can the be computed as:
+# $$La={\sigma D\over \rho\nu^2}={\sigma C^2 \over \rho b^2 D^3}$$
+# The equivalent Laplace number depends on spatial resolution as
+# illustrated in Figure \ref{laplace}.
+#
+# \begin{figure}[htbp]
+# \caption{\label{laplace}Equivalent Laplace number estimated from the
+# numerical damping of kinetic energy.}
+# \begin{center}
+# \includegraphics[width=0.8\hsize]{laplace.eps}
+# \end{center}
+# \end{figure}
+#
+# Author: St\'ephane Popinet
+# Command: sh oscillation.sh
+# Version: 1.1.1
+# Required files: oscillation.sh fit.ref
+# Running time: 3 minutes
+# Generated files: frequency.eps k.eps laplace.eps
+#
+Define Diameter 0.2
+Define Epsilon 0.05
+Define Var(T,min,max) (min + CLAMP(T,0,1)*(max - min))
+Define Rho(T) Var(T, 1., 1e-3)
+Define Radius(x,y) (Diameter/2.*(1. + Epsilon*cos (2.*atan2 (y, x))))
+
+1 0 GfsSimulation GfsBox GfsGEdge {} {
+
+ Time { end = 1 }
+
+ Refine LEVEL
+
+ VariableTracerVOF T
+ VariableFiltered T1 T 1
+ InitFraction T ({ x += 0.5; y += 0.5; return x*x + y*y - Radius(x,y)*Radius(x,y); })
+
+ PhysicalParams { alpha = 1./Rho(T1) }
+ VariableCurvature K T
+ SourceTension T 1. K
+
+ AdaptFunction { istep = 1 } {
+ cmax = 0.01
+ maxlevel = LEVEL
+ } (T > 0 && T < 1 ? 1. : fabs (Vorticity)*ftt_cell_size (cell))
+
+ OutputScalarSum { istep = 1 } k-LEVEL {
+ v = Rho(T1)*Velocity2
+ }
+
+ EventScript { start = end } {
+ cat <<EOF | gnuplot 2>&1 | awk '{if ($1 == "result:") print LEVEL,$2,$3,$4;}'
+ k(t)=a*exp(-b*t)*(cos(c*t+3.14159265359)+1.)
+ a = 3e-4
+ b = 1.5
+ c = 153
+ fit k(x) 'k-LEVEL' u 3:5 via a,b,c
+ print "result: ", a, b, c
+EOF
+ rm -f fit.log
+ }
+}
+GfsBox {
+ left = Boundary
+ bottom = Boundary
+}
diff --git a/test/oscillation/oscillation.sh b/test/oscillation/oscillation.sh
new file mode 100644
index 0000000..c02858d
--- /dev/null
+++ b/test/oscillation/oscillation.sh
@@ -0,0 +1,68 @@
+levels="5 6 7 8"
+
+if ! $donotrun; then
+ rm -f fit
+ for level in $levels; do
+ if gerris2D -D LEVEL=$level oscillation.gfs >> fit; then :
+ else
+ exit 1
+ fi
+ done
+fi
+
+rm -f fit-*
+if awk '{
+ level = $1; a = $2; b = $3; c = $4;
+ for (t = 0; t <= 1.; t += 0.005)
+ print t, 2.*a*exp(-b*t) >> "fit-" level;
+}' < fit ; then :
+else
+ exit 1
+fi
+
+if cat <<EOF | gnuplot ; then :
+ set term postscript eps color lw 3 solid 20
+
+ D = 0.2
+ n = 2.
+ sigma = 1.
+ rhol = 1.
+ rhog = 1./1000.
+ r0 = 0.1
+ omega0 = sqrt((n**3-n)*sigma/((rhol+rhog)*r0**3))
+
+ set output 'k.eps'
+ set xlabel 'Time'
+ set ylabel 'Kinetic energy'
+ set logscale y
+ plot [0:1][1e-5:]'k-8' u 3:5 t "256x256" w l, 'k-7' u 3:5 t "128x128" w l, 'k-6' u 3:5 t "64x64" w l, 'k-5' u 3:5 t "32x32" w l, 'fit-8' t "fit" w l lt 7, 'fit-7' t "" w l lt 7, 'fit-6' t "" w l lt 7, 'fit-5' t "" w l lt 7
+
+ set output 'laplace.eps'
+ set xlabel 'Diameter (grid points)'
+ set ylabel 'Equivalent Laplace number'
+ unset logscale
+ empirical_constant = 30.
+ plot 'fit' u (D*2.**(\$1)):(1./(\$3**2.*D**3.))*empirical_constant**2. t "" w lp pt 5 ps 2
+
+ set output 'frequency.eps'
+ set xlabel 'Diameter (grid points)'
+ set ylabel 'Frequency error (%)'
+ set xzeroaxis
+ plot 'fit' u (D*2.**(\$1)):(\$4/2./omega0-1.)*100. t "" w lp pt 5 ps 2
+
+
+EOF
+else
+ exit 1
+fi
+
+if cat <<EOF | python ; then :
+from check import *
+from sys import *
+if (Curve('fit',1,3) - Curve('fit.ref',1,3)).max() > 1e-2 or\
+ (Curve('fit',1,4) - Curve('fit.ref',1,4)).max() > 1e-2:
+ exit(1)
+EOF
+else
+ exit 1
+fi
diff --git a/test/template.tex b/test/template.tex
index 24f0240..e27f6e7 100644
--- a/test/template.tex
+++ b/test/template.tex
@@ -72,6 +72,7 @@ current stable branch of the version-controlled source code.
\input{capwave/capwave.tex}
\input{capwave/density/density.tex}
\input{capwave/gravity/gravity.tex}
+\input{oscillation/oscillation.tex}
\section{Shallow-water}
diff --git a/test/tests.bib b/test/tests.bib
index 9f43322..a981f71 100644
--- a/test/tests.bib
+++ b/test/tests.bib
@@ -167,6 +167,15 @@
url = {http://marine.rutgers.edu/po/index.php?model=test-problems&title=circle}
}
+ at Article{torres00,
+ author = {D. J. Torres and J. U. Brackbill},
+ title = {The Point-Set method: front-tracking without connectivity},
+ journal = {J. Comput. Phys.},
+ year = 2000,
+ volume = 165,
+ pages = {620-644}
+}
+
@Article{vola2004,
author = {D. Vola and F. Babik and J.-C Latch\'e},
title = {On a numerical strategy to compute gravity currents of non-Newtonian fluids},
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list