# Strange pressure when solving navier stokes using newton iteration

**URL:** https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560
**Category:** General Discussion
**Created:** [October 23, 2024, 12:57pm UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560 "2024-10-23T12:57:24Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![localhost](https://avatars.discourse-cdn.com/v4/letter/l/53a042/32.png) [@localhost](https://community.freefem.org/u/localhost)
#### Post date: [October 23, 2024, 12:57pm UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/1 "2024-10-23T12:57:24Z")

</div>

I solve navier stokes problem using newton method. I used the code from ([FreeFem++ users | Main \>\> Newton1](https://www.um.es/freefem/ff++/pmwiki.php?n=Main.Newton1))  
I changed the code into transient problem as follows:

```auto
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.

 ![1](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/000c49f10a147d1b74efeafdfadd628d2a8c766c.png)  
 ![2.2](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b506487e5e5b71f8b7d416109846f11a2aa22149.png)  
 ![2](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/2/28837430dffdff171313f806950ee1a4397593c4.png)

I used PETSc and SNES to solve the navier stokes problem. But the pressure is still unreasonable.

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 24, 2024, 6:04am UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/2 "2024-10-24T06:04:27Z")

</div>

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

---

<div class="post-metadata">

### Author: ![localhost](https://avatars.discourse-cdn.com/v4/letter/l/53a042/32.png) [@localhost](https://community.freefem.org/u/localhost)
#### Post date: [October 24, 2024, 12:34pm UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/3 "2024-10-24T12:34:24Z")

</div>

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

```auto
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

```

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 24, 2024, 1:30pm UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/4 "2024-10-24T13:30:29Z")

</div>

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

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

```

---

<div class="post-metadata">

### Author: ![localhost](https://avatars.discourse-cdn.com/v4/letter/l/53a042/32.png) [@localhost](https://community.freefem.org/u/localhost)
#### Post date: [October 25, 2024, 2:31pm UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/5 "2024-10-25T14:31:54Z")

</div>

You mean

```auto
+ 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).

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 25, 2024, 2:38pm UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/6 "2024-10-25T14:38:26Z")

</div>

That looks better, yes.

---

<div class="post-metadata">

### Author: ![localhost](https://avatars.discourse-cdn.com/v4/letter/l/53a042/32.png) [@localhost](https://community.freefem.org/u/localhost)
#### Post date: [October 26, 2024, 3:09am UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/7 "2024-10-26T03:09:41Z")

</div>

but the pressure is still wrong

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 26, 2024, 5:07am UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/8 "2024-10-26T05:07:02Z")

</div>

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

---

<div class="post-metadata">

### Author: ![localhost](https://avatars.discourse-cdn.com/v4/letter/l/53a042/32.png) [@localhost](https://community.freefem.org/u/localhost)
#### Post date: [October 26, 2024, 11:31am UTC](https://community.freefem.org/t/strange-pressure-when-solving-navier-stokes-using-newton-iteration/3560/9 "2024-10-26T11:31:06Z")

</div>

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?