GfsSimulation

From Gerris

Jump to: navigation, search

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 }

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

  • 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
    AdaptGradient { istep = 1 } {
	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 radius 1./12.
	#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
    InitFraction T0 (radius*radius - R2(y,z))

    VariableTracerVOF T
    VariableCurvature K T Kmax
    SourceTension T 0.00003 K
    SourceViscosity 2.*radius/Re*rho(T)
    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
    }

    # Dynamic load-balancing
    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
    AdvectionParams {
	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 {} {
	# Start with an arbitrary velocity profile, here twice
	# Bagnold's solution
	U = 2.*UB(y)
	# Start with the hydrostatic pressure profile
	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
        #define radius 0.05
    }

    Refine level

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

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

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

  • Non-linear geostrophic adjustment

    • 1 0 GfsSimulation GfsBox GfsGEdge {} {
  Time { end = 5 dtmax = 1e-2 }
  PhysicalParams { L = 1 }
  Global {
      #include <gsl/gsl_integration.h>
      @link -lgsl -lgslcblas
      
      #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	
    }
}

Retrieved from "http://gfs.sourceforge.net/wiki/index.php/GfsSimulation"
Views
Navigation
communication