#include "navier-stokes/low-mach.h"
#include "two-phase.h"
#include "tension.h"
#include "view.h"
#define DIVERGENCE 0

scalar divu[];

u.n[top] = neumann (0.);
u.t[top] = neumann (0.);
p[top] = dirichlet (0.);

u.n[right] = neumann (0.);
u.t[right] = neumann (0.);
p[right] = dirichlet (0.);

int maxlevel;
double R0 = 0.05;

int main (void) {
  rho1 = rho2 = 1;
  mu1 = mu2 = 1;

  f.sigma = 0.03;

  for (maxlevel = 4; maxlevel <= 9; maxlevel++) {
    init_grid (1 << maxlevel);
    run();
  }
}

#define circle(x,y,R) (sq(R) - sq(x) - sq(y))

event init (i = 0) {
  fraction (f, circle (x, y, R0));
}

double exact (double t) {
  return R0 + t;
}

event output (t += 0.01) {
  double R = sqrt (4./pi*statsf (f).sum);
  double err = fabs (R - exact (t)) / exact (t);
  fprintf (stderr, "maxres %d %g %g %g %g\n",
      maxlevel, t, R, exact (t), err);
}

void exact_velocity (vector u, face vector uf, scalar divu) {
  foreach_face() {
    coord o = {x,y,z};
    double mag = sqrt (sq(x) + sq(y));
    uf.x[] = o.x / (mag + 1e-10);
  }

  foreach() {
    coord o = {x,y,z};
    double mag = sqrt (sq(x) + sq(y));
    foreach_dimension()
      u.x[] = o.x / (mag + 1e-10);
  }

  foreach() {
    divu[] = 0.;
    foreach_dimension()
      divu[] += (uf.x[1] - uf.x[]);
    divu[] /= -Delta;
  }
}

void reset_velocity (vector u, face vector uf) {
  foreach()
    foreach_dimension()
      u.x[] = 0.;

  foreach_face()
    uf.x[] = 0.;
}

#if DIVERGENCE
event tracer_diffusion (i++) {
  exact_velocity (u, uf, drhodt);
  reset_velocity (u, uf);
}

event end_timestep (t = end) {
  vector ue[];
  face vector ufe[];
  exact_velocity (ue, ufe, drhodt);

  scalar e[];
  foreach()
    e[] = norm (u) - norm (ue);
  norm n = normf (e);

  fprintf (stderr, "error %d %g %g %g\n", maxlevel, n.avg, n.rms, n.max);
}
#else
event end_timestep (i++) {
  exact_velocity (u, uf, drhodt);
}
#endif

#if TREE
event adapt (i++) {
  scalar ff[];
  foreach()
    ff[] = f[];
  adapt_wavelet ({f}, (double[]){1e-3}, maxlevel);
}
#endif

event stop (t += 0.1; t <= 0.7) {
  scalar magu[];
  foreach()
    magu[] = norm (u);

#if TREE
  restriction ((scalar *){u});
#endif
  clear();
  view (tx = -0.5, ty = -0.5);
  draw_vof ("f", lw = 1.5);
  squares ("magu", spread = -1);
  vectors ("u", scale = 0.001, level = 5);
  save ("movie.mp4");
}