Speeding up cosmological hydrodynamic simulations by exploiting space-filling curve

Space filling curves (see Bader et al. @ Encyclopedia of Parallel Computing, for a cohomprensive review) are widely used in simulations for load balancing as they approximate the ideal domain decomposition (Dubinski et al. 2009). As a result, nearby particle in memory are also nearby particle in space. Can we exploit it to speedup the neighbour search? Yes.

During my PhD, our numerical group was facing some slow-downs in running large cosmological hydrodynamic simulations with Gadget3 (at the time, private repo).

We used scalasca and profiled the hydrodynamic kernels runtimes of Magneticum simulations. We found that most of the time was not even spent actual physics computations, rather, most of the time was spent in the tree walk neighbour search!

The Insight

Gadget3 orders its particles using a Hilbert space-filling curve for domain decomposition. The Hilbert curve has a key property: points nearby along the curve are also nearby in 3D space. So particles stored consecutively in memory are physically close, and physically close particles share most of their neighbours.

The fix, which I called Neighbour Recycling: instead of one tree walk per particle, group nearby particles together and perform a single tree walk for the whole group, over a slightly larger radius.

To handle Gadget3's wildly varying densities, the grouping radius is made adaptive: R = f * hsml, (where hsml is the SPH smoothing length), and the number of particles in a group is N_group = R^3.

How to find the optimal value for f parameter?

If f is too small, we have no groupings and we revert to legacy-Gadget3, if f is too large, we will create too large neighbour groups and we will take too much time in filtering them.. so there must be a compromise where the code is fastest: for which value of f?

In a ideliazied homogeneus setup - the number of neigbour of a particle scales as N_ngb ~ hsml^3, the superset of neighbour candidets as N_candidates ~ (R + hsml)^3. There it is easy to derive that f = (N_candidates/N_ngb)^1/3 - 1 = (N_group/N_ngb)^1/3.

To find f we need to compute the time spent in neighbour search and minimize it. In legacy Gadget for instance, it would be T = N_active * N_ngb, where N_active is the number of active particles in this timestep and N_ngb is the number of neighbours (tipically N_ngb = 295), given that tree walk neighbour search scales with the number of neighbours (see wikipedia).

When it comes to the neighbour reciclying the number of searches goes with the number of groups N_groups = N_active/N_group and then an additional term needed to filter out the exta neighbours candidates for each active particle (note this is a very crude and naive estimate, just to get an approximate evaluation). So we obtain:

          T = N_groups * N_ngb  + #<- tree walks
            + N_active * N_candidates #<-filtering
        

In the paper Ragagnin et al. (2016) we find the optimal constant f = 0.5 was found numerically, and I think it is astonishing how the above analytical computation provided a grouping radius that is now since long used in Gadget - succesfully - as optimal value.

Results

Final remarks

The technique doesn't rely on anything Gadget3-specific beyond the space-filling curve ordering, which most modern N-body codes already use. If your code spends more time finding neighbours than computing physics, this idea is worth a look.

This technique can be activated in the SVN Gadget3 code by using the switches AR_GREENTREE_{DENSITY,HYDRA,CONDUCTION}. And is activated by default in OpenGadget3 (see file CodeBase/greentree.hpp).

Bibtex entry:

@INPROCEEDINGS{Ragagnin2016GreenTree,
          author={{Ragagnin}, Antonio and {Tchipev}, Nikola and {Bader}, Michael and {Dolag}, Klaus and {Hammer}, Nicolay J.},
          title="{Exploiting the Space Filling Curve Ordering of Particles in the Neighbour Search of Gadget3}",
          keywords = {Astrophysics - Instrumentation and Methods for Astrophysics, Computer Science - Performance},
          booktitle = {Advances in Parallel Computing},
          year = 2016,
          month = may,
          volume={41},
          pages = {411-420},
          doi = {10.3233/978-1-61499-621-7-411},
          archivePrefix = {arXiv},
          eprint = {1810.09898},
          primaryClass = {astro-ph.IM},
          adsurl = {https://ui.adsabs.harvard.edu/abs/2016pcre.conf..411R},
          adsnote = {Provided by the SAO/NASA Astrophysics Data System}
        }