GfsSimulation
From Gerris
Examples
- B\'enard--von K\'arm\'an Vortex Street for flow around a cylinder at Re=160
- Vortex street around a "heated" cylinder
- Rayleigh-Taylor instability
- Boussinesq flow generated by a heated cylinder
- Turbulent air flow around RV Tangaroa
- Coalescence of a pair of Gaussian vortices (Gerris logo)
- Air-water flow around a Series 60 cargo ship
- Estimation of the numerical viscosity
- Estimation of the numerical viscosity with refined box
- Convergence for a simple periodic problem
- Convergence for the three-way vortex merging problem
- Lid-driven cavity at Re=1000
- Lid-driven cavity at Re=1000 (explicit scheme)
- Poiseuille flow
- Creeping Couette flow of Generalised Newtonian fluids
- Momentum conservation for large density ratios
- Hydrostatic balance with solid boundaries and viscosity
- Hydrostatic balance with quadratic pressure profile
- Convergence of a potential flow solution
- Flow through a divergent channel
- Potential flow around a thin plate
- Circular droplet in equilibrium
- Planar capillary waves
- Air-Water capillary wave
- Fluids of different densities
- Pure gravity wave
- Shape oscillation of an inviscid droplet
8 7 GfsSimulation GfsBox GfsGEdge {} {
# Stop the simulation at t = 15
Time { end = 15 }
# Use an initial refinement of 6 levels (i.e. 2^6=64x64 for each box)
Refine 6
# Insert the solid boundary defined as x*x + y*y - 0.0625*0.0625 = 0
# (i.e. a cylinder of radius 0.0625 centered on the origin)
Solid (x*x + y*y - 0.0625*0.0625)
# Add a passive tracer called T
VariableTracer {} T
# Set the initial x-component of the velocity to 1
Init {} { U = 1 }
# Adapt the mesh using the vorticity criterion at every timestep
# down to a maximum level of 6 and with a maximum tolerance of 1e-2
AdaptVorticity { istep = 1 } { maxlevel = 6 cmax = 1e-2 }
# Adapt the mesh using the gradient criterion on variable T at
# every timestep, down to a maximum level of 6 and with a maximum tolerance of 1e-2
AdaptGradient { istep = 1 } { maxlevel = 6 cmax = 1e-2 } T
# Set a viscosity source term on the velocity vector with x-component U
# The Reynolds number is Re = D*U/Nu = 0.125*1/0.00078125 = 160
# where D is the cylinder diameter (as defined in cylinder.gts)
SourceDiffusion {} U 0.00078125
SourceDiffusion {} V 0.00078125
# Writes the time and timestep every 10 timesteps on standard error
OutputTime { istep = 10 } stderr
# Writes the simulation size every 10 timesteps on standard error
OutputBalance { istep = 10 } stderr
# Writes info about the convergence of the Poisson solver on standard error
OutputProjectionStats { istep = 10 } stderr
# Pipes a bitmap PPM image representation of the vorticity field at every other timestep
# into a conversion pipeline to create a MPEG movie called vort.mpg
# Sets the minimum used for colormapping to -10 and the maximum to 10
OutputPPM { istep = 2 } { ppm2mpeg > vort.mpg } {
min = -10 max = 10 v = Vorticity
}
# Pipes a bitmap PPM image representation of the T field at every other timestep
# into a MJPEGTools conversion pipeline to create a MPEG movie called t.mpg
# Sets the minimum used for colormapping to 0 and the maximum to 1
OutputPPM { istep = 2 } { ppm2mpeg > t.mpg } {
min = 0 max = 1 v = T
}
# Pipes a bitmap PPM image representation of the vorticity field at time 15
# into the ImageMagick converter "convert" to create the corresponding EPS file
OutputPPM { start = 15 } { convert -colors 256 ppm:- vort.eps } {
min = -10 max = 10 v = Vorticity
}
# Pipes a bitmap PPM image representation of the T field at time 15
# into the ImageMagick converter "convert" to create the corresponding EPS file
OutputPPM { start = 15 } { convert -colors 256 ppm:- t.eps } {
min = 0 max = 1 v = T
}
# Outputs profiling information at the end of the simulation to standard error
OutputTiming { start = end } stderr
}
8 7 GfsSimulation GfsBox GfsGEdge {} {
# Stop the simulation at t = 15
Time { end = 15 }
# Use an initial refinement of 6 levels (i.e. 2^6=64x64 for each box)
Refine 6
# Insert the solid boundary defined as x*x + y*y - 0.0625*0.0625 = 0
# (i.e. a cylinder of radius 0.0625 centered on the origin)
Solid (x*x + y*y - 0.0625*0.0625)
# Add a passive tracer called T
VariableTracer {} T
# Add diffusion to tracer T
SourceDiffusion {} T 0.001
# Dirichlet boundary condition for T on the cylinder
SurfaceBc T Dirichlet 1
# Set the initial x-component of the velocity to 1
Init {} { U = 1 }
# Adapt the mesh using the vorticity criterion at every timestep
# down to a maximum level of 6 and with a maximum tolerance of 1e-2
AdaptVorticity { istep = 1 } { maxlevel = 6 cmax = 1e-2 }
# Adapt the mesh using the gradient criterion on variable T at
# every timestep, down to a maximum level of 6 and with a maximum tolerance of 1e-2
AdaptGradient { istep = 1 } { maxlevel = 6 cmax = 1e-2 } T
# Set a viscosity source term on the velocity vector with x-component U
# The Reynolds number is Re = D*U/Nu = 0.125*1/0.00078125 = 160
# where D is the cylinder diameter (as defined in cylinder.gts)
SourceDiffusion {} U 0.00078125
SourceDiffusion {} V 0.00078125
# Writes the time and timestep every 10 timesteps on standard error
OutputTime { istep = 10 } stderr
# Writes the simulation size every 10 timesteps on standard error
OutputBalance { istep = 10 } stderr
# Writes info about the convergence of the Poisson solver on standard error
OutputProjectionStats { istep = 10 } stderr
# Pipes a bitmap PPM image representation of the T field at every other timestep
# into a conversion pipeline to create a MPEG movie called t.mpg
# Sets the minimum used for colormapping to 0 and the maximum to 0.4
OutputPPM { istep = 2 } { ppm2mpeg > t.mpg } {
min = 0 max = 0.4 v = T
}
# Pipes a bitmap PPM image representation of the T field at time 15
# into the ImageMagick converter "convert" to create the corresponding EPS file
OutputPPM { start = 15 } { convert -colors 256 ppm:- t.eps } {
min = 0 max = 0.4 v = T
}
# Outputs profiling information at the end of the simulation to standard error
OutputTiming { start = end } stderr
}
4 3 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1 dtmax = 5e-3 }
Refine 7
# The tracer T is used to track both phases
VariableTracerVOF {} T
# The initial sinusoidal interface (translated by 0.5 along the y-axis)
InitFraction {} T (0.05*cos (2.*M_PI*x) + y) { ty = 0.5 }
AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 2e-2 }
AdaptGradient { istep = 1 } { maxlevel = 7 cmax = 1e-2 } T
# The dynamic viscosity for both phases
SourceViscosity {} 0.00313
# This defines the inverse of the density of the fluids as a
# function of T
PhysicalParams { alpha = 1./(T*1.225 + (1. - T)*0.1694) }
# We also need gravity
Source {} V -9.81
OutputTime { istep = 10 } stderr
OutputBalance { istep = 10 } stderr
OutputProjectionStats { istep = 10 } stderr
OutputDiffusionStats { istep = 10 } stderr
OutputPPM { istep = 2 } { ppm2mpeg > vort.mpg} {
min = -30 max = 30 v = Vorticity
}
OutputPPM { istep = 2 } { ppm2mpeg > t.mpg } {
min = 0 max = 1 v = T
}
OutputPPM { start = end } { convert -colors 256 ppm:- vort.eps } {
min = -30 max = 30 v = Vorticity
}
OutputPPM { start = end } { convert -colors 256 ppm:- t.eps } {
min = 0 max = 1 v = T
}
OutputTiming { start = end } stderr
OutputSimulation { step = 0.1 } stdout
EventScript { start = 0 } { echo "Save t-0.eps { format = EPS }" }
EventScript { start = 0.7 step = 0.1 } { echo "Save t-$GfsTime.eps { format = EPS }" }
}
3 2 GfsSimulation GfsBox GfsGEdge {} {
# Limit the maximum timestep to 1e-2 so that the initial diffusion
# is properly resolved
Time { end = 20 dtmax = 1e-2 }
# Use an initial refinement of 8 levels around the solid boundary
RefineSolid 8
# Insert the solid boundary defined implicitly by the
# ellipse() function
Solid (ellipse(0.,-0.15,1./16.,1./16.))
# Add a passive tracer called T
VariableTracer T
# Add diffusion to tracer T
SourceDiffusion T 0.0001
# Add a source term to the vertical velocity component equal to T
Source V T
# Dirichlet boundary condition for T on the cylinder
SurfaceBc T Dirichlet 1
# Adapt the mesh using the vorticity criterion at every timestep
# down to a maximum level of 8 if y is smaller than 1.5, 0
# otherwise. The topmost part of the domain will not be refined and
# will act as a very efficient "sponge" layer to damp any eddies
# before they exit the domain.
AdaptVorticity { istep = 1 } { maxlevel = (y > 1.5 ? 0 : 8) cmax = 1e-2 }
# Also adapt according to the tracer gradient
AdaptGradient { istep = 1 } { maxlevel = 8 cmax = 5e-2 } T
# Writes the time and timestep every 10 timesteps on standard error
OutputTime { istep = 10 } stderr
# Writes the simulation size every 10 timesteps on standard error
OutputBalance { istep = 10 } stderr
# Writes info about the convergence of the Poisson solver on standard error
OutputProjectionStats { istep = 10 } stderr
# Outputs profiling information at the end of the simulation to standard error
OutputTiming { start = end } stderr
# Outputs the simulation every 4 timesteps
OutputSimulation { istep = 4 } stdout
# Every 4 timesteps, GfsView will read the following command, after having read
# the simulation file and will output a PPM screenshot on its standard output
EventScript { istep = 4 } { echo "Save stdout { width = 256 height = 512 }" }
# At t = 19, GfsView will create the PPM file used in the doc.
EventScript { start = 19 } { echo "Save t.ppm { width = 256 height = 512 }" }
# At the end of the simulation this file is converted to EPS.
EventScript { start = end } { convert -colors 256 t.ppm t.eps ; rm -f t.ppm }
}
2 1 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 2 }
# Insert the solid boundary defined explicitly by the
# triangulated surface contained in the GTS file tangaroa.gts
Solid tangaroa.gts
Refine 5
RefineSolid 9
Init {} { U = 1. }
# Adapt only in the first GfsBox.
# The coarse resolution of the second box acts as an efficient "sponge"
# layer to dampen any eddy before it exits the domain.
AdaptVorticity { istep = 1 } { maxlevel = (x < 0.5 ? 8 : 0) cmax = 1e-2 }
OutputSolidStats {} stderr
OutputTime { istep = 1 } stderr
OutputBalance { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
# Store in SU the integral over time of U
# At the end of the simulation SU/(Total integration time) = SU/1.
# is the mean velocity
EventSum { start = 1 istep = 1 } U SU
EventSum { start = 1 istep = 1 } V SV
EventSum { start = 1 istep = 1 } W SW
# Store in SU the integral over time of U^2 (i.e. the variance)
EventSum { start = 1 istep = 1 } U*U SU2
EventSum { start = 1 istep = 1 } V*V SV2
EventSum { start = 1 istep = 1 } W*W SW2
# Output simulation on standard output (to be read and displayed by GfsView)
OutputSimulation { istep = 4 } stdout
# Sends a command to GfsView to save a 1024x768 PPM image on standard output
EventScript { istep = 4 } { echo "Save stdout { width = 1024 height = 768 }" }
EventScript { start = 1.5 } { echo "Save sections.ppm { width = 1024 height = 768 }" }
EventScript { start = end } {
convert -colors 256 sections.ppm sections.eps ; rm -f sections.ppm
}
OutputSimulation { start = end } simulation-sum {
variables = SU,SV,SW,SU2,SV2,SW2
}
OutputTiming { start = end } stderr
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 4 }
Refine 6
# Take a large domain to minimise the influence of boundaries but
# refine only in a small central disk.
Refine (sqrt(x*x + y*y) < 0.0625 ? 12 : 6)
# Initialise a vorticity field given by two gaussian distributions
InitVorticity {} {
/* We use nested functions for simplicity (this will not work on MACOSX) */
double vortex (double xc, double yc, double r) {
double r2 = (x - xc)*(x - xc) + (y - yc)*(y - yc);
return 2.*M_PI*exp (- 2.*r2/(r*r));
}
double r = 0.01, theta = 30.*M_PI/180.;
return vortex (-r*sin(theta), r*cos(theta), 0.01) +
vortex (r*sin(theta), -r*cos(theta), 0.01);
}
AdaptVorticity { istep = 1 } { cmax = 1e-2 maxlevel = 12 }
OutputTime { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
OutputSimulation { istep = 10 } stdout
OutputPPM { istep = 2 } { ppm2mpeg > logo.mpg } {
v = Vorticity
min = -0.1348 max = 6.22219
# Only generate the movie in a small box centered on the origin
box = -0.025,-0.025,0.025,0.025
}
EventScript { start = end } {
echo "Save logo.ppm { width = 1024 height = 1024 }"
sleep 5 # to wait for GfsView to finish writing the image
convert -transparent "#0000FF" logo.ppm -geometry 156x156 logo.png
montage -background white -geometry +0+0 logo.png logo.eps
rm -f logo.ppm
}
}
3 2 GfsSimulation GfsBox GfsGEdge {} {
# The wave drag on the hull has strong starting transients,
# also the mean wave field takes a relatively long time to
# establish.
Time { end = 10 }
# Nine levels is enough to get good agreement with towing tank
# data. Adding more levels will reveal finer-scale wave patterns
# (but the runs will take even longer...)
Global {
#define LEVEL 9
#define FROUDE 0.316
#define RATIO (1.2/1000.)
#define VAR(T,min,max) (min + CLAMP(T,0,1)*(max - min))
}
# Translate the model to simulate only half the domain
Solid S60-scaled.gts { ty = 0.5 }
# Refine the hull to LEVEL
RefineSolid LEVEL
# Refine the water surface to four levels
RefineSurface { return 4; } (1e-4 - z)
VariableTracerVOF T
# For high-density ratios we cannot use the volume fraction field
# directly to define the density. We need a smoother version.
VariableFiltered T1 T 1
Init {} { U = FROUDE }
# Initialise the water surface at z = 1e-4
InitFraction T (1e-4 - z)
# air/water density ratio
PhysicalParams { alpha = 1./VAR(T1,RATIO,1.) }
# Use the reduced gravity approach
VariablePosition Z T z
# g = 3, g' = 3*(rho1 - rho2)
SourceTension T -3.*(1. - RATIO) Z
# Force the horizontal component of the velocity to relax to
# 'FROUDE' (= 0.316) in a band on inflow (x <= -0.375)
Source U (x > -0.375 ? 0 : 10.*(FROUDE - U))
# Adapt the mesh using the vorticity criterion but only in the
# water side (T > 0.)
# The 'cmax' value can be lowered (e.g. to 1e-2) to increase
# the accuracy with which the weaker far-field waves are resolved.
AdaptFunction { istep = 1 } {
cmax = 1e-2
# Faster coarsening than with the default cfactor of 4 reduces
# the size of the simulation
# cfactor = 2
# Coarse 'sponge' for x >= 1.5
maxlevel = (x < 1.5 ? LEVEL : 4)
minlevel = 4
} {
return (T > 0.)*fabs (Vorticity)*ftt_cell_size (cell)/FROUDE;
}
# Pressure (i.e. wave drag) force on the hull
OutputSolidForce { istart = 1 istep = 1 } f
OutputTime { istep = 1 } stderr
OutputBalance { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
OutputTiming { istep = 10 } stderr
# Generation of animations
OutputSimulation { istep = 5 end = 4 } stdout
EventScript { istep = 5 end = 4 } { echo "Save stdout { width = 1600 height = 1200 }" }
OutputSimulation { start = 1 step = 1 } sim-%g.gfs
# Compresses the saved simulation files
EventScript { start = 1 step = 1 } { gzip -f -q sim-*.gfs }
# Graphics
EventScript { start = 10 } {
echo "Save stdout { width = 1600 height = 1200 }" | \
gfsview-batch3D sim-10.gfs.gz closeup.gfv | \
convert -colors 256 ppm:- closeup.eps
echo "Save stdout { width = 1600 height = 1200 }" | \
gfsview-batch3D sim-10.gfs.gz front.gfv | \
convert -colors 256 ppm:- front.eps
echo "Save stdout { width = 800 height = 600 }" | \
gfsview-batch3D sim-10.gfs.gz comparison.gfv | \
convert -trim ppm:- comparison.ppm
# echo "Save stdout { width = 800 height = 600 }" | \
# gfsview-batch3D sim-10.gfs.gz tank-data.gfv | \
# convert -trim -flip ppm:- tank-data.png
convert tank-data.png tank-data.ppm
montage -geometry +0+0 -tile 1x2 tank-data.ppm comparison.ppm png:- | \
convert -colors 256 png:- comparison.eps
cat <<EOF | gnuplot
set term postscript eps lw 3 solid 20 colour
set output 'f.eps'
set xlabel 'Time'
set ylabel 'Force'
plot 'f' u 1:(\$2*2.) every 10 w l t 'Drag', 'f' every 10 u 1:(\$4*2.) w l t 'Lift'
EOF
}
}
1 2 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 2 }
Refine LEVEL
Init {} {
U = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
V = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
}
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
OutputScalarNorm { istep = 1 } divLEVEL { v = Divergence }
OutputScalarSum { istep = 1 } kineticLEVEL { v = Velocity2 }
OutputScalarSum { istep = 1 } stdout { v = Velocity2 }
}
1 2 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 2 }
Refine (x > 0.25 || x < -0.25 || y > 0.25 || y < -0.25 ? LEVEL : LEVEL + 1)
Init {} {
U = (- cos (8.*M_PI*x)*sin (8.*M_PI*y))
V = (sin (8.*M_PI*x)*cos (8.*M_PI*y))
}
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
OutputScalarNorm { istep = 1 } divLEVEL { v = Divergence }
OutputScalarSum { istep = 1 } kineticLEVEL { v = Velocity2 }
OutputScalarSum { istep = 1 } stdout { v = Velocity2 }
}
1 2 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 0.5 }
AdvectionParams { cfl = 0.75 }
Refine (x < -0.25 || x > 0.25 || y < -0.25 || y > 0.25 ? LEVEL : LEVEL + BOX)
Init {} {
U = (1. - 2.*cos (2.*M_PI*x)*sin (2.*M_PI*y))
V = (1. + 2.*sin (2.*M_PI*x)*cos (2.*M_PI*y))
}
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
OutputErrorNorm { start = end } stdout { v = U } {
s = (1. - 2.*cos (2.*M_PI*(x - t))*sin (2.*M_PI*(y - t)))
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 0.25 }
AdvectionParams { cfl = 0.9 }
ApproxProjectionParams { tolerance = 1e-5 }
ProjectionParams { tolerance = 1e-5 }
Refine {
double r = sqrt(x*x + y*y);
switch (LEVEL) {
case 6: return r > 0.25 ? 4 : r > 0.15 ? 5 : 6;
case 7: return r > 0.25 ? 4 : r > 0.2 ? 5 : r > 0.15 ? 6 : 7;
case 8: return r > 0.25 ? 4 : r > 0.2 ? 5 : r > 0.175 ? 6 : r > 0.15 ? 7 : 8;
case 9: return r > 0.25 ? 4 : r > 0.2 ? 5 : r > 0.175 ? 6 : r > 0.1625 ? 7 : r > 0.15 ? 8 : 9;
}
}
InitVorticity {} {
double vortex (double xo, double yo, double s) {
double r = sqrt ((x - xo)*(x - xo) + (y - yo)*(y - yo));
return s*(1. + tanh (100.*(0.03 - r)))/2.;
}
return vortex (0., 0., -150.) +
vortex (0.09, 0., 50.) +
vortex (-0.045, 0.0779422863406, 50.) +
vortex (-0.045, -0.0779422863406, 50.);
}
AdaptVorticity { istep = 1 } { maxlevel = LEVEL cmax = 4e-3 }
OutputSimulation { start = 0.05 } stdout
EventScript { start = 0.05 } {
echo Clear
cat levels.gfv
echo Save tm_0_05.eps { format = EPS line_width = 0.1 }
echo Clear
cat vorticity.gfv
echo Save tv_0_05.eps { format = EPS line_width = 0.1 }
}
OutputSimulation { start = 0.15 } stdout
EventScript { start = 0.15 } {
echo Clear
cat levels.gfv
echo Save tm_0_15.eps { format = EPS line_width = 0.1 }
echo Clear
cat vorticity.gfv
echo Save tv_0_15.eps { format = EPS line_width = 0.1 }
}
OutputSimulation { start = 0.25 } stdout
EventScript { start = 0.25 } {
echo Clear
cat levels.gfv
echo Save tm_0_25.eps { format = EPS line_width = 0.1 }
echo Clear
cat vorticity.gfv
echo Save tv_0_25.eps { format = EPS line_width = 0.1 }
}
OutputSimulation { start = 0.25 } SIM
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
# Stop the simulation at t = 300 if convergence has not been reached before
Time { end = 300 }
# Use an initial refinement of 6 levels (i.e. 2^6=64x64)
Refine 6
# Set a viscosity source term on the velocity vector with x-component U
# The Reynolds number is Re = L*U/Nu = 1*1/1e-3 = 1000
SourceDiffusion {} U 1e-3
SourceDiffusion {} V 1e-3
# Stops the simulation if the maximum of the absolute value of the
# difference between the current U field and the U field 10 timesteps
# before is smaller than 1e-4.
#
# Stores this difference in the DU field (this can be used for
# monitoring the convergence of the simulation).
EventStop { istep = 10 } U 1e-4 DU
OutputScalarNorm { istep = 10 } du { v = DU }
# Pipes a bitmap PPM image representation of the velocity field at the end of the simulation
# into the ImageMagick converter "convert" to create the
# corresponding EPS file
OutputPPM { start = end } { convert -colors 256 ppm:- velocity.eps } {
v = Velocity
}
# At the end of the simulation, computes the values of the variables
# at the locations defined in files xprofile, yprofile and stores the
# results in files xprof, yprof
OutputLocation { start = end } xprof xprofile
OutputLocation { start = end } yprof yprofile
# At the end of the simulation calls the script generating the EPS
# files using gnuplot and files: xprof, yprof, xprof.ghia, yprof.ghia
EventScript { start = end } {
cat <<EOF | gnuplot
set term postscript eps lw 3 solid 20
set output 'xprof.eps'
set xlabel 'Y'
set ylabel 'U'
plot [-0.5:0.5]'xprof.ghia' u 1:2 title "Ghia et al." w p ps 2 pt 9, 'xprof' u 3:7 w l title "Gerris"
set output 'yprof.eps'
set xlabel 'X'
set ylabel 'V'
plot [-0.5:0.5]'yprof.ghia' u 1:2 title "Ghia et al." w p ps 2 pt 9, 'yprof' u 2:8 w l title "Gerris"
EOF
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
# Stop the simulation at t = 300 if convergence has not been reached before
Time { end = 300 }
# Use an initial refinement of 6 levels (i.e. 2^6=64x64)
Refine 6
# Set a viscosity source term on the velocity vector with x-component U
# The Reynolds number is Re = L*U/Nu = 1*1/1e-3 = 1000
SourceViscosityExplicit 1e-3
# Stops the simulation if the maximum of the absolute value of the
# difference between the current U field and the U field 10 timesteps
# before is smaller than 1e-4.
#
# Stores this difference in the DU field (this can be used for
# monitoring the convergence of the simulation).
EventStop { istep = 10 } U 1e-4 DU
OutputScalarNorm { istep = 10 } du { v = DU }
# Pipes a bitmap PPM image representation of the velocity field at the end of the simulation
# into the ImageMagick converter "convert" to create the
# corresponding EPS file
OutputPPM { start = end } { convert -colors 256 ppm:- velocity.eps } {
v = Velocity
}
# At the end of the simulation, computes the values of the variables
# at the locations defined in files xprofile, yprofile and stores the
# results in files xprof, yprof
OutputLocation { start = end } xprof ../xprofile
OutputLocation { start = end } yprof ../yprofile
# At the end of the simulation calls the script generating the EPS
# files using gnuplot and files: xprof, yprof, xprof.ghia, yprof.ghia
EventScript { start = end } {
cat <<EOF | gnuplot
set term postscript eps lw 3 solid 20
set output 'xprof.eps'
set xlabel 'Y'
set ylabel 'U'
plot [-0.5:0.5]'../xprof.ghia' u 1:2 title "Ghia et al." w p ps 2 pt 9, 'xprof' u 3:7 w l title "Gerris"
set output 'yprof.eps'
set xlabel 'X'
set ylabel 'V'
plot [-0.5:0.5]'../yprof.ghia' u 1:2 title "Ghia et al." w p ps 2 pt 9, 'yprof' u 2:8 w l title "Gerris"
EOF
}
}
1 1 GfsSimulation GfsBox GfsGEdge {} {
Refine LEVEL
# use backward Euler to avoid Crank-Nicholson oscillations in time
SourceViscosity 1. { beta = 1 }
Source U 1
# we need this so that acceleration can be balanced by viscous stress
# and yes, fixme, this is a mess...
AdvectionParams { gc = 1 }
EventStop { istep = 1 } U 1e-6 DU
ProjectionParams { tolerance = 1e-6 }
ApproxProjectionParams { tolerance = 1e-6 }
OutputErrorNorm { start = end } { awk '{print LEVEL,$5,$7,$9}' } { v = U } {
s = 1./2.*(1./4 - y*y)
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { iend = 100 dtmax = 1e-2 }
Refine 6
Solid (ellipse (0,0,0.25,0.25)) { flip = 1 scale = 1.9999 }
Solid (ellipse (0,0,0.25,0.25))
ApproxProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
SourceViscosity {} {
double mu, ty, N, mumax = 1000.;
double m;
switch (MODEL) {
case 0: /* Newtonian */
mu = 1.; ty = 0.; N = 1.; break;
case 1: /* Power-law (shear-thinning) */
mu = 0.08; ty = 0.; N = 0.5; break;
case 2: /* Herschel-Bulkley */
mu = 0.0672; ty = 0.12; N = 0.5; break;
case 3: /* Bingham */
mu = 1.; ty = 10.; N = 1.; break;
}
if (D2 > 0.)
m = ty/(2.*D2) + mu*exp ((N - 1.)*log (D2));
else {
if (ty > 0. || N < 1.) m = mumax;
else m = N == 1. ? mu : 0.;
}
return MIN (m, mumax);
} {
# Crank-Nicholson does not converge for these cases, we need backward Euler
# (beta = 0.5 -> Crank-Nicholson, beta = 1 -> backward Euler)
beta = 1
}
GfsSurfaceBc U Dirichlet (x*x + y*y > 0.140625 ? 0. : - ay)
GfsSurfaceBc V Dirichlet (x*x + y*y > 0.140625 ? 0. : ax)
EventStop { istep = 1 } U 1e-4 DU
OutputScalarNorm { istep = 1 } du-MODEL { v = DU }
OutputLocation { start = end } { awk '{if ($1 != "#") print $2,$8;}' > prof-MODEL } profile
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 0.5 }
Global {
#define var(T,min,max) (CLAMP(T,0,1)*(max - min) + min)
#define rho(T) var(T, 0.001, 1.)
#define mu(T) var(T, 1e-6, 1e-3)
#define level 7
#define radius 0.05
}
Refine level
ProjectionParams { tolerance = 1e-6 }
ApproxProjectionParams { tolerance = 1e-6 }
VariableTracerVOF T
VariableFiltered T1 T 1
InitFraction T (- ellipse(-0.3,0,radius,radius))
Init {} { U = T }
PhysicalParams { alpha = 1./rho(T1) }
SourceViscosity mu(T1)
AdaptVorticity { istep = 1 } { cmax = 0.3 maxlevel = level }
AdaptGradient { istep = 1 } { cmax = 1e-3 maxlevel = level } T
OutputScalarSum { istep = 1 } k { v = Velocity2*rho(T1) }
OutputScalarSum { istep = 1 } t { v = T }
EventScript { start = end } {
gnuplot <<EOF
set term postscript eps color lw 3 solid 20
set output 'k.eps'
set xlabel 'Time'
set ylabel 'Kinetic energy'
set grid
plot 'k' u 3:5 w l t ''
EOF
if awk '{if ($5 > 7.2e-3) exit (1);}' < k ; then
return 0;
else
return $GFS_STOP;
fi
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Refine 3
Source V -1
SourceViscosity 1e-2
Solid (ellipse(0.,0.,0.24,0.24))
Time { iend = 10 }
ApproxProjectionParams { tolerance = 1e-12 }
ProjectionParams { tolerance = 1e-12 }
OutputScalarNorm { istep = 1 } v { v = V }
EventScript { start = end } {
if awk '{if ($9 > 1e-12) exit (1);}' < v ; then
exit 0;
else
exit $GFS_STOP;
fi
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Refine 3
# This test case only works for constant refinement
# Refine (x*x + y*y < 0.2*0.2 ? 4 : 3)
# Note: it is important to use 'cy' rather than 'y' in the formula
# below so that the hydrostatic density distribution is correct
# even for 'cut cells'
Init {} { rho = (cy + 0.5) }
Source V -rho
SourceViscosity 1e-2
Solid (ellipse(0.,0.,0.24,0.24))
Time { iend = 10 }
ApproxProjectionParams { tolerance = 1e-12 }
ProjectionParams { tolerance = 1e-12 }
OutputScalarNorm { istep = 1 } v { v = V }
# Checks that the pressure profile is close to the exact solution
OutputErrorNorm { istep = 1 } p { v = P } {
s = -(cy*cy/2. + 0.5*cy)
unbiased = 1
}
EventScript { start = end } {
if awk '{if ($9 > 1e-12) exit (1);}' < v ; then :
else
exit $GFS_STOP;
fi
if awk '{if ($9 > 1e-12) exit (1);}' < p ; then :
else
exit $GFS_STOP;
fi
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { iend = 0 end = 1 }
AdvectionParams { scheme = none }
ApproxProjectionParams { tolerance = 1e-6 }
Refine LEVEL
Solid (ellipse (0.25, 0.25, 0.1, 0.1))
Solid (ellipse (-0.25, 0.125, 0.15, 0.1))
Solid (ellipse (0., -0.25, 0.2, 0.1))
Init {} { U = 1 }
OutputSimulation { start = end } sim-LEVEL {
variables = U,V,P
}
}
4 3 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1 }
AdvectionParams { cfl = 0.9 }
ProjectionParams { tolerance = 1e-6 }
ApproxProjectionParams { tolerance = 1e-6 }
Refine LEVEL
Global {
double channel (double x) {
double y1 = 0.2/4.;
double y2 = 1e-6/4.;
return x <= -0.25 ? y1 :
x < 0.25 ? y2 + 0.5*(y1 - y2)*(1. + cos (2.*M_PI*(x + 0.25))) :
y2;
}
}
Solid (0.125 - channel (x) - y) { scale = 4 tx = 1.5 }
Solid (y + 0.125 - channel (x)) { scale = 4 tx = 1.5 }
Init {} { U = 1 }
OutputSimulation { start = end } sim-LEVEL {
variables = U,V,P
}
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { iend = 30 dtmax = 1e-2 }
Refine 5
RefineSolid 6
Solid (cube(0,0,0,0.5)) { sy = 0.06251 tx = 0.031249 ty = -0.015 }
AdvectionParams { scheme = none }
Init {} { U = 1 }
OutputScalarNorm { start = end } stdout { v = Velocity }
}
1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = TMAX }
Refine LEVEL
RefineSurface {return 10;} CIRCLE
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
VariableTracerVOF T
VariableCurvature K T
SourceTension T 1 K
SourceDiffusion U MU
SourceDiffusion V MU
InitFraction T CIRCLE
Init {} { Tref = T }
AdaptGradient { istep = 1 } { cmax = 1e-6 maxlevel = LEVEL } T
EventStop { istep = 10 } T DT
OutputSimulation { start = end } stdout { depth = LEVEL }
OutputScalarNorm { istep = 1 } {
awk '{ print MU*$3/(0.8*0.8), $9*sqrt(0.8) }' > La-LAPLACE-LEVEL
} { v = Velocity }
OutputScalarNorm { istep = 1 } {
awk '{ print MU*$3/(0.8*0.8), $5, $7, $9 }' > E-LAPLACE-LEVEL
} { v = (Tref - T) }
OutputScalarStats { istep = 1 } {
awk '{ print MU*$3/(0.8*0.8), $5, $7, $9, $11 }' > K-LAPLACE-LEVEL
} { v = (K - 2.50771) }
OutputScalarNorm { istep = 1 } {
awk '{ print MU*$3/(0.8*0.8), $5, $7, $9 }' > EK-LAPLACE-LEVEL
} { v = (T > 0 && T < 1 ? K - 2.5 : 0) }
}
3 5 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 2.2426211256 }
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
# Decrease the resolution linearly down to 3 levels close to the
# bottom and top boundaries
Refine floor(LEVEL + 1 - (LEVEL - 2)*fabs(y)/1.5)
VariableTracerVOF T
VariableCurvature K T
SourceTension T 1 K
VariablePosition Y T y
SourceDiffusion U 0.0182571749236
SourceDiffusion V 0.0182571749236
InitFraction T (y - 0.01*cos (2.*M_PI*x))
OutputScalarNorm { step = 3.04290519077e-3 } {
awk '{printf ("%g %g\n", $3*11.1366559937, $9); fflush(stdout); }' > wave-LEVEL
} { v = (T > 0. && T < 1. ? Y : 0.) }
}
3 5 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.58928694288774963184 }
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
Refine floor(LEVEL + 1 - (LEVEL - 2)*fabs(y)/1.5)
VariableTracerVOF T
VariableFiltered T1 T 1
VariableCurvature K T
SourceTension T 1 K
VariablePosition Y T y
Global {
#define VAR(T,min,max) (min + CLAMP(T,0,1)*(max - min))
#define RHO(T) VAR(T, 1.2/1000., 1.)
#define MU(T) VAR(T, 1.8e-5/1.003e-3, 1.)
}
PhysicalParams { alpha = 1./RHO(T1) }
SourceViscosity 0.0182571749236*MU(T1)
InitFraction T (y - 0.01*cos (2.*M_PI*x))
OutputScalarNorm { step = 0.00198785108553814829 } {
awk '{printf ("%g %g\n", $3*15.7402, $9); fflush(stdout); }' > wave-LEVEL
} { v = (T > 0. && T < 1. ? Y : 0.) }
}
3 5 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.66481717925811447992 }
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
Refine floor(LEVEL + 1 - (LEVEL - 2)*fabs(y)/1.5)
VariableTracerVOF T
VariableCurvature K T
SourceTension T 1 K
VariablePosition Y T y
SourceDiffusion U 0.0182571749236
SourceDiffusion V 0.0182571749236
PhysicalParams { alpha = 1./(T + 0.1*(1. - T)) }
InitFraction T (y - 0.01*cos (2.*M_PI*x))
OutputScalarNorm { step = .00225584983639310905 } {
awk '{printf ("%g %g\n", $3*15.016663878457, $9); fflush(stdout); }' > wave-LEVEL
} { v = (T > 0. && T < 1. ? Y : 0.) }
}
1 1 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.66481717925811447992 }
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
Refine LEVEL
VariableTracerVOF {} T
# Line below is for direct imposition of gravity acceleration
# Source {} V 50
# It is better to use a formulation where the first-order
# hydrostatic pressure is substracted away (in particular it
# prevents the generation of "hydrostatic spurious currents")
VariablePosition {} Y T y
# acceleration of gravity is 50, the equivalent "reduced pressure"
# is 50*(1. - 0.1) = 45
SourceTension {} T 45 Y
SourceDiffusion {} U 0.0182571749236
SourceDiffusion {} V 0.0182571749236
PhysicalParams { alpha = 1./(T + 0.1*(1. - T)) }
InitFraction {} T (y - 0.01*cos (2.*M_PI*x))
OutputScalarNorm { step = .00225584983639310905 } {
awk '{printf ("%g %g\n", $3*16.032448313657, $9); fflush(stdout); }' > wave-LEVEL
} { v = (T > 1e-6 && T < 1. - 1e-6 ? Y : 0.) }
}
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
# fixme: A small amount of viscosity seems to be necessary to
# obtain a non-increasing kinetic energy at high-resolution (8
# levels).
SourceViscosityExplicit 7e-5*RHO(T1)
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
}
}
