# Grain size
Define D 0.04623

# Density
Define RHOF 0.001
Define RHO(T) (T + RHOF*(1. - T))

# Viscosity of light phase
Define MUF 0.0001

# Domain extent
Define LDOMAIN 32.

# Initial conditions
Define H0 6.26
Define R0 1.

# Maximum refinement
Define LEVEL 10

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
    }
}
GfsBox {
    top = Boundary {
	# shift the reference pressure by the hydrostatic pressure of
	# the light phase, i.e. P = 0 at the bottom boundary in the
	# light phase. I am not sure whether this makes a real
	# difference.
	BcDirichlet P -RHOF*LDOMAIN
	BcNeumann V 0
    }
    bottom = Boundary {
	# no-slip at the bottom
	BcDirichlet U 0
    }
    # lateral symmetry
    left = Boundary
    # no slip wall on the right-hand-side
    right = Boundary {
	BcDirichlet V 0
    }
}