GfsAdvection
From Gerris
Examples
- Convergence of the Godunov advection scheme
- Time-reversed VOF advection in a shear flow
- Time-reversed advection with curvature-based refinement
- 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.
VariableTracerVOF 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 1e-2 per cell length
AdaptGradient { istep = 1 } { cmax = 1e-2 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 8
# Add tracer T, using a VOF advection scheme.
# The default scheme is a Van-Leer limited, second-order upwind scheme.
VariableTracerVOF 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 } {
return T > 0. && T < 1. ? ftt_cell_size (cell)*fabs (K) : 0.;
}
OutputSimulation { start = 2.5 } stdout
EventScript { start = 2.5 } { echo "Save t-2.5.eps { format = EPS line_width = 0.5 }" }
# 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 }
}
6 12 GfsAdvection GfsBox GfsGEdge {} {
PhysicalParams { L = 2.*M_PI/4. }
MetricCubed M LEVEL
Time { end = 1 }
Refine LEVEL
VariableTracer T {
gradient = gfs_center_gradient
cfl = 1
}
Global {
#define DTR (M_PI/180.)
double bell (double x, double y, double t) {
double h0 = 1.;
double R = 1./3.;
double lc = 3.*M_PI/2. + U0*t, tc = 0.;
x *= DTR; y *= DTR;
double r = acos (sin(tc)*sin(y) + cos (tc)*cos (y)*cos (x - lc));
return r >= R ? 0. : (h0/2.)*(1. + cos (M_PI*r/R));
}
}
Init {} { T = bell(x,y,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(x,y,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 [-120:-60][-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
}
}
