# define LEVEL 6
# include "grid/octree.h"
# define DLM_Moving_particle 0
# define adaptive 1
# define NPARTICLES 1

# if adaptive
# define MAXLEVEL (LEVEL + 1)
# endif

# define diam (1.) 
# define Um 1
# define rhoval 1
# define Re 50
# define fs_density_ratio  2. // fluid solid density ratio
# define Ld_ratio 8. // box size-particle diameter ratio
# define Ldomain (Ld_ratio*diam)
# define rhosolid (fs_density_ratio*rhoval) //particle density
# define tval (Um*rhoval*diam)/Re

# define Tc (diam/Um) // caracteristic time scale
# define mydt (Tc/100) // maximum time-step (the time step is also adaptive in time but it won't exceed this value)
# define maxtime (2.)
# define tsave (Tc/200.)

# include "DLMFD_reverse_Uzawa.h"
# include "view.h"
# include "lambda2.h"

 
double deltau;
scalar un[];

int main() {		
  L0 = Ldomain;

  // set time step
  DT = mydt;
  
  /* initialize grid */
  init_grid(1 << LEVEL);

  /* left boundary */
  u.n[left] = dirichlet(Um); 
  u.r[left] = dirichlet(0.); 
  u.t[left] = dirichlet(0.); 
  
  /* right boundary */
  u.n[right] = neumann(0.);   
  u.r[right] = dirichlet(0.);
  u.t[right] = dirichlet(0.);
 
  /* top boundary */
  u.n[top] = dirichlet(0.);
  u.r[top] = dirichlet(Um);  
  u.t[top] = dirichlet(0.);  

  /* bottom boundary */
  u.n[bottom] = dirichlet(0.); 
  u.r[bottom] = dirichlet(Um);  
  u.t[bottom] = dirichlet(0.);
  
  /* front boundary */
  u.n[front] = dirichlet(0.);  
  u.r[front] = dirichlet(0.);   
  u.t[front] = dirichlet(Um);   

  /* back boundary */
  u.n[back] = dirichlet(0.);  
  u.r[back] = dirichlet(0.);  
  u.t[back] = dirichlet(Um);   
  
  p[left]    = neumann(0.);
  p[right]   = dirichlet(0.);

  /* Convergence criteria */
   TOLERANCE = 1e-3; 
  
  run();
}


event init (i = 0) {
  /* set origin */
  origin (0, 0, 0);

  /* set dynamic viscosity */
  const face vector muc[] = {tval,tval,tval};
  mu = muc;

  /* set density of the flow */ 
  const scalar rhoc[] = rhoval;
  rho = rhoc;


  /* We set the initially horizontal velocity to unity everywhere.  */
  foreach(){
    u.x[] = Um;
    un[] = u.x[];
  }

  /* initial condition: particles position */
  particle * p = particles;
  
  init_file_pointers (pdata, fdata, 0);
  
  for (int k = 0; k < NPARTICLES; k++) {
    GeomParameter gp = {0};
    gp.center.x = L0/4;
    gp.center.y = L0/2;
    gp.center.z = L0/2;
    gp.radius = diam/2;

    p[k].g = gp;
   
    /* particle id */
    p[k].pnum = k;
    
    /* density rho_s of the particle */
    p[k].rho_s = rhosolid;

    /* Volume or surface of the particle (circle of sphere) */
    p[k].Vp = 4*pi*pow(gp.radius,3)/3;

    /* total weight of the particle */
    p[k].M = rhosolid*(p[k].Vp);

#if DLM_Moving_particle
    /* The inertia tensor is: */
    /* When in doupt check https://en.wikipedia.org/wiki/List_of_moments_of_inertia */
    /* Ip[0] = Ixx */
    /* Ip[1] = Iyy */
    /* Ip[2] = Izz */
    /* Ip[3] = Ixy */
    /* Ip[4] = Ixz */
    /* Ip[5] = Iyz */
    
    /* For a solid disk: */
    p[k].Ip[0] = 2*(p[k].M)*sq(gp.radius)/5;
    p[k].Ip[1] = p[k].Ip[0];
    p[k].Ip[2] = p[k].Ip[0]; 
    p[k].Ip[3] = 0.;
    p[k].Ip[4] = 0.;
    p[k].Ip[5] = 0.;
    

    /* initial condition: particle's velocity */
    coord c;
    foreach_dimension()
      c.x = 0;
#if TRANSLATION
    p[k].U = c;
#endif
#if ROTATION
    p[k].w = c;
#endif
   p[k].wished_ratio = 0.1;
   p[k].en = 1;
   p[k].vzero = 1;
   compute_wo (p);
#endif

  }
}

event logfile (i++) {
  deltau = change (u.x, un);
  fprintf (stderr, "%d %g %d %d %g \n", i, t, mgp.i, mgu.i, deltau);
}

#if DLM_alpha_coupling
event viscous_term (i++) {
  foreach() {
    foreach_dimension()
      u.x[] += dt*DLM_lambda.x[]/(rho[]*dv());
  }
}
#endif


event adapt (i++) {
  particle * pp = particles;
  int totalcell = totalcells ();
  if (pid() == 0) {
    printf ("total cells = %d\n", totalcell);
    printf ("starting dlmfd subproblem \n");
  }
  DLMFD_subproblem (pp, i, rhoval);
   
  /* Save forces acting on particles before the adapting the mesh */
  sumLambda (pp, fdata, t, dt, flagfield, DLM_lambda, index_lambda, rhoval);

  /* Free particle structures (we dont need them anymore) */
  free_particles (pp, NPARTICLES);
   
  /* Save particles trajectories */
  particle_data (pp, t, i, pdata);
   
#if adaptive
#if DLM_Moving_particle  
  astats s = adapt_wavelet ((scalar *){flagfield_mailleur, u}, (double[]){1e-4, 0.01,0.01,0.01}, maxlevel = MAXLEVEL);
  fprintf (ferr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);

#elif DLM_Moving_particle == 0
  
  astats s = adapt_wavelet ((scalar *){flagfield, u}, (double[]){1e-4,1e-2,1e-2,1e-2}, maxlevel = MAXLEVEL);
  fprintf (ferr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
#endif
#endif
}


event output_data (t+=maxtime/50; t < maxtime) {

  /* char name[80] = "flowsphere"; */
  /* scalar * list = {p, flagfield}; */
  /* vector * vlist = {u, index_lambda}; */
  /* save_data (list, vlist, i, t, name); */
  
  view (fov = 39.6758, quat = {-0.200394,0.310041,0.0701111,0.926715}, tx = -0.0703854, ty = -0.0351928, bg = {0.3,0.4,0.6}, width = 600, height = 600, samples = 1);

  stats statsvelox;
  statsvelox = statsf (u.x);
  clear ();
  box ();
  cells (n = {0,0,1}, alpha = L0/2);
  squares ("u.x", n = {0,0,1}, alpha = L0/2, map = cool_warm, min = statsvelox.min, max = statsvelox.max);
  squares ("u.x", n = {0,1,0}, alpha = L0/2, map = cool_warm, min = statsvelox.min, max = statsvelox.max);
  cells (n = {0,1,0}, alpha = L0/2);
  scalar l2[];
  lambda2 (u, l2);
  isosurface ("l2", -0.0002);
  save ("movie.mp4");
}
 
event lastdump (t = maxtime) {
  dump (file = "dump");
}