#define SHOW_ERROR 1

#if SHOW_ERROR
#include "embed.h"
#endif
#include "axi.h"
#define ro 1.2
#define veloc(r) (1./r)
#define fpos(t) sqrt(sq(ro) +2.*t)

#include "advection.h"
#include "vof.h"
#include "curvature.h"

scalar f[];
scalar * tracers = NULL;
scalar * interfaces = {f};

int main()
{
  X0 = -0.5;
  Y0 = 1.0;
  L0 = 1.0 [0];
  DT = 1.0 [0];
  N = 64;
  run();
}

event init (i = 0) {

  fraction(f, ro-y);

#if SHOW_ERROR
  foreach()
    cs[] = 1.0;
  foreach_face()
    fs.x[] = 1.0;
#endif  

  foreach_face(y)
    u.y[] = veloc(y)*fm.y[];

  dt = dtnext (timestep (u, DT));
}

event prof_pos (i += 5; t <= 0.8) {
  scalar pos[];
  position (f, pos, {0,1.});
  fprintf (stdout, "%g %g %g \n", t, statsf(pos).max, fpos(t));
  
}