# How do I avoid divide-by-zero when solving the Poisson equation with Fourier transforms?

I wanted to try to implement part of the method in the following article using Fourier transforms.

http://www.shodor.org/media/content/jocse/student_submissions/nocito2010/nocito2010_pdf

Right now I am using the FFT to approximate the gradients at each point and using this to calculate the force on each body. I am currently getting strange results. When I add just one massive body to the system I get the following as my potential which seems incorrect. I am aware that the nature of the Fourier transform causes the wrap around which would distort the spherical shape but I would still not expect this.

Here is my code. One issue I am having is that in the paper you are supposed to divide by |k|^2. However at some points there is a zero which causes issues. That is why I clamp it to 0.1. Not sure if there is a better way to handle that.

    #include <iostream>
#include <arrayfire.h>
#include <memory>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <cstring>

const int N = 2;
const int GRID_DIM = 1024;
const double CELL_SIZE = 100.0;
const double G = 1.0;
const double DT = 0.001;

struct Vec2
{
float x, y;
Vec2(float x = 0, float y = 0) :
x(x), y(y) { }
};

int wrapIndex(int i, int i_max)
{
return ((i % i_max) + i_max) % i_max;
}

af::array makeGrid(const Vec2* p, const float* m)
{
float* grid_h = new float[GRID_DIM * GRID_DIM];
memset(grid_h, 0, sizeof(float) * GRID_DIM * GRID_DIM);

for(int i = 0; i < N; i++)
{
int x = wrapIndex(roundf(p[i].x / CELL_SIZE), GRID_DIM);
int y = wrapIndex(roundf(p[i].y / CELL_SIZE), GRID_DIM);
grid_h[y * GRID_DIM + x] += m[i];
}

af::array grid(GRID_DIM, GRID_DIM, grid_h);
delete[] grid_h;
return grid;
}

int main()
{
af::setBackend(AF_BACKEND_CUDA);
af::info();
Vec2* p = new Vec2[N];
Vec2* v = new Vec2[N];
float* m = new float[N];

std::fill(v, v + N, 0.0f);
std::fill(m, m + N, 1.0f);
//  for(int i = 0; i < N; i++)
//      p[i] = Vec2(rand() % GRID_DIM / 5 + GRID_DIM / 2, rand() % GRID_DIM / 5 + GRID_DIM / 2);

//add one massive object at the center
m[0] = 1000000.0;
p[0] = Vec2(GRID_DIM * CELL_SIZE / 2, GRID_DIM * CELL_SIZE / 2);

m[1] = 0.0;
p[1] = Vec2(GRID_DIM * CELL_SIZE / 4, GRID_DIM * CELL_SIZE / 4);

af::array rng = af::range(GRID_DIM);
af::array k2 = 2.0 * af::Pi * af::select(rng > GRID_DIM / 2, (rng - GRID_DIM) , rng) / (GRID_DIM * CELL_SIZE);
//af::array k2 = 2.0 * af::Pi * rng / (GRID_DIM * CELL_SIZE);
k2(0) = 1.0;
k2 = af::tile((k2 * k2).T(), GRID_DIM) + af::tile(k2 * k2, 1, GRID_DIM);

af::Window wnd(GRID_DIM, GRID_DIM);

af::array grid = makeGrid(p, m);
while(!wnd.close())
{
af::timer t = af::timer::start();

//fft the density grid
grid = af::fft2(grid);

//apply the formula from the paper
grid *= -G / (af::Pi * k2);
//invert the fft to get the potentials
grid = af::ifft2(grid);
grid = af::abs(grid);

//compute the gradients in the potential
af::array dx, dy;

float* dx_h = dx.host<float>();
float* dy_h = dy.host<float>();

//use the gradients to time step each of the bodies
#pragma omp parallel for simd
for(int i = 0; i < N; i++)
{
int xi = wrapIndex(roundf(p[i].x / CELL_SIZE), GRID_DIM);
int yi = wrapIndex(roundf(p[i].y / CELL_SIZE), GRID_DIM);

v[i].x += dx_h[yi * GRID_DIM + xi] * DT / m[i];
v[i].y += dy_h[yi * GRID_DIM + xi] * DT / m[i];
p[i].x += v[i].x * DT;
p[i].y += v[i].y * DT;
}
af::freeHost(dx_h);
af::freeHost(dy_h);

//display the potential

wnd.image(af::abs(grid) / af::max<float>(af::abs(grid)));
//remake the new grid
grid = makeGrid(p, m);

std::cout << "dt = " << af::timer::stop(t) << std::endl;

}

delete[] p;
delete[] v;
delete[] m;
return 0;
}


Here is another less saturated photo of the potential for a single point. I can't figure out why it is not approximately round / radially symmetric.

You are solving Poisson's equation:

$$\nabla^2 \Phi = 4\pi G \rho.$$

Notice that if $\Phi$ is a solution, then so is $\Phi+C$ for any constant $C$. Furthermore, $C$ will have no effect on your simulation since the forces depend only on derivatives of $\Phi$.

The divide-by-zero issue is just a manifestation of this degeneracy. In the Fourier-transformed Poisson equation

$$\hat{\Phi}_k = \frac{-G}{\pi |k|^2} \hat{\rho}_k,$$

the $k=0$ mode corresponds to the constant part of $\Phi$. In practice you can simply discard it.

• How do I go about discarding it? Should I just zero out the zero'th row and col in the transform? – chasep255 Apr 8 '16 at 4:22
• What you say you are doing now is also fine (I haven't checked the code, just your claims). The point is that $\hat{\Phi}_0$ will get multiplied by zero anyway when you compute the forces. – David Ketcheson Apr 8 '16 at 4:26
• Actually, I should amend my statement. If you impose a lower bound on $|k|$ you need to make sure that it is only applied to the zero mode! In your case I think you may be setting many $k$ values to your floor value, which will give meaningless results. – David Ketcheson Apr 8 '16 at 4:28
• I just updated my code. I seem to still be getting a shape with strange symmetry which looks nothing like the potential I would expect from a point source. – chasep255 Apr 8 '16 at 4:38
• Well I got it. I was setting k(0) = 1.0 before I put it into the 2D grid so it was "corrupting" the entire first row and column. Now it seems to work great. – chasep255 Apr 8 '16 at 5:32