// example 7
// Thermal convection problem 
// Navier-Stokes equations + convection-diffusion equation
// COE-tutorail 2007, Atsushi Suzuki 2007/03/08
//
// Prandtl Number for water
real Pr=7.0;
real Ra=2.0e+3;
real dt=1.0e-3;
real alpha = 1.0 / dt;
real beta=1.0 / (dt * Pr);
real eps=0.01;
int i = 0;
int timestepmax = 5000;
real ptmp;

// initial thermal distirbution with perturbation
func thini = (1.0 - y) + (1.0 - 4.0*(x-0.5)*(x-0.5))*(1-y)*y*eps;

border ba(t=0,1.0){x=2.0*t;y=0.0;label=1;};
border bb(t=0,1.0){x=2.0;y=t;label=2;};
border bc(t=0,1.0){x=2.0-2.0*t;y=1.0;label=3;};
border bd(t=0,1.0){x=0.0;y=1.0-t;label=4;};

// unstructured mesh 
mesh Th=buildmesh(ba(51)+bb(24)+bc(50)+bd(25));


plot(Th,wait=1,value=true,coef=0.1);
//fespace Xh(Th,P2),Mh(Th,P1);
fespace Xh(Th,P1b),Mh(Th,P1);
Xh u1,u2,v1,v2,uu1,uu2,up1,up2,th0;
Mh p,q;
Mh th,psi,phi;

macro d11(u1)     dx(u1) //
macro d22(u2)     dy(u2) //
macro d12(u1,u2)  (dy(u1) + dx(u2))/2.0 //

real epslon = 1.0e-6;

// heating at the bottom on(1, th=1) and cooling at the top on(3,th=0)
problem ConvectDiff (th, psi, solver=Crout, init=i) =
	int2d(Th) (
           alpha * (th * psi) +
	   (dx(th)*dx(psi) + dy(th)*dy(psi))
        )
        - int2d(Th) (alpha * convect([uu1,uu2],-dt,th0)*psi)
        + on(1,th=1)
	+ on(3,th=0)
;

// slip boundary coniditons for velocity u.n = 0
problem NS([u1,u2,p], [v1,v2,q],solver=Crout,init=i) = 
  int2d(Th)( 
           beta * (u1*v1 + u2*v2)
         + 2.0 * (d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))
	 - dx(v1)*p - dy(v2)*p - dx(u1)*q - dy(u2)*q
    	 - p*q*epslon
	 ) 
- int2d(Th)(
	   beta * convect([up1,up2],-dt,up1)*v1
         + beta * convect([up1,up2],-dt,up2)*v2
	 + Ra * th0 *v2 
	 )
	+  on(1,3,u2=0)
	+  on(2,4,u1=0)
;
	
// stream line for visualization
problem streamlines(phi,psi,solver=Crout,init=i) =
     int2d(Th)(
	dx(phi)*dx(psi)+dy(phi)*dy(psi)
     )
     - int2d(Th)(
       (dx(u2)-dy(u1))*psi)
     +on(1,2,3,4,phi=0);

th = thini;
plot(th, wait=1,value=true,coef=0.1);

real area = int2d(Th)(1.0);
real meanp;

u1 = 0;
u2 = 0;
for (i = 0; i < timestepmax; i++) {
  th0 = th;
  up1=u1;
  up2=u2;
  NS;

//  meanp = int2d(Th)(p) / area;
//  p = p - meanp;

  uu1 = u1;
  uu2 = u2;
  ConvectDiff;

// show velocity and pressure at every 25th step
  if(i%25==0) {
     plot([u1,u2],th0, value=true,coef=.1);
//     streamlines;
//     plot(phi, value=true,coef=10.0);
  }
}