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