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


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.

enter image description here

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()
    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);

    //load the grid with bodies
    af::array grid = makeGrid(p, m);
        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;
        af::grad(dx, dy, grid);

        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;

        //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.

enter image description here


1 Answer 1


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.

  • $\begingroup$ How do I go about discarding it? Should I just zero out the zero'th row and col in the transform? $\endgroup$
    – chasep255
    Apr 8, 2016 at 4:22
  • $\begingroup$ 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. $\endgroup$ Apr 8, 2016 at 4:26
  • $\begingroup$ 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. $\endgroup$ Apr 8, 2016 at 4:28
  • $\begingroup$ 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. $\endgroup$
    – chasep255
    Apr 8, 2016 at 4:38
  • $\begingroup$ 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. $\endgroup$
    – chasep255
    Apr 8, 2016 at 5:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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