Linear Algebra / Numerical Solution Of Matrix With Nullspace

Question

I have a question relating to linear algebra.

We have a fluid solver that solves the poisson equation for pressures. When there are areas of the domain that are entirely enclosed by Neumann conditions, the matrix has a nullspace and we must correct for this. We keep around a mask of the enclosed regions, and in each iteration compute the average value and subtract it off.

We've noticed that in this case we can deduce that mask from the matrix itself, rather than from the domain. In that way we have a sort of modified numerical method that operates only on the matrix $A$ and the RHS $b$, but nevertheless converges when $A$ has a nullspace (as opposed to the first paragraph where the inputs to the method are $A$, $b$ and the enclosed regions mask).

I'm wondering if this generalises a little more, or if anybody has thought about this. For example, is there a modified Jacobi method that always converges, even if the input matrix has a nullspace, and where the only inputs are $A$ and $b$? The solution we are interested in is the one where each Neumann-enclosed region has the smallest average value.

I'm interested because our domains aren't always so simple; sometimes we don't have Neumann conditions but other conditions on the pressure that nevertheless leave $A$ with a nullspace, and it seems more robust from a software point of view to be calculating the nullspace correction from $A$ itself.

(Normally the answers to a stack overflow question metioning Jacobi are why are you using Jacobi??! FYI our full solver is a multigrid pre-conditioned conjugate gradient, where we use Jacobi and Gauss-Seidel as smoothers for the multigrid (A. McAdams , E. Sifakis , J. Teran, 2010).)

Thanks!

Answer 1

This is a pretty common problem. There are two simpler ways to what you are doing:

• Enforce a Dirichlet-kind condition on one pressure node by fixing it to zero. This way, your matrix no longer has a null space. The pressure may not have mean value zero, but you can fix this once (instead of once per iteration) after you have the solution of the linear system by subtracting the mean value.

• Many iterative Krylov-space solvers will converge even with a null space if your initial residual has no component in this null space -- in other words, if the pressure component of the initial residual has mean value zero. You then do not need to do anything at all during the iterations. The solution process is not going to be stable within this null space but since you don't care about it, if you accumulate a small nonzero mean value, you can always subtract it at the end. The difficulty is that while the solver will ignore the null space, you need to have a preconditioner that, when given a vector with mean value zero pressure produces a vector with mean value zero pressure.