A bespoke multigrid approach for magnetohydrodynamics models of magnetized plasmas in PETSc

TL;DR

This paper develops a custom multigrid solver within PETSc for magnetohydrodynamics models, demonstrating robustness and scalability for complex plasma simulations with localized nonlinearities.

Contribution

It introduces a fully integrated geometric multigrid solver tailored for MHD models in PETSc, enabling fast, robust, and application-specific solutions in high-performance computing environments.

Findings

Solver is robust and scalable for complex MHD problems

Successfully handles highly localized current sheets

Demonstrates potential for production-level applications

Abstract

Fully realizing the potential of multigrid solvers often requires custom algorithms for a given application model, discretizations and even regimes of interest, despite considerable effort from the applied math community to develop fully algebraic multigrid (AMG) methods for almost 40 years. Classic geometric multigrid (GMG) has been effectively applied to challenging, non-elliptic problems in engineering and scientifically relevant codes, but application specific algorithms are generally required that do not lend themselves to deployment in numerical libraries. However, tools in libraries that support discretizations, distributed mesh management and high performance computing (HPC) can be used to develop such solvers. This report develops a magnetohydrodynamics (MHD) code in PETSc (Portable Extensible Toolkit for Scientific computing) with a fully integrated GMG solver that is designed to demonstrate the potential of our approach to providing fast and robust solvers for production applications. These applications must, however, be able to provide, in addition to the Jacobian matrix and residual of a pure AMG solver, a hierarchy of meshes and knowledge of the application's equations and discretization. An example of a 2D, two field reduced resistive MHD model, using existing tools in PETSc that is verified with a ``tilt" instability problem that is well documented in the literature is presented and is an example in the PETSc repository (\path{src/ts/tutorials/ex48.c}). Preliminary CPU-only performance data demonstrates that the solver can be robust and scalable for the model problem that is pushed into a regime with highly localized current sheets, which generates strong, localized non-linearity, that is a challenge for iterative solvers.