Dear all,
I am a beginner in finite element methods, and I am using FreeFem++ to solve parabolic equations. However, I have encountered a problem: as I refine the mesh, the error in my calculation results actually increases. Could you please advise me on where my calculation steps might have gone wrong?
Many thanks.
Regards.
This is my code:
int jie=3;
int NT=20;
real T=0.1;
real dt=T/NT;
macro div(u1,u2) (dx(u1)+dy(u2))
real LU0,HU0;
real [int] LU(jie),HU(jie),cpu(jie),hmax(jie);
real t=0;
func U=exp(-t)sin(pix)sin(piy);
func Ux=picos(pix)sin(piy)exp(-t);
func Uy=pisin(pix)cos(piy)exp(-t);
func f=(2pi^2-1)exp(-t)sin(pix)sin(piy);
func u0= sin(pix)sin(piy);
for(int n=0;n<jie;n++){
int NN=2^(n+3);
mesh Th=square(NN,NN);
hmax[n]=1./NN;
cpu[n]=clock();
fespace Vh(Th,P1);
Vh u=u0,v,uold;
solve Parabolic(u,v)
=int2d(Th)(uv/dt+dx(u)dx(v)+dy(u)dy(v) )
-int2d(Th)( uoldv/dt)
-int2d(Th)( fv )
+on(1,2,3,4,u=0)
;
ofstream ff(“Parabolic.dat”);
for(real t = dt; t < T; t += dt){
uold = u;
Parabolic;
plot(u,fill=1,value=1,cmm=“h^-1=”+1./hmax[n]);
}
Vh QiU=U;
Vh ThU=U;
LU0=int2d(Th)( (u-QiU)^2 );
HU0=int2d(Th)( (dx(u)-Ux)^2+(dy(u)-Uy)^2 );
cpu[n]=clock()-cpu[n];
LU[n]=sqrt(LU0);
HU[n]=sqrt(HU0);
error<<“h^-1=”<<1./hmax[n]<<" “<<LU[n]<<” “<<cpu[n]+“s”<<endl;
}
error<<“convergence rate”<<endl;
error<<“LU”<<” “<<“HU”<<endl;
for(int n=0;n<(jie-1);n++){
error<<log(LU[n]/LU[n+1])/log(hmax[n]/hmax[n+1])<<” “<<log(HU[n]/HU[n+1])/log(hmax[n]/hmax[n+1])<<endl;
}
for(int n=0;n<=(jie-1);n++){
if(n==0){
error.scientific<<”"+1./hmax[n]+"\\times"+1./hmax[n]+"“<<”& “<<LU[n]<<”& &“<<HU[n]<<”& \\“<<endl;
} else {
error.scientific<<”"+1./hmax[n]+"\\times"+1./hmax[n]+"“<<”& “<<LU[n]<<”& “<<log(LU[n-1]/LU[n])/log(hmax[n-1]/hmax[n])<<”& \\"<<endl;
}
}