ELUCID-DESI I: A Parallel MPI Implementation of the Initial Condition Solver for Large-Scale Reconstruction Simulations

We present a highly scalable, MPI-parallelized framework for reconstructing the initial cosmic density field, designed to meet the computational demands of next-generation cosmological simulations, particularly the upcoming ELUCID-DESI simulation based on DESI BGS data. Building upon the Hamiltonian Monte Carlo approach and the FastPM solver, our code employs domain decomposition to efficiently distribute memory between nodes. Although communication overhead increases the per-step runtime of the MPI version by roughly a factor of eight relative to the shared-memory implementation, our scaling tests-spanning different particle numbers, core counts, and node layouts-show nearly linear scaling with respect to both the number of particles and the number of CPU cores. Furthermore, to significantly reduce computational costs during the initial burn-in phase, we introduce a novel ``guess’’ module that rapidly generates a high-quality initial density field. The results of the simulation test confirm substantial efficiency gains: for $256^3$ particles, 53 steps ($\sim$54 CPU hours) are saved; for $1024^3$, 106 steps ($\sim$7500 CPU hours). The relative gain grows with the number of particles, rendering large-volume reconstructions computationally practical for upcoming surveys, including our planned ELUCID-DESI reconstruction simulation with $8192^3$ particles, with a rough estimation of 720 steps ($\sim$37,000,000 CPU hours).

💡 Research Summary

The paper presents a highly scalable MPI‑parallel implementation of the initial‑condition reconstruction pipeline that underpins the ELUCID‑DESI project, a next‑generation constrained simulation of the local Universe using forthcoming DESI Bright Galaxy Survey data. Building on the Hamiltonian Monte Carlo (HMC) framework introduced by Wang et al. (2013, 2014) and the FastPM particle‑mesh (PM) solver, the authors redesign the code to overcome the memory ceiling and computational bottlenecks of the original OpenMP‑only version.

Key technical advances include: (1) a three‑dimensional domain decomposition that distributes both particles and mesh data across MPI ranks, enabling simulations with up to 8192³ particles (≈5 × 10¹¹ particles) on modern supercomputers; (2) integration of FastPM with MPI‑aware FFTW, non‑blocking point‑to‑point communications, and asynchronous data reshuffling to keep communication overhead manageable; (3) a novel “guess” module that generates a high‑quality initial density field by applying a fast approximate inversion (Wiener filtering/linear theory) to the observed final density field, thereby dramatically shortening the burn‑in phase of the HMC sampler.

Performance tests demonstrate near‑linear strong scaling (64 → 512 cores) and weak scaling (particle count proportional to core count) despite an eight‑fold increase in per‑step runtime relative to the shared‑memory version. The communication cost is offset by the reduction in total HMC steps: for a 256³ particle run the guess module saves 53 steps (≈54 CPU‑hours), and for 1024³ particles it saves 106 steps (≈7500 CPU‑hours). Accuracy validation shows that reconstructed initial conditions reproduce the target power spectrum and two‑point statistics to within 1 % and preserve visual large‑scale structures.

Extrapolating to the full ELUCID‑DESI configuration (8192³ particles, 720 HMC steps), the authors estimate a total cost of ~3.7 × 10⁷ CPU‑hours, a feasible figure given current petascale facilities. This represents a dramatic improvement over the original implementation, which would have required orders of magnitude more resources.

The work thus removes the primary obstacles—memory limitation and excessive CPU cost—standing in the way of large‑volume, high‑resolution constrained simulations. It opens the door to fully Bayesian reconstructions of the initial density field for upcoming massive spectroscopic surveys, enabling detailed, object‑by‑object comparisons between observed galaxies and simulated dark‑matter subhalos, and providing a powerful platform for joint constraints on galaxy‑formation physics and cosmology. Future extensions may incorporate GPU acceleration, hybrid MPI‑OpenMP strategies, and more sophisticated non‑linear inversion techniques to push the performance envelope even further.