TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Luca Pennatia,∗ , Ivy Penga and Stefano Markidisa a KTH Royal Institute of Technology, Stockholm, Sweden

ABSTRACT

Keywords: Particle-in-Cell Energy-conserving semi-implicit method Space plasma physics Magnetic reconnection GPU acceleration Matrix engines

Artificial heating and numerical damping affect particle energization in kinetic space plasma simulations, even when the numerical scheme remains stable. Implicit and semi implicit energy conserving Particle-in-Cell (PIC) methods eliminate these errors in the discrete energy balance while relaxing explicit stability constraints. We present TPIC, a massively parallel implementation of the Energy Conserving Semi-Implicit Method (ECSIM) for GPU supercomputers. We extend the exact energy conservation in the ECSIM formulation to closed domains with conducting boundaries by deriving and canceling a spurious term in the discrete Poynting flux theorem. TPIC executes the complete particle and field cycle on NVIDIA and AMD GPUs, combining matrix-engine acceleration with asynchronous communication and computation. We validate TPIC in a variety of electrostatic and electrodynamic test cases relevant for space plasma physics applications, including quasineutral plasma expansion, two-stream instability, Whistler waves, Weibel instability, Whistler turbulence, and three-dimensional magnetic reconnection with conducting walls and 12.3 billion particles. We show that in all scenarios, TPIC agrees with the corresponding analytical theory or reference numerical results, and preserves the total system energy to machine precision. On 4,096 AMD GPUs, TPIC advances 1.9 trillion particles with 96.8% weak scaling efficiency.

arXiv:2611.09624v1 [cs.CE] 6 Nov 2026

ARTICLE INFO

1. Introduction

Kinetic processes in space plasmas span a broad hierarchy of spatial and temporal scales, from electron plasma oscillations and Debye shielding to ion gyromotion, magnetic reconnection, and the global dynamics of the solar wind and Earth’s magnetosphere (Lapenta, 2012; Karimabadi et al., 2014). The electron Debye length and skin depth can be orders of magnitude smaller than the characteristic dimensions of current sheets and other macroscopic plasma structures. Nevertheless, these scales are coupled through collective electromagnetic fields and particle dynamics (Lapenta et al., 2006; Kiyani et al., 2015). For instance, magnetic reconnection evolves on ion and macroscopic timescales, while electron kinetics controls important aspects of energy conversion and particle acceleration (Yamada et al., 2010; Hesse and Cassak, 2020; Burch et al., 2016). Resolving this hierarchy with conventional explicit Particle-in-Cell (PIC) methods is computationally demanding, since numerical stability and finite grid constraints impose restrictions on the time step and spatial resolution associated with the fastest and smallest plasma scales, even when the phenomena of interest evolve on much larger scales (Hockney and Eastwood, 1988; Birdsall and Langdon, 1991). Implicit PIC methods relax these numerical stability restrictions, allowing the temporal and spatial resolution to be chosen according to the physical processes of interest. However, numerical stability alone is insufficient for reliable kinetic simulations. Conventional explicit PIC schemes can exhibit artificial numerical heating associated with particle and grid coupling and finite grid instabilities (Langdon, 1970; Okuda, 1972; Barnes and Chacón, 2021). On the other hand, many implicit PIC formulations, including moment-implicit (Mason, 1981; Brackbill and Forslund, 1982) and direct-implicit (Langdon et al., 1983) schemes, introduce numerical damping that leads to artificial cooling (Cohen et al., 1989). These errors can accumulate over long simulations, modifying the particle energy distributions, plasma temperature, and energy exchange between particles and electromagnetic fields. This is particularly problematic in collisionless space plasmas, where energy conversion and particle energization are often central quantities of interest. Energy conserving implicit methods address this problem by enforcing the balance between particle kinetic energy and electromagnetic field energy at the discrete level, eliminating spurious energy gain or loss associated with their numerical coupling (Markidis and Lapenta, 2011; Chen et al., 2011; Lapenta, 2017). ∗ Corresponding author.

pennati@kth.se (L. Pennati) ORCID (s): 0009-0009-4901-1716 (L. Pennati); 0000-0003-4158-3583 (I. Peng); 0000-0003-0639-0639 (S. Markidis)

L. Pennati et al.: Preprint submitted to Elsevier

Page 1 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

The practical application of energy conserving implicit PIC methods to large scale space plasma simulations requires the method to accommodate several important issues. These includes accurate boundary conditions, efficiently handling billions of computational particles, and solving the coupled particle and field equations on massively parallel computing systems. Modern GPU supercomputers offer the computational resources required for these simulations. However, exploiting them efficiently demands algorithms and data structures adapted to accelerator architectures and distributed memory execution. Developing a production quality implementation that combines energy conservation, realistic simulation domains, and scalability across thousands of GPUs is the topic of this paper. Large-scale kinetic simulations of space and astrophysical plasmas have been enabled by massively parallel explicit PIC codes such as VPIC (Bowers et al., 2008; Bird et al., 2022), EMPIC (Wang et al., 1995; Gary et al., 2012), PSC (Germaschewski et al., 2016), OSIRIS (Fonseca et al., 2002; Lee et al., 2025), WarpX (Vay et al., 2018; Fedeli et al., 2022), and TRISTAN-MP (Buneman, 1993; Spitkovsky, 2005) with its successor Tristan v2 (Hakobyan et al., 2024). These codes have been used to study magnetic reconnection, turbulence, collisionless shocks, and kinetic instabilities in the magnetosphere, the solar wind, and high-energy astrophysical environments, with simulations reaching trillions of particles. However, the numerical constraints of explicit schemes remain expensive when the physical processes of interest evolve on ion and macroscopic scales. Implicit and semi-implicit PIC methods offer an alternative (Peng et al., 2015; Innocenti et al., 2017; Chen et al., 2020), although their implementation on GPU supercomputers is less mature. Examples include iPIC3D (Markidis et al., 2010), recently demonstrated at large scale on GPUs (Markidis et al., 2025), the semi-implicit FLEKS (Chen et al., 2023), and implicit electromagnetic solvers in WarpX based on PICNIC (Angus et al., 2023). Energy conserving formulations such as ECSIM provide an attractive foundation for long time kinetic simulations by preserving the discrete energy balance. Further advances are needed to extend energy conserving implicit PIC algorithms to realistic configurations beyond periodic domains and to execute the complete particle and field computations efficiently on distributed GPU systems. In this work, we present TilePIC (TPIC) 1 , an energy conserving semi implicit Particle-in-Cell code designed for large scale kinetic plasma simulations on GPU supercomputers. We extend the ECSIM numerical formulation to conducting boundaries by identifying and eliminating a spurious term in the discrete energy balance. TPIC executes particle and field calculations entirely on GPUs, overlapping communication and computation across distributed NVIDIA and AMD systems. We demonstrate energy conservation to near machine precision in several linear and non-linear test cases, including a three-dimensional magnetic reconnection simulation with conducting walls and 12.3 billion particles. TPIC also achieves 96.8% weak scaling efficiency on 4,096 GPU compute dies, demonstrating its suitability for production scale plasma simulations.

2. Energy-conserving semi-implicit algorithm Implicit PIC algorithms discretize implicitly in time the Maxwell equations, yielding a linear (or more generally nonlinear) formulation where the plasma response to the fields is encoded in a numerical operator. In ECSIM (Lapenta, 2017, 2023), the linear plasma response is represented by the so called mass-matrix operator, and the particle-field coupling is retained through a linear field solve. Let 𝑝 identify a computational particle, 𝑠 its species, 𝑞𝑠 (𝑚𝑠 ) the species charge (mass), 𝑛 a time step, Δ𝑡 the time step size, 𝑔 a mesh node, 𝑉𝑔 its nodal control volume, and 𝐄 (𝐁) the electric (magnetic) field. The computational particle charge and mass, 𝑞𝑝 and 𝑚𝑝 , include its statistical weight, with 𝑞𝑝 ∕𝑚𝑝 = 𝑞𝑠 ∕𝑚𝑠 . The mass-matrix 𝐌𝑔𝑔′ and 𝑛+1∕2

the current density 𝐉𝑔

are defined as

𝑞 Δ𝑡 1 ∑ ∑ 𝛽𝑠 𝑞𝑝 𝛼𝑝𝑖𝑗 𝑊𝑝𝑔 𝑊𝑝𝑔 ′ , 𝛽𝑠 = 𝑠 , 𝑉𝑔 𝑠 2𝑚𝑠 𝑝∈𝑠 ∑ 𝑛+1∕2 𝑛+1∕2 𝐉𝑔 = 𝐉̂𝑔 + 𝐌𝑔𝑔′ 𝐄𝑔′ . 𝑖𝑗 𝑀𝑔𝑔 ′ =

(1)

𝑔′

Here, 𝐌𝑔𝑔′ is a 3 × 3 block coupling the electric field at 𝑔 ′ to the current at 𝑔. The magnetic rotation tensor 𝛼𝑝𝑖𝑗 is 𝜶𝑝𝑛 =

𝐈 − [𝐛𝑝 ]× + 𝐛𝑝 𝐛𝖳𝑝 1 + |𝐛𝑝 |2

,

𝐛𝑝 =

𝛽𝑠 𝑛 𝐁 , 𝑐 𝑝

(2)

1 https://github.com/lucak17/TPIC

L. Pennati et al.: Preprint submitted to Elsevier

Page 2 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

where [𝐛]× 𝐚 = 𝐛 × 𝐚 and 𝑐 is the speed of light. The hatted current 𝐉̂ is the field independent current contribution and is 1 ∑∑ 𝐉̂𝑔 = 𝑞 𝜶 𝑛 𝐯𝑛 𝑊 . (3) 𝑉𝑔 𝑠 𝑝∈𝑠 𝑝 𝑝 𝑝 𝑝𝑔 Substituting the current response into the Maxwell equations gives the linear system [ ] 𝐈 + (𝑐𝜃Δ𝑡)2 ∇ℎ × ∇ℎ × + 4𝜋𝜃Δ𝑡 𝐌 𝐄𝑛+1∕2 = 𝐄𝑛 + 𝑐𝜃Δ𝑡 ∇ℎ × 𝐁𝑛 − 4𝜋𝜃Δ𝑡 𝐉̂, ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⏟⏞⏞⏟⏞⏞⏟ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ implicit plasma response

wave propagation

(4)

explicit Ampère’s law

where, 𝐄𝑛 and 𝐁𝑛 are the known grid fields, 𝐄𝑛+1∕2 is the unknown electric field, 𝐈 is the identity operator, ∇ℎ × is the discrete curl, and 𝜃 = 1∕2 is the time decentering factor. Thus, 𝐄𝑛+1∕2 = (𝐄𝑛+1 + 𝐄𝑛 )∕2. Once it is computed, the fields are advanced by 𝐄𝑛+1 = 2𝐄𝑛+1∕2 − 𝐄𝑛 ,

(5)

𝐁𝑛+1 = 𝐁𝑛 − 𝑐Δ𝑡 ∇ℎ × 𝐄𝑛+1∕2 , where the second equation is the discrete Faraday law. The particle mover is 𝑛+1∕2

𝐱𝑝

𝑛−1∕2

= 𝐱𝑝

+ Δ𝑡 𝐯𝑛𝑝 , 𝑛+1∕2

𝐯𝑝 = 𝜶𝑝𝑛 (𝐯𝑛𝑝 + 𝛽𝑠 𝐄𝑝

(6)

),

= 2𝐯𝑝 − 𝐯𝑛𝑝 , 𝐯𝑛+1 𝑝 ∑ ∑ 𝑛+1∕2 = 𝑉𝑔−1 𝑠 𝑝∈𝑠 𝑞𝑝 𝐯𝑝 𝑊𝑝𝑔 . The where 𝐯𝑝 = (𝐯𝑛+1 + 𝐯𝑛𝑝 )∕2 is the mean velocity defining the deposited current, 𝐉𝑔 𝑝 current response in Eq. 1 follows exactly from the discrete particle mover, assuming the particle positions and magnetic field are fixed. For 𝜃 = 1∕2 and in the absence of particle injection or loss, energy conservation follows from the discrete Poynting theorem since the particle-field energy exchange is balanced exactly, and the change in total energy is determined only by the Poynting flux through the domain boundary, 𝑛+1∕2

( +  )𝑛+1 − ( +  )𝑛 = −Δ𝑡 Φ𝜕Ω

(7)

,

where Ω is the simulation domain and the total particle kinetic, , and electromagnetic field,  , energies are 1∑ 𝑚 |𝐯𝑛 |2 , 𝑛 = 2 𝑝 𝑝 𝑝 ] [ ∑ 1 ∑ 𝑛 𝑛 2 𝑛 2  = 𝑉 |𝐄 | + 𝑉𝑐 |𝐁𝑘 | , 8𝜋 𝑔 𝑔 𝑔 𝑘

(8)

where, 𝑔 labels electric-field nodes and 𝑘 labels magnetic-field cell centers on a uniform mesh of cell volume 𝑉𝑐 . The nodal control volume 𝑉𝑔 equals 𝑉𝑐 in the interior and 𝑉𝑐 ∕2, 𝑉𝑐 ∕4, or 𝑉𝑐 ∕8 on a physical face, edge, or corner, 𝑛+1∕2

respectively. The outward Poynting flux Φ𝜕Ω is evaluated with the midpoint fields 𝐄𝑛+1∕2 and 𝐁𝑛+1∕2 , and it is the boundary remainder of the discrete summation by parts, [ ∑ ] ∑ 𝑐 𝑛+1∕2 𝑛+1∕2 𝑛+1∕2 𝑉𝑐 𝐁𝑘 ⋅ (∇ℎ × 𝐄𝑛+1∕2 )𝑘 − 𝑉𝑔 𝐄𝑔 ⋅ (∇ℎ × 𝐁𝑛+1∕2 )𝑔 . (9) Φ𝜕Ω = 4𝜋 𝑔 𝑘 𝑐 It depends only on boundary values and is the discrete counterpart of ∮𝜕Ω 4𝜋 (𝐄 × 𝐁) ⋅ d𝐀.

𝑛+1∕2

In a periodic domain, there are no boundaries, the summation by parts has no remainder, thus Φ𝜕Ω is trivially zero. At the same time, identical weights for field interpolation and current deposition make the particle and field work terms cancel, ensuring exact energy conservation. Fig. 1 summarizes the particle and field updates in one ECSIM step and highlights the mass-matrix assembly used in the field solve. L. Pennati et al.: Preprint submitted to Elsevier

Page 3 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

2

Particle → grid

Form 𝜶𝑝𝑛 from 𝐁𝑛𝑝

Mass matrix assembly 1 ∑ ∑ 𝑖𝑗 𝑞 𝛼 𝑖𝑗 𝑊 𝑊 ′ 𝛽 𝑀𝑔𝑔 ′ = 𝑉𝑔 𝑠 𝑠 𝑝∈𝑠 𝑝 𝑝 𝑝𝑔 𝑝𝑔

Deposit 𝐉̂𝑔 and 𝐌𝑔𝑔′

1

3

Particle mover

Δ𝑡

Update 𝐱𝑝 and 𝐯𝑝

Field solve

Solve for 𝐄𝑛+1∕2 Update 𝐄𝑛+1 , 𝐁𝑛+1

𝑛 → 𝑛+1

4

Grid → particle

Interpolate 𝐄 and 𝐁 Figure 1: Diagram of the ECSIM (Lapenta, 2017) PIC cycle, with mass matrix deposition highlighted.

3. Exact energy conservation with conductive walls The original ECSIM (Lapenta, 2017) conserves energy exactly only on periodic domains. Here we extend the exact conservation to perfect electric conductor (PEC) and perfect magnetic conductor (PMC) walls, which can be combined with each other and with periodic directions. TPIC uses a uniform Cartesian mesh with all components of 𝐄, 𝐉 and 𝐌 defined at the nodes 𝑔, and all components of 𝐁 at the cell centers 𝑘. The components are not staggered with respect to one another. Numerical derivatives are implemented as centered finite differences, averaged over the four parallel cell edges, e.g., for the node-to-cell curl  (𝜕𝑥ℎ 𝑓 )𝑖+ 1 ,𝑗+ 1 ,𝑘+ 1 = 2

2

2

∑ ( ) 1 𝑓 − 𝑓𝑖,𝑗+𝑎,𝑘+𝑏 , 4Δ𝑥 𝑎,𝑏∈{0,1} 𝑖+1,𝑗+𝑎,𝑘+𝑏

(10)

and analogously for the cell-to-node curl, denoted  † since it is the adjoint of , (the proof is reported in Appendix A). With the cell volume 𝑉𝑐 = Δ𝑥 Δ𝑦 Δ𝑧 and the nodal control volume 𝑉𝑔 defined in Section 2, nodal and cell inner products are ∑ ∑ 𝐚𝑘 ⋅ 𝐛 𝑘 . (11) 𝑉𝑔 𝐚𝑔 ⋅ 𝐛𝑔 , ⟨𝐚, 𝐛⟩𝐶 = 𝑉𝑐 ⟨𝐚, 𝐛⟩𝑉 = 𝑔

𝑘

This discretization retains mimetic properties of the continuum operators, ∇ℎ ⋅ (∇ℎ × 𝐀) = 0 holds exactly, and on ∑ ∑ periodic meshes 𝑘 𝑉𝑐 𝐁𝑘 ⋅ (𝐀)𝑘 = 𝑔 𝑉𝑔 𝐀𝑔 ⋅ ( † 𝐁)𝑔 , which is the discrete counterpart of ∫Ω ∇ ⋅ (𝐀 × 𝐁) 𝑑𝑉 = 𝑛+1∕2

∫Ω (𝐁 ⋅ ∇ × 𝐀 − 𝐀 ⋅ ∇ × 𝐁) 𝑑𝑉 = 0, and ensures the Poynting flux Φ𝜕Ω in Eq. 9 to vanish. In Eqs. 4 and 5, ∇ℎ × stands for  acting on 𝐄, and  † acting on 𝐁. In the discrete scheme, the Poynting flux Φ𝜕Ω is the summation-by-parts remainder of Eq. 9. In a closed domain, at a wall node,  † uses a ghost cell outside the domain, and, using the numerical differential operator defined in Eq. 10, we can express Eq. 9 as Φ𝜕Ω =

[ ] ( ) ] Δ 𝑐 [ 𝑐 ∑ ⟨𝐄, 𝐁⟩𝐶 − ⟨𝐄,  † 𝐁⟩𝑉 = 𝐴𝑔 (𝐄𝑔 × 𝐁𝑔 ) ⋅ 𝐧̂ + 𝑛 𝐄𝑔 ⋅ ∇ℎ𝑡 × [[𝐁]] 𝑔 . 4𝜋 4𝜋 𝑔∈𝜕Ω 4

(12)

Here, 𝐧̂ is the outward normal, Δ𝑛 the mesh spacing along it, 𝐴𝑔 = 2𝑉𝑔 ∕Δ𝑛 the wall area of the control volume of node 𝑔, and ∇ℎ𝑡 × the cell-to-node curl restricted to wall-tangential derivatives. 𝐁 and [[𝐁]] are the mean and the jump of 𝐁 between the first interior cell and the ghost cell, interpolated to the wall nodes, as shown in Fig. 2. We report the full derivation of Eq. 12 in Appendix A.

L. Pennati et al.: Preprint submitted to Elsevier

Page 4 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

ghost

domain 𝐄0

𝐄−1 𝐁−1∕2

𝐄1 𝐁1∕2

𝐁3∕2 𝐁3∕2

𝐁1∕2 𝐁

𝐄2

[[𝐁]]

𝐁−1∕2

Figure 2: Cell and node field values on a 1D discretized domain, including ghost points. The electric field 𝐄 lives on nodes, while the magnetic field 𝐵 is cell-centered. 𝐁 and [[𝐁]] are the mean and the jump of 𝐁 between the first interior cell and the ghost cell.

Table 1 Ghost parities 𝜎 of the wall-tangential (𝑡) and normal (𝑛) field components. Electric components marked 0 are set to zero at the wall nodes. Wall 𝐄𝑡 𝐸𝑛 𝐁𝑡 𝐵𝑛 PEC −1, 0 +1 +1 −1 PMC +1 −1, 0 −1 +1

The first term in Eq. 12 is the discrete Poynting flux. The second is a discretization artifact without a continuous counterpart, coupling 𝐸𝑛 to [[𝐁𝑡 ]] and 𝐄𝑡 to [[𝐵𝑛 ]]. In a closed domain, exact energy conservation requires both terms to vanish. TPIC cancels them with two discrete constraints. First, we set to zero at the wall nodes the electric components that vanish on the wall, 𝐄𝑡 for PEC and 𝐸𝑛 for PMC. Second, following the image field theory (Taflove and Hagness, 2005; Ricci et al., 2002), each ghost cell holds the reflection of the adjacent interior cell, 𝐁−1∕2 = 𝜎𝐁1∕2 componentwise, with parity 𝜎 = ±1, so that 𝐁=

1+𝜎 𝐁1∕2 , 2

[[𝐁]] = (1 − 𝜎) 𝐁1∕2 .

(13)

An odd parity cancels the mean, and an even parity cancels the jump. The parities in Table 1 cancel the terms that the electric constraint leaves. On a PEC wall, 𝐄𝑡 = 0 removes the flux and the [[𝐵𝑛 ]] term, and an even 𝐁𝑡 removes the [[𝐁𝑡 ]] term. On a PMC wall, 𝐸𝑛 = 0 removes the [[𝐁𝑡 ]] term, and since 𝐄𝑡 ≠ 0, an odd 𝐁𝑡 removes the flux and an even 𝐵𝑛 the [[𝐵𝑛 ]] term. The odd 𝐵𝑛 of PEC is not required by Eq. 12, but it makes 𝐵 𝑛 = 0, so that both walls satisfy the discrete counterparts of 𝐧̂ × 𝐄 = 𝐧̂ ⋅ 𝐁 = 0 (PEC) and 𝐧̂ × 𝐁 = 𝐧̂ ⋅ 𝐄 = 0 (PMC). The same closure applies to 𝐯. Nodal ghosts use the parities of 𝐄, and edge and corner ghosts take the product of the parities of the adjacent walls. With the residual of Eq. 12 cancelled, ⟨𝐄,  † 𝐁⟩𝑉 = ⟨𝐄, 𝐁⟩𝐶 , and the discrete Poynting flux through the domain is exactly zero, thus the system energy is conserved to the linear solver tolerance and round-off. With non periodic boundary conditions, particles are specularly reflected at the walls, which leaves |𝐯𝑝 | unchanged and thus preserves the energy.

4. Numerical methods in TPIC In this section we detail the numerical methods implemented in TPIC, including the accuracy of the discretization, the energy-conserving smoothing, and the matrix-free field solver. L. Pennati et al.: Preprint submitted to Elsevier

Page 5 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

Discretization and accuracy The centered differences of Eq. 10 are second-order accurate. Particles use trilinear cloud-in-cell (CIC) shape functions 𝑊𝑝𝑔 for interpolation operations, and 𝐁 is averaged from the eight neighboring cells to the nodes before the gather. TPIC spatial discretization is therefore second order accurate. For 𝜃 = 1∕2, the field update is the implicit-midpoint (Crank-Nicolson) rule, which is second-order accurate, and the particle push of Eq. 6 is time-centred in 𝐄 and 𝐱. The rotation 𝜶𝑝 uses 𝐁𝑛 instead of 𝐁𝑛+1∕2 . This half-step lag keeps Eq. 4 linear and does not affect energy conservation since the magnetic force does no work, but it introduces an error proportional to Δ𝑡 𝜕𝑡 𝐁. Thus, the coupled scheme is second-order time accurate for static magnetic fields and first-order otherwise. The implicit field solve removes the CFL and 𝜔𝑝𝑒 Δ𝑡 stability limits, and energy conservation suppresses the finite-grid instability, so Δ𝑥 does not need to resolve the Debye length.

Energy-conserving smoothing Following Lapenta (2023), we implement smoothing in TPIC to filter particle and grid noise. Critically, to retain exact energy conservation, smoothing must be applied both to particle moments and fields. Denoting the smoothing operator by , the smoothed 𝐄 field is 𝐄𝑛+1∕2 , used in the particles push. By Eq. 1, this yields the current 𝐉𝑛+1∕2 = 𝐉̂ + 𝐌𝐄𝑛+1∕2 , and Ampère’s law then uses 𝐉𝑛+1∕2 . The work done on the particles, Δ𝑡⟨𝐄𝑛+1∕2 , 𝐉𝑛+1∕2 ⟩𝑉 , equals the work done by the field, Δ𝑡⟨𝐄𝑛+1∕2 , 𝐉𝑛+1∕2 ⟩𝑉 , provided that  is self-adjoint in the nodal inner product of Eq. 11. We implement two smoothing operators in TPIC. The seven-point operator (7 𝑓 )𝑔 = 𝛼𝑓𝑔 +

1−𝛼 ∑ 𝑓 ′, 6 𝑔′ ∈ 𝑔

(14)

𝑔

where 𝑔 contains the six face neighbors of node 𝑔, has a variable center weight 0 ≤ 𝛼 < 1. The 27-point operator is the tensor product of one-dimensional binomial filters, 𝑦 𝑥 𝑧 27 = 121 121 121 ,

) ( 𝑥 (121 𝑓 )𝑖,𝑗,𝑘 = 14 𝑓𝑖−1,𝑗,𝑘 + 2𝑓𝑖,𝑗,𝑘 + 𝑓𝑖+1,𝑗,𝑘 ,

(15)

𝑦 𝑧 defined in similar way. Expanding the product,  assigns the fixed weights 8∕64 to node 𝑔 and with 121 and 121 27 4∕64, 2∕64, and 1∕64 to each of its 6 face, 12 edge, and 8 corner neighbors, respectively. Both operators have symmetric stencils, are self-adjoint, and are second-order accurate. Both operators are applied 𝑁𝑠 times.

Matrix-free GMRES With the smoothing operator  previously defined ( =  without smoothing), the Maxwell system in Eq. 4 takes the form of 𝐄𝑛+1∕2 = 𝐛, with  =  + (𝑐𝜃Δ𝑡)2  †  + 4𝜋𝜃Δ𝑡 𝐌, ̂ 𝐛 = 𝐄𝑛 + 𝑐𝜃Δ𝑡  † 𝐁𝑛 − 4𝜋𝜃Δ𝑡  𝐉.

(16)

The antisymmetric part of 𝜶𝑝 makes 𝐌, and hence , non-symmetric, so TPIC solves it with a matrix-free restarted GMRES, where the operator  is applied without assembling the corresponding matrix. The Krylov basis, the Hessenberg matrix, the Givens rotations, and the least-squares solve reside on the device. A solution is accepted only if the true residual satisfies ‖𝐛−𝐄𝑛+1∕2 ‖2 ≤ 𝜀‖𝐛‖2 , with 𝜖 the required tolerance. Boundary conditions are enforced by solving the linear system in the constrained subspace, where the right-hand side and the initial guess are projected, and the constrained rows of the operator are replaced by the identity. The mass matrix 𝐌 is diagonally dominant, making the operator  well conditioned, thus the solver generally converges in 10 − 30 iterations and no preconditioner is needed.

5. TPIC architecture TPIC is a modern C++20 code designed for heterogeneous architectures and large-scale systems. It is GPUnative, with the entire core algorithm executed on the device. This design keeps data fully device-resident throughout the simulation, minimizing expensive host-device data transfers. The host has primarily an orchestration role and is responsible for I/O operations. Parallelism is exploited at both the data and task levels. Data-level parallelism is achieved through on-node GPU execution, while task-level parallelism is implemented using C++ threads to coordinate multiple tasks concurrently across multiple device streams. For distributed-memory simulations, the computational L. Pennati et al.: Preprint submitted to Elsevier

Page 6 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

domain is partitioned into subdomains, each assigned to a distinct process, and TPIC relies on communication libraries such as MPI, NCCL, and RCCL. TPIC builds on our custom Light Particle Field Accelerator (LPFx) library, which provides performance portability across NVIDIA and AMD GPU backends and multiple communication stacks, by implementing the distributedmemory data structures and execution framework TPIC uses. In this section we describe TPIC’s design, and we discuss the key architectural choices that enable high performance simulations.

5.1. Asynchronous pipeline TPIC executes each time step as a set of tasks distributed over persistent host worker threads, implemented in LPFx with native C++ threads. Each worker owns dedicated device streams, its communication state, and the data it produces. The workers are the particle mover 𝑀, the field advance 𝐹 , one local worker 𝐿𝑠 for each species 𝑠, the particle migration 𝐶, the moment deposition 𝐷, the Maxwell solver 𝑆, the diagnostics 𝑄, and the output writer 𝑊 . A control thread submits each time step to the workers. The host threads mainly enqueue device work and drive communication, so that the device streams execute concurrently. Dependencies between workers are explicit. TPIC uses device events, recorded on the producer stream and waited on the consumer stream, to order device work without host synchronization. Host fences notify a host thread that a task has been issued, e.g., before it starts a dependent communication. The pipeline fixes the order of communication operations, ensuring all ranks issue communication calls in the same sequence, a strict requirement for NCCL and RCCL backends. Specifically, halo exchange precedes migration, species are migrated in index order, and 𝐉̂ is assembled before 𝐌. Fig. 3 schematizes one TPIC time step with two species. The step starts from the electric field computed by the previous solve, i.e. 𝐄𝑛−1∕2 . Worker 𝑀 pushes the species in order, while 𝐹 advances 𝐄 and 𝐁 with Eq. 5 and interpolates the new 𝐁 to the nodes, later used by 𝐷 to evaluate 𝜶𝑝 . Once species 𝑠 is pushed, 𝐿𝑠 separates the particles leaving the subdomain, which 𝐶 migrates on a high-priority stream, from those staying local, which 𝐿𝑠 sorts concurrently. With partial sorting (Section 5.3), 𝐷 deposits 𝐉̂ and 𝐌 on two separate streams, first from the sorted local particles and then from the arrivals, which 𝐿𝑠 appends once the migration completes. After all species are migrated, 𝐷 completes 𝐉̂ and 𝐌 assembly with halo exchange across subdomains. Worker 𝑆 starts the field solve before the deposition completes. It computes the field terms of the right-hand side as soon as the new state is available, completes it once 𝐉̂ is assembled, and starts GMRES on the interior entries of 𝐌 while their ghost contributions are still being exchanged. It then calculates the field 𝐄𝑛+1∕2 for the next particle push. Meanwhile, 𝑄 computes local diagnostics and captures the output snapshot, and it starts its collective reductions only after 𝑆 completes, so the two workers never interleave communication. The next step starts when both 𝑆 and 𝑄 finish, while 𝑊 writes the captured snapshot to disk. In this way, migration, moment assembly, diagnostics, and output overlap with computation, avoiding device idling and maximizing compute resources utilization.

5.2. Portable data structures and communication with LPFx

To provide performance portability across different computing and communication backends, we develop LPFx, a header-only C++20 library for distributed-memory simulations. LPFx provides host and device data structures, algorithms, and communication patterns commonly required by PIC methods, while separating numerical kernels from storage and transport details. It supports CUDA and HIP devices and provides transport implementations based on MPI, NCCL, and RCCL. TPIC leverages LPFx to ensure a single code source portable across multiple backends. Fig. 4 shows the resulting software stack. LPFx represents multi-component fields on uniform 3D Cartesian grids, including per-component staggering and halo regions. Field components reside at planned offsets within a contiguous allocation and may use backend-specific pitch and alignment policies. Particles are described by typed schemas and stored in a structure-of-arrays (SoA) layout, where each array holds a different particle quantity. Schemas, field locations, communication patterns, and backend choices are resolved at compile time, eliminating runtime dispatch from performance-critical paths. LPFx provides communication patterns such as copy_to_ghost and ghost_to_owner_add for fields, particle_migratio for particles, collective operations, and fixed-message exchanges. Each operation is defined at compile time and prepared using runtime topology, geometry, and capacity information. The resulting prepared plan is bound to a data owner, a transport backend, and an execution stream. This modularity allows the same data structure and numerical algorithm to be used with different communication backends. L. Pennati et al.: Preprint submitted to Elsevier

Page 7 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers step 𝑛

step 𝑛+1

𝐿0 species 0

sort local

𝐿1 species 1 𝑀 mover

sort local 𝑠0

arrivals 𝑠0

𝑠1 𝑠0

𝐶 migration

𝑠1 𝑠0 local 𝑠1 local

𝐷 deposition 𝐹 field advance

arrivals

𝑠0 arrivals

𝑠1 arrivals

𝐉̂

RHS GMRES

𝑆 Maxwell

smooth 𝐄𝑛+1∕2

reduce capture

𝑄 diagnostics 𝑊 output writer species 𝑠0

𝐌

𝐄𝑛 , 𝐁𝑛

finalize

write previous snapshot to disk species 𝑠1

migration

fields

moments

diagnostics, output

dependency

Figure 3: Asynchronous TPIC pipeline for two particle species, with one row per host worker. Dotted arrows mark the main dependencies between workers. Worker 𝐶 migrates each species once it is pushed, while 𝐿𝑠 concurrently sorts the particles that stay local. 𝐷 deposits local particles as soon as 𝐿𝑠 sorts them, then deposits the arrivals, which 𝐿𝑠 finalizes once the migration completes. The diagnostics worker 𝑄 overlaps its work with 𝐷 and 𝑆, and the next step starts once 𝑆 and 𝑄 complete. The disk writer 𝑊 writes field and particle snapshots, and checkpoints while the next step runs.

TPIC LPFx NCCL CUDA

MPI

RCCL HIP/ROCm

Figure 4: TPIC software stack. LPFx wraps the CUDA and HIP/ROCm device backends and the communication transports, NCCL for NVIDIA GPUs, RCCL for AMD GPUs, and MPI for both.

LPFx minimizes message count by aggregating data into compact communication images. Field exchange packs strided subregions into dense contiguous buffers. Particle migration uses an explicit SoA wire format containing velocities, charge, identifier, and canonical position. Both payload types are transported as uninterpreted bytes, with one message per peer. Prepared operations retain their communication images and auxiliary storage for packing, routing, and unpacking. These resources are reused across simulation steps, avoiding repeated allocation/deallocation operations. LPFx supports device-direct MPI, staged MPI through persistent pinned host buffers, NCCL, and RCCL. The transport backend is selected independently for each binding, and different transport strategies can coexist within one application. In the device-aware MPI configuration, TPIC uses direct device communication for fields and particles, and other device-resident operations. Some MPI implementations may internally stage device collectives through pageable host memory, even when device-aware MPI is enabled. Thus, for TPIC’s reductions, we explicitly select LPFx’s staged-MPI backend, which uses persistent pinned host buffers and avoids repeated staging-buffer allocation. In NCCL or RCCL configurations, the corresponding vendor backend is used for device payloads, while MPI remains responsible for bootstrap and control. Communication operations expose an explicit lifecycle through begin(), progress() or wait(), and conclude(). Each binding can use a dedicated device stream for packing, transport preparation, and receive processing. This interface allows TPIC to overlap field exchange and particle migration with independent computation. LPFx is transparent, avoiding automatic data transfer or synchronization. It provides the required stream and event mechanisms for explicit synchronization, while TPIC retains full control over pipeline ordering and the data lifetime. L. Pennati et al.: Preprint submitted to Elsevier

Page 8 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

Portability between NVIDIA and AMD accelerators is implemented through a lightweight compile-time abstraction over CUDA and HIP. Backend traits map memory allocation, copies, streams, events, and device operations directly to the corresponding native runtime. Similarly, warp-level instructions such as lpfx::warp_shuffle_up compile to __shfl_up_sync on CUDA and __shfl_up on HIP, without runtime overhead. TPIC can therefore use a single source implementation while retaining access to low-level operations required by performance-critical kernels.

5.3. Particle sorting In TPIC, particles of each species are sorted independently by mesh cell within each subdomain. In the SoA layout, this groups particles belonging to the same cell contiguously in memory within each particle-attribute array. Particle sorting is key to achieving high performance, as it enables optimized cell-based kernels that leverage coalesced memory accesses, cache field values reused across multiple particles, and accumulate particle quantities within each cell with atomic-free deposition. The sorting algorithm is implemented in the LPFx library as a direct counting sort using the integer cell index stored with each particle. The entire sorting pipeline runs on GPU. Let 𝑁𝑝 denote the number of particles of one species to sort and 𝑁𝑐 the number of cells in the local index space. The sorting operation first builds a histogram ℎ𝑐 by atomically counting particles in cell 𝑐. An exclusive prefix sum then determines the beginning of each cell’s range. Let 𝑜𝑐 be the 𝑐 cell’s offset, 𝑜0 = 0,

𝑜𝑐 =

𝑐−1 ∑

ℎ𝑏

(1 ≤ 𝑐 ≤ 𝑁𝑐 ),

𝑜𝑁𝑐 = 𝑁𝑝 .

(17)

𝑏=0

Cell 𝑐 occupies the half-open interval [𝑜𝑐 , 𝑜𝑐+1 ), whose length is exactly its particle population. Empty cells have coincident offsets. LPFx computes the offsets with a hierarchical Blelloch scan (Blelloch, 1993) implemented with GPU kernels. The scan requires 𝑂(𝑁𝑐 ) work and 𝑂(log 𝑁𝑐 ) parallel depth. For each cell 𝑐, a separate index is initialized to 𝑜𝑐 . Particles atomically advance their cell’s index to obtain unique destination indices, defining a permutation that is applied to every particle attribute. Particles’ attributes are then scattered into new buffers according to the permutation to obtain the sorted storage. Thus, the complete sort requires 𝑂(𝑁𝑝 + 𝑁𝑐 ) work. The attributes permutation uses a rotating storage scheme to reduce the overall memory footprint. Each attribute is scattered in turn into a new available buffer, and its previous buffer becomes the destination for the next attribute. The attribute traversal order accounts for their different storage widths. Thus, despite the scattering method transfers the same amount of data as a conventional double buffer approach, in TPIC the additional particle storage is only one buffer sized for the widest attribute at the allocated particle capacity. TPIC supports both full and partial sorting pipelines. In full sorting, particle migration across subdomains completes before all locally owned particles, including the new incoming ones, are sorted. All subsequent stages then process particles using the optimized cell-based kernels. In partial sorting, particles that remain on the local rank are compacted into a resident prefix and sorted while migration of outgoing/incoming particles concurrently progresses. As soon as local sorting completes, next stages (current and mass-matrix deposition) can consume this sorted prefix while remote particles are still in transit. Sorted particles are processed with optimized cell-based kernels, while arriving particles are appended as an unsorted tail and processed separately with unsorted kernels. Partial sorting therefore overlaps resident sorting and deposition with communication, at the cost of processing the received tail without cell grouping during that step. Sorting is repeated at each cycle since particle push updates particle position, invalidating the previous cell ranges. In partial sorting, particles in the unsorted tail at time step 𝑛 become sorted particles at time step 𝑛 + 1, provided they remain within the local subdomain. As a result, the number of unsorted particles remains approximately constant throughout the simulation and is equal to the number of migrated particles, thus being only a small fraction of the total particle population.

5.4. Mass-matrix deposition on tensor-cores Particle-to-grid scatter operations are among the most expensive operations in PIC algorithms since they involve 𝑖𝑗 fine-grained reductions that hinder parallelization. In ECSIM, each mass-matrix block 𝑀𝑔𝑔 ′ couples field component ′ 𝑗 at node 𝑔 to current component 𝑖 at node 𝑔. Its assembly is therefore more expensive than standard particle-tonode deposition, since each particle contributes to every node pair in its interpolation stencil. For three-dimensional cloud-in-cell (CIC) shape functions, this yields 8 × 8 × 9 = 576 scalar contributions per particle, or 36 × 9 = 324 L. Pennati et al.: Preprint submitted to Elsevier

Page 9 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

distinct coefficients after exploiting spatial symmetry. Since multiple particles update the same global entries, a naive GPU implementation parallelized over particles requires extensive atomic operations in main memory. The resulting fine-grained atomic traffic makes mass-matrix assembly a major computational bottleneck. Following Pennati and Markidis (2026), we express the mass-matrix assembly as a tensor contraction over the particle index, and use particle sorting to decompose it into dense cell-local products. Pennati and Markidis (2026) derive this decomposition for general shape functions. With CIC, all particles in a cell couple to the same eight corners, so each cell contributes one dense 8×8 product. Let 𝑐 contain the particles in cell 𝑐, 𝑤𝑎𝑝 be the CIC weight of particle 𝑝 at the local corner 𝑎 ∈ {0, … , 7}, and 𝜅𝑝𝑖𝑗 = 𝛽𝑠(𝑝) 𝑞𝑝 𝛼𝑝𝑖𝑗 be the particle response coefficient from Eq. 1, with 𝑠(𝑝) the species of 𝑝. Particles are processed in batches 𝓁 of four, 𝑝𝓁,𝑘 with 𝑘 = 0, … , 3, which form the contraction dimension of one matrix-engine instruction. The cell-local mass matrix is then ∑ ∑( ) 𝑖𝑗,𝑐 𝑀𝑎𝑏 = 𝑤𝑎𝑝 𝜅𝑝𝑖𝑗 𝑤𝑏𝑝 = 𝐴𝑖𝑗,𝓁 𝐵 𝓁 𝑎𝑏 , 𝑝∈𝑐 𝓁 (18) 𝑖𝑗,𝓁 𝑖𝑗 𝓁 𝐴𝑎𝑘 = 𝑤𝑎𝑝𝓁,𝑘 𝜅𝑝 , 𝐵𝑘𝑏 = 𝑤𝑏𝑝𝓁,𝑘 , 𝓁,𝑘

where 𝑘 is the contraction index over particles. Each batch performs one update 𝐷𝑖𝑗 ← 𝐷𝑖𝑗 + 𝐴𝑖𝑗,𝓁 𝐵 𝓁 of an 8 × 8 𝑖𝑗,𝑐 , with 𝑔𝑐 (𝑎) the global node of corner 𝑎. accumulator, and the deposition adds 𝑀𝑔𝑖𝑗 (𝑎),𝑔 (𝑏) += 𝑀𝑎𝑏 𝑐

𝑐

𝑖𝑗,𝑐 𝑖𝑗,𝑐 For fixed (𝑖, 𝑗), the mass-matrix is symmetric, 𝑀𝑎𝑏 = 𝑀𝑏𝑎 , thus only the 36 canonical entries of each local 8 × 8 matrix are written. Contributions from adjacent cells can target the same global entry and are therefore combined atomically, but accumulation within a cell is completed in the registers before this flush. Thus, the number of global atomic updates is bounded by 36 × 9 per occupied cell and does not scale with its particle population. On NVIDIA GPUs, TPIC follows Pennati and Markidis (2026) and computes each 8 × 8 update with one FP64 WMMA instruction with tile (𝑀𝑡 , 𝑁𝑡 , 𝐾𝑡 ) = (8, 8, 4). Here, we improve upon Pennati and Markidis (2026) by adding AMD support, so TPIC executes mass-matrix assembly on the matrix engines of both NVIDIA and AMD GPUs. AMD GPUs do not have an instruction equivalent to the NVIDIA FP64 8 × 8 × 4 tensor core operation, but the HIP FP64 MFMA builtin __builtin_amdgcn_mfma_f64_4x4x4f64 can perform four independent products with tiles (𝑀𝑡 , 𝑁𝑡 , 𝐾𝑡 ) = (4, 4, 4) across a 64-lane wavefront in one single instruction. Thus, on HIP, TPIC partitions the eight cell corners into two sets, 𝐿 = {0, 1, 2, 3} and 𝐻 = {4, 5, 6, 7}, which split the rows of 𝐴𝑖𝑗,𝓁 into 𝐴𝐿 and 𝐴𝐻 and the columns of 𝐵 𝓁 into 𝐵𝐿 and 𝐵𝐻 . The 8 × 8 cell-local mass matrix is thus divided in four quadrants (𝐿, 𝐿), (𝐿, 𝐻), (𝐻, 𝐿), and (𝐻, 𝐻), computed by the four products 𝐴𝐿 𝐵𝐿 , 𝐴𝐿 𝐵𝐻 , 𝐴𝐻 𝐵𝐿 , and 𝐴𝐻 𝐵𝐻 , as illustrated in Fig. 5. One MFMA instruction executes the four 4 × 4 matrix products concurrently, and therefore covers the same 8 × 8 logical result as one NVIDIA MMA instruction. The hardware blocks are interleaved across wavefront lanes, and a shuffle obtains the partner-corner weight for each lane’s second operand. One wavefront owns a cell, with four wavefronts in a 256-thread block. Each lane loads one particle per 64-particle chunk. The wavefront processes that chunk in four-particle batches and keeps one FP64 accumulator in register per lane and tensor component, for nine persistent values. The final batch is zero padded, and all 64 lanes participate in each MFMA instruction. The operands 𝐴𝑖𝑗,𝓁 and 𝐵 𝓁 of Eq. 18 have 32 entries each, but the four blocks read 64 values of each operand, one per lane. Since 𝐴𝐿 enters the (𝐿, 𝐿) and (𝐿, 𝐻) blocks, 𝐴𝐻 the (𝐻, 𝐿) and (𝐻, 𝐻) blocks, and likewise for 𝐵𝐿 and 𝐵𝐻 , each entry is held by two lanes. The register accumulators persist during the entire cell traverse, then quantities are written in main memory with an atomic operation.

5.5. Mixed precision TPIC separates, for each quantity, the storage type used in memory from the compute type of the arithmetic that ̂ and 𝐌, produces it. All kernels are templates on storage and compute types for particle positions and velocities, 𝐄, 𝐁, 𝐉, and every conversion between the two is explicit. When storage and compute type coincide, no type casting instruction is issued. At runtime, the pair is selected from those compiled, by default, it matches the five levels in Table 2. Particle charges, GMRES inner products and norms, and diagnostic reductions are always computed in FP64. Single-precision storage halves the memory footprint and traffic for the corresponding arrays, which is most relevant for particles. Particle positions are stored as cell-relative coordinates in [0, 1), together with an integer cell index. FP32 positions thus have a resolution of at least 2−24 of the cell size, independently of the domain size. Since positions enter the energy balance only through the interpolation weights, which deposition and push compute in FP64 from the same stored coordinates, single-precision positions preserve exact energy conservation (Section 7.7). Lower precision for velocities, fields, or moments instead introduces single-precision round-off in the energy balance. L. Pennati et al.: Preprint submitted to Elsevier

Page 10 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers 𝐷𝑖𝑗 ← 𝐷𝑖𝑗 + 𝐴𝑖𝑗,𝓁 𝐵 𝓁

(one MFMA, four 4 × 4 × 4 blocks)

𝐴𝑖𝑗,𝓁

𝐷𝑖𝑗 𝐵𝓁

𝐴𝐿 0

×

1

𝐵𝐿

𝐵𝐻

2 3

1

2

𝛽=1 𝐴𝐿 𝐵𝐻

𝛽=2 𝐴𝐻 𝐵𝐿

𝛽=3 𝐴 𝐻 𝐵𝐻

⟶

𝐴𝐻

0

𝛽=0 𝐴𝐿 𝐵𝐿

3

8 corners × 4 particles 𝑘

8 × 8 corner pairs

4 particles 𝑘 × 8 corners

Figure 5: Cell-local FP64 MFMA update on AMD GPUs for one tensor component (𝑖, 𝑗) and four-particle batch 𝓁. Each column of 𝐴𝑖𝑗,𝓁 contains a particle’s eight CIC weights multiplied by its response coefficient, and the matching row of 𝐵 𝓁 contains the unscaled weights. The corner sets 𝐿 = {0, 1, 2, 3} and 𝐻 = {4, 5, 6, 7} split the product into four independent 4 × 4 × 4 blocks 𝛽 of one MFMA instruction, each updating one quadrant of the 8 × 8 accumulator 𝐷𝑖𝑗 .

Table 2 TPIC precision levels, with the floating-point width (64 or 32 bits) of particle positions 𝐱𝑝 and velocities 𝐯𝑝 , fields, current, and mass matrix. Storage and compute types coincide, except for the mass matrix of the fourth level, stored in FP32 and computed in FP64 (32/64). Particle charges and all reductions are always FP64. Level

𝐱𝑝

𝐯𝑝

𝐄, 𝐁

𝐉̂

𝐌

double single-position single-position-velocity single-particles-moments-fields-double-mass single-particles-moments-fields

64 32 32 32 32

64 64 32 32 32

64 64 64 32 32

64 64 64 32 32

64 64 64 32/64 32

In GMRES, the Krylov basis and the operator applications use the compute type of 𝐄, while the solution, the right-hand side, and the residual are kept in FP64. Thus, with FP32 fields, each restart cycle is an iterative-refinement step, and we always test convergence on the FP64 residual 𝐛 − 𝐄𝑛+1∕2 . Since GMRES reads 𝐌 at every iteration, its storage type determines a large part of the memory traffic of the solve. The fourth level of Table 2 therefore stores 𝐌 in FP32 but computes it in FP64, which enables the FP64 matrix-engine deposition of Section 5.4, required by this kernel. The matrix-engine kernel accumulates the contributions of each cell in FP64 registers and rounds them only once, at the atomic flush.

5.6. Normalization of physical quantities In TPIC, physical quantities are expressed in dimensionless plasma units based on the Gaussian (CGS) system. Let 𝑛0 be the reference particle density, 𝑚𝑖 the reference ion mass, 𝑒 the elementary charge, and 𝑐phys the physical speed of light. The reference frequency (the ion plasma frequency 𝜔𝑝𝑖 ), the reference length (the ion skin depth 𝑑𝑖 ), and the reference electromagnetic-field strength 𝐹⋆ are √ 𝑐phys 𝑚𝑖 𝑐phys 𝜔𝑝𝑖 4𝜋𝑛0 𝑒2 𝜔𝑝𝑖 = , 𝑑𝑖 = , 𝐹⋆ = . (19) 𝑚𝑖 𝜔𝑝𝑖 𝑒 Using the superscript “phys” for dimensional quantities, the normalized variables are 𝐱phys , 𝑑𝑖

𝑡 = 𝜔𝑝𝑖 𝑡phys ,

𝐱=

𝐄phys 𝐄= , 𝐹⋆

𝐁phys 𝐁= . 𝐹⋆

𝐯=

𝐯phys , 𝑐phys

L. Pennati et al.: Preprint submitted to Elsevier

(20)

Page 11 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

Thus, the normalized speed of light is 𝑐 = 1. Species charge-to-mass ratios are expressed in units of 𝑒∕𝑚𝑖 . Singly charged ions √ have 𝑞∕𝑚 = +1, while electrons have 𝑞∕𝑚 = −𝜇, with 𝜇 = 𝑚𝑖 ∕𝑚𝑒 . At the reference density, 𝜔𝑝𝑒 ∕𝜔𝑝𝑖 = 𝜇 and 𝑑𝑒 ∕𝑑𝑖 = 𝜇−1∕2 . The normalized magnetic-field magnitude equals Ω𝑐𝑖 ∕𝜔𝑝𝑖 . Charge and current densities retain the Gaussian convention of the Maxwell equations: 𝜌=

𝜌phys , 4𝜋𝑒𝑛0

𝐉=

𝐉phys . 4𝜋𝑒𝑛0 𝑐phys

(21)

With these scalings, the Maxwell equations take the form ∇ ⋅ 𝐄 = 4𝜋𝜌,

𝜕𝐄 = ∇ × 𝐁 − 4𝜋𝐉, 𝜕𝑡

𝜕𝐁 = −∇ × 𝐄, 𝜕𝑡

∇ ⋅ 𝐁 = 0,

(22)

and the particle equation of motion 𝑞 𝑑𝐯 = (𝐄 + 𝐯 × 𝐁) . 𝑑𝑡 𝑚

(23)

We choose the factor 4𝜋 in Eq. 21 so the sources enter the field equations as 4𝜋𝜌 and 4𝜋𝐉, as in the Gaussian system. Let 𝑍𝑠 𝑒 be the charge of a particle of species 𝑠 and let 𝜌̂𝑠 =

|𝑍𝑠 | 𝑛𝑠 𝑛0

(24)

be the charge density of species 𝑠 in units of 𝑒 𝑛0 . For singly charged species 𝜌̂𝑠 = 𝑛𝑠 ∕𝑛0 . The normalized charge density is then 𝜌𝑠 = sign(𝑞∕𝑚)𝑠 𝜌̂𝑠 ∕(4𝜋), where sign(𝑞∕𝑚)𝑠 = ±1 is the sign of the charge of species 𝑠. The chargeto-mass ratio in units of 𝑒∕𝑚𝑖 is (𝑞∕𝑚)𝑠 = ±𝑍𝑠 𝑚𝑖 ∕𝑚𝑠 , and the plasma frequency of species 𝑠 in normalized units is 𝜔2𝑝𝑠 = 4𝜋𝜌𝑠 (𝑞∕𝑚)𝑠 = 𝜌̂𝑠 |𝑞∕𝑚|𝑠 . For the reference ions (𝜌̂ = 1, 𝑞∕𝑚 = 1) this gives 𝜔𝑝𝑖 = 1. Pressure and energy density are measured in units of 2 , 𝑃⋆ = 𝐹⋆2 = 4𝜋𝑛0 𝑚𝑖 𝑐phys

⋆ = 𝑃⋆ 𝑑𝑖3 ,

(25)

where ⋆ is the unit of integrated energy. With these scales, the field energy density is (|𝐄|2 + |𝐁|2 )∕(8𝜋) and the ∑ particle kinetic energy is 𝑝 𝑚𝑝 |𝐯𝑝 |2 ∕2, with the computational particle mass 𝑚𝑝 of Section 2 in units of 4𝜋𝑛0 𝑚𝑖 𝑑𝑖3 . √ Thermal velocities denote one-dimensional standard deviations: 𝑣th,𝑠 = 𝑘B 𝑇𝑠 ∕𝑚𝑠 ∕𝑐phys .

6. Experimental setup We run the experiments on a single node system equipped with an NVIDIA Grace Hopper 200 GPU, where the host is an ARM Neoverse-V2 72 cores CPU. The device has 96 GB of HBM, while the host has 586 GB of RAM. Large-scale simulations are performed on the accelerated partition of LUMI (LUMI-G), where each compute node has one 64-core AMD EPYC 7A53 CPU with 512 GB of DDR4 memory and four AMD MI250x accelerators. Each MI250x comprises two independent Graphics Compute Dies (GCDs), each with 64 GB of HBM. Thus, each node provides eight GCDs, and a total of 512 GB of GPU memory. The accelerators and CPU are interconnected through AMD Infinity Fabric, while inter-node communication is provided by four 200 Gbit/s HPE Cray Slingshot-11 network interfaces. Table 3 lists the software used to build and run TPIC on the two systems. Both builds use CMake Release mode. On LUMI-G, we trace MPI with CrayPat, GPU kernels with rocprofv3, and the full timeline with rocprof-sys. Kernel metrics come from Nsight Compute on GH200 and rocprof-compute on MI250x. If not reported otherwise, we run all experiments with TPIC partial sorting, matrix engines mass matrix assembly, and the single-position setup. Thus, all quantities are in FP64, except particle positions, which are in FP32. We use device-aware MPI as the communication backend.

L. Pennati et al.: Preprint submitted to Elsevier

Page 12 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Table 3 Software used to build and run TPIC. LUMI-G (MI250x)

nvc++ (NVHPC 26.5) nvcc, GCC 14.2.1 host CUDA 13.2 sm_90 Open MPI 5.0.10rc2 (HPC-X 2.50) – Nsight Compute

hipcc, AMD Clang 19.0.0 AMD Clang 19.0.0 ROCm/HIP 6.4.4 gfx90a Cray MPICH 9.0.1 2.22.3, aws-ofi-nccl 1.21.1 CrayPat 25.09, rocprofv3, rocprof-sys, rocprof-compute 100

|Etot(t)

− Etot(0)| / Etot(0)

10−4

10−2

10−6 10−8 10−10

TPIC, energy error explicit PIC, energy error TPIC, D/ x

10−12

λD = Δx

λ Δ

explicit PIC,

10−4

λD / Δx

C++ compiler Device compiler GPU toolkit Target MPI RCCL Profilers

GH200

λD/Δx

10−6

10−14 10−8

10−16 10−18

0

200

400

ωpet

600

800

1000

Figure 6: Finite-grid instability test. Relative error of the total energy (solid, left axis) and electron Debye length in units of the cell size (dashed, right axis) in TPIC and in the explicit PIC code. TPIC values are at round-off level; a zero energy error is drawn at 10−17 . The dotted line marks 𝜆𝐷𝑒 = Δ𝑥.

7. Numerical results 7.1. Finite grid instability removal We use a cold drifting plasma to verify that TPIC is free of the finite-grid instability (Lapenta, 2017). In explicit PIC, a plasma whose Debye length 𝜆𝐷𝑒 is not resolved by the grid is heated numerically by the aliasing of the particle density on the grid. We initialize cold electrons and ions of density 𝑛0 , with mass ratio 𝜇 = 2000, both drifting along 𝑥 with 𝑢0 = 0.2. The plasma carries no current and is in equilibrium, i.e., a plasma at rest seen from a moving frame, with no physical instability. We use a periodic domain of length 𝐿𝑥 = 128 𝑑𝑒 = 2.86 𝑑𝑖 , discretized with 128 × 1 × 1 cells, Δ𝑥 = 𝑑𝑒 = 2.24 × 10−2 𝑑𝑖 , so that the grid does not resolve 𝜆𝐷𝑒 = 0. We initialize 256 particles per cell per species, with electrons and ions at the same position, so that the initial charge density vanishes, and 𝐄 = 𝐁 = 0. We run the −1 −1 simulation for 40000 steps with 𝜃 = 1∕2 and Δ𝑡 = 5.59 × 10−4 𝜔−1 𝑝𝑖 = 0.025 𝜔𝑝𝑒 , up to 𝑇 = 1000 𝜔𝑝𝑒 . We compare TPIC with a 1D explicit, momentum-conserving electrostatic PIC code, with CIC shape functions, a second-order finite difference Poisson solver and a leapfrog mover, which starts from the same particles and uses the same Δ𝑥 and Δ𝑡. Fig. 6 shows the relative error of the total energy and the electron Debye length 𝜆𝐷𝑒 = 𝑣th,𝑒 ∕𝜔𝑝𝑒 , with 𝑣th,𝑒 the velocity spread about the drift. TPIC removes the finite-grid instability and keeps the initial equilibrium, with a GMRES solver tolerance of 1 × 10−16 , the total energy is conserved to machine precision, the electric energy stays at 10−31 of the initial electron kinetic energy, and the electrons remain cold. In the explicit code, the round-off perturbations of the charge density grow exponentially at 𝛾 ≃ 0.21 𝜔𝑝𝑒 up to 𝜔𝑝𝑒 𝑡 ≃ 180, and heat the plasma. The heating continues until the Debye length is comparable to the grid spacing, 𝜆𝐷𝑒 grows from zero and saturates at 𝜆𝐷𝑒 ≃ 0.36 Δ𝑥. Figs. 7 and 8 show the electron velocity distribution and phase space at 𝑡 = 0 and 𝑡 = 𝑇 . In TPIC, the electrons remain a cold beam at 𝑣𝑥 = 𝑢0 . In the explicit code, the distribution broadens, and the numerical heating accelerates electrons to unphysical superluminal velocities. Electrons exceed |𝑣𝑥 | = 𝑐 after 𝜔𝑝𝑒 𝑡 ≃ 500, and at 𝑇 , 1.2% of them are superluminal, with tails up to 𝑣𝑥 ≃ 3.8 𝑐.

L. Pennati et al.: Preprint submitted to Elsevier

Page 13 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

ωpet = 0 101

ωpet = 1000

TPIC explicit PIC

f(vx)

100 10−1 10−2 10−3

−4

−2

0

2

vx/c

−4

4

−2

0

2

vx/c

4

Figure 7: Finite-grid instability test. Electron velocity distribution in TPIC and in the explicit PIC code at 𝑡 = 0 (left) and 𝑡 = 𝑇 (right).

ωpet = 0

ωpet = 1000 explicit PIC

3

TPIC

2

vx/c

1 0

−1 −2 −3 0

20

40

60

80

x/de

100

120

0

20

40

60

80

100

120

x/de

Figure 8: Finite-grid instability test. Electron phase space (𝑥, 𝑣𝑥 ) of TPIC (dark blue) and the explicit PIC code (light purple), at 𝑡 = 0 (left) and 𝑡 = 𝑇 (right).

7.2. Two-stream instability We use the two-stream instability test to verify the electrostatic response of TPIC (Lapenta, 2017). In this experiment, we analyze the linear growth rate of a seeded mode and the exchange of energy between particles and the 𝐄 field. We initialize the system with two cold electron beams, each of density 𝑛0 ∕2, drifting along 𝑥 with velocities ±𝑢𝑏 , 𝑢𝑏 = 0.2, with a cold ion background of density 𝑛0 , with mass ratio 𝜇 = 2000. We use a periodic domain of length 𝐿𝑥 = 2𝜋 𝑑𝑒 = 0.1405 𝑑𝑖 , discretized with 128 × 1 × 1 cells, Δ = 1.10 × 10−3 𝑑𝑖 = 4.9 × 10−2 𝑑𝑒 . We initialize 256 particles per cell per species. We seed the mode 𝑚 = 3, which has 𝑘 = 2𝜋𝑚∕𝐿𝑥 = 3 𝑑𝑒−1 introducing a perturbation in particle densities 𝑛 = 𝑛0 (1 + 𝛿 cos(𝑘𝑥)), with 𝛿 = 10−4 . The electric field is initialized from Gauss’s −1 law, 𝐸𝑥 = −(𝛿∕𝑘) sin(𝑘𝑥). We run the simulation for 2000 steps with 𝜃 = 1∕2 and Δ𝑡 = 5.59 × 10−4 𝜔−1 𝑝𝑖 = 0.025 𝜔𝑝𝑒 , up to 𝑇 = 50 𝜔−1 𝑝𝑒 . For cold beams and immobile ions, a perturbation ∝ exp[𝑖(𝑘𝑥 − 𝜔𝑡)] follows the dispersion relation 1=

𝜔2𝑝𝑒 ∕2 (𝜔 − 𝑘𝑢𝑏 )2

+

𝜔2𝑝𝑒 ∕2 (𝜔 + 𝑘𝑢𝑏 )2

(26)

,

whose solutions are 𝜔2 ∕𝜔2𝑝𝑒 = 𝑎2 + 21 ± 12 growth rate 𝛾 = 𝜔𝑝𝑒

( √ 1 2

1 + 8𝑎2 − 𝑎2 − 21

)1∕2

√

.

1 + 8𝑎2 , with 𝑎 = 𝑘𝑢𝑏 ∕𝜔𝑝𝑒 . For 𝑎 < 1, we have a growing mode with

(27)

In the experiment, the seeded mode has 𝑎 = 0.6, thus 𝛾 = 0.3534 𝜔𝑝𝑒 . Fig. 9 (left) shows the amplitude |𝐸̂ 𝑥 (𝑚 = 3)| of Fourier transform of the seeded mode, which grows exponentially. A fit of ln |𝐸̂ 𝑥 | over 10 ≤ 𝜔𝑝𝑒 𝑡 ≤ 25 gives 𝛾 = 0.3539 𝜔𝑝𝑒 , 0.13% above the theoretical growth rate. L. Pennati et al.: Preprint submitted to Elsevier

Page 14 of 36

100 10 2 10 4 10 6 10 8 10 10 10 12

TPIC Cold theory

100 10 2

E(t)/Etot(0)

e |Ex(m = 3)|/(mic pi)

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

10 4 10 6 10 8 0

10

20

pet

30

40

50

Kinetic Electric Total

0

10

20

pet

30

40

50

Figure 9: Two-stream instability. Left: amplitude of the seeded Fourier mode 𝑚 = 3 of 𝐸𝑥 in TPIC and exponential growth at the cold rate 𝛾 of Eq. 27. Right: kinetic, electric and total energy, in units of the initial total energy.

ωpet = 0

0.50 0.25

/

vx c

ωpet = 30

= +0.2 c u0 = −0.2 c

ωpet = 50

u0

0.00 −0.25 −0.50 0

1

2

3

/

x de

4

5

6 0

1

2

3

/

x de

4

5

6 0

1

2

3

/

x de

4

5

6

Figure 10: Two-stream instability. Electron phase space (𝑥, 𝑣𝑥 ) of the two beams at 𝜔𝑝𝑒 𝑡 = 0, 30 and 50.

The growth saturates at 𝜔𝑝𝑒 𝑡 ≃ 30, due to the trapping of the beam electrons. Fig. 10 shows the electron phase space. At 𝜔𝑝𝑒 𝑡 = 30, the two beams reach 𝑣𝑥 = 0 and are trapped by the wave, and at 𝑇 the trapped electrons form phase-space vortices. Fig. 9 (right) shows the kinetic, electric, and total energy, normalized to the initial total system energy. The electric energy grows from 1.4 × 10−8 at the rate 2𝛾 and peaks at 36% of the total energy at 𝜔𝑝𝑒 𝑡 = 31.4, where the kinetic energy has dropped by the same amount. Then, the two exchange energy periodically. The magnetic energy remains below 10−28 of the total. With a GMRES solver tolerance of 1 × 10−16 , the total system energy is conserved to a relative error of 5.6 × 10−16 over the run, including the nonlinear phase.

7.3. Quasineutral plasma expansion We use the expansion of a quasineutral plasma into vacuum to verify TPIC solver in the electrostatic limit, in three dimensions and in a closed domain. In particular, this nonlinear test investigates the energy exchange between electrons and ions through the ambipolar electric field. We initialize electrons and ions with the same Gaussian density 𝑛 = 𝑛0 exp(−𝑟2 ∕2𝜎02 ), where 𝑟 is the distance from the domain center and 𝜎0 = 0.5 𝑑𝑖 . The electrons are Maxwellian with 𝑣th,𝑒 = 0.05, the ions are cold, and the mass ratio is 𝜇 = 100. Thus, Λ = 𝜎0 ∕𝜆𝐷𝑒0 = 100, with 𝜆𝐷𝑒0 = 𝑣th,𝑒 ∕𝜔𝑝𝑒 = 5 × 10−3 𝑑𝑖 the Debye length at the center. The domain of 14 × 14 × 14 𝑑𝑖 is discretized with 256 × 256 × 256 cells, Δ𝑥 = Δ𝑦 = Δ𝑧 = 5.47 × 10−2 𝑑𝑖 = 10.9 𝜆𝐷𝑒0 , with PEC walls and reflecting particle boundaries on all faces. We initialize both species at the same positions, 8×8×4 per cell within 𝑟 ≤ 4.3 𝜎0 , for 1.30×108 particles in total. Thus, the initial charge density vanishes, and we start from 𝐄 = 𝐁 = 0. We apply six passes of the filter of Eq. 14 with 𝛼 = 0.2, and run the simulation for 1200 steps with 𝜃 = 1∕2 and Δ𝑡 = 0.25 𝜔−1 𝑝𝑖 (𝜔𝑝𝑒 Δ𝑡 = 2.5, 𝑐Δ𝑡∕Δ𝑥 = 4.57), −1 up to 𝑇 = 300 𝜔𝑝𝑖 . We use 8 AMD MI250x GCDs of LUMI-G with a 2 × 2 × 2 domain decomposition.

L. Pennati et al.: Preprint submitted to Elsevier

Page 15 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

t/τ = 1.49

9 6 3 0 0

9 6 3 0 0

9 6 3 0 0

6

9

x [di]

6

14 0

9

14

3 y [di]

3

6

9

x [di]

6

14 0

9

14

3 y [di]

10−1 ni/ni(0, 0)

14

z [di]

14

3

100

t/τ = 2.99

14

z [di]

z [di]

t/τ = 0.00

10−2

3

6

9

x [di]

6

14 0

9

14

10−3

3 y [di]

10−4

Figure 11: Quasineutral plasma expansion. Ion density at 𝑡∕𝜏 = 0, 1.49, and 2.99, in units of the initial central density. The octant facing the viewer is cut away. The white circles on the cut faces mark the analytic 𝑟 = 𝜎 (solid) and 𝑟 = 2𝜎 (dashed) of Eq. 28.

For Λ ≫ 1, the expansion is quasineutral, and the Vlasov equation admits an exact self-similar solution (Dorozhkina and Semenov, 1998b,a). The density remains Gaussian, the flow is linear in 𝑟, and the electrons cool adiabatically: 𝑛(𝑟, 𝑡) = 𝑛0

𝜎03 𝜎3

2

2

𝑒−𝑟 ∕2𝜎 ,

𝑡 𝐮= 𝐫, 2 𝑡 + 𝜏2

( ) 𝑡2 , 𝜎 2 = 𝜎02 1 + 𝜏2 𝜎2 𝑇𝑒 = 0, 𝑇𝑒0 𝜎2

(28)

with 𝜏 = 𝜎0 ∕𝑐𝑠 , and 𝑐𝑠 = 𝑣th,𝑒 (1 + 𝜇)−1∕2 the initial ion sound speed, which includes the electron inertia. Here, 𝑐𝑠 = 4.98 × 10−3 and 𝜏 = 100.4 𝜔−1 𝑝𝑖 , so the run ends at 𝑡 = 2.99 𝜏. The ambipolar field is radial, so 𝐁 remains zero, and the electron thermal energy is converted into ion kinetic energy, 𝐾𝑖 ∕𝐸0 = 𝜇∕(1 + 𝜇) 𝑡2 ∕(𝑡2 + 𝜏 2 ) and 𝐾𝑒 = 𝐸0 − 𝐾𝑖 , with 𝐸0 the initial energy, which is only kinetic since 𝐄 = 𝐁 = 0. Fig. 11 shows the ion density at three times. The plasma remains spherical, the rms widths of the ion density along 𝑥, 𝑦, and 𝑧 differ by less than 0.06%. Fig. 12 (left) compares the width, averaged over the three axes, and the ion density at the center with Eq. 28. The width agrees within 0.22%, and the center density within 3.0%. Fig. 12 (right) shows the radial ion flow, averaged over spherical shells, which follows the linear profile of Eq. 28 up to 𝑟 = 3 𝜎. For 𝑡 > 𝜏∕4, the fitted flow slope is within 0.6% of the analytic one. Fig. 13 (left) shows the kinetic energy of the two species. At 𝑡 = 2.99 𝜏, the ions carry 88.7% of 𝐸0 , 0.3% below the analytic value. The electric field energy, neglected in the analytic solution, reaches at most 0.78% of 𝐸0 , at 𝑡 = 0.13 𝜏, and the magnetic energy at most 1.5 × 10−5 . Fig. 13 (right) shows the relative change of the total energy and the rms residual of Gauss’s law. With a GMRES solver tolerance of 1 × 10−15 , the total system energy is conserved to a relative error of 6.3 × 10−15 over the run, with PEC walls on all faces. Since ECSIM does not conserve charge exactly, the Gauss’s law residual is not zero, but it remains below 4 × 10−5 .

L. Pennati et al.: Preprint submitted to Elsevier

Page 16 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

0.6 2.0

0.4

1.5 1.0 0.0

0.5

1.0

1.5

2.0

2.5

t/τ (τ = 100.4)

3.0

t/τ = 1.25 t/τ = 1.99

2.0

ui, r/cs

0.8

t/τ = 0.25

2.5

ni(0, t)/ni(0, 0)

Centre density Analytic

2.5

σ(t)/σ0

1.0

Radius σ

3.0

t/τ = 2.99 Analytic

1.5 1.0

0.2

0.5

0.0

0.0 0.0

0.5

1.0

1.5

r/σ(t)

2.0

2.5

3.0

Figure 12: Quasineutral plasma expansion in TPIC (dotted lines) and Eq. 28 (wide light lines). Left: width of the ion density 𝜎, averaged over the three axes, and ion density at the center. Right: radial ion flow, averaged over spherical shells, in units of the initial ion sound speed 𝑐𝑠 , at four times.

1.0

0

1e−15

1e−5

0.4

−2 −4

0.2

Energy Gauss law

−6

0.0 0.0

0.5

1.0

1.5

2.0

t/τ (τ = 100.4)

2.5

3.0

0.0

0.5

1.0

1.5

2.0

/ (τ = 100.4)

tτ

4.0 3.5 3.0 2.5 2.0 1.5 1.0

rms(∇ ⋅ E − 4πρ)

Electrons Ions Analytic

0.6

(E − E0)/E0

Ks/E0

0.8

2.5

3.0

Figure 13: Quasineutral plasma expansion. Left: kinetic energy of electrons and ions in units of the initial energy 𝐸0 , in TPIC (dotted lines) and analytic (wide light lines). Right: relative change of the total energy (left axis) and rms residual of Gauss’s law (right axis).

7.4. Whistler waves Whistler waves are low frequency waves propagating in magnetized plasmas (Yu et al., 2018; Ballerini, Giulio et al., 2024). We use them to verify the dispersion of electromagnetic waves propagating along a uniform magnetic field, with a setup adapted from the whistler test of Gonzalez-Herrero et al. (2018) to a mass ratio 𝜇 = 1000. We use a periodic 2D domain of 20 𝑑𝑖 × 1.25 𝑑𝑖 (632 𝑑𝑒 × 39.5 𝑑𝑒 ), discretized with 1024 × 64 cells, Δ𝑥 = Δ𝑦 = 1.95 × 10−2 𝑑𝑖 = 0.62 𝑑𝑒 . The background field is 𝐁0 = (0.015, 0, 0), which gives |Ω𝑐𝑒 | = 15 𝜔𝑝𝑖 , Ω𝑐𝑖 = 0.015 𝜔𝑝𝑖 and 𝜔𝑝𝑒 ∕|Ω𝑐𝑒 | = 2.11. The electrons are bi-Maxwellian, with thermal velocities 𝑣th,𝑒∥ = 0.01 along 𝐁0 and 𝑣th,𝑒⟂ = 0.03 across it (𝑇𝑒⟂ ∕𝑇𝑒∥ = 9). The ions are Maxwellian with 𝑣th,𝑖 = 9.49 × 10−4 , i.e., 𝑇𝑖 = 𝑇𝑒⟂ . We initialize 256 particles per cell per species. We advance the simulation for 6000 steps with 𝜃 = 1∕2 and Δ𝑡 = 0.015 𝜔−1 𝑝𝑖 (|Ω𝑐𝑒 |Δ𝑡 = 0.225), up −1 −1 to 𝑇 = 90 𝜔𝑝𝑖 = 1350 |Ω𝑐𝑒 | . We do not seed any wave, the modes are excited by the particle noise. For 𝐤 ∥ 𝐁0 , the cold-plasma dispersion relation separates into a right-handed and a left-handed relation, with the handedness being the one of the field rotation about 𝐁0 : 𝜔2𝑝𝑖 𝜔2𝑝𝑒 𝑘2 𝑐 2 =𝑅=1− − , 𝜔 (𝜔 − |Ω𝑐𝑒 |) 𝜔 (𝜔 + Ω𝑐𝑒 ) 𝜔2 𝜔2𝑝𝑖 𝜔2𝑝𝑒 𝑘2 𝑐 2 =𝐿=1− − , 𝜔 (𝜔 + |Ω𝑐𝑒 |) 𝜔 (𝜔 − Ω𝑐𝑖 ) 𝜔2

(29) (30)

where Ω𝑐𝑒 (Ω𝑐𝑖 ) is the electron (ion) gyrofrequency and 𝜔𝑝𝑒 (𝜔𝑝𝑖 ) the electron (ion) plasma frequency. Each relation has two propagating solutions. Eq. 29 gives the whistler wave, with 0 < 𝜔 < |Ω𝑐𝑒 |, and a high-frequency electromagnetic L. Pennati et al.: Preprint submitted to Elsevier

Page 17 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers 3

Light wave (R)

10

Ion-cyclotron Light wave (L)

2

0

−2

̂ − iB̂ z|2

ω/Ωce

5

ω/Ωce

−1

Whistler

1

−3

0

−4

−1

−5

−5

−2

−10

−6

−3 −4

−2

0

k ∥ de

2

4

log10 |By

Cold relation

−2

−1

−7 0

k ∥ de

1

2

Figure 14: Whistler waves test. Power spectrum log10 |𝐵̂ 𝑦 − 𝑖𝐵̂ 𝑧 |2 of the transverse magnetic field. Positive (negative) frequencies are right-handed (left-handed) waves. Left panel: all resolved wavenumbers, with frequencies up to the Nyquist frequency 𝜋∕Δ𝑡. The box indicates the region enlarged on the right. The dashed lines represent the solutions of the cold dispersion relations Eq. 29 and Eq. 30, drawn for 𝑘 > 0 only.

mode named the R light wave. Eq. 30 gives the ion-cyclotron wave, with 0 < 𝜔 < Ω𝑐𝑖 , and the high-frequency L light wave. Fig. 14 shows the power spectrum |𝐵̂ 𝑦 − 𝑖𝐵̂ 𝑧 |2 , where the hat denotes the Fourier transform, in space 𝑥 of the field averaged over 𝑦, and in time, as shown in Matsumoto and Omura (1985). With this combination, positive (negative) frequency corresponds to a right-handed (left-handed) wave, and the sign of 𝜔∕𝑘 gives its direction of propagation. The spectral resolution is Δ𝑘 = 9.9 × 10−3 𝑑𝑒−1 and Δ𝜔 = 4.7 × 10−3 |Ω𝑐𝑒 |. The Nyquist limits are |𝑘|𝑑𝑒 = 5.09 and 𝜋∕Δ𝑡 = 13.96 |Ω𝑐𝑒 |. The TPIC whistler mode agrees with Eq. 29. Over the 200 modes with 0.5 ≤ |𝑘|𝑑𝑒 ≤ 1.5, the spectral peak lies, on average, 1.9% above the cold dispersion relation and, at most, 4.4%. This discrepancy is explained by the discretization and the finite electron temperature. To follow the propagation of the whistlers, we keep in the spectrum of 𝐵̂ 𝑦 − 𝑖𝐵̂ 𝑧 the 20 modes with 0.6 ≤ 𝑘𝑑𝑒 ≤ 0.8 and the frequencies within 30% of Eq. 29, and transform back to 𝑥 and 𝑡. This isolates the whistlers traveling towards +𝑥, which form into wave packets, as shown in Fig. 15, that travel at the group speed. The cold relation, Eq. 29 without its ion term, gives a group speed of 0.608 𝑑𝑒 |Ω𝑐𝑒 |, while the tracked packet, in orange in Fig. 15, moves at 0.613 𝑑𝑒 |Ω𝑐𝑒 |. The light waves have frequencies below the cold relation prediction, with the discrepancy increasing as |𝑘| grows. The difference is 1-6% for |𝑘|𝑑𝑒 < 1 and up to 46% at the largest wavenumber. This is due to the numerical dispersion of the ECSIM scheme, which, for 𝜃 = 1∕2, does not remove unresolved high-frequency modes, but it shifts their frequencies with the map 𝜔𝑛𝑢𝑚 = (2∕Δ𝑡) arctan(𝜔Δ𝑡∕2) < 𝜋∕Δ𝑡. For 𝜃 > 1∕2 the scheme is dissipative and damps these modes (Gonzalez-Herrero et al., 2018). The ion-cyclotron mode is not resolved in frequency. The run lasts 0.21 ion gyroperiods, so that Δ𝜔 = 4.7 Ω𝑐𝑖 , and the wave collapses onto the line 𝜔 = 0. Resolving it requires a run about two orders of magnitude longer. With a GMRES solver tolerance of 1 × 10−16 , the total system energy is conserved to a relative error of 9.7 × 10−15 over the run.

7.5. Weibel instability The Weibel instability generates magnetic fields from the temperature anisotropy of an unmagnetized plasma (Weibel, 1959). We use it to verify the kinetic electromagnetic response of TPIC. We initialize bi-Maxwellian electrons with thermal velocities 𝑣th,𝑒,𝑥 = 𝑣th,𝑒,𝑧 = 0.05 and 𝑣th,𝑒,𝑦 = 0.1, so that 𝑇𝑒,𝑦 ∕𝑇𝑒,𝑥 = 4, and cold ions of the same density, with mass ratio 𝜇 = 100. We use a periodic 2D domain of 4 𝑑𝑖 × 4 𝑑𝑖 (40 𝑑𝑒 × 40 𝑑𝑒 ), discretized with 448 × 448 cells, Δ𝑥 = Δ𝑦 = 8.93 × 10−3 𝑑𝑖 = 8.93 × 10−2 𝑑𝑒 . We initialize 256 particles per cell per species. We seed the mode 𝑚 = 4, with 𝑘 = 2𝜋𝑚∕𝐿𝑥 = 0.628 𝑑𝑒−1 , through the magnetic field 𝐵𝑧 = 𝛿𝐵 cos(𝑘𝑥), with 𝛿𝐵 = 10−5 and 𝐄 = 0. We run the −1 −1 simulation for 12000 steps with 𝜃 = 1∕2 and Δ𝑡 = 2.5 × 10−3 𝜔−1 𝑝𝑖 = 0.025 𝜔𝑝𝑒 , up to 𝑇 = 300 𝜔𝑝𝑒 .

L. Pennati et al.: Preprint submitted to Elsevier

Page 18 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Packet speed 0.613 deΩce

1.0

1200

Envelope of By / (10−4B0)

0.8 1000

0.6

Ωcet

800 600

0.4

400

0.2 200 0

0

200

400

600 0

x/de

200

400

600

0.0

x/de

Figure 15: Whistler waves test. Whistler wave packets with 0.6 ≤ 𝑘𝑑𝑒 ≤ 0.8 travelling towards +𝑥. Left: 𝐵𝑦 (light blue) and its envelope (dark blue) at seven different times. The orange line marks one packet. The dots lie on the line fitted through its peak. Right: envelope evolution over 𝑥 and 𝑡, the dashed line is the path of the tracked packet.

For 𝐤 along the cold direction 𝑥, the transverse mode with components 𝐸𝑦 and 𝐵𝑧 follows the kinetic dispersion relation (Davidson et al., 1972) 𝜔2 − 𝑘2 𝑐 2 − 𝜔2𝑝𝑒 − 𝜔2𝑝𝑖 + 𝜔2𝑝𝑒

𝑇𝑒,𝑦 𝑇𝑒,𝑥

[1 + 𝜉𝑍(𝜉)] = 0,

𝜔 𝜉=√ , 2 𝑘 𝑣th,𝑒,𝑥

(31)

where 𝑍 is the plasma dispersion function, ∞

2

1 𝑒−𝑠 𝑍(𝜉) = √ 𝑑𝑠, 𝜋 ∫−∞ 𝑠 − 𝜉

Im 𝜉 > 0.

(32)

𝑍 describes the electron response. The instability requires |𝜉| ≲ 1, where the thermal motion along 𝐤 couples the √ anisotropy to the mode. For a purely growing mode, 𝜔 = 𝑖𝛾, 𝜉 = 𝑖𝜂 with 𝜂 = 𝛾∕( 2 𝑘 𝑣th,𝑒,𝑥 ), and 𝑍(𝑖𝜂) = √ 2 𝑖 𝜋 𝑒𝜂 erf c(𝜂), so that Eq. 31 becomes [ ] √ 2 𝛾 2 + 𝑘2 𝑐 2 + 𝜔2𝑝𝑖 = 𝜔2𝑝𝑒 𝐴 − (𝐴 + 1) 𝜋 𝜂 𝑒𝜂 erf c(𝜂) , (33) with 𝐴 = 𝑇𝑒,𝑦 ∕𝑇𝑒,𝑥 − 1. Modes with 𝑘𝑑𝑒 < (𝐴 − 1∕𝜇)1∕2 = 1.73 are unstable. The seeded mode grows at 𝛾 = 0.0307 𝜔𝑝𝑒 , close to the fastest-growing mode. Fig. 16 (left) shows the amplitude of the seeded mode |𝐵̂ 𝑧 (𝑚 = 4)| averaged over 𝑦. The magnetic seed also excites a light wave at 𝜔 = 1.19 𝜔𝑝𝑒 , undamped for 𝜃 = 1∕2, which superposes with the mode and remains clearly visible up to 𝜔𝑝𝑒 𝑡 ≃ 60. The electron anisotropy relaxes during the run, so we evaluate Eq. 33 using temperatures derived from the electron pressure at each time. A fit of ln |𝐵̂ 𝑧 | over 60 ≤ 𝜔𝑝𝑒 𝑡 ≤ 160 gives 𝛾 = 0.0301 𝜔𝑝𝑒 , 0.9% above the theoretical rate averaged over the same interval. Fig. 16 (right) shows the instantaneous growth rate, measured over windows of −3 10 𝜔−1 𝑝𝑒 . It follows the theory within 2.3 × 10 𝜔𝑝𝑒 for 60 ≤ 𝜔𝑝𝑒 𝑡 ≤ 240, while 𝑇𝑒,𝑦 ∕𝑇𝑒,𝑥 decreases from 3.9 to 1.9. The mode saturates at 𝜔𝑝𝑒 𝑡 ≃ 260, with amplitude 306 𝛿𝐵 and a magnetic energy of 3.3% of the total energy. With a GMRES solver tolerance of 1 × 10−15 , the total system energy is conserved to a relative error of 5.4 × 10−14 over the run. Fig. 17 shows the evolution of 𝐵𝑧 averaged over 𝑦 in space and time. The seeded mode holds more than 80% of the power of the averaged field from 𝜔𝑝𝑒 𝑡 ≃ 80 and 99% at saturation. The mode grows without propagating, as expected for a purely unstable mode with 𝜔 = 𝑖𝛾.

L. Pennati et al.: Preprint submitted to Elsevier

Page 19 of 36

| ̂ z(

0.04

TPIC Kinetic theory

102

0.03

γ/ωpe

B m = 4)| / δB

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

101

100

0.02 TPIC Kinetic theory Kinetic theory, = 0

0.01

t

0

50

100

150

200

ωpet

0.00

250

0

50

100

150

ωpet

200

250

Figure 16: Weibel instability. Left: amplitude of the seeded mode 𝑚 = 4 of 𝐵𝑧 , in units of the seed 𝛿𝐵, and growth at the kinetic rate of Eq. 33 with the measured electron temperatures. Right: instantaneous growth rate of the seeded mode and Eq. 33 with the measured and with the initial temperatures.

250

1.5

325 δB

1.0

200

181 δB

150

62 δB

0.0

100

19 δB

−0.5

0.5

−1.0

6 δB

50

−1.5

2 δB

00

10

20

/

x de

30

40

̄ /rmsx B̄ z

ωpet

2.0

307 δB

Bz

300

0

10

20

/

x de

30

40

−2.0

Figure 17: Weibel instability: 𝐵𝑧 averaged over 𝑦, with the Fourier modes 𝑚 ≤ 16. Left: 𝐵𝑧 at seven different times, each profile rescaled to its own maximum |𝐵𝑧 |; the right axis gives this maximum in units of the seed 𝛿𝐵. Right: evolution over 𝑥 and 𝑡, with 𝐵𝑧 at each time divided by its root mean square over 𝑥. The dashed line follows a crest of the seeded mode, the dots mark it at the times of the left panel.

7.6. Whistler turbulence We run a three-dimensional whistler turbulence simulation to investigate the energy cascade from the seeded modes to small electron scales. We use a setup adapted from Gary et al. (2012); Saito et al. (2008). The √ plasma is homogeneous, with mass ratio 𝜇 = 1836 and a uniform magnetic field 𝐁0 = 𝐵0 𝐳̂ , with 𝜔𝑝𝑒 ∕|Ω𝑐𝑒 | = 5. Electrons and ions are Maxwellian, with 𝑇𝑖 = 𝑇𝑒 and 𝛽𝑒 = 𝛽𝑖 = 0.1, so that 𝑣th,𝑒 = 0.1 and 𝜆𝐷𝑒 = 0.1 𝑑𝑒 . The periodic cubic domain of side 𝐿 = 1.40 𝑑𝑖 = 60 𝑑𝑒 , is discretized with 300 × 300 × 300 cells, Δ𝑥 = Δ𝑦 = Δ𝑧 = 0.2 𝑑𝑒 = 2 𝜆𝐷𝑒 . We seed 150 helical whistler modes of equal magnetic energy, with random phases and the cold-plasma frequency, at the wavevectors 𝐤 = 2𝜋(𝑛𝑥 , 𝑛𝑦 , 𝑛𝑧 )∕𝐿 with 1 ≤ |𝑛𝑧 | ≤ 3 and 𝐤⟂ along 𝑥, 𝑦, or the diagonals, with |𝑛𝑥 |, |𝑛𝑦 | ≤ 3, so that 𝑘∥ 𝑑𝑒 ≤ 0.31 and 𝑘⟂ 𝑑𝑒 ≤ 0.44. The electric field follows from Faraday’s law and the electron flow from Ampère’s law, while the ions are at rest. The initial fluctuation energy is 𝜖0 = ⟨|𝛿𝐁|2 ⟩∕𝐵02 = 0.1. We initialize 256 particles per cell per species, 1.38 × 1010 particles in total. We apply one pass of the filter of Eq. 14 with 𝛼 = 0.5,√and run the simulation −1 for 4000 steps with 𝜃 = 1∕2 and Δ𝑡 = 5.83 × 10−3 𝜔−1 𝑝𝑖 = 0.25 𝜔𝑝𝑒 (|Ω𝑐𝑒 |Δ𝑡 = 0.11, 𝑐Δ𝑡 3∕Δ𝑥 = 2.17), up to −1 −1 −1 𝑇 = 23.3 𝜔𝑝𝑖 = 1000 𝜔𝑝𝑒 = 447 |Ω𝑐𝑒 | . We use single-precision particle positions and 32 AMD MI250x GCDs on 4 LUMI-G nodes, with a 4 × 4 × 2 domain decomposition. The 4000 steps take 1.6 h of wall-clock time, including the periodic output of fields, particles, and checkpoints, and the median time step takes 779 ms. Fig. 18 shows the electron density fluctuations, absent in the initial state, which develop during the cascade. Fig. 19 shows the reduced magnetic energy spectra. At 𝑡 = 0, the energy is only in the seeded modes. The nonlinear interaction of the whistlers transfers energy to the unseeded modes at shorter wavelengths, which hold up to 10% of the initial L. Pennati et al.: Preprint submitted to Elsevier

Page 20 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers ωpet = 0

ωpet = 300

ωpet = 1000 0.3

0.1 0.0

δne/n0 (arrow: B0)

0.2

−0.1 −0.2 −0.3

Figure 18: Whistler turbulence. Volume rendering of the electron density fluctuation 𝛿𝑛𝑒 ∕𝑛0 at 𝜔𝑝𝑒 𝑡 = 0, 300, and 1000. The arrow marks the direction of 𝐁0 .

10−1

10−2

10−2

10−3

10−3

S(k ∥ ) [B02de]

S(k ⟂ ) [B02de]

10−1

10−4 10−5 10−6 10−7 10−8

ωpet = 0 ωpet = 125 ωpet = 300 ωpet = 1000 largest seeded k

10−4 10−5 10−6

4.13 4.90 fit: k − , k− , kbde = 1.04 ⟂ ⟂

10−7

largest seeded k

0.3

0.5

1.0

2.0

k ⟂ de

4.0

8.0

10−8

0.3

0.5

1.0

2.0

4.0

8.0

k ∥ de

Figure 19: Whistler turbulence. Reduced magnetic energy spectra as a function of 𝑘⟂ (left) and 𝑘∥ (right) at four times. The dotted lines mark the largest seeded wavenumber, the dashed line is the broken power-law fit at 𝜔𝑝𝑒 𝑡 = 300.

fluctuation energy at 𝜔𝑝𝑒 𝑡 ≃ 125. The spectra extend to 𝑘𝑑𝑒 ≃ 3, where they reach the particle noise level. The cascade is anisotropic, with more energy at 𝑘⟂ than at 𝑘∥ . At 𝜔𝑝𝑒 𝑡 = 300 (|Ω𝑐𝑒 |𝑡 = 134), 𝑆(𝑘⟂ ) follows a broken power law, with slopes −4.1 and −4.9 below and above 𝑘⟂ 𝑑𝑒 = 1.04, close to the electron-scale break at 𝑘⟂ 𝑑𝑒 ≃ 1 found by Gary et al. (2012). Fig. 20 (left) shows that the electrons heat anisotropically. The parallel electron temperature, computed without the local bulk flow, grows to 1.96 𝑇0 and the perpendicular one to 1.25 𝑇0 , with 𝑇0 = 𝑚𝑒 𝑣2th,𝑒 , so that 𝑇∥ ∕𝑇⟂ = 1.57 at 𝑇 , close to the 1.6 found by Gary et al. (2012). Fig. 20 (right) details the energy exchange. At 𝑇 , the magnetic fluctuation energy has decayed to 31% of its initial value. The electrons absorb 64% of the initial fluctuation energy and the ions 5%, while the electric energy is 1.3%. With a GMRES solver tolerance of 1×10−15 , the total system energy is conserved to a relative error of 1.3 × 10−15 over the run.

7.7. Magnetic reconnection We perform a three-dimensional magnetic reconnection simulation with TPIC, following the GEM challenge setup (Birn et al., 2001), and we compare the results against the well-established iPIC3D (Markidis et al., 2010) code. The system is a single Harris current sheet of thickness 𝛿 = 0.5 𝑑𝑖 , where the magnetic field is 𝐵𝑥 = 𝐵0 tanh[(𝑦 − 𝐿𝑦 ∕2)∕𝛿], with 𝐵0 = 0.0195, so that Ω𝑐𝑖 = 0.0195 𝜔𝑝𝑖 and 𝑣𝐴 = 0.0195 𝑐. The plasma in the current sheet has density 𝑛0 sech2 [(𝑦−𝐿𝑦 ∕2)∕𝛿], while the uniform background plasma has density 0.1 𝑛0 . We use Maxwellian electrons and ions with mass ratio 𝜇 = 64 and thermal velocities 𝑣th,𝑒 = 0.045 and 𝑣th,𝑖 = 0.0126, so that 𝑇𝑖 ∕𝑇𝑒 = 5. A perturbation of the magnetic flux function, localized at the center of the sheet with width 𝛿 and uniform along 𝑧, seeds the reconnection at the X-line 𝑥 = 𝐿𝑥 ∕2, 𝑦 = 𝐿𝑦 ∕2. The drift of the sheet species along 𝑧, computed from L. Pennati et al.: Preprint submitted to Elsevier

Page 21 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers 2.0

T ∥ /T0 T ⟂ /T0

1.8

10−2

1.5

T ∥ /T ⟂

10−4

T ∥ /T ⟂

Te/T0

1.6 1.3 1.4

1.2

Energy/E0

1.4

10−6

UδB UE

10−8

Ke Ki

10−10

|Etot

− E0|

10−12

1.2

1.1

1.0

1.0 0

200

400

ωpet

600

800

1000

10−14 10−16 0

200

400

ωpet

600

800

1000

Figure 20: Whistler turbulence. Left: parallel and perpendicular electron temperatures in units of the initial temperature 𝑇0 , and their ratio (right axis). Right: magnetic fluctuation energy 𝑈𝛿𝐵 , electric energy 𝑈𝐸 , kinetic energy of electrons 𝐾𝑒 and ions 𝐾𝑖 , and change of the total energy, in units of the initial total energy 𝐸0 .

∇ × 𝐁, carries the current of the initial field. The domain of 30 × 20 × 20 𝑑𝑖 is discretized with 300 × 200 × 200 cells, Δ𝑥 = Δ𝑦 = Δ𝑧 = 0.1 𝑑𝑖 = 0.8 𝑑𝑒 . It is periodic in 𝑥 and 𝑧, with PEC walls and reflecting particle boundaries along 𝑦. 10 We initialize 256 particles per cell for each of the four species (sheet and background √ particles), 1.23 × 10 particles −3 in total. We use Δ𝑡 = 0.1 𝜔−1 𝑝𝑖 (Ω𝑐𝑖 Δ𝑡 = 1.95 × 10 , 𝜔𝑝𝑒 Δ𝑡 = 0.8), so that 𝑐Δ𝑡 3∕Δ𝑥 = 1.73 exceeds the 3D CFL limit of explicit PIC. Both TPIC and iPIC3D apply six passes of the filter of Eq. 14 with 𝛼 = 0.2, and use a GMRES tolerance of 10−15 . TPIC uses 𝜃 = 1∕2 and single-precision particle positions, while iPIC3D, which implements the implicit moment method (Brackbill and Forslund, 1982), uses 𝜃 = 1. We run both codes on 32 AMD MI250x GCDs on 4 LUMI-G nodes, with a 4 × 4 × 2 domain decomposition for 8000 steps, up to Ω𝑐𝑖 𝑡 = 15.6. We measure the reconnected flux along the X-line, [ ] 𝐿𝑥 𝐿𝑥 ∕2 1 𝜓(𝑧, 𝑡) = 𝐵 𝑑𝑥 − 𝐵𝑦 𝑑𝑥 , (34) ∫0 2 ∫𝐿𝑥 ∕2 𝑦 𝑦=𝐿𝑦 ∕2

and compute the reconnection rate 𝜓∕(𝐵 ̇ 0 𝑣𝐴 ) from field snapshots taken every 200 steps. In Fig. 21, we report its mean and standard deviation along the X-line in 𝑧, comparing TPIC and iPIC3D results. The left panel reports the reconnection rate. In both codes the rate is initially negative, as the initial perturbation decays and part of the seeded flux is lost. The rate becomes positive at Ω𝑐𝑖 𝑡 ≃ 1.4, remains at 0.02-0.03 for 2 ≲ Ω𝑐𝑖 𝑡 ≲ 4.7, and then increases to a first maximum, 0.047 at Ω𝑐𝑖 𝑡 = 5.85 in TPIC and 0.067 at Ω𝑐𝑖 𝑡 = 6.63 in iPIC3D. From Ω𝑐𝑖 𝑡 ≃ 9, the reconnection accelerates in both codes. Over the whole run, the mean rates differ by at most 0.059, at Ω𝑐𝑖 𝑡 = 12.09, comparable to the standard deviation along the X-line. The rate peaks at 0.366 at Ω𝑐𝑖 𝑡 = 14.04 in TPIC and at 0.383 at Ω𝑐𝑖 𝑡 = 13.26 in iPIC3D, and it decreases to 0.33 in both codes at Ω𝑐𝑖 𝑡 = 15.21. The standard deviation along the X-line is below 0.01 up to Ω𝑐𝑖 𝑡 ≃ 3.5. Then it grows, as the lower-hybrid drift instability develops on the flanks of the current sheet and the reconnection becomes non-uniform along the X-line, and at the end, when the drift kink instability develops, it is comparable to the mean rate. The right panel in Fig. 21 shows the relative change of the total energy. With a GMRES solver tolerance of 1×10−15 , TPIC conserves the total energy to a relative error of 1.7×10−14 over the 8000 steps. The iPIC3D scheme is dissipative, it loses 5.1 × 10−4 of the initial energy within the first 50 steps, mainly through the damping of the seeded perturbation. The iPIC3D energy loss remains nearly constant up to Ω𝑐𝑖 𝑡 ≃ 3, and then it increases as instabilities develop to 1.9 × 10−3 at Ω𝑐𝑖 𝑡 = 15.5. The lower-hybrid drift instability (LHDI) develops first, driven by the diamagnetic drift on the flanks of the current sheet, where the density gradient is large. We observe it through fluctuations in 𝐸𝑧 along 𝑧, as shown in Fig. 22. On both flanks of the sheet, the wave fronts extend along 𝑥 over the whole box, and 𝐸𝑧 alternates in sign along 𝑧 with wavelength 0.71 𝑑𝑖 . Flux ropes form in the outflow regions from Ω𝑐𝑖 𝑡 ≃ 5.5, while the reconnection rate increases. They are secondary magnetic islands in the 𝑥-𝑦 planes that extend along 𝑧 (Markidis et al., 2013). On longer timescales, the current sheet L. Pennati et al.: Preprint submitted to Elsevier

Page 22 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

0.8

(E − E0)/E0

0.6

ψ̇/(B0vA)

0.0000

TPIC iPIC3D

0.4 0.2

−0.0010 −0.0015

0.0

−0.2

−0.0005

0.0

2.5

5.0

7.5

Ωcit

10.0

12.5

15.0

TPIC iPIC3D

0.0

2.5

5.0

7.5

Ωcit

10.0 12.5 15.0

Figure 21: Magnetic reconnection in TPIC and iPIC3D. Left: reconnection rate 𝜓∕(𝐵 ̇ 0 𝑣𝐴 ), with 𝜓 given by Eq. 34. The lines show the mean along the X-line (averaged along 𝑧), and the shaded bands one standard deviation along it. Right: relative change of the total energy.

rho_i 1.2

Ez/(B0 vA) 0.94

1 0.8 0.6

0

0.4 0.2 0

-0.94

Figure 22: Lower-hybrid drift instability in TPIC at Ω𝑐𝑖 𝑡 = 4.68. Isosurfaces of positive (red) and negative (blue) 𝐸𝑧 ∕(𝐵0 𝑣𝐴 ), and volume rendering of the density of the current-sheet ions, in units of 𝑛0 , for 8 ≤ 𝑦∕𝑑𝑖 ≤ 12.

kinks along 𝑧 due to the onset of the drift kink instability (Daughton, 1998; Lapenta and Brackbill, 2002). Fig. 23 shows the flux ropes on both sides of the X-line at Ω𝑐𝑖 𝑡 = 15.21, and the kink of the rope at 𝑥 ≃ 𝐿𝑥 ∕4. Fig. 24 shows the density of the sheet ions and the magnetic field lines at the end of the TPIC run, at Ω𝑐𝑖 𝑡 = 15.6. The field lines reconnect across the sheet at the X-line, and the sheet ions accumulate in the outflow regions. The current sheet is no longer uniform along 𝑧, since the drift kink instability has developed. The figure also shows the most energetic electrons, with kinetic energy 𝐸 = 𝑚𝑒 |𝐯|2 ∕2 ≥ 10.6 𝑇0 , where 𝑇0 = 𝑚𝑒 𝑣2th,𝑒 is the initial temperature of the sheet electrons, corresponding to a speed of |𝐯| ≥ 0.207 𝑐 = 4.61 𝑣th,𝑒 = 10.6 𝑣𝐴 . TPIC runs this simulation in 2.2 hours, with a median time step, including fields and particle output, of 995 ms. On the same 32 GCDs, the median iPIC3D time step is 1330 ms 1.33× slower. Part of the difference comes from the field solve. The implicit moment method of iPIC3D yields a much worse-conditioned Maxwell system, and its GMRES solver requires on average 147 iterations per step to converge to the same 1 × 10−15 tolerance, against 18 in TPIC.

L. Pennati et al.: Preprint submitted to Elsevier

Page 23 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

rho_i 1.2 1 0.75 0.5 0.25 0

Figure 23: Flux ropes and kink of the current sheet in TPIC at Ω𝑐𝑖 𝑡 = 15.21. Volume rendering of the density of the current-sheet ions, in units of 𝑛0 , and magnetic field lines on both sides of the X-line, for 8 ≤ 𝑦∕𝑑𝑖 ≤ 12.

Figure 24: Magnetic reconnection in TPIC at Ω𝑐𝑖 𝑡 = 15.6. Volume rendering of the density of the current-sheet ions, in units of 𝑛0 , reconnected magnetic field lines, and a random sample of 20000 electrons with kinetic energy above the 98th percentile, 𝐸 ≥ 10.6 𝑇0 , colored by 𝐸∕𝑇0 . The ion density is shaped by the drift kink instability.

8. Performance analysis 8.1. Strong & weak scaling We run strong and weak scaling tests using the same GEM setup described in Section 7.7 without smoothing, with 𝜃 = 1∕2, and Δ𝑡 = 0.01 𝜔−1 𝑝𝑖 . We use full double precision, partial sorting, and the matrix-engine mass-matrix deposition. We use device-aware MPI as the communication backend, and we map one rank per GPU/GCD. Table 4 lists the baseline of the scaling tests. Weak scaling keeps 643 cells per GCD, up to 10243 cells and 1.9 × 1012 particles on 4,096 GCDs, while strong scaling keeps the initial global grid fixed and increases only the computational resources. We measure the execution time for 100 time steps, excluding start-up and finalization time, and repeat each test three times. The parallel efficiency is 𝑁0 𝑇𝑁0 ∕(𝑁 𝑇𝑁 ) in strong and 𝑇𝑁0 ∕𝑇𝑁 in weak scaling tests, where 𝑇𝑁 is the time on 𝑁 GCDs and 𝑁0 the baseline. Fig. 25 (left) shows the weak scaling results. At constant work per GCD, the execution time per time step increases from 650 ms at 32 GCDs to 671 ms on 4,096 GCDs, i.e., an efficiency of 96.8%. This 3% increase is due to communication. L. Pennati et al.: Preprint submitted to Elsevier

Page 24 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

Cells

ppc

GMRES tol

Subdomains

Weak Strong

256 × 256 × 128 256 × 256 × 256

4 × 448 4 × 448

10−5 10−5

4 × 4 × 2 = 32 4 × 4 × 4 = 64

40 30 20 10 0 32

60 40 Execution time 20 Scaling efficiency 64 128 256 512 1024 2048 4096 0

Number of GPUs/GCDs

50100.0% 98.7% 97.8% 95.1% 20

100 89.1%

80.6%

Parallel Efficiency (%)

70 100 60100.0% 99.4% 99.0% 98.5% 98.1% 97.6% 97.2% 96.8% 80 50

Execution time (s)

Test

Parallel Efficiency (%)

Execution time (s)

Table 4 Baseline of the scaling tests on LUMI-G. Particles per cell (ppc) are given as species × particles per species.

80

10

66.6% 60

5

40

2 64

Execution time Scaling efficiency 128 256 512

20 1024 2048 4096 0

Number of GPUs/GCDs

Figure 25: Weak and strong scaling on LUMI-G, where each node is equipped with 4× AMD MI250x GPUs, for a total of 8 GCDs per node. Wall-clock time of 100 steps (left axis) and parallel efficiency (right axis). Left: weak scaling from 32 to 4,096 GCDs. Right: strong scaling from 64 to 4,096 GCDs, with the time axis in logarithmic scale. Plots report the average time and efficiency measured over three runs per point. Standard deviation is always less than 0.5%.

In strong scaling, reported in Fig. 25 (right), the time per step decreases from 652 ms on 64 GCDs to 15.3 ms on 4,096 GCDs, a 42.7× speedup on 64× more GCDs. The efficiency remains above 95% up to 512 GCDs, with 323 cells per GCD, and drops to 89.1%, 80.6%, and 66.6% at 1024, 2048, and 4,096 GCDs. At 4,096 GCDs, each subdomain holds only 163 cells, with four times the surface-to-volume ratio of the baseline, making the communication cost from halo exchanges, particle migration, and GMRES global reductions more significant.

8.2. Mass-matrix deposition on tensor-cores We compare the matrix-engine mass matrix deposition kernel of Section 5.4 with two FP64 baselines, the unsorted per-particle kernel with global atomics, and a sorted kernel with vector FMAs. We test varying the number of particles per cell from 64 to 1024, with a fixed 643 cell domain. The two sorted kernels thus execute the same FLOPs and move the same data. We do not include the sorting cost here. Execution time and DRAM traffic are measured on each GPU (Nsight Compute on GH200, rocprof-compute on one MI250x GCD). We always count FLOPs with Nsight Compute on GH200, using only the FP64 vector operations of active threads plus tensor-core operations. We use the same counts for MI250x, whose counters report wavefront instructions and not per-thread operations. FP64 atomic additions in the main memory are treated as data transfer and are not counted as FLOPs. Fig. 26 show the rooflines for GH200 and MI250x, Table 5 details the execution times. The unsorted kernel has an arithmetic intensity (AI) of 0.075 FLOP/B on GH200 and 0.02 FLOP/B on MI250x, more than two orders of magnitude below the ridge points (8.3 and 16.6 FLOP/B for the FP64 and tensor-core (TC) roofs on GH200, 14.5 and 29.3 FLOP/B for the FP64 and MFMA roofs on MI250x). At the same FLOPs, the lower AI on MI250x is due to higher DRAM traffic, consistent with the smaller L2 cache of a GCD (8 MB vs. 60 MB). The kernel reaches a throughput of 81 and 6 GFLOP/s on the two architectures, respectively. It is memory-bound, but it sustains only 27% and 19% of peak bandwidth because atomic throughput limits it. Sorting reduces the DRAM traffic per particle by 92-163× on GH200 and 216-548× on MI250x, to about 40 B on both GPUs at 1024 ppc, as the per-cell traffic is amortized with more particles. The full 8 × 8 tiles increase the FLOPs per particle by 2.9×, and the AI grows by more than two orders of magnitude, to 20-35 FLOP/B on GH200 and 15-35 FLOP/B on MI250x. Thus, the sorted kernels move from the memory-bound to the compute-bound region. Without matrix engines, the sorted kernel uses less than 9% of the DRAM bandwidth and is 29-34× (GH200) and 77-82× (MI250x) faster than the unsorted kernel. L. Pennati et al.: Preprint submitted to Elsevier

Page 25 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Table 5 Execution time (ms) of the mass-matrix deposition kernels for different numbers of particles per cell (ppc): unsorted (𝑡U ), sorted without (𝑡S ) and with (𝑡M ) matrix engines, and speedups. Profiled kernel durations, averaged over two runs.

104

Speedup

Device

ppc

𝑡U

𝑡S

𝑡M

𝑡U ∕𝑡S

𝑡S ∕𝑡M

𝑡U ∕𝑡M

GH200

64 256 512 1024

98.5 390.8 789.7 1587.7

3.40 11.56 24.20 48.06

1.26 4.21 8.09 15.98

29.0 33.8 32.6 33.0

2.69 2.75 2.99 3.01

78.0 92.9 97.6 99.3

MI250x (1 GCD)

64 256 512 1024

1278.4 5252.7 10446.7 21306.6

16.59 65.32 130.29 261.18

5.49 17.61 34.86 69.72

77.0 80.4 80.2 81.6

3.02 3.71 3.74 3.75

232.7 298.4 299.7 305.6

FP64 roof: 33.45 TFLOP/s TC FP64 roof: 66.91 TFLOP/s Unsorted Sorted w/o TC Sorted TC

103

102 0.02 0.05 0.1 0.2

0.5 1

2

5 10 20

Arithmetic Intensity (FLOP/byte)

PPC 64 256 512 1024 50 100

Throughput (GFLOP/s)

Throughput (GFLOP/s)

Time (ms)

104

FP64 roof: 20.31 TFLOP/s MFMA FP64 roof: 41.06 TFLOP/s Unsorted Sorted w/o MFMA Sorted MFMA

103 102 101 0.01 0.03

0.1

0.3

1

3

10

Arithmetic Intensity (FLOP/byte)

PPC 64 256 512 1024 30 100

Figure 26: Roofline analysis of the mass-matrix deposition kernels, including unsorted deposition, sorted w/o matrix engines, and sorted with matrix engines. Left panel: roofline on NVIDIA GH200 GPU; right panel: roofline on AMD MI250x single GCD. Kernels are executed in isolation on a 643 -cell domain, varying the number of particles per cell from 64 to 1024. Arithmetic intensity and throughput of the tensor core kernel combine FLOPs from both conventional and TC operations, counted with Nsight Compute and rocprof-compute.

With the same AI, matrix engines only raise the throughput, executing 83% of the FLOPs. On GH200, throughput increases from 6.8-8.0 to 18.4-23.2 TFLOP/s (up to 27-35% of the TC peak), a 2.7-3.0× speedup. On MI250x, it increases from 1.4 to 4.2-5.3 TFLOP/s (10-13% of the MFMA peak), a 3.0-3.8× speedup. The gain grows with ppc as the per-cell work outside the matrix engines is amortized. Overall, sorting and matrix engines accelerate deposition by 78-99× on GH200 and 233-306× on MI250x, with the latter suffering particularly from atomic contention.

8.3. Trace & profiling We profile the 32-GCD weak scaling baseline (Table 4) in full double precision, with GMRES tolerance 10−15 , in three separate runs on LUMI-G to obtain a time breakdown and detailed information on computation and communication patterns. We use rocprof-sys to record a 10 steps trace of every rank. We collect ROCTX ranges marking each TPIC worker, GPU kernels, and MPI calls. rocprofv3 records every kernel dispatch and memory copy. CrayPat traces every MPI call over 50 steps and reports calls, messages, bytes, and time per rank and thread. We collect only data of the simulation main loop, thus excluding simulation startup and finalize. The three runs take about 647 ms per step, similar to the 650 ms measured in the weak scaling test. The trace in Fig. 27 shows one time step of rank 0. Each row shows a host worker of Section 5.1 submitting, polling, or waiting, while its GPU work runs asynchronously on the device streams. Worker 𝑀 enqueues the four particle pushes. Worker 𝐶 routes the migrating particles while 𝐿0 and 𝐿3 sort their resident particles (TPIC:: ResidentReorder). We omit workers 𝐿1 and 𝐿2 from the pictures for clarity. Each 𝐿𝑠 then waits for 𝐷 to deposit L. Pennati et al.: Preprint submitted to Elsevier

Page 26 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

Figure 27: rocprof-sys trace of one time step on rank 0 of the 32-GCD weak scaling baseline, on four LUMI-G nodes. Each row is a host worker with its ROCTX ranges. The local workers 𝐿1 and 𝐿2 are omitted. The top row collects the MPI calls of all threads, except the polling calls (MPI_Testall), which rocprof-sys does not record. Boxes mark the particle migration MPI messages and the moment and field solver MPI exchanges. The picture does not report GPU kernels. Given the asynchronous pipeline, host worker times do not reflect the actual GPU compute time.

Table 6 GPU utilization time by kernel per time step on rank 0 of the 32-GCD weak scaling baseline (rocprofv3). Kernels of different streams run concurrently, so their raw durations overlap. Each stage is the duration of its kernels when run alone on the GPU, multiplied by the launches per step. The idle time is the step time minus the GPU busy time. Sorting and deposition differ by less than 2% between ranks, while the idle time ranges from 25 to 41 ms. Stage

Time [ms]

Share

Particle sorting (4 species) Mass-matrix deposition (matrix engines) Particle mover Current deposition Field solver, migration kernels, copies GPU idle

324 131 99 36 20 37

50.1% 20.2% 15.3% 5.6% 3.1% 5.7%

Time step

647

100%

its sorted particles before appending the arrivals (TPIC::MigrationFinalize). Worker 𝐷 stays active for the whole step. It deposits each species once it is sorted and then assembles 𝐉̂ and 𝐌, after which 𝑆 runs the field solve. In the MPI row, the migration messages overlap the particle sorting, and they end before the first species is fully sorted. The moment halos and the GMRES exchanges close the step. Table 6 breaks down one time step on rank 0 from the rocprofv3 run. The GPU is busy 95% of the step on every rank, so the kernels drive the step time. Particle sorting takes half of the step, followed by the mass-matrix deposition (20%) and the mover (15%). Each rank sends 35 MB per step to its neighbors in about 320 point-to-point messages, and 69% of these bytes belong to the 𝐉̂ and 𝐌 halos, as detailed in Table 7. Only two exchanges lie on the critical path, the moment assembly and the field solve, which run between the end of the deposition and the next particle push. The GPU idles for 31 ms per step (median over ranks and steps), almost all of it in this window. This idling mainly comes from moment halo communication and GMRES reductions, whose completion depends on the slowest rank. Particle migration, instead, is fully hidden. Worker 𝐶 routes all species for 138 ms per step while the movers and sorts keep the GPU busy. The GPU idles only 0.05 ms in this interval, and routing completes at least 190 ms before the first species finishes sorting and needs its arrivals.

L. Pennati et al.: Preprint submitted to Elsevier

Page 27 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Table 7 MPI traffic per time step of a rank with 17 neighbor ranks (CrayPat), by exchange and worker. Exchange

Worker

MB

Messages

Critical path

𝐉̂ and 𝐌 halos GMRES halos and reductions Particle migration 𝐁 halo

𝐷 𝑆 𝐶, 𝐿𝑠 𝐹

26.6 8.7 2.7 0.6

49 251 21 + 4 × 17 17

yes yes no no

Execution time (s)

67.0

MPI RCCL

66.5 66.0 65.5 65.0 64.5

32

64

128

256

Number of GPUs/GCDs

512

Figure 28: Weak scaling execution time of 100 steps on LUMI-G with device-aware MPI and RCCL communication. Error bars show the standard deviation over runs.

8.4. Communication framework comparison We compare the RCCL and device-aware MPI communication configurations of TPIC on the weak scaling test of Table 4, from 32 to 512 GCDs. With RCCL, field halos, mass-matrix exchanges, particle payloads, and all reductions use RCCL. Fig. 28 shows the execution time, averaged over three runs per point. The two backends differ by less than 1% at all scales, with RCCL slightly faster up to 64 GCDs, and MPI better at larger scales. Differences for all points are within the statistical standard deviation. From 32 to 512 GCDs, the time increases by about 2% with MPI and 2.8% with RCCL, i.e., a weak scaling efficiency of 98% and 97.3%.

9. Discussion The ECSIM discretization introduces a spurious term in the Poynting theorem (Eq. 12) that hinders exact energy conservation in the case of a closed domain. This is canceled only by combining the zeroed electric components with the ghost-cell parities prescribed by the image field theory in Table 1. With this closure, the discrete Poynting flux vanishes for any combination of periodic, PEC, and PMC boundaries, edges, and corners included. All experiments reported in Section 7 confirm that TPIC conserves the total energy to machine precision, including three-dimensional runs in a closed domain with PEC walls. Exact energy conservation also removes the finite grid instability, in a cold drifting plasma with an unresolved Debye length, TPIC keeps the initial equilibrium, while the explicit PIC algorithm heats the plasma numerically. The two-stream instability, quasineutral plasma expansion, whistler waves, and Weibel instability tests agree with the corresponding analytical linear and non-linear theories, the whistler turbulence simulation reproduces the electron-scale spectral break and the anisotropic electron heating of Gary et al. (2012), and the magnetic reconnection simulation is consistent with the iPIC3D results. TPIC obtains these results with grid spacings larger than the Debye length in all tests and, in the three-dimensional ones, with time steps beyond the explicit CFL limit, while conserving the total energy to machine precision. This reduces the computational cost and keeps long runs stable, without numerical heating or damping. Although ECSIM does not exactly preserve the charge continuity equation, yielding an error in Gauss’ law, the experiments show that the error stays bounded and does not significantly affect result accuracy. Matrix engines combined with particle sorting greatly accelerate the mass-matrix deposition. The unsorted kernel is limited by atomic throughput at an arithmetic intensity below 0.1 FLOP/B. Sorting and in-register accumulation within each cell turn the scatter into dense 8 × 8 cell-local products, reduce the DRAM traffic per particle by two orders L. Pennati et al.: Preprint submitted to Elsevier

Page 28 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

of magnitude, and move the kernel into the compute-bound region. FP64 matrix engines then execute these products 2.7-3.0× faster on GH200 and 3.0-3.8× faster on MI250x than vector FMAs, for an overall speedup of up to 99× and 306×, respectively. The same matrix engine formulation maps onto the different instructions of the two vendors, one 8 × 8 × 4 MMA on NVIDIA and the four 4 × 4 × 4 blocks of one MFMA on AMD. The measured throughput is below the matrix-engine peaks due to the work executed outside the matrix engines, i.e., the evaluation of 𝜶𝑝 and the CIC weights, the operand shuffles, and the final atomic flush. The TPIC asynchronous pipeline executes particle migration, moment assembly, diagnostics, and output as independent tasks, hiding communication and yielding quasi-optimal scaling, as shown by LUMI-G results. In weak scaling, the time per step increases by only 3% from 32 to 4,096 GCDs, where TPIC pushes 1.9 × 1012 particles on a 10243 grid. In strong scaling, the efficiency stays above 95% up to 512 GCDs, and decreases to 66.6% at 4,096 GCDs, where each subdomain holds 163 cells and a time step takes 15.3 ms. In small subdomains, the calculation-tocommunication ratio is less favorable, and latency due to halo exchanges and GMRES global reductions hinders the performance. Profiling on 32 GCDs shows that TPIC is bound by GPU kernels. The GPU is busy 95% of the time step, and particle migration is fully hidden. Particle sorting takes half of the step, resulting in the most time consuming operation of the loop. Sorting nevertheless pays off by a large margin, since it enables cell-based mover, current, and mass-matrix kernels that reuse cached cell fields, access memory coalesced, and accumulate without global atomics (Pennati and Markidis, 2026). On one MI250x GCD, the unsorted mass-matrix kernel alone takes 10.4 s per species at 512 ppc (Table 5). TPIC instead sorts one species with 448 ppc in 81 ms and deposits its mass matrix with matrix engines in 33 ms. Neither MPI nor vendor-specific communication backends are uniformly faster. RCCL and device-aware MPI differ by less than 1% up to 512 GCDs, with RCCL ahead at small scale and MPI at larger scale (Fig. 28). TPIC’s communication patterns have different requirements, from bandwidth-bound halo, mass-matrix, and particle exchanges to small, latency-bound GMRES reductions on the critical path of the field solve. Since LPFx binds the transport backend to each communication operation independently (Section 5.2), the vendor libraries (NCCL and RCCL) and MPI can be combined within one run. This granularity enables a future investigation of hybrid xCCL+MPI configurations, where the optimal backend is selected for each communication pattern and scale, particularly in the latency-dominated strong-scaling regime. The mixed-precision infrastructure separates storage and compute types for each quantity, and the precision level is selected at run time among the compiled pairs. Single-precision positions retain exact energy conservation, as shown by the reconnection run. Lower precision for velocities, fields, and moments introduces FP32 round-off in the energy balance, in exchange for halving the memory footprint and traffic of the corresponding arrays. The infrastructure therefore enables a systematic investigation of the accuracy and performance of each precision level. This also concerns the deposition, since most recent matrix engines hardware targets lower-precision formats. Computing the cell-local products in lower precision, or emulating FP64 products on lower-precision matrix engines, could extend the speedup to future architectures, provided that the accuracy required by the energy balance is preserved.

10. Conclusions In this paper, we presented TPIC, a GPU-native implicit PIC code for massive space plasma physics simulations. TPIC is based on ECSIM, and it is built on the LPFx library for portable distributed-memory particle and field data structures and communication on NVIDIA and AMD GPUs. We extended the exact energy conservation of ECSIM from periodic to closed domains, by cancelling the spurious boundary term of the discrete Poynting theorem at PEC and PMC walls. TPIC conserves the total energy to machine precision in all tests, including a three-dimensional magnetic reconnection simulation with PEC walls and 1.23×1010 particles, which is also consistent with iPIC3D results. We also verified TPIC on a wide range of space plasma conventional test cases, including Whistler waves, Weibel instability, quasi neutral plasma expansion, and wave turbulence, showing that the simulations match the analytical theory and results previously published in literature. We executed the ECSIM mass-matrix deposition on the FP64 matrix engines of NVIDIA and AMD GPUs. Combined with particle sorting, the deposition is up to 99× faster on GH200 and 306× faster on MI250x than the unsorted per-particle kernel with global atomics, with the matrix engines contributing a factor 2.7-3.8×. TPIC scales quasi-optimally on LUMI-G, with a weak scaling efficiency of 96.8% on 4,096 GCDs, and a strong scaling efficiency above 95% up to 512 GCDs.

L. Pennati et al.: Preprint submitted to Elsevier

Page 29 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

A natural algorithmic extension is to develop open boundary conditions for ECSIM with a consistent discrete energy balance. Additionally, the design of TPIC and LPFx opens three major directions for future investigation. First, particle sorting takes half of the time step, while most particles remain in their cell between steps. Incremental sorting of only the particles that change cell, or more efficient sorting algorithms, could therefore greatly reduce the simulation time. Second, LPFx binds a transport backend to each communication pattern, enabling hybrid configurations of xCCL (NCCL, RCCL) and MPI, where the optimal backend is selected for each pattern. Third, the mixed-precision infrastructure enables a systematic study of reduced precision for particles, fields, and the mass matrix, including lower-precision matrix engines for the deposition.

Acknowledgments This work is funded by the European Union. This work has received funding from the European High Performance Computing Joint Undertaking (JU) and Sweden, Finland, Germany, Greece, France, Slovenia, Spain, and the Czech Republic under grant agreement No. 101093261, Plasma-PEPSC.

L. Pennati et al.: Preprint submitted to Elsevier

Page 30 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

A. Discrete Poynting identity at a wall Here we derive Eq. 12 from the stencil of Eq. 10. We first prove that  † is the adjoint of  on a periodic mesh. Then, reducing the bounded mesh to a periodic one we obtain an exact expression for the discrete Poynting theorem at boundary, valid for arbitrary ghost values and multiple walls in any direction. Evaluating it for walls normal to one direction gives Eq. 12.

Mesh and inner products Let 𝑒 denote a direction along which the system is bounded by walls. Nodes are indexed by 𝑖 = 0, … , 𝑁𝑒 , and cell centers by 𝑖 + 12 , with 𝑖 = 0, … , 𝑁𝑒 − 1 for the interior cells. Walls are located at the node

planes 𝑖 = 0 and 𝑖 = 𝑁𝑒 , while the ghost cells are at − 12 and 𝑁𝑒 + 21 . We assume the other directions to be periodic. We denote by  the interior cells, and by  the ghost cells, i.e., the cells outside the domain with at least one vertex in ∏ it. We indicate wall nodes by 𝜕Ω. The nodal control volume is 𝑉𝑔 = 𝑉𝑐 𝑒 𝑤𝑒𝑔 , with { 1∕2, 𝑔 on a wall normal to 𝑒, 𝑤𝑒𝑔 = (35) 1, otherwise, so that 𝑉𝑔 = 𝑉𝑐 ∕2, 𝑉𝑐 ∕4, and 𝑉𝑐 ∕8 at face, edge, and corner nodes. The cell and nodal inner products are ∑ ∑ 𝑉𝑔 𝐮𝑔 ⋅ 𝐯𝑔 , 𝐮𝑘 ⋅ 𝐯 𝑘 , ⟨𝐮, 𝐯⟩𝑁 = ⟨𝐮, 𝐯⟩𝐶 = 𝑉𝑐

(36)

𝑔

𝑘∈

where 𝑔 runs over all nodes, wall nodes included.

Operators Along each direction 𝑑 ∈ {𝑥, 𝑦, 𝑧}, with coordinate 𝜉 and mesh spacing Δ𝑑 , we define the centraldifference and mean operators [ ( ] ( ( ( Δ ) Δ ) Δ ) Δ ) (𝜇𝑑 𝑓 )(𝜉) = 12 𝑓 𝜉 + 2𝑑 + 𝑓 𝜉 − 2𝑑 , (37) (𝛿𝑑 𝑓 )(𝜉) = 𝑓 𝜉 + 2𝑑 − 𝑓 𝜉 − 2𝑑 , which map node values to cell values along 𝑑 and vice versa. Operators along different directions commute, and Eq. 10 can be written as 𝜕𝑥ℎ =

𝛿𝑥 𝜇𝑦 𝜇𝑧 Δ𝑥

𝜕𝑦ℎ =

,

𝜇𝑥 𝛿𝑦 𝜇𝑧 Δ𝑦

,

𝜕𝑧ℎ =

𝜇𝑥 𝜇𝑦 𝛿𝑧 Δ𝑧

(38)

,

both for the node-to-cell derivatives in  and for the cell-to-node derivatives in  † . Thus, it holds ∑ ∑ (𝐄)𝜆 = 𝜖𝜆𝑑𝜈 𝜕𝑑ℎ 𝐸𝜈 , ( † 𝐁)𝜈 = 𝜖𝜈𝑑𝜆 𝜕𝑑ℎ 𝐵𝜆 , 𝑑,𝜈

(39)

𝑑,𝜆

with 𝜖𝜆𝑑𝜈 the Levi-Civita symbol. At the wall nodes,  † uses the ghost cells.

Adjointness on a periodic mesh On a periodic mesh, every direction has as many nodes as cells, and 𝑉𝑔 = 𝑉𝑐 . Here we show that operators  and  † then commute, i.e. ⟨𝐄, 𝐁⟩𝐶 = ⟨𝐄,  † 𝐁⟩𝑁

for all 𝐄, 𝐁.

(40)

Along one periodic direction, let 𝑎 be node values and 𝑏 cell values, and let 𝐹𝑖 = 𝑎𝑖 𝑏𝑖− 1 at the nodes. Expanding 2

Eq. 37 gives, at every node 𝑖, 𝑏𝑖+ 1 (𝛿𝑎)𝑖+ 1 + 𝑎𝑖 (𝛿𝑏)𝑖 = 𝐹𝑖+1 − 𝐹𝑖 , 2 2 ( ) 𝑏𝑖+ 1 (𝜇𝑎)𝑖+ 1 − 𝑎𝑖 (𝜇𝑏)𝑖 = 12 𝐹𝑖+1 − 𝐹𝑖 . 2

(41)

2

Summing over one period direction, the right-hand sides telescope to zero, and ∑ ∑ ∑ ∑ 𝑏𝑖+ 1 (𝛿𝑎)𝑖+ 1 = − 𝑎𝑖 (𝛿𝑏)𝑖 , 𝑏𝑖+ 1 (𝜇𝑎)𝑖+ 1 = 𝑎𝑖 (𝜇𝑏)𝑖 , 𝑖

2

2

𝑖

L. Pennati et al.: Preprint submitted to Elsevier

𝑖

2

2

(42)

𝑖

Page 31 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

so 𝛿𝑑 is skew-adjoint and 𝜇𝑑 is self-adjoint. These identities hold with the indices along the other directions fixed, whatever their staggering. Applying them ∑ ∑ successively along 𝑥, 𝑦, and 𝑧 and using that the operators commute gives 𝑘 𝑏𝑘 (𝜕𝑑ℎ 𝑎)𝑘 = − 𝑔 𝑎𝑔 (𝜕𝑑ℎ 𝑏)𝑔 for each 𝑑. With Eq. 39 and 𝜖𝜆𝑑𝜈 = −𝜖𝜈𝑑𝜆 , ∑∑ ⟨𝐄, 𝐁⟩𝐶 = 𝑉𝑐 𝜖𝜆𝑑𝜈 𝐵𝜆,𝑘 (𝜕𝑑ℎ 𝐸𝜈 )𝑘 𝑘 𝜆,𝑑,𝜈

= 𝑉𝑐

∑∑

𝜖𝜈𝑑𝜆 𝐸𝜈,𝑔 (𝜕𝑑ℎ 𝐵𝜆 )𝑔 = ⟨𝐄,  † 𝐁⟩𝑁 .

𝑔 𝜆,𝑑,𝜈

Remainder for arbitrary walls We embed the domain in a periodic mesh with 𝑁𝑒 +2 cells along each wall direction ̃ equal 𝐄 at the nodes of the domain and vanish 𝑒, so that the ghost layers at opposite walls are distinct. On this mesh, let 𝐄 ̃ elsewhere, and let 𝐁 equal 𝐁, ghost values included, on  ∪  and vanish elsewhere. These extensions truncate the periodic sums but, by locality of the stencils, leave the curls unchanged where the inner products use them. Specifically, ̃ †𝐁 ̃ =  † 𝐁. Similarly, ( † 𝐁)𝑔 involves the cells with vertex 𝑔 in the domain, i.e. cells in ∪, where by definition of 𝐁, ̃ ̃ (𝐄)𝑘 involves the vertices of 𝑘, which for 𝑘 ∈  are nodes of the domain, so  𝐄 = 𝐄 on . On a ghost cell,  𝐄 involves 𝐄 only at the vertices of the cell that lie on the wall. ̃ and 𝐁, ̃ the cell sum reduces to  ∪  and the node sum, with uniform weights 𝑉𝑐 , to the nodes Applying Eq. 40 to 𝐄 of the domain: ∑ ∑ ̃ 𝑘 = ⟨𝐄,  † 𝐁⟩ ̃ 𝑁 = 𝑉𝑐 ̃ 𝐁⟩𝐶 = ⟨𝐄, 𝐁⟩𝐶 + 𝑉𝑐 𝐄𝑔 ⋅ ( † 𝐁)𝑔 . (43) 𝐁𝑘 ⋅ ( 𝐄) ⟨ 𝐄, 𝑔

𝑘∈

̃ 𝑁 defined on the periodic domain differs from ⟨𝐄,  † 𝐁⟩𝑁 only through the weights of the The inner product ⟨𝐄,  † 𝐁⟩ wall nodes, so ∑ ̃ 𝑁 − ⟨𝐄,  † 𝐁⟩𝑁 = ⟨𝐄,  † 𝐁⟩ (𝑉𝑐 − 𝑉𝑔 ) 𝐄𝑔 ⋅ ( † 𝐁)𝑔 . (44) 𝑔∈𝜕Ω

We then express ⟨𝐄, 𝐁⟩𝐶 as ̃ 𝐁⟩𝐶 − 𝑉𝑐 ⟨𝐄, 𝐁⟩𝐶 = ⟨ 𝐄,

∑

̃ 𝑘 = ⟨𝐄,  † 𝐁⟩ ̃ 𝑁 − 𝑉𝑐 𝐁𝑘 ⋅ ( 𝐄)

∑

̃ 𝑘. 𝐁𝑘 ⋅ ( 𝐄)

Combining Eq. 44 and Eq. 45 we obtain ∑ ∑ ̃ 𝑘. ⟨𝐄, 𝐁⟩𝐶 − ⟨𝐄,  † 𝐁⟩𝑉 = (𝑉𝑐 − 𝑉𝑔 ) 𝐄𝑔 ⋅ ( † 𝐁)𝑔 − 𝑉𝑐 𝐁𝑘 ⋅ ( 𝐄) 𝑔∈𝜕Ω

(45)

𝑘∈

𝑘∈

(46)

𝑘∈

The first term is the part of the nodal control volume outside the domain, i.e. 𝑉𝑐 − 𝑉𝑔 = 𝑉𝑐 ∕2, 3𝑉𝑐 ∕4, and 7𝑉𝑐 ∕8 at face, edge, and corner nodes, while the second is the ghost-cell contribution, which the uniform cell sum counts and ⟨⋅, ⋅⟩𝐶 does not. No assumption on the ghost values or on the arrangement of walls has been made. On a periodic mesh, 𝜕Ω and  are empty and Eq. 46 reduces to Eq. 40. Eq. 46 is the discrete counterpart of ∫ (∇ × 𝐄) ⋅ 𝐁 − ∫ 𝐄 ⋅ (∇ × 𝐁) = ∮ (𝐄 × 𝐁) ⋅ 𝐧̂ 𝑑𝑆.

Evaluation at a wall Consider walls normal to 𝑥, with 𝑛̂ 𝑥 = −1 at 𝑖 = 0, 𝑛̂ 𝑥 = +1 at 𝑖 = 𝑁𝑥 , 𝐧̂ = 𝑛̂ 𝑥 𝐞𝑥 , Δ𝑛 = Δ𝑥, and 𝐴𝑔 = Δ𝑦 Δ𝑧. Then  consists of the two ghost layers and 𝑉𝑐 − 𝑉𝑔 = 𝑉𝑐 ∕2 = 𝐴𝑔 Δ𝑛 ∕2 at the wall nodes. Let 𝐼 and 𝐺 denote the first interior and the ghost layer, (𝐼, 𝐺) = ( 21 , − 12 ) at 𝑖 = 0 and (𝐼, 𝐺) = (𝑁𝑥 − 12 , 𝑁𝑥 + 21 ) at 𝑖 = 𝑁𝑥 . For a cell quantity ℎ, the mean and the jump across the wall, ( ) ℎ = 12 ℎ𝐼 + ℎ𝐺 , [[ℎ]] = ℎ𝐼 − ℎ𝐺 , (47) are defined at the centers of the wall faces, with the jump taken as interior minus ghost at both walls. Symbols ℎ𝑔 , [[ℎ]]𝑔 , ℎ = 𝛿 𝜇 ∕Δ𝑦 and ℎ𝐺,𝑔 denote ℎ, [[ℎ]], and ℎ𝐺 interpolated to the wall node 𝑔 with 𝜇𝑦 𝜇𝑧 . The tangential derivatives 𝜕𝑡,𝑦 𝑦 𝑧 ∑ ∑ ℎ = 𝜇 𝛿 ∕Δ𝑧 map face values to wall nodes and back, and (∇ℎ × 𝐚) = ℎ 𝑎 . Since the wall and 𝜕𝑡,𝑧 𝜖 𝜕 𝑦 𝑧 𝜈 𝜈𝑑𝜆 𝜆 𝑑∈{𝑦,𝑧} 𝜆 𝑡 𝑡,𝑑 L. Pennati et al.: Preprint submitted to Elsevier

Page 32 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers

is periodic, the same argument leads to Eq. 40, in two dimensions, for face-centred 𝐚 and wall-node 𝐮, with both sums over one wall, gives ∑ ∑ ∑ ∑ 𝐚 ⋅ (𝜇𝑦 𝜇𝑧 𝐮) = (𝜇𝑦 𝜇𝑧 𝐚)𝑔 ⋅ 𝐮𝑔 , 𝐚 ⋅ (∇ℎ𝑡 × 𝐮) = 𝐮𝑔 ⋅ (∇ℎ𝑡 × 𝐚)𝑔 . (48) 𝑔

faces

faces

𝑔

̃ vanishes at the outer nodes, At a wall node, 𝛿𝑥 𝐁 = −𝑛̂ 𝑥 [[𝐁]] and 𝜇𝑥 𝐁 = 𝐁. On the adjacent ghost layer, where 𝐄 1 ̃ ̃ 𝛿𝑥 𝐄 = −𝑛̂ 𝑥 𝐄𝑔 and 𝜇𝑥 𝐄 = 2 𝐄𝑔 . The stencils in Eq. 46 therefore reduce to 1 𝐧̂ × [[𝐁]]𝑔 + (∇ℎ𝑡 × 𝐁)𝑔 , Δ𝑛 ̃ 𝐺 = − 1 𝐧̂ × (𝜇𝑦 𝜇𝑧 𝐄) + 1 ∇ℎ × 𝐄. ( 𝐄) Δ𝑛 2 𝑡

( † 𝐁)𝑔 = −

(49)

̂ and move the tangential operators We insert them into Eq. 46, use 𝐄 ⋅ (𝐧̂ × 𝐚) = −(𝐄 × 𝐚) ⋅ 𝐧̂ and 𝐚 ⋅ (𝐧̂ × 𝐄) = (𝐄 × 𝐚) ⋅ 𝐧, of the ghost sum from 𝐄 to 𝐁𝐺 with Eq. 48: [( ∑ ) 𝐴𝑔 𝐄𝑔 × ( 1 [[𝐁]]𝑔 + 𝐁𝐺,𝑔 ) ⋅ 𝐧̂ ⟨𝐄, 𝐁⟩𝐶 − ⟨𝐄,  † 𝐁⟩𝑉 = 2

𝑔∈𝜕Ω

+

( ) ] Δ𝑛 𝐄𝑔 ⋅ ∇ℎ𝑡 × (𝐁 − 𝐁𝐺 ) 𝑔 . 2

Since 21 [[𝐁]] + 𝐁𝐺 = 𝐁 and 𝐁 − 𝐁𝐺 = 21 [[𝐁]], ⟨𝐄, 𝐁⟩𝐶 − ⟨𝐄,  † 𝐁⟩𝑉 =

[ ) ] ( Δ 𝐴𝑔 (𝐄𝑔 × 𝐁𝑔 ) ⋅ 𝐧̂ + 𝑛 𝐄𝑔 ⋅ ∇ℎ𝑡 × [[𝐁]] 𝑔 , 4 𝑔∈𝜕Ω ∑

which is Eq. 12. With walls in several directions, Eq. 46 remains exact, but Eq. 48 leaves terms on the boundary of each face, which add to those of the edge and corner nodes and ghosts and do not vanish for arbitrary ghost values. For the TPIC closure, the whole right-hand side of Eq. 46 vanishes, so its explicit form is not needed.

Consistency Assuming a smooth 𝐁 field across the wall, with a Taylor expansion we find 𝐁 = 𝐁|𝜕Ω + 𝑂(Δ2𝑛 ) and

[[𝐁]] = −Δ𝑛 𝜕𝑛 𝐁 + 𝑂(Δ3𝑛 ). The first term of Eq. 12 then converges with second order accuracy to the continuous flux Δ2

̂ while the second is an 𝑂(Δ2𝑛 ) discretization artifact, equal to − 𝑛 𝐄 ⋅ (∇𝑡 × 𝜕𝑛 𝐁) to leading order. density (𝐄 × 𝐁) ⋅ 𝐧, 4

TPIC closure At a wall normal to 𝑒, let 𝑝𝑒 = +1 for PEC and 𝑝𝑒 = −1 for PMC, and { 𝜎𝑒 (𝐵𝜆 ) =

−𝑝𝑒 , 𝑝𝑒 ,

𝜆 = 𝑒, 𝜆 ≠ 𝑒,

𝜏𝑒 (𝐸𝜈 ) = −𝜎𝑒 (𝐵𝜈 ).

(50)

The closure of Table 1 sets each ghost value of 𝐵𝜆 to the mirror image of the corresponding interior value times 𝜎𝑒 (𝐵𝜆 ), with edge and corner ghosts carrying the product of the parities of the adjacent walls, and sets to zero, at the wall nodes, the components 𝐸𝜈 with 𝜏𝑒 (𝐸𝜈 ) = −1.

References Angus, J.R., Farmer, W., Friedman, A., Ghosh, D., Grote, D., Larson, D., Link, A., 2023. An implicit particle code with exact energy and charge conservation for electromagnetic studies of dense plasmas. Journal of Computational Physics 491, 112383. URL: https://www. sciencedirect.com/science/article/pii/S0021999123004783, doi:https://doi.org/10.1016/j.jcp.2023.112383. Ballerini, Giulio, Lavorenti, Federico, Califano, Francesco, Henri, Pierre, 2024. Whistler-mode waves in the tail of mercury´s magnetosphere: A numerical study. A&A 687, A204. URL: https://doi.org/10.1051/0004-6361/202349093, doi:10.1051/0004-6361/202349093. Barnes, D., Chacón, L., 2021. Finite spatial-grid effects in energy-conserving particle-in-cell algorithms. Computer Physics Communications 258, 107560. URL: https://www.sciencedirect.com/science/article/pii/S001046552030268X, doi:https://doi.org/10.1016/j. cpc.2020.107560.

L. Pennati et al.: Preprint submitted to Elsevier

Page 33 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Bird, R., Tan, N., Luedtke, S.V., Harrell, S.L., Taufer, M., Albright, B., 2022. Vpic 2.0: Next generation particle-in-cell simulations. IEEE Transactions on Parallel and Distributed Systems 33, 952–963. doi:10.1109/TPDS.2021.3084795. Birdsall, C.K., Langdon, A.B., 1991. Plasma Physics via Computer Simulation. Birn, J., Drake, J.F., Shay, M.A., Rogers, B.N., Denton, R.E., Hesse, M., Kuznetsova, M., Ma, Z.W., Bhattacharjee, A., Otto, A., Pritchett, P.L., 2001. Geospace environmental modeling (GEM) magnetic reconnection challenge. Journal of Geophysical Research: Space Physics 106, 3715–3719. URL: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/1999JA900449, doi:10.1029/1999JA900449. Blelloch, G.E., 1993. Prefix sums and their applications, in: Reif, J.H. (Ed.), Synthesis of Parallel Algorithms. Morgan Kaufmann. Bowers, K.J., Albright, B.J., Yin, L., Bergen, B., Kwan, T.J.T., 2008. Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation. Physics of Plasmas 15, 055703. URL: https://doi.org/10.1063/1.2840133, doi:10.1063/1.2840133. Brackbill, J., Forslund, D., 1982. An implicit method for electromagnetic plasma simulation in two dimensions. Journal of Computational Physics 46, 271–308. URL: https://www.sciencedirect.com/science/article/pii/002199918290016X, doi:https://doi.org/10.1016/ 0021-9991(82)90016-X. Buneman, O., 1993. TRISTAN: The 3-D electromagnetic particle code, in: Matsumoto, H., Omura, Y. (Eds.), Computer Space Plasma Physics: Simulation Techniques and Software. Terra Scientific, Tokyo, pp. 67–84. URL: https://www.terrapub.co.jp/e-library/cspp/. Burch, J.L., Torbert, R.B., Phan, T.D., Chen, L.J., Moore, T.E., Ergun, R.E., Eastwood, J.P., Gershman, D.J., Cassak, P.A., Argall, M.R., Wang, S., Hesse, M., Pollock, C.J., Giles, B.L., Nakamura, R., Mauk, B.H., Fuselier, S.A., Russell, C.T., Strangeway, R.J., Drake, J.F., Shay, M.A., Khotyaintsev, Y.V., Lindqvist, P.A., Marklund, G., Wilder, F.D., Young, D.T., Torkar, K., Goldstein, J., Dorelli, J.C., Avanov, L.A., Oka, M., Baker, D.N., Jaynes, A.N., Goodrich, K.A., Cohen, I.J., Turner, D.L., Fennell, J.F., Blake, J.B., Clemmons, J., Goldman, M., Newman, D., Petrinec, S.M., Trattner, K.J., Lavraud, B., Reiff, P.H., Baumjohann, W., Magnes, W., Steller, M., Lewis, W., Saito, Y., Coffey, V., Chandler, M., 2016. Electron-scale measurements of magnetic reconnection in space. Science 352, aaf2939. URL: https://www.science.org/doi/abs/ 10.1126/science.aaf2939, doi:10.1126/science.aaf2939. Chen, G., Chacón, L., Barnes, D., 2011. An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm. Journal of Computational Physics 230, 7018–7036. URL: https://www.sciencedirect.com/science/article/pii/S0021999111003421, doi:https://doi.org/10.1016/j.jcp.2011.05.031. Chen, Y., Tóth, G., Hietala, H., Vines, S.K., Zou, Y., Nishimura, Y., Silveira, M.V.D., Guo, Z., Lin, Y., Markidis, S., 2020. Magnetohydrodynamic with embedded particle-in-cell simulation of the Geospace Environment Modeling dayside kinetic processes challenge event. Earth and Space Science 7, e2020EA001331. URL: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2020EA001331, doi:10.1029/ 2020EA001331. Chen, Y., Tóth, G., Zhou, H., Wang, X., 2023. Fleks: A flexible particle-in-cell code for multi-scale plasma simulations. Computer Physics Communications 287, 108714. URL: https://www.sciencedirect.com/science/article/pii/S0010465523000590, doi:https: //doi.org/10.1016/j.cpc.2023.108714. Cohen, B.I., Langdon, A., Hewett, D.W., Procassini, R.J., 1989. Performance and optimization of direct implicit particle simulation. Journal of Computational Physics 81, 151–168. URL: https://www.sciencedirect.com/science/article/pii/0021999189900685, doi:https://doi.org/10.1016/0021-9991(89)90068-5. Daughton, W., 1998. Kinetic theory of the drift kink instability in a current sheet. Journal of Geophysical Research: Space Physics 103, 29429–29443. doi:10.1029/1998JA900028. Davidson, R.C., Hammer, D.A., Haber, I., Wagner, C.E., 1972. Nonlinear development of electromagnetic instabilities in anisotropic plasmas. The Physics of Fluids 15, 317–333. URL: https://doi.org/10.1063/1.1693910, doi:10.1063/1.1693910. Dorozhkina, D.S., Semenov, V.E., 1998a. Exact solution of the problem of quasineutral expansion into vacuum of a localized collisionless plasma with cold ions. Journal of Experimental and Theoretical Physics Letters 67, 573–578. URL: https://doi.org/10.1134/1.567728, doi:10.1134/1.567728. Dorozhkina, D.S., Semenov, V.E., 1998b. Exact solution of vlasov equations for quasineutral expansion of plasma bunch into vacuum. Phys. Rev. Lett. 81, 2691–2694. URL: https://link.aps.org/doi/10.1103/PhysRevLett.81.2691, doi:10.1103/PhysRevLett.81.2691. Fedeli, L., Huebl, A., Boillod-Cerneux, F., Clark, T., Gott, K., Hillairet, C., Jaure, S., Leblanc, A., Lehe, R., Myers, A., Piechurski, C., Sato, M., Zaim, N., Zhang, W., Vay, J.L., Vincenti, H., 2022. Pushing the frontier in the design of laser-based electron accelerators with groundbreaking meshrefined particle-in-cell simulations on exascale-class supercomputers, in: SC22: International Conference for High Performance Computing, Networking, Storage and Analysis, pp. 1–12. doi:10.1109/SC41404.2022.00008. Fonseca, R.A., Silva, L.O., Tsung, F.S., Decyk, V.K., Lu, W., Ren, C., Mori, W.B., Deng, S., Lee, S., Katsouleas, T., Adam, J.C., 2002. Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators, in: Sloot, P.M.A., Hoekstra, A.G., Tan, C.J.K., Dongarra, J.J. (Eds.), Computational Science — ICCS 2002, Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 342–351. Gary, S.P., Chang, O., Wang, J., 2012. Forward cascade of whistler turbulence: Three-dimensional particle-in-cell simulations. The Astrophysical Journal 755, 142. URL: https://doi.org/10.1088/0004-637X/755/2/142, doi:10.1088/0004-637X/755/2/142. Germaschewski, K., Fox, W., Abbott, S., Ahmadi, N., Maynard, K., Wang, L., Ruhl, H., Bhattacharjee, A., 2016. The plasma simulation code: A modern particle-in-cell code with patch-based load-balancing. Journal of Computational Physics 318, 305–326. URL: https: //www.sciencedirect.com/science/article/pii/S0021999116301413, doi:https://doi.org/10.1016/j.jcp.2016.05.013. Gonzalez-Herrero, D., Boella, E., Lapenta, G., 2018. Performance analysis and implementation details of the energy conserving semi-implicit method code (ecsim). Computer Physics Communications 229, 162–169. URL: https://www.sciencedirect.com/science/article/ pii/S0010465518300924, doi:https://doi.org/10.1016/j.cpc.2018.03.020. Hakobyan, H., Spitkovsky, A., Chernoglazov, A., Philippov, A., Groselj, D., Mahlmann, J., Bacchini, F., Vanthieghem, A., 2024. Princetonuniversity/tristan-mp-v2: v2.9. URL: https://doi.org/10.5281/zenodo.14564000, doi:10.5281/zenodo.14564000. Hesse, M., Cassak, P.A., 2020. Magnetic reconnection in the space sciences: Past, present, and future. Journal of Geophysical Research: Space Physics 125, e2018JA025935. doi:10.1029/2018JA025935.

L. Pennati et al.: Preprint submitted to Elsevier

Page 34 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Hockney, R.W., Eastwood, J.W., 1988. Computer Simulation Using Particles. CRC Press, Boca Raton, FL. doi:10.1201/9781439822050. often cited by original publication year; DOI corresponds to CRC/T&F electronic edition. Innocenti, M.E., Johnson, A., Markidis, S., Amaya, J., Deca, J., Olshevsky, V., Lapenta, G., 2017. Progress towards physics-based space weather forecasting with exascale computing. Advances in Engineering Software 111, 3–17. Karimabadi, H., Roytershteyn, V., Vu, H.X., Omelchenko, Y.A., Scudder, J., Daughton, W., Dimmock, A., Nykyri, K., Wan, M., Sibeck, D., Tatineni, M., Majumdar, A., Loring, B., Geveci, B., 2014. The link between shocks, turbulence, and magnetic reconnection in collisionless plasmas. Physics of Plasmas 21, 062308. URL: https://doi.org/10.1063/1.4882875, doi:10.1063/1.4882875. Kiyani, K.H., Osman, K.T., Chapman, S.C., 2015. Dissipation and heating in solar wind turbulence: from the macro to the micro and back again. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 373, 20140155. URL: https: //doi.org/10.1098/rsta.2014.0155, doi:10.1098/rsta.2014.0155. Langdon, A., 1970. Effects of the spatial grid in simulation plasmas. Journal of Computational Physics 6, 247–267. URL: https://www. sciencedirect.com/science/article/pii/0021999170900240, doi:https://doi.org/10.1016/0021-9991(70)90024-0. Langdon, A., Cohen, B.I., Friedman, A., 1983. Direct implicit large time-step particle simulation of plasmas. Journal of Computational Physics 51, 107–138. URL: https://www.sciencedirect.com/science/article/pii/0021999183900839, doi:https://doi.org/10.1016/ 0021-9991(83)90083-9. Lapenta, G., 2012. Particle simulations of space weather. Journal of Computational Physics 231, 795–821. URL: https:// www.sciencedirect.com/science/article/pii/S0021999111001860, doi:https://doi.org/10.1016/j.jcp.2011.03.035. special Issue: Computational Plasma Physics. Lapenta, G., 2017. Exactly energy conserving semi-implicit particle in cell formulation. Journal of Computational Physics 334, 349– 366. URL: https://www.sciencedirect.com/science/article/pii/S0021999117300128, doi:https://doi.org/10.1016/j. jcp.2017.01.002. Lapenta, G., 2023. Advances in the implementation of the exactly energy conserving semi-implicit (ecsim) particle-in-cell method. Physics 5, 72–89. URL: https://www.mdpi.com/2624-8174/5/1/7, doi:10.3390/physics5010007. Lapenta, G., Brackbill, J.U., 2002. Nonlinear evolution of the lower hybrid drift instability: Current sheet thinning and kinking. Physics of Plasmas 9, 1544–1554. URL: https://doi.org/10.1063/1.1464149, doi:10.1063/1.1464149. Lapenta, G., Brackbill, J.U., Ricci, P., 2006. Kinetic approach to microscopic-macroscopic coupling in space and laboratory plasmas. Physics of Plasmas 13, 055904. URL: https://doi.org/10.1063/1.2173623, doi:10.1063/1.2173623. Lee, R.P., Pierce, J.R., Miller, K.G., Almanza, M., Tableman, A., Decyk, V.K., Fonseca, R.A., Alves, E.P., Mori, W.B., 2025. Acceleration of the particle-in-cell code osiris with graphics processing units. Journal of Plasma Physics 91, E8. doi:10.1017/S0022377824001569. Markidis, S., Henri, P., Lapenta, G., Divin, A., Goldman, M., Newman, D., Laure, E., 2013. Kinetic simulations of plasmoid chain dynamics. Physics of Plasmas 20, 082105. URL: https://doi.org/10.1063/1.4817286, doi:10.1063/1.4817286. Markidis, S., Hu, A., Peng, I., Pennati, L., Lumsden, I., Yokelson, D., Brink, S., Pearce, O., Scogland, T.R., de Supinski, B.R., et al., 2025. Exascale implicit kinetic plasma simulations on el˜ capitan for solving the micro-macro coupling in magnetospheric physics. arXiv preprint arXiv:2507.20719 . Markidis, S., Lapenta, G., 2011. The energy conserving particle-in-cell method. Journal of Computational Physics 230, 7037–7052. URL: https: //www.sciencedirect.com/science/article/pii/S0021999111003445, doi:https://doi.org/10.1016/j.jcp.2011.05.033. Markidis, S., Lapenta, G., Rizwan-uddin, 2010. Multi-scale simulations of plasma with ipic3d. Mathematics and Computers in Simulation 80, 1509– 1519. URL: https://www.sciencedirect.com/science/article/pii/S0378475409002444, doi:https://doi.org/10.1016/j. matcom.2009.08.038. multiscale modeling of moving interfaces in materials. Mason, R.J., 1981. Implicit moment particle simulation of plasmas. Journal of Computational Physics 41, 233–244. URL: https://www. sciencedirect.com/science/article/pii/0021999181900942, doi:https://doi.org/10.1016/0021-9991(81)90094-2. Matsumoto, H., Omura, Y., 1985. Particle simulation of electromagnetic waves and its application to space plasmas doi:10.1007/ 978-94-009-5321-5_2. Okuda, H., 1972. Nonphysical noises and instabilities in plasma simulation due to a spatial grid. Journal of Computational Physics 10, 475–486. URL: https://www.sciencedirect.com/science/article/pii/0021999172900484, doi:https://doi.org/10.1016/ 0021-9991(72)90048-4. Peng, I.B., Markidis, S., Vaivads, A., Vencels, J., Amaya, J., Divin, A., Laure, E., Lapenta, G., 2015. The formation of a magnetosphere with implicit particle-in-cell simulations. Procedia Computer Science 51, 1178–1187. URL: https://www.sciencedirect.com/science/article/ pii/S1877050915010960, doi:https://doi.org/10.1016/j.procs.2015.05.288. international Conference On Computational Science, ICCS 2015. Pennati, L., Markidis, S., 2026. Mass matrix assembly on tensor cores for implicit particle-in-cell methods. Journal of Computational Physics , 115442URL: https://www.sciencedirect.com/science/article/pii/S0021999126007916, doi:https://doi.org/10.1016/j. jcp.2026.115442. Ricci, P., Lapenta, G., Brackbill, J., 2002. A simplified implicit maxwell solver. Journal of Computational Physics 183, 117–141. URL: https: //www.sciencedirect.com/science/article/pii/S0021999102971706, doi:https://doi.org/10.1006/jcph.2002.7170. Saito, S., Gary, S.P., Li, H., Narita, Y., 2008. Whistler turbulence: Particle-in-cell simulations. Physics of Plasmas 15, 102305. URL: https://doi.org/10.1063/1.2997339, doi:10.1063/1.2997339. Spitkovsky, A., 2005. Simulations of relativistic collisionless shocks: shock structure and particle acceleration. AIP Conference Proceedings 801, 345–350. URL: https://doi.org/10.1063/1.2141897, doi:10.1063/1.2141897. Taflove, A., Hagness, S., 2005. Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd edition.. volume 2062. doi:10.1016/B978-012170960-0/50046-3. Vay, J.L., Acciarri, M., Almgren, A., Amorim, L.D., Andriyash, I., Angus, J.R., Antoun, T., Belkin, D., Bizzozero, D., Blelly, A., Clark, S.E., Dammak, E., Fedeli, L., Formenti, A., Garten, M., Ge, L., Giacomel, L., Gott, K., Groenewald, R.E., Grote, D., Gu, J., Harrison, C., Haseeb, M.,

L. Pennati et al.: Preprint submitted to Elsevier

Page 35 of 36

TPIC: Energy-Conserving Space Plasma Simulations on GPU Supercomputers Huebl, A., Jambunathan, R., Kara-Mostefa, I., Klion, H., Kumar, P., Lehe, R., Loring, B., Marks, T., Miller, P., Myers, A., Pech, J., Rheaume, E., Richardson, G., Rowan, M.E., Sandberg, R.T., Scherpelz, P., Sinn, A., Shapoval, O., Terzani, D., Thévenet, M., Weichman, K., Yang, E., Zaim, N., Zhang, W., Zhao, Y., Zhu, K.Z., Zoni, E., 2018. WarpX: An advanced Particle-In-Cell code. https://blast-warpx.github.io. URL: https://doi.org/10.5281/zenodo.4571577, doi:10.5281/zenodo.4571577. https://doi.org/10.5281/zenodo.4571577. Wang, J., Liewer, P., Decyk, V., 1995. 3d electromagnetic plasma particle simulations on a mimd parallel computer. Computer Physics Communications 87, 35–53. URL: https://www.sciencedirect.com/science/article/pii/001046559400167Z, doi:https://doi. org/10.1016/0010-4655(94)00167-Z. particle Simulation Methods. Weibel, E.S., 1959. Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution. Phys. Rev. Lett. 2, 83–84. URL: https://link.aps.org/doi/10.1103/PhysRevLett.2.83, doi:10.1103/PhysRevLett.2.83. Yamada, M., Kulsrud, R., Ji, H., 2010. Magnetic reconnection. Rev. Mod. Phys. 82, 603–664. URL: https://link.aps.org/doi/10.1103/ RevModPhys.82.603, doi:10.1103/RevModPhys.82.603. Yu, Y., Delzanno, G.L., Jordanova, V., Peng, I.B., Markidis, S., 2018. Pic simulations of wave-particle interactions with an initial electron velocity distribution from a kinetic ring current model. Journal of Atmospheric and Solar-Terrestrial Physics 177, 169–178. URL: https://www.sciencedirect.com/science/article/pii/S1364682617302195, doi:https://doi.org/10.1016/j.jastp.2017. 07.004. dynamics of the Sun-Earth System: Recent Observations and Predictions.

L. Pennati et al.: Preprint submitted to Elsevier

Page 36 of 36