# Pressure projection methods **URL:** \ **Category:** General Discussion\ **Created:** [May 14, 2024, 8:49am UTC](https://community.freefem.org/t/pressure-projection-methods/3247 "2024-05-14T08:49:31Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 14, 2024, 8:49am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/1 "2024-05-14T08:49:31Z")
Hello, everyone. I’m a beginner in FreeFEM. Currently, I’m attempting to decouple the velocity field and pressure field of the Navier-Stokes equations using the pressure projection method. In the first step, I aim to calculate the convergence order to validate the accuracy of my code. However, I’ve encountered some issues. Here is my code[test\_pressure\_projection.edp](https://community.freefem.org/uploads/short-url/nvH9665JSoIgBU87KmPFFgOfsHg.edp) (3.7 KB). I’m not sure whether it’s an issue with my understanding of the projection method or simply a coding problem. The format I’m using is: ![截屏2024-05-14 16.42.45](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/2/264f64aa2482fda4a073537642260344d3f630c3.png) ---
**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, 3:55pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/2 "2024-05-14T15:55:04Z")
Hello, Your algorithm is correct (but you need to take care of the boundary condition on p^{n+1}, this is a sensitive point). You forgot to code the convection. Moreover your computation of right-hand side f was mistaken. The way you evaluate the order of convergence is not good: you cannot evaluate the convergence in time by diminishing the timestep at fixed mesh. You need to reduce both simultaneously (I did not modify that). Francois. [test\_pressure\_projection.edp](https://community.freefem.org/uploads/short-url/6RoyH7RPJhxbEZh8plVJrjKpFX5.edp) (3.9 KB) ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 14, 2024, 4:33pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/3 "2024-05-14T16:33:06Z")
Thank you very much for responding to my questions. I looked up some references on the pressure projection method. In the second variational problem of updating pressure, the Neumann boundary condition of pressure is required. In practice, this condition is responsible for the errors this method shows close to the boundary of the domain since the real pressure (i.e., the pressure in the exact solution of the Navier-Stokes equations) does not satisfy such boundary conditions. But if the velocity can satisfy `u·n=0` boundary condition, the Neumann boundary conditions of the pressure will naturally be obtained. I see that what you added in the code is `on(1,2,3,4,pnp1=pe(t))`Dirichlet boundary condition. Is it correct? ```auto Convergence order in L2 norm of u: 0, 0.158643, 1.05926, 0.962136, 0.74141, 0.496186. Convergence order in H1 norm of u: 0, 0.153184, 0.967829, 0.871569, 0.691116, 0.541268. convergence order in L2 norm of p: 0, 1.72902, -1.21749, 0.519209, 0.767257, 0.801716. ``` Is this the result of me running the code, or is it incorrect and is it caused by the above issue? In addition, for the convection term, do I have to use `convect`? Can’t I define `macro UgradV(ux,uy,vx,vy) [ux*dx(vx) + uy*dy(vx), ux*dx(vy) + uy*dy(vy)] //` myself? ---
**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, 5:19pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/4 "2024-05-14T17:19:29Z")
I put Dirichelt BC on p with the exact solution on the boundary in order to be sure that you get the right solution. This is of course nor possible when you do not know the exact solution. Thus in general you need to do something else, and I am not at all expert on that (it seems that you know more than me on that). For the convection you can do it also with UgradV implicitly, but for that you need enough viscosity for it to be stable. ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 15, 2024, 4:25am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/5 "2024-05-15T04:25:31Z")
Thank you very much for your reply. I have two more questions. In 3D, the convection term is expressed as `convect([u1,u2,u3],−dt,cold)` ？ In addition to the convection, are there other expressions affected by physical parameters, such as the diffusion second-order curl term `curl(curl B)` and the convection first-order curl `curl(B \times u)$` in the magnetic fluid equation? Can I define curl and cross product operators directly using macro? ```auto macro curl(Bx,By,Bz) [dy(Bz)-dz(By), dz(Bx)-dx(Bz), dx(By)-dy(Bx)] // macro crossProduct(ux,uy,uz,vx,vy,vz) [uy*vz-uz*vy, uz*vx-ux*vz, ux*vy-uy*vx] // ``` ![截屏2024-05-15 12.16.57](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b9be367485d865d28f2642cd0fe513171a48cd57.png) ---
**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 15, 2024, 8:56am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/6 "2024-05-15T08:56:22Z")
I have not tried convection in 3d, but your formula looks correct. Your macros are correct also, except that you must put some () for cross product: `macro crossProduct(ux,uy,uz,vx,vy,vz) [(uy)*(vz)-(uz)*(vy), (uz)*(vx)-(ux)*(vz), (ux)*(vy)-(uy)*(vx)] //` This is because `macro` performs a formal replacement with the arguments (the arguments have to be understood as strings, not as numbers). In case the argument ux is “a+b” the macro without parenthese would produce terms `-a+b*vz` and `a+b*vy`, missing parentheses! ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 16, 2024, 2:33am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/7 "2024-05-16T02:33:14Z")
Thank you ,Professor. I will correct the code. ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 19, 2024, 8:36am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/8 "2024-05-19T08:36:41Z")
Hello dear professor. I looked up some more information about the pressure projection method. Now, regardless of the accuracy of this method, I want to add the following boundary conditions to the second step update pressure and the third step update speed, how do I enter the code? The boundary condition of the pressure in the second step is derived from `u·n = 0 ( or = ue·n )` ![截屏2024-05-19 16.32.56](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/5c6d97e94ae70759fb34133f9f6f20779c3200a0.png) Here’s my updated theoretical algorithm ![截屏2024-05-19 16.32.43](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/f/ff1142f4129698998b06104f229144709fca2050.png) ---
**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 19, 2024, 3:50pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/9 "2024-05-19T15:50:02Z")
The equation you have to solve for p^{n+1} is a Laplace equation with non homogeneous Neumann boundary condition. You can look at this post > [@Uzawa scheme for stokes bingham flow](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/16): > > In that paper the flow is taken incompressible, whereas in your code you use compressible. This is why the solution is different (but looks correct). You can modify easily your code to do the incompressible case, by using the extra variables p,q as solve bing([ux,uy,p],[vx,vy,q]) and adding -int2d(Th)(p\*Divergence(vx, vy)+q\*Divergence(ux, uy)) But then it happens that the scheme converges extremely slowly (you need thousands of iterations). In order to improve it you should use the Augmente… From “We are thus led to the problem…” ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 19, 2024, 4:17pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/10 "2024-05-19T16:17:35Z")
Yes, professor. But why do not you enter `f` and `g` in this code? What is `dyg2*(dy(ux) - dx(uy))`? I want to know the weak form by which the code is written. ![截屏2024-05-20 00.13.18](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/1/16b5bd7a167509fb6e3db4207a7e62fcc01fa870.png) ```auto fespace Rh(Th, P2); Rh dyg,dyg2; solve streamlines (dyg, dyg2,solver=UMFPACK) = int2d(Th)( dx(dyg)*dx(dyg2) + dy(dyg)*dy(dyg2) ) + int2d(Th)( - dyg2*(dy(ux) - dx(uy)) ) + int1d(Th)( - dyg2*(uy*N.x-ux*N.y) ) + int1d(Th)(1e-8*dyg*dyg2) ; ``` ---
**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 19, 2024, 5:14pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/11 "2024-05-19T17:14:28Z")
replace dy(ux) - dx(uy) by f replace uy_N.x-ux_N.y by g ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 20, 2024, 2:17am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/12 "2024-05-20T02:17:12Z")
Okay, Professor. I understand what you mean. ![截屏2024-05-20 09.52.47](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/00922e05a771fb1c7a998d2b50e82a8d5e6ea3d7.png) Pressure projection method first step is to calculate \boldsymbol{\tilde{u}}, BC is \boldsymbol{\tilde{u}}|\_{\partial\Omega} = 0, but for examples with analytical solutions \boldsymbol{u}\_e, I think BC is \boldsymbol{\tilde{u}}|\_{\partial\Omega} = \boldsymbol{u}\_e. The third step is to calculate \boldsymbol{u}^{n+1}, which cannot use Dirichlet BC, so take \boldsymbol{u}^{n+1}\cdot\boldsymbol{n}|\_{\partial\Omega} = 0(theoretical analysis), but For examples with analytical solutions, the boundary condition should be \boldsymbol{u}^{n+1}\cdot\boldsymbol{n}|\_{\partial\Omega} = \boldsymbol{u}\_e\cdot\boldsymbol{n}. Can I analyze it this way, professor? So, I can get \nabla (p^{n+1} - p^n) \cdot \boldsymbol{n} = - \frac{\boldsymbol{u}^{n+1} -\boldsymbol{\tilde{u}}}{\Delta t}\cdot\boldsymbol{n} = 0. What should I enter for pressure BC: ```auto solve pb4up(pnp1,pp) =...(omit this part) + int1d(Th)( (dx(pn)*N.x + dy(pn)*N.y)*pp ); ``` Is it right, professor? ---
**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 20, 2024, 10:05am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/13 "2024-05-20T10:05:12Z")
First remark, if the original problem you want to solve is `u=0` on boundary, they you have that `u_e=0` on boundary. Then to solve the homogeneous Neumann BC for p^{n+1} or p^{n+1}-p^n, you use the formulation I sent with` g=0`. Thus the only `int1d` that remains is the one with `1e-8`. You do not have to put the term dx(pn)\*N.x + dy(pn)\*N.y. ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 24, 2024, 8:18am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/14 "2024-05-24T08:18:36Z")
Okay, Professor. I understand what you mean. The third step: update \boldsymbol{u}^{n+1} with \frac{\boldsymbol{u}^{n+1} -\boldsymbol{\tilde{u}}}{\Delta t} + \nabla p^{n+1} - \nabla p^n = 0. Should I use the variational problem and add the boundary condition \boldsymbol{u}^{n+1}\cdot\boldsymbol{n} = 0 or \boldsymbol{u}^{n+1}\cdot\boldsymbol{n} = \boldsymbol{u}\_e\cdot\boldsymbol{n}, like this ```auto solve pb4u(uxnp1, uynp1, vx, vy) = int2d(Th)([uxnp1, uynp1]'*[vx, vy]) - int2d(Th)([uxtilde, uytilde]'*[vx, vy] - dt(i)*grad(pnp1)'*[vx, vy] + dt(i)*grad(p0)'*[vx, vy]) + int1d(Th,1,2,3,4)(1.e10*(uxnp1*N.x + uynp1*N.y)*(vx*N.x + vy*N.y) ) - int1d(Th,1,2,3,4)(1.e10*(uxe(t)*N.x + uye(t)*N.y)*(vx*N.x + vy*N.y) ); ``` This is my second step. ```auto solve pb4p(pnp1, pp) = int2d(Th)(grad(pnp1)'*grad(pp)) - int2d(Th)(grad(p0)'*grad(pp) - div(uxtilde, uytilde)*pp/dt(i)) + int1d(Th,1,2,3,4)(1e-8*pnp1*pp) ; ``` ---
**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 24, 2024, 8:48am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/15 "2024-05-24T08:48:09Z")
The third step is an L^2 projection thus there is no boundary condition. You should not put any int1d. Alternatively you can just set ```auto uxnp1 = uxtilde - dt(i)*dx(pnp1)+ dt(i)*dx(p0); uynp1 = uytilde - dt(i)*dy(pnp1)+ dt(i)*dy(p0); ``` This performs an interpolation instead of an L^2 projection. ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 25, 2024, 6:42am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/16 "2024-05-25T06:42:01Z")
Thank you very much, professor. I changed my code based on your guidance. [pressure\_projection\_BDF2.edp](https://community.freefem.org/uploads/short-url/ay5znQYSMBzrC5UPBXo3tNmV5Bv.edp) (6.0 KB) This is the standard incremental pressure-projection method. The time format is BDF2, and the first step (n=0) is backward Euler. Linearized convection term: n = 0, (\boldsymbol{u}^n\cdot\nabla)\boldsymbol{u}^n; n \>= 1, ((2\*\boldsymbol{u}^n - \boldsymbol{u}^{n-1})\cdot\nabla)(2\*\boldsymbol{u}^n - \boldsymbol{u}^{n-1}). This is my output: ```auto L2 norm error of u: 8.90828e-05, 5.76243e-05, 1.52496e-05, 3.91707e-06, 9.99719e-07, 2.78619e-07. H1 norm error of u: 0.00129628, 0.000582098, 0.000196062, 0.000120427, 0.00011179, 0.000111002. L2 norm error of p: 0.000535824, 0.000129077, 3.14533e-05, 8.38204e-06, 2.28884e-06, 6.33111e-07. Convergence order in L2 norm of u: 0, 0.62847, 1.9179, 1.96093, 1.97018, 1.84323. Convergence order in H1 norm of u: 0, 1.15505, 1.56995, 0.703146, 0.107374, 0.0102053. Convergence order in L2 norm of p: 0, 2.05353, 2.03695, 1.90784, 1.87268, 1.85409. ``` I really don’t understand why the H1 norm convergence order of velocity is so bad? I test time accuracy. I refer to page 4 of this paper (taking Stokes equation as an example): [An overview of projection methods for incompressible flows.pdf](https://community.freefem.org/uploads/short-url/hH1ESAZTNDGpYHk7LhM8zBuQ9Vp.pdf) (943.4 KB) ![截屏2024-05-25 14.32.08](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/3412551e223725b35646543d2ada439de84d483a.png) ![截屏2024-05-25 14.32.25](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/4/48ad22d8183739b9cc9aeaa17cf1771f27e22390.png) ---
**Author:** ![emmawilson](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/emmawilson/32/2346_2.png) [@emmawilson](https://community.freefem.org/u/emmawilson)\ **Post date:** [May 25, 2024, 7:45am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/17 "2024-05-25T07:45:05Z")
Could you please guide me why my post is hidden? ![Screenshot (197)](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/dd72a4dc39be89a2a337dace1649be2a876860f2.png) ---
**Author:** ![LXX](https://avatars.discourse-cdn.com/v4/letter/l/258eb7/32.png) [@LXX](https://community.freefem.org/u/LXX)\ **Post date:** [May 25, 2024, 8:58am UTC](https://community.freefem.org/t/pressure-projection-methods/3247/18 "2024-05-25T08:58:09Z")
So sorry, I do not know why. ---
**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 25, 2024, 5:41pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/19 "2024-05-25T17:41:45Z")
What this paper says is also what I am thinking of: the Neumann boundary condition that is imposed on p is not necessarily satisfied by the exact solution. As a consequence it could be that the pressure well converges to the exact solution far from the boundary, but does not converge close to the boundary (probably within a length of the order of the mesh size). This is what is called a boundary layer. The consequence is that the convergence rate on the pressure could be lower than expected. The H1 norm of u has usually a behavior similar to that of the pressure. I would expect a rate of 2 for L2u, and a rate of 1 (or lower) for H1u and L2p. But this behavior is for when you refine also the mesh, not only the timestep as you do. However this is not what you observe: you get a low rate for the H1 norm, and a surprisingly high rate for the pressure. This is a bit mysterious. This could be due to the fact that you refine only in time, or to the BDF2 formulas that you apply to u but not to p. I don’t know. I am not sure also of the explicit treatment of the convection. You could try it implicit as for example `UgradV(uxbar, uybar, uxtilde, uytilde)` ---
**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 27, 2024, 3:25pm UTC](https://community.freefem.org/t/pressure-projection-methods/3247/20 "2024-05-27T15:25:06Z")
Several comments: 1. After some tries, I find that taking the convection implicit is slightly better. 2. It is also slightly better in the step of computing utilde, to take the pressure guess as pbar= 2\*pn-pnm1 instead of pn. Then it is also pbar in the step of computing pnp1 by the Laplace equation and ```auto uxnp1 = uxtilde - dt(i)*dx(pnp1)+ dt(i)*dx(pbar); uynp1 = uytilde - dt(i)*dy(pnp1)+ dt(i)*dy(pbar); ``` 1. It seems that there is no boundary layer, since if we replace the Neumann boundary condition on p by Dirichlet p=p\_exact, nothing changes. 2. The apparent lack of convergence in the H1 norm is due to the fixed mesh 100\*100. Indeed an idea of the H1 error associated to the mesh can be obtained setting just after the definition of initial `ux0,uy0` ```auto real errdiv0=sqrt(int2d(Th)(abs(dx(ux0)+dy(uy0))^2)); cout << "errdiv0 " << errdiv0 << endl; ``` The result is 1.4e-4. This means that a typical H1 error induced by projection onto the finite element space is of this order. It follows that even for high order time integration and small timestep, the H1 norm cannot get below this value (without refining the mesh). This is the reason why once the H1 error attains this value, it will no longer decay. [Next page](https://community.freefem.org/t/pressure-projection-methods/3247.md?page=2)