arXiv CSopen access

TokaGLINT: A Scalable GPU-Tailored Implicit Solver for Full 3D Tokamak Electromagnetic Simulations

TokaGLINT: A Scalable GPU-Tailored Implicit Solver for Full 3D Tokamak Electromagnetic Simulations Zifan Yang1,2 , Haoyuan Zhang1 , Jialin Li3 , Wu Yuan1 , Xiazhen Liu1 , Jian Zhang1,* , Jianyuan Xiao4,* , Shan Liang1,*

arXiv:2609.21366v1 [cs.DC] 18 Sep 2026

1

Computer Network Information Center, Chinese Academy of Sciences, Beijing, China 2 University of Chinese Academy of Sciences, Beijing, China 3 Tsinghua University, Beijing, China 4 University of Science and Technology of China, Hefei, China

Abstract—We introduce TokaGLINT, a GPU-accelerated implicit solver for electromagnetic field computations in full 3D tokamak simulations, aimed at efficient large-scale parallel GPU computing. Its central innovation lies in the co-design of hierarchical domain decomposition and a fast exact local solver, where hierarchical partitioning is tailored to match fine-grained intracard subdomains and exploit the tensor-based solver dedicated to curvilinear-coordinate symplectic CN-FDTD-discretized 3D Maxwell equations. Backed by automated operator fusion and batching customized for the intra-card multi-subdomain structure, the solver decouples unknowns through discrete transformations and leverages tensor-structured computations to achieve high hardware utilization, while preserving the long-time stability characteristic of symplectic discretizations. TokaGLINT scales the electromagnetic field solve beyond 10,000 GPUs, achieving 90.1% weak and 53.9% strong scaling efficiency, while delivering a 2.67× single-node speedup over an unpreconditioned BiCGStab baseline (HIP-enabled HYPRE). It is validated in EAST tokamak simulations within the SymPIC plasma simulation code, enabling high-fidelity long-duration modeling. Index Terms—High performance computing, Linear systems, Partitioning algorithms, Plasma simulation.

I. I NTRODUCTION When studying wave heating, current drive, and energeticparticle interactions with the background plasma in tokamaks [1]–[3], an electromagnetic (EM) fully kinetic model is desirable since it can capture the non-Maxwellian velocity distributions, finite Larmor radius effects, or wave-particle resonances that are central to these processes. The particle-in-cell (PIC) method is the standard computational approach for simulating EM fully kinetic plasmas. In this method, macro-particles sample the 6D distribution function, Maxwell’s equations govern the evolution of EM fields, which in turn exert Lorentz forces on charged plasma particles to update their motion state, while the spatial distribution and motion of plasma particles (i.e., charge and current densities) serve as the source terms of Maxwell’s equations, forming a closed self-consistent loop. * Corresponding authors ([email protected]; [email protected]; [email protected]).

Fig. 1. 3D view of the electric field excited by a point source simulated with the CN-FDTD scheme accelerated by TokaGLINT. The simulation domain size and wave frequency are referenced from the lower hybrid wave injection experiment conducted on the EAST Tokamak.

The wave-heating, current-drive, and energetic-particle processes of interest evolve on time scales far longer than the fundamental periods that limit explicit PIC time steps, i.e., the inverse plasma frequency, the inverse cyclotron frequency, and the EM wave crossing time across a grid cell (the EM Courant-Friedrichs-Lewy). As a result, simulations routinely span millions of time steps. Under these conditions, cumulative discretization errors in physical invariants (energy, charge, or more general, symplectic structure) must remain bounded, a requirement that places strict demands on the numerical scheme. In response, a family of structure-preserving geometric PIC algorithms has been developed [4]–[9], whose spacetime discretizations preserve key physical symmetries exactly in a discrete sense. The SymPIC [10]–[12] code, a member of this family, features an explicit symplectic schemes and has been applied to full 3D tokamak simulations at extreme scale [13]. However, explicit time integration is constrained by the Courant-Friedrichs-Lewy (CFL) condition. In the low-density plasma regimes targeted here, or in simulations employing a

This paper has been accepted at the International Conference for High Performance Computing, Networking, Storage, and Analysis (SC 2026).

reduced ion-to-electron mass ratio, the EM wave CFL is the binding constraint, impractically small time steps are required, making long-time scale simulations computationally inefficient and prohibitively costly. Tokamak geometry is toroidal, naturally calling for a cylindrical or toroidal coordinate mesh. SymPIC already operates on a cylindrical mesh, the challenge is to upgrade its EM solver to an implicit formulation that preserves both the symplectic structure and compatibility with the curvilinear discretization. The Crank-Nicolson (CN) time discretization is the most natural choice for this upgrade. Within SymPIC’s discrete action principle, the EM part of the Lagrangian action integral is discretized in time. Applying the mid-point rule to this integral yields, via the discrete Euler-Lagrange equations, exactly the CN-FDTD formulation of Maxwell’s equations [14]–[16]. This choice is not unique, other implicit timeintegration methods can also remove the EM CFL, but each carries practical drawbacks. The widely used alternatingdirection implicit (ADI) FDTD method [17], [18] achieves unconditional stability by splitting each time step into directional substeps, but it sacrifices the discrete symplectic structure and suffers from more severe numerical dispersion than the CN scheme at large time steps [19], [20]. A symplectic variant of the backward Euler method can be constructed, yet it is only first-order accurate, losing the second-order precision essential for long-time scale integrations. The mid-point/CN scheme is the simplest extension that simultaneously removes the EM CFL, preserves the symplectic structure exactly, retains charge conservation, and integrates seamlessly with SymPIC’s existing spatial discretization and current deposition framework, all at the cost of a single linear solve per time step. In this work, we upgrade SymPIC to the implicit CNFDTD framework, achieving both high numerical fidelity and improved computational efficiency while retaining excellent long-term stability. A brief derivation of charge conservation for the field-implicit SymPIC scheme will be given in the supplement material. The implementation of CN-FDTD entails the large-scale parallel solution of curl-curl linear systems. The primary challenge arises from the frequent neighbor and global communication operations inherent in iterative methods for sparse linear systems, which severely hinder parallel scalability. Domain decomposition preconditioning serves as an effective solution to address this scalability challenge. A well-designed preconditioner can substantially reduce the number of iterations while lowering the proportion of communication overhead, thereby improving parallel scalability. For the overall performance, the cost associated with preconditioning operations must be offset by the reduction in iteration count and the enhancement in scalability. In general, high-precision preconditioners are more effective in reducing the number of iterations but typically incur higher computational overhead. Thus, there exists a tradeoff between the cost and strength of the preconditioner. The exceptional computing power of GPUs offers ample room for the design and optimization of high-strength preconditioners. Motivated by these facts, we propose TokaGLINT (GPU

Linear Iterative Solver with Numerical Transforms for Tokamak), a composite preconditioner based on multi-level domain decomposition and a tensor-form fast solver for subdomains, which is specifically tailored for the large-scale GPU-parallel solution of CN-FDTD. A key feature of this composite preconditioner is its ability to accommodate multiple overlapping subdomains within a single accelerator card. This capability allows us to employ a high-computational-intensity tensorform solver for each subdomain, while the overlap between subdomains can be efficiently implemented using High Bandwidth Memory (HBM). When combined with the upper-level traditional overlapping domain decomposition across different cards, this composite preconditioner substantially reduces the number of outer iterative method iterations compared to conventional preconditioners. Specifically, the superior convergence performance of the proposed scheme is quantitatively validated in Section III-D, while its enhanced computational efficiency is demonstrated in Section V. A further key advantage of the composite preconditioner is that the subdomain solver is specifically designed to preserve the tensor structure of the curl-curl differential operator and the EM field. We develop a fast algorithm that achieves partial decoupling of the unknowns via discrete transformations derived from the differential operator and its eigenvalue/singularvalue decomposition. While this approach has been demonstrated in the literature [21]–[24], extending it to efficiently handle the complex curl-curl operator with variable coefficients introduced by cylindrical coordinates poses a highly nontrivial challenge. Nevertheless, we successfully generalize and adapt this approach to the tokamak simulation regime with substantial and practically impactful speedup, owing not only to the co-design of hierarchical domain decomposition framework and fast exact subdomain solver but also to the careful engineering of operator fusion and batching across multiple intra-card subdomains. Our primary contributions are outlined as follows. • We propose a novel GPU-tailored Hierarchical Additive Schwarz Method (HASM) preconditioner. Built around the co-design of intra-card domain decomposition and physics-driven fast subdomain solvers, this architecture is natively optimized to deliver efficient large-scale parallel GPU solutions for full 3D tokamak EM simulations. • We develop a fast exact local solver tailored to the linear systems arising from the curvilinear-coordinate symplectic CN-FDTD-discretized 3D Maxwell equations. It leverages discrete transformations to decouple unknowns, exploits tensor-structured operations for accelerated computation, and is seamlessly integrated with the intracard multi-subdomain framework to enable automated operator fusion and batching, consistently achieving high hardware utilization across diverse simulation scenarios. • The proposed TokaGLINT delivers outstanding performance, sustaining 90.1% weak and 53.9% strong scaling efficiency across more than 10,000 GPUs, together with a 2.67× single-node speedup over an unpreconditioned BiCGStab baseline (HIP-enabled HYPRE). Integrated

into SymPIC, TokaGLINT has been validated in full-3D EAST tokamak EM simulations, confirming its practicality for real-world fusion plasma simulation. II. R ELATED WORK A. EM solver in PIC PIC codes are widely recognized as essential tools for plasma simulation, and numerous research teams worldwide have dedicated efforts to advancing this technology. A variety of high-quality full-kinetic PIC software packages have been developed, with representative examples including SMILEI [25], PICADOR [26], EPOCH [27], VPIC 2.0 [28], WarpX [29], iPIC3D [30], [31], PIConGPU [32], [33], OSIRIS [34], and SymPIC [13]. iPIC3D is an implicit PIC code with its field solver running on CPUs. All other packages adopt explicit numerical schemes and exhibit remarkable scalability when deployed on large-scale computing platforms. For instance, VPIC 2.0 has successfully completed weak scaling tests involving over 17,000 GPUs [28], while WarpX has been optimized for exascale supercomputers such as Frontier, Fugaku, and Summit, leveraging thousands of AMD and NVIDIA GPUs [29]. SymPIC [13] scaled up to over 40,000,000 cores on the new Sunway supercomputer. These explicit PIC codes have been widely applied in plasma physics research, covering areas including laser-plasma interaction, magnetic reconnection, and inertial confinement fusion. Their excellent scalability originates from the local nature of particle updates. However, they are tightly constrained by the CFL condition when plasma density is low or a reduced ion-electron mass-ratio is used. A complementary strategy is implicit PIC [35]–[38], which treats both field and particle advances implicitly. The direct implicit approach [35], [36] linearizes the coupled system about the previous time step, requiring only a single linear solve per time step and avoiding particle-orbit Jacobian evaluations. In contrast, the fully implicit Newton-Krylov approach [37], [38] solves the coupled nonlinear system via Newton-Krylov iteration, incurring multiple linear solves and particle-orbit Jacobian evaluations per time step, with irregular communication patterns that severely degrade GPU scalability. Implicit PIC becomes necessary when the plasma-frequency CFL (ωpe ∆t < 2) is the binding constraint. In the lowdensity or reduced mass-ratio regimes targeted in this work, however, the EM CFL dominates, and the CN-FDTD fieldimplicit strategy pursued here captures the majority of the practical speedup at a fraction of the algorithmic complexity. B. Multilevel domain decomposition Domain decomposition methods [39]–[42] are widely recognized as effective solution approaches suitable for parallel computations. Among these techniques, the two-level overlapping Schwarz domain decomposition methods [42]–[44] incorporate a coarse space correction mechanism to mitigate the convergence degradation of one-level methods when dealing with a large number of subdomains. Although they exhibit superior convergence speed compared to one-level schemes, they often require an accurate global solution at the coarse level,

which ultimately becomes a computational bottleneck. Distinct from the two-level approaches, hierarchical partitioning here in this work serves solely to distinguish between inter-card and intra-card interactions, while the second level is primarily designed to cooperate with batched subdomain solvers based on tensor operations and resolve the data exchange issues among these subdomains. Concurrent with our research, Ichitaro Yamazaki et al. [45] also adopt the strategy of assigning multiple subdomains to each GPU to improve accelerator utilization. Their approach achieves parallel computation of subdomains on the same GPU through NVIDIA MPS. As noted in [45], running multiple MPI processes on a single GPU may not be the most optimal approach; however, achieving the same decomposition by having multiple subdomains per MPI process would necessitate substantial innovations in algorithms and extensive software development efforts. TokaGLINT follows this latter route: one MPI process manages multiple subdomains on a GPU. This design exposes opportunities for fusing and batching subdomain solves on the same GPU. Detailed discussions are provided in Sec. IV.B. At the subdomain solver level, the family of batched solvers developed by Anzt et al. [46]–[50] is primarily designed for scenarios requiring the batch inversion and factorization of thousands of small matrices, such as the general purpose block-Jacobi and ISAI preconditioning, where the involved matrix dimensions are typically modest (≤ 100 × 100). In this study, we integrate the batched solver into the overlapping additive Schwarz method (ASM) preconditioning framework, whereas the equivalent matrix size is considerably larger (> 3000×3000) but the number of matrices is drastically fewer (typically only a few dozen). By exploiting the symmetry of the curl-curl differential operator, we design a discrete transform to partially decouple EM field unknowns, allowing each solver to efficiently handle matrices of this scale. C. Discrete transform based solver A compact representation of the difference operator was first proposed by Frederic in [51] to achieve storage compression. The authors of [21] further extended the compact representation to three-dimensional problems via a tensor-based formulation and subsequently proposed fast exact solvers for partial differential equations. The approach has since been generalized to a broader class of PDEs, as demonstrated in [22], [23], [52], among others. Within this line of research, the method was adopted in [24] for solving complex curl-curl systems in Cartesian coordinates. These discrete-transformbased direct solvers achieve high efficiency via unknown decoupling. FlashMP [24] leverages the symmetry of the curlcurl operator defined on Yee grids to decouple unknowns over separate grid points. However, cylindrical meshes break this symmetry, making FlashMP inapplicable in this setting. In the present work, the method is extended to curvilinear coordinates for the first time, and is efficiently coupled with the aforementioned HASM preconditioner, which adopts the co-design of numerical algorithms and GPU implementations

to achieve a favorable balance among convergence, computational complexity and GPU throughput, enabling efficient solving on more than 10,000 GPUs.

For each point c, define qc = rc /r0 and the local diagonal blocks Pc = diag(∆x, ∆y, ∆z), Rc = diag(1, qc , 1), Wc = diag(qc , 1, qc ),

III. A LGORITHM

(5)

Qc = diag(qc−1 , qc , qc−1 ) = Rc Wc−1 .

A. Background SymPIC [10]–[13] constructs explicit charge-conservative symplectic PIC algorithms on cylindrical meshes. In the present work, we revise the temporal discretization of the EM field action integral to a midpoint, equivalently CrankNicolson, form. Coupled with SymPIC’s existing particle push and charge-conservative current deposition, which remain unchanged from the explicit formulation. This yields a fieldimplicit, charge-conservative, symplectic PIC scheme on the cylindrical mesh. Charge conservation follows directly from the discrete gauge invariance of the action [10], [11], and a brief derivation is provided in the supplement material. The linear system associated with this field-implicit scheme is solved via our TokaGLINT framework. We start from the EM field equations (Faraday’s and Ampère’s laws)  ∂E   = ∇ × B − J, ∂t ∂B   = −∇ × E, ∂t

(1)

where E, B, and J denote the electric field, magnetic field, and electric current density, respectively. Applying the CN-FDTD scheme to the above system and eliminating B gives the vector-field curl-curl system   2 fw I + ∆t4 curlbw E n+1 = b, (2) d curld where b contains known fields, source terms, boundary contributions, and particle-current contributions. The operators bw curlfw d and curld denote the forward and backward discrete curl operators associated with staggered discrete field layout. The SymPIC field representation is geometric: the electric unknown is a discrete 1-form [11], [12] rather than the vector E itself. Accordingly, the vector-field system above must be rewritten in the 1-form unknown used by SymPIC. The conversion from vector components to 1-form components is determined by the metric of the coordinate system. For the tokamak geometry considered in this work, we adopt the cylindrical-coordinate formulation used in SymPIC [12], [13]. The standard cylindrical coordinates (r, θ, z) are mapped to the straightened coordinates (x, y, z) as x = r − r0 ,

y = r0 θ,

z = z,

(3)

where r0 is a fixed radial reference length and the axial coordinate is unchanged. For typical tokamak applications, r0 is the major radius, and θ denotes the toroidal angle. the line-element vector is ∆s =  Under this mapping,  ∆x, ∆y rr0 , ∆z . It relates the conventional vector field to the 1-form field E1-form = P RE. (4)

The global diagonal matrices P , R, W , and Q are assembled from these local blocks using the point-wise ordering (x, y, z)c1 , (x, y, z)c2 , . . .. The linear system obtained from symplectic CN-FDTD discretization is then given by 2 fw n+1
n+1 R−1 P −1 E1-form E1-form + ∆t4 P R curlbw = P Rb. (6) d curld

To expose the tensor-product stencil used by the fast subdomain solver, we separate the logical-coordinate spacings from the coordinate metric weights and introduce the Cartesian-like ˜ d . Its forward form is logical-grid curl curl  Fy,i,j,k+1 − Fy,i,j,k Fz,i,j+1,k − Fz,i,j,k −   ∆y ∆z      Fz,i+1,j,k − Fz,i,j,k  Fx,i,j,k+1 − Fx,i,j,k fw ˜  . (7)  curld Fi,j,k =  −  ∆z ∆x      Fy,i+1,j,k − Fy,i,j,k Fx,i,j+1,k − Fx,i,j,k  − ∆x ∆y 

Its relation to the original curvilinear discrete curl curld can be written as
˜ d F = W curld R−1 F . curl (8) We further introduce the scaled electric unknown Ẽ = P −1 E1-form .

(9)

Since P , R, W and Q are diagonal under the pointwise ordering, their products and inverses reduce to local component-wise scalings within each point. Combining (8) with (9) and applying straightforward algebra, (6) becomes   2 ˜ fw Ẽ n+1 = Rb. ˜ bw Q curl (10) Ẽ n+1 + ∆t4 Q curl d d This metric-weighted curl-curl system is the linear system addressed by TokaGLINT. B. Pipelined BiCGStab with HASM preconditioning The outer Krylov subspace iteration in TokaGLINT is realized using a preconditioned pipelined BiCGStab algorithm [53], whose full procedure is described in Algorithm 1. In contrast to standard preconditioning setups, we introduce a customized HASM preconditioner. Its key merit lies in a GPU-oriented hierarchical decomposition, which allows a single GPU to handle multiple overlapping subdomains simultaneously, combined with strong, computationally intensive subdomain solvers. The detailed derivation and formulas of the subdomain solvers are presented in Sec. III.C, while the hierarchical domain decomposition and its efficient GPU implementation are described in Sec. IV. In the pipelined BiCGStab framework, the global reductions are overlapped with SpMV and preconditioner applications. This strategy is most effective when the preconditioner application is compute-bound and requires only limited neighbor

Algorithm 1 Preconditioned Pipelined BiCGStab

can be written in the compact tensor form

Input: Matrix, RHS, Preconditioner Output: Solution 1: Initialization 2: loop 3: computation local axpby, dot-product partial sums 4: Begin global sum allreduce 1 5: Apply HASM preconditioner 6: Apply SpMV with halo-exchange overlap 7: End global sum allreduce 1 8: computation local axpby, dot-product partial sums 9: Begin global sum allreduce 2 10: Apply HASM preconditioner 11: Apply SpMV with halo-exchange overlap 12: End global sum allreduce 2 13: computation local axpby 14: end loop

x ∆x G = Df w ⃝G, where

communication [37], both of which match the design of our HASM preconditioner: the subdomains are compute-intensive, and inter-subdomain communication within a single GPU is realized via HBM. With this framework, the bottleneck caused by global reductions and overlapping-region communications in large-scale parallel execution can be substantially mitigated. This effect is evaluated experimentally in Sec. V. Moreover, the intra-GPU domain decomposition also plays an important role in complexity control and supports operator fusion across subdomains, which we discuss in detail in Sec. IV. C. Fast subdomain solver We adopt the compact operator notation introduced in [21]. For any matrix T = {Tij } ∈ Rn×n and 3D field 3 x ⃝, y ⃝ z are q = {qijk } ∈ Rn , a set of compact operators ⃝, defined as follows, x = Tim qmjk , T ⃝q y = Tjm qimk , T ⃝q

(11)

z = Tkm qijm . T ⃝q It is easy to verify that for any T1 , T2 ∈ Rn×n , x 2 ⃝q x = (T1 T2 )⃝q, x T1 ⃝T x + T2 ⃝q x = (T1 + T2 )⃝q, x T1 ⃝q x 2 ⃝q y = T2 ⃝T y 1 ⃝q. x T1 ⃝T Recall (10) and notice that both curl operators (7) are composed of one-dimensional forward and backward finite differences. The forward difference of field G

−1



1

  Df w =  

−1

 ..

.

..

.

   ∈ Rn×n , 

1 −1

and Dbw = −DfTw . The backward differences, as well as differences in all other coordinate directions can be represented similarly. In (5), qc = rc /r0 varies only along the x direction. We −1 −1 define H = diag(q111 , q211 , . . . , qn−1 ) ∈ Rnx ×nx , where nx , x 11 ny , and nz denote the lengths in the x, y, and z directions, respectively. For simplicity, we set ∆x = ∆y = ∆z = 1. Unless otherwise specified, the 3D field vectors in this subsection are stored in field-major order, whereas Sec. III.A uses point-major ordering; no additional notation is introduced to distinguish these two orderings. From (5) and (7), we then obtain 

z y y z − H ⃝D x f w ⃝e x f w ⃝e H ⃝D



  −1 ˜ fw  Qcurl x z x f w ⃝e z x − H −1 ⃝D x f w ⃝e ⃝D d E = H ,

(12)

y x x y − H ⃝D x f w ⃝e x f w ⃝e H ⃝D

where E = (ex , ey , ez )T . This gives the alternative form (15) of (10), where the unknown is E. For simplicity, both sides have been multiplied by β = 4/∆t2 , and the resulting scalar factor and the right-handside matrix R have been absorbed into b = (bx , by , bz )T . Let U, S, V be the singular value decomposition of Df w , Df w = U SV T and Dbw = −V SU T . Exploring the symmetry of the double-curl system, we can define discrete transforms in tensor form to decouple the unknowns in (10). However, in cylindrical coordinates, this decoupling is only partial, as the metric Q appears between the two curl operators. Define a forward transformation (mapping the electric field vector and right-hand side vector to the feature space) y T ⃝b z x hx = V T ⃝V y T ⃝b z y hy = U T ⃝V T

(13)

T

y z z hz = V ⃝U ⃝b The corresponding backward transformation (mapping from the feature space back to the physical space) is defined as: y ⃝f z x ex = V ⃝V y ⃝f z y ey = U ⃝V

(14) y ⃝f z z ez = V ⃝U where fx , fy , fz are the vectors to be solved in the feature (∆x G)i,j,k = Gi+1,j,k − Gi,j,k space. The transforms essentially correspond to dense matrixmatrix multiplications in the size of ni × ni and ni × (nj nk ).
2 x bw Df w ⃝ y − Dbw Df w ⃝ z ex + H 2 Df w ⃝D x bw ⃝e y y + Dbw ⃝D z f w ⃝e x z = bx , βI − H ⃝D
−1 −1 x f w ⃝e y x + βI − Dbw Df w ⃝ z − H Dbw HDf w ⃝ x ey + Dbw ⃝D z f w ⃝e y z = by , (15) H Dbw H ⃝D
x f w ⃝e z x + H 2 ⃝D x bw ⃝D y f w ⃝e z y + βI − HDbw H −1 Df w ⃝ x − H 2 ⃝D x bw Df w ⃝ y ez = bz . HDbw H −1 ⃝D

in these two directions. For example, expanding the first block row in (17) at point (i, j, k) yields   r2 r2 r2 β + s2j r02 + s2k fx,i,j,k − sj r02 fy,i+1,j,k + sj r02 fy,i,j,k i i i −sk fz,i+1,j,k + sk fz,i,j,k = hx,i,j,k . The coupling in the remaining direction is represented by the terms involving H1 , H2 and Df w . Since H is diagonal, the upper-left block of the field-major x 2⃝ y + S 2 ⃝, z is diagonal. matrix M , namely βI + H 2 ⃝S Representing M in block-matrix form, (17) becomes      D11 A12 A13 fx hx  A21 A22 A23  fy  = hy  . (18) A31 A32 A33 fz hz We introduce the following auxiliary variables: fx′ = fx , Fig. 2. Fast exact subdomain solver. The original physical-space system (a) is transformed into the feature-space system (b), where the discrete transforms partially decouple the unknowns. Using the block structure, the system is algebraically reduced to two equivalent subsystems: a diagonal subsystem (c) and a two-component subsystem coupled along x-lines (d). After permuting (d) from field-major to point-major ordering, the coupled subsystem becomes block diagonal (e). Solving (c) and (e), followed by recombination and the backward transform, gives the solution (f) of the original system.

We first illustrate the forward transformation using the ey term in the first equation of (15) as an example. This term x bw ⃝e y y . According to (13), the corresponding is H 2 Df w ⃝D y T⃝ z for the first equaforward transform operator is V T ⃝V tion. This simplifies to y T ⃝H z 2 Df w ⃝D x bw ⃝e y y V T ⃝V T

T

2

(16)

Then (18) becomes  D11  A21 A31

D11 A22 A−1 12 D11 A32 A−1 12 D11

 ′   D11 fx hx −1 A23 A13 D11  fy′  = hy  . (20) hz fz′ A33 A−1 13 D11

Next, we introduce other auxiliary variables: fx′′ = fx′ + fy′ + fz′ , fy′′ = D11 fy′ ,

(21)

fz′′ = D11 fz′ .

x ⃝ z −Df w ⃝S fw z − H 1 Df w ⃝ x y ⃝ z βI + S 2 ⃝ −S ⃝S x ⃝S y ⃝ z x + H 2 ⃝S x 2⃝ y −H 2 ⃝S βI − H2 Df w ⃝

y T ⃝b z x, h′x = hx = V T ⃝V −1 x −1 U T ⃝V y T ⃝b z y, h′y = A−1 21 hy = H1 ⃝S

By applying the same derivation and simplification to all terms, the full system (15) is transformed into      fx hx  M  fy  = hy  , (17) fz hz

x ⃝ y H1 ⃝S x ⃝ z H2 ⃝S

(19)

−1 fz′ = D11 A13 fz .

and new forward transformation T

y z x SU ⃝U y ⃝V y ⃝f z y = − V ⃝V ⃝H Df w ⃝V x ⃝f y y. = − H 2 Df w ⃝S

where M takes the following form  βI + H 2 ⃝S x 2⃝ y + S2 ⃝ z x ⃝ y −H 2 D ⃝S

−1 fy′ = D11 A12 fy ,



and H1 = H −1 Dbw H, H2 = HDbw H −1 . Note that H x direction; If nx , ny , and H −1 are applied only to the ⃝ and nz differ, Df w and Dbw should be replaced by directiondependent matrices of the corresponding sizes, with U , S and V adjusted accordingly. Per Definition (11), the compact operator only generates node dependencies on the field q along the operating direction. When diagonal matrices T are applied, it introduces no crossnode dependencies. Notice that in M , the linear system in the feature space, only the diagonal singular value matrix y and ⃝ z directions. From the above S is applied to the ⃝ observation, we can easily verify that the system is decoupled

(22)

−1 x T ⃝S y −1 U T ⃝b z z. h′z = A−1 31 hz = H2 ⃝V

Then, by left-multiplying the second and third equations in −1 (20) by A−1 21 and A31 , respectively, we obtain D11 fx′′ = h′x ,  −1 −1 A21 A22 A−1 12 − D11 −1 −1 −1 A A32 A12 − D11  31′  hy − fx′′ = . ′ ′′ hz − fx

−1 −1 A−1 21 A23 A13 − D11 −1 −1 −1 A31 A33 A13 − D11

(23a)   ′′  fy fz′′ (23b)

Since D11 is diagonal, fx′′ in (23a) can be solved directly. Moreover, all A∗ matrices and their inverses are decoupled in the y and z directions. Therefore, after permutation, the 2 × 2 block matrix of (23b) becomes a block-diagonal matrix with ny nz dense blocks, each of size (2nx ) × (2nx ). By precomputing the inverses of these (2nx ) × (2nx ) blocks during initialization, fy′′ and fz′′ can be obtained efficiently. The corresponding storage cost is 4n2x ny nz , which scales as O(N 4/3 ) for balanced subdomains with nx ∼ ny ∼ nz and N = nx ny nz . By combining (14), (19), and (21) with the definitions of A12 and A13 , the solution of the linear system is obtained

x ⃝, y and ⃝ z operations in this through (24). The composed ⃝, expression define the new backward transformation.
−1 ′′ −1 ′′ y ⃝ z fx′′ − D11 ex = V ⃝V fy − D11 fz , −2 y ⃝f z y′′ , x S −1 ⃝V ey = −Df−1 ⃝U wH

(24)

x ⃝U y S −1 ⃝f z z′′ . ez = −Df−1 w ⃝V Figure 2 illustrates the workflow of this fast exact subdomain solver. D. Preconditioning effect In the following tables, L1 and L2 denote the first-level inter-GPU and second-level intra-GPU ASM subdomains, respectively. The notation “count” gives the number of subdomains, and “size” gives the grid size of each subdomain; both are written as a × b × c in the x, y, and z directions. The parameter l denotes the overlap width used in HASM. TABLE I C OMPARISON OF B I CGS TAB I TERATIVE C ONVERGENCE S TEPS U NDER D IFFERENT P RECONDITIONING S ETTINGS . PARAMETERS :tolr = 10−12 , ∆x = 1.1, ∆y = 1.4, ∆z = 1.0, r0 = 192, L1 COUNT =2 × 2 × 2, L1 SIZE = 32 × 32 × 32; T OKAGLINT: L2 COUNT =2 × 2 × 2. MAX: R EACHED MAXIMUM ITERATIONS WITHOUT CONVERGENCE . DIV: S OLUTION DIVERGED DURING COMPUTATION . ∆t

1.0

2.0

4.0

8.0

16.0

Baseline NON PRE

16

34

60

102

167

JACOBI (PETSc) SOR (PETSc) ISAI (HYPRE) ILUT (HYPRE) AMG (HYPRE) TokaGLINT’s HASM(l = 1) TokaGLINT’s HASM(l = 2) TokaGLINT’s HASM(l = 3) TokaGLINT’s HASM(l = 4)

16 12 41 15 5 3 2 2 2

34 25 MAX 31 9 5 3 3 2

61 51 24 115 18 8 6 4 4

96 82 46 MAX 37 12 9 7 7

174 148 100 DIV 63 19 14 11 10

To assess the extent to which preconditioning improves convergence speed, we conduct comparative evaluations against conventional preconditioning techniques. All general-purpose preconditioners are realized through interfaces to the PETSc and HYPRE libraries, with their default parameter settings. Table I reports the BiCGStab iteration counts under different ∆t. TokaGLINT’s HASM requires fewer iterations than the general-purpose preconditioners, especially for large time steps. The results also show improved convergence as the overlap size increases. However, enlarging the overlap region also raises computational overhead and data movement costs (both inter-GPU and intra-GPU). In terms of overall runtime, this creates a trade-off between convergence and per-iteration cost. In most cases, an overlap size of 3 delivers the optimal overall performance. Table II reports the weak-scaling convergence as the total problem size grows proportionally with the number of accelerators. As the parallel scale is increased, the number of BiCGStab iterations remains nearly constant at approximately 7-8 steps. The results demonstrate robust and scalable convergence behavior for large-scale parallel simulations, which

provides a solid basis for applying the proposed method to ultra-large-scale EM problems. TABLE II W EAK - SCALING CONVERGENCE OF T OKAGLINT. PARAMETERS : tolr = 10−12 , ∆x = 1.1, ∆y = 1.4, ∆z = 1.0, ∆t = 8.0, r0 = 1920, L1 SIZE = 128 × 128 × 128, l = 3, L2 COUNT = 4 × 4 × 4. ACCELERATORS EQUALS THE PRODUCT OF THE L1 COUNT.

Accelerators

L1 count

BiCGStab Steps

8 64 512 4096

2×2×2 4×4×4 8×8×8 16 × 16 × 16

7.40 7.90 7.75 8.00

IV. T OKAGLINT FRAMEWORK While the outer Krylov iteration of the TokaGLINT solver is detailed in Section III-B, this section focuses on the HASM preconditioner and the batched multi-subdomain solver implemented on GPU accelerators. A. Hierarchical ASM Figure 3(a) illustrates the two-level hierarchical overlapping domain decomposition. While the actual simulation is fully 3D, a 2D schematic is shown for illustrative purposes, and the detailed overlap regions between neighboring subdomains are omitted for clarity. A sector of the global domain is first partitioned into first-level (L1) overlapping subdomains, each assigned to a GPU. Each first-level subdomain is then further decomposed into smaller overlapping sub-subdomains (L2). The L2 subdomains at same radial layers in the cylindrical coordinate system share the same metric tensors i.e., the transforms defined in (22) and (24), and the diagonal solve (23a) and (23b), and they are grouped into one bin to facilitate operator fusion and batching. The streaming between different bins is illustrated in Figure 3(b). As derived in Sec. III.C and illustrated in Figure 3(c), two types of expensive operations in the subdomain solve can be identified: the block-diagonal solve (23b) with precomputed inverse blocks, which has memory footprint and computational complexity both O(n4 ), where n (assuming nx = ny = nz = n) denotes the subdomain side length; the forward/backward transforms (22) and (24), with computational complexity O(n4 ). In a single subdomain solve, the transforms are in the form of transposed matrix multiplication, and the block diagonal solve is matrix-vector multiplication. The former is batched within each bin, and the latter is fused into matrix multiplication. The fusion of the block diagonal solve (23b) is the most desirable feature for the L2 domain decomposition and geometry-driven binning. Within a bin, the block solves at the same diagonal-block position can be fused into one GEMM. Consequently, the block-diagonal solves step across L2 subdomains in this bin can be executed as batched GEMM, which is much more efficient than launching multiple separate matrixvector multiplications. Furthermore, this fusion requires a

Fig. 3. Schematic diagram of HASM for parallel EM simulation. (a) Two-level domain decomposition: L1 subdomains follow the inter-GPU decomposition, with each assigned to one MPI process and one GPU; each L1 subdomain is further split into L2 subdomains and automatically grouped into bins sharing geometry-dependency of the subdomain solver. (b) Stream-level execution of an L1 solve: The bins adopt identical color scheme as in (a). (c) Main operators and data layouts: The green geometry-dependent transformation and block diagonal matrices are shared across all L2 subdomains within each bin. Data exclusive to the L2 subdomains in one bin are color-coded with alternative colors. Left: transforms on field-major variables are batched across subdomains. Right: block-diagonal solves on point-major vectors, fused together and realized by batched GEMM. Note that the fusion operation relies on the custom data layout illustrated at the bottom, where lij refers to line segments within each L2 subdomain, with i and j corresponding to indices along the y and z direction, respectively.

customized data layout for L2 subdomains within each bin. All of the above designs are illustrated in the right half of Figure 3(c). The transforms (22) and (24) natively adopt GEMM form. Nevertheless, small subdomain size n limits single-subdomain performance. By binning subdomains together, we execute them as batched GEMM operations, which coarsens the operation granularity. The reason these operations are not fused together, even though they share identical metric-dependent transform matrices, is that successive transformations along different directions require transposed matrices multiplication. Fusing multiple such operations would demand highly customized low-level implementations, and the performance gains cannot offset the associated development overhead. Binning further enables fusion of additional operations including permutation and L1 communication buffer packing, and so on. Most of them are adopted in our implementation, but since they contribute minor performance gains, we do not elaborate further here. To support this grouping strategy efficiently, the solver integrates an automated mechanism for meta-data reorganization and mask mapping, handling memory management and boundary condition adaptation. B. Discussion We further discuss alternative strategies corresponding to different features of TokaGLINT. Large single-domain solve per GPU. Different from the large single-domain solve per GPU adopted by many conventional methods, including FlashMP [24], our work introduces

L2 domain decomposition and geometry-driven binning strategies, effectively reducing the overall computational complexity and fuses, batches, and pipelines matrix transforms, blockdiagonal solves, and inter-subdomain data movement across multiple sub-solvers, which large single-domain solves cannot achieve. We construct a representative example, whose results are presented in Table III. When solving with a single large subdomain yields poor efficiency (the bottom row), the L2 domain decomposition substantially boosts solving efficiency, with iteration counts staying consistent. The one iteration discrepancy observed in the case for relative tolerance 10−10 can be attributed to different convergence paths taken to reach the prescribed tolerance. TABLE III R EPRESENTATIVE B ENCHMARKS FOR L2 D OMAIN D ECOMPOSITION . E ACH ENTRY IS REPORTED AS T /I , WHERE T IS THE TOTAL SOLVE TIME IN MILLISECONDS AND I IS THE NUMBER OF ITERATIONS TO CONVERGENCE . PARAMETERS :∆x = ∆y = ∆z = 1.0, ∆t = 10.0, r0 = 1152, l = 3, L1 SIZE = 64 × 64 × 64, L1 COUNT = 2 × 2 × 2. Relative tolerance L2 size

L2 count

10−8

10−10

10−12

16 × 16 × 16 32 × 32 × 32 64 × 64 × 64

4×4×4 2×2×2 1×1×1

49.8/6 63.0/6 209.0/6

63.8/8 71.7/7 238.9/7

74.6/9 92.1/9 302.7/9

Sparse direct solver. Regarding computational cost, the transformation and block-diagonal solve, both exhibit O(n4 ) complexity, which is comparable to the solve-phase cost of sparse direct solvers such as a general purpose sparse

LU solve. However, sparse triangular solves exhibit inherent data dependencies, and suffer from indirect memory accesses, irregular execution, and unstructured memory access patterns. By contrast, empowered by operator fusion realized through geometry-driven binning of L2 subdomains, the dominant computation within TokaGLINT’s solving stage is restructured into batches of concurrent, uniform-sized GEMM operations. Why skip the MPS-based approach in [45]? A key design advantage of TokaGLINT comes from fusing the O(n4 ) blockdiagonal solves across multiple L2 subdomains, converting many independent matrix-vector multiplications into batched GEMMs. This scheme relies on geometry-driven binning and custom-designed data layouts, which is difficult to achieve without explicit inter-process control of bin partitioning and synchronization. For this reason, the single-process-per-GPU deployment becomes a natural choice.

TABLE IV P ERFORMANCE COMPARISON WITH HIP- ENABLED HYPRE. E ACH ENTRY IS REPORTED AS T /I , WHERE T IS THE TOTAL SOLVE TIME IN MILLISECONDS AND I IS THE NUMBER OF ITERATIONS TO CONVERGENCE . PARAMETERS : tolr = 10−12 , ∆x = ∆y = ∆z = 1.0, ∆t = 8.0, r0 = 1920, L1 SIZE = 64 × 64 × 64; T OKAGLINT: l = 3, L2 COUNT = 2 × 2 × 2.

HYPRE TokaGLINT Speedup

np=8

np=64

np=256

np=512

212.38/125 79.67/8

224.02/122 81.96/8

235.10/122 85.09/8

257.04/122 84.79/8

2.67×

2.73×

2.76×

3.03×

V. E XPERIMENTS All experiments are conducted on China’s latest heterogeneous supercomputer supporting both scale-up and scale-out capabilities. Each compute node is equipped with 8 GPUs and two 64-bit CPUs. Each CPU operates at 2.4 GHz with 64 cores, adopts a NUMA-based memory organization, supports eight-channel DDR5-6400 memory, and connects to accelerators via PCIe Gen5. Each GPU integrates 320 SIMD units, with a theoretical peak double-precision (FP64) performance of 32.7 TFLOPS. The GPU is further equipped with 64 GB HBM and supports a theoretical peak memory bandwidth of 1.8 TB/s. Within each node, accelerators are interconnected by a high-speed intra-node accelerator link. Inter-node cluster networking relies on 4×400 Gbps InfiniBand-like, RDMAcapable links. The implementation uses HIP-compatible GPU kernels and GPU-aware MPI for inter-GPU communication. The software stack consists of a GPGPU programming environment compatible with mainstream GPGPU API standards, Clang 17.0.0, and Open MPI 5.0.3. In the experiments presented in this section, we use overlap layers l = 3 for the reason that preliminary tuning showed the best time-to-solution in most cases. The L1 sizes 643 –1283 per GPU, are representative of PIC workloads, consistent with VPIC 2.0 using 1003 cells per GPU V100 [28] and WarpX reporting tens of millions of cells per GPU in FOM tests [29]. These workloads provide substantial local computation to amortize communication and expose GPU batching efficiency, which should be considered when interpreting the reported scaling efficiency. A. Performance We compare TokaGLINT against the BiCGStab solver in HIP-enabled HYPRE 2.32. We use unpreconditioned BiCGStab in HYPRE as the baseline. As shown in Table IV, TokaGLINT achieves a speedup of 2.67× on a single node (8 accelerators), and the speedup ratio further increases as the problem scale expands. This scaling trend is fully consistent with the expected high parallel scalability of the TokaGLINT solver.

Fig. 4. Weak-scaling runtime breakdown of TokaGLINT. Stacked bars show the average time of one linear-system solve to convergence. The two numbers above each bar denote the average solve time in seconds and the parallel efficiency, respectively. The dashed line marks the single-GPU reference case: L1 domain decomposition is absent, so no cross-GPU communication takes place.

B. Parallel Scalability To assess the scalability of TokaGLINT, we carry out weakand strong-scaling studies. All reported times in Fig. 4 and Fig. 5 are averaged over 100 complete linear-system solves to convergence. For the weak scaling test (Table V), we keep the number of unknowns per accelerator fixed. We use the 16-accelerator case (2 nodes) as the scaling baseline because a single-GPU run has no L1 decomposition or inter-GPU communication, and a single-node run mainly exercises direct GPU communication rather than representative multi-node behavior. As shown in Fig. 4, TokaGLINT achieves an excellent parallel efficiency of 90.1% when scaling from 16 to 10,000 accelerators (625×). These results compellingly confirm TokaGLINT’s robust scalability for ultra-large-scale EM simulations. For the strong-scaling test, the total problem size is fixed, as listed in Table VI. TokaGLINT reaches the relative tolerance in 6 iterations for all tested parallel configurations. As shown in Fig. 5, when scaling from 672 to 10,752 accelerators, TokaGLINT maintains a parallel efficiency of 53.9%. Thanks to the communication-hiding BiCGStab solver, global Allreduce operations are overlapped with preconditioning and SpMV, effectively reducing communication overhead. As shown in Fig. 5 (right), the fractions of SpMV and ASM

Fig. 5. Strong-scaling runtime breakdown of TokaGLINT. The left panel shows the absolute average time of one linear-system solve to convergence, with the two numbers above each bar denoting the average solve time in seconds and the parallel efficiency. The right panel shows the normalized time distribution, with the number above each bar denoting the communication percentage of the total solve time.

TABLE V W EAK - SCALING TEST SETUP. PARAMETERS : tolr = 10−8 , ∆x = ∆y = ∆z = 1.0, ∆t = 8.0, r0 = 1920, l = 3. Accelerators

Global grid

L1 count

L2 count

16 128 512 2048 8192 10000

512 × 256 × 256 1024 × 512 × 512 1024 × 1024 × 1024 2048 × 2048 × 1024 4096 × 2048 × 2048 3200 × 2560 × 2560

4×2×2 8×4×4 8×8×8 16 × 16 × 8 32 × 16 × 16 25 × 20 × 20

4×4×4 4×4×4 4×4×4 4×4×4 4×4×4 4×4×4

TABLE VI S TRONG - SCALING TEST SETUP. PARAMETERS : tolr = 10−7 , ∆x = ∆y = ∆z = 1.0, ∆t = 8.0, r0 = 1920, l = 3. Accelerators

Global grid

L1 count

L2 count

672 1344 2688 5376 10752

2048 × 1536 × 1792 2048 × 1536 × 1792 2048 × 1536 × 1792 2048 × 1536 × 1792 2048 × 1536 × 1792

8 × 6 × 14 8 × 12 × 14 16 × 12 × 14 16 × 12 × 28 16 × 24 × 28

8×8×4 8×4×4 4×4×4 8×8×4 8×4×4

and a long-time energy-conservation test. We first test the correctness of the field solver using a wave-excitation case. The simulation domain is nx = ny /4 = nz = 256, r0 = 384∆l, ∆x = ∆z = ∆l, ∆y = 2πr0 ∆l/ny . PEC boundaries are adopted for x and z directions, while periodic boundary condition is chosen for y direction. Wave source is located at x = nx ∆x, y = 102∆y, z = nz ∆z/2 and the frequency is ω = 0.2c/∆l. Time step is set to ∆t = 2∆l/c and total number of time steps is nt = 1600. Evolution of electric field is shown in Fig. 6, which clearly shows the propagation process of the EM wave inside the toroidal domain. For visualization, results computed in the logical cylindrical coordinate system (x, y, z) are mapped to Cartesian coordinates (X, Y, Z) and plotted on the slice Z = 128. The coordinate transformation is defined as follows: X = (i∆x + r0 ) cos(j∆y/r0 ) ,

halo exchange rise with parallel scaling, while Allreduce increases more significantly, though the overall communication ratio remains below 50%. Owing to the high computational cost of subdomain solvers, the algorithm maintains good strong-scaling efficiency. The observed efficiency loss mainly comes from the growing Allreduce proportion and reduced local computation time. At large core counts, lower peraccelerator computational intensity limits the batched local solver from reaching peak performance. Allreduce suffers the largest performance degradation, followed by ASM halo communication, whose data volume is three times that of SpMV halo exchange. C. Simulation We next validate the integration of TokaGLINT into the SymPIC EM field solve through a wave-propagation test

Y = (i∆x + r0 ) sin(j∆y/r0 ) ,

Z=z .

A key merit of adopting symplectic structure-preserving PIC scheme is that the truncation errors of fundamental system invariants such as the total energy remain uniformly bounded over long simulation time intervals. To verify this property, we performed a three-dimensional magnetized electron toroidal plasma simulation with the parameters specified below. λd = 1.85 × 10−3 ∆l, ωpe,0 ∆l/c = 2.82 × 10−1 ,

ωce,0 ∆l/c = 1.69 × 10−1 , B = B0⃗j,

x0 = 1920∆l, ∆x = 1.0∆l,

nx = 2ny = nz = 128,

∆y = 1.2∆l, ∆tc/∆l = 2.5,

∆z = 1.3∆l,

Fig. 6. Wave-propagation validation of the CN-FDTD field solver in SymPIC using TokaGLINT. Wave source is located at (X ≈ 519∆l, Y ≈ 375∆l). The panels show time snapshots of the Ey value on the z-midplane, visualized in Cartesian coordinates after mapping the logical cylindrical grid (x, y, z) to (X, Y, Z).

q 2 n0 q e qe B0 where ωpe,0 = ε0 me , vte , λd = vte /ωpe and ωce = me are the plasma frequency, thermal speed, Debye length and cyclone frequency of electrons, n0 is the reference density, density distribution of electrons is 2

2

ne,i (x, y, z) = n0 e−R /R0

p where R = (x − 64)2 + (z − 64)2 , R0 = 7.1, c is the speed of light in the vacuum. The total number of simulation time steps is nt = 1.0 × 106 , which means ∆tnt ωpe = 7.05 × 105 . where me , qe , vte are mass, charge, thermal speed of electrons, respectively. The external magnetic field is B = B0⃗j, and the number of sampling particles per grid at R = 0 is set to 56. The total energy is recorded every 2000 time steps, which is shown in Fig. 7. As the initial field is far from equilibrium, the energy at the second output step (nt = 2000) is taken as the baseline E0 . The internal energy E is found to vary within ±1‰ of E0 for the rest of the simulation. It is clear that the total energy is well conserved, which verifies the advantage of the present field-implicit symplectic PIC scheme for cylindrical coordinates.

computational bottlenecks of the CN-FDTD discretization of Maxwell’s equations in curvilinear coordinates on modern heterogeneous GPU platforms, as part of the upgrade to the EM solver of the large-scale parallel PIC simulation software SymPIC. Built upon a hierarchical domain decomposition framework coupled with tensor-structured fast subdomain solvers tailored for curvilinear coordinate systems, this composite preconditioner effectively overcomes the scalability challenges of linear system solution on GPU clusters, delivering substantial speedup and excellent parallel scalability while preserving rigorous numerical fidelity for large-scale tokamak simulations. Extensive numerical experiments thoroughly verify the effectiveness and competitiveness of the proposed method. Compared with the widely adopted general-purpose solver library HYPRE, TokaGLINT achieves a performance improvement of 2.67× on a single node. At massively parallel scales, the solver maintains a weak scaling efficiency of 90.1% and a strong scaling efficiency of 53.9% across more than 10,000 GPUs. Moreover, its practical deployment within the SymPIC code demonstrates its capability to enable highly efficient implicit EM field solutions under realistic physical scenarios. Future work will aim to broaden the applicability to more diverse physical scenarios and multiphysics simulation workflows. In particular, we will investigate fully implicit formulations that treat both electromagnetic field evolution and particle kinetics within a consistent implicit framework. ACKNOWLEDGMENTS

Fig. 7. Long-time total-energy variation of a magnetized toroidal electron plasma simulated by the field-implicit symplectic PIC scheme in SymPIC using TokaGLINT.

VI. C ONCLUSION In this work, we propose TokaGLINT, a highly scalable, hardware-aware linear solver designed to address the

We sincerely thank our anonymous SC reviewers for their valuable comments. This work was supported by the National Key Research and Development Program of China (Grant No. 2025YFB3003403) and the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB0500101).

R EFERENCES [1] N. J. Fisch, “Theory of current drive in plasmas,” Reviews of Modern Physics, vol. 59, no. 1, pp. 175–234, 1987. https://doi.org/10.1103/ RevModPhys.59.175 [2] R. Prater, “Heating and current drive by electron cyclotron waves,” Physics of Plasmas, vol. 11, no. 5, pp. 2349–2376, 2004. https://doi. org/10.1063/1.1690762 [3] A. Fasoli, C. Gormenzano, H. L. Berk, B. Breizman, S. Briguglio, D. S. Darrow, N. Gorelenkov, W. W. Heidbrink, A. Jaun, S. V. Konovalov, et al., “Chapter 5: Physics of energetic ions,” Nuclear Fusion, vol. 47, no. 6, pp. S264–S284, 2007. https://doi.org/10.1088/0029-5515/47/6/S05 [4] J. Squire, H. Qin and W. M. Tang, “Geometric integration of the VlasovMaxwell system with a variational particle-in-cell scheme,” Physics of Plasmas, vol. 19, no. 8, pp. 084501, 2012. https://doi.org/10.1063/1. 4742985 [5] Y. He, H. Qin, Y.-J. Sun, J.-Y. Xiao, R.-L. Zhang and J. Liu, “Hamiltonian time integrators for Vlasov-Maxwell equations,” Physics of Plasmas, vol. 22, no. 12, pp. 124503, 2015. https://doi.org/10.1063/1. 4938034 [6] H. Qin, J. Liu, J.-Y. Xiao, R.-L. Zhang, Y. He, Y.-L. Wang, Y.-J. Sun, J. W. Burby, L. Ellison and Y. Zhou, “Canonical symplectic particle-incell method for long-term large-scale simulations of the Vlasov-Maxwell equations,” Nuclear Fusion, vol. 56, no. 1, pp. 014001, 2016. https: //doi.org/10.1088/0029-5515/56/1/014001 [7] Y. He, Y.-J. Sun, H. Qin and J. Liu, “Hamiltonian particle-in-cell methods for Vlasov-Maxwell equations,” Physics of Plasmas, vol. 23, no. 9, pp. 092108, 2016. https://doi.org/10.1063/1.4962573 [8] M. Kraus, K. Kormann, P. J. Morrison and E. Sonnendrücker, “GEMPIC: geometric electromagnetic particle-in-cell methods,” Journal of Plasma Physics, vol. 83, no. 4, pp. 905830401, 2017. https://doi.org/10.1017/ s002237781700040x [9] P. J. Morrison, “Structure and structure-preserving algorithms for plasma physics,” Physics of Plasmas, vol. 24, no. 5, pp. 055502, 2017. https: //doi.org/10.1063/1.4982054 [10] J.-Y. Xiao, H. Qin, J. Liu, Y. He, R.-L. Zhang and Y.-J. Sun, “Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems, ” Physics of Plasmas, vol. 22, Nov 2015, https://doi.org/10.1063/1.4935904. [11] J.-Y. Xiao, H. Qin and J. Liu, “Structure-preserving geometric particlein-cell methods for Vlasov-Maxwell systems,” Plasma Science and Technology, vol. 20, no. 11, Sep 2018, https://doi.org/10.1088/20586272/aac3d1. [12] J.-Y. Xiao and H. Qin, “Explicit structure-preserving geometric particlein-cell algorithm in curvilinear orthogonal coordinate systems and its applications to whole-device 6D kinetic simulations of tokamak physics,” Plasma Science and Technology, vol. 23, no. 5, pp. 055102, 2021, https://doi.org/10.1088/2058-6272/abf125. [13] J.-Y. Xiao et al., “Symplectic Structure-Preserving Particle-in-Cell Whole-Volume Simulation of Tokamak Plasmas to 111.3 Trillion Particles and 25.7 Billion Grids,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’21), St. Louis, MO, USA, 2021, pp. 1–13, https://doi.org/10.1145/3458817.3487398. [14] Y. Yang, R. S. Chen, Edward K. N. Yung, “The unconditionally stable Crank Nicolson FDTD method for three-dimensional Maxwell’s equations,” Microwave and Optical Technology Letters, vol. 48, pp. 1619–1622, 2006. https://doi.org/10.1002/mop.21684 [15] D. M. Sullivan, Electromagnetic Simulation Using the FDTD Method. IEEE Press, 2000, https://doi.org/10.1002/9781118646700 [16] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302–307, May 1966. http://dx.doi.org/10.1109/TAP.1966.1138693 [17] T. Namiki, “A new FDTD algorithm based on alternating-direction implicit method,” IEEE Transactions on Microwave Theory and Techniques, vol. 47, no. 10, pp. 2003–2007, 1999. https://doi.org/10.1109/ 22.795075 [18] F. Zheng, Z. Chen and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 9, pp. 1550–1558, 2000. https://doi.org/10.1109/22.869007 [19] S. G. Garcia, T.-W. Lee and S. C. Hagness, “On the accuracy of the ADI-FDTD method,” IEEE Antennas and Wireless Propagation Letters, vol. 1, pp. 31–34, 2002. https://doi.org/10.1109/LAWP.2002.802583

[20] J. Shibayama, M. Muraki, R. Takahashi, J. Yamauchi and H. Nakano, “Performance evaluation of several implicit FDTD methods for optical waveguide analyses,” Journal of Lightwave Technology, vol. 24, no. 6, pp. 2465–2472, 2006. https://doi.org/10.1109/JLT.2006.874570 [21] Q. Nie, F. Y. M. Wan, Y.-T. Zhang and X.-F. Liu, “Compact integration factor methods in high spatial dimensions,” Journal of Computational Physics, vol. 227, no. 10, pp. 5238–5255, 2008, https://doi.org/10.1016/ j.jcp.2008.01.050. [22] L. Ju, J. Zhang, L. Zhu, and Q. Du, “Fast explicit integration factor methods for semilinear parabolic equations,” Journal of Scientific Computing, vol. 62, pp. 431–455, 2015, https://doi.org/10.1007/s10915-014-9862-9. [23] L. Ju, J. Zhang, and Q. Du, “Fast and accurate algorithms for simulating coarsening dynamics of Cahn-Hilliard equations,”Computational Materials Science, vol. 108, pp. 272–282, 2015, https://doi.org/10.1016/j.commatsci.2015.04.046. [24] H.-Y. Zhang, Y.-Q. Gao, X.-X. Zhang, J.-L. Li, R.-F. Jin, Y.-D. Chen, F. Zhang, W. Yuan, W.-P. Ma, S. Liang, J. Zhang and Z.-H. Lu, “FlashMP: Fast Discrete Transform-Based Solver for Preconditioning Maxwell’s Equations on GPU”, in The 43rd IEEE International Conference on Computer Design (ICCD’2025), Dallas, USA. 2025, https://doi.org/10.1109/ICCD65941.2025.00118. [25] J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, et al., “SMILEI: A collaborative, open-source, multi-purpose particle-incell code for plasma simulation,” Computer Physics Communications, vol. 222, pp. 351–373, Jan. 2018. https://doi.org/10.1016/j.cpc.2017.09. 024 [26] S. Bastrakov, R. Donchenko, A. Gonoskov, E. Efimenko, A. Malyshev, I. Meyerov, and I. Surmin, “Particle-in-cell plasma simulation on heterogeneous cluster systems,” Journal of Computational Science, vol. 3, pp. 474–479, 2012, https://doi.org/10.1016/j.jocs.2012.08.012. [27] T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, “Contemporary particle-in-cell approach to laserplasma modelling,” Plasma Physics and Controlled Fusion, vol. 57, no. 11, p. 113001, 2015, https://doi.org/10.1088/0741-3335/57/11/113001. [28] R. Bird, N. Tan, S. V. Luedtke, S. L. Harrell, M. Taufer, and B. Albright, “VPIC 2.0: Next Generation Particle-in-Cell Simulations,” IEEE Transactions on Parallel and Distributed Systems, vol. 33, no. 4, pp. 952–963, 2022, https://doi.org/10.1109/TPDS.2021.3084795. [29] L. Fedeli, A. Huebl, F. Boillod-Cerneux, T. Clark, K. Gott, C. Hillairet, S. Jaure, A. Leblanc, R. Lehe, A. Myers, C. Piechurski, M. Sato, N. Zaim, W. Zhang, J.-L. Vay, and H. Vincenti, “Pushing the Frontier in the Design of Laser-Based Electron Accelerators with Groundbreaking Mesh-Refined Particle-In-Cell Simulations on Exascale-Class Supercomputers,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’22), Dallas, TX, USA, 2022, https://doi.org/10.1109/SC41404.2022.00008. [30] S. Markidis, et al., “Multi-scale simulations of plasma with iPIC3D,” Mathematics and Computers in Simulation, vol. 80, no. 7, pp. 1509– 1519, 2010, https://doi.org/10.1016/j.matcom.2009.08.038. [31] J. J. Williams, D. Medeiros, I. B. Peng, and S. Markidis, “Characterizing the Performance of the Implicit Massively Parallel Particle-in-Cell iPIC3D Code,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis (SC’23), Denver, CO, USA, 2023. https://doi.org/10.48550/arXiv.2408.01983 [32] H. Burau, R. Widera, W. Hönig, G. Juckeland, A. Debus, T. Kluge, et al., “PIConGPU: A Fully Relativistic Particle-in-Cell Code for a GPU Cluster,” IEEE Transactions on Plasma Science, vol. 38, no. 10, pp. 2831–2839, 2010. https://doi.org/10.1109/TPS.2010.2064310 [33] M. Bussmann, H. Burau, T. E. Cowan, A. Debus, A. Huebl, et al., “Radiative Signatures of the Relativistic Kelvin-Helmholtz Instability,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’13), Nov. 17–21, 2013, Denver, CO, USA, https://doi.org/10.1145/2503210.2504564. [34] R.A. Fonseca, L.O. Silva, F.S. Tsung, V.K. Decyk, W. Lu, C. Ren, W.B. Mori, S. Deng, S. Lee, T. Katsouleas, and J.C. Adam, “OSIRIS: A Three-Dimensional, Fully Relativistic Particle in Cell Code for Modeling Plasma Based Accelerators,” in Computational Science – ICCS 2002, LNCS 2331, pp. 342–351, Springer, 2002, https://doi.org/10.1007/3-54047789-6 36. [35] B. I. Cohen, A. B. Langdon and A. Friedman, “Implicit time integration for plasma simulation,” Journal of Computational Physics, vol. 46, no. 1, pp. 15–38, 1982. https://doi.org/10.1016/0021-9991(82)90002-X

[36] A. B. Langdon, B. I. Cohen and A. Friedman, “Direct implicit large time-step particle simulation of plasmas,” Journal of Computational Physics, vol. 51, no. 1, pp. 107–138, 1983. https://doi.org/10.1016/ 0021-9991(83)90083-9 [37] G. Chen, L. Chacón and D. C. Barnes, “An energy- and chargeconserving, implicit, electrostatic particle-in-cell algorithm,” Journal of Computational Physics, vol. 230, no. 18, pp. 7018–7036, 2011. https://doi.org/10.1016/j.jcp.2011.05.031 [38] L. Chacón, G. Chen and D. C. Barnes, “A charge- and energyconserving implicit, electrostatic particle-in-cell algorithm on mapped computational meshes,” Journal of Computational Physics, vol. 233, pp. 1–9, 2013. https://doi.org/10.1016/j.jcp.2012.07.042 [39] D. E. Keyes and W. D. Gropp, “A comparison of domain decomposition techniques for elliptic partial differential equations and their parallel implementation,” SIAM Journal on Scientific and Statistical Computing, vol. 8, no. 2, pp. s166–s202, 1987. https://doi.org/10.1137/0908020 [40] D. E. Keyes, “Domain decomposition: A bridge between nature and parallel computers,” Adaptive, Multilevel, and Hierarchical Computational Strategies (A. K. Noor, ed.), vol. 157 of AMD, pp. 293–334, ASME, 1992. [41] B. F. Smith, P. E. Bjørstad, and W. D. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations. Cambridge University Press, 1996. [42] V. Dolean, P. Jolivet, F. Nataf, An introduction to domain decomposition methods - algorithms, theory, and parallel implementation. Society for Industrial and Applied Mathematics, Philadelphia, 2015. https://doi.org/ 10.1137/1.9781611974065 [43] N. Spillane, V. Dolean, P. Hauret, F. Nataf, C. Pechstein, R. Scheichl, “A robust two-level domain decomposition preconditioner for systems of PDEs,” Comptes Rendus. Mathématique, vol. 349, no. 23-24, pp. 1255–1259, 2011. http://dx.doi.org/10.1016/j.crma.2011.10.021. [44] J.-M. Gratien, “A robust and scalable multi-level domain decomposition preconditioner for multi-core architecture with large number of cores,” Journal of Computational and Applied Mathematics, vol. 373, 2020, https://doi.org/10.1016/j.cam.2019.112614. [45] I. Yamazaki, A. Heinlein and S. Rajamanickam, “An Experimental Study of Two-level Schwarz Domain-Decomposition Preconditioners on GPUs,” in IEEE International Parallel and Distributed Processing Symposium (IPDPS’2023), St. Petersburg, FL, USA, pp. 680–689, 2023, https://doi.org/10.1109/IPDPS54959.2023.00073.

[46] H. Anzt, E. Chow, T. Huckle, and J. Dongarra, “Batched Generation of Incomplete Sparse Approximate Inverses on GPUs,” in Proceedings of the 7th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (ScalA’16), pp. 49–56, 2016 , http://dx.doi.org/10.1109/ScalA.2016.011. [47] H. Anzt, J. Dongarra, G. Flegar, and E. S. Quintana-Ortı́, “Batched Gauss-Jordan Elimination for Block-Jacobi Preconditioner Generation on GPUs,” in Proceedings of the 8th International Workshop on Programming Models and Applications for Multicores and Manycores (PMAM’17), pp. 1–10, New York, NY, USA: ACM, 2017, https://doi. org/10.1145/3026937.3026940. [48] H. Anzt, J. Dongarra, G. Flegar, and E. S. Quintana-Ortı́, “VariableSize Batched LU for Small Matrices and Its Integration into BlockJacobi Preconditioning,” in 46th International Conference on Parallel Processing (ICPP’2017), pp. 91–100, 2017, https://doi.org/10.1109/ ICPP.2017.18 [49] H. Anzt, J. Dongarra, G. Flegar, E. S. Quintana-Ortı́, and A. E. Toms, “Variable-Size Batched Gauss-Huard for Block-Jacobi Preconditioning,” Procedia Computer Science, vol. 108, pp. 1783–1792, 2017, (International Conference on Computational Science, ICCS 2017), https: //doi.org/10.1016/j.procs.2017.05.186. [50] H. Anzt, J. Dongarra, G. Flegar, and E. S. Quintana-Ortı́, “Variable-size batched Gauss-Jordan elimination for block-Jacobi preconditioning on graphics processors,” Parallel Computing, vol. 81, pp. 131–146, 2019. https://doi.org/10.1016/j.parco.2017.12.006 [51] F. Y. M. Wan, “An In-Core finite difference method for separable boundary value problems on a rectangle,” Studies in Applied Mathematics, vol. 52, pp. 103–113, June, 1973, https://doi.org/10.1002/sapm1973522103 [52] J. Zhang, C.-B. Zhou, Y.-G. Wang, L.-L. Ju, Q. Du, X.-B. Chi, D.-S. Xu, D.-X. Chen, Y. Liu, and Z. Liu, “Extreme-Scale Phase Field Simulations of Coarsening Dynamics on the Sunway TaihuLight Supercomputer.” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’16), Salt Lake City, UT, USA, 2016, https://doi.org/10.1109/SC.2016.3. [53] S. Cools and W. Vanroose, “The communication-hiding pipelined BiCGstab method for the parallel solution of large unsymmetric linear systems,” Parallel Computing, vol. 65, pp. 1–20, 2017, https://doi.org/ 10.1016/j.parco.2017.04.005