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

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.