# Particle Mesh Method - FFT with Green Function

1. Feb 25, 2014

### fab13

Hello,

For generating initial velocities on a N-Body code (modeling galaxy dynamics), I need to solve the Poisson equation with Green's function by applying Fast Fourier Transform.

The Fourier analysis being more straightforwardly performed in a periodic system, I use the "zero padding trick" (Hockney & Eastwood 1981) for a N x N grid. Let's consider the 2D case :

1. I compute the Green's function in real space, in the N x N grid

2. I duplicate and mirror the Green's function in the 3 "ghost" quadrants, so that the Green's function is now periodic in the Ng = 2N x 2N grid (see below figure)

3. Set the density field to zero in the 3 ghosts quadrants. (see below figure)

4. Perform the Fourier analysis on the 2N x 2N grid : multiply the Fourier transform of Green's function by the Fourier transform of density field ( product wich is equal to the Fourier transform of the gravitational potential)

5. Perform the inverse Fourier transform of the Fourier potential to get the potential in space.

In the original code on wich I am working, the author uses DCT (discrete cosine transform) of Green's function and density field. The DCT is just a particular case of the Discrete Fourier transform, i.e when the function has an even symmetry and is periodic.

So, why the DCT is used for computing the Fourier transform of density field (which is not periodic and even symmetric) ?

Actually, the grid has 3D dimension. The author of this code is printing into a data file the values of the potential on the grid at z=0. Once the inversion Fourier transform of the potential has been got, here the sample code (with Ng = 2 * N) for printing values :

Code (Text):

for (i = Ng/4; i < 3*Ng/4; i++) {
for (j = Ng/4; j < 3*Ng/4; j++) {
// fprintf left bottom on right top
fprintf(file,"%lf %lf %le \n",(double) ((i-Ng/2)*space_x), (double) ((j-Ng/2)*space_y), gal->potential[i][j][Ng/2]);
}
}

I don't understand why the author takes the interval [Ng/4, 3*Ng/4] for x (i index) and y (j index). Is the inverse Fourier transform of potential symmetric ? I don't think so because the Fourier transform of density field is not symmetric ( the 3 ghost quadrants are set to zero).

Moreover, why the z = 0 coordinate is got by taking gal->potential[j][Ng/2] (k index is set to Ng/2 for having z = 0) ?

If anyone has already implemented this method, could you help me please.

Last edited: Feb 25, 2014