Fast solvers for Tokamak fluid models with PETSC

TL;DR

This paper introduces a multigrid solver for MHD tokamak models in the M3D-C1 code, improving convergence and robustness over existing methods for complex fusion energy simulations.

Contribution

It develops the first multigrid solver tailored for M3D-C1 tokamak models, enhancing solver performance and robustness in fusion energy applications.

Findings

Multigrid solver demonstrates competitive performance on runaway electron models.

The solver shows superior robustness on stellarator models where block Jacobi fails.

Improves convergence for implicit momentum equation solves in M3D-C1.

Abstract

Multigrid (MG) is widely recognized as a highly effective solver for the model problem, the Laplacian, but textbook MG fails on most problems of interest. MG methods have been applied to complex, real-world applications with careful consideration of the physical model and discretization. This work develops the first step in applying MG methods to science and engineering relevant magnetohydrodynamics (MHD) tokamak models in the \textit{M3D-C1} https://m3dc1.pppl.gov fusion energy science code. The semi-implicit time integrator in \textit{M3D-C1} is composed of many linear solves. The implicit advance of the momentum equation is the most challenging and is the focus of this work. The current production solver in \textit{M3D-C1} is a block Jacobi (BJ) preconditioner within a Krylov solver, where blocks group degrees of freedom on planes of constant toroidal coordinate. BJ convergence degrades as the number of planes increases due to the spectral properties of the matrix preconditioned with BJ. The partially magnetic field-aligned, regular toroidal grid structure in \textit{M3D-C1} is amenable to semi-coarsening geometric MG in the toroidal direction. This paper develops such a solver and demonstrates competitive performance on a runaway electron model of a SPARC https://cfs.energy/technology/sparc disruption, and superior robustness on a stellarator model on which the BJ solver fails to converge.