Plasma-PEPSC D1.11 Updated Assessment of Developments and Revision of Scientific Challenge, Optimising Magnetically Confined Plasmas with GENE

HORIZON JU Research and Innovation Actions HORIZON-EUROHPC-JU-2021-COE-01-01 European High-Performance Computing Joint Undertaking Plasma Exascale-Performance Simulations CoE 101093261 D1.11 Updated assessment of developments and revision of scientific challenge - Optimising Magnetically Confined Plasmas with GENE WP1: Plasma Simulations - Codes and Grand Challenges Date of preparation (latest version): 30/06/2025 Copyright©2023 – 2027 The Plasma-PEPSC Consortium D1.11: Optimising Magnetically Confined Plasmas with GENE 2 DOCUMENT INFORMATION Deliverable Number D1.11 Deliverable Name Updated assessment of developments and revision of scientific challenge - Optimising Magnetically Confined Plasmas with GENE Due Date 30/06/2025 Deliverable lead MPG Authors Tilman Dannert (MPG), Jordy Trilaksono (MPG), Srikanth Satyanarayana (MPG) and Mou Lin (MPG) Responsible Author Tilman Dannert (MPG) E-mail: [email protected] Keywords GENE, GENE-X “profile prediction” and ITER WP/Task WP1/Task D1.11 Nature R Dissemination Level PU Final Version Date 30/06/2025 Reviewed by David Tskhakaya (IPP CAS), Etienne Renault (SiPearl) and C´eline Scetbun (SiPearl) D1.11: Optimising Magnetically Confined Plasmas with GENE 3 DOCUMENT HISTORY Partner Date Comment Version MPG 05/05/2025 Skeleton version 0.1 MPG 02/06/2025 First draft 0.2 KTH 04/06/2025 Final draft updated for internal review 0.3 MPG 24/06/2023 Revised draft after internal review 0.4 KTH 24/06/2023 Final cleanup for submission 1.0 D1.11: Optimising Magnetically Confined Plasmas with GENE 4 Executive Summary In this deliverable, we describe the progress made for the gyrokinetic fusion turbulence code GENE-X since the last deliverable. The advancements are mainly threefold. First, we progressed on the GPU port of the code so that is can now be run fully on the GPU in a performant way. Second, the employed field solver library PARALLAX[1] is now also fully ported to GPUs via the backend PAccX[2]. Additional work has been invested in making the fieldsolver available on AMD GPUs and to employ a library solution instead of the self-written FGMRES-Multigrid solver. Third, the implementation of the blockstructured grids reached a level of usability on the CPU which makes it possible to reduce the number of gridpoints by a factor of 4-8. It is not yet possible to combine all of the achievements, but with using the GPU code for GENE-X and the fieldsolver, we arrive at a FOM of around 20000 on 32 MI300A GPUs. The roadmap has been refined to further close the gap to the envisaged FOM of 6000. D1.11: Optimising Magnetically Confined Plasmas with GENE 5 Contents 1 Introduction 6 2 Current Status of Scientific Challenge 6 3 Code Developments 7 3.1 GPU offloading and optimization efforts in GENE-X . . . . . . . . . . . 7 3.1.1 General GPU porting progress . . . . . . . . . . . . . . . . . . . . 7 3.1.2 AMDGPUsupport.......................... 10 3.1.3 Optimization: GPU memory usage . . . . . . . . . . . . . . . . . 11 3.1.4 Optimization: GENE-X support for PARALLAX’s mesh reordering 12 3.2 Block-structured grids for GENE-X . . . . . . . . . . . . . . . . . . . . . 13 3.3 Field solver optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3.1 AMDGPUPorting.......................... 17 3.3.2 Matrix reordering optimization . . . . . . . . . . . . . . . . . . . 18 3.3.3 Multi-GPU parallelization . . . . . . . . . . . . . . . . . . . . . . 19 4 Current Status of Gap Analysis 19 4.1 Targettestcase................................ 20 4.2 Reducedtestcase............................... 20 4.3 Extrapolation................................. 21 4.4 Estimation of Compute Resources for the Target run . . . . . . . . . . . 22 5 Revision of Roadmap 23 5.1 Partial MPI-Parallelization of the Fieldsolver . . . . . . . . . . . . . . . . 23 5.2 Overall Optimization of the GPU code of GENE-X . . . . . . . . . . . . 24 5.3 Porting BSG of GENE-X to GPU . . . . . . . . . . . . . . . . . . . . . . 24 6 Conclusion 24 D1.11: Optimising Magnetically Confined Plasmas with GENE 6 1 Introduction GENE and GENE-X are grid-based microturbulence simulation codes solving the gyrokinetic Vlasov-Maxwell system of integro-differential equations. Both codes work on a 5D physical grid consisting of three spatial coordinates and two velocity space coordinates, the velocity along the field lines v∥and the magnetic moment µ. Whereas GENE uses a field aligned coordinate system with the spatial coordinates x(radial), y(binormal) and z(along the field line), GENE-X employs a fieldline-coordinate independent (FCI) approach. The latter uses the cylindrical coordinate system R,Zand ϕin real space and the same velocity coordinates as GENE. Both codes have two main parts, first the computation of the right-hand side of the Vlasov equations which mainly consists of stencil operations in different directions and the computation of the nonlinearity (which involves 1D FFTs in the case of GENE). The second important part is the solution of the field equations which are linear equations solved in all poloidal planes independently. The input to these equations are the zeroth (density) and first (parallel current) moment of the distribution function, which are computed by a reduction (MPI Allreduce) operation over the velocity space. In previous deliverables (D1.2 and D1.6), we have shown the shortcomings of the two codes and first solutions with respect to the target challenge. One of the biggest tasks is to port the GENE-X code to work on GPUs and this task has been progressed a lot in the last year. A second task has been the optimization of the fieldsolver scaling in the velocity space directions. As a predecessor step to this scaling optimization, we focused on running the fieldsolver on GPUs on all available GPU architectures. A highly optimized fieldsolver takes away some pressure from the scaling needs. But still, improvements of the velocity space scaling will be a main focus in the remaining 18 months of the project. The third point shown in D1.2 was the problem of strongly varying background temperature which leads either to too many points or to an underresolution in velocity space. To overcome this problem, we originally wanted to develop a library to be used by both GENE codes, but finally refrained from this solution and focused on a specific implementation for the GENE-X code. The GENE code already supports block-structured grids (at least on the CPU). As described in the former two deliverables, GENE-X is the most time-consuming part of the simulation, so we focus in this deliverable on the progress of this code towards the scientific challenge. 2 Current Status of Scientific Challenge The Scientific Challenge is not changed. We still head for an ab-initio prediction of density and temperature profiles in the ITER tokamak using the three codes GENE for the core, TANGO for the transport evolution and GENE-X for the edge to the wall part. We refer for a detailed description of the physical and numerical setup of the scientific target case to Deliverable D1.2. D1.11: Optimising Magnetically Confined Plasmas with GENE 7 3 Code Developments 3.1 GPU offloading and optimization efforts in GENE-X 3.1.1 General GPU porting progress The main code development results have been published in a recent paper titled ”OpenACC and OpenMP-Accelerated Fortran/C++ Gyrokinetic Fusion Code GENE-X for Heterogeneous Architectures” [3]. It has been accepted to the Platform for Advanced Scientific Computing Conference (PASC) 2025. This paper covers mainly two points. First, the interoperation efforts between the main Fortran layer and the auxiliary C++ layer of GENE-X and second, the GPU offloading efforts of GENE-X via OpenACC and OpenMP offload. Both directive-based approaches use very similar concepts, just the syntax is different. The largest difference is the portability aspect, as OpenMP is much more widely portable across different GPU vendors, e.g., NVIDIA, AMD, Intel, than OpenACC. Both directives are implemented in GENE-X, but we will focus in the project on the more portable OpenMP offload set of directives. This paper excludes topics related to field solvers that are powered by PARALLAX and its PAccX libraries. It also excludes comparisons with kernels running on Fortran OpenMP on the CPU. The excluded field solver regions are relevant parts of the code, but even without it the results are important as the native regions of the application side take most of the computing time on the CPU. All benchmarks and performance comparisons in the paper are done strictly within the auxiliary C++ layer. All benchmark cases used in the paper are based on the production run of TCV-X21 tokamak published in [4] and are run on MareNostrum 5 ACC partition with 4 NVIDIA H100 per node. OMP ACC OMPX 0 2 4 6 8 10 Runtime per timestep per MPI process [s] Op Vlasov static Op Vlasov dynamic Op coll BGK Misc. kernels isend & irecv allreduce allreduce sync Figure 1: Runtime comparison published in [3]. One benchmark case used in the paper has a size of 5.4 billion grid points and is spread across 80 MPI processes and 80 GPUs to get an overview of the performance gain from GPU offloading. This makes up a quarter of the total grid size in [4] but the size per MPI process is eight times larger. In summary, we achieve a speedup of 6.9×with OpenMP Target offload directives in the elapsed time of the timestep as shown in Figure 1. The three main operators or compute kernels, i.e., Vlasov static,Vlasov dynamic and BGK collision, all together gain a significant speedup of 21×. In the paper, a different case with 212,387 grid points on a single GPU with an MPI D1.11: Optimising Magnetically Confined Plasmas with GENE 8 10 210 1100101102 Arithmetic Intensity [FLOP/Byte] 102 103 104 105 Performance [GFLOP/s] HBM: 1631 GB/s FP64 (non-tensor) Peak: 28.1 TFLOP/s OpenACC OpenMP Offload Op Vlasov static Op Vlasov Dynamic Op coll BGK Op mom Maxwell's eqs Op mom Ohm's law Op mom coll Op linear combination Op axpy Figure 2: Roofline model of OpenACC and OpenMP offload kernels in GENE-X from [3] provided by NVIDIA Nsight Compute on MareNostrum5. process provides us with the roofline model as shown in Figure 2. This figure shows that all GENE-X compute kernels are memory-bound by the 1631 GB/s bandwidth of GPU high-bandwidth memory (HBM). The roofline model also shows that the default configuration of OpenMP offload performs slightly better than OpenACC. After we see some of the metrics related to the L1 and L2 caches provided by Nsight Compute, we conclude that OpenMP offload performs better for memory-bound kernels with coalesced access and high cache reuse, as the default configuration results in better memory access pattern and warp occupancy. 1 2 4 8 16 32 64 Number of nodes 0 10 20 30 Runtime per timestep per MPI process [s] OpenMP (CPU) 1 2 4 8 16 32 64 Number of nodes 0.0 0.5 1.0 1.5 Runtime per timestep per MPI process [s] OpenACC 1 2 4 8 16 32 64 Number of nodes 0.0 0.5 1.0 1.5 Runtime per timestep per MPI process [s] OpenMP Offload allreduce sync allreduce isend & irecv Misc. kernels Op coll BGK Op Vlasov dynamic Op Vlasov static Total speedup 1 2 4 8 16 32 64 Number of nodes 0 10 20 30 40 Speedup w.r.t. CPU OpenACC 1 2 4 8 16 32 64 Number of nodes 0 10 20 30 40 Speedup w.r.t. CPU OpenMP Offload 4 8 16 32 64 128 256 Number of MPI processes 4 8 16 32 64 128 256 Number of GPUs and MPI processes 4 8 16 32 64 128 256 Number of GPUs and MPI processes 4 8 16 32 64 128 256 Number of GPUs and MPI processes 4 8 16 32 64 128 256 Number of GPUs and MPI processes Figure 3: Strong scaling of the runtime and GPU speedup factor published in [3]. Figure 3 shows the strong scaling of OpenMP on CPUs, OpenACC and OpenMP offload on GPUs. Information on the grid size and the MPI decomposition is provided in Table 1. The same case as the overview case in Figure 1 cannot be incorporated in this strong scaling analysis due to the 64 GB GPU memory limit when we try to run the same size within one node. Therefore, these cases were chosen so that the one-node case can fit into GPU memory. This also motivates the GPU memory optimization D1.11: Optimising Magnetically Confined Plasmas with GENE 9 efforts reported in the next two sections. The compute kernels on GPUs show similar performance improvements compared to the ones on CPUs as the number of GPUs used is increased. However, ghost exchange and synchronization time become more prominent as bottlenecks when more GPUs are used. The speedup factors from GPU offloading are diminished as a result of the increasing contribution of the bottlenecks. However, the speedup factor of the Vlasov static operator is slightly increasing. Table 1: Problem size and MPI decomposition for the strong scaling analysis in Figure 3 [3]. Dim. Num. points MPI processes for each case RZ 212387 1 1 1 1 1 1 1 φ8 4 4 4 4 4 4 8 v∥24 1 1 1 1 1 1 1 µ12 1 1 2 4 8 16 16 species α2 1 2 2 2 2 2 2 Number of MPI processes 4 8 16 32 64 128 256 Number of GPUs 4 8 16 32 64 128 256 Number of nodes 1 2 4 8 16 32 64 In Figure 1 and Figure 3, we observe a significant runtime contribution from the communication regions, especially the synchronization time of MPI Allreduce. Focusing on MPI Allreduce, the overwhelming synchronization time compared to the pure communication time may also indicate a similar phenomenon for the ghost exchange region. However, we cannot isolate the pure communication time easily with non-blocking MPI communications. This indicates some load imbalances occurring during the timestep, but the damage is observed mainly in the region of blocking communication, e.g. MPI Allreduce. Due to this, the total performance is undermined by load imbalances that become stronger as the compute kernels become faster by the GPU offloading. We find that the synchronization time of a central collective communication routine (MPI Allreduce) takes more than 50% of the time per timestep in the OpenMP offload cases for the cases with 16 or 32 GPUs in Figure 3. This bottleneck is then followed by the ghost exchange which takes approximately 16%. Outside the scope of the paper, we have performed a preliminary investigation to find the source of the load imbalances in GENE-X. We have extended our built-in profiler in GENE-X to support tracing of the runtime. GENE-X built-in profiler only stores cumulative performance data such as total runtime, number of calls, and time per call for each profiling region. We recently extended the capability of our built-in profiler, to store all the timestamps for each start and stop instances of all profiling regions. This allows us to construct a trace view of the profiling regions in GENE-X. The trace file in csv format not only can be visualized our homemade visualization script, which can show the profiler regions of each MPI process, but can also be converted to json format according to the Perfetto scheme and visualized by the Perfetto-UI web interface. Figure 4 is the trace visualization generated by our script and it demonstrates the load imbalance caused by one of our diagnostic I/O, i.e., diagnose em fields. The mitigation of the D1.11: Optimising Magnetically Confined Plasmas with GENE 16 1 2 3 4 Resolution 10 −4 10 −3 10 −2 10 −1 || f − f MMS || p || f MMS || p Convergence with L 2 norm Second order reference ions electrons es_pot A_par E_par 1 2 3 4 Resolution 10 −3 10 −2 10 −1 10 0 || f − f MMS || p || f MMS || p Con ergence with L ∞ norm Second order reference ions electrons es_pot A_par E_par Figure 8: Convergence study with MMS: Observed order of accuracy for distribution function and field derivatives, compared to second-order reference lines across four progressively finer resolutions. Resolution1 Resolution2 Resolution3 Resolution4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Slowdown Slowdown Comparison Across Resolutions Number of Blocks = 2 Number of Blocks = 3 Figure 9: Code slowdown: Comparison of GENE-X code slowdown for different number of blocks. D1.11: Optimising Magnetically Confined Plasmas with GENE 17 A secondary consideration when reducing the velocity space grid points is the corresponding performance scaling. Figure 9 depicts the slowdown of the GENE-X code which is basically the ratio of the runtime with and without BSG. This quantity is compared for different number of blocks. Usually the number of blocks for the target case will be limited to 3-5. A first observation from the plot is that the slowdown is generally greater than 1 which is due to the additional interpolation costs performed at the block boundaries. Further, we also observe that with increasing number of points the slowdown saturates between 1.2 - 1.4. However, we expect this number to go down since the targeted ITER case is supposed to be much finer than the test resolutions used here while keeping the number of blocks fairly constant. Therefore the ratio between the number of points in the poloidal plane that require interpolation and the total number of RZ points reduces. 3.3 Field solver optimization The field solver in GENE-X, PARALLAX, handles the quasi-neutrality equations, the generalized Amp`ere’s law, and Ohm’s law on two-dimensional drift planes. These equations can be reformulated as elliptic Helmholtz-type equations of the form: λϕ −ξ∇ · (c∇⊥ϕ) = ρ, (2) where the perpendicular gradient is defined as ∇⊥=∇ − b(b· ∇), with bdenoting the unit vector along the magnetic field. The coefficients λ,ξ, and care known functions, and ρrepresents the charge density. These equations are discretized using finite difference methods and are solved using the flexible generalized minimal residual method (fGMRES), preconditioned with a geometric multigrid scheme. 3.3.1 AMD GPU Porting In the previous development update, PARALLAX was successfully ported to NVIDIA GPUs by implementing a new CUDA C++ backend of the solver, PAccX, achieving a significant socket-to-socket speedup (approximately 20×) compared to the original Fortranbased CPU solver. In this update, substantial efforts have been dedicated to extending the portability of PARALLAX/PAccX to AMD GPUs. To accomplish this, we adopt the hipifly approach [7], which enables porting the CUDA backend to AMD GPUs by redefining CUDA function calls as their HIP counterparts through a header file, e.g., #define cudaMalloc hipMalloc #define cudaFree hipFree #define cudaMemcpy hipMemcpy #define cudaMemcpyHostToDevice hipMemcpyHostToDevice ... We shall distinguish hipifly from the HIPIFY [8] approach, which modifies the CUDA code to the HIP code with hipify-clang and hipify-perl provided by the ROCm library. By using the header-based hipifly approach, we achieve multi-platform portability of the solver while maintaining a single code base, with CUDA as the primary backend. Any future development on the CUDA backend will automatically propagate to the HIP backend, provided that corresponding functions are available in ROCm/HIP. In cases D1.11: Optimising Magnetically Confined Plasmas with GENE 18 where no equivalent functionality exists, we will rely on third-party libraries that are compatible with both backends. 3.3.2 Matrix reordering optimization While the updated field solver runs on both NVIDIA and AMD GPUs, we have observed a significant performance gap on AMD GPUs, despite comparable hardware specifications. For example, in a large circular test case with 37,679,076 grid points, the solving time for the CUDA backend on an NVIDIA A100 GPU is 2.09 s, whereas the HIP backend on an AMD MI250X takes 3.02 s. By profiling the field solver, we identified that the performance difference between the CUDA and HIP backends is primarily caused by the Gauss–Seidel red-black kernel in the multigrid preconditioner, which accounts for approximately 80% of the total computation time. The original mesh ordering results in a staggered memory access pattern for the red-black kernel, which is suboptimal for GPU acceleration. The AMD MI250X is more severely affected by this discontinuous memory access than the NVIDIA A100, due to its smaller L2 cache size. To address this issue, we reorder the matrix according to the red-black partition, ensuring that grid points of the same color are stored contiguously in GPU memory (Figure 10). Figure 10: Sparsity parttern of the original mesh ordering and red-black ordering. Table 4 presents the performance comparison of the A100 and MI250X GPUs, with and without red-black reordering. As shown, adopting the GPU memory-friendly redblack ordering yields a greater speedup on the HIP-MI250X backend (×1.65) than on the CUDA-A100 backend (×1.30), thereby significantly narrowing the performance gap between HIP and CUDA from 44.5% to 13.6%. The HIP-MI250X backend benefits more from the matrix reordering than the CUDA-A100 backend due to the fact that MI250X has smaller L2 Cache (16 MB) than A100 (40 MB). As a result, the MI250X relies more heavily on coalesced and continuous memory accesses to maximize cache efficiency and minimize memory latency. D1.11: Optimising Magnetically Confined Plasmas with GENE 19 Table 4: Performance comparison with and without red-black reordering Original (s) Red-black reordering (s) Speedup NVIDIA A100 2.09 1.61 1.30 AMD MI250X 3.02 1.83 1.65 3.3.3 Multi-GPU parallelization To further improve the scalability of the field solver, we explore the possibility of solving the field equations in parallel using multiple GPUs. Preliminary efforts have been conducted on AMD hardware using rocALUTION [9], a sparse linear algebra library designed for fine-grained parallelism on the ROCm platform. We developed a standalone unit test program that reads the raw matrix and right-hand side (RHS) data from an ITER-scale test case (approximately 11 million grid points) and solves it using rocALUTION. The matrix is initially scaled by its diagonal entries and then solved using the BiCGStab linear solver, preconditioned with the classical Ruge–St¨uben algebraic multigrid method. The strong scaling results of this test case on 1 to 8 MI250X GPUs are shown in Figure 11. Figure 11: Multi-GPU strong scaling with rocALUTION for the ITER-size test problem. As shown in Figure 11, the general out-of-the-box parallelism provided by rocALUTION still falls short of ideal scaling. To address this, we plan to enhance the multi-GPU performance by developing a customized communication pattern tailored to the specific matrix structure of our problem. 4 Current Status of Gap Analysis This section deals with the progress of the code GENE-X towards the target challenge. All of the described code improvements of the last section help in reducing the runtime of the code and makes it possible to further close the gap to the target FOM, although not all of them can be used together at the moment. The Figure of Merit (FOM) has D1.11: Optimising Magnetically Confined Plasmas with GENE 20 been defined in the proposal and in previous deliverables to be FOM = Time per time step ×total stream BW FP size ×total number of grid points full where the ”total number of grid points full” stand for the number of grid points that would be necessary without algorithmic changes like block-structured grids, which substantially reduce the number of needed grid points. 4.1 Target test case The envisaged target test case definition did not change. We will simulate the turbulence in the edge and the Scrape Off Layer (SOL) of the ITER tokamak. GENE-X uses a Cartesian grid in the poloidal RZ-planes and a discrete number of poloidal planes along the toroidal direction. For the planned ITER runs, we estimate that at least 48 millions points are needed in the poloidal RZ plane and 32 of these poloidal planes are needed to discretize along the toroidal direction. For the velocity space we plan to use ≈80 points in v∥(without BSG first, for BSG advantages see later) direction and ≈20 in µ direction. With two species as in the GENE case, this multiplies to a total of roughly 5 trillion points (4.9·1012). We have to run the simulation long enough to have good statistics in the turbulent phase. For this we expect 50k to 100k time steps. As all timesteps are nearly equivalent in runtime, we only run 10-100 timesteps for the performance assessment runs in this section. 4.2 Reduced test case For assessing the performance of the code, we have to run a reduced test case, that still contains all the relevant physics (at least the same code paths as with the target case), but is much less resource-intensive. As in all the previous deliverables D1.2 and D1.6, we are using the following reduced case for GENE-X, which is roughly a factor 800-1000 smaller than the target case. We use a grid spacing in RZ of 9 ·10−4Lref (with Lref = 6.41327m) which leads to 339411 RZ points. Instead of 32 poloidal planes we use only 16 and also in v∥and µdirection we reduce by a factor of 2 each to 40 v∥points and 10 µpoints. We still use 2 species to cover the same physics. We tried to weak scale the test case up to the target values for the dimensions φ,v∥ and µbut keep the RZ points constant. Starting with 8 nodes of the MPCDF’s VIPERGPU machine (each node is equipped with two AMD MI300A APUs) we scale up to 64 nodes and get the results shown in Table 5. The scaling from 8 to 16 nodes in the φdirection is fine, but then a big jump happens for larger cases. This is a clear indication that something is not working correctly in this architecture. Doing the same runs on the Raven system at MPCDF, which is equipped with four A100 for each node and running the fieldsolver on the CPU, we get a decent weak scaling as shown in figure 13. D1.11: Optimising Magnetically Confined Plasmas with GENE 21 Table 5: Weak scaling. RZ φ v∥µspecies MPI nodes Time per timestep 330661/1 16/2 40/2 10/2 2/2 16 8 12.66 330661/1 32/440/2 10/2 2/2 32 16 11.19 330661/1 32/4 80/410/2 2/2 64 32 73.47 330661/1 32/4 80/4 20/42/2 128 64 63.60 0 20 40 60 80 100 120 140 160 8 16 32 64 exchange rhs static collisions Ohm Maxwell time per timestep Number of VIPER-GPU nodes (with 2 MI300A) Weak scaling of GENEX on VIPER-GPU Figure 12: Breakdown of a weak scaling on VIPER-GPU for the reduced test case. The fieldsolver is executed on the GPU and the data is passed via device pointers from GENE-X to PARALLAX/- PAccX. In the last column of Table 6, you can find the computed Figure of Merit of the runs. It is around 20000 which is still a factor of 3 too high compared to the target FOM of 6000. As we can see from the graph, a big part of the runtime goes into the fieldsolver, which has been ported to GPU already, but the combination of the two codes does not yet work fine. 4.3 Extrapolation Extrapolation of the results from the reduced test case to the target case on an Exascale system is a difficult task. Especially because we do not know how the scaling develops on a larger scale. We start from the reduced test case on 16 VIPER-GPU nodes which shows a time per timestep of 11.2s. One VIPER-GPU node has a performance of roughly 100TFlops/s. For reaching Exascale, one would need a total of 10000 of these kind of nodes with 2 Table 6: A breakdown of different code parts with the scaling is shown in figure 12. nodes timestep exchange static coll Ohm Maxwell Allreduce Barrier FOM 8 11.180 3.192 0.872 1.127 0.739 5.624 0.929 0.518 19895 16 11.989 3.375 0.872 1.509 0.832 5.798 1.401 0.657 21334 32 12.338 3.382 0.834 1.631 0.883 6.012 1.394 0.824 21955 64 12.152 3.386 0.834 1.910 1.017 5.416 1.809 0.898 21625 D1.11: Optimising Magnetically Confined Plasmas with GENE 22 0 2 4 6 8 10 12 14 8 16 32 64 exchange rhs static collisions Ohm Maxwell time per timestep Number of RAVEN GPU nodes (with 4 A100) Weak scaling of GENEX on RAVEN Fieldsolver on CPU Figure 13: Weak scaling on Raven with 4 A100 per node. The fieldsolver is executed on the CPU. MI300A each, hence a total of 20000 MI300A. On the problem side we have to scale up a factor of 400 to come to the target case. With the (naive) assumption of ideal scaling, we would come to a time per timestep of the target case on an Exascale system of 11.2s·400 10000/16 = 7.2s. The FOM does not change by ideal scaling. 4.4 Estimation of Compute Resources for the Target run In the last year of the project, we are planning for the simulation runs for the target case. To get some more realistic estimation of the compute resources needed for these runs, we try to get these number from the FOM measured and the one we aim for. Taking the results of the last subsection 4.3, we can estimate the needed compute resources to do the target run on a suitable system. We need 100000 timesteps of 7.2s on 20000 MI300A (or equivalent), this is a total of 200 hours on 20000 MI300A or 4 million GPU hours. This large amount of computing resources will decrease in the same way as the FOM decreases in the follow-up on this project. Hence, for the target case with a FOM<6000, we should end up with around 1.2 million GPU hours. To become more concrete, we estimate the needed resources on LUMI-G, MareNostrum 5 and Jupiter, assuming ideal scaling (which is not given at the moment). So all of the numbers are just rough estimates. D1.11: Optimising Magnetically Confined Plasmas with GENE 23 VIPER-GPU LUMI-G MN5 Jupiter FP64 peak perf. 1 node 100T 128T 157T 167T mem bw per node 7.2T 10T 10T 12T # nodes full machine 16 2978 1120 6000 time per timestep (red.) 11.2 time per timestep (full) 4480 17 46 7.2 node-hours needed (current,ideal) 1.4M 1.4M 1.2M node-hours needed (target,FOM) 420k 420k 360k available in recent call 1.8M 580k 4.7M But we can already with these numbers see that one production run of the target case needs a considerable amount of the available computing time on the largest EuroHPC machines. And it becomes even more clear that we have to reach at least the FOM promised. Ideally we will go below the FOM of 6000 substantially. 5 Revision of Roadmap From the weak scaling runs we can see that we need further development to overcome the still existing gap. The most obvious problem is the bad scaling of the field solver with the velocity space processes, which is an inherent problem of (gyro-)kinetic plasma simulations and comes from the 3D nature of the field equations. Next is that the GPU code has still potential to be optimized, our kernels do not yet reach the roofline of the GPUs used and also the occupancy metric is quite low. Finally, by porting the block-structured grid approach to GPU, we can use it for the target case. This would decrease the number of velocity space points and speed up the time to solution. 5.1 Partial MPI-Parallelization of the Fieldsolver Although the problem might become less prominent due to optimizations on the fieldsolver itself, we still need a moderate parallelization of the update and solve step over the (unsused) velocity space processes. To come to a large number of GPUs used, as necessary for Exascale runs, we have to highly parallelize GENE-X in the velocity space directions. Unfortunately, this parallelization is usually across nodes, because the innermost φ-direction is parallelized first. A full parallelization is not advised as then the local problem might become too small for one GPU, so we aim at a moderate parallization level of 4-8 MPI processes. However, the final setup has to be determined by performance measurements of the optimized code. So, the next steps would be 1. Further assess the library solution with rocALUTION on AMD GPUs and with the ginkgo[10] library as a portable solution. 2. Use the MPI parallelization in these libraries to parallelize over the velocity space processes from GENE-X 3. Assess performance of these runs and find the “sweet spot” of parallelization level D1.11: Optimising Magnetically Confined Plasmas with GENE 24 5.2 Overall Optimization of the GPU code of GENE-X We have a well performing GPU code now, but it is not yet very reliable on different architectures. Also the usage of the GPU ported PARALLAX/PAccX library is not yet always working. In addition the op rhs static kernel is quite large and does still have optimization potential for faster execution on the GPUs. So the next steps here are 1. Stabilize the GPU code on the two main architectures for NVIDIA and AMD GPUs. Target systems here are MPG’s computers RAVEN and VIPER-GPU but also LUMI-G and MareNostrum 5. 2. Measure kernel performance of the most time-consuming kernels on these architectures and optimize the kernels 3. Integrate the MPI-parallelized PARALLAX solver 4. Integrate the GPU ported block-structured grids 5.3 Porting BSG of GENE-X to GPU With regard to the upcoming/remaining work, we plan to target two major milestones. First, we would like to perform code validation of GENE-X using BSG. For this, we intend to perform simulations with TCV-X21 configuration which has been well validated using GENE-X [4]. Finally, we plan to port the BSG implementation to GPU. This makes it then possible to combine all of the achievements in the project in the target run. We expect a reduction in points for the ITER case of 4-6 (we will not get the full factor of 10, because in general a higher velocity space resolution is advisable in these runs, but were not possible without BSG). This would come with a reduction in the runtime by a factor 3-4×and simultaneously reduce the FOM by the same factor. 6 Conclusion We presented the recent code developments and how they aim at reaching the target simulations. Making use of available GPUs is the main part which has been followed in the GPU porting of GENE-X and PARALLAX and successfully reached for both of them independently. On the numerical side, we finalized the implementation of a CPU version of the block-structured grids for GENE-X which promised a reduction in runtime of a factor 3-4. Validation in a larger run is now necessary and then a porting of these changes also to GPU and integration with the rest of the GENE-X GPU code. In the project we reached now a state where all developments on their own did quite well progress and in the remaining time we have to focus on optimization and integration of all of these developments into a common simulation. Only then, we will be able to reach the promised Figure of Merit. D1.11: Optimising Magnetically Confined Plasmas with GENE 25 References [1] phoenix-public: PARALLAX, 2025. https://gitlab.mpcdf.mpg.de/ phoenix-public/parallax. [2] phoenix-public: PAccX, 2025. https://gitlab.mpcdf.mpg.de/phoenix-public/ paccx. [3] J. Trilaksono, P. Ulbl, J. J. Williams, C.-M. Pfeiler, M. Finkbeiner, T. Dannert, E. Laure, S. Markidis, and F. Jenko. Openacc and openmp-accelerated fortran/c++ gyrokinetic fusion code gene-x for heterogeneous architectures. In Proceedings of the Platform for Advanced Scientific Computing Conference, PASC ’25, New York, NY, USA, 2025. Association for Computing Machinery. [4] Philipp Ulbl, Thomas Body, Wladimir Zholobenko, Andreas Stegmeir, Jan Pfennig, and Frank Jenko. Influence of collisions on the validation of global gyrokinetic simulations in the edge and scrape-off layer of TCV. Physics of Plasmas, 30(5), 2023. [5] D. Jarema, H. J. Bungartz, T. Goerler, F. Jenko, T. Neckel, and D. Told. Blockstructured grids for Eulerian gyrokinetic simulations. Computer Physics Communications, 198:105–117, JAN 2016. [6] D. Jarema, H. J. Bungartz, T. Goerler, F. Jenko, T. Neckel, and D. Told. Blockstructured grids in full velocity space for Eulerian gyrokinetic simulations. Computer Physics Communications, 215:49–62, JUN 2017. [7] hipifly, 2025. https://github.com/amd/HPCTrainingExamples/tree/main/ hipifly. [8] hipify, 2025. https://github.com/ROCm/HIPIFY. [9] rocALUTION, 2025. https://github.com/ROCm/rocALUTION. [10] Hartwig Anzt, Terry Cojean, Goran Flegar, Fritz G¨obel, Thomas Gr¨utzmacher, Pratik Nayak, Tobias Ribizel, Yuhsiang Mike Tsai, and Enrique S. Quintana-Ort´ı. Ginkgo: A Modern Linear Operator Algebra Framework for High Performance Computing. ACM Transactions on Mathematical Software, 48(1):2:1–2:33, February 2022.