# Problem in Navier-Stokes 3D formulation

**URL:** <https://community.freefem.org/t/problem-in-navier-stokes-3d-formulation/2630>\
**Category:** General Discussion\
**Created:** [August 1, 2023, 1:47pm UTC](https://community.freefem.org/t/problem-in-navier-stokes-3d-formulation/2630 "2023-08-01T13:47:32Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![MariannaDellOtto](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mariannadellotto/32/1839_2.png) [@MariannaDellOtto](https://community.freefem.org/u/MariannaDellOtto)\
**Post date:** [August 1, 2023, 1:47pm UTC](https://community.freefem.org/t/problem-in-navier-stokes-3d-formulation/2630/1 "2023-08-01T13:47:32Z")

</div>

Hello,  
I am trying to write Freefem++ code for the steady Navier Stokes equations in 3D using Newton’s method taking my cue from the 2D model in the FreeFem documentation. The code I have written is as follows:

load “msh3”  
load “medit”

//Data  
real mu=1.825\*10^(-5); // viscosity  
real g=9.81; //gravitational acceleration  
real uMax = 10. ;// Inlet boundary velocity (for Dirichlet condition)

//Mesh Generation  
real L = 1;  
mesh3 Th = cube(20, 20, 20, [L_x,L_y,L\*z]);  
medit(“Mesh”, Th);  
// {2,4,5,6} = wall; 1 = inflow boundary ; 3 = outflow boundary

//Finite Elements Spaces:  
fespace Vh(Th, [P2, P2, P2]);  
Vh [ux,uy,uz] , [vx,vy,vz] ; // Fluid velocity, test function  
Vh [ux1,uy1,uz1] , [dux,duy,duz] ;

fespace Qh(Th, P1);  
Qh p , q ; // pressure e test function  
Qh p1 ;  
Qh rhoa = 1.293; /\* Just to define the macro, in this case I am assuming rhoa constant but in further implementations I shall need to define more complex functions for rhoa, so I want to keep the code as general as possible \*/

//Macro  
macro grad(v) [dx(v), dy(v), dz(v)] //  
macro div(vx,vy,vz) ( dx(vx) + dy(vy) + dz(vz) ) //  
macro DivRho(vx,vy,vz,rho) ( grad(rho)’ \* [vx,vy,vz] + rho \* div(vx,vy,vz) ) //  
macro ugradv(ux,uy,uz,vx,vy,vz) [[ux,uy,uz]‘\*grad(vx) , [ux,uy,uz]’\*grad(vy) , [ux,uy,uz]'\*grad(vz) ] //

// Macro for solving the steady nonlinear problem using Newton method:  
real toll = 1e-09; // Tolerance (for fluid velocity u)  
int kmax = 15; // max number of iterations

macro NavierStokes(){  
Stokes; // I start solving Stokes

```
for(int k = 0 ; k < kmax ; k++ ){

    linNS; // I compute u_(k+1) using linNS, with u_k taken from the previous iteration

    // I compute the error:
    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) ) );

    // I update the fluid velocity and the pressure:
    ux[] = ux1[];
    uy[] = uy1[];
    uz[] = uz1[];
    p[] = p1[];

    if( err < toll ){
	    // I print the error at the last iteration:
	    cout << "Error after " << k+1 << " iterations: " << err << endl ;
	    break;
    }

    if( k == kmax - 1 )
	    cout << "WARNING: Newton method does not converge in " << kmax << " iterations!" << endl;
}

```

} //

func h = 0; // Neumann boundary datum  
func uIn = uMax; // Dirichlet bounday datum (both are constant functions)

// Stokes problem (cut out the nonlinear term):  
problem Stokes( [ux,uy,uz,p] , [vx,vy,vz,q] )

```
= int3d(Th)( mu * ( grad(ux)'*grad(vx) + grad(uy)'*grad(vy) + grad(uz)'*grad(vz) ) +
	     mu * div(ux,uy,uz) * div(vx,vy,vz) -
         div(vx,vy,vz) * p +
	     DivRho(ux, uy, uz, rhoa) * q ) // l.h.s. della weak formulation

- int3d(Th)( rhoa * (-g) * vz ) // r.h.s. della weak formulation

+ on(2,4,5,6, ux = 0 , uy = 0 , uz = 0 ) // Dirichlet BCs at the wall
+ on(1, ux = 0 , uy = uIn , uz = 0 ) //Dirichlet BCs at the inflow
+ int2d(Th,3)( h*vy ); // Neumann BCs at the outflow

```

// Navier Stokes linarized with semi-implicit expression of the nonlinear term:  
problem linNS( [ux1,uy1,uz1,p1] , [vx,vy,vz,q] )

```
= int3d(Th)( 

	    rhoa * ugradv(ux,uy,uz,ux1,uy1,uz1)' * [vx,vy,vz] +

	     mu * ( grad(ux1)'*grad(vx) + grad(uy1)'*grad(vy) + grad(uz1)'*grad(vz) ) +
 	     mu * div(ux1,uy1,uz1) * div(vx,vy,vz) -
	     div(vx,vy,vz) * p1 +
	     DivRho(ux1, uy1, uz1, rhoa) * q )

- int3d(Th)( rhoa * (-g) * vz )

+ on(2,4,5,6, ux1 = 0 , uy1 = 0 , uz1 = 0 )
+ on(1, ux1 = 0 , uy1 = uIn , uz1 = 0 )
+ int2d(Th,3)( h*vy );

```

// Finally I solve the problem:  
NavierStokes;

// Plot:  
plot( p , ps=“Power.ps”);  
plot( [ux,uy,uz] , ps=“FluidVelocity.ps”);

Nevertheless, the program does not compile and I do not understand what mistake I am making since I followed step by step the model found in the documentation (obviously adapting it to my case).  
Thank you in advance for your help.

Best regards,  
Marianna

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 1, 2023, 3:08pm UTC](https://community.freefem.org/t/problem-in-navier-stokes-3d-formulation/2630/2 "2023-08-01T15:08:27Z")

</div>

Can you upload the code as one file? Copy from browser is a hassle.  
And while on the topic, curious if there is anyway to set up a github  
repository for code like this? That is, original posted code and a  
revision history that may be useful to archive.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 1, 2023, 5:32pm UTC](https://community.freefem.org/t/problem-in-navier-stokes-3d-formulation/2630/4 "2023-08-01T17:32:00Z")

</div>

You can’t put a comment in a macro, that is what terminates the macro.  
So you get the error on line 44 as it tries to compile that instead of just make  
iit part of the macro…

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 1, 2023, 10:32pm UTC](https://community.freefem.org/t/problem-in-navier-stokes-3d-formulation/2630/5 "2023-08-01T22:32:46Z")

</div>

If you want to comment in a macro use `/*this syntax*/`