# Navier Stokes with Newton method

Dear all,

I am trying to solve the Navier Stokes problem (a channel flow with a parabolic inflow) using the Newton method. I have the following code (inspired from what I have found on the FreeFEM tutorial Newton Method for the Steady Navier-Stokes equations).

``````verbosity = 0;

{
string Createrep="mkdir -p ./SOL_VTU ";
exec(Createrep);
}

mesh Th = square(30,30);

func int labelbord(){
return (1*(x==0)+2*(x==1)+3*(y==0)+4*(y==1)+10);} //function to put a label on each boundary edge of the mesh;

Th = change(Th, flabel=labelbord());

fespace Xh(Th, P2);
Xh u1, u2;
Xh v1,v2;
Xh du1,du2;

fespace Mh(Th,P1);
Mh p;
Mh q;
Mh dp;

// Macro
macro Grad(u1,u2) [dx(u1), dy(u1), dx(u2),dy(u2)] //
macro div(u1,u2) (dx(u1) + dy(u2)) //

// Initialization
u1 = 0;
u2 = 0;
real nu=1;
real err=1;
real eps = 1e-6;

for (int n = 0; n < 100; n++){

solve Oseen ([du1, du2, dp], [v1, v2, q])
+ UgradV(u1, u2,du1,du2)' * [v1,v2] - div(du1,du2) * q - div(v1,v2) * dp
- 1e-8*dp*q)
- int2d(Th) ( nu * (Grad(u1,u2)' * Grad(v1,v2)) + UgradV(u1,u2, u1, u2)' * [v1,v2] - div(u1,u2) * q - div(v1,v2) * p)
+on(11, du1=y*(1-y), du2=0) + on(12,du1=0, du2=0) + on(13,du1=0, du2=0) + on(14, du1=0 , du2=0) ; //+ on(10, du1=0, du2=0);

u1[]-=du1[];
u2[]-=du2[];
p[]-=dp[];

if (n>1){
real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
}

if(err < eps) break; //converge

if( n>3 && err > 10.) break; //blowup

cout << err << endl;

int[int] fforder=[1, 1];
savevtk("./SOL_VTU/SOL_REF.vtu", Th, [u1, u2], p, dataname = "velocity pressure", order=fforder);
}

``````

However, the results I obtain are false (the velocity is too high and I do not obtain a poiseuille as expected).

Could someone help me to understand what is wrong with this code ?

Best regards,

Loïc,

This is obviously wrong.

``````    u1[]-=du1[];
u2[]-=du2[];
p[]-=dp[];
``````

Why ?
and how to update the solution otherwise ?

Loïc,

I think this is the boundary condition that is wrong, it should be
`du1=u1-y*(1-y)`

Because `u1[]`, `u2[]`, and `p[]` are all the same. So you are basically updating the solution with the Newton step three times. My bad, you are not using a vectorial `fespace`, so I was wrong.
There is a steady-state NS solver available out-of-the-box at FreeFem-sources/examples/hpddm/navier-stokes-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub.

But yes, see the above answer: at convergence, your linear form will never evaluate to 0 because of the `on`, so it is not consistent.

Thank you ! Change the boundary conditions solves the problem.

Best regards,

Loïc,