I am a beginner in finite element methods and I am using FreeFem++to solve parabolic equations. However, I encountered a problem: when I solved the Posion equation without a time term, the refinement grid error and error order were the same as the theory, but when I added a time term to the equation to solve the parabolic equation, I refined the grid, and the error and error order were clearly independent of the refinement grid. I asked this question a while ago because the exam did not continue, but this issue is very important to me and I cannot solve it. Therefore, I can only seek help from teachers or classmates to understand where the mistake lies.
Thank you very much for your help.
This is the error result I obtained
This is the code I received after being reminded by the teacher
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 error = cerr;
func U=exp(-t)*sin(pi*x)*sin(pi*y);
func Ux=pi*cos(pi*x)*sin(pi*y)*exp(-t);
func Uy=pi*sin(pi*x)*cos(pi*y)*exp(-t);
func f=(2*pi^2-1)*exp(-t)*sin(pi*x)*sin(pi*y);
func u0= sin(pi*x)*sin(pi*y);
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)(u*v/dt+dx(u)*dx(v)+dy(u)*dy(v) )
-int2d(Th)( uold*v/dt)
-int2d(Th)( f*v )
+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;
}
}