Gravity Grid Parallelization

Whereas we tend to divide up our list of particles among processors to parallelize hydro routines, we are forced to divide up our gravity grids among processors spatially. For its parallel FFT routines, the FFTW requires that all grids be divided equally along the first dimension of the grid, but it is VITAL to note that the FFTW is written in C. Since C and Fortran use opposite notations for indexing vectors, C stepping along the last index and Fortran along the first, this means that in the Fortran part of the code we must break up the grids passed to the FFTW in the z-direction. Thus, all grids in the gravity routine have dimensions NNGRAVxNNGRAVxNNP, where ${\tt NNP}={\tt NNGRAV}/{\tt NPROCS}$. This is straightforward when inputting values onto the grids, but harder when outputting values back to the particles. We keep careful track of the z-coordinate of vertices whose values contribute to the potential or force of a given particle, only summing those that lie on a given processor, before using an MPI_ALLREDUCE to gather the sums properly.

The sectioning of the 3-d grids sent to the parallel 3-d FFT algorithm is computed in the subroutine ``GRAVQUANT'', which is invoked from all of the initialization subroutines. See Sec. 4 on initial conditions for more on the topic.