convect_dervieux.edp

bool fast=true; 
load "mat_dervieux";  // external module in C++ must be loaded 
border a(t=0, 2*pi){ x = cos(t); y = sin(t);  } 
mesh th = buildmesh(a(100)); 
fespace Vh(th,P1); 

Vh vh,vold,u1 = y, u2 = -x; 
Vh v = exp(-10*((x-0.3)^2 +(y-0.3)^2)), vWall=0, rhs =0; 

real dt = 0.025; 
// qf1pTlump means mass lumping is used 
problem  FVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt) 
                  - int2d(th,qft=qf1pTlump)(vold*vh/dt) 
      + int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh) 
+ rhs[] ; 


matrix A; 
MatUpWind1(A,th,vold,[u1,u2]); 
if(fast) 
  { 
    varf  vFVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt) 
      - int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh)      ; 
    real[int] rhs0=vFVM(0,Vh); 
    matrix M=vFVM(Vh,Vh,solver=CG); 
    A=-A+M; 

    for ( real t=0; t< pi ; t+=dt) 
      { 
	vold=v; 
        rhs[]=rhs0; 
	rhs[] += A * vold[] ; 
	v[]= M^-1*rhs[]; 
	plot(v,wait=0); 
      } 
  } 
else 
for ( real t=0; t< pi ; t+=dt){ 
  vold=v; 
  rhs[] = A * vold[] ; 
  FVM; 
  plot(v,wait=0); 
};