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

double tau_y = 0.2;

int main() {
  periodic (top);
  size (0.5);
  DT = HUGE [0];
  
  stokes = true;
  TOLERANCE = 1e-5;
  
  u.t[left] = dirichlet(0);

  const face vector g[] = {0.,1.};
  a = g;
  const face vector muc[] = {1.,1.};
  mu = muc;
  const face vector tauc[] = {tau_y,tau_y};
  tau_0 = tauc;

  for (N = 8; N <= 64; N *= 2)
    run();
}

double analytical (double x)
{
  double h = 0.5, dpdy = 1., mu = 1.;
  double U = sq(tau_y - h*dpdy)/(2.*mu*dpdy);
  double X = h - tau_y/dpdy;
  return (x < X ? dpdy*(h*x - sq(x)/2.) - tau_y*x : U);
}

scalar un[];

event logfile (t += 0.1; i <= 1000) {
  double du = change (u.y, un);
  if (i > 0 && du < 1e-7)
    return 1; /* stop */
  //  fprintf (stderr, "%g %d %g %d %d\n", t, i, du, mgp.i, mgu.i);
}

event profile (t = end) {
  if (N == 16) {
    foreach()
      fprintf (stdout, "%g %g %g %g %g %g %g\n", x, y, u.x[], u.y[], p[],
	       x - Delta/2., (u.y[] - u.y[-1])/Delta);
    printf ("\n");
  }
  scalar e[];
  foreach()
    e[] = u.y[] - analytical(x);
  norm n = normf (e);
  fprintf (ferr, "%d %g %g %g\n", N, n.avg, n.rms, n.max);
}