Mercury-Ar$χ$es: a high-performance n-body code for planet formation studies

TL;DR

Mercury-Arχes addresses the computational challenge of simulating planet formation by integrating a high-performance n-body engine (Mercury+) with a comprehensive multi-physics layer (Arχes) that models protoplanetary disks, gas drag, disk gravity, planet growth, and migration. The approach combines a hybrid symplectic integrator with OpenMP-based parallelism and vectorization, and introduces detailed HPC optimizations (pre-computation, static scheduling, and guarded critical regions) to sustain performance during large, complex simulations. Key contributions include a complete description of the HPC implementation, refined gas-disk physics (including drag, η corrections, and disk-gravity effects), two-phase growth and migration models, an updated ejection criterion, and new snapshot-based outputs for efficient post-processing. Performance assessments on Intel Performance Hybrid architectures and the Leonardo pre-exascale system demonstrate strong parallel speedups (with 85–90% parallel fraction in large runs) and identify critical cost centers (notably gas-disk interactions), establishing Mercury-Arχes as a scalable, physics-rich platform for exoplanet and Solar System formation studies. The work also outlines future GPU-port developments (Mercury-OPAL) and emphasizes compatibility with existing n-body libraries, facilitating broad adoption in the planetary science community.

Abstract

Forming planetary systems are populated by large numbers of gravitationally interacting planetary bodies, spanning from massive giant planets to small planetesimals akin to present-day asteroids and comets. All these planetary bodies are embedded in the gaseous embrace of their native protoplanetary disks, and their interactions with the disk gas play a central role in shaping their dynamical evolution and the outcomes of planet formation. These factors make realistic planet formation simulations extremely computationally demanding, which in turn means that accurately modeling the formation of planetary systems requires the use of high-performance methods. The planet formation code Mercury-Ar$χ$es was developed to address these challenges and, since its first implementation, has been used in multiple exoplanetary and Solar System studies. Mercury-Ar$χ$es is a parallel n-body code that builds on the widely used Mercury code and is capable of modeling the growth and migration of forming planets, the interactions between planetary bodies and the disk gas, as well as the evolving impact flux of planetesimals on forming planets across the different stages of their formation process. In this work we provide the up-to-date overview of its physical modeling capabilities and the first detailed description of its high-performance implementation based on the OpenMP directive-based parallelism for shared memory environments, to harness the multi-thread and vectorization features of modern processor architectures.

Figures (4)

• Figure 1: Flat profile produced by Gprof for the serial run of Mercury-Ar$\chi$es described in Sect. \ref{['sec:use_case1']}. The simulated planetary system hosted two growing and migrating planets and about 18000 massless particle interacting with the disk gas. The flat profile shows the ten most computationally expensive subroutines.
• Figure 2: Roofline analysis produced by Intel Advisor for the parallel run of Mercury-Ar$\chi$es described in Sect. \ref{['sec:use_case1']} using the same planetary system as Fig. \ref{['fig:gprof-profile']}. The red symbol is the most computationally intensive subroutine, arxes_gas, while the two yellow ones are mco_x2el and mce_cent (see discussion in Sect. \ref{['sec:use_case1']} for details). The roofline plot highlights how most of the computational load of the simulation resides in the compute and memory bound region.
• Figure 3: Comparison of the runtime values of Mercury-Ar$\chi$es in the tests of Sect. \ref{['sec:use_case1']} when running serially, in parallel with no thread binding (i.e. leaving the management of the workload schedule to the processor and the operating system) and in parallel with thread binding. All runtime values are normalized to the duration of the serial run. The theoretical dashed curve shows the expected runtime based on Amdahl's law when 90% of the workload is parallelized. The vertical dashed line marks the number of P-cores available on the processor.
• Figure 4: Comparison of the runtime of the parallel simulations of Mercury-Ar$\chi$es in Sect. \ref{['sec:use_case2']}. The simulations are run on one node of Leonardo's DCGP module using 14 threads and considering planetesimal disks extending by about 10 AU and populated by 1000, 2000, 5000 and 10000 particles/au. The green curve with circle symbols shows the evolution of the wallclock runtime, the golden curve with square symbols shows the average runtime per particle. The parallel simulations are run with the version of Mercury-Ar$\chi$es containerized with Singularity adopted by the OPAL project (Polychroni et al., this issue). Mercury-Ar$\chi$es is compiled with Intel OneAPI using the compilation flags described in Sect. \ref{['sec:use_case2']}.