• 1.Introduction
• 2.Solution of 3D Poisson equation using Green's function
• 3.Discretization : Green's function
• 4.Discretization : Matter Distribution - "Particle Mesh" method
• 5.Implementation using FFT/GPU/OpenCL
• 6.Final solution
• 7.Validation

1.Introduction :

We are interested here in solving numerically the 3D Poisson equation in the gravitational domain. The aim is to compute the gravitational potential of a disk galaxy from an initial distribution of particles positions representing stars; This computation is achieved into BmyOGG simulation code in order to generate initial velocities of these stars.

The method used is based on FFT (Fast Fourier Transform) as they are implemented in GPU/OpenCL API. One starts by express Poisson equation into Fourier space using Green's function, in continuous and discret case. Then one has to best describe the distribution of matter on the chosen mesh : this is performed with PM (Particle Mesh) method. Finally, we apply FFT synthesis to get the gravitational potential of the galaxy.

Source code of this project is in the following archive :

• Solving_Poisson_FFT_GPU.tar.gz

  (tested with AMD APP SDK 3.0 and clFFT 2.8)

2.Solution of 3D Poisson equation using Green's function :

3D Poisson equation in gravitatonal domain is expressed as :

where $\Phi$ is the gravitational potential and $\rho$ the matter density. Green's function is defined as the solution of :

or equivalently in cartesian coordinates :

We can demonsrate that Green's function, solution of above equation, is written :

Then, we get the following solution of Poisson equation :

where $*$ is the convolution product. If we apply a forward Fourier Transform (FT) to \eqref{sol1}, we get :

For continuous case, we can easily calculate $\hat{G}(\vec{k})$ :

So for $\hat{\Phi}(\vec{k})$ :

Finally, we need to perform a backward FT to get $\Phi$ in spatial domain, i.e $\Phi(\vec{r})$.

We perform the numerical solving of Poisson equation so that we need to work in discret case (replacing FT by FFT). So, our code has the following steps :

• Read a data file of initial positions for particles (here stars in a disk galaxy)
• Define a cartesian mesh $(x,y,z)$
• Discretization of the Green's function $G$ on the previously defined mesh
• Discretization of matter distribution $\rho$ on the mesh with interpolation "Particle-Mesh" method
• Apply a forward FFT on Green's function $G$
• Apply a forward FFT on Matter distribution $\rho$
• Apply a backward FFT on the product of the two previous FFT to get potential

3.Discretization : Green's function

Green's function in spatial domain is defined by \eqref{green-start}. We notice there is a singularity for $|\vec{r}|=0$. This can be circumvented by setting $G(x=0,y=0,z=0)=1$. Below the Green's function (in absolute value) represented for $z=0$ on a cartesian mesh with size $N=256$. We have indices equals to [i][j][k] = [-128,127][-128,127][-128,127] :

4.Discretization : Matter Distribution - "Particle Mesh" method

In order to compute matter distribution, we use a data set of $(x,y,z)$ stars coordinates in spiral galaxy. These data have been generated thanks to Rejection-Sampling method used in BmyOGG code. We represent the matter distribution of this spiral galaxy (10240 stars) on the following WebGL canvas; you can zoom with mouse wheel or rotate the scene by dragging. To check if your browser supports WebGL or if you want to activate it, go to Enabling WebGL.

"Particle Mesh" method :

"Particle Mesh" is a numerical method for computing forces in particles system. Particles being atoms, stars or components of a fluid, this method is used in different fields like plasma physics or astrophysics. The main idea is that particles system can be converted into a mesh of density values. Then, the potential is solved on this mesh.

There are a lot of methods to convert the particles system into grid. The simpliest is to assign mass particle to its closest point on the grid : this is the "Nearest-Grid-Point" (NGP) method. Another is to model each particle like a cube of constant value and the mass is decomposed along the neighboring cells : it is called the "Cloud-In-Cell" method (CIC). We have implemented the 2 methods into code but the best interpolation is done by CIC method.

Below a schema representing in 2D case the distribution of particle mass located at $(x,y)$ point on neighboring cells (CIC method).

Now, we are going to see how mass repartition is discretized, taking :

case 1 - If $\Delta x > \dfrac{H_{x}}{2}$ and $\Delta y > \dfrac{H_{y}}{2}$ :

case 2 - If $\Delta x < \dfrac{H_{x}}{2}$ and $\Delta y > \dfrac{H_{y}}{2}$ :

case 3 - If $\Delta x > \dfrac{H_{x}}{2}$ and $\Delta y < \dfrac{H_{y}}{2}$ :

case 4 - If $\Delta x < \dfrac{H_{x}}{2}$ and $\Delta y < \dfrac{H_{y}}{2}$ :

We have also to handle cases where we are on edges of the mesh : for example, if we have $i=1$ and $\Delta x<\dfrac{H_{x}}{2}$ or $i=N$ and $\Delta x>\dfrac{H_{x}}{2}$, $N$ being the number of points on X axis. In these special cases, we must ignore external mesh cells.

Results :

We show here results with a 256x256x256 mesh size. Density function is computed with data set Positions_Disk.dat and the numerical interpolation method is plotted on the following figure :

Density values on this figure are not significant since a factor has been added to remain in range of float type. This factor is deleted in the final calculation of potential.

5.Implementation using FFT/GPU/OpenCL :

FFT computation is carried out using clMathLibraries/clFFT. Maximum size allowed for float type (single precision) is 2^24. For double type (double precision), it's 2^22. We made the choice in our code to work with float type so that 3D arrays have 256x256x256 dimensions.

After reading the input file (containing $(x,y,z)$ coordinates of stars), we produce a cartesian mesh that we automatically adapt to the max,min of positions. All coordinates are in kiloparsec unit (1kpc = 3,086e+21 cm). Then, we have to numerically implement Green's function as it is indicated in (3) part. We do the same to compute density function. Finally, we do forward FFT synthesis for each of these two 3D arrays, make their product and apply backward FFT to get the gravitational potential on 3D mesh. All this computing chain is done in CGS units : so gravitational potential will be expressed in $cm^{2}.s^{-2}$ units.

Concerning the Fourier transforms, we must not forget to include right factors when we compute a backward FT on the product of 2 forward FT. Indeed, each forward FT makes appear a $(space_{x}\,space_{y}\,space_{z})$ factor and backward FT includes a $1/(N^{3}\,space_{x}\,space_{y}\,space_{z})$ factor. So, we have to multiply final solution by $(space_{x}\,space_{y}\,space_{z})$.

6.Final solution :

Since it is difficult to visualize potential values on a 3D space (volume rendering), we choose to do a "slice plot", i.e we show values in a 3D space cut by (Oxz) and (Oyz) plans with a colormap scale.

Below the final solution on a 256x256x256 mesh plotted with the Matlab script plot_potential_slice.m :

7.Validation :

Validation of final solution is performed in 2 ways; the first one is done from the definition of gravitational potential :

We call this way the "direct computation" method. The second way is to use multi-dimensional Matlab/FFT. The script compute_potential_Matlab.m implements these 2 options.

Here the direct computation on a 256x256x256 mesh :

As with C language version, we start by loading input positions, then discretize Green's and density function. Their FFT are applied with "fftn" Matlab function for 3D case. On this figure the potential got with Matlab :

Finally, the results of these 2 validation tests are similar to those obtained with FFT/GPU. Let's compare the average value and the standard deviation between FFT/GPU, direct computation and FFT/Matlab methods :

We see that the values remain in the same range for the 3 methods. While there may be differences locally due to numerical artifacts, the method we used in our code is validated.