# Inquiries about solving the N-S equation, projection algorithm and parallel

**URL:** <https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006>\
**Category:** General Discussion\
**Created:** [February 28, 2024, 6:08pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006 "2024-02-28T18:08:38Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [February 28, 2024, 6:08pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/1 "2024-02-28T18:08:38Z")

</div>

Hello everyone，

I am not very familiar with FEM and have limited knowledge about its implementation for the Navier-Stokes equation. Currently, I am trying to understand the principles of FEM and parallel computing, but my learning progress is slow.

For now, my research involves solving 2D/3D N-S equation and temperature equation under a large range of Ra/Re number(10^5\<=Ra\<=10^10). Consequently, the computational cost increases significantly due to the Kolmogorov length scale. To study a 2D convective heat transfer problem, I tried to adapt the script [Rayleigh-Benard-time.edp](https://community.freefem.org/t/non-linear-time-dependent-problem/856/2) by [frederichecht](https://community.freefem.org/u/frederichecht), it works well but, slower and slower when Ra increases.  
My serial script uploaded:[R-B\_obstacles.edp](https://community.freefem.org/uploads/short-url/queNHIdzTKrX7S08Jg7ga0BPJUE.edp) (7.8 KB)

My supervisor recommends the projection algorithm and parallel computing.  
The projection method avoids solving a large system at a time, replaced by 2 or 3 decoupled equations for velocity and a Poisson-like equation for pressure. A implementation [NSprojection.edp](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/examples/NSprojection.edp) uses only **P1** elements, but some common practices using the pressure space 1 order lower than the velocity space such as **P2-P1**.  
I watched the tutorial [Pierre Jolivet introduction to parallel FreeFEM part 1](https://www.youtube.com/watch?v=-Aw2O46V2bo&list=PLYaPP8GBlaGf7AtAW8_NpFzJxEYvNLs5N), I noticed that the equation for pressure `pb4p(q, w, solver=CG, init=n, eps=epsp) = int2d(Th)(dx(q)*dx(w) + dy(q)*dy(w)) - int2d(Th)((dx(u) + dy(v))*w/dt)` could be solved in parallel just as the video did.

To find some power law through a series of simulations under a large range of Ra, I need to ensure my numerical results reflect the properties of the equation derived from true physical world correctly.

So the problems are:

1. Can NSprojection.edp handle high Reynolds numbers and be easily adapted in parallel?
2. If so, how to deal with convective term:

```auto
    = int2d(Th)(
      u*w/dt + nu*(dx(u)*dx(w) + dy(u)*dy(w))
    )
    - int2d(Th)(
      (convect([uold, vold], -dt, uold)/dt - dx(p))*w
    )

```

1. Is it reasonable to use the same element for u and p in projection algorithm? Is it better to use P2-P1 element?

Best regards

---

<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:** [February 29, 2024, 12:42pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/2 "2024-02-29T12:42:22Z")

</div>

Hello,  
I am not an expert on parallel issues, but I can give you some hints on the other parts.  
2. You should not use convect, which is relevant when there is no viscosity. Since you have viscosity you can solve the convection implicitly as

```auto
= int2d(Th)(u*w/dt + nu*(dx(u)*dx(w) + dy(u)*dy(w)))
  - int2d(Th)(uold*w/dt)
  + int2d(Th)(dx(p)*w)
  + int2d(Th)((uold*dx(u)+vold*dy(u))*w)

```

1. In the projection method it is not mandatory to satisfy an inf-sup condition (thus you can use P1/P1), but you will get a more accurate solution if the inf-sup condition is satisfied (as for P2/P1).

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [February 29, 2024, 5:58pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/3 "2024-02-29T17:58:27Z")

</div>

Thanks a lot! Your answer has addressed many of my concerns, I will try to implement a simple script using such an implicit scheme for convection.

Best regards,  
Haocheng Weng.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [March 1, 2024, 2:10pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/4 "2024-03-01T14:10:37Z")

</div>

> [@fb77](#):
>
> 1. You should not use convect, which is relevant when there is no viscosity.

Sorry but I am a bit confused. Why the built-in function `convect` relevant to inviscid flow?  
I have no idea because the convection term `u'*grad(u)` is without nu.  
Also, those serial examples related to N-S equation like [NSUzawaCahouetChabart.edp](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/tutorial/NSUzawaCahouetChabart.edp) using “convect” and they are indeed viscous flow. Could you please explain this further or if I have some misunderstanding?  
Thank you very much for your help,  
H. Weng

---

<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:** [March 1, 2024, 3:14pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/5 "2024-03-01T15:14:51Z")

</div>

The fact is that when there is no viscosity the problem is hyperbolic. Thus you should not use the implicit formulation because it could be unstable, and give weak or strong oscillations. You should use `convect`, which is adapted to the hyperbolic structure.  
With viscosity you have elliptic regularization of the gradient of the solution, that stabilizes the implicit resolution of convection (intuitively second-order derivatives dominate the first-order ones involved in the convection). This holds with a moderate Reynolds number `Re=L*u/nu`. But you are right in this case you can also use `convect` even if it is of really different nature than the spirit of finite elements (but it is truly relevant if you consider small viscosity i.e. large Reynolds number).  
Indeed since there is space discretization the real key value is the discrete Reynolds number `Re_h=h*u/nu` with `h` the size of the mesh cell. If `Re_h` is not too large then the implicit formulation is possible.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [March 1, 2024, 5:32pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/6 "2024-03-01T17:32:47Z")

</div>

Thank you very much for such detailed response. I need some time to digest those different stability properties of PDEs.

So, if the discrete Reynolds number Re\_h is not too large, computational stability can be guaranteed, right? From the definition `Re_h=h*u/nu`, it appears to imply that the size h is proportional to `nu/u`, when I take `Re_h` fixed. Before, when simulating some 3D pipe flow I just build a initial 2D mesh, and for different `Re=1/nu`, and do splitting to the initial mesh to ensure `h_local*Re^3/4=const`. My problem definition, similar to NSprojection.edp, using `convect`, but with an extra dimension z, I noticed that sometimes, during the transition to turbulence, the computation blows up when timestep too large.

Is it better, for the stability and accuracy of computing, to use something like `adaptmesh` to let local mesh size `h = const*nu/u`? Here, the exponent of `Re_h` is 1, whereas that of the Kolmogorov scale is `3/4`. I would like to know if maintaining the condition of an appropriate Re\_h is more strict than that of an appropriate Kolmogorov scale to resolve the details. Further more, what are the criteria for selecting an appropriate mesh size in DNS within FEM?

H. Weng

---

<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:** [March 1, 2024, 7:23pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/7 "2024-03-01T19:23:26Z")

</div>

The condition say `Re_h<1` does not mean that you get an accurate solution, it just means that your computation is stable with respect to space discretisation if you use the implicit method for advection.  
If you use `convect` this condition is not required.  
If you have other criterion on `h` to achieve accuracy, you have to respect them additionally.  
Concerning time stability, probably you need to respect that `nu*dt/h^2` and `u*dt/h` of order of unity at most.  
The mesh size h should be dictacted a priori by your data: it should be the cutoff length under which you do not wish to compute the details of the solution. But as you say it should be small enough in relation with nu in order to resolve your Kolmogorov scale.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [March 2, 2024, 2:11am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/8 "2024-03-02T02:11:31Z")

</div>

Thank you once again! Your comprehensive answer has provided me with a further understanding of the issue. Now it’s much more clear for me.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [April 29, 2024, 6:00am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/9 "2024-04-29T06:00:13Z")

</div>

I have written a FreeFEM code to compute vorticity formulation of incompressible flow using `problem`, which is a bit slow and can’t be parallelized.  
So I tried to write a matrix form but met some problem.  
I found an algorithm to solve this problem(JIAN-GUO LIU AND WEINAN E, 2000"SIMPLE FINITE ELEMENT METHOD IN VORTICITY FORMULATION FOR INCOMPRESSIBLE FLOWS"), which decouples the vorticity \omega and stream function \psi :

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/54e1b4489d994dc7171f25a9b666c6dc70c60ac5.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c99aea639b829af70869078cc6ff85e21389aada.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/61892598ad560b6ab0712874d5b9f108542675d2.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/2/25669e6ebb8b626db0256b21aa5a696332dd1d5f.png)  
The X\_{0,h}^k is the subspace of X\_h^k with zero boundary value.  
To compute an auxiliary term \bar{\omega}^{n+1}, I need to get a submatrix (A\_{0,0}A\_{0,b}) from the stiffness matrix A to compute right hand side of (I’), and compute the nonlinear term N.  
In my code using `problem`, the explicit nonlinear term (\nabla \varphi,\omega \nabla^{\perp}\psi) is:

```auto
grad(v1)'*[-dy(psiold),dx(psiold)]*omegaold

```

How to write a `varf` of nonlinear term as (2.6)?  
I know that I can use `varf Stiffness(u,v)=int2d(Th)(grad(u)'*grad(v));` and set `tgv=-10` to build a matrix of Stiffness with zero rows of label 1 , but I have no idea to build A\_{0,b}.  
I found the method to get list of DoF boundary in [Frédéric Hecht](https://community.freefem.org/u/frederichecht)’s reply:

> [@How do I obtain the matrix ordering of the boundaries?](https://community.freefem.org/t/how-do-i-obtain-the-matrix-ordering-of-the-boundaries/2834/2):
>
> I think , you need the list of DoF boundary. border C1(t=0,2\*pi) {x=cos(t); y=sin(t);label=1;} mesh Th=buildmesh(C1(6)); fespace Vh(Th,P1); Vh phi, w; varf vBord(phi,w) = int1d(Th)(w); Vh on1=vBord(0,Vh); for(int i=0; i\< Vh.ndof;++i) if (on1[][i]\>0) cout \<\< i \<\< " " \<\< endl;

If I get the list of DoF on boundary, should I permutate the rows and columns of `Stiffness` to get A？Or is there a better way to implement such algorithm?  
Could you give me some advice?

---

<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:** [April 29, 2024, 9:02am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/10 "2024-04-29T09:02:46Z")

</div>

I think that in order to compute your nonlinear term in (I’ ) you do not really need the matrix, because u and w are known. You can do something like

```auto
fespace Uh(Th, P1);
fespace Vh(Th, [P1dc,P1dc]);
Uh w,psi;
Vh wu;
w=...;
psi=...;
wu=w*[-dy(psi),dx(psi)];
varf WU(unused,v1)=int2d(Th)(grad(v1)'*wu);

```

On the question of extracting a submatrix, you have to make the list of indices corresponding to each block, then you can extract the submatrices. An example is  
[matrix-extract-bdry.edp](https://community.freefem.org/uploads/short-url/sdxqn9leYvAMNaUz3yI3S3dXRJC.edp) (1.3 KB)

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [May 14, 2024, 9:02am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/11 "2024-05-14T09:02:20Z")

</div>

Dear Prof. Bouchut,  
The vorticity formulation gives somewhat abnormal results, which conflicted with previous ones. I don’t have much time left to complete my paper, so I go back to solve the original N-S equation, but in parallel. I am working on a convection problem, unlike Boussinesq approximation, the density is modelled nonlinearly: \rho=\rho\_0(1-\beta|T-T\_M|^q),where β \> 0 is the thermal expansion coefficient and ρ\_0 is the density at T = T\_M.  
Conventionally, q=2, and the dimensionless equation set reads

\left\{ \begin{aligned} &\frac{\partial \textbf{u}}{\partial t}+\textbf{u}\cdot\nabla\textbf{u}=-\nabla p+\sqrt{\frac{Pr}{Ra}}\nabla^2\textbf{u}+(T-T\_M)^2\textbf{e}\_z \\ &\frac{\partial T}{\partial t}+\nabla\cdot (\textbf{u} T)=\frac{1}{\sqrt{PrRa}}\nabla^2T \end{aligned} \right.

I modified [cmd](https://community.freefem.org/u/cmd)’s parallel code. I added the temperature equation and a linear buoyancy term `dubT * vy` referring to T\textbf{e}\_z.

> [@Non homogeneous boundary and initial conditions](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/21):
>
> Dear Chris: I added a componet T for temperature and modified vR and vJ, but I am uncertain about the derivation of residual and Jacobian. the equations are: [image] Here I put the code: macro def(u)[u, u#y, u#T, u#p]// macro grad(u) [dx(u),dy(u)]//EOM macro dot(v, u) ( v \* u + v#y \* u#y + v#T \* u#T ) //EOM macro div(u) ( dx(u) + dy(u#y) ) // velocity divergence macro ugradu(v, U, u) ( v \* (U \* dx(u) + U#y \* dy(u)) + v#y \* (U \* dx(u#y) + U#y \* dy(u#y)) ) // scaled convection term macro visc…

When I turned to the nonlinear buoyancy term: (T-T\_M)^2\textbf{e}\_z, it reported an error with

```auto
varf vJ([u1,u2,uT,up],[v1,v2,vT,vp])=int2d(Th)(...-(uT-TM)^2*v2)

```

,other terms are omitted. I found in **varf** , a term should be either bilinear or linear, so the `uT^2*v2` referred to trilinear term, is illegal.  
I tried to substitute `uT^2*v2` to `uT*uTb*v2`, uTb is temperature at previous step, and I can run code. If it is a right doing? I also tried to use fully explicit term `(uTb-TM)^2`, but instability won’t be observed: the system remains stable at a high Rayleigh number, and the growing of time step does not have an effect.

---

<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:** [May 14, 2024, 10:24am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/12 "2024-05-14T10:24:40Z")

</div>

The Jacobian operator is the bilinear form you need and is the linearization of the nonlinear residual operator.

---

<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:** [May 14, 2024, 1:26pm UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/13 "2024-05-14T13:26:53Z")

</div>

Your two tries `uT*uTb*v2` and `(uTb-TM)^2` could be possible, numerical stability analysis would be necessary in order to know what is the best. In practice you can just try, see, and choose what you prefer. Next if the system does not behave numerically as you expect (unstable whereas you expect stable, or stable whereas you expect unstable), this is a hard problem of numerical analysis!

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [May 30, 2024, 11:01am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/14 "2024-05-30T11:01:26Z")

</div>

If you are NOT working on developing a FEM code in FreeFEM from scratch anf your aim is solve problems using FEM over the mid-range problem, I suggest checking out NEK5000. You can define your problem and run DNS as well as LES for wide-range of Re using high-accuracy and faster computations. It is a well-established software.

NOTE: Just a suggestion, if your aim is to study a problem and not develop code from scratch.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [May 30, 2024, 11:16am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/15 "2024-05-30T11:16:43Z")

</div>

Thanks for the suggestion. I know a bit about Nektar++, but not a skilled user. For the problem I am working on, it is hard for me to handle the complex mesh, and the multi-physical domain in Nektar++. Nek5000 or Nektar should have the competence to solve my problem, but at this stage, I don’t have much time to adapt my FreeFEM++ code to Nektar++.

Also, it takes time to be accustomed to the way Nektar++ solving a problem: to configure a .xml file and run Nektar++ session.

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [May 30, 2024, 11:38am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/16 "2024-05-30T11:38:16Z")

</div>

I am currently switching from a IBM based DNS code to NEK5000. It worked well for me. In NEK5000 or Nektar, you don’t have to code the NS equation as it is already there. You just have to tune some parameters. So once you get a hang of it, it is way faster than writing the code from scratch and optimizing it.

---

<div class="post-metadata">

**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\
**Post date:** [May 30, 2024, 11:49am UTC](https://community.freefem.org/t/inquiries-about-solving-the-n-s-equation-projection-algorithm-and-parallel/3006/17 "2024-05-30T11:49:00Z")

</div>

Nektar is a very good software, no doubt. For a “pure” fluid DNS, I will choose it.