GfsAxi

From Gerris

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
}