/**
# Bursting bubble

With this code, we want to simulate a symplified version of a bursting
bubble in a water bath. We based this simulation on the work of
Duchemin et al. published in 2002. */

#include "axi.h"
#include "navier-stokes/centered.h"
// #define mu(f) ((mu1*mu2)/(clamp(f,0,1)*(mu2-mu1)+mu1))
#include "two-phase.h"
#include "tension.h"
#include "tag.h"

/**
Definition of the geometrical parameters. The 2 radii correspond
respectively to the radius of the bubble and to the radius of a smooth
interface. */

#define LEVEL 14
#define L0 10
#define R1 1.

/**
The parameters L1, L2 and L3 are used to define the hole between the upper 
interface and the spherical droplet. Once the hole is done, we smooth the 
interface by adding a circle between the upper interface and the drop. This 
circle as to be continuous with the upper interface but also with the bubble 
interface. $2\times L2$ defines the distance between the top of the bubble and 
the interface. L1 defines the radius of the smoothing circle, but also the 
distance between the center of this cylinder with the upper interface.*/

#define L2 0.05
#define L3 0.15
#define L1 (-(sq(R1+2.*L2)-sq(R1)+sq(L3))/(4.*(L2+R1)))

#define Sigma 1.

uf.n[left] = 0.;
uf.n[right] = 0.;
uf.n[top] = 0.;
uf.n[bottom] = 0.;

double uemax = 0.3;

int main(int argc, char *argv[]) {
  size (L0);
  origin (-L0/2., 0.);
  init_grid (1 << (9));

  /**
  By default, the value of the Laplace number is 10 000. But, we can
  change that with an initial input argument. */
  double La = 10000.;
  if (argc >= 2){
      La = atof (argv[1]);
    }

  /**
  The Ohnesorge number is defined such that: $Oh = \frac{1}{\sqrt{La}} */
  
  double Oh = sqrt(1./La);

  /**
  We then define the properties of the fluids. In this case, the fluid
  1 correspond to water, and the fluid 2 correspond to air.*/

  rho1 = 1., mu1 = Oh;
  rho2 = 1./998., mu2 = Oh/55;
  
  /**
  The surface tension is defined with the Ohnesorge number*/

  f.sigma = Sigma;

  /**
  We print the characteristics of the fluid in the logfile. */
  fprintf (stderr, "props %f %f %f %f %f\n", rho1, mu1, La, R1, Sigma);
  
  run();
}

double geometry(double x, double y) {

  /**
  The geometry need the definition of 2 circles (C*) and 3 lines
  (D*). Then we need to proceed correctly to the intersection and the
  union of those parameters.*/
  
  double C1 = sq(x + R1 + 2*L2) + sq(y) - sq(R1);
  double C2 = sq(x - L1) + sq(y - L3) - sq(L1);
  double D1 = - x - 1e-8;
  double D2 = y - L3;
  double D3 = - x + (2*L1);

  /**
  Our geometry is:
  $$geom = C1 \cup (D3 \cap (D1 \cup D2))$$ */

  double D1D2 = min(D1, D2);
  double D1D2D3 = max(D1D2, D3);
  double D1D2D3C1 = min(D1D2D3, C1);
  double D1D2D3C1C2 = min(-D1D2D3C1, C2);
  return -D1D2D3C1C2;
} 

/**
We initialise the geometry of the interface.

![Initial interface of the simulation. The length of the rim is half of the 
size of the hole](bursting/initial.png)
*/
  
double intemax = 0.005;

event init (t = 0) {
  double iteration = 0;
  do {
    iteration++;
    fraction(f, geometry(x,y))
  } while( adapt_wavelet({f,u}, (double []){intemax,uemax,uemax,uemax},
     maxlevel = LEVEL, 5).nf != 0 && iteration <= 10);

  static FILE * fp = fopen("initial.png", "w");
  output_ppm(f, fp, 400, min = 0, max = 1, box = {{-2.5,0},{0.5,3}});
  static FILE * fp2 = fopen("initial.gfs", "w");
  output_gfs(fp2);
}


/**
We use an adaptive mesh. Our refinement criteria are based on the
interface and the velocity. */

event adapt (i++) {
  adapt_wavelet ({f,u}, (double []){intemax,uemax,uemax,uemax},
		 maxlevel = LEVEL, 5);
}

/**
We output the shape of the interface at regular intervals. */


/**
This gives the following image

~~~gnuplot Evolution of the interface
unset key
set size ratio -1
plot [-2:2]'interface' u 2:1 w l, 'interface' u (-$2):1 w l lt 1
~~~

We output the position of the interface of the fluid allong the axis of
symmetry at every timestep. We track this position with the height
function. We also compute the corresponding velocity. */

static double tprev = -1, xprev = -1,xDropPrev = -1;
static double compt = 0, toto = 0;

event droplets (i ++; t<=2) {
  scalar m[];
  foreach()
    m[] = f[] > 1.5e-1;
  int n = tag(m);

  /**
  We want to observe the formation of the first drop. We detect when the tag 
  function return more than 1 phase of the internal fluid. We then output 15 
  steps in gfs file to observe the behaviour of the fluid at this moment.*/

  if (n>1 && t>0.15) {
    compt++;
    if (compt>=1 && compt<=15) {
      char rupture[80];
      sprintf(rupture, "rupture-%07ld.gfs",i);
      FILE* fprup = fopen (rupture, "w");
      output_gfs(fprup);
    }
  }  

  /**
  We tag the droplet, when there is some.*/

  double v[n];
  coord b[n];
  double dropVelo[n];

  for (int j = 0; j < n; j++)
    dropVelo[j] = v[j] = b[j].x = b[j].y = b[j].z = 0.;

  foreach_leaf()
    if (m[] > 0) {
      int j = m[] - 1;
      v[j] += dv()*f[];
      dropVelo[j] += dv()*f[]*u.x[];
      coord p = {x,y,z};
      foreach_dimension()
      b[j].x += dv()*f[]*p.x;
    }

  /**
  We compute the droplets velocity. We need at least 2 steps with a droplet to 
  have a correct velocity. We add for that an extra variable, which is 
  counting the number of steps with the presence of a drop. We also compute 
  the volume of the first drop.*/

  double xDrop, veloDrop, volume, radius;

  for (int j = 1; j < n; j++) {

    /**
    Computation of the drop velocity, and the drop position.*/

    xDrop = b[j].x/v[j];
    veloDrop = dropVelo[j]/v[j];

    /**
    Computation of the drop radius.*/ 

    volume = v[j];
    radius = sqrt(2*volume/M_PI);

    fprintf (fout, "%d %g %d %g %g %g %g %g %d\n", i, t,
       j, v[j], xDrop, veloDrop, volume, radius, n);
  }
  fflush (fout); 

/**
We also compute the velocity of the top of the jet. We use for that the height 
function.*/

  vector h[];
  double xMax = -HUGE;;
  heights (f, h);
  foreach()
    if (y <= Delta && h.x[] != nodata && m[]==0) {
      double xi = x + height(h.x[])*Delta;
      if (xi > xMax)
        xMax = xi;
    }
  double velo = tprev >= 0 ? (xMax - xprev)/(t - tprev) : 0.;
  double deltaT = t - tprev;
  printf ("%d %g %g %g %g\n", i, t, deltaT, xMax, velo);
  fflush (stdout);

  if (xMax > 0. && xprev < 0.) {
    fprintf (stderr, "end %g %g\n", t, velo);
    char name[80];
    sprintf(name, "jetX-0.gfs");
    FILE* fp = fopen (name, "w");
    output_gfs (fp);
  }  

  /**
  The stop condition of the simulation is based on 2 conditions. First, we 
  need to overpass the axi x=0. Secondly, we need to observe at least 2 steps 
  consecutive steps with a drops.*/

  if (xMax>=0. && n>1) {
    fprintf(stderr, "firstDrop %g %g %g %g\n", t, veloDrop, volume, radius);
    return 1;
  }

  tprev = t, xprev = xMax, xDropPrev = xDrop;
}

/**
This gives the following figure

~~~gnuplot Evolution of the velocity of the centerline
set xlabel 'Time'
set ylabel 'Velocity'
set grid
plot [0.01:]'out' u 1:4 w l t ''
~~~
*/

event interface(i+= 10) {
  static FILE * fp = popen ("ppm2mpeg > tracer.mpeg", "w");
  static FILE * fp2 = popen ("ppm2mpeg > grid.mpeg", "w");
  output_ppm(f, fp, 512, min=0, max=1, box = {{-2.5,0},{0.5,3}});
  scalar l[];
  foreach()
    l[] = level;
  output_ppm(l,fp2, min=5, max=LEVEL, box = {{-2.5,0},{0.5,3}});
}

/**
We can follow the evolution of the simulation with a video

![The bursting bubble](bursting/tracer.mpeg)
*/

#if 0
event gfsview (i += 10) {
  static FILE * fp = popen ("gfsview2D bursting.gfv", "w");
  output_gfs (fp);
}
#endif

event logfile (i += 10) {
  fprintf (stderr, "%d %g %d %d %d\n", i, t, mgp.i, mgpf.i, mgu.i);
}

/**
We output the final state of the simulation. In the final state, we include 
the tag function, to see where basilisk identify the different phase.*/

event finalState (t = end) {
  scalar m[];
  foreach()
    m[] = f[] > 1.5e-1;
  int n = tag(m);
  output_gfs (file = "final.gfs");
  static FILE * fp = fopen("final.ppm", "w");
  output_ppm(f, fp, 512, min=0, max=1, box = {{-2.5,0},{0.5,3}});
  dump (file = "dump");
}

/**
We can control the refinement parameter by loocking at the evolution 
of the mesh during the simulation.

![Evolution of the mesh](bursting/grid.gif) */