Hello,

I'm trying to solve the pressure Poisson equation with FULL Neumann boundary condition. Unfortunately, it leads to the sparse linear system equation Ax = b, where A is singular (because the Neumann BC in all boundaries). Well, I've tried Gauss-Seidel and it works very well, but it is very very slow. Now, I'm trying to solve this equation with GMRES or BiCGStab with the library gmm++, but it seems that all Krylov methods get blocked with singular matrixes and do not converge. I'd like to know how can I solve the pressure Poisson equation with GMRES? Is there any library which implements these methods with possibility to solve singular linear systems?

You need to remove the null space. The best way is to tell the Krylov solver about this (see KSPSetNullSpace(...) in PETSc) but you can also fix the pressure at one point. That is, pick any row of your matrix and set it to a row of the identity and put 0 on the right hand side. This should also make the matrix nonsingular, but it disrupts the spectrum so may effect convergence rates (highly preconditioner dependent).

What can you talk about the difference of PETSc and AztecOO in Trilinos package? It seems that PETSc has more options when solving a linear system.

I'm not very familiar with AztecOO but there should be a way to remove the null space. I appreciate the PETSc philosophy of making everything available from the command line. Avoid at all costs writing code that you have to recompile to change solvers or preconditioners. I think PETSc is somewhat more comprehensive and it provides many third-party preconditioners (including ML which is arguably the best part of Trilinos, and Hypre from LLNL) and direct solvers all of which are tunable from the command line. Also PETSc cleanly separates the preconditioning matrix from the Jacobian. This is especially useful for me because I have a fast way to apply the Jacobian matrix-free, but use different (sparser) matrices for preconditioning. This is possible with AztecOO, but less clean in my opinion. The matrix and vector interface is certainly much richer in PETSc which makes it possible to write solvers which are more generic.

my $.02

Allan

I would image that the PPE formulation would have neumann BC for closed cell problems. but what sorts of problems are you designing this code for where you anticipate all the pressure boundaries to be gradient (Neumann)? Also, I haven't experience with the CG solver you describe but shouldn't it solve for N BC if the Neumann condition is upheld?

curious.

Regards, D

The case is that in the development of the numerical procedure for solving the incompressible Navier-Stokes equation, you end up with the Neumann BC for all boundaries in a Poisson equation for a variable phi (which is not the pressure, but is related to it).

laplacian(phi) = divergence(V*)/dt, interior nabla(phi) = 0, boundary

If you wish, I can cite to you the paper where I got it.

I've found it in the PETSc's manual:

4.5 Solving Singular Systems

Sometimes one is required to solver linear systems that are singular. That is systems with the matrix has a null space. For example, the discretization of the Laplacian operator with Neumann boundary conditions as a null space of the constant functions. PETSc has tools to help solve these systems. First, one must know what the null space is and store it using an orthonormal basis in an array of PETSc Vecs. (The constant functions can be handled separately, since they are such a common case). Create a MatNullSpace object with the command MatNullSpaceCreate(MPI Comm,PetscTruth hasconstants,int dim,Vec *basis,MatNullSpace *nsp); Here dim is the number of vectors in basis and hasconstants indicates if the null space contains the constant functions. (If the null space contains the constant functions you do not need to include it in the basis vectors you provide). One then tells the KSP object you are using what the null space is with the call KSPSetNullSpace(KSP ksp,MatNullSpace nsp); The PETSc solvers will now handle the null space during the solution process. But if one chooses a direct solver (or an incomplete factorization) it may still detect a zero pivot. You can run with the additional options -pc_factor_shift_nonzero <dampingfactor> or -pc_ factor_shift_nonzero <dampingfactor> to prevent the zero pivot. A good choice for the damping factor is 1.e-10.

Hi Allan

But if you have an outflow boundary, you would normally set the PPE boundary as fixed Dirichlet. Gradient BC on vixed velocity faces derive P,n = mu Delta(V.n). Thus you would automatically satisfy the Neumann condition. Essentially, by analogy, the diffusion of pressure (via the PPE equation) has access to this fixed boundary thereby eliminating the non-uniqueness I think you are talking about.

Actually BiCGstab should still work without fixing the unknown constant in the pressure - in practise it actually converges faster when the constant isn't fixed.

Hello Q,

In my case, if there is an outflow BC, I still have to set up the Neumann BC in all boundaries for the Pressure Poisson equation. This is because of an enforcement of the BC on the Momentum equation (without pressure). The treatment of the outflow goes in other place.

1º) Calculate u* and v* (without pressure [p] in the momentum eq.) with Crank-Nicholson for viscous term and Adams-Bashforth for advective term. The boundary conditions are u^(n+1) and v^(n+1).

2º) Calculate phi^(n+1) with laplacian(phi) = divergence(u*, v*) / dt. Boundary conditions must be Neumann in all boundaries.

3º) Calculate u^(n+1) and v^(n+1) with the projection eq. V^(n+1) = V* - dt*Grad(phi)

You can find all that I've written here in: http://citeseerx.ist.psu.edu/viewdoc...=10.1.1.7.5949

This is an excellent paper.

[shyamdsundar]

Hi Allan,
I am wondering if you were able to solve the problem you had mentioned in the earlier posts. I am trying to use GMRES based solver for solution of pressure poisson equation.

Please let me know if you managed to solve the problem.

Thanks and Regards,
Shyam

[camaosos]

The problem maybe is that the "pressure" that enters is not the same as in the boundaries gets out. If you have for example a Laplace equation (laplacian(P) = 0) it says that "not any pressure inside is being generated". The same as a thermic wall (L=1) with no generation (the heat that enters in one side SHOULD BE THE SAME as in the other side). And if you have 0 flux in the left side and 1 in the right, the only Poisson 1D equation that can be solved is d^2 P /dx^2 = 1 (because we have to generate 1 "heat" to receive 1 heat in the right side).

Whe i was simulating that "singular system" in 2D i always got Planes wherever the f (right side) was. And the GMRES converged only when the generation and the Neumman BC's was in "equilibrium".