I've been encountering some convergence problems with the technique described in section 19.4 for inhomogenous Dirichlet conditions (which essentially amounts to taking the known values over to the RHS).

The rms error initially goes as h^2 as I refine the grid, but at about 32x32 drops to h^0.5, and the residual actually increases (h^-1.5). (This is with the FAS routine that automatically terminates when it thinks truncation error is reached; it terminates after only 1 V-cycle, by the way.) If I add enough V-cycles then I get good convergence, so it would appear that the algebraic problem is correct, but not of high enough order.

I've also tried extra relaxation sweeps (up to 4 pre- and 4-post relaxations) and W-cycles, but none of it really seems to help much (it basically puts off the slowdown in convergence, but it eventually happens). If I simply apply the boundary conditions at every level (and don't modify the RHS) then everything works fine, even with plain old V-cycles and 1 pre- and post-relaxation sweep.

A plot of the residual clearly shows that it is smoothed successfully in the interior, but not on the boundaries.

I've been over and over my code, and can't find anything wrong. (In fact, just to be sure, I created the u_B array explicitly, then applied the operator to it and subtracted the resulting array from the RHS, which is inefficient, but fairly idiot proof.)

Does anyone have any idea what might be wrong? Is it my code, or is the technique described in 19.4 just not suited to multigrid? If it's my code, any ideas for what to try?

The test problem I am using, by the way, is just
del^2 u = rho
with the solution u =sin(pi x) * sin(2 pi y) + 9
(and boundary condition u_B = 9, obviously).