Dear all,
I am trying to compute the L2 velocity rate for time dependent Navier-Stokes equations using P2-P1 finite element space. when I take Reynolds number =1, then it works well, but when I take Reynolds number= 10 or above 10, then the L2 error convergence rate for velocity decreases and it is not coming 3, according to theory. What should I do to recover the convergence rate?
Here I have attached the code for Reynolds number =10.
/*
- Incompressible Navier Stokes with Taylor-Hood Finite element
- NSE use Newton method for solution 1
- Sol x y
*/
//load “Element_P3”
verbosity = 0;
// Parameters
real nu = 1./10.;
int maxiter = 6; // Maximum iteration
real[int] L2errorV(maxiter);
real[int] H1errorV(maxiter);
real[int] L2errorP(maxiter);
real[int] cputime(maxiter); // for store cpu time
// Start mesh loop
for (int i=1; i<= maxiter; i++ ) {
real cpu = clock();
// mesh size
int NN = 2^(i); // mesh size
real dt = 1./NN^3;
real T = 1; // final time
real TN = T/dt; // Total number of time iteration
real eps = 1e-6; // newton error
real epsilon = 1e-6; // // penalty parameter
// Mesh
mesh Th = square(NN, NN, flags=0);
// Fespace
fespace Xh(Th, P2);
Xh u1, u2, du1, du2, v1, v2;
Xh up1, up2, uu1, uu2, psi, phi;
fespace Mh(Th, P1);
Mh p, q, pp, pp1, dp;
// Macro
macro div(u1, u2) (dx(u1) + dy(u2)) //
//macro grad(u1, u2) [dx(u1), dy(u2)] //
macro ugrad(u1, u2, v) (u1*dx(v) + u2*dy(v)) //
macro Ugrad(u1, u2, v1, v2) [ugrad(u1, u2, v1), ugrad(u1, u2, v2)] //
// Define exact solution as a func of x y and t
func real u1exact(real t)
{
return exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1);
}
func real u2exact(real t)
{
return -exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2;
}
func real pexact(real t)
{
return exp(t)*2*(x-y);
}
// Define right hand side for steady Stokes
Xh fL1 = - nu* (4*y*(2*y - 1)*(x - 1)^2*(y - 1) + 4*x^2*y*(2*y - 1)*(y - 1)
+ 8*x*y*(2*x - 2)*(2*y - 1)*(y - 1) + 8*x^2*(x - 1)^2*(y - 1) + 4*x^2*(2*y - 1)*(x - 1)^2 + 8*x^2*y*(x - 1)^2)
+ 2;
Xh fL2 = - nu* (- 8*y^2*(x - 1)*(y - 1)^2 - 4*y^2*(2*x - 1)*(y - 1)^2 - 8*x*y^2*(y - 1)^2
- 4*x*(2*x - 1)*(x - 1)*(y - 1)^2 - 4*x*y^2*(2*x - 1)*(x - 1) - 8*x*y*(2*x - 1)*(2*y - 2)*(x - 1))
-2;
// Right hand side for unsteady NSE
func real f1(real t)
{
return exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1) - nu* exp(t)*(4*y*(2*y - 1)*(x - 1)^2*(y - 1) + 4*x^2*y*(2*y - 1)*(y - 1)
+ 8*x*y*(2*x - 2)*(2*y - 1)*(y - 1) + 8*x^2*(x - 1)^2*(y - 1) + 4*x^2*(2*y - 1)*(x - 1)^2 + 8*x^2*y*(x - 1)^2)
+ 2*exp(t) + exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1) *exp(t)*4*x*(x-1)*(2*x-1)*y*(y-1)*(2*y-1)
-exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2* exp(t)* 2*x^2*(x-1)^2*((2*y-1)^2+2*y*(y-1));
}
func real f2(real t)
{
return -exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2 - nu* exp(t)*(- 8*y^2*(x - 1)*(y - 1)^2 - 4*y^2*(2*x - 1)*(y - 1)^2 - 8*x*y^2*(y - 1)^2
- 4*x*(2*x - 1)*(x - 1)*(y - 1)^2 - 4*x*y^2*(2*x - 1)*(x - 1) - 8*x*y*(2*x - 1)*(2*y - 2)*(x - 1))
-2*exp(t) + exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1) * (-exp(t)*2*y^2*(y-1)^2*((2*x-1)^2+2*x*(x-1)))
-exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2 * (-exp(t)*4*x*(x-1)*(2*x-1)*y*(y-1)*(2*y-1));
}
// Problem Stokes (with solve)
solve Stokes ([u1, u2, p], [v1, v2, q], solver=sparsesolver)
= int2d(Th)
(
nu*( dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
- epsilon* p * q // Penalty term
- p*div(v1, v2) - div(u1, u2)*q )
- int2d(Th)([fL1,fL2]'*[v1,v2])
+ on(1, 2, 3, 4, u1=0, u2=0);
// Update sol for NSE
uu1 = u1;
uu2 = u2;
pp1 = p;
//plot(coef=0.2, cmm="Stokes [u1, u2] and p" ,p, [uu1, uu2], wait=1);
//cout << "Stokes L2 Error= " << sqrt(int2d(Th)((u1-u1exact)^2 + (u2-u2exact)^2)) << endl;
// Navier Stokes Solve
// Time loop
for (int j = 1; j <= TN; j++) {
real err=0;
// Update right hand side at each time step
real t = j*dt;
// Newton method
for (int n = 0; n <= 20; n++) {
// Solve
solve NSE ([du1, du2, dp], [v1, v2, q], solver=sparsesolver)
= int2d(Th)(
1./dt*(du1*v1 + du2*v2)
+nu * (
dx(du1)*dx(v1) + dy(du1)*dy(v1)
+ dx(du2)*dx(v2) + dy(du2)*dy(v2)
)
+ Ugrad(du1, du2, u1, u2)'*[v1, v2]
+ Ugrad(u1, u2, du1, du2)'*[v1, v2]
- dp*div(v1, v2) - div(du1, du2)*q
- epsilon*dp * q
)
+ int2d(Th)(
1./dt*(u1*v1 + u2*v2)
+ nu * (
dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2)
)
+ Ugrad(u1, u2, u1, u2)'*[v1, v2]
- p*div(v1, v2) - div(u1, u2)*q
)
- int2d(Th)([f1(t),f2(t)]'*[v1,v2])
- int2d(Th) (
1./dt*(uu1*v1 + uu2*v2)
)
+ on(1, 2, 3, 4, du1=0, du2=0);
// Update
u1[] += du1[];
u2[] += du2[];
p[] += dp[];
real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
cout << "Mesh = " << NN << ", Time = " << j*dt << ", Newton iteration = " << n+1 << ", err = " << err << ", eps = " << eps << endl;
if(err < eps) break; //converge
}
// update NSE solution for next step
uu1 = u1;
uu2 = u2;
// plot
plot(coef=0.2, cmm="Time= "+t+", [u1, u2] and p" ,p, [u1, u2], wait=0);
//plot(coef=0.2, cmm="u1" ,u1, wait=1);
//plot(coef=0.2, cmm="u2" ,u2, wait=1);
//plot(coef=0.2, cmm="p" ,p, wait=1);
//plot(coef=0.2, cmm="pexact" ,pexact, p, wait=1);
}
// Exact sol at final time
Xh u1ex = u1exact(T);
Xh u2ex = u2exact(T);
Mh pex = pexact(T);
// error
L2errorV[i-1] = sqrt(int2d(Th)((u1-u1ex)^2 + (u2-u2ex)^2));
H1errorV[i-1] = sqrt(int2d(Th)( (dx(u1)-dx(u1ex))^2+(dy(u1)-dy(u1ex))^2 + (dx(u2)-dx(u2ex))^2+(dy(u2)-dy(u2ex))^2));
L2errorP[i-1] = sqrt(int2d(Th)((p-pex)^2));
//cout << "L2 Error= " << L2error << endl;
// cpu time
cputime[i-1] = clock()-cpu;
}
// Print error
//
cout << "L2 Error= " << L2errorV << endl;
cout << "L2 Error= " << H1errorV << endl;
cout << "L2 ErrorP= " << L2errorP << endl;
// calculte rate
real[int] rcvl2(maxiter-1);
real[int] rcvh1(maxiter-1);
real[int] rcpl2(maxiter-1);
for(int j=1; j<maxiter; j++)
{
rcvl2[j-1] = log(L2errorV[j-1]/L2errorV[j])/log(2.);
rcvh1[j-1] = log(H1errorV[j-1]/H1errorV[j])/log(2.);
rcpl2[j-1] = log(L2errorP[j-1]/L2errorP[j])/log(2.);
}
cout << "L2 Error rate vel= " << rcvl2 << endl;
cout << "H1 Error rate vel= " << rcvh1 << endl;
cout << "L2 Error rate pre= " << rcpl2 << endl;
cout << " CPU time = " << cputime << endl;