//  run with MPI:  ff-mpirun -np 4 script.edp
// time on 4 proc times: compile 0.013393s, execution 2.28842s,  mpirank:0
// time on 1 proc times: compile 0.012316s, execution 4.92302s,  mpirank:0
// NBPROC 4

load "MUMPS_mpi"

verbosity=0; 
real ttgv=1e30;

/*
  Incompressible Navier Stokes 
    with Taylor-Hood Finite element
    No linearity : Newton methode 
    continuation on Reynols Number
    Mesh adaptation 
*/
real  reymax = 9000; // ok < 125000 
mesh Th=square(8,8);
Th=change(Th,fregion=nuTriangle%mpisize); 
fespace Xh(Th,P2);
fespace Mh(Th,P1);
fespace XXMh(Th,[P2,P2,P1]);
XXMh [u1,u2,p];
XXMh [v1,v2,q]; 

macro div(u1,u2) (dx(u1)+dy(u2))//
macro grad(u1,u2) [dx(u1),dy(u2)]//
macro ugrad(u1,u2,v) (u1*dx(v)+u2*dy(v)) //
macro Ugrad(u1,u2,v1,v2) [ugrad(u1,u2,v1),ugrad(u1,u2,v2)]//

solve Stokes ([u1,u2,p],[v1,v2,q],solver=sparsesolver,tgv=ttgv,master=-1) =
    int2d(Th,mpirank)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            - p*div(v1,v2)-q*div(u1,u2)
           )
  + on(3,u1=4*x*(1-x),u2=0) 
  + on(1,2,4,u1=0,u2=0);

 Xh uu1=u1,uu2=u2;  
if(mpirank==0)
plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[uu1,uu2],wait=0);

Xh psi,phi;


solve streamlines(psi,phi,master=-1) = 
      int2d(Th,mpirank)( dx(psi)*dx(phi) + dy(psi)*dy(phi))
   +  int2d(Th,mpirank)( -phi*(dy(u1)-dx(u2)))
   +  on(1,2,3,4,psi=0);
if(mpirank==0)
plot(psi,wait=0);
int i=0;
real  nu=1./100.;
real dt=0.1;
real alpha=1/dt;



/* NL 
 varf   vNS ([u1,u2,p],[v1,v2,q],solver=Crout,init=i) =
    int2d(Th)(
             alpha*( u1*v1 + u2*v2) 
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            + p*dx(v1)+ p*dy(v2)
            + dx(u1)*q+ dy(u2)*q
            + Ugrad(u1,u2,u1,u2)'*[v1,v2]
           )   
  + on(3,u1=1,u2=0) 
  + on(1,2,4,u1=0,u2=0) 
*/

XXMh [up1,up2,pp];
varf   vDNS ([u1,u2,p],[v1,v2,q]) =
    int2d(Th,mpirank)(
            
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            + p*dx(v1)+ p*dy(v2)
            + dx(u1)*q+ dy(u2)*q
            + Ugrad(u1,u2,up1,up2)'*[v1,v2]
            + Ugrad(up1,up2,u1,u2)'*[v1,v2]
           )
  + on(1,2,3,4,u1=0,u2=0) 
;


varf   vNS ([u1,u2,p],[v1,v2,q]) =
    int2d(Th,mpirank)(
          
            + nu * ( dx(up1)*dx(v1) + dy(up1)*dy(v1)
            +  dx(up2)*dx(v2) + dy(up2)*dy(v2) )
            + pp*q*(0.000001) 
            + pp*dx(v1)+ pp*dy(v2)
            + dx(up1)*q+ dy(up2)*q
            + Ugrad(up1,up2,up1,up2)'*[v1,v2]//'
	      )
  + on(1,2,3,4,u1=0,u2=0) 
  ;

for(real re=100;re<=reymax;re *=2)
  { 
    
    real lerr=0.04;
    
    if(re>8000) lerr=0.01;
    if(re>10000) lerr=0.005; 
    for(int step=0;step<2;step++)
      {
		  
	   if( mpirank==0) {
		    Th=change(Th,fregion=1);
			Th=adaptmesh(Th,[u1,u2],p,err=lerr,nbvx=100000);
		    Th=change(Th,fregion=nuTriangle%mpisize);			
		   } 
	    broadcast(processor(0),Th);
      
	//plot(Th,wait=0);
	[u1,u2,p]=[u1,u2,p];
	[up1,up2,pp]=[up1,up2,pp];
	
	for (i=0;i<=20;i++)
	  {
	    nu =1./re;
	    up1[]=u1[];
	    real[int] b = vNS(0,XXMh);
	    matrix Ans=vDNS(XXMh,XXMh,tgv=ttgv);
	    set(Ans,solver=sparsesolver,master=-1);
	    real[int] w = Ans^-1*b;
	    u1[] -= w;
	    if( mpirank==0)  cout << step << " iter = "<< i << "  " << w.l2 <<  " rey = " << re << endl;
	    if(w.l2<1e-6) break; 
	    // uu1=u1;uu2=u2;
	    //plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[uu1,uu2]);  
	    
	  } ;
      }
    uu1=u1;uu2=u2;
    streamlines;
    if( mpirank==0) 
    plot(coef=0.2,cmm="rey="+re+" [u1,u2] et psi  ",psi,[uu1,uu2],wait=0,nbiso=20,ps="cavity-"+re+".ps");  
  }
  
