//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "MUMPS_mpi"
////////////////////
//  parameters

real ttgv=1e30;
//string ssparams="nprow=1, npcol="+mpisize;
////////////////////
real nu=0.01,dt=0.3;
real alpha=1./dt,alpha2=sqrt(alpha);

int nn=7;

mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);
Vh2 ux,uz,p2;
int[int] rup=[0,2],  rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;

mesh3 Th=buildlayers(Th2,nn,
  zbound=[zmin,zmax],
  // region=r1, 
  labelmid=rmid, 
  reffaceup = rup,
  reffacelow = rdown);

fespace VVh(Th,[P23d,P23d,P23d,P13d]);
fespace Vh(Th,P23d);
fespace Ph(Th,P13d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM
  
cout << "initilisation" << endl;

real t1,t2,t3,t4;

t1=clock();
varf vStokes([u1,u2,u3,p],[v1,v2,v3,q]) = 
  int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3) //' for emacs
             - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
 + on(2,u1=1.,u2=0,u3=0)
 + on(1,u1=0,u2=0,u3=0)
 ;
matrix A=vStokes( VVh, VVh, tgv=ttgv );
t1=clock()-t1;

t4=clock();
set(A,solver=sparsesolver );
t4=clock()-t4;

t2=clock();
real[int] b= vStokes(0,VVh);
t2=clock()-t2;

VVh [u1,u2,u3,p];
VVh [X1,X2,X3,Xp];
VVh [x1,x2,x3,xp]=[x,y,z,0];

t3=clock();
u1[]= A^-1 * b;
t3=clock()-t3;

ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5",wait=1);
macro XX1() (x-u1*dt)//
macro XX2() (y-u2*dt)//
macro XX3() (z-u3*dt)//

  varf vNS([uu1,uu2,uu3,p],[v1,v2,v3,q]) = 
  int3d(Th)( alpha*(uu1*v1+uu2*v2+uu3*v3) + nu*(Grad(uu1)'*Grad(v1) +  Grad(uu2)'*Grad(v2) +  Grad(uu3)'*Grad(v3)) //' for emacs
  - div(uu1,uu2,uu3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
  + on(2,uu1=1,uu2=0,uu3=0)
  + on(1,uu1=0,uu2=0,uu3=0)
 
  +  int3d(Th,optimize=1,qforder=4)(   alpha*(  convect([u1,u2,u3],-dt,u1)*v1  +   convect([u1,u2,u3],-dt,u2)*v2  +   convect([u1,u2,u3],-dt,u3)*v3 )  ) ;
  //   +  int3d(Th,optimize=1)(   alpha*(  u1(X1,X2,X3)*v1  +  u2(X1,X2,X3)*v2  +  u3(X1,X2,X3)*v3 )  ) ;
//  +  int3d(Th,optimize=1)(   alpha*(  u1(XX1,XX2,XX3)*v1  +  u2(XX1,XX2,XX3)*v2  +  u3(XX1,XX2,XX3)*v3 )  ) ;
//+  int3d(Th,optimize=1)(   alpha*(  u1(x,y,z)*v1  +  u2(x,y,z)*v2  +  u3(x,y,z)*v3 )  ) ;
//+  int3d(Th,optimize=1)(   alpha*(  u1*v1  +  u2*v2  +  u3*v3 )  ) ;

real time1,time2=0.,time3=0.,time4,time5;

time1=clock();
A = vNS( VVh, VVh, tgv=ttgv);
time1=clock()-time1;

time5=clock();
set(A,solver=sparsesolver,verb=1);
time5=clock()-time5;

time4=clock();
real t=0;
int nbtemps=5;
for(int i=0;i<nbtemps;++i)
  {
    real time2tmp,time3tmp;
    t += dt;
    cout << " iteration " << i << " t = " << t << endl;
    X1[]=x1[]+u1[]*(-dt);
    //    verbosity=1000;

    time2tmp=clock();
    b=vNS(0,VVh);
    time2=time2+(clock()-time2tmp);

    time3tmp=clock();
    u1[]= A^-1 * b;
    time3=time3+(clock()-time3tmp);

    ux= u1(x,0.5,y);
    uz= u3(x,0.5,y);
    p2= p(x,0.5,y);
    plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=0);
    if(i%5==6)
    {
      //exec("mkdir dd");
      string prefu="dd/pastix-nn-"+nn+"+u-"+(100+i);
      string prefp="dd/pastix-nn-"+nn+"p-"+(100+i);
      savemesh(Th,prefu+".mesh");
      savemesh(Th,prefp+".mesh");
     
      ofstream file(prefu+".bb"); 
      ofstream filep(prefp+".bb"); 
      Ph up1=u1,up2=u2,up3=u3,pp=p;
      file << "3 1 3 "<< up1[].n << " 2 \n";
      filep << "3 1 1 "<< pp[].n << " 2 \n";
      for (int j=0;j<up1[].n ; j++)  
	{
	  file << up1[][j] <<" " <<up2[][j] <<" "<< up3[][j] <<"\n";
	  filep << pp[][j] <<  endl; 
	}  
    }
  }
time4=clock()-time4;

//plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=1);
cout << "============ CPU TIME ==========================" << endl;
cout << "= Intialisation Stokes :::         =" << endl;
cout << "= matrix                           =" << t1<< endl;
cout << "= Factorization                    =" << t4 << endl; 
cout << "= second member                    =" << t2 << endl;
cout << "= solving                          =" << t3<< endl;
cout << "= all                              =" << t1+t2+t3+t4 << endl; 
cout << "============ CPU TIME ==========================" << endl;
cout << "= Navier Stokes :::                =" << endl;
cout << "= matrix                           =" << time1<< endl;
cout << "= Factorization                    =" << time5 << endl; 
cout << "= second member by step (dt)       =" << time2/nbtemps << endl;
cout << "= solving by step (dt)             =" << time3/nbtemps<< endl;
cout << "= all step solving                 =" << time4 << endl;
cout << "= solving Navier Stokes            =" << time1+time2+time3+time4+time5 << endl; 
cout << "============ CPU TIME ==========================" << endl;
