Convergence rate of FFT Poisson solver

Question

What is the theoretical convergence rate for an FFT Poison solver?

I am solving a Poisson equation: $$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$ with $$n(x, y, z) = {3\over\pi} ((x-1)^2 + (y-1)^2 + (z-1)^2 - 1)$$ on the domain $[0, 2] \times [0, 2] \times [0, 2]$ with periodic boundary condition. This charge density is net neutral. The solution is given by: $$ V_H({\bf x}) = \int {n({\bf y})\over |{\bf x}-{\bf y}|} d^3 y $$ where ${\bf x}=(x, y, z)$. In reciprocal space $$ V_H({\bf G}) = 4\pi {n({\bf G})\over G^2} $$ where ${\bf G}$ are the reciprocal space vectors. I am interested in the Hartree energy: $$ E_H = {1\over 2} \int {n({\bf x}) n({\bf y})\over |{\bf x}-{\bf y}|} d^3 x d^3 y = {1\over 2} \int V_H({\bf x}) n({\bf x}) d^3 x $$ In reciprocal space this becomes (after discretization): $$ E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2} $$ The ${\bf G}=0$ term is omitted, which effectively makes the charge density net neutral (and since it is already neutral, then everything is consistent).

For the test problem above, this can be evaluated analytically and one gets: $$ E_H = {128\over 35\pi} = 1.16410... $$ How fast should this energy converge?

Here is a program using NumPy that does the calculation.

from numpy import empty, pi, meshgrid, linspace, sum
from numpy.fft import fftn, fftfreq
E_exact = 128/(35*pi)
print "Hartree Energy (exact):      %.15f" % E_exact
f = open("conv.txt", "w")
for N in range(3, 384, 10):
    print "N =", N
    L = 2.
    x1d = linspace(0, L, N)
    x, y, z = meshgrid(x1d, x1d, x1d)

    nr = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
    ng = fftn(nr) / N**3

    G1d = N * fftfreq(N) * 2*pi/L
    kx, ky, kz = meshgrid(G1d, G1d, G1d)
    G2 = kx**2+ky**2+kz**2
    G2[0, 0, 0] = 1  # omit the G=0 term

    tmp = 2*pi*abs(ng)**2 / G2
    tmp[0, 0, 0] = 0  # omit the G=0 term
    E = sum(tmp) * L**3
    print "Hartree Energy (calculated): %.15f" % E
    f.write("%d %.15f\n" % (N, E))
f.close()

And here is a convergence graph (just plotting the conv.txt from the above script, here is a notebook that does it if you want to play with this yourself):

As you can see, the convergence is linear, which was a surprise for me, I thought that FFT converges much faster than that.

Update:

The solution has a cusp at the boundary (I didn't realize this before). In order for FFT to converge fast, the solution must have all derivatives smooth. So I also tried the following right hand side:

nr = 3*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)/4

Which you can just put into the script above (updated script). The exact solution is $V_H=\sin(\pi x)\sin(\pi y)\sin(\pi z)$, which should be infinitely differentiable. The exact integral in this case is $E_H = \frac{3\pi}{8}$. Yet the FFT solver still converges only linearly towards this exact solution, as can be checked by running the script above and plotting the convergence (updated notebook with plots).

Does anyone know any benchmark in 3D so that I can see faster convergence than linear?

Answer 1

Let me first answer all the questions:

What is the theoretical convergence rate for an FFT Poison solver?

The theoretical convergence is exponential as long as the solution is sufficiently smooth.

How fast should this energy converge?

The Hartree energy $E_H$ should converge exponentially for a sufficiently smooth solution. If the solution is less smooth, then the convergence is slower.

Does anyone know any benchmark in 3D so that I can see faster convergence than linear?

Any right hand side that produces a solution that is periodic and infinitely differentiable (including across the periodic boundary) should converge exponentially.

---

In the code above there happens to be a bug, that causes the convergence to be slower than exponential. Using the smooth density code (https://gist.github.com/certik/5549650/), the following patch fixes the bug:

@@ -6,7 +6,7 @@ f = open("conv.txt", "w")
 for N in range(3, 180, 10):
     print "N =", N
     L = 2.
-    x1d = linspace(0, L, N)
+    x1d = linspace(0, L, N+1)[:-1]
     x, y, z = meshgrid(x1d, x1d, x1d)

     nr = 3*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)/4

The problem was that the real space sampling cannot repeat the first and last point (which has the same value due to periodic boundary condition). In other words, the problem was in setting up the initial sampling.

After this fix, the density converges in one iteration, as Matt said above. So I didn't even plot the convergence graphs.

However, one can try a more difficult density, for example:

     nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4

then the convergence is exponential, as expected. Here are the convergence graphs for this density: