Strange pressure when solving navier stokes using newton iteration

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.

You mention that you are using PETSc, but the code you copy/pasted does not use PETSc. What solver are you using?

Both, but the reslut was almost the same. Here is the code:

load "PETSc"
macro dimension()2//
include "macro_ddm.idp"

verbosity=0;
real dt = 0.1;
real T = 40;
int saveEach = 0.5/dt;
int[int] orderOut = [1, 1, 1, 1, 1];
string outputFileName = "output/output-2d.vtu";



macro div(u) (dx(u#x) + dy(u#y))//
macro grad(u) [dx(u), dy(u)]//
macro Grad(u) [dx(u#x),dy(u#x),dx(u#y),dy(u#y)]// EOM
macro UgradV(u, v)[[u#x, u#y]' * [dx(v#x), dy(v#x)],
                   [u#x, u#y]' * [dx(v#y), dy(v#y)]]//
				   
real cpu=clock();
//mesh Mesh = readmesh(meshFileName);
real h=0.04;
real L=3.2;
real H=5*h;
real l=4*H;
int[int] ll=[3,2,5,1];
// label 1  in
//       2  out 
//       3  down wall
//       4  step h
//       5  up wall
//       6  step v

mesh Mesh=square(400,60,[x*(l+L),y*H],label=ll);
Mesh=trunc(Mesh, (x>l) | (y >(H-h)),label=4); 

func uin=(H-y)*(y-(H-h))/square((H-(H-h))/2.);
real re=5000;
real nu = (H-h)/re;

if (mpirank == 0)
  cout << "Number of Elements: " + Mesh.nt << endl;

func PkVector = [P2, P2, P1];

buildDmesh(Mesh);
Mat J;
{
  macro def(i)[i, i#B, i#C]//
  macro init(i)[i, i, i]//
  createMat(Mesh, J, PkVector)
}

fespace SpaceVector(Mesh, PkVector);
SpaceVector [ucx, ucy, pc] = [1, 0, 0];
SpaceVector [uckx, ucky, pck];
SpaceVector [ucxold, ucyold, pcold];

if (mpirank == 0)
  cout << "Finite Element DOF (in each partition): " + SpaceVector.ndof << endl;

varf vRes([ux, uy, p], [vx, vy, q]) = int2d(Mesh)(
      (1/dt*[uckx, ucky]' * [vx, vy]) - (1/dt*[ucxold, ucyold]' * [vx, vy]) +
      nu * (Grad(uc)'*Grad(v))
    + UgradV(uc, uc)' * [vx, vy]
    - pc * div(v) - div(uc) * q)
    + on(3,4,5, ux = ucx-0, uy = ucy-0)
    + on(1, ux = 0, uy = 0);

varf vJ([ux, uy, p], [vx, vy, q]) = int2d(Mesh)(
      (1/dt*[ux, uy]' * [vx, vy]) +
      (UgradV(uc, u) + UgradV(u, uc))' * [vx, vy] +
      nu * (Grad(u)'*Grad(v))
    - p * div(v) - div(u) * q)
    + on(3,4,5, ux = ucx-0, uy = ucy-0)
    + on(1, ux = 0, uy = 0);

set(J, sparams = "-pc_type lu");

func real[int] funcRes(real[int]& inPETSc) {
    ChangeNumbering(J, ucx[], inPETSc, inverse = true, exchange = true);
    [uckx, ucky, pck]=[ucx, ucy, pc];
    real[int] out(SpaceVector.ndof);
    out = vRes(0, SpaceVector, tgv = -1);
    ChangeNumbering(J, out, inPETSc);
    return inPETSc;
}

func int funcJ(real[int]& inPETSc) {
    ChangeNumbering(J, ucx[], inPETSc, inverse = true, exchange = true);
    J = vJ(SpaceVector, SpaceVector, tgv = -1);
    return 0;
}

for(int i = 0; i < T/dt; i++)
{
  if (mpirank == 0)
    cout << "time step: " << i << endl;

    [uckx, ucky, pck] = [ucx, ucy, pc];

  real[int] xPETSc;
  ChangeNumbering(J, ucx[], xPETSc);
  SNESSolve(J, funcJ, funcRes, xPETSc, sparams = "-snes_monitor ");
  ChangeNumbering(J, ucx[], xPETSc, inverse = true, exchange = true);

  [ucxold, ucyold, pcold]=[ucx, ucy, pc];

  if (i % saveEach == 0)
    savevtk(outputFileName, Mesh, ucx, ucy,[ucx,ucy,0], pc, dataname="u v velocity p", order=orderOut,
            append = i ? true : false);
}
cout << " CPU time = " << clock()-cpu << endl;//CPU time = 1335.78

Your linearisation is incorrect, for example, this can’t be right in vRes, as it does not depend on the linearisation point:

    + on(1, ux = 0, uy = 0);

You mean

+ on(1, ux = ucx-uin, uy = ucy-0);

I initialized the velocity and pressure with [uin, 0,0], so ucx-uin would be 0 at boundary 1(inlet).

That looks better, yes.

but the pressure is still wrong

Do you get the proper pressure without SNESSolve()? Your linearization is likely wrong, but it’s difficult to debug without a reference working code.

I seem to have figured out what the problem is, I reduced the Reynolds number, and then the pressure was normal.
But why does the velocity seems normal when the Re number is very large, while the pressure doesn’t?