// resolution de NS avec THERMIQUE contre une paroi PC 3
// PYL  jv 07 dec 17
// inspire de l'exemple NSP1P2.edp et transmission.edp
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC3/FF++_PC3/poindarretkimarche.edp
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC3/index.html
//
// run
// freefem++  poindarretkimarche.edp
//
//

exec("echo poindarret kimarche2");

real s0=clock();
real h0=.33;   //hauteur domaine
real e=.2;   //epaisseur plaque
real L0=1.;   //longueur
int n=25;    //nbre de points
    
// definition des cotes Maillage 
border b(t=0,1) { x= t*L0;    y = 0      ; };
border d(t=0,1) { x= L0;      y = h0 * t ; };
border h(t=1,0) { x= L0*(t);  y = h0     ; };
border g(t=1,0) { x= 0;       y = h0 * t ; };

mesh Th= buildmesh(b(n)+d(n)+h(n)+g(n));
plot(Th,wait=0,ps="pla.eps");

border hs(t=0,1) { x= L0*(1-t);   y =0;          };
border gs(t=0,1) { x= 0;          y = e * (-t);  };
border bs(t=0,1) { x= t*L0;       y = -e      ;  };
border ds(t=0,1) { x= L0;         y = -e*(1- t) ;} ;

mesh Sh=buildmesh(bs(n)+ds(n*.2)+d(n/3)+h(n)+g(n/3)+gs(n*.2)+hs(-n));
plot(Sh,wait=0,ps="pla.eps");
cout << " mesh "<< endl;

fespace Vh2(Th,P2);
fespace Vh(Th,P2);

//fespace W0(Th,P0);
//fespace Vh2s(Sh,P2);
fespace Vhs(Sh,P2);
fespace W0s(Sh,P0);

Vh2 u2,v2,up1,up2;
Vh2 u1,v1; 
Vhs T,TT,Told;
Vh2 Ue,Ve;
Vh p=0,q; 


real reynolds=500; 
real ks=0.0001*1000;  //
real dt = 0.1;// 0.0005;
real alpha=1/dt;
real U0=1;
real dd1=0.6479/sqrt(reynolds);

Ue=U0*(x-L0/2);
Ve=-(y-dd1)*U0;
u1=Ue ;
u2=Ve;
up1=u1;
up2=u2;  
T=0;
Told=0;

real [int] viso=[ 0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1];
int mur=Sh(L0/2,-e/2).region,fluid=Sh(L0/2,h0/2).region;
W0s  nu=(region==fluid)/reynolds  +  ks*(region==mur);  

int i=0,iter=0, ii=1;
plot(nu,fill=1,wait=0);



problem  NS (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
             )
 
  + int2d(Th)( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 ) 
  + on(b,u1=0,u2=0)
  + on(h,u1=Ue,u2=Ve) 
  + on(g,p=0) 
  + on(d,p=0) 
   ;
                      

Vh2 uu1=u1*(region==fluid), uu2=u2*(region==fluid);     

problem Ture (T,TT,solver=Crout) =
    int2d(Sh)( T*TT/dt +
             nu* ( dx(T)*dx(TT) + dy(T)*dy(TT))  )
 + int2d(Sh) ( -(1./dt)*convect([uu1,uu2],-dt,Told)*TT  
            )  
 + on(h,T=0)
 + on(bs,T=1)
 ;
      
/*--------------------------------------------------------------------------------*/

int nbiter = 100000,it=0;
real t=0;
real errl2=1;  

 
real bi=(e/ks/sqrt(reynolds));
cout << "Reynolds= "<<  reynolds <<" ks = "<<  ks << " b " << bi << " Tp=" << 1/(1+b)<< endl;
 

 alpha=1./dt; 
  
  
   	
while((errl2>5e-5)) //for (iter=0;iter<=nbiter;iter++)
{
  for (i=0;i<=4;i++)
   {
     up1=u1;
     up2=u2; 
     uu1=u1*(region==fluid); 
     uu2=u2*(region==fluid);     

     Told=T;
     NS;
     
      
     Ture;
     t=t+dt;
     
      errl2 = sqrt(int2d(Th)( (u1-up1)^2 + (u2-up2)^2) );
      errl2=errl2+sqrt(int2d(Sh)( (T -Told)^2) );
      cout << "t = "<< t <<" ---------------- "<< iter <<" +++++++++++++++ errl2 " << errl2<< endl;

     
     if ( !(i % 4))  {
     	  it++; 
     plot(T,cmm="T  t="+t,T,fill=1 ,ps="FILM/T"+it+".eps");
   //  plot(u1,cmm="U  t="+t,fill=1 );

    
real x,y;
int nn=100;

 { ofstream gnu("NplotT.gp"); 
gnu << "# vals" << endl; 
for (int i=0;i<nn;i++) 
{  x = L0/2;  y=-e+i*(h0+e)/nn; 
gnu<< y <<" "<<  T(x,y)    << endl; }
 }  
 
{ ofstream gnu("NplotT"+ii+".gp"); 
gnu << "# vals" << endl; 
for (int i=0;i<nn;i++) 
{  x = L0/2;  y=-e+i*(h0+e)/nn; 
gnu<< y <<" "<<  T(x,y)    << endl; }
 } 
 
 
{ ofstream gnu("NplotUT.gp"); 
gnu << "# vals x y u v T" << endl; 
for (int i=0;i<nn;i++) 
{ for (int j=0;j<nn;j++)
  {x =(i*L0)/nn;  
   y =-e+j*(h0+e)/nn; 
gnu<< (x-L0/2) <<"  " << y << " " << uu1(x,y) <<" "<< uu2(x,y)<< " " << T(x,y)  << endl; }
gnu<< endl;}}
}
   


// plot(coef=.02,cmm=" [u1,u2] et T  t="+t,fill=1,T,[uu1,uu2],viso=viso,ps="FILM/u1u2T"+iter+".eps");
   
cout << "CPU " << clock()-s0 << "s "<< endl; 
cout << "Reynolds= "<<  reynolds <<" ks = "<<  ks << " b " << b << " Tp=" << 1/(1+b)<< endl;
cout << "~~~~~~~~~~~~~~~~~~   R=" << reynolds << " " <<i << endl;     
   } 

     
//  Erreur 
 
if ((t>1.5) &&((errl2 ) < (1e-8))) break;
 

 }   
 
 
 
cout << "CPU " << clock()-s0 << "s " << endl;     

// run
// freefem++  poindarretkimarche.edp


/*
Pour gnuplot....


set xlabel "x"
set ylabel "y"
UU=4
p[:][0:]'NplotUT.gp' u ($1):2:($3/UU):($4/UU) every 4:4:4:4 t 'vitesse' w vect ,0 not w l linec -1 , 1/500**.5 t'1/sqrt(Re)'

#p[:.5][]'NplotUT.gp' u 1:2:($3/UU):($4*0/UU)  w vect  , 1/500**.5,0 not w l -1

#p[][]'NplotUT.gp' u ($2):4,0 not w l -1 , 1/500**.5


set xlabel "u(x,y)/x"
set ylabel "y"
#p[0:2][0:.5]'NplotUT.gp' u ($1*0):2:($3/$1):($4*0) every 2:2:2:2 t 'vitesse auto'w vect ,\
   0 not w l linec -1 , 1/500**.5 t'1/sqrt(Re)'


set xlabel "T(x,y)"
set ylabel "y"
p[0:][:.5]'NplotUT.gp' u 5:2 every 2:2:2:2 t 'temperature auto' ,\
   0 not w l linec -1 , 1/500**.5 t'1/sqrt(Re)'

*/