string dir="/tmp/xxx";
exec("mkdir "+dir);
int iter = 0; // depart si  diff=0  => depart fichier iter 
real dt=0.1,tfinal=1;
real tolNewton =1e-3; 
real T0 = 0;
real dTh = 0.05; 
real Tr = -dTh , Tl = 1; 
func  Tinit = (x < 0.0001 ) ? Tl : Tr ; // Solid .. 
real nus = 1e8, nuf = 1;

real Ra = 3.27e5;
real Pr = 56.2; 
real Ste = 0.045; 
real kf =1; 
real ks = 1; 

real ktanh = 50; //  parametre de regularisation de S 50 => couche 1/50 en temp.
real tanhshift  = 0  - T0*ktanh; 
real T0s = tanhshift/ktanh;  
real  errh = 0.2; //  Niveau erreur adapatation
real hmax = 0.1; // 
int n = 20;
real xl=0,xr=1;
mesh Th=square(n,n,[xl+x^1.3*(xr-xl),y]);
plot(Th,wait=1);
macro Save(dir,v,Th,i,temps)
{
    ofstream fv(dir+"/v-"+i);
     ofstream ff(dir+"/t-"+i);
    savemesh(Th,dir+"/Th-"+i);
    ff << temps << endl; 
    fv <<  v[0][];
    
}
//
macro Read(dir,v,Th,i,temps)
{
    ifstream fv(dir+"/v-"+i);
        ifstream ff(dir+"/t-"+i);
    Th=readmesh(dir+"/Th-"+i);    
    v = [0,0,0,0]; /* resize v */
    ff >> temps ;
    fv >>  v[0][];
    cout << " -- Read " << " temps = " << temps << " itera =" << i << endl; 
}
// 


fespace Wh(Th,[P2,P2,P1,P1]);
fespace Vh(Th,P2);
fespace Ph(Th,P1); 
real[int] viso(20);
for(int i=0;i<5;i++)
  viso(i) = Tr +i* (T0-Tr)/4.;
  
for(int i=6;i<viso.n;i++) 
  viso(i) = T0 + i*(Tl - T0) / (viso.n-7);
  viso(5)=T0-dTh/10;
  viso(6)=T0+dTh/10;
cout << " viso =" << viso << endl;    
real eps=1e-9;

Wh [u1,u2,p,T];
Wh [u1w,u2w,pw,Tw];// increment ...  
Wh [u1p,u2p,pp,Tp];

Wh [v1,v2,q,TT]; 



real cdt=1./dt;
 


real TCF=T0,DTMm=dTh*2;
real TM=TCF+DTMm/2,Tm=TCF-DTMm/2;
macro Phi(T) ((min(TM,max(Tm,T))-Tm) /DTMm)// 
macro dPhi(t) ( (t < Tm)  ? 0. : ( (t > TM) ? 0. : 1./DTMm ))  //
cout << Phi(1.) << " " << Phi(-1.) << " "  << endl ;


macro S(T) ((1+tanh(T*ktanh+tanhshift))/Ste/2.)//
macro dS(T) ((cosh(T*ktanh+tanhshift)^-2)*ktanh/Ste/2.)//



Vh kT ;
[u1p,u2p,pp,Tp]=[0,0,0,Tinit];

macro grad(u)[dx(u),dy(u)]//
macro Grad(u1,u2)[dx(u1),dy(u1),dx(u2),dy(u2)]//
macro div(u1,u2) (dx(u1)+dy(u2)) //
macro UgradV(u1,u2,v1,v2,T) [ [u1,u2]'*[dx(v1),dy(v1)] , [u1,u2]'*[dx(v2),dy(v2)], [u1,u2]'*[dx(T),dy(T)] ]// 


real fT = 0; 
real cmT = - Ra/Pr; 


func nuT =  Phi(T)*nuf + (1-Phi(T))* nus ;
func dnuT =  dPhi(T)*nuf  -dPhi(T)* nus ;



Vh nu , dnu;

problem BoussinesqNL([u1w,u2w,pw,Tw],[v1,v2,q,TT])
= int2d(Th) (
     [u1w,u2w,Tw]'*[v1,v2,TT]*cdt
     + UgradV(u1,u2,u1w,u2w,Tw)' * [v1,v2,TT]
     + UgradV(u1w,u2w,u1,u2,T)' * [v1,v2,TT]
     + ( Grad(u1w,u2w)'*Grad(v1,v2)) * nu 
     + ( Grad(u1,u2)'*Grad(v1,v2)) * dnu* Tw   
     + cmT*Tw*v2 
     + grad(Tw)'*grad(TT)*kT  
     - div(u1w,u2w)*q -div(v1,v2)*pw - eps*pw*q 
     + dS(T)*Tw*TT*cdt 
     )
   - int2d(Th)(
     [u1,u2,T]'*[v1,v2,TT]*cdt
     + UgradV(u1,u2,u1,u2,T)' * [v1,v2,TT]
     + ( Grad(u1,u2)'*Grad(v1,v2)) * nu 
     + cmT*T*v2 
     +  grad(T)'*grad(TT)*kT  
     - div(u1,u2)*q -div(v1,v2)*p 
     - eps*p*q  
     +  S(T)*TT*cdt 
     - [u1p,u2p,Tp]' *[v1,v2,TT]*cdt 
     -  S(Tp)*cdt*TT
      ) 
   + on(1,2,3,4, u1w=0,u2w=0) + on(2,Tw=0) + on(4,Tw=0) ;

u1[]=u1p[];
kT  =  1. / Pr ; 
real temps = 0;
if(iter>0) 
{
 Read(dir,[u1,u2,p,T],Th,iter,temps)
 real Umax=2; 
  [u1w,u2w,pw,Tw]=[u1w,u2w,pw,Tw];// increment ...  
 [u1p,u2p,pp,Tp]=[u1p,u2p,pp,Tp];
 [v1,v2,q,TT]=[v1,v2,q,TT]; 

 plot([u1,u2], Tp, wait=1,coef=0.1/Umax,value=0,viso=viso,cmm=iter+" U max = " + Umax+ " T= "+ temps, dim=2 );
 iter++; 
}
 
for( iter; iter<1000; iter++ )
{
 temps += dt;  
 
 u1p[]=u1[]; 


 { 
  Ph ph = S(T),pph=S(Tp);
  Th= adaptmesh(Th,T,Tp,ph,pph,[u1,u2],err=errh,hmax=hmax,hmin=hmax/100,ratio = 1.2);
//  plot(Th,phase,wait=1);
 [u1,u2,p,T]=[u1,u2,p,T];
 [u1w,u2w,pw,Tw]=[u1w,u2w,pw,Tw];// increment ...  
 [u1p,u2p,pp,Tp]=[u1p,u2p,pp,Tp];
 [v1,v2,q,TT]=[v1,v2,q,TT]; 
 } 
 real err=1e100,errp ;
 for(int kk=0;kk<2;++kk)
 {
   nu = nuT;
   dnu = 0; //dnuT;
   err = 1e100;
   
 for(int niter=0;niter<10; ++ niter)
 {
   errp = err;
   BoussinesqNL; 
   cout << niter << " err NL " << u1w[].linfty << endl; 
   u1[] -= u1w[];	
    err = u1w[].linfty;
   if( err > errp*100) break; 
   if(u1w[].linfty< tolNewton) break; 
  
  // plot([u1,u2], Tp); 
 }
 if( err > tolNewton  ) exit(1); 

}
 Vh U=sqrt(u1*u1+u2*u2); 
 real Umax = U[].max;
 U = U *( Phi(T) <0.0001) ;
 real Umaxs = U[].max;

 cout << "temps =" << temps << " u Max = " << Umax << " u max Sol " << Umaxs << endl;

if(iter%1==0) 
plot([u1,u2], Tp, wait=0,coef=0.05/Umax,value=0,viso=viso,cmm=iter+" U max = " + Umax+ " T= "+ temps, dim=2 );
if(iter%20==0) Save(dir,[u1,u2,p,T],Th,iter,temps)
if( temps > 80) break; 
}