# Problem solving Navier Stokes equations with Newton method

**URL:** <https://community.freefem.org/t/problem-solving-navier-stokes-equations-with-newton-method/3538>\
**Category:** General Discussion\
**Created:** [October 8, 2024, 10:10am UTC](https://community.freefem.org/t/problem-solving-navier-stokes-equations-with-newton-method/3538 "2024-10-08T10:10:05Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![SusannaM](https://avatars.discourse-cdn.com/v4/letter/s/f04885/32.png) [@SusannaM](https://community.freefem.org/u/SusannaM)\
**Post date:** [October 8, 2024, 10:10am UTC](https://community.freefem.org/t/problem-solving-navier-stokes-equations-with-newton-method/3538/1 "2024-10-08T10:10:06Z")

</div>

Hi, i am trying to solve NS equation with Newton method in a 3D cubic domain with non zero value of the y-component of the velocity at the inflow boundary.  
The problem I have is that the velocity is almost not updating at each iteration and each component of the final velocty has exactly the same values.

which could be the reason of this behaviour?  
here is the code I ame running:  
real L = 1; // Characteristic length of the cell (m)  
real mu = 1e-1; // Air viscosity at 300K  
real U = 0.6; // Characteristic velocity and inlet boundary velocity (Dirichlet condition NS)  
real absu = 0.; // Module of the velocity  
real Re = rhoa_U_L/mu;  
load “msh3”;  
load “iovtk”;  
load “medit”; // Per visualizzazione 3D con medit  
//Mesh size  
int meshSize = 20;

// Mesh generation  
mesh3 Th = cube(meshSize, meshSize, meshSize, [L_x,L_y,L\*z]);  
// Visualize the mesh using the plot command  
plot(Th,cmm = “3D Mesh Visualization”);

fespace Vh(Th, [P2, P2, P2]);  
Vh [ux,uy,uz] , [vx,vy,vz]; //Fluid velocity, test function  
Vh [ux1,uy1,uz1] , [dux,duy,duz]; //Fluid Velocity at the next iteration, variation between one iteration and another  
fespace Qh(Th, P1);  
Qh p, q; // pressure e test function  
Qh p1; // pressure at the next iteration  
Qh Q;  
// Stabilization term definitions  
Qh hK = hTriangle; // stabilization parameter based on the element size (hTriangle is usually a built-in function in FreeFem++ that returns the size of the element).  
Qh tauK; //will be used to store stabilization coefficients.  
//Macro NB: ‘\* rappresenta il prodotto scalare  
macro grad(v) [dx(v), dy(v), dz(v)] //  
macro div(vx,vy,vz) ( dx(vx) + dy(vy) + dz(vz) ) //  
macro ugradv(ux,uy,uz,vx,vy,vz) [[ux,uy,uz]’\*grad(vx) , [ux,uy,uz]‘\*grad(vy) , [ux,uy,uz]’\*grad(vz) ] //  
func uIn = U; // Dirichlet BC

// Stokes (linear) problem, used to initialized Newton’s method:  
problem Stokes( [ux,uy,uz,p] , [vx,vy,vz,q], solver=sparsesolver )  
= int3d(Th)( 1/Re\*(grad(ux)‘\*grad(vx) + grad(uy)’\*grad(vy) + grad(uz)'_grad(vz) )  
- div(vx,vy,vz) \* p  
+ div(ux, uy, uz) \* q ) /_ l.h.s. della weak formulation _/  
+ on(2,4,5,6, ux = 0 , uy = 0 , uz = 0 ) /_ Dirichlet BCs on the walls _/  
+ on(1, ux = 0 , uy = uIn , uz = 0 ); /_ Dirichlet BCs at the inflow _/  
/_ Neumann BCs at the outflow (automatically set by FreeFem) \*/

// Linearized Navier Stokes:  
problem linNS( [ux1,uy1,uz1,p1] , [vx,vy,vz,q], solver=sparsesolver )

```
= int3d(Th)(ugradv(ux,uy,uz,ux1,uy1,uz1)' * [vx,vy,vz] +
            ugradv(ux1,uy1,uz1,uz,uy,uz)' * [vx,vy,vz] +
		    1/Re * ( grad(ux1)'*grad(vx) + grad(uy1)'*grad(vy) + grad(uz1)'*grad(vz) ) -
	        div(vx,vy,vz) * p1 +
	        div(ux1, uy1, uz1) * q)
- int3d(Th)(ugradv(ux,uy,uz,ux,uy,uz)' * [vx,vy,vz])
+ on(2,4,5,6, ux1 = 0 , uy1 = 0 , uz1 = 0 )
+ on(1, ux1 = 0 , uy1 = uIn , uz1 = 0 );

```

real toll = 1e-06; // Error variation tolerance (of fluid velocity u)  
int kmax = 30; // Maximum number of iteration for Newton to diverge  
// Macro containing Newton’s iterations:

macro NavierStokes(){  
Stokes;  
for(int k = 0; k \< kmax; k++){  
linNS;  
dux = ux1 - ux;  
duy = uy1 - uy;  
duz = uz1 - uz;  
real err = sqrt( int3d(Th) ( grad(dux)‘\*grad(dux) + grad(duy)’\*grad(duy) + grad(duz)‘\*grad(duz) ) ) /  
sqrt( int3d(Th) ( grad(ux)’\*grad(ux) + grad(uy)‘\*grad(uy) + grad(uz)’\*grad(uz) ) );  
ux = ux1;  
uy = uy1;  
uz = uz1;  
p = p1;  
if( err \< toll ){  
cout \<\< "Value of the error after " \<\< k+1 \<\< " iterations: " \<\< err \<\< endl ;  
break;  
}  
if( k == kmax - 1 )  
cout \<\< “WARNING: Newton method does not converge in " \<\< kmax \<\< " iterations!” \<\< endl;

```
}

```

} //

// Non-linear problem solution with Newton’s method:  
NavierStokes;

absu = sqrt(int3d(Th)([ux,uy,uz]'_[ux,uy,uz]))+ 1e-8; // norma L2 del campo di velocità su tutto il dominio, che è una misura dell’energia cinetica totale del flusso.  
tauK = hK[].max/(2_absu);  
cout \<\< "|u|= " \<\< absu \<\< endl;

ofstream fileOut(“Pressure\_velocity.txt”); // Create a text file for Temperature

// Loop through the mesh points and save the scalar values  
for (int i = 0; i \< Th.nv; i++) {  
// Access the x, y, z coordinates and corresponding scalar value  
real x = Th(i).x;  
real y = Th(i).y;  
real z = Th(i).z;  
real valueP = p[i]; // Access the value of the scalar field at each node  
real valueux = ux[i]; // Access the value of the scalar field at each node  
real valueuy = uy[i]; // Access the value of the scalar field at each node  
real valueuz = uz[i]; // Access the value of the scalar field at each node

```
// Write coordinates and value to the file
fileOut << x << " " << y << " " << z << " " << valueP << " " << valueux << " " <<valueuy<< " " <<valueuz << endl;

```

};

plot(p, cmm = “pressure Solution”, value = true, fill=true, nbiso = 30);

Thank you in advance!

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [October 8, 2024, 12:42pm UTC](https://community.freefem.org/t/problem-solving-navier-stokes-equations-with-newton-method/3538/2 "2024-10-08T12:42:14Z")

</div>

Hello,  
Here is your file after a few corrections:  
[navierstokesnewton.edp](https://community.freefem.org/uploads/short-url/vxj7bGPnRiFI7HipzaIb5yUNxvK.edp) (4.2 KB)  
I have done the following changes:

```auto
5c5
< real Re = rhoa*U*L/mu;
---
> real Re = U*L/mu;
10c10
< int meshSize = 20;
---
> int meshSize = 10;
45c45
< ugradv(ux1,uy1,uz1,uz,uy,uz)' * [vx,vy,vz] +
---
> ugradv(ux1,uy1,uz1,ux,uy,uz)' * [vx,vy,vz] +
61,63c61
< dux = ux1 - ux;
< duy = uy1 - uy;
< duz = uz1 - uz;
---
> [dux,duy,duz] = [ux1,uy1,uz1] - [ux,uy,uz];
66,68c64
< ux = ux1;
< uy = uy1;
< uz = uz1;
---
> [ux,uy,uz] = [ux1,uy1,uz1];
70d65
< if( err < toll ){
71a67
> if( err < toll ){
92,98c88,94
< real x = Th(i).x;
< real y = Th(i).y;
< real z = Th(i).z;
< real valueP = p[i]; // Access the value of the scalar field at each node
< real valueux = ux[i]; // Access the value of the scalar field at each node
< real valueuy = uy[i]; // Access the value of the scalar field at each node
< real valueuz = uz[i]; // Access the value of the scalar field at each node
---
> real posx = Th(i).x;
> real posy = Th(i).y;
> real posz = Th(i).z;
> real valueP = p(posx,posy,posz); // Access the value of the scalar field at each node
> real valueux = ux(posx,posy,posz); // Access the value of the scalar field at each node
> real valueuy = uy(posx,posy,posz); // Access the value of the scalar field at each node
> real valueuz = uz(posx,posy,posz); // Access the value of the scalar field at each node
100c96
< fileOut << x << " " << y << " " << z << " " << valueP << " " << valueux << " " <<valueuy<< " " <<valueuz << endl;
---
> fileOut << posx << " " << posy << " " << posz << " " << valueP << " " << valueux << " " <<valueuy<< " " <<valueuz << endl;

```