// NBPROC 4 
// PARAM   -cas 2    -n 20 -ndt 50
// usage :
//  ff-mpirun -np 8 NSNewtonCyl-100-mpi.edp -cas 2    -n 20 -ndt 50 
/* args can are:
   -ns                no output of script 
   -cas 1|2|3         case of benchmark  or Reynorl integer)
   -n nbseg           nb of seg on input edge ..  
   -ndt  ndt         nb time step par period 
   -nNewton nNx      number of of iteration in Newtow loop.    
   -gp      gnuplot ..
   -T final time 
 */
// Author: F. Hecht  .. 
// jan 2012 Stationnary imcompressible Navier Stokes Equation with Newton method.
//  a round a 2d Cylinder 
//  Benchmark  
// http://www.mathematik.tu-dortmund.de/lsiii/cms/papers/SchaeferTurek1996.pdf
//  
//  Thank to Macarena Gómez Mármol <macarena@us.es>
// need version 3.20-2 ****  for pipe pluging 
// otherwise remove line with  pgnuplot 
// corresponding to gnuplot visualization
// --------------------------------------------------
include "getARGV.idp"
include "getARGV.idp"
load "MUMPS_mpi"
load "shell"
load "scotch"
load "pipe"

int cas = getARGV("-cas",2);; // 1, 2 or 3 ..
int n=getARGV("-n",15); //  mesh size ...
int ndt=getARGV("-ndt",25); 
int schema = 2; //  order 1  
bool gplt = usedARGV("-gp")>=0  && (mpirank==0) && (cas!=1) ; // 
real Tend = getARGV("-T",2);;  
bool save=0; 
real cpu0=mpiWtime() ; 
string data="./dd";//"2D-"+cas+"-"+n;
if(mpirank==0) mkdir(data);
string fgnup =data+"/fld2-"+cas+".gp";
string fdata = data+"/sol-cas-";
string fTh = data+"/Th-"+cas+".msh";

pstream  pgnuplot(  gplt  ? "gnuplot": "cat" );
 
//   parameter ... 
real D=0.1;
real H=0.41;
int ncurl=20; 

real[int] colorhsv=[  // color hsv model
  3./6., 1 , 0.5, // dark cyan
  3./6., 1 , 1, //  cyan
  4./6., 1 , 1, //  blue  
  5./6., 1 , 1, //  magenta
  1, 1. , 1, //  red
  1, 0.5 , 1 // light red 
   ];

real[int] vcurl( ncurl*2+1);
{ int n = vcurl.n;
   real d = 100;
   real c = 1.3 ;
    real ccc=0,cc = d/(1- c^(ncurl+1))/ (1-c) ; 
   vcurl[ncurl]=0; 
    for( int i =1 ; i <= ncurl;++i)
      {
        ccc += cc;
        cc *= c ; 
      //  cout << i << " " << ncurl-i << " " << ncurl+i << endl;
        vcurl[ncurl-i] = -ccc;
        vcurl[ncurl+i] = ccc;
      }
   cout << vcurl << endl; 
}
real cx0 = 0.2, cy0 = 0.2; // center of cyl. 
real xa = 0.15, ya=0.2, xe = 0.25,ye =0.2;// point for pressure..
//Definición del dominio 
border fr1(t=0,2.2){x=t; y=0; label=1;}
border fr2(t=0,H){x=2.2; y=t; label=2;}
border fr3(t=2.2,0){x=t; y=H; label=1;}
border fr4(t=H,0){x=0; y=t; label=1;}
border fr5(t=2*pi,0){x=cx0+D*sin(t)/2; y=cy0+D*cos(t)/2; label=3;}

//plot(fr1(n)+fr2(n)+fr3(n)+fr4(n)+fr5(n));
mesh Th;

if(mpirank==0)
{
  Th=buildmesh(fr1(5*n)+fr2(n)+fr3(5*n)+fr4(n)+fr5(-n*3));
  int[int] nupart(Th.nt);
  nupart=0; 
  if(mpisize>1)
    scotch(nupart, Th, mpisize);
  Th=change(Th,fregion= nupart[nuTriangle]);
  savemesh(Th,fTh);
  plot(Th,wait=0);

}
broadcast(processor(0),Th); 
// partition of the mesh ...
// for computation of the Strouhal number
include "func-max.idp"
real hsize = D*pi/n; // size of the mesh ???? 

// bounding box for the plot 
func bb=[[0,H/4],[H*2,3*H/4]];

//  operator 
macro Grad(u1,u2) [ dx(u1),dy(u1), dx(u2),dy(u2)]// 
macro Eps(u1,u2) [ dx(u1),(dy(u1)+dx(u2))*0.5,(dy(u1)+dx(u2))*0.5,dy(u2)]// 
macro UgradV(u1,u2,v1,v2) [ [u1,u2]'*[dx(v1),dy(v1)] , [u1,u2]'*[dx(v2),dy(v2)] ]// 
macro div(u1,u2)  (dx(u1)+dy(u2))//
//  FE Space 
fespace Xh(Th,P2);fespace Mh(Th,P1);
fespace Wh(Th,[P2,P2,P1]);
fespace Rh(Th,P1dc);
Wh [u1,u2,p];
Wh [up1,up2,pp];
Wh [upp1,upp2,ppp];

Wh [vx1,vx2,vxp];
Wh [vy1,vy2,vyp];
varf von3x([u1,u2,p],[v1,v2,q]) = on(3 ,u1= 1); 
varf von3y([u1,u2,p],[v1,v2,q]) = on(3 ,u2= 1); 
vx1[]= von3x(0,Wh,tgv=1);
vy1[]= von3y(0,Wh,tgv=1);
real coefv =1, vtheta=0;  
real Um= 1.5;// max volicite case 2,3 (Rey 100) 
if( cas>3) 
  {
    
    Um = cas*1.5/100.;
    cas =2;
  }

func  Ub = Um*2./3.; 
real nu = 1e-3; 
real mu = 2*nu; //  Formulation in Eps : Eps 
func Rey = Ub*D/nu;
func ccdrag = 2./ square(Ub) /D ;
func ccfreq = D/Ub; 
real freq = 0.3/ccfreq; // frequence theorique
// initial guess..
u1[]=0; 

real T=0; // current time 

// velocity BC .. 
func coefV = (cas<3)? 1 : sin(pi*min(T/8.,1.));
func U1 = 4.*Um*y*(H-y)/(H*H) * coefv ;
func U2 = 0 ;

// stop test for Newton 
real eps=1e-8*Um;
//  choise of time step ..
int  mdt = rint(Tend*freq*ndt);
real dt = Tend/mdt;//  ntd  pas de temps pas periode 
real CFL = Um*dt/hsize; 
if(mpirank==0) cout << " dt =" << dt << " CFL = " << CFL <<  endl; 

verbosity=0;
// BDF2 schema  order 2 in time ... 
real[int] beta=[1.5, -2., +0.5];  
if(schema ==1) 
  beta=[1,-1,0]; // Euler implicit + ruse OP 
  
real[int] ab=1./dt*beta  ;

if(cas ==1) 
{ // static case
  Um = 0.3;
  Tend=0;
  ab=0; // no time .. 
 
}
// init..
int iter=0;
if(mpirank==0) {
        ofstream f(fgnup);// clean the file ..   
 }

varf GStokesl([du1,du2,dp],[v1,v2,q]) =
		   int2d(Th,mpirank) ( 
		                ab[0]*[du1,du2]'*[v1,v2]
		              + mu*(Eps(du1,du2)'*Eps(v1,v2) )
		              - div(du1,du2)*q - div(v1,v2)*dp 
		              - 1e-8*dp*q // stabilization term 
		             )
	    + on(1,3,du1=0,du2=0)
	    + on(2,du2=0) 
	    
;
varf GStokesl0([du1,du2,dp],[v1,v2,q]) =
		   int2d(Th,mpirank) ( 
		                ab[0]*[du1,du2]'*[v1,v2]
		              + mu*(Eps(du1,du2)'*Eps(v1,v2) )
		              - div(du1,du2)*q - div(v1,v2)*dp 
		              - 1e-8*dp*q // stabilization term 
		             )
	   ;
varf Vconvect([du1,du2,dp],[v1,v2,q]) =
  -int2d(Th,mpirank) ( 
    ab[1]* ( [v1,v2]'* [
      convect([u1,u2],-dt,up1),
      convect([u1,u2],-dt,up2)
    ] )
   +  ab[2]* ( [v1,v2]'* [
      convect([u1,u2],-dt*2,upp1),
      convect([u1,u2],-dt*2,upp2)
    ] )
     ); 

varf Vconvect1([du1,du2,dp],[v1,v2,q]) =
  -int2d(Th,mpirank) ( 
    ab[1]* ( [v1,v2]'* [
      convect([u1,u2],-dt,up1),
      convect([u1,u2],-dt,up2)
    ] )
); 
    
  
varf VBC([u1,u2,p],[v1,v2,q]) = on(1,u1=U1,u2=U2); 
real[int] bcl=VBC(0,Wh);
real[int] bl(Wh.ndof),b(Wh.ndof),bcv(Wh.ndof); 
matrix A,A0;
 
 A=GStokesl(Wh,Wh,solver=GMRES); //sparsesolver,master=-1);
//mpiAllReduce(Al,A,mpiCommWorld,mpiSUM);
A0=GStokesl0(Wh,Wh,solver=GMRES);
//mpiAllReduce(A0,A,mpiCommWorld,mpiSUM);

set(A,solver=sparsesolver,master=-1);
Wh [w1,w2,wp]; 
// loop in time if need ..
cout << "Stort\n";
real CPU0=mpiWtime();

while(1) 
{	
    real err=0;
    T += dt;
    iter++;
    real co = coefv;
    coefv  =  coefV;// coef in velocty ..
    if( co != coefv ) bcl=VBC(0,Wh);
    upp1[]=up1[];
    up1[]=u1[];
    if(schema>1)
     { u1[] *=2; u1[] -= upp1[]; } 
     if( ab[2])    
      bcv = Vconvect(0,Wh);
    else 
      bcv = Vconvect1(0,Wh);
    bl = bcv;
    bl += bcl;
       	
	{
	   //mpiAllReduce(bl,b,mpiCommWorld,mpiSUM);
        w1[]=A^-1*bl;
	    u1[] =  w1[];
	    if(mpirank==0)  
	      cout << T << "\t " << n << " rey  =" << Rey << " U1 max = " << u1[].max <<" " ;
	}
	//if( schema==1)
	 //{ u1[] *= 2; u1[] -= up1[]; }   
   if(cas==1) 
    {
        plot([u1,u2],p,  coef = 0.02/u1[].max, bb=bb,wait=1, cmm=" cas 1 "); 
    }
     b  = A0*w1[];
     b += bcv;
     real cdrag = -(b'*vx1[]) ; 
     real clift = b'*vy1[] ; 
		             
    real drag,lift,Cd,Cl,Ta;
    mpiAllReduce(clift,lift,mpiCommWorld,mpiSUM);  
    mpiAllReduce(cdrag,drag,mpiCommWorld,mpiSUM);  
    if ( iter > 2)
    {
    Cd = ccdrag*drag;
    Cl = ccdrag*lift; 
    }
    // end of compute the Drag and lift Coef 
    
    Ta = T/ccfreq; // The adimensional time ..
    real Deltap = p(xa,ya) - p(xe,ye) ; // the Deltap
    real strouhal = AddStrouhal(iter,Ta,Cl);
    AddMaxO2(Cd,iter,mxdrag,mxdragi);   
    AddMaxO2(Cl,iter,mxlift,mxlifti);   
    real Cdx = mxdrag.max, Clx = mxlift.max; 
   if(mpirank==0) {
        {
        ofstream f(fgnup,append);
        f << T << " " << drag << " " << lift << " " 
          << Ta << " " << Cd << " " << Cl    << " " << Deltap <<"   " << mpiWtime()-cpu0 << " " 
          << iter <<  "  " << strouhal << endl;
        }
     real Td = max((int(Ta)-20),0); //  period  = 3 en Ta 
     if(cas==3) Td=0; 
     if(gplt) {
        pgnuplot << "set title \"Rey = " +Rey + ", Strouhal= "+strouhal+",  T = "
                    + T + ", cas = " + cas << ", Cara CFL : "+CFL+"\";\n";
        pgnuplot << "plot ["+Td+":] \""+fgnup+"\" u 4:5 w l t \"Cd\", \""+fgnup
                    +"\"  u 4:6 w l t \"Cl\","+Cdx+",   "+ Clx << endl;
        flush(pgnuplot); 
     }
 
    }
    
    // plot curl of the flow preetty ...
    Rh curlu = -dy(u1)+ dx(u2);
    plot(curlu,viso=vcurl,  fill=1,hsv=colorhsv,
    // bb=bb,
      coef=0.01/Um,cmm = " T = " + T+" / "+Tend+" s + "); 
    
    upp1[] = up1[] -u1[];
    upp2[] = up2[] -u2[];
     real CPU = mpiWtime()-CPU0;
     CPU0=mpiWtime();
    real errt = upp1[].linfty + upp2[].linfty; 
    if(mpirank==0)  cout << " errt = " << errt<< " T = " << T <<  " Strouhal  " << strouhal  <<" Cdx= " << Cdx<< " Clx =  " << Clx << " " << CPU << "s" <<endl; 
    if(mpirank ==0 && save)
    {
        ofstream f(fdata+"-"+iter+".txt");
        f << T << endl;
        f << u1[] << endl;
        f << u2[] << endl;
        f << p[] << endl; 
    }
    if( errt < 1e-5) break; // stop if steady state 
    if(T>Tend*1.0000001) break; // stop after final time ..
}


