convect-apt.edp

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

// This a the rotating hill problem with one turn.
// First 1/2 turn is a convection equation and second 1/2 a convection diffusion
int kt=6;
int kloop=5;
int nbadap=5;

bool inq=1;
real tol=0.05;
real tol2=1e-4;

border a(t=0, 2*pi)     {    x = cos(t);    y = sin(t);  }; // the unit circle
mesh th = buildmesh(a(70));                                 // triangulates the disk
fespace Vh(th,P1);

Vh u1 = y, u2 = -x;           // rotation velocity
Vh m11=0,m22=0,m12=0;         // to store the metric field
Vh vT;                        // to save the initial in big loop
int i;
Vh vv,vo,vp; // working  Finite element function 



func rhill=sqrt((x-0.3)^2 +(y-0.3)^2);
func  hill = 1-tanh(30*(rhill-0.2));
//func real hill(real r2)  { return 1-tanh(100*(r2-0.04));}

Vh v = hill;//((x-0.3)^2 +(y-0.3)^2);     // initial condition

plot(v);
vp=0;

for (int i=0;i<nbadap;i++)
 {
  th=adaptmesh(th,v,err=tol,inquire=inq);
  v = hill;//((x-0.3)^2 +(y-0.3)^2);
        real errl2= sqrt(int2d(th)(square(vp-v))) ;
  vp=vp-v;      
  cout << " interation adp =" << i << "  err l2 = " << errl2 
       << " diff min = " << vp[].min << " max =" << vp[].max << endl
       << " --------------- " << endl; 
  vp=v;
  if (errl2  < tol2) break;
 } 
plot(th,wait=1);

real dt = 0.17,t=0;                                                 // time step

vT[]=v[];

real error;
for ( i=0; i< 20/kt ; i++) {
    real T=t; 
    vp=0;
    for (int k=0;k<kloop;k++)
     {
       t=T; //  restart
       v=vT; // interpolation
       m11=0;m22=0;m12=0; // reset metric   
       adaptmesh(th,v,err=tol,metric=[m11[],m12[],m22[]],nomeshgeneration= true );  // warning change the order in version 1.28
       
       cout << "m11  = " << m11[].min << " " << m11[].max << endl;
       cout << "m22  = " << m22[].min << " " << m22[].max << endl;
       cout << "m12  = " << m12[].min << " " << m12[].max << endl;
        
       for (int j=0;j<kt;j++)
        {
          t = t+dt;
          vo[]=v[];
          v=convect([u1,u2],-dt,vo);                          // convect v by u1,u2, dt seconds, results in f
          plot(v,cmm="convection: t="+t + ", min=" + v[].min + ", max=" +  v[].max,wait=0);
          adaptmesh(th,v,err=tol,metric=[m11[],m12[],m22[]],nomeshgeneration= true);
        }
        th=adaptmesh(th,v,err=tol,metric=[m11[],m12[],m22[]]);
        vo=0; //
        plot(th,wait=1,cmm="k=" + k + ", t= " + t + " ,i= " + i );  
        real errl2= sqrt(int2d(th)(square(vp-v))) ;
        cout << " interation " << k << "  err l2 = " << errl2 << " --------------- " <<endl; 
        vp=v;
	error=errl2;
        if (errl2  < tol2) break;
     }
     vT=v; 
};


plot(th,v,wait=1);