Pressure Poisson Equation (PPE)
 

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.

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.

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.

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.cfd-online.com/Forums/mai...-equation.html

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.

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.

Well to my knowledge this is how we could solve all Neumann kind of problem since there is no unique solution.


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.

Think about it.

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.

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.

Dear arjun,

Thank you very much for the reference, I was following for such book.

Ghasemi.

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

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.

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

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.

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.

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.

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.
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.

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 ?

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.

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.

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.


It seem so now to me too.



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