GfsAxi

From Gerris

Jump to: navigation, search

The GfsAxi solver is the axisymmetric (cylindrical coordinates) version of GfsSimulation; y is the radial coordinate and x the along-axis coordinate. Note that the default position of the center of the reference GfsBox is shifted to (0,0.5) so that the axis of symmetry (y = 0) is aligned with the bottom boundary.

Note that it is important to explicitly set a boundary condition on the bottom boundary of GfsBoxes i.e.

bottom = Boundary

By default the azimuthal velocity (swirl) is zero. The W variable can be used to add a non-zero azimuthal velocity. The corresponding terms must be added to the parameter file. See the boundary layer on a rotating disk test case for an example.

Examples

  • Potential flow around a sphere
  • 1 0 GfsAxi GfsBox GfsGEdge {} {
        Time { end = 0 }
        PhysicalParams { L = 50 }
        AdvectionParams { scheme = none }
        ApproxProjectionParams { tolerance = 1e-10 }
        Refine 4
        Refine (LEVEL + 1./50.*(x*x + y*y)*(4. - LEVEL))
        Global {
    	#define A0 0.5
    	#define U0 1.
        }
        Solid (ellipse (0., 0., A0, A0))
        Init {} {
    	U = U0
    	Phi = {
    	    double r = sqrt (cx*cx + cy*cy);
    	    return U0*A0*A0*A0*cx/(2.*r*r*r);
    	}
        }
        OutputErrorNorm { start = end } { awk '{ print LEVEL " " $7 " " $9}' } { v = U } {
     	s = (dx("Phi") + 1.)
        }
        OutputSimulation { start = end } sim-LEVEL.gfs
    }
    

  • Viscous flow past a sphere
  • 1 0 GfsAxi GfsBox GfsGEdge {} {
        Time { end = 100 }
        PhysicalParams { L = 50 }
        Refine 4
        Refine (LEVEL + 1./50.*(x*x + y*y)*(4. - LEVEL))
        Solid (ellipse (0., 0., A0, A0))
        SourceViscosity 1./RE
        Init {} { U = U0 }
        AdaptGradient { istep = 1 } { cmax = 5e-2 maxlevel = LEVEL } U
        AdaptGradient { istep = 1 } { cmax = 5e-2 maxlevel = LEVEL } V
        AdaptFunction { istep = 1 } { cmax = 1e-2 maxlevel = LEVEL } {
    	return (fabs(dx("U"))+fabs(dy("U")))/fabs(U)*ftt_cell_size (cell);
        }
        EventStop { step = 0.1 } U 1e-3 DU
    
    #    OutputTime { step = 1 } stderr
    #    OutputScalarNorm { step = 1 } stderr { v = DU }
        OutputSimulation { start = end } end-LEVEL-RE.gfs
        OutputLocation { step = 0.1 } {
    	awk 'BEGIN { t = 2.; oldl = -1.; oldt = 0.; } {
              if ($1 != t) { t = $1; x1 = $2; u1 = $7; }
              else {
                x2 = $2; u2 = $7;
                if (u1 <= 0. && u2 > 0.) {
                  l = (u1*x2 - u2*x1)/(u1 - u2) - A0;
                  dl = (l - oldl)/(t - oldt);
                  print t, l, dl;
                  fflush (stdout);
                  oldl = l;
                  oldt = t;
                }
                x1 = x2; u1 = u2;
              }
            }' > l-LEVEL-RE
        } axis
    }
    

  • Mass conservation
  • 1 0 GfsAxi GfsBox GfsGEdge {} {
        Time { end = 1.3 }
        Refine 6
        Init {} { U = 1 }
        VariableTracerVOF T
        VariableTracer T1
        InitFraction T (- ellipse (0, 0.3, 0.1, 0.1))
        InitFraction T1 (- ellipse (0, 0.3, 0.1, 0.1))
        AdaptGradient { istep = 1 } { cmax = 1e-3 minlevel = 4 maxlevel = (x < 0.25 ? 6 : 7) } T1
        OutputScalarSum { istep = 1 end = 0.8 } srt { v = T }
        OutputScalarSum { istep = 1 end = 0.8 } srt1 { v = T1 }
        OutputSimulation { step = 0.2 } stdout
        EventScript { step = 0.2 } {
    	echo "Save stdout { format = Gnuplot }"
        }
        EventScript { start = 1.2 } {
    	echo "Clear"
    	cat vectors.gfv
    	echo "Save vectors.gnu { format = Gnuplot }"
        }
    }
    

  • Mass conservation with solid boundary
  • 1 0 GfsAxi GfsBox GfsGEdge {} {
        Time { end = 0.55 }
        Refine (x < 0. ? 7 : 8)
        Init {} { U = 1 }
        ApproxProjectionParams { tolerance = 1e-6 }
        ProjectionParams { tolerance = 1e-6 }
        Solid (ellipse (0., 0., 0.05, 0.05))
        VariableTracerVOF T
        VariableTracer T1
        InitFraction T (- ellipse (-0.2, 0., 0.05, 0.05))
        InitFraction T1 (- ellipse (-0.2, 0., 0.05, 0.05))
        AdaptGradient { istep = 1 } { cmax = 1e-3 minlevel = 4 maxlevel = (x < 0. ? 7 : 8) } T1
        OutputScalarSum { istep = 1 } srt { v = T }
        OutputScalarSum { istep = 1 } srt1 { v = T1 }
    }
    

  • Boundary layer on a rotating disk
  • 4 3 GfsAxi GfsBox GfsGEdge { x = 0.5 } {
        PhysicalParams { L = 3 }
        SourceViscosity 0.2
    
        # the block below adds the azimuthal velocity (swirl) and
        # associated terms
        VariableTracer W
        SourceDiffusion W 0.2
        Init { istep = 1 } { W = W*(1. - dt*V/y) }
        Source V W*W/y
    
        Refine 5
        Time { end = 12 dtmax = 2e-2 }
        EventStop { step = 0.1 } U 1e-3 DU
    #    OutputTime { istep = 10 } stderr
    #    OutputProjectionStats { istep = 10 } stderr
    #    OutputSimulation { istep = 10 } stdout
        OutputSimulation { start = end } end.txt { format = text variables = U,W }
        OutputSimulation { start = end } end.gfs
    
        EventScript { start = end } {
    	awk 'BEGIN{ nu = 0.2 } {
                   if ($2 != 0. && $1 != "#" && ($1*$1 + $2*$2) < 8.0)
                     print $1/sqrt(nu),-$4/sqrt(nu),$5/$2;
                 }' < end.txt > nu
    	if gnuplot <<EOF ; then :
        set term postscript eps lw 3 color solid 20 enhanced
        set output 'profiles.eps'
        set xlabel '{/Symbol z}'
        set key center right
        plot [0:6]'analytical' u 1:2 w l t '-F({/Symbol z})', 'nu' u 1:2 t 'Gerris', \
                  'analytical' u 1:3 w l t 'G({/Symbol z})', 'nu' u 1:3 t 'Gerris'
    EOF
    	else
    	    exit $GFS_STOP;
    	fi
    
    	if python <<EOF ; then :
    from check import *
    from sys import *
    if (Curve('nu',1,2) - Curve('analytical',1,2)).normi() > 8e-3 or \
       (Curve('nu',1,3) - Curve('analytical',1,3)).normi() > 8e-3 :
        print (Curve('nu',1,2) - Curve('analytical',1,2)).normi()
        print (Curve('nu',1,3) - Curve('analytical',1,3)).normi()
        exit(1)
    EOF
    	else
    	    exit $GFS_STOP;
    	fi
        }
    }
    

  • Axisymmetric spherical droplet in equilibrium
  • 1 0 GfsAxi 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); fflush (stdout); }' > La-LAPLACE-LEVEL
      } { v = Velocity }
      OutputScalarNorm { istep = 1 } {
        awk '{ print MU*$3/(0.8*0.8), $5, $7, $9; fflush (stdout); }' > E-LAPLACE-LEVEL
      } { v = (Tref - T) }
      OutputScalarStats { istep = 1 } {
        awk '{ print MU*$3/(0.8*0.8), $5, $7, $9, $11; fflush (stdout); }' > K-LAPLACE-LEVEL
      } { v = (K - 2.50771) }
      OutputScalarNorm { istep = 1 } {
        awk '{ print MU*$3/(0.8*0.8), $5, $7, $9; fflush (stdout); }' > EK-LAPLACE-LEVEL
      } { v = (T > 0 && T < 1 ? (K - 5.)/2. : 0) }
    }
    

  • Scalings for Plateau--Rayleigh pinchoff
  • 1 0 GfsAxi GfsBox GfsGEdge {} {
       Refine 5
       VariableTracerVOFHeight T
       PhysicalParams { alpha = 1./(T + 1e-2*(1. - T)) }
       VariableCurvature K T Kmax
       VariablePosition Y T y
       InitFraction T (0.2*(1. + 0.1*sin(M_PI*x)) - y)
       SourceTension T 1. K
    #   SourceViscosity 1./20.*(T + 1e-2*(1. - T))
    
       AdaptGradient { istep = 1 } { maxlevel = 6 cmax = 0 } T
       AdaptFunction { istep = 1 } { 
           # adapt down to 14 levels to capture breakup asymptotics and
           # only 10 levels after breakup
           maxlevel = (t < 0.7455 ? 14 : 10)
           cmax = 0.2 
       } (T > 0 && T < 1)*Kmax*dL
    
       Time { end = 0.8 }
    
       OutputScalarStats { istep = 1 } p { v = P }
       # we need more precision for the time output
       OutputScalarStats { istep = 1 } y { v = Y format = "%10.6e" }
       OutputScalarStats { istep = 1 } k { v = K }
       OutputBalance { istep = 1 } balance
       OutputScalarSum { istep = 1 } ke { v = Velocity2*T }
       OutputScalarNorm { istep = 1 } u { v = U }
    #   OutputSimulation { istep = 10 } stdout
      
       GModule gfsview
       OutputView { start = 0.6 istep = 5 end = 0.7455 } { 
           ppm2mpeg -s 640x480 > plateau.mpg 
       } { width = 1280 height = 960 } zoom.gfv
       OutputView { step = 0.2 } plateau-%g.eps { format = EPS } plateau.gfv
       OutputView { start = 0.7451 } plateau-t0.eps { format = EPS } plateau.gfv
       OutputView { start = 0.7451 } zoom-t0.eps { format = EPS } zoom.gfv
    
       EventScript { start = end } {
           # set t0 (breakup time) as the time for which the minimum
           # radius is smaller than the size of the finest cell: 1/2^14
           t0=`awk '{ if ($5 < 1./2.**14) { print $3; exit(0);} }' < y`
           rm -f fit.log
           gnuplot <<EOF
           set term postscript eps lw 3 20 color enhanced
           set key spacing 1.5
           set grid
           set output 'y.eps'
           set xlabel 't_0 - t'
           set ylabel 'r_{min}'
           set logscale
           fit [1e-6:1e-3] a*x**(2./3.) 'y' u ($t0 - \$3):5 via a
           plot [1e-6:][1./2**14:]'y' u ($t0 - \$3):5 ps 0.5 t '', a*x**(2./3.) t 'x^{2/3}'
    EOF
           status=0
           if awk '{if ($1 == "rms" && (NF < 8 || $8 > 4e-5)) {
                      print "rmsy: " $8 > "/dev/stderr"; exit (1); 
                   }}' < fit.log; then :
           else
    	   status=$GFS_STOP;
           fi
           rm -f fit.log
           gnuplot <<EOF
           set term postscript eps lw 3 20 color enhanced
           set key spacing 1.5
           set grid
           set output 'u.eps'
           set xlabel 't_0 - t'
           set ylabel 'u_{max}'
           set logscale
           fit [1e-6:1e-3] a*x**(-1./3.) 'u' u ($t0 - \$3):9 via a
           plot [1e-6:]'u' u ($t0 - \$3):9 ps 0.5 t '', a*x**(-1./3.) t 'x^{-1/3}'
    EOF
           if awk '{if ($1 == "rms" && (NF < 8 || $8 > 8.)) {
                      print "rmsu: " $8 > "/dev/stderr"; exit (1); 
                   }}' < fit.log; then :
           else
    	   status=$GFS_STOP;
           fi
           exit $status;
       }
    }
    GfsBox {
        bottom = Boundary
        left = Boundary
        right = Boundary
    }
    

Personal tools
communication