Practical implementation of the discrete compatibility condtion

Question

I've recently started looking into writing a finite-differences-based solver for the Poisson equation of the form $$\nabla\left(\varepsilon\nabla\varphi\right)=\rho$$ in 2D for arbitrary geometries (some additional details can be found here). However, I've soon discovered that upon imposing pure Neumann boundaries the system no longer can be solved in a straight forward manner, as the problem becomes ill-posed. First I thought the failure to solve the system arises perhaps due to some kind of a limitation of the solver I use (successive relaxation method), but from answer here I learned that the compatibility conditions are not satisfied in my problem (in particular the discrete compatibility condition), on top of the fact that the all Neumann boundaries make it an ill-posed problem (Edit: not sure about the nomenclature here, but the lack of satisfied compatibility condition could be classified as an ill-conditioned problem).

Now, how I understood one can impose the compatibility condition is just to modify the right-hand side (RHS) of the Poisson equation (including both the original $\rho$ and the Neumann boundary conditions that are stored in the RHS vector) by subtracting its mean as: $$\rho_{h}=\rho-mean(\rho),$$ where mean is calculated over all points in the vector $\rho$. This should be trivial enough, right? Well, my implementation of this in C++ did not yield any results and the convergence issues continued. In fact, the solvers I tried out with this (SOR, biCGstat, CG, CGsquare) all failed to converge as their respective residuals either diverged or just get stuck at a particular value and looped around it until the iteration reached maximal allowed iteration count.

I have recently encountered this post, where a user complains about a similar issue. He, luckily, managed to resolve it, but did not go too much into the details (he only explained that something was incorrect with the counter for the mean value calculation/subtraction).

Thus I wanted to ask if somebody can explain the practical implementation details of the discrete compatibility condition. What needs to be taken into account when calculating the average of the RHS vector, and what needs to be taken into account when subtracting the said average from the original RHS?

Answer 1

It is actually not that complicated, try to think less analytically.

Assume you have already given the ill-posed & ill-conditioned system $A u=f$ of the original Poisson Neumann problem ($\varphi\equiv u,~\rho\equiv f, \epsilon =1).$ All you have to do is to introduce one single constraint

\begin{align} \frac{1}{b-a}\int_a^b u \,\mathrm{d}x&=0, \\ \int_a^b u \,\mathrm{d}x&=0, \end{align}

which can be discretisized with the most simple quadrature rule (rectangle rule)

\begin{align} \sum_{i=1}^N h u_i &=0, \\ u_1 +u_2 + \dots + u_N &= 0 \equiv \lambda. \end{align}

Now add this single line to your system

$$\begin{bmatrix} A \\ 1^T \end{bmatrix} \begin{bmatrix} u \end{bmatrix} = \begin{bmatrix} f \\ 0\end{bmatrix}.$$

At this point your system is still irregular because you have more lines than columns. In order to fix this add a "clever" zero using $\lambda$ resulting in

$$\begin{bmatrix} A & 1 \\ 1^T & 0 \end{bmatrix} \begin{bmatrix} u \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}$$

which is invertible and solvable by most standard solvers.

Edit: Credits to @Will P. which however is not general at all.

The Poisson equation with Neumann BCs is given by

\begin{align} \Delta u(x) &= f(x),\qquad x\in\Omega\\ \frac{\partial u(x)}{\partial\mathbf n} &= g(x),\qquad x\in\partial\Omega. \end{align}

To actually be well posed the divergence theorem has to be fullfilled

\begin{align} \int_\Omega f(x)\,dV = \int_\Omega \operatorname{div}\nabla u(x) \,dV = \oint_{\partial\Omega}\nabla u(x)\cdot \mathbf n(x)\,dS = \oint_{\partial\Omega} g(x)\,dS, \end{align}

which ultimately results in following condition

\begin{align} \int_\Omega f(x)\,dV = \oint_{\partial\Omega} g(x)\,dS. \end{align}

Note that the answer by @Will P. only covers the case $\oint_{\partial\Omega} g(x) = 0$.

Answer 2

The problem $$ \nabla (\varepsilon \nabla u) = \rho $$ with pure Neumann conditions admits a family of solutions that differ by a constant value (i.e. $u + c$ is also a solution, for any constant $c$).

For most numerical methods, after discretizing, you will obtain a linear system $$ L \boldsymbol{u} = \boldsymbol{\rho}, $$ where $L$ is a symmetric matrix with a one-dimensional nullspace spanned by the constants, i.e. $$ L \boldsymbol{1} = \boldsymbol{0}, $$ where $\boldsymbol{1}$ is the vector with 1 in every entry.

In order for the discrete problem to be well posed, you need $\boldsymbol{\rho}$ to be in the range of $L$. Note that for any $\boldsymbol{y} \in \operatorname{Range}(L)$ (i.e. with $\boldsymbol{y} = L \boldsymbol{x}$, we have $$ \boldsymbol{1}^T \boldsymbol{y} = 0, $$ since, by symmetry of $L$, $$ \boldsymbol{1}^T \boldsymbol{y} = \boldsymbol{1}^T L \boldsymbol{x} = \boldsymbol{1}^T L^T \boldsymbol{x} = \boldsymbol{0}^T \boldsymbol{x} = 0. $$

Therefore, the discrete compatibility condition is that the vector $\boldsymbol{1}$ should be orthogonal to your right-hand side.

Assume you are given $\tilde{\boldsymbol{\rho}}$, which may not satisfy this condition. You can set $$ \boldsymbol{\rho} = \tilde{\boldsymbol{\rho}} - \left( \frac{\boldsymbol{1}^T \tilde{\boldsymbol{\rho}}}{n} \right) \boldsymbol{1}, $$ where $n$ is the dimension of the discrete space (i.e. $L \in \mathbb{R}^{n \times n}$). This is what is meant by subtracting off the mean of $\tilde{\boldsymbol{\rho}}$. It is easy to check that this modified vector $\boldsymbol{\rho}$ satisfies the discrete compatibility condition. Then, the linear system $$ L \boldsymbol{u} = \boldsymbol{\rho} $$ will have a solution.

The conjugate gradient method applied to this problem will converge to a solution without performing any special modifications (see this answer).

Note that if you are using preconditioned conjugate gradient, you may need to ensure that the output of the preconditioner remains orthogonal to $\boldsymbol{1}$. This can be achieved by performing an orthogonalization step every preconditioner application (i.e, if $M$ is the preconditioner, first compute $\tilde{\boldsymbol{z}} = M^{-1} \boldsymbol{r}$, and then set $\boldsymbol{z} = \tilde{\boldsymbol{z}} - \left( \frac{\boldsymbol{1}^T \tilde{\boldsymbol{z}}}{n} \right) \boldsymbol{1}$). If you are not using a preconditioner, this is not necessary.