/**
# Collapse of granular columns with basal friction
We propose an implementation of the Jop Pouliquen Forterre µ(I) rheology with friction at the base. This is the counterpart  of the [no slip](http://basilisk.fr/sandbox/M1EMN/Exemples/granular_column.c) case.
 
# Code 
Includes and definitions
*/

#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "vof.h" 
// Domain extent
#define LDOMAIN 5.
// heap definition
double  H0,R0,D;
//film size
double  Lz;
//passive fluid
#define RHOF 1e-3
#define mug  1e-4
#define SIGMA 0.01
// Maximum refinement
#define LEVEL 8
char s[120];
FILE * fpf;
scalar f[];
scalar * interfaces = {f}; 
face vector alphav[];
face vector muv[];
scalar rhov[];
scalar eta[];
scalar etadudy[];
scalar muIp[];
vector u1[];
double tauw;
// note the v


// Boundary conditions
p[top] = dirichlet(-RHOF*LDOMAIN);
pf[top] = dirichlet(-RHOF*LDOMAIN);
u.n[top] = neumann(0);
u1.t[bottom] = neumann(0);
/**
We may impose no slip 
*/
//u.t[bottom] = dirichlet(0);
/** 
but here we impose a friction at the wall
$$\tau = \mu_w p$$
so that we impose a Neumann BC: 
$$\frac{\partial u}{\partial y}|_0 = \frac{\mu_w p}{\eta}$$
where $\eta$ is computed by the $\mu(I)$ law in the bulk. Note the minus sign (the normal os $-y$)), and note that we do some regularisation as we impose the stress in the reverse direction of the velocity
which may be zero (two possible regularizations). The solution is dependant on the regularisation parameter which should be not too small 'eps'.
*/
double eps=3e-3,mu_w=0.48;
//u.t[bottom] = neumann(mu.y[] ? -(mu_c*p[]/mu.y[])*u1.x[]/sqrt( sq(u1.x[]) + sq(eps)) : 0.);
u.t[bottom] = neumann(mu.y[] ? -(mu_w*p[]/mu.y[] * u1.x[]/(fabs(u1.x[] + eps)))  : 0.);
/**
see [couette with stress](http://basilisk.fr/sandbox/M1EMN/Exemples/couette_muw.c) case.*/

int main() {
  L0 = LDOMAIN;
  // number of grid points
  N = 1 << LEVEL;
  // maximum timestep
  DT = 2e-2;
  TOLERANCE = 1e-3;
// Initial conditions for the heap
  H0=1;
  R0=1.4286;
  // Grain size
  D=1./200;
  const face vector g[] = {0.,-1.};
  a = g;
  alpha = alphav;
  mu = muv;
  rho = rhov;
  Lz = LDOMAIN;
  run(); 
}

event init (t = 0) {
  fraction (f, min(H0 - y, R0 - x));
/*
initialisation of hydrostatic pressure for granular phase  
*/
  foreach()
    p[] = f[]*(H0-y);
}

#define rho(f) ((f) + RHOF*(1. - (f)))

event properties (i++) {
  trash ({alphav}); 
  foreach() {
    eta[] = mug;
    if (p[] > 0.) {
      double D2 = 0.;
      foreach_dimension() {
	double dxx = u.x[1,0] - u.x[-1,0];
	double dxy = (u.x[0,1] - u.x[0,-1] + u.y[1,0] - u.y[-1,0])/2.;
	D2 += sq(dxx) + sq(dxy);
      }
      if (D2 > 0.) {
	D2 = sqrt(2.*D2)/(2.*Delta); // this D2 is sqrt(2) D2
	double In = D2*D/sqrt(p[]);
	double muI = .48 + .25*In/(.279 + In);
	double etamin = sqrt(D*D*D);
	eta[] = max((muI*p[])/D2, etamin);
	eta[] = min(eta[],100);      
	muIp[] = muI*p[]; }
    }
  }
  boundary ({eta,muIp});
  scalar fa[];
  foreach()
    fa[] = (4.*f[] + 
	    2.*(f[-1,0] + f[1,0] + f[0,-1] + f[0,1]) +
	    f[1,1] + f[-1,1] + f[1,-1] + f[-1,-1])/16.;
  boundary ({fa});
  foreach_face() {
    double fm = (fa[] + fa[-1,0])/2.;
    muv.x[] = (fm*(eta[] + eta[-1,0])/2. + (1. - fm)*mug);
    // mu.x[] = 1./(2.*fm/(eta[] + eta[-1,0]) + (1. - fm)/mug);
    alphav.x[] = 1./rho(fm);
  }
  foreach()  
  rhov[] = rho(fa[]);
 // boundary_normal ({mu,alpha});
 boundary ({muv,alphav,rhov});
}
/**
save stress and save old field
*/
event sylvain (i++) { 
  foreach()
    etadudy[] = muv.y[]*(u.x[0,0] - u.x[0,-1])/(Delta);
  boundary({etadudy});

  foreach()  
     u1.x[] = u.x[];
 boundary ({u1.x});
}

void mg_print (mgstats mg)
{
   if (mg.i > 0 && mg.resa > 0.)
        fprintf (stderr, "#   %d %g %g %g\n", mg.i, mg.resb, mg.resa,
                 mg.resb > 0 ? exp (log (mg.resb/mg.resa)/mg.i) : 0);
}

event logfile (i++) {
  stats s = statsf (f);
  fprintf (stderr, "%g %d %g %g %g %g\n", t, i, dt, s.sum, s.min, s.max - 1.);
  mg_print (mgp);
  mg_print (mgpf);
  mg_print (mgu);
  fflush (stderr);
}

event interface (t = {0,  0.8371, 1.5068, 1.6742, 2.5113 }) {
  //output_facets (f, fpf); 
  char s[120];
  sprintf (s, "field-%g.txt", t);
  FILE * fp = fopen (s, "w");
  output_field ({f,p,u,uf,pf,etadudy,muIp}, fp, linear = true);
  fclose (fp);
}

#if 0
event movie (t += 0.05) {
  static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
  scalar l[];
  foreach()
    l[] = level;
  output_ppm (l, fp1, min = 0, max = LEVEL, 
	      n = 2048, box = {{0,0},{Lz,Lz}});

  foreach()
    l[] = f[]*(1+sqrt(sq(u.x[]) + sq(u.y[])));
  boundary ({l});
  static FILE * fp2 = popen ("ppm2mpeg > velo.mpg", "w");
  output_ppm (l, fp2, min = 0, max = 1.5, linear = true, 
	      n = 1024, box = {{0,-1},{3,1.5}});

  static FILE * fp3 = popen ("ppm2mpeg > f.mpg", "w");
  foreach()
    l[] = f[]*p[];
  output_ppm (l, fp3, min = 0, linear = true,
	      n = 2048, box = {{0,0},{Lz,Lz}});
}
 event pictures (t==3) {
  output_ppm (f, file = "f.png", min=0, max = 2,  spread = 2, n = 2048, linear = true,
  box = {{0,0},{2,2}});
}
#endif
/**

#Plots
 

Check that final $ p$ at the base multiplied by $\mu_s=0.48$ (chosen value) is the computed $\mu(I)p$

~~~gnuplot 
set xlabel "x"
  p[:3]'field-2.5113.txt'u 1:($3>.95? $2:NaN) t 'heap','' u 1:($2==0.00976562? 0.48*$4:NaN) t' mu p', ''u 1:($2==0.00976562? $11:NaN) t'mu(I)p'
~~~

 
Check that   $\mu(I)p$ is $\eta \partial u/\partial y$. this is not  when there is regularisation (which appears here when the flow is very slow)


~~~gnuplot 
set xlabel "x"
  p[:3]'field-2.5113.txt' u 1:($3>.95? $2:NaN)t 'heap',''u 1:($2==0.00976562? $10:NaN) t'eta du/dy' ,''u 1:($2==0.00976562? $11:NaN) t'mu(I)p'
~~~
 

*/