#include "grid/multigrid3D.h"
#include "embed.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/perfs.h"
#include "vof.h"
#include "tension.h"
#include "two-phase.h"
#include "view.h"
#include "maxruntime.h"
#include <math.h>
#include "common.h"
#include "reduced.h"

#define mu(f) (1./(clamp(f, 0, 1)*(1./mu1 - 1./mu2) + 1./mu2))

// Simulation properties to change
const double MAXTIME = 2.0;
const int LEVEL = 12;

// Geometrical properties of experiment
const double PI = 3.1415926;
const double DIAMETER = 0.1953, WIDTH = 0.2 [1], HEIGHT = 8.0; 
const int NX = 40; //this should be equal to height/width
const int NY = 1, NZ = 1;

const double VLiqIn = 1.028949;
const double VAirIn = 2.312438;
const double Vrad = 0.2;
const double Dinj = 0.014;
#define NHOLES 32 // makes 2x this number of holes, two rings of them
const int Nholes = NHOLES;
double angles[NHOLES];

// Global variables
bool master;
scalar f0[], vInjY[], vInjZ[];

// Boundary conditions
u.n[top] = dirichlet(0.0);
u.n[bottom] = dirichlet(0.0);
u.n[front] = dirichlet(0.0);
u.n[back] = dirichlet(0.0);
u.n[left] = dirichlet(f0[] > 1e-3 ? VAirIn:VLiqIn);
u.n[right] = neumann(0.0);

u.t[top] = dirichlet(0.0);
u.t[bottom] = dirichlet(0.0);
u.t[front] = dirichlet(0.0);
u.t[back] = dirichlet(0.0);
u.t[left] = dirichlet(vInjY[]);

u.r[top] = dirichlet(0.0);
u.r[bottom] = dirichlet(0.0);
u.r[front] = dirichlet(0.0);
u.r[back] = dirichlet(0.0);
u.r[left] = dirichlet(vInjZ[]);

u.n[embed] = dirichlet(0.0);
u.t[embed] = dirichlet(0.0);
u.r[embed] = dirichlet(0.0);

f[left] = (f0[]);
p[left] = neumann(0);
pf[left] = neumann(0);

f[right] = neumann(0);
p[right] = dirichlet(0.);
pf[right] = dirichlet(0.);

void updateAperatures()
{
  vertex scalar phi[];

  foreach_vertex()
  {
    phi[] = sq(DIAMETER/2.0) - sq(y) - sq(z);
  }

  fractions(phi, cs, fs);
  fractions_cleanup(cs, fs);

  // Define angles for injection holes

  for (int i = 0; i < Nholes; i++) {
    angles[i] = i*PI/(Nholes/2);
  }
}

int main (int argc, char * argv[]) {
  maxruntime(&argc, argv);
  #if !_MPI
  master = true;
  #else
  master = (mpi_rank == 0);
  #endif
  
  // Air-water properties at 0.25 MPa operating conditions
  
  rho1 = 2.35 [0]; // air density (average) kg/m^3, 2.79 at inlet, 1.91 at outlet
  rho2 = 997.561 [0]; // water density kg/m^3
  mu1 = 1.85508e-5; // air viscocity Pa-s
  mu2 = 8.8871e-4; // water viscocity Pa-s
  f.sigma = 0.072; // surface tension N/m
  G.x = -9.81; // gravity in x-direction m/s^2

  // Simulation domain size setup
  
  size(HEIGHT);
  dimensions(nx = NX, ny = NY, nz = NZ);
  origin (0., -WIDTH/2.0, -WIDTH/2.0);
  init_grid(pow(2, LEVEL));

  TOLERANCE = 1e-3 [*];
  CFL = 0.5;

  run();
}

event init (t = 0)
{ 
  updateAperatures();

  // Initialize the injection holes at the inlet with a void of 1 
  // and the corresponding radial velocity

  scalar f00[], f01[];

  for(int i = 0; i < Nholes; i++) 
  {
    // Make first set of holes, 32 around circumference right around the edge
    fraction(f00, sq(Dinj/2.0) - sq(y - (DIAMETER-Dinj)/2.0*cos(angles[i])) 
                                - sq(z - (DIAMETER-Dinj)/2.0*sin(angles[i])));
    foreach(){
    f0[] = f0[] + f00[];
    vInjY[] = vInjY[] - Vrad*cos(angles[i])*f00[];
    vInjZ[] = vInjZ[] - Vrad*sin(angles[i])*f00[];
    }

    // Make second set of holes, 32 around circumference right around the first set
    fraction(f01, sq(Dinj/2.0) - sq(y - (DIAMETER-3.0*Dinj)/2.0*cos(angles[i] + PI/(Nholes))) 
                            - sq(z - (DIAMETER-3.0*Dinj)/2.0*sin(angles[i] + PI/(Nholes))));
    foreach(){
    f0[] = f0[] + f01[];
    vInjY[] = vInjY[] - Vrad*cos(angles[i])*f01[];
    vInjZ[] = vInjZ[] - Vrad*sin(angles[i])*f01[];
    }
  }

  boundary((scalar*){u});

  foreach()
  {
      // Initialize velocity everywhere
      u.x[] = (cs[] > 0.) ? VLiqIn:0.;
    }

}

event endImage (t = MAXTIME) 
{
  view(camera = "front", fov = 15., width = 800, ty = 0., tx = -0.5,
    height = 600, bg = {1,1,1});
  box();
  squares("u", alpha = 0.0001, min = 0, max = 3);
  colorbar(horizontal = true, label = "V (m/s)", min = 0, max = 3);
  save("u.ppm");
  clear();
}