# Pressure Poisson Equation (PPE)

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 19, 2011, 07:09 Pressure Poisson Equation (PPE) #1 New Member   Tom Join Date: Mar 2011 Posts: 12 Rep Power: 6 Hello all, I have been searching the forum and obviously there is a lot of information regarding the pressure poisson equation. However, I have not come across anyone who mentions the problem I have so here goes. Discretizing the PPE leads to a system of linear algebraic equations to be solved at each time-step and as the boundary conditions for the PPE are Neumann type the coefficient matrix is singular. I am using the Biconjugate Gradient (BCG) method, as detailed in NUMERICAL RECIPES in FORTRAN, as my solver. In order to get a solution a reference-pressure-point has to be set with the value of pressure at that node fixed. Regardless of the location of this reference point the solution from the BCG method produces a discontinuity in the pressure field surrounding the location of the reference point. This means that there are false pressure gradients which in turn produce false fluid flow. Has anyone else come across this problem when solving system of linear equations and if so do you know how to resolve the issue? Thanks for reading, Tom.

 July 19, 2011, 07:53 #2 Member   HouKen Join Date: Jul 2011 Posts: 66 Rep Power: 6 I think maybe you can first get rid of BCG and use a simple way like Gauss-Seidel to test.If even GS didn't work,that means maybe there is sth. wrong in your coefficient in PPE. If GS is ok then maybe you should turn to check the BCG. hope it can help u.

July 19, 2011, 08:20
#3
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by TheBoyce Regardless of the location of this reference point the solution from the BCG method produces a discontinuity in the pressure field surrounding the location of the reference point. This means that there are false pressure gradients which in turn produce false fluid flow.
This is not supposed to happen. Pressure field should be continuous.

I think you already asked this before too. I can tell you that there is something wrong with the way you are doing and it is not much to do with Poisson equation if Poisson equation was solved properly.

Are you sure that your equations are numbered properly??

One easy way to check this is to create a very small test case. Then export the matrix in some text file. Solve it to machine precision by your method and solve it by say matlab and then compare the solution. If everything is right then solutions will match, if they do not match then you are not solving it correctly.

 July 31, 2011, 08:24 #4 New Member   Amir Join Date: Oct 2010 Posts: 16 Rep Power: 6 Dear Tom, Your solver is probably true. Fixing the pressure of a point induce some errors that results in the issue you mentioned. There is some information in following link Singular Pressure Poisson Equation Ghasemi

July 31, 2011, 10:46
#5
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by Amir_Ghasemi Dear Tom, Your solver is probably true. Fixing the pressure of a point induce some errors that results in the issue you mentioned. There is some information in following link Singular Pressure Poisson Equation Ghasemi

NO it does not. When there is outflow boundary condition, pressure equation is all Neumann. This all Neuman equation is solved in Fluent,StarCCM, CFX etc etc. If it was true these all softwares would have problems just as he has.

 July 31, 2011, 23:04 #6 New Member   Amir Join Date: Oct 2010 Posts: 16 Rep Power: 6 Dear Arjun, As we have Neumann kind on all boundaries, the solution can be determined up to a constant (matrix is singular with a constant vector in null space). Fixing the pressure of a one point in the domain is not a correct approach to solve such a matrix and induce some errors around the fixed point. This is also mentioned in the link appeared in my last post. Ghasemi.

August 1, 2011, 00:45
#7
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by Amir_Ghasemi As we have Neumann kind on all boundaries, the solution can be determined up to a constant (matrix is singular with a constant vector in null space).
Well to my knowledge this is how we could solve all Neumann kind of problem since there is no unique solution.

Quote:
 Originally Posted by Amir_Ghasemi Fixing the pressure of a one point in the domain is not a correct approach to solve such a matrix and induce some errors around the fixed point. This is also mentioned in the link appeared in my last post.
Why it is not a correct way of solving it.

If upon solution , all the solutions only differ by a constant then no matter what way you solve it, if it is a solution it would differ by constant to any solution you would calculate by the method in the link you provided.

So this means if I get a solution by fixing pressure at any point in domain and if solution converged than it is guaranteed to differ by solution you would get by other method by a constant.

Further, there will be no errors anywhere in domain if the solution converged to machine precision because this is what we mean by convergence.

 August 1, 2011, 01:09 #8 New Member   Amir Join Date: Oct 2010 Posts: 16 Rep Power: 6 Dear arjun, I do not know very well about technical details so I can't convince you. But I am sure that such error may occur around fixed points as I have experienced it by myself and I have heard similar issue from my friend and also it is mentioned in other threads in CFD-online. But probably such problem may not occur for all matrices and this is why me and you have different experience on this problem. Ghasemi.

August 1, 2011, 01:33
#9
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by Amir_Ghasemi Dear arjun, I do not know very well about technical details so I can't convince you. But I am sure that such error may occur around fixed points as I have experienced it by myself and I have heard similar issue from my friend and also it is mentioned in other threads in CFD-online. But probably such problem may not occur for all matrices and this is why me and you have different experience on this problem. Ghasemi.

I have worked with this issue a lot.
First time I wrote additive corrective multigrid, I was not aware that all neumann problem could be pain in A. And faced lack of convergence.

And this is the reason that null space based method is called correct way of solving it. Because that way it converges and may be converges faster than normal methods. Like for example in all neumann problem additive corrective multigrid converges very very poorly and many times if problem is large it does not even converge.
For this reason I switched to smoothed aggregation multigrid and now it converges very fast even for all neuman problem where I fix value at one point in domain.

I still use BiCGstab method preconditioned with smoothed aggregation multigrid and it works. I have used this with millions of points in domain. Works like charm.

About the issue of velocity shoot up or down. (that you think is due to fixing a point pressure).

In navier stokes equations or in CFD that could happen even in well poissed poisson problem like problems with Dirichlet BCs of pressure.

But this has not much to do with lack of convergence but rather NON SMOOTHNESS of solution of pressure equation.
Why? Because velocity is corrected by pressure correction gradients. If solution is non smooth gradients vary very much and that cause velocity shoot up in flow domain.

To avoid this and make solver more stable, I apply gradient limiter on pressure correction in iNavier. This keeps solver much more stable than it normally is.

 August 1, 2011, 01:43 #10 Senior Member   Arjun Join Date: Mar 2009 Location: Nurenberg, Germany Posts: 369 Rep Power: 10 By the way, if someone is looking for reading something about how to solve singular system, he should look for book http://www.amazon.com/Multilevel-Blo...2177265&sr=8-1 In this book, please refer chapter 8 Preconditioning Nonsymmetric and Indefinite Matrices. This chapter lists most of the popular methods for solving such system , including the null space orthogonalization based methods. they are given in good detail.

 August 1, 2011, 11:33 #11 New Member   Amir Join Date: Oct 2010 Posts: 16 Rep Power: 6 Dear arjun, Thank you very much for the reference, I was following for such book. Ghasemi.

 August 4, 2011, 01:31 #12 New Member   Amir Join Date: Oct 2010 Posts: 16 Rep Power: 6 Hello all, I have just found what the problem was. My unsymmetric matrix was singular with constant vector in null space and its determinate was zero and it had emerged from discretizing PPE with Neumann boundaries. One should be aware that these are not enough to decide to remove one equation (or fixing the pressure of one point). you can do this, if the system of equation has infinite solutions. Ax = b has infinite solution if det(A) = 0 and b is in the range of A, if these condition exist you can remove one equation and be sure that the result of new system also satisfy the original one and in may system this was not the case. Ghasemi

August 7, 2011, 12:34
#13
Senior Member

Join Date: Aug 2011
Posts: 251
Rep Power: 6
Quote:
 Originally Posted by Amir_Ghasemi Dear Arjun, As we have Neumann kind on all boundaries, the solution can be determined up to a constant (matrix is singular with a constant vector in null space). Fixing the pressure of a one point in the domain is not a correct approach to solve such a matrix and induce some errors around the fixed point. This is also mentioned in the link appeared in my last post. Ghasemi.
I'm fully agree with Amir. Fixing one pressure point in case of singular linear system equations provided by a Poisson equation with Neumann boundary conditions on all boundaries does not work despite what many people say!
Or if it works it's not exactly what they do...
But anyway thanks Arjun for the reference, I will have a look...

If one wants to benchmark its own solution you can use fishpack:
http://www.cisl.ucar.edu/css/software/fishpack/
which is a package software to solve poisson equation with many different types of boundary condtions including full Neumann boundary conditions. I don't know how they treat the problem especially because they also use direct solver but it can provide good solution.
You can then compare your solution while fixing one pressure point and the solution provided by fishpack.

 August 7, 2011, 12:48 #14 New Member   Amir Join Date: Oct 2010 Posts: 16 Rep Power: 6 Dear leflix, if your matrix is symmetric, check that if the summation of your source terms is zero, if this is the case then your matrix has infinite solution and fixing the pressure of one point would not induce a problem. consider following discussion, Ax = b C is a constant vector we know that A*C = 0, in case of all Neumann boundary if A is symmetry then C'*A is also zero, then you can conclude that C'*A = C'*b = 0, it means that summation over b entries must be zero to have a solution (up to a constant) for the system. above discussion is valid for symmetric A. in case of non-symmetric there should also be a such condition which i am not aware of that! Ghasemi

August 7, 2011, 18:59
#15
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by leflix I'm fully agree with Amir. Fixing one pressure point in case of singular linear system equations provided by a Poisson equation with Neumann boundary conditions on all boundaries does not work despite what many people say!

Do you have a case in mind for which you think it would not work. I would love to solve it with iNavier and show to you that it works.

As far as I know these are the softwares which fix one point value and all of them are working just fine: Fluent, CFX, StarCD, StarCCM+, Polyflow.

 August 8, 2011, 04:52 #16 Senior Member   Join Date: Aug 2011 Posts: 251 Rep Power: 6 Ok Arjun could you recall exactly the procedure you follow with your method based on fixing one reference point ? Let's imagine you want to solve LAP(P)=f with dp/dn =0 on all boundaries. LAP is the Laplacian orperator and f any given function. let's P_ref the value you want to fix at one given node and P* the resulting pressure field obtained from solving this new linear system: (LAP(P)=f + dp/dn=0 + P=P_ref at one node). what do you do after? As you kow it does not exist analytical solution for Poisson equation with full Neumann BC. so what you could do is to solve LAP(P)=f for any given f function. Then from the field P obtained recompute Grad(P) and then DIV(Grad(P)) to obtain LAP(P). Then compare your solution with f. If you got it then it's correct...

August 8, 2011, 05:05
#17
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by leflix Ok Arjun could you recall exactly the procedure you follow with your method based on fixing one reference point ? Let's imagine you want to solve LAP(P)=f with dp/dn =0 on all boundaries. LAP is the Laplacian orperator and f any given function. let's P_ref the value you want to fix at one given node and P* the resulting pressure field obtained from solving this new linear system: (LAP(P)=f + dp/dn=0 + P=P_ref at one node). what do you do after? As you kow it does not exist analytical solution for Poisson equation with full Neumann BC. so what you could do is to solve LAP(P)=f for any given f function. Then from the field P obtained recompute Grad(P) and then DIV(Grad(P)) to obtain LAP(P). Then compare your solution with f. If you got it then it's correct...

I am in hurry at the moment so here is quick outline of the approach that i took.

If you descretize a Poisson problem with all neumann you get this:

Ap = |sum_neigh(Al)|
That is diagonal term equal to sum of off diagonal terms.

Numerically solution to exist, Ap > |sum_neigh(Al)| for AT LEAST ONE equation in system.

Now imagine that you pick a cell or control volume that has a boundary face. Remember we picked just 1 such cell. Lets fix the value of pressure or pressure correction at that face center , I set it to 0.

That means the term that would get into source will be +Al * 0

and the value that goes inside the diagonal term is +Al . So for this cell, now
Ap > |sum_neigh(Al)|

Hence the solution exist and since we apply Dirichlet condition on that face this is equivalent to solving system with 1 face with fixed value.

This is difficult to solve but it does work in Navier Stokes equation solver. One has to be careful about keeping that face pressure to 0 or whatever value you set during whole Navier Stokes solution process.

 August 8, 2011, 07:49 #19 Senior Member   Arjun Join Date: Mar 2009 Location: Nurenberg, Germany Posts: 369 Rep Power: 10 i am too busy at the moment so when i get time in a day or two, i will construct a small test say 2D mesh with all neumann boundary condition and will try to solve it by fixing a value at a point and see how it works out. I did transform the equations but this is exactly what will happen if you FIX a value at a point.

August 8, 2011, 10:05
#20
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
Quote:
 Originally Posted by leflix What you did is to transform a singular linear system in a non singular one which gives you the possibility to obtain a solution. It seems that your criteria to judge that "it works " or not is that you get a solution from your Navier-Stokes solver or not.
I did little home work. I should admit you were correct about it. Fixing a value at one point at BC, does produce solution to that matrix but in the end gradient across that point and cell center is not 0. Thus it is not true all neumann problem.

So you are correct about "your criteria to judge that "it works " or not is that you get a solution from your Navier-Stokes solver or not."

Because I use pressure correction method, this issue did not bother me much because eventually solution converges to a pressure field that satisfies Navier Stokes.

Admittedly I never much bothered about this in detail too because I know that I am supposed to fix a value somewhere to make it work.

Quote:
 Originally Posted by leflix But which coupling algorithm in your Navier-Stokes solver do you use ? Projection method with a real poisson equation to solve? or SIMPLE like algorithm with a Poisson like equation? The consequences won't be the same.
I use SIMPLE so at the moment it is not an issue.
But I am exploring possibilities of applying monolithic schemes where saddle point problem is solved for u and p. Plus exploring Braess Sarazin type methods too.
In this issue this issue might be big trouble, I would need to look at it again.

Quote:
 Originally Posted by leflix How behaves the pressure solution in the cell where you fixed the value on the boundary face compared to the adjacent cells ? I'm not convinced that your Navier-Stokes solver providing a solution is a sufficient criteria to tell that your poisson solver with this procedure is valid.
It seem so now to me too.

Quote:
 Originally Posted by leflix Another question, which iterative linear system solver do you use too ?

I use smoothed aggregation multigrid. Or BiCGStab preconditioned by smoothed aggregation multigrid.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post dontknow Main CFD Forum 2 August 30, 2010 10:12 Sas CFX 15 July 13, 2010 08:56 Miguel Baritto CFX 4 August 31, 2006 12:02 abilge Main CFD Forum 12 March 29, 2005 12:34 frederic felten Main CFD Forum 1 November 17, 2000 04:17

All times are GMT -4. The time now is 04:38.

 Contact Us - CFD Online - Top