FFT Convolution

There are SPH techniques which can be used to solve for the gravitational force on each particle, but we have found it easiest to use FFT convolution to calculate the gravitational potential and force for each particle. Note that according to Newton's law, the gravitational potential $\Phi$ satisfies an equation

$$\Phi(\vec{r})= \int \rho({\vec{r}}') \frac{d^3 {\vec{r}}'}{\... ...r}}'\vert} = \sum_i \frac{m_i}{\vert\vec{r}-\vec{r}_i\vert} .$$ (10)

The direct sum requires $N^2$ operations, and is unfeasible to implement. The integral, on the other hand, is in the form of a convolution, and we find

$$\Phi(r)={\rm FFT}^{-1}\times \left[{\rm FFT}(\rho)*{\rm FFT}\left(\frac{1}{r}\right)\right].$$ (11)

To use FFT convolution techniques, we have to import the data to a 3-dimensional grid. Setting the boundaries for this grid is a topic worthy of some explanation, and can have significant effects on both the performance and accuracy of the code. Essentially, if the grid covers too large a spatial volume, the matter becomes concentrated in fewer cells, and we lose the sharpness of our solution for the gravitational force and potential. If the grid is too small, though, particles can be lost off the edges, which in some cases can be a problem. Resizing the grid allows for flexibility, but requires doing a third FFT (that of $1/r$) every timestep in addition to that of the matter and the inverse FFT of the potential, adding up to $50\%$ to the time required to compute the gravitational field. If the grid spacing is known to be constant, we can calculate the FFT of $1/r$, save the result, and reuse it over and over again. In all cases, the dimensions of the grid are fixed in advance, by the parameter ``NGRAV'' in the header file, representing the size of the grid devoted to the matter. The actual size of the grid is ${\tt NNGRAV}^3\equiv (2 \times{\tt NGRAV})^3$, since convolution techniques require double the size in every dimension to properly calculate the potential, a process known as ``zero-padding''. If you will, the convolution has to account for the fact that one particle can be located across the grid all the way to the right, or all the way to the left, so the range of values of $\vec{x}_i-\vec{x}_j$ spans a length twice that of the size of the grid. It will come as no surprise that the FFT works dramatically better if NGRAV is chosen to be a power of two, or at least contains nothing but small prime factors. The FFTW can handle multiples of 5 as well as other small primes, but at a huge loss of computational speed.

When calculating the density at each grid vertex, we use a ``cloud-in-cell'' method, whereby the mass contributions from any given particle are attributed to the eight vertices surrounding by means of a weighting function which contains the product of three terms, each of which are linear in a dimension, with the total weighting normalized to unity. If a particle is respectively $25\%, 40\%$, and $70\%$ of the way across the cell from the lower front left corner in the three dimensions, we attribute $0.75\times 0.6\times 0.3$ of the particle's mass to that corner.

After calculating the gravitational potential at each vertex, we finite difference to find the gravitational force at each vertex. The potential and gravitational force are interpolated for each particle using the ``cloud-in-cell'' technique described above, with the same weightings used on the way out as were used on the way in.

Particles not on the grid are flagged, and a warning will be printed saying ``QGRAV: NPOUT=... AMIN='', where NPOUT is the number of particles off the grid, and AMIN is the total mass of particles located within the grid. For particles off the grid, gravitational terms are calculated between them and particles on the grid (but not among particles off the grid) via a simple monopole approximation.