GfsSimulation

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 15:12, 10 October 2010
Popinet (Talk | contribs)

← Previous diff
Current revision
GeordieMcBain (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.
<examples/> <examples/>

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

Personal tools
communication