MPI parallel mesh gravity - hashmap free

Implements the mesh gravity calculation in an MPI-distributed fashion. This allows to run with much larger meshes and reduce the memory footprint.

Base implementation is similar to !1045 (closed).

• Each rank accumulates it's own contributions to the density field. These would need to be stored in some kind of sparse array to avoid making any assumptions about how the domain decomposition is done. Unlike !1045 (closed), we do not use a hashmap but rather a simple array since we can actually predict what we need and construct nice keys.
• Each rank allocates a slab of the full mesh as required by FFTW
• The mesh contributions are sent to whichever rank stores the corresponding slab. This is reasonably straightforward because the coordinates of a cell indicate which rank it's stored on.
• The FFT/Green function/CIC deconvolution/inverse FFT are carried out to make a slab-distributed mesh containing the potential
• Each rank calculates which cells it needs to calculate the potential gradient for its particles and requests them from whichever node they're on. This is slightly trickier than constructing the mesh because we need several cells around each particle to evaluate the gradient.
• We can then evaluate the potential and acceleration on the particles as usual.

Another difference with !1045 (closed) is that I am using a hand-written bucket sort instead of the qsort() that was used originally. Since we only need a crude sort with a fixed small number of bins this is much much faster. (But still the current bottleneck)

The code needs to be configures with --enable-mpi-mesh-gravity and the runtime parameter Gravity:distributed_mesh must be set to 1.

Implements #524 (closed). Bypasses #716 (closed). Likely supersedes !1045 (closed).

Possible improvements:

• Use the threadpool to speed-up the three bucket sorts. At least to construct the counts.
• Construct the array of send/recv counts at the same time as the bucket sort. Or just from the bucket counts?
• Check whether the potential patches can be smaller than currently and only cover the particle extent.

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -292,6 +292,7 @@ Simulations using periodic boundary conditions use additional parameters for the
 Particle-Mesh part of the calculation. The last five are optional:
 
 * The number cells along each axis of the mesh :math:`N`: ``mesh_side_length``,
+* Whether or not to use a distributed mesh when running over MPI: ``distributed_mesh`` (default: ``0``),
 * The mesh smoothing scale in units of the mesh cell-size :math:`a_{\rm
   smooth}`: ``a_smooth`` (default: ``1.25``),
 * The scale above which the short-range forces are assumed to be 0 (in units of
@@ -305,6 +306,18 @@ For most runs, the default values can be used. Only the number of cells along
 each axis needs to be specified. The remaining three values are best described
 in the context of the full set of equations in the theory documents.
 
+By default, SWIFT will replicate the mesh on each MPI rank. This means that a
+single MPI reduction is used to ensure all ranks have a full copy of the density
+field. Each node then solves for the potential in Fourier space independently of
+the others. This is a fast option for small meshes. This technique is limited to
+mesh with sizes :math:`N<1291` due to the limitations of MPI. Larger meshes need
+to use the distributed version of the algorithm. The code then also needs to be
+compiled with ``--enable-mpi-mesh-gravity``. That algorithm is slower for small
+meshes but has no limits on the size of the mesh and truly huge Fourier
+transforms can be performed without any problems. The only limitation is the
+amount of memory on each node. The algorithm will use ``N^3 * 8 * 2 / M`` bytes
+on each of the ``M`` MPI ranks.
+
 As a summary, here are the values used for the EAGLE :math:`100^3~{\rm Mpc}^3`
 simulation:
 
@@ -316,7 +329,8 @@ simulation:
      MAC:                    adaptive
      theta_cr:               0.6
      epsilon_fmm:            0.001
-     mesh_side_length:       512
+     mesh_side_length:       2048
+     distributed_mesh:       0
      comoving_DM_softening:         0.0026994  # 0.7 proper kpc at z=2.8.
      max_physical_DM_softening:     0.0007     # 0.7 proper kpc
      comoving_baryon_softening:     0.0026994  # 0.7 proper kpc at z=2.8.