// resolution de NS avec THERMIQUE
// PYL oct 2005 / fev 07 / dec 17
// inspire de l'exemple NSP1P2.edp
//
// freefem++ thermiquekimarcheRB.pde 
// 
// http://www.lmm.jussieu.fr/%7Elagree/COURS/ENSTA/web/CLib/RayleighBenard/index.html
// 
// mettre l'oppose du nombre de Rayleight (-Rayleigh) dans le fichier JJ2.txt
// par exemple si -1800 est dans JJ2.txt, la configuration est stable 
//
//
exec("echo kimarche");


real s0=clock();
real h0=1;
real L0=5; 
int n=8;


 
// longueur L0, epaisseur h0, ici 5 sur 1
// 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=0,1) { x= 0;          y = h0 * (1-t);    };


// Oh le beau maillage
mesh Th = buildmesh(b(n*5)+d(n)+h(n*5)+g(n));  
mesh Thc= buildmesh(b(n*5/2)+d(n/2)+h(n*5/2)+g(n/2));  
plot(Th,wait=0,ps="cuve.eps");


//espaces EF 
fespace Vh2(Th,P2);
fespace Vh(Th,P1);
fespace Vhc(Th,P2);

// declaration des champs
Vh2 u2,v2,up1,up2;
Vh2 u1,v1; 
Vhc u1c,u2c,Tc;
Vh2 T,TT,Told;

// Reynolds est unitaire par construction des echelles de tems pet de vitesse
// l'information sera dan Rayleigh
real reynolds=1;

// initialisation
u1=0 ;
u2=0 ;
up1=u1;
up2=u2; 
Told=0;
T=0;

// initialisation
Vh p=0,q;
real J1=0,J2=-1;
real dt10;
real dt =0.005;
dt10=10*dt;
real alpha=1/dt;
real  nu=1./reynolds;   
int i=0,iter=0;

// mettre (-Rayleigh) dans JJ2.txt, le fichier est ensuite relu, on peut donc faire varier le Rayleigh en direct
{ ifstream f("JJ2.txt");   f>>J2;} 
 
// Formultaion variationnelle de Navier Stokes
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)(J2*Told*v2)
  + 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=0,u2=0) 
  + on(g,u1=0,u2=0) 
  + on(d,u1=0,u2=0)  ;


// Formultaion variationnelle de l'équation de la chaleur
problem NST (T,TT,solver=Crout,init=i) =
    int2d(Th)( T*TT/dt +
             nu * ( dx(T)*dx(TT) + dy(T)*dy(TT)) 
             )
 + int2d(Th) ( -(1./dt)*convect([u1,u2],-dt,Told)*TT 
            ) 
 + on(h,T=0)
 + on(b,T=1)
 ;


// pour tracer les lignes de courant 
Vh2 psi,phi; 

problem streamlines(psi,phi) =
  int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi)) 
+ int2d(Th)( -phi*(dy(u1)-dx(u2))) 
+ on(b,psi=0)
+ on(g,psi=0)
+ on(h,psi=0)
+ on(g,psi=0);
//-------------------------------------------------------

 
int nbiter = 100;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter) , cut=0.01;
real tol=0.5,coeftol = 0.25^(1./nbiter);
real t=0;
  

// boucle de calcul
for (iter=0;iter<nbiter;iter++)
{ //if(t>1)dt=dt10;;
  real errl2 = sqrt(int2d(Th)( (u1-up1)^2 + (u2-up2)^2  +(T -Told)^2) );
  cout << "================" <<
        "t = "<< t <<" ---------------- "<< iter <<" +++++++++++++++ errl2 " << errl2<< endl;
  alpha=1./dt;
  nu=1./reynolds; 
  for (i=0;i<100;i++)
   {
     up1=u1;
     up2=u2;
     NS;
     NST; 
     Told=T;
     t=t+dt;
     if ( !(i % 10))  {
     	real maxu=u1[].max+.00001;
     	Tc=T;u1c=u1/maxu;u2c=u2/maxu;
     plot(coef=.05,cmm=" [u1,u2] et T  t="+t,Tc,[u1c,u2c]);
     
// sauvegarde gnuplot     
real x,y;
n=100;
{ ofstream gnu("Nplot.gp"); 
gnu << "# vals" << endl; 
for (int i=0;i<=n;i++) 
{ x =(i*L0)/n*1.0;  y=h0/2; 
gnu<< x <<" "<<  T(x,y)  << " "<< u1(x,y)<<" "<< u2(x,y)<< endl; }
 } }
     cout << "CPU " << clock()-s0 << "s " << "  J2=" << J2 << " " <<i << endl;     
   } 


 
streamlines;
//plot(wait=0,psi,ps="psi.eps");


 //  Erreur 
 // adaptation de maillage
 /* Th=adaptmesh(
  Th,[dx(u1),dy(u1),dx(u1),dy(u2)],abserror=0,cutoff=cut,err=tol, 
  inquire=0,ratio=1.5,hmin=1./1000);
  plot(Th);
  dt = dt*coefdt;
  tol = tol *coeftol;
  cut = cut *coefcut;
*/


// relecture du Rayleigh 
// cela permet de faire un chauffage dynamique a la main
// faire apparaitr ou disparaitre les rouleaux
//{ ifstream f("JJ1.txt");   f>>J1;}
{ ifstream f("JJ2.txt");   f>>J2;}
 
 
 if ((t>15) &&((errl2 ) < (5e-5))) break;

 }  
streamlines;
plot(wait=0,psi,ps="psi.eps");
cout << "t = "<< t << endl; 
cout << "CPU " << clock()-s0 << "s " << endl;     



// 
// plus de details
// http://www.lmm.jussieu.fr/%7Elagree/COURS/ENSTA/web/CLib/RayleighBenard/index.html
//
//