GfsAdvection

From Gerris

Revision as of 02:26, 15 May 2009; view current revision
←Older revision | Newer revision→
Jump to: navigation, search

Examples

  • Convergence of the Godunov advection scheme
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
      Time { end = 0.785398 }
      Refine LEVEL
      VariableTracer T { gradient = gfs_center_gradient }
      Init {} {
        T = {
          double r2 = x*x + y*y; 
          double coeff = 20. + 20000.*r2*r2*r2*r2;
          return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
        }
      }
      VariableStreamFunction Psi -4.*(x*x + y*y)
      OutputErrorNorm { start = end } { awk '{ print LEVEL " " $5 " " $7 " " $9}' } { v = T } {
        s = {
          double r2 = x*x + y*y; 
          double coeff = 20. + 20000.*r2*r2*r2*r2;
          return (1. + cos(20.*x)*cos(20.*y))*exp(-coeff*r2)/2.;
        }
      }
      OutputScalarSum { istep = 1 } t { v = T }
    }
    

  • Time-reversed VOF advection in a shear flow
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
        Time { end = 5 }
        Refine 8
    
        # Add tracer T, using a VOF advection scheme.
        # The default scheme is a Van-Leer limited, second-order upwind scheme.
        VariableTracerVOFHeight T
    
        InitFraction T (ellipse (0, -.236338, 0.2, 0.2))
    
        # Maintain U and V with the vortical shear flow field defined by
        # its streamfunction
        VariableStreamFunction {
    	# make sure that the entire field is reinitialised at t = 2.5
    	step = 2.5 
        } Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
        
        # Adapt the mesh dynamically so that at any time the maximum of the gradient
        # of T is less than or equal to zero per cell length
        AdaptGradient { istart = 1 istep = 1 } { cmax = 0 maxlevel = 8 } T
        
        OutputPPM { start = 0 } { convert ppm:- t-0.eps } { v = T }
        OutputPPM { start = 2.5 } { convert ppm:- t-2.5.eps } { v = T }
        OutputPPM { start = 5 } { convert ppm:- t-5.eps } { v = T }
    
        # Add a new variable 
        Variable Tref
    
        # Initialize Tref with the initial shape
        InitFraction { start = end } Tref (ellipse (0, -.236338, 0.2, 0.2))
    
        # Output the norms of the difference between T and Tref, stores that into
        # new variable DT
        OutputErrorNorm { start = end } norms { v = T } {
            s = Tref v = DT
        }
    
        OutputPPM { start = end } { convert ppm:- dt-5.eps } { v = DT }
    
        OutputScalarSum { istep = 1 } { awk '{ print $3,$5-8.743441e-01 }' > t } { v = T }
    }
    

  • Time-reversed advection of a VOF concentration
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
        Time { end = 5 }
        Refine LEVEL
        VariableTracerVOFHeight T
        VariableVOFConcentration C T
        VariableTracer G
        InitFraction T CIRCLE
        Init {} {
    	C = GAUSSIAN*T
    	G = GAUSSIAN
        }
        VariableStreamFunction {
    	step = 2.5 
        } Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
        AdaptGradient { istart = 1 istep = 1 } { cmax = 0 maxlevel = LEVEL } T
        AdaptError { istart = 1 istep = 1 } { cmax = 1e-3 maxlevel = LEVEL } C
        AdaptError { istart = 1 istep = 1 } { cmax = 1e-3 maxlevel = LEVEL } G
        OutputScalarSum { istep = 1 } t { v = T }
        OutputScalarSum { istep = 1 } t1 { v = C }
        OutputScalarSum { istep = 1 } t2 { v = G }
    #    OutputSimulation { istep = 10 } stdout
        
        OutputSimulation { start = 2.5 } half-LEVEL.gfs    
    
        EventScript { start = end } {
    	conservation() 
            {
    	    if awk -v tolerance=$1 'BEGIN { min = 1.; max = -1.; }{ 
                  if ($5 > max) max = $5; 
                  if ($5 < min) min = $5; 
                }END{ if (max - min > tolerance) { print max - min > "/dev/stderr"; exit (1); } }'; then
    		:
    	    else
    		exit $GFS_STOP;
    	    fi
    	}
            conservation 1e-4 < t
            conservation 2e-5 < t1
            conservation 0. < t2
        }
    
        VariableTracer Tref
        InitFraction { start = end } Tref CIRCLE
        OutputErrorNorm { start = end } dt-LEVEL { v = T } {
    	s = Tref
    	v = DT
        }
        OutputErrorNorm { start = end } dt1-LEVEL { v = C } {
    	s = Tref*GAUSSIAN
    	v = DC
        }
        OutputErrorNorm { start = end } dt2-LEVEL { v = G } {
    	s = GAUSSIAN
    	v = DG
        }
        OutputSimulation { start = end } end-LEVEL.gfs    
    }
    

  • Time-reversed advection with curvature-based refinement
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
        Time { end = 5 }
        Refine 8
    
        # Add tracer T, using a VOF advection scheme.
        # The default scheme is a Van-Leer limited, second-order upwind scheme.
        VariableTracerVOFHeight T
    
        # Curvature K of the interface defined by T
        VariableCurvature K T
    
        InitFraction T (ellipse (0, -.236338, 0.2, 0.2))
    
        # Maintain U and V with the vortical shear flow field defined by
        # its streamfunction
        VariableStreamFunction {
    	# make sure that the entire field is reinitialised at t = 2.5
    	step = 2.5 
        } Psi (t < 2.5 ? 1. : -1.)*sin((x + 0.5)*M_PI)*sin((y + 0.5)*M_PI)/M_PI
        
        # Adapt the mesh dynamically using a custom cost function returning
        # the cell size times the local curvature (only for cells cut by
        # the interface).
        # The maximum cost is set to 0.1 i.e. the radius of curvature must be
        # resolved with at least 10 cells.
        AdaptFunction { istep = 1 } {
        	cmax = 0.1 maxlevel = 8
        	cfactor = 2. 
        } (T > 0. && T < 1.)*dL*fabs(K)
        AdaptThickness { istep = 1 } { maxlevel = 8 } T
    
        GModule gfsview
        OutputView { start = 2.5 } t-2.5.eps { format = EPS line_width = 0.5 } curvature.gfv
        
        # Add a new variable 
        Variable Tref
    
        # Initialize Tref with the initial shape
        InitFraction { start = end } Tref (ellipse (0, -.236338, 0.2, 0.2))
    
        # Output the norms of the difference between T and Tref, stores that into
        # new variable DT
        OutputErrorNorm { start = end } norms { v = T } {
            s = Tref v = DT
        }
    
        OutputPPM { start = end } { convert ppm:- dt-5.eps } { v = DT }
    
        OutputScalarSum { istep = 1 } { awk '{ print $3,$5-0.8743665 }' > t } { v = T }
    }
    

  • Rotation of a straight interface
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
        PhysicalParams { L = 2. }
        Refine 3
        VariableTracerVOFHeight T
        Init {} { U = y }
        # fixme: 1e-9 to circumvent a bug in GfsView
        InitFraction T (x - 1e-9)
        Variable Tref
        InitFraction { istep = 1 } Tref (x - t*y)
        OutputErrorNorm { istep = 1 } error { v = T } { s = Tref }
        Time { end = 5 }
        GModule gfsview
        OutputScalarSum { istep = 1 } t { v = T }
        OutputView { step = 1 } rotate-%g.gnu { format = Gnuplot } rotate.gfv
        OutputView { start = end } cells.gnu { format = Gnuplot } cells.gfv
        EventScript { start = end } {
    	gnuplot <<EOF
    set term postscript eps color lw 3 18
    set output 'error.eps'
    set ylabel 'Maximum volume fraction error'
    set xlabel 'Time'
    set key top left
    set logscale y
    set grid y
    plot [][1e-3:]'error' u 3:9 w l t 'n=4', 'error.n1' u 3:9 w l t 'n=1'
    EOF
            for t in 0 1 2 5; do
    	    gnuplot <<EOF
    set term postscript eps color lw 2 solid
    set output 'rotate-$t.eps'
    unset border
    unset xtics
    unset ytics
    unset key
    set size ratio -1
    plot 'cells.gnu' w l lc 0, 'rotate-$t.gnu' w l lw 2 lc 1, 'n1-$t.gnu' w l lc 2 lw 2
    EOF
            done
    
            status=0;
            if awk '{if ($5 != 2.) { print $5 > "/dev/stderr"; exit (1); }}' < t ; then :
            else
                status=$GFS_STOP;
            fi
    
            paste error error.ref > tmp
    	if awk '{
              if ($9 > $20) {
                print $3,$9,$20 > "/dev/stderr"
                exit (1);
              }
            }' < tmp; then :
    	else
    	    status=$GFS_STOP;
            fi
    	
    	exit $status
        }
    }
    GfsBox {
        left = Boundary { 
    	BcDirichlet U y
    	BcAngle T (90. + atan2(1.,t)*180./M_PI)
        }
        right = Boundary { 
    	BcDirichlet U y
    	BcAngle T (90. - atan2(1.,t)*180./M_PI)	
        }
        top = Boundary {
    	BcDirichlet U y
    	BcAngle T (180. - atan2(1.,t)*180./M_PI)
        }
        bottom = Boundary {
    	BcDirichlet U y
    	BcAngle T atan2(1.,t)*180./M_PI
        }
    }
    

  • Comparison between the explicit and implicit diffusion schemes
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
        VariableTracer T { scheme = none }
    #    SourceDiffusion T 1 { beta = 1 }
        SourceDiffusionExplicit T 1
        Solid (-sphere(0,0,0,0.4))
        Refine 5
        SurfaceBc T Dirichlet 1
    
        # Just checks that it works with hypre
        GModule hypre
    
        Time { end = 0.02 dtmax = 9e-5 }
    #    OutputTime { istep = 1 } stderr
        OutputSimulation { start = 0.01 } end.gfs
        OutputSimulation { start = 2.5e-3 step = 2.5e-3 } {
    	awk '{
               if ($1 != "#")
                 print sqrt($1*$1+$2*$2),$8;
            }' > prof
        } { format = text }
    }
    

  • Comparison between explicit and implicit diffusion schemes on concentration tracer
  • 1 0 GfsAdvection GfsBox GfsGEdge {} {
        VariableTracerVOF T
        VariableTracer T1 { scheme = none }
        VariableVOFConcentration C T
        Refine 5
        InitFraction T (-difference(sphere(0,0,0,0.4),sphere(0,0,0,0.1)))
        Init {} {
    	C = GAUSSIAN*T
    	T1 = C
        }
    
    #    GModule hypre
        SourceDiffusion T1 T/(T+1.e9*(1.-T)*(1.-T)) {beta = 1}
        SourceDiffusionExplicit C T/(T+1.e9*(1.-T)*(1.-T))
    #    SourceDiffusion C T/(T+1.e9*(1.-T)*(1.-T))
    
    
        Time { end = 0.02 dtmax = 9e-5 }
    
        OutputSimulation { start = 0.01 } end.gfs
        OutputScalarSum {istep = 10}  T1vol { v = T1 }
        OutputScalarSum {istep = 10}  Cvol { v = C }
        EventScript { start = end } {
    	conservation() 
            {
    	    if awk -v tolerance=$1 'BEGIN { min = 1.; max = -1.; }{ 
                  if ($5 > max) max = $5; 
                  if ($5 < min) min = $5; 
                }END{ if (max - min > tolerance) { print max - min > "/dev/stderr"; exit (1); } }'; then
    		:
    	    else
    		exit $GFS_STOP;
    	    fi
    	}
            conservation 2.1e-7 < T1vol
            conservation 2.1e-7 < Cvol
        }
        OutputSimulation {start = 0.005} {
    	awk '{
               if ($1 != "#" && $8 > 0.)
                 print sqrt($1*$1+$2*$2),$12,$13;
            }' > prof
        } { format = text}
    }
    

  • Numerical viscosity for vorticity/streamfunction formulation
  • 1 2 GfsAdvection GfsBox GfsGEdge {} {
      Time { end = 2 }
      Refine LEVEL
      VariableTracer Omega
      Init {} {
        U0 = (- cos (2.*M_PI*x)*sin (2.*M_PI*y))
        V0 = (sin (2.*M_PI*x)*cos (2.*M_PI*y))
        # we need (U,V) to compute the vorticity
        U = U0
        V = V0
      }
      # In the general case the vorticity cannot be consistently
      # defined directly (because it needs to verify discrete properties
      # on the circulation)
      # Initialising from the velocity field (U,V) as done below should be safe
      Init {} { Omega = Vorticity }
      VariablePoisson Psi -Omega
      VariableStreamFunction { istep = 1 } Psi1 Psi
      OutputScalarNorm { istep = 1 } divLEVEL { v = (dx("U") + dy("V")) }
      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
      }
    }
    

  • Advection of a cosine bell around the sphere
  • 6 12 GfsAdvection GfsBox GfsGEdge {} {
      PhysicalParams { L = 2.*M_PI/4. }
      MetricCubed M LEVEL
    
      # Use alternative implementation
      # MetricCubed1 M
      # MapFunction {
      #     x = atan2 (X, Z)*180./M_PI
      #     y = asin (CLAMP(Y,-1.,1.))*180./M_PI
      # }
    
      Time { end = 1 }
      Refine LEVEL
      VariableTracer T { 
          gradient = gfs_center_gradient 
          cfl = 1
      }
      Global {
          #define DTR (M_PI/180.)
          double lambda1 (double lambda, double theta, double lambdap, double thetap) {
    	  /* eq. (8) */
    	  return atan2 (cos (theta)*sin (lambda - lambdap),
    	                cos (theta)*sin (thetap)*cos (lambda - lambdap) - 
                            cos (thetap)* sin (theta));
          }
          double theta1 (double lambda, double theta, double lambdap, double thetap) {
    	  /* eq. (9) */
    	  return asin (sin (theta)*sin (thetap) + cos (theta)*cos (thetap)*cos (lambda - lambdap));
          }
          double bell (double lambda, double theta, double lc, double tc) {
    	  double h0 = 1.;
    	  double R = 1./3.;
    	  double r = acos (sin(tc)*sin(theta) + cos (tc)*cos (theta)*cos (lambda - lc));
    	  return r >= R ? 0. : (h0/2.)*(1. + cos (M_PI*r/R));
          }
          double bell_moving (double lambda, double theta, double t) {
    	  double lambda2 = lambda1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
    	  double theta2 = theta1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
    	  double lc = 3.*M_PI/2. - U0*t, tc = 0.;
    	  return bell (lambda2, theta2, lc, tc);
          }      
      }
      Init {} { T = bell_moving(x*DTR,y*DTR,0) }
      VariableStreamFunction Psi -U0*(sin (y*DTR)*cos (DTR*ALPHA) - cos (x*DTR)*cos (y*DTR)*sin (DTR*ALPHA))
      AdaptGradient { istep = 1 } { cmax = 1e-4 maxlevel = LEVEL } T
      OutputTime { istep = 10 } stderr
    #  OutputSimulation { istep = 10 } stdout
      OutputSimulation { start = end } end-LEVEL-ALPHA.gfs
      OutputErrorNorm { istep = 1 } { awk '{ print LEVEL,$3,$5,$7,$9}' > error-LEVEL-ALPHA } { v = T } {
          s = bell_moving(x*DTR,y*DTR,t)
          v = E
          relative = 1
      }
      OutputScalarSum { istep = 1 } t-LEVEL-ALPHA { v = T }
      OutputScalarSum { istep = 1 } area-LEVEL-ALPHA { v = 1 }
      EventScript { start = end } {
          ( cat isolines.gfv
            echo "Save isolines.gnu { format = Gnuplot }"
            echo "Clear"
            cat reference.gfv
            echo "Save reference.gnu { format = Gnuplot }"
    	echo "Clear"
            cat zero.gfv
            echo "Save zero.gnu { format = Gnuplot }"
          ) | gfsview-batch2D end-LEVEL-ALPHA.gfs
          cat <<EOF | gnuplot
            set term postscript eps lw 2 18 color
            set output 'isolines-LEVEL-ALPHA.eps'
            set size ratio -1
            set xlabel 'Longitude'
            set ylabel 'Latitude'
            unset key
            plot [60:120][-30:30]'isolines.gnu' w l, 'reference.gnu' w l, 'zero.gnu' w l
    EOF
          fixbb isolines-LEVEL-ALPHA.eps
          rm -f isolines.gnu reference.gnu zero.gnu
    
          # check mass conservation
          if awk '
            BEGIN { min = 1000.; max = -1000.; }{ 
              if ($5 < min) min = $5; 
              if ($5 > max) max = $5; 
            }
            END {
              if (max - min != 0.)
                exit (1);
            }' < t-LEVEL-ALPHA; then
    	  exit 0
          else
    	  exit $GFS_STOP
          fi
      }
    }
    

Personal tools
communication