Particle Mesh Method - FFT with Green Function

[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 which 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 which 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 :

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.

 

[AndyW]

Thank you for sharing your question about using Fast Fourier Transform to solve the Poisson equation in a N-Body code. I understand the importance of accurately simulating galaxy dynamics and appreciate your dedication to finding the best approach for generating initial velocities.

To answer your first question about why the DCT is used for computing the Fourier transform of the density field, let's first clarify the difference between the Discrete Fourier Transform (DFT) and the Discrete Cosine Transform (DCT). Both of these transforms are used to convert a signal from its original domain (e.g. time or space) to the frequency domain. The DFT is used for signals that are periodic and symmetric, while the DCT is used for signals that are periodic but not symmetric. In your case, the density field is periodic along the x and y dimensions, but not along the z dimension. Therefore, the DCT is the appropriate choice for computing the Fourier transform of the density field.

Moving on to your second question about why the author takes the interval [Ng/4, 3*Ng/4] for the x and y indices when printing the values of the potential, let's first consider the concept of symmetry in Fourier transforms. The Fourier transform of a real-valued function is complex and can be divided into a real and imaginary part. For a symmetric function, the real part of the Fourier transform is even and the imaginary part is odd. Similarly, for an anti-symmetric function, the real part is odd and the imaginary part is even. In your case, since the density field is symmetric, the real part of the Fourier transform of the density field is even, and the imaginary part is odd. Therefore, by setting the values in the 3 ghost quadrants to zero, the author is effectively making the density field anti-symmetric, which ensures that the real part of the Fourier transform is odd. This is why the author takes the interval [Ng/4, 3*Ng/4] for the x and y indices, as it ensures that the real part of the Fourier transform of the potential is also odd.

Lastly, regarding your question about why the z = 0 coordinate is obtained by taking gal->potential[j][Ng/2], let's first consider the concept of the Nyquist frequency. In Fourier analysis, the Nyquist frequency is the highest frequency that can be represented in a discrete Fourier transform. In your