I solve navier stokes problem using newton method. I used the code from (FreeFem++ users | Main >> Newton1)
I changed the code into transient problem as follows:
load "iovtk"
real h=0.04;
real L=3.2;
real H=5*h;
real l=4*H;
bool dplot=1; // debug plot
bool adapt=0;
real T=40.,dt=0.1;
real alpha=1./dt;
string FileName;
int[int] ll=[3,2,5,1];
// label 1 in
// 2 out
// 3 down wall
// 4 step h
// 5 up wall
// 6 step v
border A1(t=0,1){x=0; y=H-t*h; label=1;};
border A2(t=0,1){x=t*l; y=H-h; label=4;};
border A3(t=0,1){x=l; y=(H-h)*(1-t); label=6;};
border A4(t=0,1){x=l+6*(H-h)*t; y=0; label=3;};
border A5(t=0,1){x=l+6*(H-h)+(L-6*(H-h))*t; y=0; label=3;};
border S1(t=0,1){x=l+L; y=t*H; label=2;};
border S2(t=0,1){x=l+6*(H-h)+(L-6*(H-h))*(1-t); y=H; label=5;};
border S3(t=0,1){x=l+6*(H-h)*(1-t); y=H; label=5;};
border S4(t=0,1){x=l*(1-t); y=H; label=5;};
mesh Th=buildmesh( A1(8)+A2(80)+A3(40)+A4(80)+A5(160)+S1(40)+S2(200)+S3(100)+S4(100) );
plot(Th);
func uin=(H-y)*(y-(H-h))/square((H-(H-h))/2.);
//
// Finite element spaces
//
fespace Xh(Th,P2);
fespace Mh(Th,P1);
fespace XXMh(Th,[P2,P2,P1]);
XXMh [u1,u2,p];
XXMh [v1,v2,q];
XXMh [up1,up2,pp]; //temp
XXMh [uc1,uc2,pc]; //previous time
//
// macros for the variational formulation
//
macro Grad(u) [dx(u#1),dy(u#1),dx(u#2),dy(u#2)]// EOM
macro div(u) (dx(u#1)+dy(u#2)) //EOM
macro UgradV(u,v) [ [u#1,u#2]'*[dx(v#1),dy(v#1)] , [u#1,u#2]'*[dx(v#2),dy(v#2)] ]//EOM
//real nu=1./100.;
real re=1000;
real nu = (H-h)/re;
Xh ms=hTriangle;
Xh CFL;
//cout << " mesh size = " << ms[].min << " " << ms[].max << " "<< ms[].sum/ms[].n << endl;
//
// Stokes problem
//
solve Stokes ([u1,u2,p],[v1,v2,q],solver=UMFPACK) =
int2d(Th)( nu*(Grad(u)'*Grad(v))
+ p*q*(0.00000001)
- p*div(v)-q*div(u)
)
+ on(1,u1=uin,u2=0)
+ on(3,4,5,6,u1=0,u2=0)
+ on(2,p=0)
;
CFL=u1/ms*dt;
cout << " CFL max = " << CFL[].max<<endl;
int[int] Order = [1,1,1,1];
string DataName = "u1 u2 vector p";
savevtk("stokes.vtu", Th, u1,u2,[u1,u2,0],p, dataname=DataName, order=Order);
int i=0;
//
// To build the matrix for the linear system in Newton method
//
// this is the matrix that multiply the increment, dF
//
varf vDNS ([u1,u2,p],[v1,v2,q]) =
int2d(Th)(alpha*(u1*v1+u2*v2)
+ nu * (Grad(u)'*Grad(v))
+ p*q*1e-8// stabilization term
- p*div(v)
- div(u)*q
+ UgradV(u,up)'*[v1,v2]
+ UgradV(up,u)'*[v1,v2]
)
+ on(1,3,4,5,6,u1=0,u2=0)
;
//
// To build right hand side for the linear system in Newton method
//
varf vNS ([u1,u2,p],[v1,v2,q]) =
int2d(Th)(-alpha*(up1*v1+up2*v2-uc1*v1-uc2*v2)
-nu * (Grad(up)'*Grad(v))
+ pp*q*1e-8// stabilization term
+ pp*div(v)
+ div(up)*q
- UgradV(up,up)'*[v1,v2]
)
+ on(1,3,4,5,6,u1=0,u2=0)
;
real lerr=0.01;
real finalerror;
for(real time=0;time<T;time+=dt)
{
//restore values of previous time step
uc1[]=u1[];
up1[]=u1[];
for(int step=0;step<(adapt?2:1) ;step++)
{
if(adapt)
{
Th=adaptmesh(Th,[u1,u2],p,abserror=1,cutoff=1e-5,err=lerr,nbvx=200000,hmin=0.002);
if(dplot) plot(Th,wait=1);
}
[u1,u2,p]=[u1,u2,p];
[up1,up2,pp]=[up1,up2,pp];
for (i=0;i<=20;i++)// Newton iterations until convergence for a given Rey
{
//up1[]=u1[];
real[int] b = vNS(0,XXMh); // Right hand side for the linear system
matrix Ans=vDNS(XXMh,XXMh);// Matrix for the linear system
set(Ans,solver=UMFPACK);// Set solver to matrix
real[int] w = Ans^-1*b; // Solve the system
up1[] += w; // Perform the update of the increment in both variables at the same time
cout << " iter = "<< i << " " << w.l2 << " rey = " << re << endl;
if(w.l2<1e-6)
{
finalerror=w.l2;
break;
}
// uu1=u1;uu2=u2;
};
CFL=u1/ms*dt;
//cout << " CFL max = " << CFL[].max<<endl;
if(dplot) plot(coef=0.2,cmm="Time="+time+" final error="+finalerror+"CFL[].max="+CFL[].max+ " [u1,u2] et p ",pp,[up1,up2]);
u1[]=up1[];
}
}
savevtk("final.vtu", Th, u1,u2,[u1,u2,0],p, dataname=DataName, order=Order);
Newton iteration converged. The velocity seems good.
Everything goes fine until I check the pressure. In the first time step the pressure is reasonable but in the second time step the pressure goes wrong.
I used PETSc and SNES to solve the navier stokes problem. But the pressure is still unreasonable.