Assessing Performance and Porting Strategies for Gravitational N-Body Simulations on the RISC-V-Based Tenstorrent Wormhole™ Jenny Lynn Almerola,∗, Elisabetta Boellaa,∗∗, Mario Sperab , Daniele Gregoria a E4 Computer Engineering SpA, Scandiano, 42019, Italy
arXiv:2605.02744v1 [cs.DC] 4 May 2026
b Scuola Internazionale Superiore di Studi Avanzati (SISSA), Trieste, 34136, Italy
Abstract While RISC-V-based accelerators were initially designed with artificial intelligence applications in mind, they are increasingly being recognized as promising platforms for high performance scientific computing. In this work, we present three strategies for scaling an N-body code across multiple Tenstorrent Wormhole accelerators based on the RISC-V architecture. We assess the performance of these approaches by measuring both the execution time and the energy consumption required to complete a representative simulation, ultimately identifying the configuration that offers the most favorable balance between efficiency and performance. Keywords: RISC-V accelerator, Tenstorrent Wormhole, N-body simulations, Computational astrophysics, Energy efficiency
1. Introduction Interest in the RISC-V Instruction Set Architecture (ISA) has grown rapidly across the High Performance Computing (HPC) community, driven by its open, royalty-free standard and its inherently modular and extensible framework (Patterson and Waterman, 2017). Though conceived for embedded systems, RISC-V has matured through the addition of 64-bit cores, vector processing units, and domain-specific accelerators, making it increasingly suitable for high-end computational tasks (Venieri et al., 2025; Mahale et al., 2023; Bartolini et al., 2022) This evolution is further fueled by the rise of artificial intelligence (AI) and machine learning (ML), whose computational demands strongly overlap with those of traditional scientific workloads (Brown and Barton, 2024; Brown et al., 2025). In scientific computing, the increasing demand for performance and energy efficiency is driving the exploration of new architectures. Astronomy and astrophysics, in particular, stand at the forefront of data-intensive science, with next-generation surveys and simulations requiring unprecedented computational throughput (Shukla et al., 2025b; Suarez et al., 2025; Shukla et al., 2025a). To meet these requirements, researchers are increasingly adopting heterogeneous and ∗ [email protected] ∗∗ [email protected]
Preprint submitted to Astronomy & Computing
May 5, 2026
accelerator-based systems that can deliver high performance per Watt while scaling to massive problem sizes. Among the most computationally demanding workloads are N-body simulations, which model the dynamical evolution of astrophysical systems. Their reliance on large-scale linear algebra and particle–particle interactions makes them particularly well suited to hardware acceleration (Khan et al., 2021; Rexroth et al., 2020), and therefore ideal candidates for evaluating the potential of novel RISC-V–based accelerators. A noteworthy example of a RISC-V-based accelerator targeting the HPC domain is the Tenstorrent Wormhole n300 (Tenstorrent, 2025g). Designed to decouple computation from data movement, it offers high parallel efficiency and energy performance at a comparatively modest cost (Brown et al., 2025). These characteristics render it a compelling testbed for exploring the potential of RISC-V accelerators in astrophysical contexts. However, adopting new architectures for scientific computing presents a significant challenge, primarily due to the extensive and complex codebases developed for large-scale scientific numerical experiments, which often comprise millions of lines of code. Porting such codes is justified only when substantial performance or energy efficiency improvements can be achieved. Motivated by this, we developed an N-body simulation code in C++ and offloaded the most computationally intensive calculations using the Tenstorrent TT-Metalium programming interface (Tenstorrent, 2025e). To the best of our knowledge, this represents the first example of an N-body application ported to Tenstorrent accelerators. Thus, our application serves as a concrete test case for assessing the capability of RISC-V in demanding astrophysical workloads. The accelerated code was evaluated on the Tenstorrent Wormhole n300 card. Performance measurements show that it achieves over a 2× speedup and roughly 2× higher energy efficiency compared to a highly optimized baseline version executed on x86 64 (Almerol et al., 2025). In this work, we enhance the code parallelism by leveraging the Message Passing Interface (MPI) paradigm to use multiple accelerators simultaneously. In addition to this distributed memory approach, we implement strategies to exploit card level parallelism, enabling the utilisation of both chips present on the n300 accelerator. As far as we are aware, this is the first scientific code capable of leveraging multiple Tenstorrent cards and multiple chips for large-scale computations. Thus, the main contributions of this paper are: • We explore multiple approaches to enhance the parallelization of our N-body code on the RISC-V-based Tenstorrent accelerators; • We evaluate the performance, scalability, and energy efficiency of the different proposed strategies. This paper is structured as follows. In Section 2, we introduce the direct N-body simulation model and the Tenstorrent Wormhole accelerator. Section 3 outlines the porting and optimization strategy. In Section 4 validation and benchmarking procedures are presented and performance results, including single- and multi-chip scalability as well as energy efficiency analysis, are discussed. Finally, Section 5 illustrates concluding remarks and future perspectives. 2. Background 2.1. Direct N-Body Simulation Direct N-body simulations numerically solve the equations of motion for a system of particles by directly summing the gravitational forces exerted by all other particles in the system, 2
expressed as: Fi =
N X
G
j=1 j,i
mi m j ri3j
(r j − ri ),
(1)
where mi and m j are the particle masses, ri and r j their position vectors, ri j = |r j − ri | is the inter-particle distance, and G the gravitational constant. Direct summation is the most straightforward approach to the N-body problem and also the most accurate, as it does not rely on approximations (Spurzem, 1999). However, this accuracy comes at a steep computational cost, since evaluating all pairwise interactions scales as O(N 2 ). Therefore, developing accurate, efficient, and scalable direct N-body codes is essential, particularly in advancing the astrophysical interpretation of gravitational-wave observations from current detectors such as LIGO, Virgo, and KAGRA (Abbott et al., 2020), as well as future observatories such as the Einstein Telescope (Einstein Telescope Italy, 2025). In this work, we employ a high-order time integration scheme, the sixth-order Hermite integrator (Nitadori and Makino, 2008), and port its implementation to the RISC-V–based Tenstorrent Wormhole accelerator. The Hermite scheme comprises three iterative stages: prediction, evaluation, and correction (Spera, 2014). In the prediction step, the positions, velocities, and accelerations of all particles are estimated from their previously known values. In the evaluation step, accelerations and their first time derivatives (jerks) are computed using these predicted quantities. Finally, in the correction step, the predicted positions and velocities are refined using the newly evaluated accelerations and jerks, achieving sixth-order time integration accuracy. In this work, we optimize the computationally intensive all-to-all interaction of particles which involves calculating particle accelerations and jerks to exploit the parallel capabilities of the Wormhole accelerator effectively. 2.2. The Tenstorrent Wormhole™ Accelerator The Wormhole is an Application-Specific Integrated Circuit (ASIC) consisting of 64 programmable processing elements, known as Tensix cores (Tenstorrent, 2025f,g). Each Tensix core contains five “Baby” RISC-V CPUs, 1.5 MB of local SRAM (L1), two network routers interfacing with Networks-on-Chip (NoCs), a math unit (Floating Point Unit, FPU) for low-precision matrix arithmetic supporting up to bfloat16 (BF16) operations, and a wide SIMD vector engine (Scalar Floating Point Unit, SFPU)supporting up to single-precision (FP32) operations (Tenstorrent, 2025f; Corsix, 2024), as illustrated in Fig. 1. The “Baby” RISC-V processors are 32-bit, single-issue cores operating at 1 GHz (Corsix, 2024). Two of these cores (RISC-V “NC” and “B”) are dedicated to data movement, managing transfers between the Tensix core and off-chip DRAM. The remaining three are compute cores that coordinate arithmetic and logic operations by issuing instructions to the various coprocessor units, including the tensor FPU and SFPU. Each compute core drives a specific stage of the compute pipeline (Tenstorrent, 2025f). For instance, RISC-V T0 (UNPACK) issues instructions to the unpacker module to load data from SRAM into two source registers, srcA and srcB for the matrix engine or the dst register for the vector engine. RISC-V T1 (MATH) then issues instructions to the SFPU and FPU to perform operations on these registers. Finally, RISC-V T2 (PACK) coordinates the movement of results from the destination register, dst, back to SRAM (Tenstorrent, 2025e). This data flow is illustrated in Fig. 1. 3
Figure 1: Simplified schematic overview of a Tensix core in the Tenstorrent Wormhole AI accelerator. Blue arrows correspond to instruction dispatch, while black arrows indicate data movement. Adapted from Chang (2025).
This hardware partitioning is mirrored in the software execution model, where we implement three custom kernels, read, compute, and write, each corresponding to a distinct stage of the processing pipeline. These kernels are executed across data movement and compute cores in a dataflow-driven manner, communicating via software-managed circular buffers (CBs). Kernel interactions follow a producer–consumer paradigm, with data transfers implemented as First-In, First-Out (FIFO) streams to ensure ordered and efficient communication between stages (Tenstorrent, 2025a). Synchronization between kernels is managed via TT-Metalium buffer control primitives: cb wait front, cb pop front, and cb reserve back (Tenstorrent, 2025b,e). The first two functions manage consumer-side synchronization, ensuring that a kernel proceeds only when the required input data are available and that data are consumed in the correct order. The cb reserve back primitive manages producer-side synchronization by suspending write operations until adequate buffer space becomes available, thereby preventing data overwrites and enforcing back-pressure. Once space is reserved, the producer writes data and invokes cb push back to finalize the insertion. 3. Exploiting Parallelism on Tenstorrent Wormhole 3.1. Porting the N-Body Simulation to the Tenstorrent Wormhole Building on the execution model described in Section 2.1, we port the direct N-body simulation to the Tenstorrent Wormhole architecture using three custom kernels: read, compute, and write, each responsible for a distinct stage of the computation pipeline. Data flow across these kernels is organized in a tiled layout. A tile is a 2D grid of values the architecture operates on, has a size of 32 × 32, and can thus contain up to 1024 elements of the same data type. Since the Tenstorrent Wormhole natively supports up to FP32, a mixed-precision strategy is adopted for our N-body code. Calculations performed during the evaluation phase, including the acceleration and jerk computations on the Wormhole accelerator, are executed with FP32 precision, whereas all remaining computations on the CPU cores are performed with FP64. While in our case, this configuration is dictated by the hardware capabilities of the accelerator, previous studies have consistently shown that mixed-precision schemes deliver substantial performance and 4
Algorithm 1 Asynchronous Data Loading for Particle Tiles 1: for each source tile i do 2: Reserve buffer space for positions (r x , ry , rz ) and velocities (v x , vy , vz ) 3: Load positions and velocities of the i-th particle from memory 4: Wait for data to become available 5: Push loaded data to consumer buffers 6: for each target tile j do 7: Reserve buffer space for positions (r x , ry , rz ), velocities (v x , vy , vz ), and weights (pw ) 8: Load positions, velocities, and weights of the j-th particle from memory 9: Wait for data to become available 10: Push loaded data to consumer buffers 11: end for 12: end for
energy efficiency gains without compromising numerical accuracy (Micikevicius et al., 2018; Schäfer et al., 2020; Rexroth et al., 2020). To manage data movement across the memory hierarchy, we allocate buffers in off-chip DRAM and asynchronously copy data from the CPU to the device using non-blocking EnqueueWriteBuffer operations. The kernels then consume this data directly from DRAM. In addition, we allocate on-chip SRAM circular buffers to stage data between the read, compute, and write kernels. The read kernel loads particle data from off-chip DRAM, which corresponds to one of the “Baby” RISC-V cores issuing instruction to load data to the SRAM. The compute kernel waits for data availability from the reader and then performs the core force calculations (see Appendix A). Within this kernel, the copy tile operation loads input tiles into the SFPU input register, and element-wise tile operations are executed using compute API functions such as sub binary tile(), square tile(), rsqrt tile(), and related primitives (Tenstorrent, 2025b). The results are then stored in output circular buffers via pack tile, enabling the write kernel to stream the computed results to the output buffers (Tenstorrent, 2025b,e). Finally, the output data is copied back to the host using EnqueueReadBuffer operations. To further optimize execution, we also implement custom ternary SFPU functions tailored to the computational patterns of the force calculation kernel, including the squared-distance calculation between particles and the mul-add operation needed for acceleration and jerk calculation. Each Tensix core destination (dst) register can accommodate up to 8 tiles of FP32 data format. Given this limited register space, frequently reused intermediate values within the force computation, such as the displacement vector components (r x,i j , ry,i j , rz,i j ), are staged in on-chip SRAM using circular buffers (CBs). This staging strategy is essential for efficient data reuse and minimizing off-chip memory traffic. Consequently, the compute kernel performs multiple acquire–release operations on the dst registers to accommodate intermediate computations within the available register space. These acquire–release operations are necessary because the kernel must reuse more intermediate values than the 8-tile dst registers can hold, temporarily staging data in on-chip SRAM to avoid overwriting values needed for subsequent particle interactions. To enable parallel execution, we employ a Single Program, Multiple Data (SPMD) model. Each Tensix core executes the same program but operates on a distinct subset of the dataset to compute particle accelerations and jerks. Workload distribution is managed using Metalium split work to cores function (Tenstorrent, 2025e). However, to compute the net gravitational 5
Algorithm 2 Asynchronous N-Body Acceleration and Jerk Computation and Accumulation 1: for each source tile i do 2: Wait for input positions (r x , ry , rz ) and velocities (v x , vy , vz ) of the i-th particle 3: Reserve buffer space for acceleration and jerk 4: Initialize local acceleration (atx , aty , atz ) and jerk (ȧtx , ȧty , ȧtz ) tiles to zero 5: for each target tile j do 6: Wait for input positions (r x , ry , rz ), velocities (v x , vy , vz ) and weights (pw ) of the j-th particle 7: Compute and update pairwise accelerations and jerks (see Appendix A) 8: end for 9: Pack accumulated (atx , aty , atz ) and (ȧtx , ȧty , ȧtz ) into output buffers 10: Push finalized acceleration and jerk to consumer buffers 11: end for
force on a particle, each core must access the complete set of particle data. To accommodate this and to match the architecture tile-based data-access pattern, particle data are replicated and organized into tiles of 1024 elements, with each tile containing the same information for a single particle attribute (e.g., position or velocity component). In other words, each tile contains 1024 identical copies of one scalar quantity for a given particle. This replication strategy scales with adding more computational resources such as using multiple Wormhole cards and chips, since each chip stores its own copy of the replicated data in local DRAM. The data loading process implemented in the read kernel is illustrated in Algorithm 1, while the data-parallel organization is illustrated in Fig. 2: column tiles are distributed across Tensix cores, while rows correspond to parallel computations. Acceleration and jerk values are accumulated along the indicated directions and, once complete, written to the output CBs. Fig. 2 also shows the pipelined execution of read, compute, and write kernels, allowing data movement to overlap with computation. The data consumption and computation process implemented in Algorithm 2 follows the same asynchronous streaming pattern as the data production stage shown in Algorithm 1. An alternative approach to parallel execution on Wormhole devices is the mesh programming model (Tenstorrent, 2025c). In this model, a MeshDevice acts as a virtual abstraction that unifies multiple physical devices into a single logical mesh, in principle allowing efficient distribution of computations across a multi-device topology. When operations are issued to a MeshDevice, command queues automatically distribute workloads to all participating devices concurrently, ideally reducing dispatch overhead compared to sequential, device-by-device execution. To scale the simulation across multiple devices, we utilize up to four Wormhole n300 cards in isolated configuration (e.g. without any internal interconnects between cards). Each card hosts two chips: the L-chip, which connects directly to the host via PCIe, and the R-chip, which connects to the L-chip over an Ethernet link (Tenstorrent, 2025e). The simulation is deployed in a multi-host, distributed configuration using tt-run, the TT-Metalium distributed process launcher, which provides a declarative YAML-based interface for MPI-based orchestration (Tenstorrent, 2025c). During the evaluation step, MPI is used to distribute the particle data across the Wormhole devices. Furthermore, OpenMP is employed to parallelize both the prediction and correction steps of the N-body calculation within each MPI process, improving overall scalability and computational throughput. The full implementation of the code can be found in the 6
Figure 2: (Top) Tile-based parallel force calculation. Gray row tiles correspond to the replicated data. Tiles within the vertical lines show the first batch of calculations when, for example, 2 Tensix cores are used. When using more Tensix cores and/or using more Wormhole devices, more tiles are processed in parallel. (Bottom) Pipelined kernel execution illustrating the overlap between data movement and computation.
Figure 3: Schematic of three parallel execution strategies for leveraging Multi-Host and Multi-Chip Wormhole architectures, illustrated for the case of two Wormhole cards: (a) Multi-Host and Single-Chip configuration, where only one chip per card is utilized; (b) Multi-Host and Multi-Chip configuration, where both chips per card are used; and (c) MeshBased configuration, where both chips are organized as MeshDevice instances.
Github repository1 . Scaling Configurations for Multi-Card Execution. As shown in Fig. 3, we consider three approaches to scale the N-body simulation across multiple Wormhole cards and chips: 1 https://github.com/jlalmerol/N-Body-Code.git
7
1. Multi-Host and Single-Chip Configuration: In this configuration, only the L-chip directly connected to the host via PCIe is utilized from each Wormhole card, while the corresponding R-chip remains idle. Parallel execution is achieved by assigning each active chip to a dedicated MPI process. Each MPI process independently initializes its device using the CreateDevices function and allocates interleaved buffers to store the subset of particle data assigned to that process. Global data decomposition is performed using MPI, and the computational workload within each chip is further distributed across Tensix cores using the split work to cores function. 2. Multi-Host and Multi-Chip Configuration: This configuration extends the Multi-Host approach by employing both chips per card within a single MPI process. Separate program instances are defined for each device, including explicit buffer creation and program enqueueing. As a result, each MPI process iterates over its local devices to manage data transfers and kernel execution. This includes explicit data partitioning between the two chips as well as further distribution of the workload across Tensix cores using the split work to cores function. 3. Mesh-Based Configuration: In this configuration, parallel execution across multiple chips is managed by employing the MeshDevice abstraction. Unlike the previous approaches, command queues automatically distribute operations to all devices in the mesh, reducing per-device dispatch overhead and avoiding redundant per-device compilation. In our implementation, a (1, 2) mesh topology is instantiated per MPI process, corresponding to the two chips available on each Wormhole card. Data exchange between devices is managed through sharded buffers for domain-decomposed data and replicated buffers for globally shared particle data. In addition to data decomposition across Wormhole cards using MPI, sharded buffers distribute data across the devices within each card (i.e., the two chips), such that each device operates independently on its assigned subset. Replicated buffers, on the other hand, maintain a complete copy of the data on every device, ensuring that globally shared values are immediately accessible without requiring inter-device communication. 4. Benchmarking of the Wormhole N-Body Application 4.1. Validation To verify the accuracy of the ported N-body application, the calculated acceleration and jerk values are compared against a “golden reference”, namely the naive, double-precision brute-force implementation of the N-body algorithm executed serially on a conventional CPU. Discrepancies are confirmed to be within acceptable tolerances with each acceleration and jerk component deviating by no more than 0.05% and 0.2%, respectively, relative to the golden reference, thereby validating the correctness of the TT-Metalium implementation and the integrity of its hardware operations. A comparison between the “golden reference” and the accelerated codes is shown in Fig. 4, where the energy distributions of the particles at the end of a typical simulation are reported for both cases. The accelerated simulation employs Approach 1; however, comparable distributions are observed when using the other two approaches. 4.2. Performance evaluation Setup. For performance evaluation, the Wormhole-based N-body simulations are benchmarked against an optimized CPU-based reference implementation. The reference code performs the 8
Golden Reference Accelerated
Probability Density
10 3 10 4 10 5 10 6 10 7 0
2000
4000 6000 Particle Kinetic Energy, E
8000
Figure 4: Energy distribution of particles obtained from the Wormhole-accelerated simulation (orange) and the golden reference (gray) after time cycle t = 3. The two distributions illustrate the consistency between the accelerated and the golden reference dataset. All quantities are expressed in normalized units.
evaluation step entirely on the CPU in parallel using MPI and OpenMP, and adopts mixedprecision arithmetic to ensure consistency with the accelerated implementation. Additionally, the computation of acceleration and jerk is further optimized through the use of AVX-512 intrinsics, thereby maximizing single-core performance. We note that a detailed evaluation of the efficiency of the mixed-precision approach relative to a full double-precision implementation is beyond the scope of this work. The reference code is compiled using CMake with the GNU Compiler Collection (GCC) version 13.3.0, conforming to the C++20 standard. Distributed parallelism is enabled through Open MPI version 4.0.6. The build process applies several optimization and performanceoriented flags, including -O3, -fopenmp, -march=native, and -mavx512f. The accelerated version uses the same software environment, except that Open MPI 5.0.7-ULFM replaces version 4.0.6. This version is required as a dependency of the TT-Metal library (v0.62.2), the Tenstorrent Software Development Kit used to operate and program the cards, although our application does not employ any MPI fault-tolerance features. The accelerated code is compiled with the flags -O3, -fopenmp, and -march=x86-64-v3. The host system used to run the N-body simulations for the experimental campaign is equipped with the dual-socket AMD EPYC 9124 processor, offering a total of 64 hardware threads (2 sockets × 16 cores × 2 threads per core) and a maximum clock frequency of 3.71 GHz. The system is equipped with 1.5 TB of DDR5 memory and runs Ubuntu 24.04.2 LTS (64-bit) with Linux kernel version 6.8.0. Four Wormhole n300 devices are connected to the processor via PCIe Gen 4. Performance results. We evaluate the three proposed approaches by measuring the time required to complete a representative simulation (time-to-solution, in seconds) and the energy consumed by the processor and accelerators during its execution (energy-to-solution, in Joule). In this study, we consider simulations modeling 409600 particles evolving over three time steps. This low number of time steps has been chosen to keep the total execution time manageable. While this restricts the total number of simulated time steps, it does not affect the validity of the performance characterization, since the initialization phase of our N-body code represents only a minor fraction of the overall execution time and does not dominate the measured behavior. We execute multiple simulations (typically 20) to obtain sufficient statistical data, launching 9
them in batches. To ensure the system returns to idle conditions between successive runs, a 120-second sleep period is inserted both before and after each simulation. The time-to-solution is measured using hardcoded MPI_Wtime() calls at the beginning and end of the simulation, excluding any time spent in sleep. Table 1 reports the average timeto-solution and the corresponding standard deviation for simulations performed using the three approaches described in Section 3. For reference, in the case of the Multi-Host Single-Chip scenario, we also include the execution time obtained when employing two chips (and thus two cards). All the simulations are launched using tt-run, with each MPI task bound to one device. Additionally the Open MPI flag --bind-to core is used and OMP_NUM_THREADS is set to one. In other words, we use one MPI rank per n300 card, with each MPI rank pinned to a dedicated physical CPU core and exclusively controlling a single n300 card. According to the results in Table 1, the Multi-Host Single-Chip strategy delivers the best performance. Simulations performed on a single chip of an n300 card require, on average, 1459 ± 0.47 s to complete, corroborating the ≈ 2× speedup relative to CPU-only executions reported in Almerol et al. (2025). For the CPU implementation, only the time-to-solution obtained using a shared-memory paradigm is reported, since this approach was verified to achieve superior performance on a single node compared to pure MPI or hybrid MPI+OpenMP configurations. When both chips on the same card are utilized, the execution time increases by approximately 3.6%. Even longer times-to-solution (about 6.58× higher than the Single-Chip configuration) are observed when running the code on two chips and one card in the Mesh-Based configuration. This degradation can be attributed to the additional communication overhead introduced by the Ethernet interconnect between accelerators, which has lower bandwidth and higher latency compared to the PCIe link. Interestingly, when two chips located on different cards are employed through the Multi-Host Single-Chip approach, the simulation completes faster than when using a single card with two active chips, confirming that communication topology plays a crucial role in overall performance. In the Mesh-Based approach, data is distributed across devices using sharded buffers, with subsets of the global particle data assigned to each device. Unlike the first two approaches, which explicitly map data to Tensix cores, the mesh abstraction relies on runtime-managed buffer distribution and command scheduling. This simplifies multi-device programming but introduces overheads from buffer management, shard metadata handling, and synchronized command dispatch. In our experiments, we used device-level sharding to scale the N-body code across multiple devices. While this enables multi-chip execution, it does not fully exploit intra-device parallelism (Tenstorrent, 2025d). We are currently investigating strategies to optimize communication and data movement in the Multi-Host Multi-Chip configuration, as well as ways to better levarege the intra-device parallelism in the Mesh-Based approach to mitigate performance bottlenecks. Scalability. We perform strong scaling tests for the Multi-Host Single-Chip and Multi-Host Multi-Chip approaches, which currently appear more promising than the Mesh-Based Configuration. The results are presented in Fig. 5, which shows the decrease in time-to-solution with increasing number of MPI tasks (left panel) and the corresponding speedup (right panel). Reported values represent averages over at least 16 simulations. For the Multi-Host Single-Chip setup, the number of MPI tasks corresponds to the number of cards and chips used; for the Multi-Host Multi-Chip setup, the number of MPI tasks corresponds to the number of cards, while the number of chips is doubled. Shorter execution times are consistently achieved with the Multi-Host Single-Chip approach, which is on average ≈ 1.04× faster for the same number of MPI tasks. Scalability appears comparable for both strategies, though far from ideal, with 10
Approach
# cards
# chips
Time-to-solution [s]
EDP [kJ s]
1
1
1459.46 ± 0.47
563.10 ± 6.47
2
2
1318.54 ± 2.72
479.43 ± 5.41
1
2
1511.67 ± 0.95
636.84 ± 0.82
1
2
9614.54 ± 1.14
27822.11 ± 52.17
0
0
2875.39 ± 5.30
1927.70 ± 7.30
Multi-Host Single-Chip Multi-Host Single-Chip Multi-Host Multi-Chip Mesh-Based Configuration CPU 32 OMP threads + AVX-512
4
1500 Speedup
Time-to-solution [s]
Table 1: Execution time and EDP for simulations modeling 409600 particles over three time steps. Reported values correspond to the average time-to-solution and EDP over approximately 20 runs for the three approaches described in Section 3. The table also specifies the number of accelerator cards (second column) and processor chips (third column) used in each simulation, and includes, in the last row, the average execution time and EDP for CPU-only simulations. Note that simulations using a single n300 card employ one MPI rank, regardless of the number of chips utilized per card, whereas simulations using two cards and one chip per card in the Multi-Host Single-Chip approach employ two MPI ranks.
1400 1300 1200
2
1 1
2 MPI tasks
4
1
2 MPI tasks
4
Figure 5: Time-to-solution (left panel) and strong scaling parallel speedup (right panel) as functions of the number of MPI tasks for simulations performed using the Multi-Host Single-Chip (green dashed line) and Multi-Host Multi-Chip (orange solid line) configurations. The black dash-dotted line represents the ideal speedup.
speedups and parallel efficiencies of 1.10× and 55%, and 1.16× and 29% at two and four MPI tasks, respectively, for the Multi-Host Single-Chip configuration. These observations indicate that, for both configurations, further optimization of the communication and workload distribution strategies is likely necessary to achieve improved parallel performance. We emphasize that the limited scalability is not caused by inter-card communication overhead, as Wormhole cards do not communicate directly with one another. Rather, the observed scaling limitations are likely due to overheads in the TT-Metalium multi-device orchestration layer, reflecting the still limited maturity of the software stack for multi-card execution. Although our scalability tests were conducted on a single node hosting four cards, the same experiments could be performed on two or four host systems with two or one card each, respectively, without requiring code modifications. However, these configurations would likely exhibit lower scalability due to the higher latency of inter-node communication compared with the single-node setup used in our experiments. Energy-to-solution values are obtained by combining raw and post-processed data sampled in user space at a frequency of ≈ 1 Hz over the entire duration of the job, including both simulation 11
450 Peak power [W]
Energy-to-solution [kJ]
450 425 400 375 350 1
2 MPI tasks
4
400 350 300 250 1
2 MPI tasks
4
Figure 6: Energy-to-solution (left panel) and peak power (right panel) as functions of the number of MPI tasks for simulations performed using the Multi-Host Single-Chip (green dashed line) and Multi-Host Multi-Chip (orange solid line) configurations.
and sleep periods (Amati et al., 2025). Following the methodology outlined in Almerol et al. (2025), we record the timeseries of the power drawn by the eight chips on the four n300 cards, as well as the energy consumed by the two CPUs. For measuring the chip power, we use the Tenstorrent system management interface tt-smi, while CPU energy is obtained via perf stat -a -e with a one-second sampling interval. The energy dissipated by a single n300 chip during a simulation is calculated as the discrete integral of its power over the simulation duration, excluding the sleep periods. These values are then summed across all eight chips to obtain the total energy consumption of the four cards. Similarly, the energy consumed by the dual-socket processor is computed by summing the values recorded with perf, considering only the active simulation time. Conversely, given the energy recorded for each 1 s interval, the instantaneous power absorbed by the CPUs is obtained as the discrete derivative of the energy. The total energy per simulation (energy-to-solution) is then calculated as the sum of the energy consumed by both the four n300 cards and the dual CPUs. Given the energy-to-solution and the time-to-solution, it is also possible to compute the so-called Energy–Delay Product (EDP), defined as the product of the energy-to-solution and the time-to-solution (Amati et al., 2025). Table 1 reports the EDP values for simulations performed using the three implemented approaches, as well as for the reference CPU simulation. The results show that the Multi-Host Single-Chip approach not only minimizes the time-to-solution, but also achieves the lowest EDP. Finally, to determine the overall power, we sum the power values from the n300 cards and the CPUs, performing linear interpolation when measurements are not aligned in time. The left panel of Fig. 6 presents the energy-to-solution as a function of the number of MPI tasks for the Multi-Host Single-Chip and Multi-Host Multi-Chip configurations. As for the timeto-solution, these values represent averages over multiple simulation runs. Unlike the time-tosolution, however, the energy-to-solution does not exhibit a monotonic trend, displaying instead a minimum when two MPI tasks are employed. In particular, using two MPI tasks yields an energy reduction of approximately 10% compared to the configuration achieving the shortest time-to-solution in both cases. The EDP results minimized in this case. This indicates that running with two MPI tasks provides a good compromise between execution time and overall energy efficiency. The right panel of Fig. 6 reports the maximum system peak power observed (excluding memory, storage, networking, cooling system, and other contributions) for simulations using one, two, and four MPI ranks. Employing two chips per card increases the peak power by up to ≈ 8.7% relative to configurations using a single chip per card. 12
5. Summary and perspectives Recently, we developed an N-body astrophysical simulation code and accelerated its force calculation using the Tenstorrent TT-Metal programming interface. By executing the code on the RISC-V-based Tenstorrent Wormhole n300 accelerator, we demonstrated that this class of devices can offer a competitive and energy-efficient solution for next-generation astrophysical computations. In the present work, we further investigated methodologies to enhance code parallelism. In particular, we enabled execution across multiple cards simultaneously by exploiting MPI and card-level parallelism. To the best of our knowledge, this represents the first example of a scientific code capable of running across multiple Tenstorrent devices and cards. Three approaches were implemented to scale simulations across cards and devices: the MultiHost Single-Chip Configuration, the Multi-Host Multi-Chip Configuration, and the Mesh-Based Configuration. Among these, the first approach delivered the best performance in terms of both time-to-solution and energy-to-solution, while the Mesh-Based Configuration yielded the lowest performance. In terms of scalability, the first and second approaches showed comparable behavior, both exhibiting relatively low parallel efficiency at four MPI tasks, suggesting that further optimization is required. Our results also show that while the time-to-solution decreases monotonically with the number of MPI tasks, the energy-to-solution and the enrgy delay product display a minimum when using two MPI tasks, indicating that this configuration provides an effective trade-off between execution time and energy efficiency. As an immediate next step, we plan to benchmark the Tenstorrent results against conventional NVIDIA and AMD GPUs, as well as Xilinx Alveo FPGAs. We have already started developing the force kernel in CUDA. Preliminary results obtained on a system featuring an AMD Genoa CPU and an NVIDIA H100 GPU, using the same parameters described in Section 4, show a speedup of 2.6× with respect to the Multi-Host Single-Chip approach. These results indicate a performance gap with respect to more established accelerator technologies, which can reasonably be attributed to the relative maturity of the current RISC-V–based accelerator platform and its software stack. Further evaluation on newer hardware is therefore required to better assess the evolution of performance, and for this reason we plan to extend our experiments to the Tenstorrent p150 (codename Blackhole™). Since the p150 is built upon the same underlying architecture as the Wormhole series, we anticipate that the strategies developed in this work will transfer effectively to this platform. Acknowledgements We thank E. Duffy, R. Friedman, R. Ganisetti, and F. LeClair (Tenstorrent) for valuable discussions. This research is supported by the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations). E.B. and D.G. acknowledge support from the SPACE project, funded by the European Union. This project has received funding from the European High Performance Computing Joint Undertaking (JU) and from Belgium, the Czech Republic, France, Germany, Greece, Italy, Norway, and Spain under grant agreement No. 101093441. 13
Algorithm 3 Tile-based pairwise acceleration and jerk evaluation within the Compute Kernel 1: Acquire lock on dst register 2: Compute relative displacement ri j = r j − ri 3: Store components r x,i j , ry,i j , rz,i j in intermediate circular buffers (CBs) for reuse 4: Acquire lock on dst register q 5: Compute distance metrics ri2j = |ri j |2 , ri−1 j = 1/
−1 3 ri2j , and ri−3 j = (ri j )
−2 6: Compute scaling coefficients t j = pw, j ri−3 j and Ai j = −3 ri j 7: Compute velocity difference vi j = v j − vi 8: Compute radial velocity component vr,i j = ri j · vi j 9: Compute scaling coefficient qi j = Ai j vr,i j 10: Store t j and qi j in CBs for reuse 11: Acquire lock on dst register
12: Update temporary acceleratio ait += t j ri j 13: Acquire lock on dst register 14: Update temporary jerk ȧit += t j (vi j + qi j ri j )
Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work we used ChatGPT and Grammarly in order to improve the clarity and correctness of the English language. The content of the paper remains entirely original. After using these tools, we reviewed and edited the content as needed. We take full responsibility for the content of the published article. Appendix A. Acceleration and Jerk Calculations within the Compute Kernel Algorithm 3 outlines the tile-based pairwise acceleration and jerk computation implemented within the compute kernel of the Tenstorrent Wormhole accelerator. Intermediate quantities such as the relative displacement, distance metrics, and scaling coefficients are staged into circular buffers (CBs) to enable data reuse and minimize redundant computation across tiles. Synchronization primitives are used to acquire and release locks on destination registers, ensuring safe sequencing between computation and data movement stages. Releasing the destination register before staging data to a CB allows other compute or DMA operations to proceed concurrently, improving pipeline utilization. A small softening parameter, ϵ = 1.0 × 10−7 , is included in the distance computation to avoid numerical singularities when particles are close. References Abbott, B., Abbott, R., Abbott, T., et al., 2020. Prospects for observing and localizing gravitational-wave transients with advanced ligo, advanced virgo and kagra. Living reviews in relativity 23, 1–69. doi:10.1007/ s41114-020-00026-9. Almerol, J.L., Boella, E., Spera, M., Gregori, D., 2025. Accelerating gravitational n-body simulations using the risc-vbased tenstorrent wormhole. arXiv preprint doi:10.48550/arXiv.2509.19294, arXiv:2509.19294. Amati, G., Turisini, M., Monterubbiano, A., Paladino, M., Boella, E., Gregori, D., Croce, D., 2025. Experience on clock rate adjustment for energy-efficient gpu-accelerated real-world codes, in: ISC proceedings in the Springer LNCS series.
14
Bartolini, A., Ficarelli, F., Parisi, E., Beneventi, F., Barchi, F., Gregori, D., Magugliani, F., Cicala, M., Gianfreda, C., Cesarini, D., Acquaviva, A., Benini, L., 2022. Monte cimone: Paving the road for the first generation of riscv high-performance computers, in: 2022 IEEE 35th International System-on-Chip Conference (SOCC), pp. 1–6. doi:10.1109/SOCC56010.2022.9908096. Brown, N., Barton, R., 2024. Accelerating stencils on the tenstorrent grayskull risc-v accelerator. arXiv preprint doi:10. 48550/arXiv.2409.18835, arXiv:2409.18835. Brown, N., Davies, J., LeClair, F., 2025. Exploring fast fourier transforms on the tenstorrent wormhole. arXiv preprint doi:10.48550/arXiv.2506.15437, arXiv:2506.15437. Chang, M., 2025. Programming tenstorrent processors. URL: https://clehaxze.tw/gemlog/2025/ 04-21-programming-tensotrrent-processors.gmi. Corsix, 2024. Tenstorrent wormhole series (part 1: “physicalities” - part 7: “bits of the matmul”). URL: https: //www.corsix.org/content/tt-wh-part1. Einstein Telescope Italy, 2025. Einstein telescope. URL: https://www.einstein-telescope.it/en/ einstein-telescope-en/. Khan, M.A., Al-Mouhamed, M.A., Mohammad, N., 2021. Adaptive tiling for parallel n-body simulations on many core. Astronomy and Computing 36, 100466. Mahale, G., Limbasiya, T., Aleem, M.A., Plana, L., Duricic, A., Monemi, A., Abancens, X., Cervero, T., Davis, J.D., 2023. Optimizations for very long and sparse vector operations on a risc-v vpu: A work-in-progress, in: Bienz, A., Weiland, M., Baboulin, M., Kruse, C. (Eds.), High Performance Computing, Springer Nature Switzerland, Cham. pp. 472–485. Micikevicius, P., Narang, S., Alben, J., Diamos, G., Elsen, E., Garcia, D., Ginsburg, B., Houston, M., Kuchaiev, O., Venkatesh, G., Wu, H., 2018. Mixed precision training. arXiv preprint arXiv:1710.03740 doi:10.48550/arXiv. 1710.03740. Nitadori, K., Makino, J., 2008. Sixth-and eighth-order hermite integrator for n-body simulations. New Astronomy 13, 498–507. doi:https://doi.org/10.1016/j.newast.2008.01.010. Patterson, D., Waterman, A., 2017. The RISC-V Reader: An Open Architecture Atlas. 1st ed., Strawberry Canyon. Rexroth, M., Schäfer, C., Fourestey, G., Kneib, J.P., 2020. High performance computing for gravitational lens modeling: Single vs double precision on gpus and cpus. Astronomy and Computing 30, 100340. doi:https://doi.org/10. 1016/j.ascom.2019.100340. Schäfer, C., Fourestey, G., Kneib, J.P., 2020. Lenstool-hpc: A high performance computing based mass modelling tool for cluster-scale gravitational lenses. Astronomy and Computing 30, 100360. Shukla, N., Romeo, A., Caravita, C., Redenti, M., Vavrik, R., Riha, L., Mignone, A., Rossazza, M., Truzzi, S., Tornatore, L., Ragagnin, A., Castro, T., Karademir, G.S., Dolag, K., Deka, P.J., Bacchini, F., Wilhelm, R.P., Gregori, D., Boella, E., 2025a. Towards exascale computing for astrophysical simulation leveraging the leonardo eurohpc system. Procedia Computer Science 267, 112–123. doi:https://doi.org/10.1016/j.procs.2025.08.238. proceedings of the Third EuroHPC user day. Shukla, N., Romeo, A., Caravita, C., et al., 2025b. Eurohpc space coe: Redesigning scalable parallel astrophysical codes for exascale. invited paper, in: Proceedings of the 22nd ACM International Conference on Computing Frontiers: Workshops and Special Sessions, Association for Computing Machinery, New York, NY, USA. p. 177–184. doi:10. 1145/3706594.3728892. Spera, M., 2014. High Precision, High Performance Simulations of Astrophysical Stellar Systems. Phd thesis. Università degli Studi di Roma ”La Sapienza”. Spurzem, R., 1999. Direct n-body simulations. Journal of Computational and Applied Mathematics 109, 407–432. Suarez, E., Amaya, J., Frank, M., Freyermuth, O., Girone, M., Kostrzewa, B., Pfalzner, S., 2025. Energy efficiency trends in hpc: what high-energy and astrophysicists need to know. Frontiers in Physics Volume 13 - 2025. doi:10. 3389/fphy.2025.1542474. Tenstorrent, 2025a. Dram loopback. URL: https://docs.tenstorrent.com/tt-metal/latest/tt-metalium/ tt_metal/examples/dram_loopback.html. Tenstorrent, 2025b. Kernel apis. URL: https://docs.tenstorrent.com/tt-metal/latest/tt-metalium/tt_ metal/apis/kernel_apis.html. Tenstorrent, 2025c. Programming mesh of devices with tt-nn. URL: https://github.com/tenstorrent/ tt-metal/blob/main/tech_reports/Programming_Mesh_of_Devices/Programming_Mesh_of_ Devices_with_TT-NN.md#4-single-program-multiple-device. Tenstorrent, 2025d. Tensor sharding. https://github.com/tenstorrent/tt-metal/blob/main/tech_ reports/tensor_sharding/tensor_sharding.md. Accessed: 2025-12-15. Tenstorrent, 2025e. Tt architecture and metalium guide. URL: https://github.com/tenstorrent/tt-metal/ blob/main/METALIUM_GUIDE.md. Tenstorrent, 2025f. tt-isa-documentation: Tenstorrent isa documentation. URL: https://github.com/ tenstorrent/tt-isa-documentation.
15
Tenstorrent, 2025g. Wormhole™. URL: https://tenstorrent.com/hardware/wormhole. Venieri, E., Manoni, S., Ceccolini, G., Madella, G., Ficarelli, F., Gregori, D., Cesarini, D., Benini, L., Bartolini, A., 2025. Monte cimone v2: Down the road of risc-v high-performance computers. arXiv preprint doi:10.48550/ arXiv.2503.18543, arXiv:2503.18543.
16