What could be the reason that the solution of a Poisson equation is smooth when obtained by an iterative solver, only if the maximum residual is set to a high value (e.g. 0.1)? When the maximum residual is set as a typical lower value (e.g. 1e-6), then the solution oscillates - often (but not exclusively) near Neumann boundary conditions. The solution with the high residual is underpredicted, and with the low residual it oscillates around expected values.
The matrix is non-symmetric sparse, weakly diagonally-dominant (a diagonal entry equals the negative sum of other row entries). The system includes both Dirichlet and Neumann boundary conditions. Similar outcomes happen using various iterative solvers: BiCGStab, CGS, QMR and GMRES. The preconditioner is always a diagonal preconditioner, since I have to use matrix-free solving approach on the GPU. Using no preconditioner yields similar (though a bit smoother) results.