Memory- and compute-optimized geometric multigrid GMGPolar for curvilinear coordinate representations -- Applications to fusion plasma

TL;DR

This work tackles efficient solution of the gyrokinetic Poisson equation on curvilinear tokamak cross sections using a matrix-free geometric multigrid GMGPolar. It introduces a fully refactored, object-oriented GMGPolar with two matrix-free implementations (Give and Take), higher-order implicit extrapolation, and FMG/W/F-cycle features to achieve linear complexity and elevated accuracy. The numerical results demonstrate substantial solver speedups (up to 16–18× for the solver, 25–37× as a preconditioner) and reduced memory footprints across multiple geometries, confirming strong and weak scalability. The approach enables fast, memory-efficient high-fidelity plasma simulations on HPC systems, with potential GPU porting and multi-patch extensions for more complex tokamak geometries.

Abstract

Tokamak fusion reactors are actively studied as a means of realizing energy production from plasma fusion. However, due to the substantial cost and time required to construct fusion reactors and run physical experiments, numerical experiments are indispensable for understanding plasma physics inside tokamaks, supporting the design and engineering phase, and optimizing future reactor designs. Geometric multigrid methods are optimal solvers for many problems that arise from the discretization of partial differential equations. It has been shown that the multigrid solver GMGPolar solves the 2D gyrokinetic Poisson equation in linear complexity and with only small memory requirements compared to other state-of-the-art solvers. In this paper, we present a completely refactored and object-oriented version of GMGPolar which offers two different matrix-free implementations. Among other things, we leverage the Sherman-Morrison formula to solve cyclic tridiagonal systems from circular line solvers without additional fill-in and we apply reordering to optimize cache access of circular and radial smoothing operations. With the Give approach, memory requirements are further reduced and speedups of four to seven are obtained for usual test cases. For the Take approach, speedups of 16 to 18 can be attained. In an additionally experimental setup of using GMGPolar as a preconditioner for conjugate gradients, this speedup could even be increased to factors between 25 and 37.

Figures (12)

• Figure 1: Particular solution of \ref{['eq:poisson']} on a 3D tokamak with Culham geometry. The figure shows a 3D tokamak geometry spanning over 220 degrees in toroidal direction and clipped over the remaining 140 degrees. It furthermore show several 2D cross sections in the clipped part of the geometry.
• Figure 2: Visualization of the Shafranov and Czarny cross section geometries of a tokamak. As both geometries can be described in curvilinear coordinates, mappings from a rectangular grid $[r_1,1.3]\times[0,2\pi]$ (center) to the considered geometries are indicated. The Shafranov geometry (left) is given by the mapping $F_S:[r_1,1.3]\times[0,2\pi]\rightarrow \mathbb{R}^2$ and the Czarny geometry (right) is given by the mapping $F_C:[r_1,1.3]\times[0,2\pi]\rightarrow \mathbb{R}^2$; see \ref{['eq:shafranov']} and \ref{['eq:czarny']}. For the existence of the invertible mappings $F_S^{-1}$ and $F_C^{-1}$, we need $r_1>0$. The depicted grids are symbolically refined at 2/3 of the generalized radius.
• Figure 3: Refactored class layout highlighting the modular structure and object-oriented design approach.
• Figure 4: Combined circle-radial smoother and indexing for the smoothing operations. Circle and radial lines colored black and white on a grid of dimension $16 \times 32$ with eight circle lines of 32 nodes and 16 radial lines of eight nodes; for visualization simplification, the curvilinear lines are shown without the nodes (left). Optimized grid indexing for a periodic grid of dimension $8 \times 4$ with four circle lines of four nodes and four radial lines of four nodes. Vertical lines of the same color represent circular smoothers, while horizontal lines of the same color correspond to radial smoothers (right).
• Figure 5: Dependency graph for the application of the system matrix $A^{(i)}$, of multigrid level $i\in\{0,\ldots,L-1\}$, using the Give implementation. Each vertex represents a task corresponding to a line of nodes. The edges visualize the dependencies between these tasks.

Theorems & Definitions (1)

• Remark 1