Solving the Poisson equation with Neumann Boundary Conditions - Finite Difference, BiCGSTAB

Question

I am trying to solve the Poisson equation in a rectangular domain using a finite difference scheme with a rectangular mesh.

I have happily generated the matrix system of equations Ax = b which is required to be solved, but when I try to impose Neumann boundary conditions (using ghost cells) of zero gradient I am running in to problems, I believe because the matrix is singular.

I know that there are several similar threads on this topic: Discrete Poisson Equation with Pure Neumann Boundary Conditions , Writing the Poisson equation finite-difference matrix with Neumann boundary conditions , Solving a linear equation system with pure Neumann condition but I cannot find how to fix my problem. Some people suggest that I should impose a Dirichlet condition at one point on the boundary, but when I try this I am getting a spike at this point.

I was using a Gauss-Seidel procedure to solve the system, and at each iteration I would re-set the solution to have an average of 0. However, as this is very slow I moved to using the Bicongugate Gradient Stabilised Method which is much quicker. However, I cannot see a way in which to do something similar in this case, so any help or suggestions would be welcome.

One reason I think this may be the case is that the test problems I am using may be unphysical. For the case where I had b being a sinusoidal forcing with an average of zero over the domain (-2 pi^2 cos(pi x)cos(pi y) from 0-1 in x and y), setting one corner to a Dirichlet condition appeared to work fine. However, for a point source in the centre of the domain, setting a Dirichlet condition at one corner caused a spike there. In my head/reasoning physically a point source in the centre of the domain with Neumann conditions is just going to keep 'pushing' and the solution will float off to infinity with an ever growing arbitrary constant and cannot be anchored. So there is some kind of condition that says I need zero net forcing by b on my closed domain must be zero i.e. the integral of the source term across the domain must be zero. Is this correct?

If it helps, the application of this is solving a Poisson equation for pressure in a fluid dynamics pressure correction scheme. Here I cannot guarantee that b will have an average of zero, though it should hopefully be close.

I am happy to provide more information as necessary.

Answer 1

One of the standard ways to handle singular problems like Poisson with Neumann boundary conditions which lead to a sinuglar problem $Ax=b$ is to make the obtained solution orthogonal to the kernel (which is for simple domains only the constant vector) if you use something like a Krylov subspace method.

For instance, you can solve your system with conjugated gradients (CG) and then make the solution orthogonal to the kernel, i.e., subtract the average from your solution. But before solving the linear system make sure that your righthand side vector is also orthogonal to the kernel (again, subtract the average of its values), which is the compatibility condition for singular systems.

The approach with 1-point Dirichlet condition is known to be bad because you spoil the underlying boundary value problem (which is singular, but nice) and one should usually avoid this way.