/**
# Shape oscillation of an inviscid droplet with density ratio of 1/10

A two-dimensional elliptical droplet (density ratio 1/10) is
released in a large domain. Under the effect of surface-tension the
shape of the droplet oscillates around its (circular) equilibrium
shape. The fluids inside and outside the droplet are inviscid so
ideally no damping of the oscillations should occur. As illustrated on
the figures some damping occurs in the simulation due to numerical
dissipation.

This simulation is also a stringent test case of the accuracy of the
surface tension representation as no explicit viscosity can damp
eventual parasitic currents. 

We use a the Navier--Stokes solver with VOF interface tracking and
surface tension using FastP* method for solving poisson equation. */

#include "fast-navierstokes.h"
#include "vof.h"
#include "tension.h"

/**
The interface is represented by the volume fraction field *c*. */

scalar c[], * interfaces = {c};

/**
The diameter of the droplet is 0.2. The density inside the droplet is
one and outside 10^-1^. */

#define D 0.2
#define rho(c) (clamp(c,0,1)*(1. - 1e-1) + 1e-1)

FILE * fp = NULL;

int LEVEL;

int main() {
  
  /**
  The surface tension is unity. We cleanup existing files and vary the
  level of refinement. */      
  c.sigma = 1;
  remove ("error");
  remove ("laplace");
/** 
We choose FastP* method for solving poisson equation*/
  poisson_solver=2;
/** 
We set rho0 parameter as smaller density */
  rho0=1e-1;
  /** 
We reducing the *TOLERANCE* to imporve resutls*/
TOLERANCE=1e-6;

  /**
We will vary the level of refinement to study convergence. */

  for (LEVEL = 5; LEVEL <= 8; LEVEL++) {
    N = 1 << LEVEL;
    
    /**
    We open a file indexed by the level to store the time evolution of
    the kinetic energy. */

    char name[80];
    sprintf (name, "k-%d", LEVEL);
    fp = fopen (name, "w");
    run();
    fclose (fp);
  }

  /**
  We use *grep* to filter the lines generated by gnuplot containing
  the results of the fits (see below). */

  system ("grep ^fit out >> log");
}

/**
We have the option of using some "smearing" of the density jump. */

#if 0
#define cf c
#else
scalar cf[];
#endif

/**
The density is variable. We allocate a new field to store its
inverse. */

face vector alphav[];

event init (i = 0) {

  /**
  We initialise the shape of the interface, a slightly elliptic droplet. */

  fraction (c, D/2.*(1. + 0.05*cos(2.*atan2(y,x))) - sqrt(sq(x) + sq(y)));
#ifndef cf
  foreach()
    cf[] = c[];
#endif

  /**
  The density is variable, and defined by the function below. */

  alpha = alphav;
}

/**
The density is defined at each timestep via the *properties()* event
declared by the Navier--Stokes solver. */

event properties (i++) {

  /**
  When using smearing of the density jump, we initialise *cf* with the
  vertex-average of *c*. */

#ifndef cf
  foreach()
    cf[] = (4.*c[] + 
	    2.*(c[0,1] + c[0,-1] + c[1,0] + c[-1,0]) +
	    c[-1,-1] + c[1,-1] + c[1,1] + c[-1,1])/16.;
  boundary ({cf});
#endif

  /**
  The inverse of the density $\alpha$ is then given by the
  face-averaged value of *cf* and the arithmetic average of density
  defined by *rho()*. */

  foreach_face() {
    double cm = (cf[] + cf[-1])/2.;
    alphav.x[] = 1./rho(cm);
  }
  boundary ((scalar *){alphav});  
}

/**
By default tracers are defined at $t-\Delta t/2$. We use the *first*
keyword to move VOF advection before the *logfile* output i.e. at
$t+\Delta/2$. This improves the results. */

event vof (i++, first);

/**
At each timestep we output the kinetic energy. */

event logfile (i++; t <= 1) {
  double ke = 0.;
  foreach (reduction(+:ke))
    ke += sq(Delta)*(sq(u.x[]) + sq(u.y[]))*rho(cf[]);
  fprintf (fp, "%g %g %d\n", t, ke, mgp.i);
  fflush (fp);
}

/**
At the end of the simulation, we use gnuplot to fit a function of the form
$$
k(t) = ae^{-bt}(1-\cos(ct))
$$
to the kinetic energy. This gives estimates of the oscillation
pulsation *c* and of the damping *b*.*/

event fit (t = end) {
  FILE * fp = popen ("gnuplot 2>&1", "w");
  fprintf (fp, 
           "k(t)=a*exp(-b*t)*(1.-cos(c*t))\n"
           "a = 3e-4\n"
           "b = 1.5\n"
	   "\n"
           "D = %g\n"
           "n = 2.\n"
           "sigma = 1.\n"
           "rhol = 1.\n"
           "rhog = 1./10.\n"
           "r0 = D/2.\n"
           "omega0 = sqrt((n**3-n)*sigma/((rhol+rhog)*r0**3))\n"
	   "\n"
           "c = 2.*omega0\n"
           "fit k(x) 'k-%d' via a,b,c\n"
	   "level = %d\n"
	   "res = D/%g*2.**level\n"
           "print \"fit \", res, a, b, c, D\n"
	   "\n"
	   "set table 'fit-%d'\n"
	   "plot [0:1] 2.*a*exp(-b*x)\n"
	   "unset table\n"
	   "\n"
	   "set print 'error' append\n"
	   "print res, c/2./omega0-1., D\n"
	   "\n"
	   "set print 'laplace' append\n"
	   "empirical_constant = 30.\n"
	   "print res, (1./(b**2.*D**3.))*empirical_constant**2, D\n"
	   "\n",
	   D, LEVEL, LEVEL, L0, LEVEL);
  pclose (fp);
}


#if TREE
event adapt (i++) {
  adapt_wavelet ({c,u}, (double[]){5e-3,1e-3,1e-3}, LEVEL);
  event ("properties");
}
#endif

/**
## Results

~~~gnuplot Evolution of the kinetic energy as a function of time for the spatial resolutions (number of grid points per diameter) indicated in the legend. The black lines are fitted decreasing exponential functions.
set output 'k.png'
set xlabel 'Time'
set ylabel 'Kinetic energy'
set logscale y
plot [0:1][8e-5:1e-3]'k-8' t "51.2" w l, 'k-7' t "25.6" w l,               \
  'k-6' t "12.8" w l, 'k-5' t "6.4" w l,			       \
  'fit-8' t "" w l lt 7, 'fit-7' t "" w l lt 7, 'fit-6' t "" w l lt 7, \
  'fit-5' t "" w l lt 7
~~~

 An equivalent Laplace number can then be computed as:
$$La=\frac{\sigma D}{\rho\nu^2}=\frac{\sigma C^2}{\rho b^2 D^3}$$
The equivalent Laplace number depends on spatial resolution as
illustrated below.

~~~gnuplot Equivalent Laplace number estimated from the numerical damping of kinetic energy.
set output 'laplace.png'
set xlabel 'Diameter (grid points)'
set ylabel 'Equivalent Laplace number'
set grid
plot 'laplace' t "" w p pt 5 ps 2
~~~
*/