convectchacon-3d.edp

load "msh3" 
load "mat_psi" 
border a(t=0, 2*pi)     {    x = cos(t);    y = sin(t);  }; 
mesh th2 = buildmesh(a(70)); 
mesh3 th=  buildlayers(th2,5,zbound=[0,1]); 
fespace Vh(th,P1); 

Vh vh,vo,u1 = y, u2 = -x, v = exp(-10*((x-0.3)^2 +(y-0.3)^2)), u3=0, rhs =0; 
Vh wh=0,Bwh; 
real dt = 0.001,t=0, tmax=3.14; 
int kk=0; 
varf aa(v,vh) = int3d(th)(v*vh/dt) ; 
matrix  AA= aa(Vh,Vh,solver=UMFPACK); 

problem  A(v,vh) = int3d(th)(v*vh/dt)   - int3d(th)(vo*vh/dt) + rhs[]; 

real[int] viso=[-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,10,100,1000]; 
verbosity = 4; 
for ( t=0; t< tmax ; t+=dt) 
{ 
  vo=v; 
  matrix B; 
  MatUpWind0(B,th,vo,[u1,u2,u3]); 
//  cout << B << endl; 
  rhs[] = AA* vo[]; 
  rhs[] += B* vo[] ; 
  v[] = AA^-1*rhs[]; 
 // wh[][50]=1; 
  //Bwh[] = B*wh[]; 
 // plot(wh,Bwh,wait=1); 
  plot(v,fill=1,cmm="convection: t="+t + ", min=" + v[].min + ", max=" +  v[].max,viso=viso,wait=0); 
}; 
plot(v,wait=1,value=1,fill=1);