I have written Poisson solvers using two different methods:

A classic Jacobi scheme and one using the multigrid solver Hypre. I made up a couple of test cases ensuring the validity of those solvers.

For both cases the domain is defined as $(x,y) \in [-1,1]^2$ with periodic boundary conditions. Also, the grids first and last point are the same:

$$p(0,y) = p(N_x-1,y)$$

$$p(x,0) = p(x,N_y-1)$$

Test case 1

  • $f_{rhs} = -8 \pi^2 \sin(2\pi x) \sin(2\pi y)$
  • $p_{exact} = \sin(2\pi x) \sin(2\pi y)$

For both solvers, the solution is 2$^{\mathrm{nd}}$ order accurate is space. No problems here so far.

Test Case 2

  • $p_{rhs} = e^{-10 (x^2 + y^2)}$
  • $p_{exact}:$ No analytical solution, and therefore the numerical solution is differentiated using a high-order compact scheme and compared to $p_{rhs}$

Note that in this case, $\int_V p_{rhs}dV \neq 0$ and therefore the problem is ill defined. Therefore, the $rhs$ must be modified:

$p_{rhs} = e^{-10 (x^2 + y^2)} - \dfrac{\int_D e^{-10 (x^2 + y^2)} dx dy}{V}$

where $V$ is the domain volume. The integral is computed using the Trapezoidal rule.

This is where things get tricky. No matter how fine my grid is, $\left(p_{num}\right)_{xx} + \left(p_{num}\right)_{yy}$ never converge to $p_{rhs}$. When the grid is fine enough, the solvers converge, but the solution is off by ~20% while the overall profiles are relatively correct. When the grid is coarse, the Hypre solver simply diverges.


Have I missed anything? Is my approach inconsistent/wrong?

  • 2
    You say that the solution is unique up to a constant, but I think there exist no solution. But if you correct it to be zero on average, it should exist (and be determined up to a constant). – Vladimir F Aug 9 '17 at 15:14
  • @VladimirF Yes you are right. Fixed. – solalito Aug 9 '17 at 15:40
  • 1
    I am not sure you can use the trapezoidal integration to correct the p_rhs. You need the RHS of the discretized linear system to be zero on average. Is it? – Vladimir F Aug 10 '17 at 13:36
  • 1
    Fix one node, could be anywhere, for example on corner, and set value at that node to your exact solution, you will get well posed problem. – likask Aug 10 '17 at 19:51
  • Please skip the trapezoidal rule and compute the integral by yourself to the exact value. There is no need to introduce an additional approximation error. Hint: the integral is the product of two one-dimensional integrals. – shuhalo Aug 11 '17 at 17:23

In your Test Case 1, the solution is periodic. Hence, you had no problems when you imposed periodic boundary conditions.

For the Test Case 2, the solution may not be periodic (I suspect that it is not). By enforcing the periodic conditions, you get some solution which does not actually corresponding to the source term. That is why your solution is off by 20%. This should be due to incompatibility between the actual solution and the obtained solution at the boundaries. If you check carefully, you may find that the maximum error occurs at the boundaries.

You can verify your code with another benchmark example from FEniCS. https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/periodic/python/documentation.html

  • Why shouldn't there exist a periodic solution? – shuhalo Aug 11 '17 at 17:15
  • 1
    In order to answer further, I need to know the value you impose on the boundary for the Test Case 2. – Chenna K Aug 11 '17 at 21:56
  • There are periodic boundary conditions. Nothing to impose on the boundary. – shuhalo Aug 12 '17 at 4:07
  • Even for periodic boundaries, you have to specify the boundary values on one of the boundaries in each direction. However, I don't think it really matters whether you use periodic boundary conditions or not. As the Jacobi Poisson solver is used here, it is assumed that the values of all the boundary points are known. If you don't specify anything, then the value is assumed to be zero. – Chenna K Aug 12 '17 at 14:08
  • I tested these two problems using Jacobi solver. For 40x40 grid, Case 1 requires 167 iterations, while Case 2 takes 2674 iterations (~16 times more iterations). I suspect that this is due to incompatible boundary conditions. Even if the periodic solution exists for Case 2, it might not vanish on the boundaries of the domain considered for this problem. So, by enforcing periodic bcs with zero value everywhere on the boundary, you are asking the solver to solve for some unnatural (to the problem) conditions. – Chenna K Aug 12 '17 at 14:12

So after much testing, I have found a solution to this problem. However, I am not quite sure of the mathematical justification for it.

Basically, when using purely periodic BC, iterative solvers would not converge and the solution will keep growing. By fixing Dirichlet BC on all sides of the domain to a fixed constant $C \in \mathbb{R}$, on top of using periodicity, the RHS term can be obtained by differentiating the solution with an order of accuracy of 2 (note that high-order scheme are using for numerical differentiation to ensure the leading error term comes from the Poisson solver).

Furthermore, this got me thinking that there was no need for periodic BC to begin with since only the gradient of the solution is of interest to me:$\; \dfrac{\partial p}{\partial x_i}$. By just enforcing dirichlet BC on all sides, this gradient can be retrieved with the desire order of accuracy. I will be happy to hear more of why that is though.

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.