I am trying to write a multigrid solver for Poisson's equation, $-\Delta u=f$, on the unit square, $\Omega=(0,1)^2$ with periodic boundaries. My primary source has been Multigrid by Trottenberg, Oosterlee and Schuller, 2001. Based on their discussion I'm using the standard 2nd order central difference for the derivatives: $$-\Delta_h u(x_i,y_j) \approx \frac{1}{h^2} \left[ 4u(x_i,y_j) -u(x_{i+1},y_j) - u(x_{i-1},y_j)-u(x_i,y_{j+1})-u(x_i,y_{j-1}) \right]$$ where $$x_i=\frac{i}{h} \text{ } i =1,2,\dots ,n \text{ and } y_j=\frac{j}{h} \text{ } j=1,2,\dots , n$$ with the following boundary conditions: $$ u(x_0, y_j) = u(x_n, y_j)$$ $$ u(x_1, y_j) = u(x_{n+1}, y_j) $$ $$ u(x_i, y_0) = u(x_i,y_n)$$ $$ u(x_i, y_1) = u(x_i,y_{n+1})$$
For existence, I'm enforcing the compatibility condition by altering the right hand side as follows: $$f_h(x_i,y_j)=f(x_i,y_j)-h^2\sum_{i,j}f(x_i,y_j).$$ For uniqueness, I'm enforcing the global constraint: $$\sum_{i,j}u_h(x_i,y_j)=0.$$
Now, here is where my question comes in:
Before writing the multigrid solver, I wanted to test my implementation of the smoother (in this case I'm using Gauss-Seidel, although eventually I plan to switch to Red-Black Gauss Seidel). This gives me the following iterative scheme: $$u^{m+1}_{i,j}=\frac{h^2}{4}f_h(x_i,y_j)+\frac{1}{4}\left[ u^m_{i+1,j}+ u^{m+1}_{i-1,j} +u^m_{i,j+1}+u^{m+1}_{i,j-1} \right] \text{ }i,j=1,2,\dots n$$ where I have used the shorthand notation: $u_{i,j}=u(x_i,y_j)$. My first test case was: $$u(x,y)=\cos(2\pi x)\cos(2 \pi y),$$ giving us the source term: $$f(x,y)=8\pi^2\cos(2\pi x)\cos(2 \pi y).$$
The solution is periodic on $(0,1)^2$ and also satisfies the continuous compatibility condition:
$$\int_0^1\int_0^1 f(x,y) dx dy=0.$$
However, it doesn't satisfy the discrete compatibility condition It does saitsisfy the discrete compatibility condition (see edit):
$$ \sum_{i,j}f(x_i,y_j)=0,$$
and so I have imposed the condition above so that $\sum_{i,j}f_h(x_i,y_j)=0$.
What is strange to me is that when I solve the system $-\Delta_h u_{i,j}=f_h(x_i,y_j)$, the residual does not converge to zero (see edit). After noticing this, I removed the enforcement of the discrete compatibility condition and solved the system $-\Delta_h u_{i,j}=f(x_i,y_j)$. In this case the residual converges to zero. Since the correction to the RHS is $O(h^2)$ in both cases my numerical solution converges to the analytic solution as $n\to\infty$.
Since I plan to use a Multigrid method, I really want the residual to go to zero. For a practical problem I will only have data on the source term at discrete points, and thus, I will only know if the discrete compatibility condition is satisfied. Now for the question, based on my experiments it seems like my residual, $R_{i,j}$, should be: $$R_{i,j}=f(x_i,y_j)+\Delta_h u_{i,j}$$ where $f(x_i,y_j)$ is the source term without the compatibility correction.
This confuses me when I look at the mathematics, though. I am solving the system $-\Delta_h u_{i,j}=f_h(x_i,y_j)$ and thus, to me it seems like I should be checking how well my iterative scheme satisfies this equation. In other words, it seems like the residual should be: $$\bar{R}_{i,j}=f_h(x_i,y_j)+\Delta_h u_{i,j}.$$
Can anyone find the flaw in my reasoning?
EDIT: Problem solved. I was evaluating the compatibility condition incorrectly. The correct compatibility condition is: $$ \sum_{i=1}^n \sum_{j=1}^n f(x_i, y_j)=0 $$ NOT $$ \sum_{i=0}^n \sum_{j=0}^n f(x_i, y_j)=0. $$ This counts the boundary terms twice.
The correct residual does converge to zero and is given by: $$\bar{R}_{i,j}=f_h(x_i,y_j)+\Delta_h u_{i,j}.$$
