Mass Matrix Assembly on Tensor Cores for Implicit Particle-In-Cell Methods Luca Pennati, Stefano Markidis
arXiv:2604.19286v1 [cs.CE] 21 Apr 2026
KTH Royal Institute of Technology, Stockholm, 114 28, Sweden
Abstract Matrix-multiply-accumulate (MMA) units, or tensor cores, are now widespread across modern computing architectures. Yet, their use for particle-grid operators remains limited. In implicit particle methods, mass-matrix assembly is a reduction-dominated kernel in which weighted outer products of interpolation weights are accumulated over particle support. We show that this operation can be reformulated exactly, cell by cell, as a sequence of matrix products matched to hardware MMA tiles. The formulation is general with respect to interpolation order and hardware platform, and applies to both scalar mass matrices and the tensorial block mass matrix arising in implicit in the Energy-Conserving Semi-Implicit Method (ECSIM) for Particle-in-Cell simulations. We introduce particle batching and a support-group decomposition for higher-order shape functions whose stencil extends beyond a single cell, specialize the method to first- and second-order B-spline interpolation, and implement it on NVIDIA tensor cores. The resulting kernels achieve up to 3× over optimized conventional implementations and reduce end-to-end ECSIM runtime by ∼ 15%. Keywords: Mass Matrix, Particle-In-Cell, Matrix Engines, Kinetic Plasma Simulation, ECSIM, GPUs 1. Introduction Modern computing architectures increasingly own a substantial fraction of their computational throughput to hardware Matrix-Multiply-Accumulate (MMA) units, commonly referred to as tensor cores or matrix engines. Originally introduced to accelerate machine-learning workloads to perform an
MMA operation per clock cycle [1, 2, 3, 4], these units are now available across essentially all major accelerator and processor families. However, their effective use in scientific computing, still depends on the ability to recast an application kernel into a sequence of small dense matrix products with sufficient regularity to match the underlying hardware tiles. MMA units have found broader application across additional scientific domains. In particular, tensor cores have been employed in linear algebra [5] and linear solvers [6], especially in the context of low and mixed-precision algorithms. Other studies have investigated how to recast finite element methods [7] and general stencil operations [8] into MMA-friendly forms. Tensor cores have also been successfully applied to domain-specific problems, including signal processing [9], molecular docking [10], and quantum molecular dynamics [11]. However, tensor core applicability to irregular particle-grid operators remains much less explored and it is the topic of this paper. An important example of such an operator is the mass-matrix assembly arising in semi-implicit particle methods. In the Energy-Conserving SemiImplicit Method (ECSIM) [12] for Particle-In-Cel (PIC) plasma simulations, the mass matrix represents the linear response of the plasma to the electric field and enters the field solve as a grid-defined operator assembled from particle information. Its construction requires, for each particle, the accumulation over all pairs of support nodes of a weighted outer product of interpolation values, scaled in ECSIM by a particle-dependent response tensor. Although mathematically simple, this operation is computationally demanding. It is dominated by fine grained reductions and irregular scatter patterns that map poorly to conventional Single Instruction Multiple Data (SIMD) and Single Instruction Multiple Threads (SIMT) execution. In practice, it becomes one of the most expensive stages of the ECSIM cycle. The same algebraic structure also appears more broadly in mass-matrix-based particle-grid formulations, including scalar variants in related methods, most notably in the Material Point Method (MPM) for continuum mechanics [13, 14, 15], where a scalar mass matrix couples grid-node momenta. In this work, we show that the mass-matrix assembly can be reformulated exactly in a form that is naturally matched to tensor cores. We express the weighted outer-product accumulation as a tensor contraction over particles and decompose that contraction cell by cell. In this formulation, the local assembly reduces to a sequence of batched matrix products whose inner dimension coincides with the contraction dimension of MMA tiles. This approach leads to a general mapping strategy, independent of interpolation 2
order and kind of matrix engine, and applies to both the scalar mass matrix and the tensorial block structure arising in ECSIM. In this work, we focus on ECSIM as the primary and most demanding case study. The main contributions of this work are as follows. First, we present a general reformulation of the mass matrix assembly as a tensor contraction that factors into a matrix product Mab = Aak Bkb , valid for arbitrary spatial dimensions and interpolation orders, applicable to both scalar and tensorial mass matrices. Second, we introduce a formal particle-batching and support-group decomposition strategy that maps the cell-local contraction onto fixed-size MMA tiles, with particular attention to higher-order shape functions whose stencil extends beyond a single cell. Third, we describe a reference implementation on NVIDIA tensor cores for first-order (CIC) and second-order (TSC) B-spline interpolation in three dimensions, with performance benchmarks against optimized conventional GPU kernels that demonstrate an end-to-end simulation speedup. The paper is organized as follows. In Section 2 we introduce the ECSIM mass matrix, reformulate its assembly as a tensor contraction, and recall the B-spline shape functions used in this work. We also describe the hardware tensor cores abstraction used in the paper. In Section 3 we develop the general mass matrix assembly strategy, including the cell-local outer-product decomposition, particle batching, support-group decomposition and sparse stencil deposition. In Section 4 we present numerical results for our implementation on NVIDIA tensor cores, including comparisons against conventional GPU kernels and end-to-end speedup in a production kinetic plasma simulation. Section 5 discusses the implications and limitations of the proposed approach, and Section 6 is dedicated to conclusions. 2. Preliminaries In the Particle-In-Cell (PIC) method [16, 17], particles evolve in Lagrangian coordinates while field quantities live on a discrete Eulerian grid, and interpolation functions mediate all field–particle interactions. In ECSIM [12], the plasma medium response is represented by the mass matrix operator, defined on the grid, which takes the form of a 3 × 3 tensor block for every pair of grid nodes within the support of a particle shape function, thus requiring the assembly of nine components per node pair. Its construction involves accumulating, for each particle, the outer product of the correspond3
Wgp support nodes g ∈ Np p
Figure 1: Two-dimensional particle-grid coupling for the first-order (CIC) case. The particle p in the cell deposits to the four corner nodes in its support Np .
ing shape-function weight vector with itself, scaled by a particle-dependent coefficient tensor. This procedure is traditionally implemented as a particleby-particle scatter loop with fine-grained reductions, a computational pattern that maps poorly onto conventional accelerated architectures. As a result, mass-matrix assembly fails to fully exploit the computational capabilities of contemporary hardware and becomes the most time-consuming stage of the PIC cycle. 2.1. Mass matrix in the ECSIM PIC method We briefly recall the origin of the mass matrix in the ECSIM formulation by Lapenta [12], the reader is referred to that work for a complete derivation. (s) Consider a plasma described by Ns species, each represented by Np computational particles on a d-dimensional Cartesian grid G with nodes indexed by g. Each particle at position xp interacts with the grid through a compactly supported shape (or weight) function W (xp − xg ). We write Wpg ≡ W (xp − xg ),
(1)
and denote by Np the compact support of particle p, i.e. the set of grid nodes with nonzero weights: Np := { g ∈ G | Wpg ̸= 0 }.
(2)
Figure 1 illustrates the particle-grid coupling in the two-dimensional case, for a first order Cloud-In-Cell (CIC) interpolation. Figure 2 describes one cycle of the ECSIM algorithm. In the ECSIM method, the implicit coupling between the unknown electric field at time 4
n + 1/2, E n+1/2 , and the plasma current gives rise to a linear system for the field, defined on the grid nodes, of the form ! X L+ Ms E n+1/2 = b, (3) s
where L is a discrete curl-curl operator and b a right-hand-side vector depending on known quantities, such as current J, at time n. Once the electric field is known, the magnetic field is advanced via the discrete Faraday’s law. The operator Ms in Eq. 3 is the mass matrix for species s, its entries couple grid-node pairs (g, g ′ ) through a 3 × 3 tensor block: ij βs X qp αpij Wpg Wpg′ , Ms gg′ = c Vg p∈s
i, j ∈ {1, 2, 3},
(4)
where qp is the particle charge, Vg is the cell volume, βs = qs ∆t/(2ms ) is a species-dependent time-step parameter, and αpij is the (i, j) component of the particle rotation-response tensor αp . The tensor αp is defined as αp =
1 I − C(ωp ) + ωp ωpT , 2 1 + ∥ωp ∥
(5)
where ωp = βs B(xp )/c is the dimensionless magnetization vector, with B(xp ) being the magnetic field interpolated to the particle position, and C(ω) u ≡ ω × u is the skew-symmetric cross-product operator. Its assembly is the dominant cost of the implicit field solve, since for every particle one must evaluate and accumulate N 2 products of shape-function values, each multiplied by nine tensor components. Dropping the species index for notational simplicity, the mass matrix can be written in the general form ij Mgg ′ = σ
Np X
s̃ij p Wpg Wpg ′ ,
(6)
p=1 ij where σ absorbs all constant prefactors and s̃ij p ≡ qp αp is the per-particle coefficient tensor. In the scalar case (e.g., the mass matrix in MPM [13]), ij one has s̃ij p = qp δ and only a single component per node pair.
5
Mass Matrix Assembly X ij Mgg qp αpij Wpg Wpg′ ′ =
2. Particle → Grid Deposit J ; compute αp , assemble M 1. Particle Mover Update xp , vp
p
Reformulated as Matrix Engine MMA operation
3. Field Solver P L + s Ms E n+1/2 = b
∆t
B n+1 = B n −c∆t ∇g ×E n+1/2
4. Grid → Particle Interpolate E n+1/2 , B n+1 to particles
Figure 2: Diagram of the PIC cycle for the ECSIM [12] algorithm. The mass matrix calculation, the topic of this work, is highlighted in red.
2.2. Mass matrix as a tensor contraction It has been recognized in the literature that the mass matrix computation, Eq. (6), admits a natural interpretation as a tensor contraction over the particle index [18, 19]. Here we formalize this observation. Consider a spatial grid with Ng nodes and a system of Np particles. The weight matrix W ∈ RNg ×Np collects all shape-function evaluations: Wgp ≡ W (xp − xg ),
g ∈ {0, . . . , Ng −1}, p ∈ {0, . . . , Np −1}.
(7)
For each tensor component (i, j), we define the diagonal coefficient matrix ij Np ×Np Sij = diag s̃ij . (8) 0 , . . . , s̃Np −1 ∈ R With Einstein summation, the mass matrix assembly in Eq. (6) can be recast in a tensor contraction ij ij (9) Mgg ′ = Wgp s̃p Wg ′ p , where the three factors are contracted over the repeated particle index p. In matrix notation, this reads Mij = W Sij WT ∈ RNg ×Ng .
(10)
Although Mij ∈ RNg ×Ng formally, it is extremely sparse since each particle contributes to at most N × N entries (with N = (n+1)d ), thus every row 6
contains at most (2n+1)d nonzero entries. Additionally, the mass matrix is symmetric, for each fixed tensor component (i, j), ij ij Mgg ′ = Mg ′ g ,
(11)
since the product Wgp Wg′ p in Eq. (9) is invariant under interchange of g and g ′ . Each Ng × Ng block Mij is therefore a real symmetric matrix. These symmetries, combined with the compact sparsity pattern, reduce the storage from Ng2 entries to O(Ng ) per component, indexed by canonical stencil offsets. 2.3. Shape functions The shape function W (xp − xg ) in Eq. (1) assigns to each particlenode pair a non-negative interpolation weight that determines the coupling strength. On a uniform Cartesian grid with spacing ∆xµ along coordinate µ ∈ {1, . . . , d}, the standard choice is a product of one-dimensional B-splines of order n: µ d Y xp − xµg (n) ϕ W (xp − xg ) = , (12) ∆xµ µ=1 where ϕ(n) is compactly supported on [−(n+1)/2, (n+1)/2] and the support d of each particle includes N = (n+1) P grid nodes. The shape functions satisfy the partition-of-unity property g Wpg = 1 and are non-negative, ensuring conservative interpolation [16, 20]. We consider the two cases of main practical relevance: first-order CloudIn-Cell (CIC), with N = 2d support nodes per particle (N = 8 in 3D), and second-order Triangular-Shaped Cloud (TSC), with N = 3d support nodes (N = 27 in 3D). For CIC, the stencil always coincides with the 2d corner nodes of the cell containing the particle, so all particles in a cell share the same support. For TSC, the identity of the support nodes depends on the particle position within the cell. In each dimension µ, the stencil is centered on the nearest grid node, so that for fractional coordinate ξpµ ∈ [0, 1) a particle with ξpµ < 1/2 uses support nodes {j−1, j, j+1} (base offset bµ = −1), while ξpµ ≥ 1/2 gives {j, j+1, j+2} (bµ = 0). 2.4. Tensor Core Architectures For the scope of this work, we abstract a tensor core as any hardware unit that realizes the following operations. Given fixed positive integers Mt , Nt , Kt , which define the tile shape, the engine accepts operand tiles A ∈ FinMt ×Kt 7
Nt
Kt
×
Kt
+
Mt
Mt
Nt
Figure 3: Abstract matrix-engine MMA update written in tile form as D ← D + AB. The accumulator has size Mt × Nt , the left operand has size Mt × Kt , and the right operand has size Kt × Nt .
and B ∈ FinKt ×Nt , where Fin is the floating-point format of the input operands, Mt ×Nt and updates an accumulator tile D ∈ Facc , stored in a (generally wider) accumulation format Facc ⊇ Fin , according to the MMA rule Dab ← Dab + Aak Bkb ,
a ∈ {0, . . . , Mt −1}, b ∈ {0, . . . , Nt −1},
(13)
where summation over the repeated index k ∈ {0, . . . , Kt −1} is implied. Figure 3 summarizes this abstraction in the compact form of the MMA update D ← D+AB. Each invocation of Eq. (13) represents a tiled matrix multiplyaccumulate involving 2 Mt Nt Kt floating-point operations.This operation is exposed as a single MMA instruction in the programming model and is executed on specialized tensor-core hardware, yielding higher throughput than an equivalent implementation built from scalar FMA instructions [21]. Additionally, the use of a reduced-precision input format Fin (such as FP16, TF32, or BF16) with a wider accumulation format Facc (such as FP32 or FP64) allows the hardware to maximize throughput for the multiply stage while preserving numerical accuracy in the accumulation. All major contemporary computer architectures provide such tensor cores: NVIDIA tensor cores [21], AMD matrix cores [22], Intel Advanced Matrix Extensions (AMX) [23], and Google Tensor Processing Units (TPUs) [1]. These implementations differ in their supported tile shapes (Mt , Nt , Kt ), and precision formats (Fin , Facc ), but all conform to the abstract MMA interface defined by Eq. (13). 3. Methodology In Section 2.2 we recall that the mass matrix assembly is a tensor conij ij traction Mgg ′ = Wgp s̃p Wg ′ p over the particle index p (Eq. (9)), and that 8
p0
p1
pN
Wg0 p0 s̃ij p0
Wg0 p1 s̃ij p1
···
Wg0 pN s̃ij pN
p0
Wg0 p0
Wg1 p0
···
WgNg −1 p0
Wg1 p0 s̃ij p0 .. .
Wg1 p1 s̃ij p1 .. .
··· .. .
Wg1 pN s̃ij pN .. .
p1
Wg0 p1 .. .
Wg1 p1 .. .
··· .. .
WgNg −1 p1 .. .
WgNg −1 p0 s̃ij p0
WgNg −1 p1 s̃ij p1
pN
W g0 pN
W g1 pN
· · · WgNg −1 pN
· · · WgNg −1 pN s̃ij pN
×
Figure 4: Schematic of the factorization Mij = Aij B. Each column of Aij corresponds to one particle and contains that particle’s interpolation weights scaled by s̃ij p . The matching row of B contains the same particle weights unscaled. The matrix product contracts over the shared particle index and sums the per-particle outer products into the mass matrix.
the compact support of the shape functions restricts each particle’s contribution to a small block of Mij . In this section, we show how we leverage the inherent sparsity of Mij to decompose its calculation at the cell level and map the tensor contraction onto the fixed-size tile operations provided by hardware matrix engines. The derivation is carried out in full generality, independent of the interpolation order, the number of spatial dimensions, the scalar or tensorial nature of the per-particle coefficient, and the tile shapes of the matrix engine. 3.1. Cell-local tensor contraction Firstly, wo factor the three-tensor product in Eq. (9) into a two-operand matrix multiply by absorbing the diagonal coefficient into the left factor: ij ij Mgg ′ = Agp Bpg ′ ,
Aij = W Sij ∈ RNg ×Np , B = WT ∈ RNp ×Ng .
(14)
Figure 4 visualizes this factorization: column p of Aij stores the interpolation weights of particle p scaled by its coefficient s̃ij p , while row p of B stores the same particle weights without the scaling. The product Aij B therefore contracts over the shared particle index and accumulates the outer-product contribution of each particle into the grid-grid matrix. Given the compact support of the shape function W (xp −xg ), Eq. 12, each column (row) p in the global matrix Aij (B) would have at most (n + 1)d non zero entries, leading to an highly sparse matrix-matrix multiplication. We can therefore leverage the regularity and compactness of the shape functions support to decompose the mass matrix assembly in a series of cell-local tensor contractions rather than a single large global operation. 9
We assume that particles have been sorted by cell, which is typically the case for production-level simulations [24]. Let a cell contain P particles and let N = {n0 , . . . , nN −1 } be a fixed set of N grid nodes such that the support of every particle in the cell is contained in N . Restricting the global weight tensor (Eq. (7)) and coefficient tensor (Eq. (8)) to these P particles and N nodes yields cell-local matrices Ŵ ∈ RN ×P and Ŝij ∈ RP ×P . The global tensor contraction (Eq. (9)) then applies block-wise and the cell-local mass matrix is ij M̂ab = Ŵap ŝij (15) p Ŵbp . 3.2. Particle batching for MMA tiles The cell-local tensor contraction in Eq. (15) sums over all P particles in a cell, but the hardware MMA tile contracts over a fixed inner dimension Kt . Thus, we partition the P particles into ⌈P/Kt ⌉ consecutive batches of Kt particles each, with the last batch zero-padded if Kt ∤ P . For batch β containing particles {pβ,0 , . . . , pβ,Kt −1 }, we define the batch weight matrix (β) Wak ≡ Ŵa,pβ,k , W (β) ∈ RN ×Kt , (16) and the corresponding diagonal coefficient slice ij Kt ×Kt Sij,(β) = diag ŝij . pβ,0 , . . . , ŝpβ,K −1 ∈ R t
(17)
The MMA operands for batch β are ij,(β)
Aak
(β)
(β)
≡ Wak ŝij pβ,k ,
(β)
Bkb ≡ Wbk ,
(18)
with Aij,(β) ∈ RN ×Kt and B(β) ∈ RKt ×N . The full cell-local mass matrix is the sum of per-batch products ⌈P/Kt ⌉−1 ij M̂ab =
X
ij,(β)
Aak
(β)
Bkb ,
(19)
β=0
crucially, this summation maps exactly onto the hardware MMA instruction ij defined in Eq. (13). Initializing the accumulator tile to Dab ← 0, each batch β triggers the in-place update ij,(β)
ij ij Dab ← Dab + Aak
10
(β)
Bkb ,
(20)
thus, after all ⌈P/Kt ⌉ MMA calls the accumulator holds the exact cell-local ij ij . The mathematical sum over particle batches = M̂ab mass matrix: Dab in Eq. (19) is therefore realized by a loop of hardware MMA instructions that accumulate in place, requiring no intermediate storage and no explicit reduction step. When N > Mt (or N > Nt ), the N × N accumulator is covered by ⌈N/Mt ⌉ × ⌈N/Nt ⌉ tiles, each executing an independent MMA instruction per batch. The weight matrix rows are then padded to the nearest multiple of Mt to fill incomplete tiles. 3.3. Support-group decomposition The cell-local contraction in Eq. (15) assumes that all particles contributing in a given product share the same support nodes. For shape functions of order n ≥ 2, the stencil placement depends on the particle position within the cell, so two particles in the same cell may touch different (though overlapping) subsets of (n+1)d grid nodes. Let Pc = {p0 , . . . , pP −1 } be the particles in cell c. We partition Pc into G groups Π0 , . . . , ΠG−1 by grouping particles that share identical support nodes: G−1 G Pc = Πγ , ∀ p ∈ Πγ : Np = Nγ , (21) γ=0
where Nγ denotes the common support of group Πγ and |Πγ | = Pγ . Figure 5 illustrates this decomposition for second-order interpolation in two dimensions, showing the four possible 3 × 3 TSC nodal supports inside one cell. The matrix product Eq. (15) applies independently within each group. We define the per-group weight matrix W(γ) ∈ RN ×Pγ restricted to the particles in Πγ and the nodes in Nγ . The per-group mass matrix is ij,(γ)
Mab
(γ)
(γ) ij = Wap s̃p Wbp ,
(22)
and the full cell contribution to the global mass matrix is the sum over all groups: G−1 X ij,(γ) ij M̂ab += Mab , (23) γ=0
where each Mij,(γ) is deposited to the (generally distinct) global nodes in Nγ . 11
An alternative is to embed all particles’ weights into a single vector of S length | γ Nγ |, padding with zeros for unsupported nodes, and form one large outer product. This is mathematically correct since zero weights remove cross-terms, but computationally wasteful, as the accumulator grows from N 2 S 2 to | γ Nγ | with many structurally zero entries. The number of support groups G depends on the interpolation order n: • First-order (CIC): G = 1. (n+1)d = 2d nodes.
All particles in a cell share the same
• Second-order (TSC): G ≤ nd = 2d in d dimensions. Each particle’s stencil can be shifted by one node per dimension relative to the cell corner. • In general, for order n in d dimensions: G ≤ nd . The complete mass matrix assembly at the cell level thus has a two-level structure: 1. Support groups (γ = 0, . . . , G−1): partition particles by the set of N grid nodes they touch. 2. Batches of Kt : within each support group, particles are further partitioned into batches of size Kt for the MMA tile operation. After each support group is processed, the accumulated tile(s) are deposited to the global mass matrix at the addresses determined by the group’s node set Nγ . 3.4. Mass matrix sparse stencil deposition ij Once the cell-local (or per-group) mass matrix M̂ab has been accumulated, it must be scattered into the global mass matrix, which is stored in a compact sparse format indexed by canonical stencil offsets. For an order-n shape function in d dimensions, the displacement between any two nodes in a particle’s support ranges over {−n, . . . , +n}d , giving ij ji (2n+1)d possible offsets. By exploiting the symmetry Mgg ′ = Mg ′ g , which reduces to Mgg′ = Mg′ g in the scalar case, only the "forward half" plus the diagonal need be stored. Concretely, in the case of first and second order interpolation functions, we have: • CIC (n = 1): displacements in {−1, 0, +1}d , giving 3d offsets, of which (3d +1)/2 are canonical. 12
p
p
(a) Π0 (lower left)
(b) Π1 (lower right)
p
p
(c) Π2 (upper left)
(d) Π3 (upper right)
Figure 5: Four possible two-dimensional TSC supports inside a fixed cell. The orange cell is the particle-containing cell, the dashed lines mark its half-cell boundaries, and the blue dots are the 3 × 3 interpolation nodes used for the particle position shown in each panel. The four panels correspond to the lower-left, lower-right, upper-left, and upperright particle classes within the cell. Particles that fall in the same panel share the same node set and therefore belong to the same group.
• TSC (n = 2): displacements in {−2, . . . , +2}d , giving 5d offsets, of which (5d +1)/2 are canonical. When support groups are present (G > 1), different groups deposit to different, but overlapping, sets of global nodes. The canonical stencil index for a given local entry (a, b) therefore depends on the group’s base offset. A precomputed lookup table can be used to map each local node pair to the corresponding canonical stencil index and global node address, ensuring correct assembly irrespective of the number of groups. 3.5. Algorithmic summary Algorithm 1 describes the general cell-level procedure for assembling the mass matrix using hardware matrix engines with generic MMA tiles of shape 13
(Mt , Nt , Kt ), with G support groups per cell and Nc tensor components per node pair. For first-order shape functions (n = 1), only a single support group exists (G = 1) and the outer loop is trivial. For higher-order shape functions (n ≥ 2), the support-group loop executes up to G = nd iterations. Within each group, the batch loop processes all assigned particles in chunks of Kt . Algorithm 1 Cell-local mass matrix assembly via tensor cores. 1: Input: Cell particle list {p0 , . . . , pP −1 }, grid geometry, interpolation or-
der n, tile shape (Mt , Nt , Kt ).
2: Output: Contributions accumulated into global mass matrix Mij .
▷ Nodes per particle support ▷ Padded node count ▷ Tile grid
3: N ← (n+1)d 4: Npad ← ⌈N/Mt ⌉ · Mt 5: Tr ← Npad /Mt , Tc ← Npad /Nt
6: for each support group γ ∈ {0, . . . , G−1} do 7: Determine node set Nγ and deposit addresses. ij ← 0, r ∈ 8: Initialize accumulator tiles: Drc
{0, . . . , Tr −1}, c ∈
{0, . . . , Tc −1}. for each particle p ∈ Πγ do Compute shape-function weights Wap , a ∈ {0, . . . , N −1}. Compute per-particle coefficients s̃ij p. ij Buffer Wap and s̃p into current batch. if batch full (Kt particles accumulated) then (β) ij,(β) 14: Form Aak = Wa,pβ,k s̃ij pβ,k and Bkb = Wb,pβ,k . 15: for each tile (r, c) and each component (i, j) do ij,(β) (β) ij ij 16: MMA: Drc ← Drc + Ar Bc 17: end for 18: end if 19: end for 20: Process remaining partial batch with zero-padded MMA. ij 21: Deposit Drc to global Mij via stencil lookup. 22: end for
9: 10: 11: 12: 13:
14
4. Numerical results To demonstrate the practical viability of the mass matrix matrix-product reformulation, we specialize the general framework of Section 3 to first-order (CIC) and second-order (TSC) B-spline interpolation in three dimensions and implement it on NVIDIA GPUs using the Warp Matrix Multiply-Accumulate (WMMA) intrinsics provided by the CUDA programming model. We target two tile formats: • FP64 tile (Mt , Nt , Kt ) = (8, 8, 4): all operands in double precision. The 8 × 8 accumulator matches the CIC support size exactly (N = 8), so a single tile covers the full outer product with batch size Kt = 4. • TF32 tile (Mt , Nt , Kt ) = (16, 16, 8): inputs in TF32 (10-bit mantissa, 8-bit exponent) with FP32 accumulation. For TSC (N = 27), the weight matrix is padded to 32 rows and the 32 × 32 accumulator is covered by 2 × 2 = 4 tiles, with batch size Kt = 8. Exploiting the spatial symmetry M̂ab = M̂ba (Eq. (11)), tile (1, 0) is skipped and only the three upper-triangle tiles are computed. Both kernels assign one warp per cell with a grid-stride loop over cells. Particles are processed in bounded chunks (up to Pmax = 64) to control register pressure, and accumulator fragments persist across chunks to reduce atomic reductions in main memory. For CIC, 8 × Kt -sized batches are assembled via warp shuffles. For TSC, fragment data is staged through shared memory to assemble the 16 × Kt sub-tiles. Table 1 summarizes the mapping from interpolation order to WMMA tile parameters in the experiments. Firstly, we assess the performance of MMA mass matrix assembly in isolation, comparing the execution times of the WMMA kernels with those of conventional, highly optimized GPU implementations. Then, we assess the benefit provided by MMA in a production PIC simulation, using the ECSIM algorithm. We run all the isolation experiments on a single node machine, equipped with an AMD EPYC 7302P 32-core CPU, and an NVIDIA A100 GPU. The production PIC simulations are run on a multi-node cluster equipped with 2x AMD Rome 7H12 CPUs and 4x NVIDIA A100 GPUs per node. 4.1. WMMA comparison against optimized conventional GPU kernels The isolation experiments measure the mass matrix assembly kernel on a 3D domain of Nx × Ny × Nz cells with a single species and a uniform 15
Table 1: Mapping of cell-local mass matrix dimensions to MMA tile parameters for the NVIDIA tensor core implementation. N is the number of nodes in a particle’s support, Mt , Nt , and Kt are the tile size, with Kt corresponding to the batch size, and G is the number of support groups per cell. P added is the matrix size after rounding up, while T iles represents the number of tiles required to cover the matrix of size P added.
Interpolation
Precision
N
Padded
Tile (Mt , Nt , Kt )
Tiles
G
CIC (n = 1) TSC (n = 2)
FP64 TF32/FP32
8 27
8 32
(8, 8, 4) (16, 16, 8)
1 3†
1 ≤8
†
Upper-triangle tiles only: (0,0), (0,1), (1,1). Tile (1,0) is skipped by symmetry.
number of particles per cell (ppc), pre-sorted by cell. We test both CIC and TSC interpolation for both scalar and 3 × 3 ECSIM tensorial mass matrices (Eq. 4). CIC experiments use FP64 data format, while TSC experiments use FP32 data with TF32 inputs and FP32 accumulation. Two parameter investigations are performed: (i) varying ppc with a fixed 16 × 16 × 16 grid, and (ii) varying the grid size with a fixed 128 ppc. Figure 6 reports the CIC scalar mass matrix results in FP64. The WMMA kernel incurs no penalty even at 1 ppc, and its speedup grows monotonically with particle density: from 1.2× at 13 ppc to 3.7× at 1024 ppc. Larger domains also benefit more, with speedups exceeding 2× across all tested grid sizes. Figure 7 shows the corresponding CIC results for the 3 × 3 ECSIM tensorial mass matrix. The trend is similar, with speedups exceeding 2× above 64 ppc and reaching 2.6× at 1024 ppc. The lower acceleration relative to the scalar case is due to the additional non-MMA work in the tensorial kernel, namely loading magnetic field values and precomputing the rotation tensor αp (Eq. (5)). These stages share the same implementation in both kernels and dilute the tensor-core advantage. Figures 8 and 9 report the corresponding TSC results in FP32. The overall trend mirrors the CIC case: tensor cores are increasingly beneficial at higher particle densities, with peak speedups at 1024 ppc of 2.2× (scalar) and ∼ 2× (tensorial). The speedups are lower than in the CIC case for two reasons. First, as in the tensorial CIC case, non-MMA pre-computation stages reduce the advantage. Second, particles are sorted by cell but not by support group. We process four of the eight TSC groups per pass, requiring two full passes and thus loading each particle from main memory twice. This overhead is common to both WMMA and non-WMMA kernels but 16
83 123 Grid size
Particles per cell
1.00x
NVIDIA w/o TC 1 NVIDIA-TC 1.20x 13 1.85x 32 2.38x 64 2.81x 128 3.12x 256 3.49x 512 3.70x 1024 10 1 100 Time [ms]
NVIDIA w/o TC NVIDIA-TC
2.12x 2.44x 2.82x
163
3.32x
243
3.55x
323
3.59x
483 10 1
(a)
100
Time [ms]
(b)
Figure 6: Tensor cores kernel speedup with respect to a conventional GPU kernel for calculating the scalar mass matrix with CIC interpolation in FP64 precision. a) Performance varying the number of particles per cell, with a fixed grid of 16 × 16 × 16 cells; b) Performance varying the grid size with a fixed number of 128 particles per cell.
NVIDIA w/o TC NVIDIA-TC
83 123 Grid size
Particles per cell
1.00x
1 1.08x 13 1.53x 32 2.10x 64 2.36x 128 256 512 1024 10 1
2.49x
2.05x 2.38x
163
2.49x
243
2.40x
323
2.59x
2.44x
483
2.63x
Time [ms]
NVIDIA w/o TC NVIDIA-TC
2.09x
100
10 1
(a)
Time [ms]
100
(b)
Figure 7: Tensor cores kernel speedup with respect to a conventional GPU kernel for calculating the 3 × 3 ECSIM tensorial mass matrix with CIC interpolation in FP64 precision. a) Performance varying the number of particles per cell, with a fixed grid of 16 × 16 × 16 cells; b) Performance varying the grid size with a fixed number of 128 particles per cell.
inflates the total cost, reducing the relative gain from tensor cores (peak scalar speedup 2.2× vs. 3.7× for CIC, tensorial ∼ 2× vs. 2.6×). With TSC, WMMA shows a slight disadvantage at low particle density (< 32 ppc) or in small domains. The 2.6× speedup at 8 × 8 × 8 in Figure 8 b) is an outlier caused by the conventional kernel underperforming, with such a small domain 17
0.46x
NVIDIA w/o TC NVIDIA-TC
0.85x
83
1.39x 1.73x 1.98x
NVIDIA w/o TC NVIDIA-TC
2.60x
123
1.09x
Grid size
Particles per cell
1 13 32 64 128 256 512 1024
1.87x 1.71x
163
1.70x
243
1.75x
323
2.15x
1.77x
483
2.23x
101
100 Time [ms]
100
(a)
Time [ms]
101
(b)
1 13 32 64 128 256 512 1024
0.43x
NVIDIA w/o TC NVIDIA-TC
0.66x
123
0.99x 1.27x 1.55x 1.81x
1.55x 1.61x
243
1.62x 1.61x
483
2.08x
101 Time [ms]
163
1.46x
323
1.97x
100
NVIDIA w/o TC NVIDIA-TC
83 0.75x Grid size
Particles per cell
Figure 8: Tensor cores kernel speedup with respect to a conventional GPU kernel for calculating the scalar mass matrix with TSC interpolation in FP32 precision. a) Performance varying the number of particles per cell, with a fixed grid of 16 × 16 × 16 cells; b) Performance varying the grid size with a fixed number of 128 particles per cell.
102
101
(a)
Time [ms]
102
(b)
Figure 9: Tensor cores kernel speedup with respect to a conventional GPU kernel for calculating the 3 × 3 ECSIM tensorial mass matrix with TSC interpolation in FP32 precision. a) Performance varying the number of particles per cell, with a fixed grid of 16 × 16 × 16 cells; b) Performance varying the grid size with a fixed number of 128 particles per cell.
the GPU is heavily underutilized, and the WMMA kernel better hides the control-flow latency of TSC deposition.
18
4.2. End-to-end acceleration of a kinetic plasma simulation We assess the end-to-end impact of tensor cores in a production 3D double Harris current sheath magnetic reconnection simulation (Figure 11) using the ECSIM algorithm with full periodic boundary conditions. The domain is 28 × 14 × 14 di discretized on a 160 × 80 × 80 grid, with two species (ions q/m = 1, electrons q/m = −256) at 768 ppc per species (∼ 1.5×109 particles total). Since the methodology requires particles sorted by cell, we compare three ECSIM pipeline variants to isolate the contributions of sorting and tensor cores: (i) unsorted particles with a naive atomic-based mass matrix kernel; (ii) sorted particles with an optimized conventional GPU kernel; (iii) sorted particles with WMMA mass matrix assembly. In variants (ii) and (iii), sorted particle layout is exploited in all particle-related kernels. The only difference between them is the use of tensor cores in mass matrix deposition. All three variants use CIC interpolation with FP64 precision and FP64 (8, 8, 4) tiles, and run with four MPI processes on four NVIDIA A100 GPUs. Figure 10 reports the per-cycle time breakdown, averaged over 300 steps. Because the ECSIM pipeline overlaps MPI communication, host-device transfers, and computation via task-based parallelism, we group the cycle into three non-overlapping stages: deposition (mass matrix assembly and moment deposition), sort & communicate (particle sorting and MPI particle exchange), and other (field solver, particle mover, diagnostics). I/O is excluded. In variant (i), each step takes ∼ 2, 800 ms, with deposition accounting for more than 90%. Introducing sorting (ii) reduces deposition from ∼ 2, 500 ms to ∼ 200 ms at the cost of increasing sort & communicate from ∼ 150 ms to ∼ 240 ms, a strongly favorable trade-off. Tensor cores (iii) further reduce the deposition time by ∼ 40% relative to variant (ii). Profiling with NVIDIA Nsight Systems shows that the mass matrix kernel alone runs in ∼ 24 ms with WMMA versus ∼ 61 ms without, a ∼ 2.5× speedup consistent with the isolation results at 768 ppc (Figure 7). Overall, combining sorting with tensor cores yields a ∼ 5.8× end-to-end speedup over the unsorted baseline. To verify physical correctness, we run the WMMA variant for 600 ωpi (Figure 11) and compare the total-energy evolution across all three variants over 200 ωpi (Figure 12). The plot shows the signed difference E(t) − Eu (0), where Eu (0) is the initial energy of the unsorted variant. All three implementations preserve energy to machine precision, confirming their physical equivalence. 19
2521.4 ms 90.5%
Unsorted Sorted w/o NVIDIA-TC Sorted NVIDIA-TC 0
2787.3 ms
198.2 ms 36.2% 547.8 ms
Deposition Sort & Communicate Other Total avg +/- std 1000 1500 2000 2500 3000 Average Time Per Cycle [ms]
123.5 ms 25.9% 476.9 ms
500
Figure 10: Execution time breakdown of a single PIC cycle, averaged over 300 time steps, for a naive ECSIM implementation without particle sorting, a partially optimized implementation with particle sorting without WMMA mass matrix deposition, and a fully optimized implementation that leverages particle sorting and WMMA mass matrix deposition. Colors identify the time spent in the current deposition and mass matrix assembly (orange); particle sorting and communication (blue); particle mover, field solver and diagnostics (yellow). The simulations are run with four MPI processes on 4x NVIDIA A100 GPUs.
5. Discussion Kernel-level speedup. Matrix engines consistently accelerate mass matrix assembly for both CIC and TSC interpolation, and for both scalar and tensorial mass matrices. The magnitude of the speedup depends on the fraction of kernel instructions that can execute on tensor cores: • Lighter kernels (e.g., CIC scalar) benefit most, as particle data are loaded once and no preliminary computation is required, yielding up to 3.7× acceleration. • Heavier kernels (e.g., TSC tensorial) include non-MMA stages, such as magnetic field loading and rotation-tensor precomputation, that reduce the tensor-core advantage; nevertheless, speedups of ≳ 1.5× are still achieved. End-to-end impact. The ∼ 2.5× kernel-level acceleration observed in isolation is reproduced in the production simulation at 768 ppc. However, the end-to-end speedup depends on the fraction of the PIC cycle spent in mass matrix assembly. In our 3D multi-GPU ECSIM simulation, the optimized 20
Figure 11: Charge density of ion species and reconnecting field lines in a 3D double Harris current sheath magnetic reconnection simulation, at t = 600 ωpi . The simulation is run accelerating the ECSIM algorithm with tensor cores for the mass matrix assembly.
1.50
1e 16
E(t) Eu(0)
1.25 1.00 0.75 0.50 0.25 0.00
0
Unsorted Sorted w/o NVIDIA-TC Sorted NVIDIA-TC 50 100 Time [ pit]
150
200
Figure 12: Exact energy conservation to machine precision for the three ECSIM implementations (unsorted, sorted w/o NVIDIA TC, sorted with NVIDIA TC). The plot shows the signed difference between the total system energy E(t) and the total system energy at time t = 0 measured in the unsorted pipeline, Eu (0).
pipeline (without tensor cores) already reduces the deposition stage to ∼ 36% of the cycle, limiting the tensor-core benefit to ∼ 15% end-to-end. In configurations where deposition dominates, such as 2D simulations without domain decomposition (not reported here), we observed end-to-end speedups of up 21
to 30%. Sorting particles by cell is a prerequisite for the MMA-based assembly and incurs a cost of up to ∼ 20% of each cycle, but this is more than compensated by the deposition speedup. Exactness and precision. The reformulation is mathematically exact, no approximations are introduced beyond the floating-point rounding inherent in a given tile precision. In our experiments, no rounding error was observed with the FP64 (8, 8, 4) tile, while a relative error of ∼ 10−4 was measured with the TF32/FP32 (16, 16, 8) tile, consistent with the reduced TF32 mantissa. Mixed-precision workflows are natively supported, since data can be cast on the fly to match the MMA input format. In the production simulation of Section 4.2, for instance, field values and particle positions are stored in FP32, while velocities, statistical weights, and the mass matrix use FP64. Tile-support matching and portability. The approach is most effective when the tile dimensions match the interpolation support. The (8, 8, 4) tile fits CIC (N = 8) exactly, and the (16, 16, 8) tile is well suited to TSC (N = 27 padded to 32). Higher-order interpolation functions generally benefit from larger tiles; when substantial padding is required, correctness is preserved but throughput is reduced. The methodology is formulated in terms of a generic MMA abstraction and is directly portable to any platform exposing tile-level MMA instructions. For example, AMD matrix cores offer FP64 (16, 16, 4) and FP32 (32, 32, 2) tiles, the latter covers the 27-node TSC support in a single tile without padding. 6. Conclusion In this work, we showed that the mass matrix assembly arising in implicit PIC methods can be reformulated exactly as a tensor contraction that maps onto hardware matrix-multiply-accumulate units. The key observation is that the mass matrix is inherently sparse due to the compact support of the interpolation functions, and the weighted outer-product accumulation over particles decomposes, cell by cell, into a sequence of matrix products whose inner dimension coincides with the contraction dimension of hardware MMA tiles. We developed a complete algorithmic framework that includes: i) a celllocal factorization of the global tensor contraction into two-operand matrix products, ii) a particle-batching scheme that partitions particles into groups 22
of size Kt matching the MMA tile inner dimension, and iii) a supportgroup decomposition that handles the position-dependent stencil placement of higher-order shape functions. The resulting formulation is general with respect to interpolation order, spatial dimension, and the scalar or tensorial nature of the mass matrix, while remaining independent of the specific hardware platform. We specialized the framework to first-order (CIC) and second-order (TSC) B-spline interpolation in three dimensions and implemented on NVIDIA tensor cores using WMMA intrinsics. Isolation benchmarks on an NVIDIA A100 GPU demonstrated speedups of up to 3.7× for CIC scalar and 2.6× for CIC mass matrices, and up to 2.2× and ∼ 2× for the corresponding TSC cases, with the acceleration being more appreciable at high particle densities. In a production 3D magnetic reconnection simulation using the ECSIM algorithm with CIC interpolation, the tensor-core-accelerated mass matrix kernel achieved a ∼ 2.5× speedup over the optimized conventional GPU kernel, translating into a ∼ 15% reduction of the end-to-end wall-clock time per PIC cycle. Importantly, the proposed reformulation is exact and introduces no approximations beyond the floating-point rounding inherent in a given tile precision. The methodology is expressed in terms of a generic MMA abstraction and is therefore directly portable to AMD matrix cores, Intel AMX, Google TPUs, and any future architecture exposing a tile-level MMA interface. More broadly, the same tensor-contraction approach extends to any particleto-grid scatter operation, where one MMA operand encodes the deposited quantities and the other encodes the interpolation weights. The effectiveness depends on tile occupancy, standard charge and current deposition (four quantities per particle in 3D) would leave an (8, 8, 4) tile half-occupied, whereas algorithms that scatter higher-order moments, such as the Implicit Moment Method [25, 26, 27] (ten quantities per particle) or high-order moment diagnostic calculations, can fill the tile more efficiently. 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. 23
References [1] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, A. Borchers, R. Boyle, P.-l. Cantin, C. Chao, C. Clark, J. Coriell, M. Daley, M. Dau, J. Dean, B. Gelb, T. V. Ghaemmaghami, R. Gottipati, W. Gulland, R. Hagmann, C. R. Ho, D. Hogberg, J. Hu, R. Hundt, D. Hurt, J. Ibarz, A. Jaffey, A. Jaworski, A. Kaplan, H. Khaitan, D. Killebrew, A. Koch, N. Kumar, S. Lacy, J. Laudon, J. Law, D. Le, C. Leary, Z. Liu, K. Lucke, A. Lundin, G. MacKean, A. Maggiore, M. Mahony, K. Miller, R. Nagarajan, R. Narayanaswami, R. Ni, K. Nix, T. Norrie, M. Omernick, N. Penukonda, A. Phelps, J. Ross, M. Ross, A. Salek, E. Samadiani, C. Severn, G. Sizikov, M. Snelham, J. Souter, D. Steinberg, A. Swing, M. Tan, G. Thorson, B. Tian, H. Toma, E. Tuttle, V. Vasudevan, R. Walter, W. Wang, E. Wilcox, D. H. Yoon, In-datacenter performance analysis of a tensor processing unit, SIGARCH Comput. Archit. News 45 (2) (2017) 1–12. doi:10.1145/3140659.3080246. URL https://doi.org/10.1145/3140659.3080246 [2] P. Micikevicius, S. Narang, J. Alben, G. Diamos, E. Elsen, D. Garcia, B. Ginsburg, M. Houston, O. Kuchaev, G. Venkatesh, H. Wu, Mixed precision training (10 2017). doi:10.48550/arXiv.1710.03740. [3] A. Reuther, P. Michaleas, M. Jones, V. Gadepally, S. Samsi, J. Kepner, Survey of machine learning accelerators, in: 2020 IEEE High Performance Extreme Computing Conference (HPEC), 2020, pp. 1–12. doi:10.1109/HPEC43674.2020.9286149. [4] N. Jouppi, G. Kurian, S. Li, P. Ma, R. Nagarajan, L. Nai, N. Patil, S. Subramanian, A. Swing, B. Towles, C. Young, X. Zhou, Z. Zhou, D. Patterson, Tpu v4: An optically reconfigurable supercomputer for machine learning with hardware support for embeddings (04 2023). doi: 10.48550/arXiv.2304.01433. [5] N. J. Higham, T. Mary, Mixed precision algorithms in numerical linear algebra, Acta Numerica 31 (2022) 347–414. doi:10.1017/ S0962492922000022. [6] A. Haidar, S. Tomov, J. Dongarra, N. J. Higham, Harnessing gpu tensor cores for fast fp16 arithmetic to speed up mixed-precision iterative 24
refinement solvers, in: Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis, SC ’18, IEEE Press, 2019. doi:10.1109/SC.2018.00050. URL https://doi.org/10.1109/SC.2018.00050 [7] C. Cui, Acceleration of tensor-product operations with tensor cores, ACM Trans. Parallel Comput. 11 (4) (Nov. 2024). doi:10.1145/ 3695466. URL https://doi.org/10.1145/3695466 [8] X. Liu, Y. Liu, H. Yang, J. Liao, M. Li, Z. Luan, D. Qian, Toward accelerated stencil computation by adapting tensor core unit on gpu, in: Proceedings of the 36th ACM International Conference on Supercomputing, ICS ’22, Association for Computing Machinery, New York, NY, USA, 2022. doi:10.1145/3524059.3532392. URL https://doi.org/10.1145/3524059.3532392 [9] L. Oostrum, B. Veenboer, R. Rook, M. Brown, P. Kruizinga, J. W. Romein, The Tensor-Core Beamformer: A High-Speed Signal-Processing Library for Multidisciplinary Use , in: 2025 IEEE International Parallel and Distributed Processing Symposium (IPDPS), IEEE Computer Society, Los Alamitos, CA, USA, 2025, pp. 582–592. doi:10.1109/IPDPS64566.2025.00058. URL https://doi.ieeecomputersociety.org/10.1109/ IPDPS64566.2025.00058 [10] G. Schieffer, I. Peng, Accelerating drug discovery in autodock-gpu with tensor cores, in: J. Cano, M. D. Dikaiakos, G. A. Papadopoulos, M. Pericàs, R. Sakellariou (Eds.), Euro-Par 2023: Parallel Processing, Springer Nature Switzerland, Cham, 2023, pp. 608–622. [11] J. Finkelstein, J. S. Smith, S. M. Mniszewski, K. Barros, C. F. A. Negre, E. H. Rubensson, A. M. N. Niklasson, Quantum-based molecular dynamics simulations using tensor cores, Journal of Chemical Theory and Computation 17 (10) (2021) 6180–6192. doi:10.1021/acs.jctc. 1c00726. URL https://doi.org/10.1021/acs.jctc.1c00726 [12] G. Lapenta, Exactly energy conserving semi-implicit particle in cell formulation, Journal of Computational Physics 334 (2017) 349–366. 25
doi:https://doi.org/10.1016/j.jcp.2017.01.002. URL https://www.sciencedirect.com/science/article/pii/ S0021999117300128 [13] D. Burgess, D. Sulsky, J. Brackbill, Mass matrix formulation of the flip particle-in-cell method, Journal of Computational Physics 103 (1) (1992) 1–15. doi:https://doi.org/10.1016/0021-9991(92)90323-Q. URL https://www.sciencedirect.com/science/article/pii/ 002199919290323Q [14] D. Sulsky, Z. Chen, H. Schreyer, A particle method for historydependent materials, Computer Methods in Applied Mechanics and Engineering 118 (1) (1994) 179–196. doi:https: //doi.org/10.1016/0045-7825(94)90112-0. URL https://www.sciencedirect.com/science/article/pii/ 0045782594901120 [15] D. Sulsky, S.-J. Zhou, H. L. Schreyer, Application of a particlein-cell method to solid mechanics, Computer Physics Communications 87 (1) (1995) 236–252, particle Simulation Methods. doi:https://doi.org/10.1016/0010-4655(94)00170-7. URL https://www.sciencedirect.com/science/article/pii/ 0010465594001707 [16] C. K. Birdsall, A. B. Langdon, Plasma Physics via Computer Simulation, 1991. [17] R. Hockney, Computer Simulation Using Particles, CRC Press, 1988. URL https://books.google.se/books?id=SVslEAAAQBAJ [18] T. Montoya, D. W. Zingg, A unifying algebraic framework for discontinuous galerkin and flux reconstruction methods based on the summationby-parts property, Journal of Scientific Computing 92 (3) (2022) 87. doi:10.1007/s10915-022-01935-3. URL https://doi.org/10.1007/s10915-022-01935-3 [19] B. Perse, K. Kormann, E. Sonnendrücker, Geometric particle-in-cell simulations of the vlasov–maxwell system in curvilinear coordinates, SIAM Journal on Scientific Computing 43 (1) (2021) B194–B218. arXiv:
26
https://doi.org/10.1137/20M1311934, doi:10.1137/20M1311934. URL https://doi.org/10.1137/20M1311934 [20] J. Monaghan, Particle methods for hydrodynamics, Computer Physics Reports 3 (2) (1985) 71–124. doi:https: //doi.org/10.1016/0167-7977(85)90010-3. URL https://www.sciencedirect.com/science/article/pii/ 0167797785900103 [21] S. Markidis, S. W. D. Chien, E. Laure, I. B. Peng, J. S. Vetter, Nvidia tensor core programmability, performance & precision, in: 2018 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), 2018, pp. 522–531. doi:10.1109/IPDPSW.2018. 00091. [22] G. Schieffer, D. Medeiros, J. Faj, A. Marathe, I. Peng, On the rise of amd matrix cores: Performance, power efficiency, and programmability, 2024, pp. 132–143. doi:10.1109/ISPASS61541.2024.00022. [23] H. Kim, G. Ye, N. Wang, A. Yazdanbakhsh, N. S. Kim, Exploiting intel advanced matrix extensions (amx) for large language model inference, IEEE Comput. Archit. Lett. 23 (1) (2024) 117–120. doi:10.1109/LCA. 2024.3397747. URL https://doi.org/10.1109/LCA.2024.3397747 [24] K. Bowers, Accelerating a particle-in-cell simulation using a hybrid counting sort, Journal of Computational Physics 173 (2) (2001) 393– 411. [25] J. Brackbill, D. Forslund, An implicit method for electromagnetic plasma simulation in two dimensions, Journal of Computational Physics 46 (2) (1982) 271–308. doi:https: //doi.org/10.1016/0021-9991(82)90016-X. URL https://www.sciencedirect.com/science/article/pii/ 002199918290016X [26] S. Markidis, G. Lapenta, Rizwan-uddin, Multi-scale simulations of plasma with ipic3d, Mathematics and Computers in Simulation 80 (7) (2010) 1509–1519, multiscale modeling of moving interfaces in materials. doi:https://doi.org/10.1016/j.matcom.2009.08.038. 27
URL https://www.sciencedirect.com/science/article/pii/ S0378475409002444 [27] S. Markidis, P. Henri, G. Lapenta, K. Rönnmark, M. Hamrin, Z. Meliani, E. Laure, The fluid-kinetic particle-in-cell method for plasma simulations, Journal of Computational Physics 271 (2014) 415–429.
28