
[Sponsors] 
July 19, 2011, 07:09 
Pressure Poisson Equation (PPE)

#1 
New Member
Tom
Join Date: Mar 2011
Posts: 12
Rep Power: 8 
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 timestep 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 referencepressurepoint 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: 8 
I think maybe you can first get rid of BCG and use a simple way like GaussSeidel 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: 742
Rep Power: 19 
Quote:
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: 9 
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 http://www.cfdonline.com/Forums/mai...equation.html Ghasemi 

July 31, 2011, 10:46 

#5  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 742
Rep Power: 19 
Quote:
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: 9 
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: 742
Rep Power: 19 
Quote:
Quote:
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. Think about it. 

August 1, 2011, 01:09 

#8 
New Member
Amir
Join Date: Oct 2010
Posts: 16
Rep Power: 9 
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 CFDonline. 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: 742
Rep Power: 19 
Quote:
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: 742
Rep Power: 19 
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/MultilevelBlo...2177265&sr=81 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: 9 
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: 9 
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: 271
Rep Power: 9 
Quote:
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: 9 
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 nonsymmetric 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: 742
Rep Power: 19 
Quote:
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: 271
Rep Power: 9 
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: 742
Rep Power: 19 
Quote:
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, 05:42 

#18 
Senior Member
Join Date: Aug 2011
Posts: 271
Rep Power: 9 
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 NavierStokes solver or not.
But which coupling algorithm in your NavierStokes 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. 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 NavierStokes solver providing a solution is a sufficient criteria to tell that your poisson solver with this procedure is valid. So when you will get time, to figure out this issue and in order to focus only on the real poisson equation please perform the test I mentionned previously. Choose a f function, discretize your poisson equation, apply dp/dn=0 every where, follow your procedure to fixe a P_ref. From the P field then obtained from your poisson solver compare the result of div(grad(P)) with f on each grid node (and especially on the node where you fixed P_ref on boundary face instead of dp/dn=0) Another question, which iterative linear system solver do you use too ? 

August 8, 2011, 07:49 

#19 
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 742
Rep Power: 19 
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: 742
Rep Power: 19 
Quote:
So you are correct about "your criteria to judge that "it works " or not is that you get a solution from your NavierStokes 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:
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:
Quote:
I use smoothed aggregation multigrid. Or BiCGStab preconditioned by smoothed aggregation multigrid. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Physical meaning of pressure in pressure eqn  dontknow  Main CFD Forum  2  August 30, 2010 10:12 
Constant velocity of the material  Sas  CFX  15  July 13, 2010 08:56 
TwoPhase Buoyant Flow Issue  Miguel Baritto  CFX  4  August 31, 2006 12:02 
pressure poisson equation on nonstaggered grids  abilge  Main CFD Forum  12  March 29, 2005 12:34 
pressure poisson equation!!  frederic felten  Main CFD Forum  1  November 17, 2000 04:17 