15
$\begingroup$

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):

FFT convergence graph

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?

$\endgroup$
  • $\begingroup$ Ondrej, isn't the Fourier transform of your smooth density a delta function? I admit to being too lazy to run it, but it should get the exact answer on the first try. $\endgroup$ – Matt Knepley May 9 '13 at 22:25
  • $\begingroup$ I think it is. But it doesn't converge in one iteration, as can be seen from the notebook plots. I have no idea what is going on. $\endgroup$ – Ondřej Čertík May 10 '13 at 16:34
  • $\begingroup$ Ondrej, are you sure your implementation is correct? I remember trying to implement spectral solvers for a homework assignment in grad school and totally flubbing the constants. I notice that you are measuring error by looking at the absolute distance between the computed and exact energy. What does your convergence look like to the problem's actual solution? This should be easy to compute and even plot over a 2-D slice of the problem. $\endgroup$ – Aron Ahmadia May 11 '13 at 23:41
  • $\begingroup$ Aron --- I checked my implementation against some other code, but I was checking it for my wrong initial sampling, so I had the same bug in both codes. Matt was right, now it converges on the first try. See my answer below. $\endgroup$ – Ondřej Čertík May 12 '13 at 5:05
9
$\begingroup$

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: enter image description here enter image description here

$\endgroup$
  • $\begingroup$ Awesome. Sorry I wasn't more help! $\endgroup$ – Aron Ahmadia May 12 '13 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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