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 UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)], [u1,u2]'*[dx(v2),dy(v2)]] //
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])
= int2d(Th)(nu * (Grad(du1,du2)' * Grad(v1,v2))+ UgradV(du1,du2, u1, u2)' * [v1,v2]
+ 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;
load "iovtk"
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 ?

Thank you in advance,

Best regards,

Loïc,