# GfsSimulation

(Difference between revisions)
 Revision as of 15:12, 10 October 2010Popinet (Talk | contribs)← Previous diff Current revisionGeordieMcBain (Talk | contribs) (described the two numbers beginning the syntax) Line 1: Line 1: - GfsSimulation solves the incompressible Euler equations by default. If viscosity is added (e.g. using [[GfsSourceViscosity]]) the equations become incompressible Navier-Stokes. A variable density (the default density is unity) can be defined using [[GfsPhysicalParams]]. + GfsSimulation solves the incompressible Euler equations by default, or, if [[GfsSourceViscosity]] is added, the incompressible Navier–Stokes equations. A variable density (the default density is unity) can be defined using [[GfsPhysicalParams]]. + + The GfsSimulation syntax begins with two numbers, being the number of [[GfsBox]]es containing the domain, and the number of connections between them via edges (or faces in three dimensions), either via adjacency or periodicity.

## Current revision

GfsSimulation solves the incompressible Euler equations by default, or, if GfsSourceViscosity is added, the incompressible Navier–Stokes equations. A variable density (the default density is unity) can be defined using GfsPhysicalParams.

The GfsSimulation syntax begins with two numbers, being the number of GfsBoxes containing the domain, and the number of connections between them via edges (or faces in three dimensions), either via adjacency or periodicity.

### Examples

• B\'enard--von K\'arm\'an Vortex Street for flow around a cylinder at Re=160
• ```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

}
```

• Vortex street around a "heated" cylinder
• ```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

}
```

• Parallel simulation on four processors
• ```8 7 GfsSimulation GfsBox GfsGEdge {} {

# Stop the simulation at t = 15
Time { end = 15 }

# 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)

# Use an initial refinement of 6 levels (i.e. 2^6=64x64 for each
# box) only around the solid boundary
RefineSolid 6

# 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
# 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)
SourceViscosity 0.00078125

# Balance the number of elements across parallel subdomains at every
# timestep if the imbalance is larger than 0.1 (i.e. 10% difference
# between the largest and smallest subdomains).
EventBalance { istep = 1 } 0.1

# Writes the time and timestep every 10 timesteps on standard error
OutputTime { istep = 10 } stderr

# Writes the time and simulation balance every timestep in 'balance'
OutputTime { istep = 1 } balance
OutputBalance { istep = 1 } balance

# Writes info about the convergence of the Poisson solver on standard error
OutputProjectionStats { istep = 10 } stderr

# Save MPEG movie using GfsView module
GModule gfsview
OutputView { step = 0.05 } {
ppm2mpeg -s 800x100 > pid.mpg
} { width = 1600 height = 200 } pid.gfv

# Outputs profiling information at the end of the simulation to standard error
OutputTiming { start = end } stderr

# Generate graphics
OutputSimulation { start = end } end.gfs
EventScript { start = end } {
echo "Save pid.eps { format = EPS width = 800 height = 100 line_width = 0.2 }" | \
gfsview-batch2D pid.gfv end.gfs
awk '{
if (\$1 == "step:")
t = \$4;
else if (\$1 == "domain")
print t, 100.*(\$9/\$3 - 1.), \$3, \$5, \$9;
}' < balance > balance1
cat <<EOF | gnuplot
set term postscript eps lw 3 solid 20 colour
set output 'balance.eps'
set xlabel 'Time'
set ylabel 'Number of elements per processor'
set key bottom right
plot 'balance1' u 1:3 w l t 'minimum', '' u 1:4 w l t 'average', '' u 1:5 w l t 'maximum'
EOF
}
}
```

• Rayleigh-Taylor instability
• ```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 }" }
}
```

• Boussinesq flow generated by a heated cylinder
• ```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 }

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 }
}
```

• Coalescence of a pair of Gaussian vortices (Gerris logo)
• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 4 }
Refine 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 minlevel = 6 }
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. We also need to make sure that box size is a multiple
# of 1./64. so that the PPM image size stays constant (ffmpeg
# crashes on variable image sizes).
condition = (Level < 6 ||
(x >= -3./128. && x <= 3./128. && y >= -3./128. && y <= 3./128.))
}
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
}
}
```

• Collapse of a column of grains
• ```1 0 GfsSimulation GfsBox GfsGEdge {
# shift origin of the domain
x = 0.5 y = 0.5
} {
PhysicalParams { L = LDOMAIN }

Time { end = 5.8 dtmax = 1e-2 }

# We need to tune the solver
ApproxProjectionParams { tolerance = 1e-4 }
ProjectionParams { tolerance = 1e-4 }

# VOF tracer and interface positions
VariableTracerVOF T
VariablePosition X T x
VariablePosition Y T y

# "internal" tracer
VariableTracerVOF Ti
InitFraction Ti ({
double diametre = 5e-3;
double r0 = 0.677/6.26;
double centre = x - 4.*diametre/r0;
double top = y - H0 + 9.*diametre/r0;
double side = x - R0 + 5.*diametre/r0;
return union (union (-top, -side), centre);
})

# mu(I) granular rheology
SourceViscosity {} {
double Eta = MUF;
if (P > 0. && D2 > 0.) {
double In = sqrt(2.)*D*D2/sqrt(P);
double muI = .32 + (.28)*In/(.4 + In);
double Etamin = sqrt(pow(D,3));
Eta = MAX((muI*P)/(sqrt(2.)*D2), Etamin);
Eta = MIN(Eta,10);
}
// Classic trick: use harmonic mean for the dynamic viscosity
return 1./(T/Eta + (1. - T)/MUF);
} {
beta = 1
tolerance = 1e-4
}

# Track a "band" around the interface to resolve surface gradients
# properly
cmax = 0
maxlevel = LEVEL
} T

# Use constant resolution inside the granular material
AdaptFunction { istep = 1 } {
cmax = 0
maxlevel = LEVEL
} T

# density
PhysicalParams { alpha = 1./RHO(T) }

# gravity
Source V -1

# initial conditions
Refine 6
InitFraction T (union(H0 - y, R0 - x))

OutputTime { istep = 10 } stderr
OutputProjectionStats { istep = 10 } stderr
OutputDiffusionStats { istep = 10 } stderr

EventSum { istep = 1 } U SU
EventSum { istep = 1 } V SV

# remove ejected droplets (just in case)
RemoveDroplets { istep = 1 } T -1

# stop when the acceleration of any cell full of granular material
# is less than 1e-2/0.1
# Init { istep = 1 } { UT = U*(T == 1.) }
# EventStop { start = 0.1 step = 0.1 } UT 1e-2 DU
# OutputScalarNorm { istep = 10 } du-LEVEL { v = DU }

# check mass conservation
OutputScalarSum { istep = 10 } t-LEVEL { v = T }

OutputSimulation { istep = 10 } stdout

# generate profiles
OutputSimulation { start = 0 } snapshot-%g.gfs
OutputSimulation { start = 1.65132 } snapshot-%g.gfs
OutputSimulation { start = 2.3769 } snapshot-%g.gfs
OutputSimulation { start = 3.10248 } snapshot-%g.gfs
OutputSimulation { start = 3.80304 } snapshot-%g.gfs
OutputSimulation { start = 5.70456 } snapshot-%g.gfs

# interface positions
OutputScalarNorm { istep = 10 } X-H0-LEVEL { v = (T > 0.1 ? X : G_MAXDOUBLE) }
OutputScalarNorm { istep = 10 } Y-H0-LEVEL { v = (T > 0.1 ? Y : G_MAXDOUBLE) }
# position of the "center" of the column
OutputScalarNorm { istep = 10 } Yc-H0-LEVEL { v = (x < 0.1 ? Y : G_MAXDOUBLE) }

# center of mass
OutputScalarSum { istep = 10 } {
awk '{
print \$3,\$5/(H0*R0);
fflush (stdout);
}' > xg-H0-LEVEL
} { v = T*x }
OutputScalarSum { istep = 10 } {
awk '{
print \$3,\$5/(H0*R0);
fflush (stdout);
}' > yg-H0-LEVEL
} { v = T*y }

# movie
GModule gfsview
OutputView { step = 5e-2 } { ppm2theora -s 640x240 > movie.ogv } {
width = 1280 height = 480
} movie.gfv

# generate figures
EventScript { start = end } {
for t in 0 1.65132 2.3769 3.10248 3.80304 5.70456; do
echo "Save snapshot-\$t.gnu { format = Gnuplot }" | \
gfsview-batch2D column.gfv snapshot-\$t.gfs
done
echo "Save movie.ppm { format = PPM width = 1280 height = 480 }" | \
gfsview-batch2D movie.gfv snapshot-0.gfs
convert movie.ppm -geometry 640x240 movie.png
convert movie.png movie.eps
rm -f movie.ppm
tar xzf grains.tgz
gnuplot comparison.plot
}
}
```

• Starting vortex of a NACA 2414 aerofoil
• ```3 2 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.5 }		# exit - trailing edge = (3-0.5) - 1 = 1.5
RefineSolid 9
Solid AEROFOIL.gts { rz = -INCIDENCE }
Init {} { U = 1 }
AdaptVorticity { istep = 1 } { maxlevel = 8 cmax = 5e-2 }

OutputTime { step = FRAMEPERIOD } stderr

OutputSimulation { step = FRAMEPERIOD } stdout

GModule gfsview
OutputView { step = FRAMEPERIOD } { ppm2mpeg > starting.mpg } {
width = 720 height = 240
} starting.gfv
OutputView { start = 1. } { convert ppm:- -geometry 720x240 starting.eps } {
format = PPM width = 1440 height = 480
} starting.gfv
}
```

• Viscous folding of a fluid interface
• ```4 3 GfsSimulation GfsBox GfsGEdge {} {
Global {
// Parameters
static double Re = 50;
static double ReSr = 1.2;
static double cavity_depth = 1.0;
static double floor_depth = 0.5;

// Solver parameters
static int min_level = 4;
static int max_level = 8;

// Density functions
static double density(double tracer) {
return 1.0;
}

// Velocity distributions
static double velocity_bc(double y, double t) {
if (y > floor_depth) return (0.125*(2*y - 1)*(2*y - 1)*(cos(2*M_PI*t*(ReSr / Re)) + 1));
else return 0;
}

// Initial tracer distribution
static double tracer_init(double y) {
if (y <= floor_depth) return 1.0;
else return 0.0;
}
}

# Define the depth of the cavity by filling the region below it
# with solid. The cavity will be filled with tracer on
# initialisation. Gerris will remove the cells which are
# completely solid.
Solid ({
return y - (floor_depth - cavity_depth);
})

# Specify simulation times
Time { end = 250 }

# The viscosity is nu = 1/Re.
SourceViscosity (1/Re)

# Initialise the tracer location
VariableTracerVOF T
Init {} {
T = tracer_init(y)
U = velocity_bc(y, 0)
V = 0
}

# We want adaptivity around the fluid interface
AdaptFunction { istep = 1 } {
cmax = 0
minlevel = min_level
maxlevel = max_level
} (T > 0 && T < 1)

# We also need to control the error on the velocity field
AdaptError { istep = 1 } { cmax = 1e-2 maxlevel = max_level } U
AdaptError { istep = 1 } { cmax = 1e-2 maxlevel = max_level } V

# Balance the number of elements across parallel subdomains at
# every timestep if the imbalance is larger than 0.1 (i.e. 10%
# difference between the largest and smallest subdomains).
EventBalance { istep = 1 } 0.1

# Output of solution information/data
OutputTime { istep = 10 } stderr
OutputProjectionStats { istep = 10 } stderr
OutputDiffusionStats { istep = 10 } stderr
OutputTiming { start = end } stderr

# use gfsview to generate movies and figures
GModule gfsview
OutputView { step = 0.5 } { ppm2mpeg -s 400x300 > viscmix.mpg } {
width = 1280 height = 960
} viscmix.gfv
OutputView { start = end } { convert -colors 256 ppm:- viscmix.eps } {
width = 1280 height = 960
} viscmix.gfv
}
```

• Turbulent air flow around RV Tangaroa
• ```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
}
```

• Savart--Plateau--Rayleigh instability of a water column
• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.7 }
Refine 5

VariableTracerVOF T
# Filter the volume fraction for smoother density transition for
# high density ratios: this helps the Poisson solver
VariableFiltered T1 T 1
PhysicalParams { alpha = 1./RHO(T1) }

# We need Kmax as well as K(mean) for adaptivity
VariableCurvature K T Kmax
SourceTension T 1 K
VariablePosition Y T y
VariablePosition Z T z
#    SourceViscosityExplicit 1e-2*MUR(T1)
SourceViscosity 1e-2*MUR(T1)

# Initial deformed tube (only a quarter of it)
InitFraction {} T ({
x -= 0.5;
y += 0.5; z += 0.5;
double r = RADIUS*(1. + EPSILON*cos(M_PI*x));
return r*r - y*y - z*z;
})

# Adapt according to interface curvature. Make sure we have at
# least 5 grid points per radius of curvature, up to level
# 10. Only start after 5 timesteps to ignore transients due to
# initialisation
AdaptFunction { istart = 5 istep = 10 } {
cmax = 0.2
maxlevel = 10
cfactor = 2
} (T > 0 && T < 1 ? dL*Kmax : 0)

# Remove small satellite droplets (only keep the two largest pieces)
RemoveDroplets { istep = 10 } T -2

OutputTime { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
OutputSimulation { istep = 1000 } plateau-%ld.gfs
EventScript { istep = 1000 } { gzip -f -q plateau-*.gfs }

# Generate three movies on the fly
GModule gfsview
OutputView { istep = 7 } { ppm2theora -s 480x480 > plateau.ogv } {
width = 960 height = 960
} plateau.gfv
OutputView { istep = 7 } { ppm2theora -s 480x480 > closeup.ogv } {
width = 960 height = 960
} closeup.gfv
OutputView { istep = 7 } { ppm2theora -s 480x480 > white.ogv } {
width = 960 height = 960
} white.gfv

# Generate figures
EventScript { start = end } {
for f in white plateau closeup; do
echo "Save \$f.ppm { format = PPM width = 960 height = 960 }" | \
gfsview-batch3D \$f.gfv plateau-1000.gfs.gz
convert \$f.ppm -geometry 480x480 \$f.png
convert \$f.png \$f.eps
rm -f \$f.ppm
done
cat <<EOF | gnuplot
set term postscript eps color lw 2 20
set output 'size.eps'
set xlabel 'Timestep'
set ylabel 'Total number of cells'
unset key
plot [10:]'< grep domain size' u 3 w l
EOF
}

OutputTiming { istep = 100 } stderr

OutputScalarNorm { istep = 1 } v { v = Velocity }

# Evolution of the minimum and maximum interface radii
OutputScalarStats { istep = 1 } r {
v = (T > 1e-2 && T < 1. - 1e-2 ?
(sqrt((Y + 0.5)*(Y + 0.5) + (Z + 0.5)*(Z + 0.5))/RADIUS - 1.)/EPSILON : 0)
}
OutputScalarStats { istep = 1 } k { v = K }
OutputScalarStats { istep = 1 } kmax { v = Kmax }
OutputScalarSum { istep = 1 } t { v = T }
OutputBalance { istep = 1 } size
}
```

• Atomisation of a pulsed liquid jet
• ```3 2 GfsSimulation GfsBox GfsGEdge {} {
Global {
#define length 0.025
#define level 9
#define Re 5800
#define R2(y,z) ((y)*(y) + (z)*(z))
#define rho(T) (T + 1./27.84*(1. - T))
/* Weber = rhoV^2D/sigma = 5555 */
}
Time { end = 1.6 }
# Initial refinement of the inlet
Refine (x < -0.5 + length && R2(y,z) < 2.*radius*radius ? level : 5)

# Define a static field used to enforce the boundary conditions
# for volume fraction corresponding to a cylindrical jet
Variable T0

VariableTracerVOF T
VariableCurvature K T Kmax
SourceTension T 0.00003 K
PhysicalParams { alpha = 1./rho(T) }

# Use constant (maximum) resolution on the interface
AdaptFunction { istep = 1 } {
minlevel = 0
maxlevel = level
} (T > 0 && T < 1)

# Initialise a short jet
Init {} {
T = (x < -0.5 + length ? T0 : 0)
U = T
}

EventBalance { istep = 1 } 0.1

OutputTime { istep = 1 } log
OutputBalance { istep = 1 } log
OutputProjectionStats { istep = 1 } log
OutputTiming { istep = 100 } log

# Use the gfsview module to generate movies on-the-fly and in parallel
GModule gfsview
OutputView { step = 4e-3 } { ppm2theora -s 640x480 > jet.ogv } {
format = PPM width = 1280 height = 960
} jet.gfv
OutputView { step = 4e-3 } { ppm2theora -s 640x480 > back.ogv } {
format = PPM width = 1280 height = 960
} back.gfv

# Save a (single) snapshot every 100 timesteps
EventScript { istep = 100 } { rm -f snapshot-*.gfs }
OutputSimulation { istep = 100 } snapshot-%ld.gfs { }

# Generate figures
OutputView { start = end } jet.ppm { format = PPM width = 1280 height = 960 } jet.gfv
OutputView { start = end } back.ppm { format = PPM width = 1280 height = 960 } back.gfv
EventScript { start = end } {
for f in jet back; do
convert \$f.ppm -geometry 640x480 \$f.png
convert \$f.png \$f.eps
rm -f \$f.ppm
done
awk '{if (\$1 == "step:") t = \$4;
else if (\$1 == "domain") print t,\$3,\$5,\$9;}' < log > balance
cat <<EOF | gnuplot
set term postscript eps color lw 2 20 solid
set output 'balance.eps'
set xlabel 'Time'
set ylabel 'Number of cells per processor'
plot [0:1.6]'balance' u 1:2 w l t 'Minimum', \
'balance' u 1:3 w l t 'Average', \
'balance' u 1:4 w l t 'Maximum'
EOF
}
}
```

• Air-water flow around a Series 60 cargo ship
• ```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
}
}
```

• Forced isotropic turbulence in a triply-periodic box
• ```1 3 GfsSimulation GfsBox GfsGEdge {} {

Time { end = 300 }

Refine 4
PhysicalParams{ L = 2.*M_PI }

# The initial condition is "ABC" flow. This is a laminar base flow that
# is easy to implement in both Gerris and a spectral code.
Init {} {
U = cos(y) + sin(z)
V = sin(x) + cos(z)
W = cos(x) + sin(y)
}

# Set the viscosity mu
SourceViscosity MU

# Calculate the mean flow
SpatialSum { istep = 1 } SU U
SpatialSum { istep = 1 } SV V
SpatialSum { istep = 1 } SW W
SpatialSum { istep = 1 } Un Velocity

# Adapt according to the relative error on the velocity field (with
# a 5% threshold)
AdaptError { istep = 1 } { cmax = 5e-2 maxlevel = 7 } U/Unbar
AdaptError { istep = 1 } { cmax = 5e-2 maxlevel = 7 } V/Unbar
AdaptError { istep = 1 } { cmax = 5e-2 maxlevel = 7 } W/Unbar

# Calculate -eps = mu * sum_{ij} (partial_j u_i)^2
SpatialSum { istep = 1 } Dissipation {
return MU*(dx("U")*dx("U") + dy("U")*dy("U") + dz("U")*dz("U") +
dx("V")*dx("V") + dy("V")*dy("V") + dz("V")*dz("V") +
dx("W")*dx("W") + dy("W")*dy("W") + dz("W")*dz("W"));
}

# The mean fluctuating kinetic energy
SpatialSum { istep = 1 } FluctKinEn {
return 0.5*((U - Ubar)*(U - Ubar) +
(V - Vbar)*(V - Vbar) +
(W - Wbar)*(W - Wbar));
}

# Add the linear forcing, subtracting the mean
Source U 0.1*(U - Ubar)
Source V 0.1*(V - Vbar)
Source W 0.1*(W - Wbar)

# Output
OutputTime { istep = 1 } log
OutputBalance { istep = 1 } log
OutputScalarStats { istep = 1 } log { v = Unbar }
OutputScalarStats { istep = 1 } log { v = U }
OutputScalarStats { istep = 1 } Reynolds.dat {
v = 2./3.*FluctKinEn/VOLUME/MU*sqrt(15*MU/(Dissipation/VOLUME))
}
OutputScalarStats { istep = 1 } Dissipation.dat { v = Dissipation/VOLUME }
OutputScalarStats { istep = 1 } Energy.dat { v = FluctKinEn/VOLUME }
OutputScalarStats { istep = 1 } Vorticity.dat { v = Vorticity }
EventScript { istep = 100 } { rm -f snapshot-*.gfs }
OutputSimulation { istep = 100 } snapshot-%ld.gfs
OutputSimulation { start = end } end.gfs

# Generate graphics
GModule gfsview
OutputView { step = 0.1 end = 150 } { ppm2mpeg > multiview.mpg } {
width = 512 height = 512
} multiview.gfv
OutputView { start = end } { convert ppm:- multiview.eps } {
width = 512 height = 512
} multiview.gfv

EventScript { start = end } {
gnuplot <<EOF
set term postscript eps lw 3 solid 20 colour

set output 'Energy.eps'
set xrange [0:300]
set xlabel 'Time'
set ylabel 'Kinetic energy'
set logscale y
plot 'Energy.dat' u 3:7 w l t 'Gerris', 'spectral.dat' u 1:(\\$3*3/2) w l t 'Spectral'

set output 'Reynolds.eps'
set ylabel 'Microscale Reynolds number'
plot 'Reynolds.dat' u 3:7 w l t 'Gerris', 'spectral.dat' u 1:4 w l t 'Spectral'

set output 'Dissipation.eps'
set ylabel 'Dissipation function'
plot 'Dissipation.dat' u 3:7 w l t 'Gerris', 'spectral.dat' u 1:2 w l t 'Spectral'

set output 'size.eps'
set ylabel 'Total number of grid points'
unset logscale
plot "< awk '{ if (\\$1 == \"step:\") t = \\$4; else if (\\$1 == \"domain\") print t,\\$5*8.;}' < log" w l t 'Gerris', 128**3 t 'spectral'
EOF
}
}
```

• Wingtip vortices behind a rectangular NACA 2414 wing
• ```6 7 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.5 }		# exit - trailing edge = (3-0.5) - 1 = 1.5
RefineSolid 7
Solid AEROFOIL.gts { rz = -INCIDENCE }
Init {} { U = 1 }
AdaptVorticity { istep = 1 } { maxlevel = 6 cmax = 5e-2 }

OutputTime { step = FRAMEPERIOD } stderr
OutputSimulation { step = FRAMEPERIOD } stdout

GModule gfsview
OutputView { step = FRAMEPERIOD } { ppm2mpeg > wingtip.mpg } wingtip.gfv
OutputView { start = 1 } { convert ppm:- wingtip.eps } wingtip.gfv
}
```

• Conservation of diffusive tracer
• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { iend = 100 dtmax = 2e-1 }
Refine 3
VariableTracer T
VariableTracer Te
InitFraction T (0.01 - (x*x + y*y + z*z))
InitFraction Te (0.01 - (x*x + y*y + z*z))
SourceDiffusion T 1e-4 { beta = 0.5 }
SourceDiffusionExplicit Te 1e-4
AdaptGradient { istep = 1 } { minlevel = 3 maxlevel = 5 cmax = 1e-2 } T
OutputScalarSum { istep = 1 } st { v = T }
OutputScalarSum { istep = 1 } ste { v = T }
OutputScalarStats { istep = 1 } t { v = T }
OutputScalarStats { istep = 1 } te { v = T }
OutputScalarNorm { istep = 1 } diff { v = (T - Te) }
EventScript { start = end } {
if awk '{if (\$9 > 8e-3) {
print "diff: " \$9 > "/dev/stderr"; exit (1);
}}' < diff &&
awk 'BEGIN{ s=-1 } {
if (s < 0.) s = \$5;
else if (\$5 - s != 0.) {
print "st: " \$5 - s > "/dev/stderr"; exit (1);
}
}' < st &&
awk 'BEGIN{ s=-1 } {
if (s < 0.) s = \$5;
else if (\$5 - s != 0.) {
print "ste: " \$5 - s > "/dev/stderr"; exit (1);
}
}' < ste ; then :
else
exit \$GFS_STOP;
fi
}
}
GfsBox{}
```

• Estimation of the numerical viscosity
• ```1 2 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 2 }
Refine LEVEL
Init {} {
U0 = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
V0 = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
U = U0
V = V0
}
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 }
OutputErrorNorm { start = end } { awk '{ print \$3,\$5,\$7,\$9 }' > errorLEVEL.dat } {
v = Velocity
} {
s = sqrt(U0*U0+V0*V0)
v = E
relative = 1
}
}
```

• Estimation of the numerical viscosity with refined box
• ```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 {} {
U0 = (- cos (8.*M_PI*x)*sin (8.*M_PI*y))
V0 = (sin (8.*M_PI*x)*cos (8.*M_PI*y))
U = U0
V = V0
}
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 }
OutputErrorNorm { start = end } { awk '{ print \$3,\$5,\$7,\$9 }' > errorLEVEL.dat } {
v = Velocity
} {
s = sqrt(U0*U0+V0*V0)
v = E
relative = 1
}
}
```

• Convergence for a simple periodic problem
• ```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)))
}
}
```

• Convergence for the three-way vortex merging problem
• ```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
}
```

• Flow created by a cylindrical volume source
• ```1 0 GfsSimulation GfsBox GfsGEdge { } {

Global {
#define R0 (0.05)
double sol (double x, double y) {
double r = sqrt(x*x+y*y);
return r >= R0 ? R0*R0*0.5/r : r*0.5;
}
}

Time { iend = 10 }

Refine (LEVEL - 4*pow((x*x+y*y)/0.25, 0.5))

Variable F
InitFraction F (R0*R0 - x*x - y*y)
Source P F
PhysicalParams { L = 2 }

OutputTime { istep = 1 } stderr

OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' > error-LEVEL } {
v = Velocity
} {
s = sol(x,y)
v = E
w = (x*x + y*y < 25.*R0*R0)
relative = 1
}

OutputSimulation { start = end } end-LEVEL.gfs
}
```

• Lid-driven cavity at Re=1000
• ```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

OutputSimulation { start = end } end.gfs

# 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 } {
gnuplot <<EOF
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
}
}
```

• Lid-driven cavity at Re=1000 (explicit scheme)
• ```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
}
}
```

• Lid-driven cavity with a non-uniform metric
• ```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

# Use a non-uniformly-stretched metric in the x- and y-directions
Metric M {
x = tanh(4.*rx)/tanh(4./2.)/2.
y = tanh(4.*ry)/tanh(4./2.)/2.
}

# Use hypre to accelerate convergence. It also works with the native
# solver but one needs to be careful about the tolerance for the
# projection
GModule hypre
ProjectionParams { tolerance = 1e-4 }
ApproxProjectionParams { tolerance = 1e-4 }

# 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 }

# Use gfsview module to generate a Gnuplot file with the mesh and isolines
# We use gnuplot because gfsview cannot (yet) take the metric into
# account when displaying results
GModule gfsview
OutputView { start = end } isolines.gnu { format = Gnuplot } isolines.gfv

# 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

OutputSimulation { start = end } end.gfs

# 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 } {
gnuplot <<EOF
set term postscript eps
set output 'velocity.eps'
set size ratio -1
unset border
unset key
unset xtics
unset ytics
plot 'isolines.gnu' w l
EOF
gnuplot <<EOF
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
}
}
```

• Lid-driven cavity on an anisotropic mesh
• ```2 1 GfsSimulation GfsBox GfsGEdge {
# we need to shift the origin of the reference box to (0,0.5)
y = 0.5
} {

# 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

# The mesh is stretched by a factor 0.5 (compressed) in the y direction.
MetricStretch {} { sy = 0.5 }

# 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:- -resize 128x128! 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

OutputSimulation { start = end } end.gfs

# 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 } {
gnuplot <<EOF
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
}
}
```

• Poiseuille flow
• ```1 1 GfsSimulation GfsBox GfsGEdge {} {
Refine LEVEL
# use backward Euler to avoid Crank-Nicholson oscillations in time
SourceViscosity 1. { beta = 1 }
Source U 1
# add a transverse source term to also test hydrostatic balance
Source V 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)
}
}
```

• Bagnold flow of a granular material
• ```1 1 GfsSimulation GfsBox GfsGEdge { y = 0.5 } {
Time { iend = 1000 dtmax = 5 }

# Ignore inertia
scheme = none
}
ApproxProjectionParams { tolerance = 1e-8 }
ProjectionParams { tolerance = 1e-8 }

# mu(I) granular rheology
SourceViscosity {} {
double In = sqrt(2.)*D*D2/sqrt(fabs(P));
double muI = .38 + (.26)*In/(.3 + In);
return MAX((muI*P)/(sqrt(2.)*D2), 0.);
} { beta = 1 }

# gravity
Source V -cos(ALPHA)
Source U  sin(ALPHA)

Refine LEVEL
Init {} {
# Bagnold's solution
U = 2.*UB(y)
P = (1. - y)*cos(ALPHA)
}

OutputTime { istep = 10 } stderr
OutputProjectionStats { istep = 10 } stderr
OutputDiffusionStats { istep = 10 } stderr

EventStop { istart = 10 istep = 10 } U 5e-6 DU
OutputScalarNorm { istep = 10 } du-LEVEL { v = DU }

OutputSimulation { start = end } end-LEVEL.txt { format = text }
OutputErrorNorm { start = end } { awk '{print LEVEL,\$5,\$7,\$9}' } { v = U } {
# Bagnold's solution
s = UB(y)
}
}
```

• Poiseuille flow with metric
• ```1 1 GfsSimulation GfsBox GfsGEdge {} {
# use a non-uniformly-stretched metric in the y-direction
# the grid points are denser close to the wall
Metric M {
y = tanh(4.*ry)/tanh(4./2.)/2.
}

Refine LEVEL
# use backward Euler to avoid Crank-Nicholson oscillations in time
SourceViscosity 1. { beta = 1 }
Source U 1
# add a transverse source term to also test hydrostatic balance
Source V 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)
}
}
```

• Creeping Couette flow of Generalised Newtonian fluids
• ```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
}
SurfaceBc U Dirichlet (x*x + y*y > 0.140625 ? 0. : - ay)
SurfaceBc 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
}
```

• Momentum conservation for large density ratios
• ```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
}

Refine level

ProjectionParams { tolerance = 1e-6 }
ApproxProjectionParams { tolerance = 1e-6 }

VariableTracerVOFHeight T
VariableFiltered T1 T 1
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
}
}
```

• Hydrostatic balance with solid boundaries and viscosity
• ```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 > 1.5e-12) { print \$9 > "/dev/stderr"; exit (1); }}' < v ; then
exit 0;
else
exit \$GFS_STOP;
fi
}
}
GfsBox {
bottom = Boundary
}
```

• Hydrostatic balance with quadratic pressure profile
• ```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 > 1.5e-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
}
}
```

• Coriolis formulation in 3-D
• ```1 3 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 30000  }

Refine 2

# The reference length is set to 60 m
PhysicalParams { L = 60 }

# The domain is 60 x 30 x 15 m with an anisotropic mesh
MetricStretch {} { sy = 0.5 sz = 0.25 }

# Initialisation of the velocity field
Init {} {
U = U0
V = V0
W = W0
}

# Viscosity (or eddy-viscosity)
SourceViscosity 10

# Coriolis forces
SourceCoriolis (2.*Omega) 0. { x = 0. y = 1. z = 1. }

# Output of the integral of each component of the velocity field over the domain every
# 1000 time steps
OutputScalarStats { step = 3000 } { awk '{print \$3, \$7}' > u.dat } { v = U }
OutputScalarStats { step = 3000 } { awk '{print \$3, \$7}' > v.dat } { v = V }
OutputScalarStats { step = 3000 } { awk '{print \$3, \$7}' > w.dat } { v = W }

# Output of some kind of error measurement
OutputScalarStats { istep = 10 } { awk '{print \$3,\$7}' > error.dat } {
v = sqrt((U - Usol0)*(U - Usol0) + (V - Vsol0)*(V - Vsol0) + (W - Wsol0)*(W - Wsol0))
}

# Output of the final simulation file
OutputSimulation { start = end } end.gfs
}
```

• Wind-driven lake
• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
GModule hypre
MetricStretch {} { sy = 0.1 }
Time { end = 5 dtmax = 0.1 }
SourceViscosity 1./400.
Refine 6
OutputTime { istep = 1 } stderr
OutputProjectionStats { istep = 1 } stderr
OutputDiffusionStats { istep = 1 } stderr
OutputScalarStats { istep = 1 } p { v = P }
GModule gfsview
OutputView { start = end } lake.eps { format = EPS } lake.gfv
OutputSimulation { start = end } end.gfs
}
```

• Convergence of a potential flow solution
• ```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
}
}
```

• Flow through a divergent channel
• ```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
}
}
```

• Potential flow around a thin plate
• ```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 }
}
```

• Circular droplet in equilibrium
• ```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 }

VariableTracerVOFHeight T
VariableCurvature K T
SourceTension T 1 K
SourceViscosity 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) }
}
```

• Planar capillary waves
• ```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)
VariableTracerVOFHeight 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.) }
}
```

• Air-Water capillary wave
• ```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)
VariableTracerVOFHeight 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.) }
}
```

• Fluids of different densities
• ```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)
VariableTracerVOFHeight 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.) }
}
```

• Pure gravity wave
• ```1 1 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1.66481717925811447992 }
ApproxProjectionParams { tolerance = 1e-6 }
ProjectionParams { tolerance = 1e-6 }
Refine LEVEL
VariableTracerVOFHeight 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.) }
}
```

• Shape oscillation of an inviscid droplet
• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 1 }

Refine LEVEL

VariableTracerVOFHeight 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)*dL)

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,\$5;}'
k(t)=a*exp(-b*t)*(1.-cos(c*t))
a = 3e-4
b = 1.5

D = DIAMETER
n = 2.
sigma = 1.
rhol = 1.
rhog = 1./1000.
r0 = D/2.
omega0 = sqrt((n**3-n)*sigma/((rhol+rhog)*r0**3))

c = 2.*omega0
fit k(x) 'k-LEVEL' u 3:5 via a,b,c
print "result: ", a, b, c, D
EOF
rm -f fit.log
}
}
```

• Height-Function on parallel subdomains
• ```2 1 GfsSimulation GfsBox GfsGEdge {
#     x = -0.5
x = -0.38
} {
#    Refine (x > 0 && y > 0 ? 4 : 3)
Refine 4
VariableTracerVOFHeight T
VariableCurvature K T
InitFraction T (ellipse (0, 0, 0.2, 0.3))
Time { end = 0 }
OutputSimulation { start = end } stdout
}
```

• Sessile drop
• ```1 0 GfsSimulation GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
Refine LEVEL
VariableTracerVOFHeight T
VariableCurvature K T
SourceTension T 1. K
PhysicalParams { alpha = 1./(T + 0.01*(1. - T)) }
SourceViscosity 0.2/(T + 100.*(1. - T))
InitFraction T (- ellipse (0, 0, 0.3, 0.3))
AdaptFunction { istep = 1 } { cmax = 0 maxlevel = LEVEL } (T > 0 && T < 1)
Time { end = 3 }
EventStop { istep = 10 } K 1e-5 DK
OutputScalarNorm { istep = 10 } v-ANGLE { v = Velocity }
OutputScalarStats { istep = 10 } k { v = (T > 0.05 && T < 0.95 ? K : NODATA) }
OutputScalarSum { start = end } vol { v = T }
#    OutputSimulation { istep = 10 } stdout
OutputSimulation { start = end } end-ANGLE.gfs
}
```

• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
Time { end = 5 dtmax = 1e-2 }
PhysicalParams { L = 1 }
Global {
#include <gsl/gsl_integration.h>

#define G 1.
#define FROUDE 0.1

double vtheta (double r) {
return FROUDE*(r < 0.4)*(1. + cos((r - 0.2)/0.2*M_PI))/2.;
}
double h0p (double r, void * p) {
double vt = vtheta(r);
return vt*(2.*OMEGA + vt/r)/G;
}
double h0 (double r) {
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
double result, error;
gsl_function F;
F.function = &h0p;
gsl_integration_qags (&F, 0, r, 0, 1e-7, 1000, w, &result, &error);
gsl_integration_workspace_free (w);
return result;
}
}
SourceCoriolis 2.*OMEGA
Init {} {
U = - vtheta(sqrt (x*x + y*y))*y/sqrt (x*x + y*y)
V = vtheta(sqrt (x*x + y*y))*x/sqrt (x*x + y*y)
}
ApproxProjectionParams { tolerance = 1e-9 }
ProjectionParams { tolerance = 1e-9 }
Refine 6
OutputScalarSum { istep = 1 } energy-OMEGA { v = Velocity2 }
OutputErrorNorm { istep = 1 } error-OMEGA { v = P } {
s = h0(sqrt (x*x + y*y))
v = E
unbiased = 1
relative = 1
}
OutputSimulation { start = end } end-OMEGA.gfs
}
```

• Creeping Couette flow between cylinders
• ```1 1 GfsSimulation GfsBox GfsGEdge {} {
# the polar coordinates. rx and ry vary between [-0.5:0.5]
# the radius varies between [0.25:0.5]
Metric M {
x = (rx/4. + 0.375)*cos(ry)
y = (rx/4. + 0.375)*sin(ry)
}
VariableTracer T { scheme = none }
SourceDiffusion T 1
Time { dtmax = 1e-2 }
Refine LEVEL
ApproxProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
SourceViscosity 1
EventStop { step = 1e-2 } V 1e-5 DV

OutputScalarNorm { istep = 1 } dv { v = DV }
OutputErrorNorm { start = end } { awk '{print LEVEL,\$5,\$7,\$9}' } { v = V } {
s = {
double r = rx/4. + 0.375;
// the analytical solution for the tangential velocity
return r*((0.5/r)*(0.5/r) - 1.)/((0.5/0.25)*(0.5/0.25) - 1.);
}
v = E
}
OutputSimulation { start = end } end-LEVEL.gfs
Init { start = end } { Rx = rx }
OutputSimulation { start = end } end-LEVEL.txt { format = text }
}
```

• Creeping Couette flow between eccentric cylinders
• ```1 0 GfsSimulation GfsBox GfsGEdge {} {
PhysicalParams { L = 2.5 }
# Upper bound on time, it should converge much before this
Time { end = 100 }
Refine LEVEL
Solid (- ellipse (0.,ECC,R2,R2))
Solid (ellipse (0.,0.,R1,R1))
ApproxProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
SourceViscosity 1
# Fixed outer cylinder and rotating inner cylinder (tangential velocity unity)
SurfaceBc U Dirichlet (x*x + y*y > 1.5*R1*R1 ? 0. : - ay/R1)
SurfaceBc V Dirichlet (x*x + y*y > 1.5*R1*R1 ? 0. :   ax/R1)
# Stop when stationnary
EventStop { step = 1e-2 } U 1e-4 DU

Global {
#include "wannier.c"
double psi (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return p;
}
double ux (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return u;
}
double uy (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return v;
}
}

OutputScalarNorm { istep = 1 } du { v = DU }
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' } { v = Velocity } {
s = {
double p, u, v;
psiuv (x, y - ECC, R1, R2, ECC, 1., 0., &u, &v, &p);
return sqrt (u*u + v*v);
}
v = EU
}
OutputSimulation { start = end } end-LEVEL.gfs
}
```

• Flow between eccentric cylinders using bipolar coordinates
• ```8 8 GfsSimulation GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
#  GModule hypre
Time { end = 100 }

# Bipolar coordinates. see http://en.wikipedia.org/wiki/Bipolar_coordinates
Metric M {
x = sinh(rx/2. + 1.)/(cosh(rx/2. + 1.) - cos (M_PI*ry/4.))
y = sin(M_PI*ry/4.)/(cosh(rx/2. + 1.) - cos (M_PI*ry/4.))
}
Refine LEVEL
ApproxProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
SourceViscosity 1.

Global {
// The radii R1,R2 and center positions X1,X2 below
// correspond to the extent of the computational domain
// i.e. tau=rx/2+1    in [1:1.5]
// and  sigma=pi*ry/4 in [0:2*pi]

#define R1 (1./sinh(1.5))
#define R2 (1./sinh(1.))
#define X1 (1./tanh(1.5))
#define X2 (1./tanh(1.))
#define ECC (X2 - X1)
#include "../wannier.c"
}

EventStop { step = 1e-2 } U 1e-4 DU

#  OutputProjectionStats { istep = 1 } stderr
#  OutputDiffusionStats { istep = 1 } stderr
#  OutputScalarNorm { istep = 1 } stderr { v = DU }

OutputScalarNorm { istep = 1 } du { v = DU }
# Note that the U and V components are defined in the "normalised
# basis" so that U^2+V^2 is independent from the basis (which is why
# we chose to compute the error on the norm of the velocity rather
# than on individual components).
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' } { v = Velocity } {
s = {
double p, u, v;
// we need to rotate coordinates by 90 degrees
psiuv (y, x - X2, R1, R2, ECC, 1., 0., &u, &v, &p);
return sqrt (u*u + v*v);
}
v = EU
}
OutputSimulation { start = end } end-LEVEL.gfs
}
```

• Flow between eccentric cylinders on a stretched grid
• ```2 1 GfsSimulation GfsBox GfsGEdge { y = -0.5 } {
MetricStretch {} { sy = 0.5 }
PhysicalParams { L = 2.5 }
# Upper bound on time, it should converge much before this
Time { end = 100 }
Refine LEVEL
Solid (- ellipse (0.,ECC,R2,R2))
Solid (ellipse (0.,0.,R1,R1))
ApproxProjectionParams { tolerance = 1e-6 }
AdvectionParams { scheme = none }
SourceViscosity 1
# Fixed outer cylinder and rotating inner cylinder (tangential velocity unity)
SurfaceBc U Dirichlet (x*x + y*y > 1.5*R1*R1 ? 0. : - ay/R1)
SurfaceBc V Dirichlet (x*x + y*y > 1.5*R1*R1 ? 0. :   ax/R1)
# Stop when stationnary
EventStop { step = 1e-2 } U 1e-4 DU

Global {
#include "../wannier.c"
double psi (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return p;
}
double ux (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return u;
}
double uy (double x, double y) {
double p, u, v;
psiuv (x, y, R1, R2, ECC, 1., 0., &u, &v, &p);
return v;
}
}

OutputScalarNorm { istep = 1 } du { v = DU }
OutputErrorNorm { start = end } { awk '{ print LEVEL,\$5,\$7,\$9 }' } { v = Velocity } {
s = {
double p, u, v;
psiuv (x, y - ECC, R1, R2, ECC, 1., 0., &u, &v, &p);
return sqrt (u*u + v*v);
}
v = EU
}
OutputSimulation { start = end } end-LEVEL.gfs
}
```

• Rossby--Haurwitz wave
• ```6 12 GfsSimulation GfsBox GfsGEdge {} {
Global {
#define AR 6.37122e6
#define N 4.
#define Umax 50.
#define M (Umax/(N*AR))
#define K M

#define Omega 7.292e-5
#define DTR (M_PI/180.)
#define H0 8e3
#define G 9.81

// Williamson 1992, eq. (137)
static double u0 (double lambda, double phi) {
double cosphi = cos (phi), sinphi = sin (phi);
return AR*M*cosphi + AR*K*pow (cosphi, N - 1)*
cos (N*lambda)*(N*sinphi*sinphi - cosphi*cosphi);
}

// Williamson 1992, eq. (138)
static double v0 (double lambda, double phi) {
double cosphi = cos (phi), sinphi = sin (phi);
return - AR*K*N*pow (cosphi, N - 1)*sinphi*sin (N*lambda);
}

// Williamson 1992, eq. (139)
static double vorticity0 (double lambda, double phi) {
return 2.*M*sin(phi) - K*sin(phi)*pow (cos (phi), N)*(N*N + 3.*N + 2.)*cos (N*lambda);
}

// Williamson 1992, eqs. (140-143)
static double p0 (double lambda, double phi, double t) {
double nu = (N*(3. + N)*M - 2.*Omega)/((1. + N)*(2. + N));
lambda -= nu*t;
double cosphi = cos (phi);
double Aphi = M/2.*(2.*Omega + M)*cosphi*cosphi + K*K/4.*pow (cosphi, 2.*N)*
((N + 1.)*cosphi*cosphi + 2.*N*N - N - 2. - 2.*N*N*pow(cosphi, -2.));
double Bphi = 2.*(Omega + M)*K/((N + 1.)*(N + 2.))*pow(cosphi, N)*
(N*N + 2.*N + 2. - (N + 1.)*(N + 1.)*cosphi*cosphi);
double Cphi = K*K/4.*pow(cosphi,2.*N)*((N + 1.)*cosphi*cosphi - (N + 2.));
return AR*AR*(Aphi + Bphi*cos (N*lambda) + Cphi*cos (2.*N*lambda));
}
}
PhysicalParams { L = 2.*M_PI*AR/4. }
MetricCubed M LEVEL
SourceCoriolis 2.*Omega*sin(y*DTR)

Init {} { (U,V) = (u0(x*DTR,y*DTR), v0(x*DTR,y*DTR)) }

ProjectionParams { tolerance = 1e-8 }
ApproxProjectionParams { tolerance = 1e-8 }

Refine LEVEL

# ~24 days
Time { end = 2073534 dtmax = 1e3 }
#    OutputProgress { istep = 10 } stderr
#    OutputProjectionStats { istep = 1 } stderr

OutputScalarNorm { istep = 10 } v-LEVEL { v = V }
OutputScalarSum { istep = 10 } ec-LEVEL { v = Velocity2 }
OutputScalarSum { istep = 10 } zeta-LEVEL { v = Vorticity }
OutputScalarSum { istep = 10 } p-LEVEL { v = P }
OutputErrorNorm { istart = 1 istep = 10 } eh-LEVEL { v = P/G } {
s = p0(x*DTR,y*DTR,t)/G
v = EH unbiased = 1 relative = 1
}
OutputSimulation { start = end } end-LEVEL.gfs
#    OutputSimulation { istep = 10 } stdout

GModule gfsview
OutputView { start = end } ehp-LEVEL.gnu { format = Gnuplot } ehp.gfv
OutputView { start = end } ehm-LEVEL.gnu { format = Gnuplot } ehm.gfv
OutputView { start = end } h-LEVEL.gnu { format = Gnuplot } h.gfv
OutputView { start = end } href-LEVEL.gnu { format = Gnuplot } href.gfv
}
```

• Simple example of groundwater flow following Darcy's law
• ```1 0 GfsSimulation GfsBox GfsGEdge { x = 0.5 y = -0.5 } {
Time { dtmax = 1 }
Refine LEVEL
AdvectionParams { scheme = none }
Solid (-difference(ellipse(0,0,0.7,0.7), ellipse(0,0,0.3,0.3)))
SourceCoriolis 0 1
EventStop { istep = 1 } U 1e-9

EventList { start = end } {
OutputErrorNorm {} { awk '{print LEVEL, \$5, \$7, \$9}' >> p } { v = P } {
s = atan2 (y, x)
}
OutputErrorNorm {} { awk '{print LEVEL, \$5, \$7, \$9}' >> U } { v = U } {
s = y / (x*x + y*y)
}
OutputErrorNorm {} { awk '{print LEVEL, \$5, \$7, \$9}' >> V } { v = V } {
s = -x / (x*x + y*y)
}
OutputSimulation solution.gfs
}
}
```

• Groundwater flow with piecewise constant permeability
• ```1 0 GfsSimulation GfsBox GfsGEdge { x = 0.5 y = -0.5 } {
Time { dtmax = 1 }
Refine LEVEL
AdvectionParams { scheme = none }
Solid (-difference(ellipse(0,0,0.7,0.7), ellipse(0,0,0.3,0.3)))
SourceCoriolis 0 1
EventStop { istep = 1 } U 1e-9
Variable k
Init {} { k = THETA>-M_PI/3 ? 1 : 0.5 }
PhysicalParams { alpha = k }

EventList { start = end } {
OutputErrorNorm {} { awk '{print LEVEL, \$5, \$7, \$9}' >> p } { v = P } {
s = THETA>-M_PI/3 ? THETA : 2*THETA+M_PI/3
}
OutputErrorNorm {} { awk '{print LEVEL, \$5, \$7, \$9}' >> U } { v = U } {
s = y / (x*x + y*y)
}
OutputErrorNorm {} { awk '{print LEVEL, \$5, \$7, \$9}' >> V } { v = V } {
s = -x / (x*x + y*y)
}
OutputSimulation {} solution.gfs
}
}
```