# Problem in Navier-Stokes 3D formulation

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:

//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, [Lx,Ly,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 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[];

// 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 * 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).

Best regards,
Marianna

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.

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…

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