convectchacon.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

load "mat_psi"
border a(t=0, 2*pi)     {    x = cos(t);    y = sin(t);  };
mesh th = buildmesh(a(70));
fespace Vh(th,P1);

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

problem  A(v,vh) = int2d(th)(v*vh/dt)   - int2d(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]);
//  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);