GfsAdvection
From Gerris
Examples
- Convergence of the Godunov advection scheme
- Time-reversed VOF advection in a shear flow
- Time-reversed advection of a VOF concentration
- Time-reversed advection with curvature-based refinement
- Rotation of a straight interface
- Comparison between the explicit and implicit diffusion schemes
- Comparison between explicit and implicit diffusion schemes on concentration tracer
- Numerical viscosity for vorticity/streamfunction formulation
- Advection of a cosine bell around the sphere
1 0 GfsAdvection GfsBox GfsGEdge {} {
Time { end = 0.785398 }
Refine LEVEL
VariableTracer T { gradient = gfs_center_gradient }
Init {} {
T = {
double r2 = x*x + y*y;
double coeff = 20. + 20000.*r2*r2*r2*r2;
return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
}
}
VariableStreamFunction Psi -4.*(x*x + y*y)
OutputErrorNorm { start = end } { awk '{ print LEVEL " " $5 " " $7 " " $9}' } { v = T } {
s = {
double r2 = x*x + y*y;
double coeff = 20. + 20000.*r2*r2*r2*r2;
return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
}
}
OutputScalarSum { istep = 1 } t { v = T }
}
1 0 GfsAdvection GfsBox GfsGEdge {} {
Time { end = 5 }
Refine 8
# Add tracer T, using a VOF advection scheme.
# The default scheme is a Van-Leer limited, second-order upwind scheme.
VariableTracerVOFHeight T
InitFraction T (ellipse (0, -.236338, 0.2, 0.2))
# Maintain U and V with the vortical shear flow field defined by
# its streamfunction
VariableStreamFunction {
# make sure that the entire field is reinitialised at t = 2.5
step = 2.5
} Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
# Adapt the mesh dynamically so that at any time the maximum of the gradient
# of T is less than or equal to zero per cell length
AdaptGradient { istart = 1 istep = 1 } { cmax = 0 maxlevel = 8 } T
OutputPPM { start = 0 } { convert ppm:- t-0.eps } { v = T }
OutputPPM { start = 2.5 } { convert ppm:- t-2.5.eps } { v = T }
OutputPPM { start = 5 } { convert ppm:- t-5.eps } { v = T }
# Add a new variable
Variable Tref
# Initialize Tref with the initial shape
InitFraction { start = end } Tref (ellipse (0, -.236338, 0.2, 0.2))
# Output the norms of the difference between T and Tref, stores that into
# new variable DT
OutputErrorNorm { start = end } norms { v = T } {
s = Tref v = DT
}
OutputPPM { start = end } { convert ppm:- dt-5.eps } { v = DT }
OutputScalarSum { istep = 1 } { awk '{ print $3,$5-8.743441e-01 }' > t } { v = T }
}
1 0 GfsAdvection GfsBox GfsGEdge {} {
Time { end = 5 }
Refine LEVEL
VariableTracerVOFHeight T
VariableVOFConcentration C T
VariableTracer G
InitFraction T CIRCLE
Init {} {
C = GAUSSIAN*T
G = GAUSSIAN
}
VariableStreamFunction {
step = 2.5
} Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
AdaptGradient { istart = 1 istep = 1 } { cmax = 0 maxlevel = LEVEL } T
AdaptError { istart = 1 istep = 1 } { cmax = 1e-3 maxlevel = LEVEL } C
AdaptError { istart = 1 istep = 1 } { cmax = 1e-3 maxlevel = LEVEL } G
OutputScalarSum { istep = 1 } t { v = T }
OutputScalarSum { istep = 1 } t1 { v = C }
OutputScalarSum { istep = 1 } t2 { v = G }
# OutputSimulation { istep = 10 } stdout
OutputSimulation { start = 2.5 } half-LEVEL.gfs
EventScript { start = end } {
conservation()
{
if awk -v tolerance=$1 'BEGIN { min = 1.; max = -1.; }{
if ($5 > max) max = $5;
if ($5 < min) min = $5;
}END{ if (max - min > tolerance) { print max - min > "/dev/stderr"; exit (1); } }'; then
:
else
exit $GFS_STOP;
fi
}
conservation 1e-4 < t
conservation 2e-5 < t1
conservation 0. < t2
}
VariableTracer Tref
InitFraction { start = end } Tref CIRCLE
OutputErrorNorm { start = end } dt-LEVEL { v = T } {
s = Tref
v = DT
}
OutputErrorNorm { start = end } dt1-LEVEL { v = C } {
s = Tref*GAUSSIAN
v = DC
}
OutputErrorNorm { start = end } dt2-LEVEL { v = G } {
s = GAUSSIAN
v = DG
}
OutputSimulation { start = end } end-LEVEL.gfs
}
1 0 GfsAdvection GfsBox GfsGEdge {} {
Time { end = 5 }
Refine 8
# Add tracer T, using a VOF advection scheme.
# The default scheme is a Van-Leer limited, second-order upwind scheme.
VariableTracerVOFHeight T
# Curvature K of the interface defined by T
VariableCurvature K T
InitFraction T (ellipse (0, -.236338, 0.2, 0.2))
# Maintain U and V with the vortical shear flow field defined by
# its streamfunction
VariableStreamFunction {
# make sure that the entire field is reinitialised at t = 2.5
step = 2.5
} Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
# Adapt the mesh dynamically using a custom cost function returning
# the cell size times the local curvature (only for cells cut by
# the interface).
# The maximum cost is set to 0.1 i.e. the radius of curvature must be
# resolved with at least 10 cells.
AdaptFunction { istep = 1 } {
cmax = 0.1 maxlevel = 8
cfactor = 2.
} (T > 0. && T < 1.)*dL*fabs(K)
AdaptThickness { istep = 1 } { maxlevel = 8 } T
GModule gfsview
OutputView { start = 2.5 } t-2.5.eps { format = EPS line_width = 0.5 } curvature.gfv
# Add a new variable
Variable Tref
# Initialize Tref with the initial shape
InitFraction { start = end } Tref (ellipse (0, -.236338, 0.2, 0.2))
# Output the norms of the difference between T and Tref, stores that into
# new variable DT
OutputErrorNorm { start = end } norms { v = T } {
s = Tref v = DT
}
OutputPPM { start = end } { convert ppm:- dt-5.eps } { v = DT }
OutputScalarSum { istep = 1 } { awk '{ print $3,$5-0.8743665 }' > t } { v = T }
}
1 0 GfsAdvection GfsBox GfsGEdge {} {
PhysicalParams { L = 2. }
Refine 3
VariableTracerVOFHeight T
Init {} { U = y }
# fixme: 1e-9 to circumvent a bug in GfsView
InitFraction T (x - 1e-9)
Variable Tref
InitFraction { istep = 1 } Tref (x - t*y)
OutputErrorNorm { istep = 1 } error { v = T } { s = Tref }
Time { end = 5 }
GModule gfsview
OutputScalarSum { istep = 1 } t { v = T }
OutputView { step = 1 } rotate-%g.gnu { format = Gnuplot } rotate.gfv
OutputView { start = end } cells.gnu { format = Gnuplot } cells.gfv
EventScript { start = end } {
gnuplot <<EOF
set term postscript eps color lw 3 18
set output 'error.eps'
set ylabel 'Maximum volume fraction error'
set xlabel 'Time'
set key top left
set logscale y
set grid y
plot [][1e-3:]'error' u 3:9 w l t 'n=4', 'error.n1' u 3:9 w l t 'n=1'
EOF
for t in 0 1 2 5; do
gnuplot <<EOF
set term postscript eps color lw 2 solid
set output 'rotate-$t.eps'
unset border
unset xtics
unset ytics
unset key
set size ratio -1
plot 'cells.gnu' w l lc 0, 'rotate-$t.gnu' w l lw 2 lc 1, 'n1-$t.gnu' w l lc 2 lw 2
EOF
done
status=0;
if awk '{if ($5 != 2.) { print $5 > "/dev/stderr"; exit (1); }}' < t ; then :
else
status=$GFS_STOP;
fi
paste error error.ref > tmp
if awk '{
if ($9 > $20) {
print $3,$9,$20 > "/dev/stderr"
exit (1);
}
}' < tmp; then :
else
status=$GFS_STOP;
fi
exit $status
}
}
GfsBox {
left = Boundary {
BcDirichlet U y
BcAngle T (90. + atan2(1.,t)*180./M_PI)
}
right = Boundary {
BcDirichlet U y
BcAngle T (90. - atan2(1.,t)*180./M_PI)
}
top = Boundary {
BcDirichlet U y
BcAngle T (180. - atan2(1.,t)*180./M_PI)
}
bottom = Boundary {
BcDirichlet U y
BcAngle T atan2(1.,t)*180./M_PI
}
}
1 0 GfsAdvection GfsBox GfsGEdge {} {
VariableTracer T { scheme = none }
# SourceDiffusion T 1 { beta = 1 }
SourceDiffusionExplicit T 1
Solid (-sphere(0,0,0,0.4))
Refine 5
SurfaceBc T Dirichlet 1
# Just checks that it works with hypre
GModule hypre
Time { end = 0.02 dtmax = 9e-5 }
# OutputTime { istep = 1 } stderr
OutputSimulation { start = 0.01 } end.gfs
OutputSimulation { start = 2.5e-3 step = 2.5e-3 } {
awk '{
if ($1 != "#")
print sqrt($1*$1+$2*$2),$8;
}' > prof
} { format = text }
}
1 0 GfsAdvection GfsBox GfsGEdge {} {
VariableTracerVOF T
VariableTracer T1 { scheme = none }
VariableVOFConcentration C T
Refine 5
InitFraction T (-difference(sphere(0,0,0,0.4),sphere(0,0,0,0.1)))
Init {} {
C = GAUSSIAN*T
T1 = C
}
# GModule hypre
SourceDiffusion T1 T/(T+1.e9*(1.-T)*(1.-T)) {beta = 1}
SourceDiffusionExplicit C T/(T+1.e9*(1.-T)*(1.-T))
# SourceDiffusion C T/(T+1.e9*(1.-T)*(1.-T))
Time { end = 0.02 dtmax = 9e-5 }
OutputSimulation { start = 0.01 } end.gfs
OutputScalarSum {istep = 10} T1vol { v = T1 }
OutputScalarSum {istep = 10} Cvol { v = C }
EventScript { start = end } {
conservation()
{
if awk -v tolerance=$1 'BEGIN { min = 1.; max = -1.; }{
if ($5 > max) max = $5;
if ($5 < min) min = $5;
}END{ if (max - min > tolerance) { print max - min > "/dev/stderr"; exit (1); } }'; then
:
else
exit $GFS_STOP;
fi
}
conservation 2.1e-7 < T1vol
conservation 2.1e-7 < Cvol
}
OutputSimulation {start = 0.005} {
awk '{
if ($1 != "#" && $8 > 0.)
print sqrt($1*$1+$2*$2),$12,$13;
}' > prof
} { format = text}
}
1 2 GfsAdvection GfsBox GfsGEdge {} {
Time { end = 2 }
Refine LEVEL
VariableTracer Omega
Init {} {
U0 = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
V0 = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
# we need (U,V) to compute the vorticity
U = U0
V = V0
}
# In the general case the vorticity cannot be consistently
# defined directly (because it needs to verify discrete properties
# on the circulation)
# Initialising from the velocity field (U,V) as done below should be safe
Init {} { Omega = Vorticity }
VariablePoisson Psi -Omega
VariableStreamFunction { istep = 1 } Psi1 Psi
OutputScalarNorm { istep = 1 } divLEVEL { v = (dx("U") + dy("V")) }
OutputScalarSum { istep = 1 } kineticLEVEL { v = Velocity2 }
OutputScalarSum { istep = 1 } stdout { v = Velocity2 }
OutputErrorNorm { start = end } { awk '{ print $3,$5,$7,$9 }' > errorLEVEL.dat } {
v = Velocity
} {
s = sqrt(U0*U0+V0*V0)
v = E
relative = 1
}
}
6 12 GfsAdvection GfsBox GfsGEdge {} {
PhysicalParams { L = 2.*M_PI/4. }
MetricCubed M LEVEL
# Use alternative implementation
# MetricCubed1 M
# MapFunction {
# x = atan2 (X, Z)*180./M_PI
# y = asin (CLAMP(Y,-1.,1.))*180./M_PI
# }
Time { end = 1 }
Refine LEVEL
VariableTracer T {
gradient = gfs_center_gradient
cfl = 1
}
Global {
#define DTR (M_PI/180.)
double lambda1 (double lambda, double theta, double lambdap, double thetap) {
/* eq. (8) */
return atan2 (cos (theta)*sin (lambda - lambdap),
cos (theta)*sin (thetap)*cos (lambda - lambdap) -
cos (thetap)* sin (theta));
}
double theta1 (double lambda, double theta, double lambdap, double thetap) {
/* eq. (9) */
return asin (sin (theta)*sin (thetap) + cos (theta)*cos (thetap)*cos (lambda - lambdap));
}
double bell (double lambda, double theta, double lc, double tc) {
double h0 = 1.;
double R = 1./3.;
double r = acos (sin(tc)*sin(theta) + cos (tc)*cos (theta)*cos (lambda - lc));
return r >= R ? 0. : (h0/2.)*(1. + cos (M_PI*r/R));
}
double bell_moving (double lambda, double theta, double t) {
double lambda2 = lambda1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
double theta2 = theta1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
double lc = 3.*M_PI/2. - U0*t, tc = 0.;
return bell (lambda2, theta2, lc, tc);
}
}
Init {} { T = bell_moving(x*DTR,y*DTR,0) }
VariableStreamFunction Psi -U0*(sin (y*DTR)*cos (DTR*ALPHA) - cos (x*DTR)*cos (y*DTR)*sin (DTR*ALPHA))
AdaptGradient { istep = 1 } { cmax = 1e-4 maxlevel = LEVEL } T
OutputTime { istep = 10 } stderr
# OutputSimulation { istep = 10 } stdout
OutputSimulation { start = end } end-LEVEL-ALPHA.gfs
OutputErrorNorm { istep = 1 } { awk '{ print LEVEL,$3,$5,$7,$9}' > error-LEVEL-ALPHA } { v = T } {
s = bell_moving(x*DTR,y*DTR,t)
v = E
relative = 1
}
OutputScalarSum { istep = 1 } t-LEVEL-ALPHA { v = T }
OutputScalarSum { istep = 1 } area-LEVEL-ALPHA { v = 1 }
EventScript { start = end } {
( cat isolines.gfv
echo "Save isolines.gnu { format = Gnuplot }"
echo "Clear"
cat reference.gfv
echo "Save reference.gnu { format = Gnuplot }"
echo "Clear"
cat zero.gfv
echo "Save zero.gnu { format = Gnuplot }"
) | gfsview-batch2D end-LEVEL-ALPHA.gfs
cat <<EOF | gnuplot
set term postscript eps lw 2 18 color
set output 'isolines-LEVEL-ALPHA.eps'
set size ratio -1
set xlabel 'Longitude'
set ylabel 'Latitude'
unset key
plot [60:120][-30:30]'isolines.gnu' w l, 'reference.gnu' w l, 'zero.gnu' w l
EOF
fixbb isolines-LEVEL-ALPHA.eps
rm -f isolines.gnu reference.gnu zero.gnu
# check mass conservation
if awk '
BEGIN { min = 1000.; max = -1000.; }{
if ($5 < min) min = $5;
if ($5 > max) max = $5;
}
END {
if (max - min != 0.)
exit (1);
}' < t-LEVEL-ALPHA; then
exit 0
else
exit $GFS_STOP
fi
}
}
