# Navier Stokes with Newton method

**URL:** <https://community.freefem.org/t/navier-stokes-with-newton-method/3371>\
**Category:** General Discussion\
**Created:** [July 1, 2024, 2:59pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371 "2024-07-01T14:59:51Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [July 1, 2024, 2:59pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/1 "2024-07-01T14:59:51Z")

</div>

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](https://doc.freefem.org/tutorials/navierStokesNewton.html)).

```auto
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,

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 1, 2024, 3:05pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/2 "2024-07-01T15:05:19Z")

</div>

This is obviously wrong.

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

```

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [July 1, 2024, 3:08pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/3 "2024-07-01T15:08:34Z")

</div>

Why ?  
and how to update the solution otherwise ?

Loïc,

---

<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:** [July 1, 2024, 3:11pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/4 "2024-07-01T15:11:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 1, 2024, 3:11pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/5 "2024-07-01T15:11:47Z")

</div>

~~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](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/navier-stokes-2d-PETSc.edp).

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 1, 2024, 3:14pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/6 "2024-07-01T15:14:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Loic](https://avatars.discourse-cdn.com/v4/letter/l/e9c0ed/32.png) [@Loic](https://community.freefem.org/u/Loic)\
**Post date:** [July 1, 2024, 3:22pm UTC](https://community.freefem.org/t/navier-stokes-with-newton-method/3371/7 "2024-07-01T15:22:19Z")

</div>

Thank you ! Change the boundary conditions solves the problem.

I will also see [FreeFem-sources/examples/hpddm/navier-stokes-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/navier-stokes-2d-PETSc.edp)

Best regards,

Loïc,