Scalable algorithms for 3D extended MHD.*

In the modeling of plasmas with extended MHD (XMHD), the challenge is to resolve long time scales while rendering the whole simulation manageable. In XMHD, this is particularly difficult because fast (dispersive) waves are supported, resulting in a very stiff set of PDEs. In explicit schemes, such stiffness results in stringent numerical stability time-step constraints, rendering them inefficient and algorithmically unscalable. In implicit schemes, it yields very ill-conditioned algebraic systems, which are difficult to invert. In this talk, we present recent theoretical and computational progress that demonstrate a scalable 3D XMHD solver (i.e., $CPU \sim N$, with $N$ the number of degrees of freedom). The approach is based on Newton-Krylov methods, which are preconditioned for efficiency. The preconditioning stage admits suitable approximations without compromising the quality of the overall solution. In this work, we employ optimal ($CPU \sim N$) multilevel methods on a parabolized XMHD formulation, which renders the whole algorithm scalable. The (crucial) parabolization step is required to render XMHD multilevel-friendly. Algebraically, the parabolization step can be interpreted as a Schur factorization of the Jacobian matrix, thereby providing a solid foundation for the current (and future extensions of the) approach. We will build towards 3D extended MHD\footnote{L. Chac\'on, \emph{Comput. Phys. Comm.}, {\bf 163} (3), 143-171 (2004)}$^,$\footnote{L. Chac\'on et al., {\em 33rd EPS Conf. Plasma Physics}, Rome, Italy, 2006} by discussing earlier algorithmic breakthroughs in 2D reduced MHD\footnote{L. Chac\'on et al., {\em J. Comput. Phys}. {\bf 178} (1), 15- 36 (2002)} and 2D Hall MHD.\footnote{L. Chac\'on et al., {\em J. Comput. Phys.}, {\bf 188} (2), 573-592 (2003)}