  This paper presents an accurate density computation approach for large dark
matter simulations, based on a recently introduced phase-space tessellation
technique and designed for massively parallel, heterogeneous cluster
architectures. We discuss a memory efficient construction of an oct-tree
structure to sample the mass densities with locally adaptive resolution,
according to the features of the underlying tetrahedral tessellation. We
propose an efficient GPU implementation for the computationally intensive
operation of intersecting the tetrahedra with the cubical cells of the deposit
grid, that achieves a speedup of almost an order of magnitude compared to an
optimized CPU version. We discuss two dynamic load balancing schemes - the
first exchanges particle data between cluster nodes and deposits all tetrahedra
for each block of the grid structure on single nodes, whereas the second
approach uses global reduction operations to obtain the total masses. We
demonstrate the scalability of our algorithms for up to 256 GPUs and TB-sized
simulation snapshots, resulting in tessellations with over 400 billion
tetrahedra.
