#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/perfs.h"
#include "tracer.h"

// Passive tracer
scalar s[];
scalar * tracers = {s};

// Physical and geometrical parameters

double U0;        // Jet inlet velocity
double R_d;       // Jet radius
double Re;        // Reynolds number

FILE * fpmax; 

face vector muv[];

int main() {

  // Flow parameters (Re=80 and R_d=0.03 works also well
  Re      = 53;      // Reynolds number 
  R_d     = 0.015;   // Jet radius 
  L0      = 1.;      // Box size
  U0      = 1.;      // Inlet velocity
  
  TOLERANCE = 1e-3 [*];

  
  
  // Boundary conditions
  
  // Bottom inlet: top-hat velocity profile
  u.n[bottom] = dirichlet(U0*(x > -R_d && x < R_d));
  u.t[bottom] = dirichlet(0.);

  // Top outlet
  u.n[top] = (x > -5*R_d && x < 5*R_d) ? neumann(0.) : dirichlet(0.);
  
  // Tracer injected with the jet
  s[bottom] = dirichlet(U0*(x > -R_d && x < R_d));
  s[top] = (x > -5*R_d && x < 5*R_d) ? neumann(0.) : dirichlet(0.);

  u.n[left]  = dirichlet(0.);
  u.t[left]  = dirichlet(0.);
  u.n[right] = dirichlet(0.);
  u.t[right] = dirichlet(0.);
  
  // Grid initialization
  N=256;  
  origin (-L0/2, 0);
  init_grid(N);
  
  char param_dim[80];
  sprintf (param_dim, "param_dim.txt");
  FILE * fparam = fopen(param_dim, "w");
  fprintf (fparam, "%g %g %g %d\n",L0,R_d,U0,N);
  fclose (fparam);

  mu=muv;

  fpmax =  fopen("log.dat", "w"); 

  run();
}

event properties (i++)
{
  foreach_face() {
    muv.x[] = fm.x[]*U0*2*R_d/Re;
  }
}

event init (t = 0) {
    fprintf(stderr, "Starting new simulation.\n");
}

event logfile (i++) {
  fprintf (stderr, "%d %g \n", i, t);
  fprintf (fpmax, "%d %g \n", i, t);
}

event ppm_output (t = 0; t += 0.1; t <= 100) {
  char name[80];
  sprintf (name, "uX.mp4");
  output_ppm (u.x, file = name, n = 512, min = -U0, max = U0, linear = true);
  
  char name1[80];
  sprintf (name1, "uY.mp4");
  output_ppm (u.y, file = name1, n = 512, min = -U0, max = +U0, linear = true);

  char name2[80];
  sprintf (name2, "s.mp4");
  output_ppm (s, file = name2, n = 512, min = 0., max = U0, linear = true);

  char name3[80];
  sprintf (name3, "p.mp4");
  output_ppm (p, file = name3, n = 512, min = 0, max = 0.1*U0*U0, linear = true);
}

event profile (t = end) {
  printf ("----END----\n");
}