# Problem with 3d incompressible NS equation

**URL:** <https://community.freefem.org/t/problem-with-3d-incompressible-ns-equation/3434>\
**Category:** General Discussion\
**Created:** [August 12, 2024, 1:57am UTC](https://community.freefem.org/t/problem-with-3d-incompressible-ns-equation/3434 "2024-08-12T01:57:51Z")\
**Posts on this page:** 4\
**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:** [August 12, 2024, 1:57am UTC](https://community.freefem.org/t/problem-with-3d-incompressible-ns-equation/3434/1 "2024-08-12T01:57:51Z")

</div>

Hello,  
I am trying to write Freefem++ code for the unsteady Navier Stokes equations in 3D using Newton’s method. I found a similar example in the FreeFem hpddm examples. I changed it into 3d version.

The problem is when I set the MaxLayer = 5, everything goes fine and the result looks reasonable. But when I change the MaxLayer bigger than 20, the solver didn’t converge. Something wrong with the parameter in SNESSolve?

The code is as follows:  
///////////////////////////////////////////////////////////////////////////////////  
load “PETSc”  
macro dimension()3//  
include “macro\_ddm.idp”

verbosity=0;  
real dt = 0.001;  
real T = 200;  
int saveEach = 0.005/dt;  
int[int] orderOut = [1, 1, 1, 1, 1, 1];  
string outputFileName = “output/sub-3d.vtu”;  
/\*  
int Wall = 1;  
int Inlet = 3;  
int Outlet = 2;  
int InnerWall = 4;  
\*/  
real uin = 0.3;

real nu = 0.05;

macro velu(u) [u#x, u#y,u#z]//  
macro div(u) (dx(u#x) + dy(u#y) + dz(u#z))//  
macro grad(u) [dx(u), dy(u), dz(u)]//  
macro Grad(u) [dx(u#x),dy(u#x),dz(u#x),dx(u#y),dy(u#y),dz(u#y),dx(u#z),dy(u#z),dz(u#z)]// EOM  
macro UgradV(u,v) [[u#x,u#y,u#z]‘_[dx(v#x),dy(v#x),dz(v#x)] ,  
[u#x,u#y,u#z]'_[dx(v#y),dy(v#y),dz(v#y)] ,  
[u#x,u#y,u#z]’\*[dx(v#z),dy(v#z),dz(v#z)] ]//

mesh3 Mesh;  
{  
real r=0.0475;  
real pitch=0.063;  
real gap=2\*(0.063-0.0475);  
int nn=10;

if (mpirank == 0)  
{  
border C01(t=0, pi/2){x=-pitch+r_cos(t); y=-pitch+r_sin(t); label=1;}  
border C02(t=0, gap){x=-pitch; y=-pitch+r+t; label=2;}  
border C03(t=0, pi/2){x=-pitch+r_sin(t); y=pitch-r_cos(t); label=3;}  
border C04(t=0, gap){x=-pitch+r+t; y=pitch; label=4;}  
border C05(t=0, pi/2){x=pitch-r_cos(t); y=pitch-r_sin(t); label=5;}  
border C06(t=0, gap){x=pitch; y=pitch-r-t; label=6;}  
border C07(t=0, pi/2){x=pitch-r_sin(t); y=-pitch+r_cos(t); label=7;}  
border C08(t=0, gap){x=pitch-r-t; y=-pitch; label=8;}  
mesh Th = buildmesh(C01(-20)+C02(-10)+C03(-20)+C04(-10)+C05(-20)+C06(-10)+C07(-20)+C08(-10)  
);  
int[int] rup=[0,10], rdown=[0,9], rmid=[1,1,2,2,3,3,4,4];  
int[int] r1 = [0, 11];  
func zmin = 0.;  
func zmax = 0.5;  
int MaxLayer = 5;  
Mesh = buildlayers(Th, MaxLayer, zbound=[zmin, zmax], region=r1,labelmid=rmid,labelup = rup,labeldown = rdown);  
}  
broadcast(processor(0), Mesh);  
}

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

func PkVector = [P2, P2, P2, P1];

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

fespace SpaceVector(Mesh, PkVector);  
SpaceVector [ucx, ucy, ucz, pc] = [0,0,uin, 0];  
SpaceVector [uckx, ucky, uckz, pck]= [0,0,uin, 0];  
SpaceVector [ucoldx, ucoldy, ucoldz, pcold]= [0,0,uin, 0];

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

varf vRes([ux, uy, uz, p], [vx, vy, vz, q]) = int3d(Mesh)(  
(1/dt\*velu(uck)’ \* velu(v))

- (1/dt\*velu(ucold)’ \* velu(v)) +  
nu \* Grad(uc)'\*Grad(v)
  - UgradV(uc, uc)’ \* velu(v)

  - pc \* div(v) - div(uc) \* q  
//+ 0.000000001_pc_q  
)

  - on(1, 3,5,7, ux = ucx-0, uy = ucy-0, uz = ucz-0)
  - on(2,6, ux = ucx-0)  
+on(4,8,uy=ucy-0)
  - on(9, ux = ucx-0, uy = ucy-0, uz = ucz-uin);

varf vJ([ux, uy, uz, p], [vx, vy, vz, q]) = int3d(Mesh)(  
(1/dt_velu(u)’ \* velu(v)) +  
(UgradV(uc, u) + UgradV(u, uc))’ \* velu(v) +  
nu \* Grad(u)’_ Grad(v)  
- p \* div(v) - div(u) \* q  
//+0.000000001_p_q  
)  
+ on(1, 3,5,7, ux = ucx-0, uy = ucy-0, uz = ucz-0)  
+ on(2,6, ux = ucx-0)  
+on(4,8,uy=ucy-0)  
+ on(9, ux = ucx-0, uy = ucy-0, uz = ucz-uin);

set(J, sparams = “-pc\_type lu”);

func real[int] funcRes(real[int]& inPETSc) {  
ChangeNumbering(J, ucx, inPETSc, inverse = true, exchange = true);  
[uckx, ucky, uckz, pck]=[ucx, ucy, ucz, 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, uckz, pck] = [ucx, ucy, ucz, pc];

real[int] xPETSc;  
ChangeNumbering(J, ucx, xPETSc);  
SNESSolve(J, funcJ, funcRes, xPETSc, sparams = “-snes\_monitor”);  
ChangeNumbering(J, ucx, xPETSc, inverse = true, exchange = true);

[ucoldx, ucoldy, ucoldz, pcold]=[ucx, ucy, ucz, pc];

if (i % saveEach == 0)  
savevtk(outputFileName, Mesh, ucx, ucy, ucz,[ucx,ucy,ucz], pc, dataname=“u v w velocity p”, order=orderOut,  
append = i ? true : false);  
}

---

<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:** [August 12, 2024, 3:08pm UTC](https://community.freefem.org/t/problem-with-3d-incompressible-ns-equation/3434/2 "2024-08-12T15:08:43Z")

</div>

You probably want to do continuation to make the nonlinear solver converge.

---

<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:** [August 13, 2024, 1:12am UTC](https://community.freefem.org/t/problem-with-3d-incompressible-ns-equation/3434/3 "2024-08-13T01:12:10Z")

</div>

Thx for your reply, you mean MaxLayer continuation? or Re ?

---

<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:** [August 13, 2024, 5:56am UTC](https://community.freefem.org/t/problem-with-3d-incompressible-ns-equation/3434/4 "2024-08-13T05:56:56Z")

</div>

Continuation in Reynolds.