{'fraction_non_alphanumeric': 0.03019300221350627, 'fraction_numerical': 0.013611121898178712, 'mean_word_length': 4.947228637413395, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Monolithic Algebraic Multigrid Framework for Multiphysics Applications with Examples from Resistive MHD

A multigrid framework is described for multiphysics applications. The framework allows one to construct, adapt, and tailor a monolithic multigrid methodology to different linear systems coming from discretized partial differential equations. The main idea centers on developing multigrid components in a blocked fashion where each block corresponds to separate sets of physical unknowns and equations within the larger discretization matrix. Once defined, these components are ultimately assembled into a monolithic multigrid solver for the entire system. We demonstrate the potential of the framework by applying it to representative linear solution sub-problems arising from resistive MHD.


Introduction & Motivation.
Multigrid methods are among the fastest most scalable techniques for solving the sparse linear systems that arise from the discretization of many different systems of partial differential equations (PDEs) [59]. The basic multigrid idea is to accelerate convergence to a linear system's solution by performing relaxation (i.e., simple iterative techniques) on a hierarchy of different resolution systems. Algebaic multigrid methods (AMG) are popular as they build the hierarchy automatically requiring little effort from the application developer. While algebraic multigrid methods have been successfully applied to many complex problems, further developments are needed to more robustly adapt them to multiphysics PDE systems. In fact, the original AMG research was focused on scalar PDEs, and most of the theory relies on the linear system being symmetric positive definite.
Multigrid's rapid convergence relies on constructing its components such that relaxation sweeps applied to the different fidelity versions are complementary to each other. That is, errors not easily damped by relaxation on one version can be damped effectively by relaxation on another version. However, constructing AMG components with desirable complementary properties can be quite complicated for PDE systems when the coupling between different physical unknowns is significant. AMG challenges arise in designing effective relaxation procedures as well as in designing the grid transfer operators used to project the high fidelity version of the matrix problem to lower resolutions. A discussion of AMG issues for PDE systems and some different approaches can be found in [25]. While there are important AMG themes when considering PDE systems, the best AMG choices are typically application dependent. For this reason, it is critical that an AMG framework can support the implementation and exploration of a wide range of AMG adaptations for complex multiphysics PDE systems. In this paper, we describe such a framework that has been implemented in the MueLu package [4,3], found within the Trilinos libraries [31]. We demonstrate its utility in the context of solving difficult linear systems associated with magnetohydrodynamics (MHD) simulations. In particular, we highlight how the framework allows one to adapt the method to different scenarios. In one case, a specialized block relaxation method is devised that has significant computational advantages over a more black-box relaxation technique. In another case, we show how to adapt the solver to address situations where different finite element basis functions are used to represent the different physical fields of the MHD system.
The resistive magnetohydrodynamics (MHD) model provides a base-level continuum description of the dynamics of conducting fluids in the presence of electromagnetic fields. This model is useful in the context of plasma physics systems coming from both naturally occurring dynamics (e.g. astrophysics and planetary dynamos) and man-made systems (e.g. magnetic confinement fusion) [27]. The governing partial differential equations for the resistive MHD model consist of conservation of mass, momentum and energy augmented by the low-frequency Maxwell's equations. These systems are strongly coupled, highly nonlinear and characterized by coupled physical phenomena that span a very large range of length-and time-scales. These characteristics make the scalable, robust, accurate, and efficient computational solution of these systems extremely challenging. From this point of view, fully-implicit formulations, coupled with effective robust nonlinear/linear iterative solution methods, become attractive, as they have the potential to provide stable, higher-order time-integration of these complex multiphysics systems. These methods can follow the dynamical timescales of interest, as opposed to time-scales determined by fast normal modes that severely impact numerical stability, but do not control temporal order-of-accuracy. The existence of fast normal modes in MHD make the associated algebraic systems very stiff [27], and thus difficult to solve iteratively using modern iterative techniques [9,52,34,50]. Despite these challenges, much effort has been invested in the last decade towards the development of more fully-implicit, nonlinearly coupled solution methods for MHD, with the goal of enhancing robustness, accuracy, and efficiency (see e.g. [37,58,10,11,54,8,38,9,52,53,34,50,15]). To efficiently address the unique linear systems challenge associated with MHD simulations, special solver adaptations must be considered to leverage the underlying character of the equations. The framework described in this paper is intended to facilitate a range of different solver adaptations that might be appropriate for different MHD scenarios.
In Section 2 we introduce the MHD equations and the accompanying discrete systems of equations. In Section 3 we give an overview of the algebraic multigrid method. We discuss the implementation of these algebraic multigrid methods to multiphysics PDE systems in a truly monolithic manner through the use of blocked operators in Section 4. We demonstrate the numerical and computational performance benefits of this approach on various test problems and present the results in Section 5. Finally we end with concluding remarks in Section 6.
2. MHD equations. The model of interest for this paper is the 3D resistive isothermal MHD equations including dissipative terms for the momentum and magnetic induction equations [27]. The system of equations is: Here u is the plasma center of mass velocity; ρ is the mass density; T is the (total) stress tensor; B is the magnetic induction (here after also termed the magnetic field) that is subject to the divergence-free involution ∇ · B = 0. The associated current, J, is obtained from Ampère's law as J = 1 µ0 ∇ × B. In its simplest form, the resistive MHD equations are completed with definition of the total stress tensor, T = T + T M composed of the fluid and magnetic stress tensors, T, T M that are given by, Here, p is the plasma pressure and I is the identity tensor. The physical parameters in this model are the plasma viscosity, µ, the resistivity, η, and the magnetic permeability of free space, µ 0 . Finally, for convenience of notation, a tensor induction flux, F, that is defined as Satisfying the solenoidal involution ∇ · B = 0 in the discrete representation to machine precision is a topic of considerable interest in both structured and unstructured finite-volume and unstructured finite-element contexts (see e.g. [57,18,8]). In the formulation discussed in this study, a scalar Lagrange multiplier (ψ in Eq. (5)) is introduced into the induction equation that enforces the solenoidal involution as a divergence free constraint on the magnetic field. This procedure is common in both finite volume (see e.g. [18,57,8]) and finite element methods (see e.g. [12,13,2,50]). This paper focuses on the incompressible limit of this system, i.e., ∇ · u = 0. This limit is useful in the modeling of various applications such as low-Lundquist-number liquid-metal MHD flows [42,17], and high-Lundquist-number, large-guide-field fusion plasmas [56,30,21]. This limit is characteristic of low flow-Mach-number applications for compressible systems as well, and is the most challenging algorithmically because of the presence of the elliptic incompressibility constraint. However, the stabilized FE formulation that is presented, the strongly-coupled Newton-Krylov nonlinear iterative solvers, and the fully-coupled algebraic multilevel preconditioners also work in the variable density low-Mach-number compressible case. Together the incompressibility constraint with the solenoid involution, enforced as a constraint, produces a dual saddle point structure for the systems of equations [50].
The 3D resistive MHD equations in residual form with the introduction of the scalar Lagrange multiplier and the incompressibility assumption are given by Continuity Constraint (6b) The spatial discretization of the governing equations is based on the VMS FE method [32,20]. The semi-discretized system is integrated in time with a method of lines approach based on BDF schemes. The weak form of the VMS / Stabilized FE formulation for the resistive MHD equation (see Eq. (6)) is given by whereτ i are the stabilization parameters. Here [w h , q h , C h , s h ] are the FE weighting functions for the velocity, pressure, magnetic field and the Lagrange multiplier respectively. The sum e indicates the integrals are taken only over element interiors Ω e and integration by parts is not performed. A full development and examination of this formulation is presented in [50].

Brief
Overview of Discrete Systems of Equations. To provide context to the solution methods that follow later, we present a brief discussion of the equation structure. Here, the focus is on the structure of VMS terms generated by the induction equation and in the enforcement of the solenoidal constraint through the Lagrange multiplier, ψ. Specifically, the VMS weak form of the solenoidal constraint in expanded form, is detailed in [50] as This expression includes a weak Laplacian operator acting on the Lagrange multiplier This term is the analogue of the weak pressure Laplacian L P = e Ωe ρτ m ∇Φ · ∇pdΩ appearing in the total mass continuity equation (see the general discussion for stabilized FE CFD in [20] and our previous development in [52,50] in the context of 2D MHD). These VMS operators are critical in the elimination of oscillatory modes from the null space of the resistive MHD saddle point system for both (u, p) and (B, ψ) and allow equal-order interpolation of all the unknowns (see [20] for incompressible CFD and [12,13,2] for resistive MHD and a coercivity proof of stability). The final term of the stabilized form of the magnetics induction equation (7c) involves a weak divergence-type operator (substitute r ψ ). This term adds to the stability of the VMS form [12,13,2,50] and also enhances the ability to iteratively invert the magnetic induction sub-block in the Jacobain matrix by physics-based and approximate block factorization methods [22,14].
A finite element (FE) discretization of the stabilized equations gives rise to a system of coupled, nonlinear, non-symmetric algebraic equations, the numerical solution of which can be very challenging. These equations are linearized using an inexact form of Newton's method leading to a block system with the following form The block matrix, J u , corresponds to the discrete transient, convection, diffusion and stress terms acting on the unknowns δû; the matrix, G, corresponds to the discrete gradient operator; D, the discrete representation of the continuity equation terms with velocity (note for a true incompressible flow this would be the divergence operator denoted G T ); the block matrix, J I , corresponds to the discrete transient, convection, diffusion terms acting on magnetic induction, and the matrices, L P , L ψ , are the stabilization Laplacian's mentioned above. The right hand side vectors contain the residuals for Newton's method. We solve for the solution increments δû and δp for the velocity and pressure as well as for δB and δψ which represent the magnetics field and the Lagrange multipliers. The L P and L ψ operators help facilitate the solution of the linear systems with a number of algebraic and domain decomposition type preconditioners that rely on non-pivoting ILU factorization, Jacobi relaxation or Gauss-Seidel as sub-domain solvers [49,50,51].
The difficulty of producing robust and efficient preconditioners to (10) has motivated many different types of decoupled solution methods. Often, transient schemes combine semi-implicit methods with fractional-step (operator splitting) approaches or use fully-decoupled solution strategies [58,37,1,44,35,28,47,29,55,33,46]. In these cases, the motivation is to reduce memory usage and to produce a simplified equation set for which efficient solution strategies already exist. Unfortunately, these simplifications place significant limitations on the broad applicability of these methods. A detailed presentation of the characteristics of different linear and nonlinear solution strategies is beyond our current scope. Here, we wish to highlight that our approach of fully-coupling the resistive MHD PDEs in the nonlinear solver preserves the inherently strong coupling of the physics with the goal to produce a more robust solution methodology [49,52,50]. Preservation of this strong coupling, however, places a significant burden on the linear solution procedure.
3. Multigrid methods. Multigrid methods are based on the fact that many simple iterative methods are effective at eliminating high-frequency error components relative to the "mesh" resolution used for discretization. The basic multigrid idea is to introduce a hierarchy of discrete approximations to the PDE problem employing different resolution meshes. The newly introduced coarser approximations are used to accelerate the solution process on the finest mesh as lower-frequency error components on the finest mesh will appear to be high-frequency (relative to the grid resolution) on a particular mesh within the hierarchy. That is, different error components are efficiently reduced by essentially applying a simple iterative method to the appropriate resolution approximation. Multigrid methods generally come in two varieties: geometric multigrid (GMG) and algebraic multigrid (AMG). Typically with GMG , applications supply a hierarchy of meshes, discrete PDE operators, interpolation operators to transfer solutions from coarser resolutions to finer ones, and restriction operators to transfer finer resolution residuals to coarser levels. In geometric multigrid, inter-grid transfers are based on geometric relationships, such as using linear interpolation to define a finer level approximation from a coarse solution. With AMG methods, coarse level information and inter-grid transfers are developed automatically by analyzing the supplied fine level discretization matrix. This normally involves a combination of graph heuristics to coarsen followed by some approximation algorithm to develop inter-grid transfers with the aim of more accurately transferring information associated with small magnitude eigenvalues of the discrete operator, which is often low frequency. This paper focuses on AMG as it can be more easily adapted to complex application domains by non-multigrid scientists. Figure 1 depicts what is referred to as a multigrid V-cycle to solve a linear system A u = b . Subscripts distinguish between different resolutions. P +1→ interpolates from level +1 to level . R → +1 restricts from level to level +1. A is the discrete problem on level and for coarse levels is defined by a Petrov-Galerkin projection S pre and S post denote a basic iterative scheme (e.g., Gauss-Seidel) that is applied to damp or relax some error components. The overall efficiency is governed by the interplay of the two main multigrid ingredients: inter-grid transfer operators and the smoothing methods. One can either use several sweeps with the multigrid V-cycle as a standalone solver, or one can apply a V-cycle sweep as a preconditioner within an outer Krylov based linear iterative solver. For a general overview on multigrid methods the reader is referred to [7,59] and the references therein.
For PDE and multiphysics systems, applying AMG to the entire PDE system (i.e., monolithic multigrid) can be problematic, especially when the coupling between different solution types (e.g., pressures and velocities) is strong. Classical simple iterative methods may not necessarily reduce all oscillatory error components and might even amplify some oscillatory error components. In fact, methods requiring the inversion of the matrix diagonal (e.g., the Jacobi iteration) are not even well defined when applied to incompressible fluid formulations that give rise to zeros on the matrix diagonal. Further, many standard AMG algorithms for defining inter-grid transfers might lead to transfers that do not accurately preserve smooth functions. For example, methods such as smoothed aggregation rely on a Jacobi-like step to generate smooth inter-grid transfer basis functions. This Jacobi step, however, will obviously not generate smoother basis functions when the Jacobi iteration is not well defined (or when the matrix diagonal is small in a relative sense). In even worse situations, the AMG software may not be applicable to the PDE system when different FE basis functions are used to represent different fields within the system. This is because most AMG codes employ a simple technique to address PDE systems. This technique relies on the different equations to be effectively discretized in a similar way. In particular, AMG coarsening algorithms (or for our approach aggregation algorithms) are often applied by first grouping all DoFs at each mesh node together and then applying the coarsening algorithms to the graph induced from the block matrix. An advantage to this approach is that all unknowns at mesh points are coarsened in identically the same fashion. However, the approach cannot be applied when the number of DoFs associated with different fields varies or if all DoFs are not co-located, e.g., when using quadratic basis functions to represent velocities while pressures employ only linear basis functions. In previous versions of our AMG software, the only alternative to this simple PDE system technique would be to completely ignore the multiphysics coupling and effectively treat the entire system as if it is a scalar PDE, which almost always leads to disastrous convergence rates.

Truly monolithic block multigrid for the MHD equations.
In this section we propose a multigrid method for multiphysics systems that can be represented by block matrices allowing one to adapt both the relaxation algorithms and the grid transfer construction algorithms to the structure of the system.

Block matrices and multigrid for PDE systems.
As illustrated in the previous MHD discussion, PDE systems are often represented by block matrices as in, for example, equation (10). More generally, block systems can be written as where each component of the vector x i is a field in the multiphysics PDE and the sub-matrices A ij are approximations to operators in the governing equations. One preconditioning approach to PDE systems follows a so-called physics-based strategy (see Figure 2a). These techniques can be viewed as approximate block factorizations (involving Schur complement approximations) to the block matrix equations (11). The factorizations are usually constructed based on the underlying physics. Here, different AMG V cycle sweeps are used to approximate the different sub-matrix inverses that appear within the approximate block factors. As the sub-matrices correspond to single physics or scalar PDE operators, application-specific modifications to the multigrid algorithm are often not necessary. This makes the physics-based strategy particularly easy to implement as one can leverage ready-to-use multigrid packages. Several methods that follow this popular strategy are described in [23] and the references therein. While the physics-based approach has some practical advantages, the efficiency of the preconditioner relies heavily on how well the coupling between different equations within the PDE is approximated by the block factorization.
In this paper, we instead consider a monolithic multigrid alternative to a physicsbased approach. A monolithic scheme applies a multigrid algorithm to the entire block PDE system and so the AMG scheme effectively develops a hierarchy of block PDE matrices associated with different resolutions. Figure 2 graphically illustrates the two contrasting approaches highlighting the key potential advantage to a monolithic approach. In particular, the cross-coupling defined by A 01 and A 10 is explicitly Fig. 2: Two multigrid approaches to address multiphysics applications.
(a) Physics-based approach with multiple multigrid approximations to sub-matrix inverses within an approximate block factorization. Example with 4 and 2 multigrid levels for A00 and A11.
represented on all multigrid hierarchy levels with a monolithic approach as opposed to the physics-based scheme that attempts to represent coupling within an approximate Schur complement (whose inverse might then be approximated via multigrid sweeps). That is, a physics-based approach relies heavily on being able to develop effective Schur complements, which can be non-trivial for complex applications. On the other hand, a monolithic approach introduces its own set of application-specific mathematical challenges such as the construction of relaxation procedures for monolithic systems and the development of coarsening schemes for different fields within a multiphysics system. The design of efficient multiphysics preconditioners often requires one to make use of the specific knowledge about the block structure, the mathematical models of the underlying physics, and the problem-specific coupling of the equations. The mathematical challenges of multiphysics systems are often further compounded by non-trivial software challenges. Unfortunately, most AMG packages cannot be customized to particular multi-physics scenarios without having an in-depth knowledge of the AMG software. It is for this reason that many prefer a physics-based approach.
In the next sub-sections, we propose a monolithic algorithm for the MHD equations that can be easily adopted and customized using the MueLu package within the Trilinos framework [4,3]. Though we focus on a concrete MHD case, we cannot emphasize enough the importance of the software's generality in facilitating a monolithic approach for those with limited knowledge of the multigrid package internals.
Algebraic representation of the MHD problem. In adapting a monolithic multigrid strategy to the MHD equations, a natural approach would be to interpret the system (10) as a 2 × 2 block system where the Navier-Stokes equations are separated from the Maxwell equations. That is, we treat the Navier-Stokes part and the Maxwell part as separate entities that are coupled by the off-diagonal blocks as shown in the 2 × 2 block representation of Figure 3. In this way, we can leverage existing ideas/solvers for the Navier-Stokes equations and for the Maxwell equations. In making this 2 × 2 decomposition, we are effectively emphasizing the significance of the coupling between fields within the Navier-Stokes block and within the Maxwell block as compared to the coupling between Navier-Stokes unknowns and Maxwell unknowns. This is due to the importance of the coupling constraint equations (e.g., incompressibility conditions involving velocities or contact constraints [61]) to the associated evolution equations. These constraint equations often give rise to saddlepoint like block systems. Of course, there are physical situations where the coupling between the Navier-Stokes equations and the Maxwell equations is quite significant and so a block 2 × 2 decomposition might be less appropriate. There might also be situations where coupling relationships are more complex. If, for example, the MHD equations are embedded in another larger more complex multiphysics problem, then one might need to consider a hierarchy of of coupling configurations, which might require different arrangements/blocking of multigrid ingredients. Given the problem specific nature of multiphysics preconditioning, our emphasis here is on the importance of a flexible software framework to facilitating different types of blocking within the preconditioner.
For the remainder of the paper the block notation is used representing the corresponding blocks from Figure 3. The velocity and pressure increments δû and δp are grouped in x 0 and the Maxwell information is represented by x 1 , respectively. In a similar way, the block notation for the residual vector is adopted.

Monolithic multigrid ingredients for volume-coupled problems.
The multiphysics solver that we propose is generally applicable to volume-coupled problems. Volume coupled means that the different physics blocks are defined on the same domain. Volume-coupled examples include, e.g, thermo-structure-interaction (TSI) problems (see [16]) or in our case the MHD equations. Specifically, within our MHD formulation all physics equations (i.e., Navier-Stokes and Maxwell parts) are defined throughout the entire domain. This is in contrast to interface-coupled problems such as fluid-structure interaction (FSI) applications (see [26,19,40,36]) or structural contact problems (see [61,62]) where different equation sets are valid over distinct domains that are only coupled through a common interface. From a multigrid perspective special interface coarsening methods are necessary for interface-coupled applications, which is not the focus of this paper.
Two multigrid ingredients must be specified to fully define the monolithic solver: the inter-grid transfer operators and the relaxation or smoother procedures.

4.3.1.
Inter-grid transfers for the MHD system. Following the 2 × 2 decomposition of the MHD system from equation (12), we consider rectangular block diagonal inter-grid transfer operators (13) R i→i+1 = R 00 R 11 and P i+1→i = P 00 P 11 for restriction and prolongation between multigrid levels respectively. The basic idea uses the MueLu multigrid package to produce grid transfers for the Navier-Stokes equations and the Maxwell equations and then leverages MueLu's flexibility to combine these grid transfer operators into a composite block diagonal operator. The block perspective allows us to use completely separate invocations of the multigrid package to build individual components (e.g., P 00 and P 11 ) and then combine or compose them together. As discussed earlier, we rely on underlying core AMG kernels that are applicable to either a scalar PDE or a PDE system with co-located unknowns. Mixed finite element schemes may not normally satisfy this co-located requirement, so the ability to separately invoke these core components to produce P 00 and P 11 alleviates this restriction, allowing us to apply monolithic AMG to a wider class of PDE systems. That is, a mixed basis function discretization can be approached without having to erroneously treat the entire system as a scalar PDE. We demonstrate this capability at the end of Section 5 through a mixed formulation utilizing Q2/Q2 VMS for the Navier-Stokes degrees of freedom and Q1/Q1 VMS for the Maxwell degrees of freedom.
In the case where the unknowns between blocks are co-located (e.g., Q1/Q1 VMS for both Navier-Stokes and Maxwell) we have the option to correlate the grid transfer construction by having the multigrid invocations share the same aggregates (or coarsening definition) as depicted in Figure 4. In this way, we guarantee that there is a one-to-one relationship of the coarse Maxwell degrees of freedoms and the associated Navier-Stokes degrees of freedom. That is, we obtain the same coarsening rate for both, and so the ratio between Navier-Stokes and Maxwell degrees of freedom is constant on all multigrid levels. It should be noted that this is relatively straight-forward when all equations are defined on the same mesh using a first-order nodal finite element discretization method. Thus, there are 8 degrees of freedom at each mesh node (four associated with fluid flow and four associated with electromagnetics).
Within our aggregation approach, we first build standard aggregates using the graph of sub-block A 00 on level i as input. Next we build separate aggregates using the graph of sub-block A 11 on level i as input. However, if the degrees of freedom are co-located, instead of building new aggregates we have the option to reconstruct the corresponding aggregates for the Maxwell equations by cloning the Navier-Stokes aggregation information as depicted in Figure 4. Specifically, mesh vertices on a given level are assigned to aggregates A i such that where N denotes the number of mesh vertices on level . Each aggregate A i on level gives rise to one node on level + 1. The A i are formed by applying greedy algorithms to the graph associated with a matrix discretization. Typically, one wants aggregates to be approximately the same size and roughly spherical in shape (for isotropic problems). Piecewise-constant interpolation can then be defined over each aggregate for each solution component. Node discretization of current level.
As discussed later, this type of multigrid adaptation is relatively straight forward with MueLu using an XML input file. In addition to consistently coarsening both fluid and electromagnetic variables, the reuse or cloning of aggregates saves a modest amount of time within the multigrid setup phase. Additional reductions in the setup cost can be achieved if one defines R 11 = R 00 and P 11 = P 00 . Again, this is relatively easy to do within the MueLu framework, though not necessary. In our experiments, piecewise constant interpolation is the basis for all grid transfers. This simple grid transfer choice is more robust for highly convective flows.

4.3.2.
Block smoother for the MHD system. A block Gauss-Seidel (BGS) iteration forms the basis of the multigrid smoother using the notion of blocks already introduced via the 2 × 2 decomposition. A standard block Gauss-Seidel iteration would solve for an entire block of unknowns simultaneously, recompute residuals, and then solve for the other block of equations simultaneously. This corresponds to alternating between a Navier-Stokes sub-block solve and a Maxwell sub-block solve while performing residual updates after each sub-solve. Effectively, the 2 × 2 block perspective emphasizes the coupling within the Navier-Stokes block and within the Maxwell part as each of these blocks correspond to a saddle-point system due to the presence of constraint equations. That is, the difficulties associated with saddlepoint systems will be addressed by an approximate sub-block solve. Further, the block Gauss-Seidel iteration considers the Navier-Stokes block and the Maxwell block equally important as the iteration simply alternates equally between the two subsolves. As the exact solution of each sub-system is computationally expensive, domain decomposition with ILU(0) is used to generate approximate sub-solutions associating one domain with each computing core. While not necessarily an optimal smoother, we have found that ILU(0) is often effective for the saddle-point systems associated with incompressible flow. While ILU(0) could be applied to the entire 2 × 2 matrix, this involves significantly more computation and memory during the setup and apply phases of the smoother. This more expensive process is generally not needed as the coupling between the block sub-systems is often much less important than the coupling within the sub-blocks, though there could be MHD situations with significant crosscoupling between Navier-Stokes and Maxwell that would warrant an ILU smoother applied to the whole system. As one can see from Algorithm 1 we only need approximate inverses of the diagonal blocks A 00 and A 11 . A flexible implementation allows one to choose appropriate local smoothing methods to approximately invert the blocks A 00 and A 11 . The coupling is guaranteed by the off-diagonal blocks A 01 and A 10 in the residual calculations in lines 4 and 8.

4.4.
Software. Even though this paper demonstrates the multigrid approach for certain MHD formulations, it is clear that the concept is more general. Due to the problem-dependent nature of preconditioning multiphysics problems, it is logical to provide a general software framework that helps design application-specific preconditioning methods. The proposed methods in this paper are implemented using the next-generation multigrid framework MueLu from the Trilinos software libraries. In contrast to other publications like [60] which share the same core idea of a general software framework, MueLu is publicly available through the Trilinos library. Furthermore, it is fully embedded in the Trilinos software stack and has full native access to all features provided by the other Trilinos packages. It aims at next-generation HPC platforms and automatically benefits from all performance improvements in the underlying core linear algebra packages. The core design concept of MueLu is based on building blocks that can be combined to construct complex preconditioner layouts. Figure 5 shows an algorithmic layout to build standard aggregation-based transfer operators P and R, the level smoother S and a new coarse level operator A +1 . Further coarse levels follow by recursively applying a multigrid setup procedure to the coarse level operators. Each building block, denoted by a rectangular block in Figure 5, processes certain input data and produces output data which serves as input for the downstream building blocks. For application-specific adaptions it is usually sufficient to replace specific algorithmic components while the majority of building blocks can be reused. Figure 5 represents the process of defining coarse multigrid levels during the setup phase for a single-field multigrid method. Figure 6 shows the algorithmic design of the setup phase to create a 2 × 2 block transfer operator for a volume-coupled problem as discussed earlier. Here, aggregates built from block A 00 are re-used instead of re-building them from build block A 11 5. Experimental results.
Here, B 0 is the strength of the field and δ is a measure of the transition length-scale for application of the field. The inlet, outlet and sides are perfect conductors with B ·n = 0 and E ×n = 0, wheren is the outward facing unit normal vector. Zero Dirichlet boundary conditions are applied on all surfaces for the Lagrange multiplier. The problem is defined by three non-dimensional parameters: the Reynolds number Re = ρuL/µ, the magnetic Reynolds number Re m = µ 0 uL/η, and the Hartmann number Here, u is the maximum x-direction velocity. The parameters in this problem are taken to be u = 1.0, ρ = 1, B 0 = 3.354, µ 0 = 1, η = 1, x on = 4.0, x off = 6.0, and δ = 0.5. The linear solver used was non-restarted GMRES to reach a relative tolerance of 10 −3 . As a reference preconditioner we use a fully-coupled or nonblocked AMG method (FC-AMG), where the relaxation method, Additive Schwarz with overlap one, is applied to the entire system that includes the off-diagonal coupling between the fluids and magnetics (see [41]). For the blocked variant presented in this manuscript (Section 4), we use BGS as a relaxation method (AMG(BGS). The approximate solves for the sub-blocks are handled with Additive Schwarz with overlap one. The coarse grid for both the FC-AMG and AMG(BGS) is solved directly.
In the first study, we investigate the number of iterations and solution time as a function of the block Gauss-Seidel (BGS) damping parameter ω for various viscosities µ ∈ {0.006, 0.007, 0.008, 0.01, 0.02, 0.04}, which effectively sets the range for the nondimensional parameters: 25 ≤ Re ≤ 167, 17 ≤ Ha ≤ 43, and Re m = 1.      Figure 7b displays the accumulated solver timings (setup and iteration phase) of the corresponding preconditioning variants. The timings are averaged over 5 simulations. As one can see from Figure 7, the right choice of the BGS damping parameter is crucial for smaller viscosities. Even though the optimal damping parameter depends on the problem, in practice a damping parameter ω near 0.5 seems to work well. Figure 8 shows the accumulated linear iterations and the solver times for the MHD generator example on a finer 480 × 32 × 32 mesh using 256 processors (Intel Broadwell E5-2695 (2.1 GHz), 8 cluster nodes with 128 GB RAM each, Intel Omni-Path high speed interconnect). One can see in Figure 8a that a higher number of BGS iterations reduces the overall number of linear iterations, but not enough to compensate for the higher costs per iteration (see Figure 8b). For reference, the fully coupled approach is also shown in the back row (labeled as FC-AMG ILU(0)). This option does not require the multiphysics framework, but it generally takes more time and iterations than the multiphysics solver, varying somewhere between 2x and 3x slower.   It should be noted that the scaling of the iteration counts with respect to mesh refinement can be improved by increasing the linear solver tolerance for both the FC-AMG and AMG(BGS) methods. In experiments not shown here using the same solver settings but tightening the tolerance to 10 −6 , the FC-AMG total number of linear iterations increased by less than 10% while the AMG(BGS) iterations decreased for the same two meshes considered in Figure 7 and 8. While the scaling is better, the overall solution time with these tighter tolerances is longer and therefore it is not considered further in this paper since for engineering applications our primary goal is to minimize the solver time and not the iteration count.

Hydromagnetic Kevin-Helmholtz (HMKH).
A hydromagnetic Kelvin-Helmholtz unstable shear layer problem is a configuration used to study magnetic reconnection [50] and is posed in a domain of [0, 4] × [−2, 2] × [0, 2]. It is described by an initial condition defined by two counter flowing conducting fluid streams with constant velocities u(x, y > 0, z, 0) = (5, 0, 0) and u(x, y < 0, z, 0) = (−5, 0, 0) and a Harris sheet sheared magnetic field configuration given by (14) B(x, y, z, 0) = (0, B 0 tanh(y/δ), 0) . For these non-dimensional parameters, the shear layer is Kelvin-Helmholtz unstable and forms a vortex sheet that evolves with time and undergoes thin current sheet formation, vortex rollup and merging. Figure 9 shows the unstable shear layer evolving from smaller vortices to a larger vortex after about half the total runtime of the simulation.
In this numerical study, we perform transient simulations of the problem with CFL max = 0.25 and CFL max = 0.5, and study the behavior of the nonlinear-and linear solver. In particular, we compare the number of iterations and timings of the non-restarted GMRES solver using blocked AMG(BGS) as a preconditioner. The fully-coupled or non-blocked AMG method (FC-AMG) is used as a reference preconditioner. The linear solver tolerance is set to 10 −3 . Figures 10 and 11 show the solver performance for the different preconditioning strategies over a time sequence for the HMKH problem with a maximum CFL number of 0.25 and 0.5 respectively. We ran the problem both on a 64 × 32 × 16 mesh on 32 processors (Intel Broadwell E5-2695 (2.1 GHz), 1 cluster node with 128 GB RAM) and a 128 × 64 × 32 mesh on 256 processors (Intel Broadwell E5-2695 (2.1 GHz), 8 cluster nodes with 128 GB RAM each, Intel Omni-Path high speed interconnect). Generally, the FC-AMG method needs the least number of iterations, whereas the blocked AMG(BGS) variant with only 1 coupling iteration needs the highest number of iterations. Increasing the number of BGS coupling iterations reduces the number of linear iterations getting closer to the reference method. However, looking at the linear solver time (setup and iteration phase), we see the opposite picture. The FC-AMG method is the slowest, and the blocked AMG(1 BGS) is the most time efficient method. Comparing the solver behavior for the different meshes there is a slight increase in the linear iteration count for the finer meshes. So, while the weak scaling is not optimal, the iteration growth with problem size is mild.
Numerical results are further summarized in Tables 1 and 2. The first column denotes the average number of nonlinear iterations per timestep for the full simulation. One can see that the number of BGS coupling iterations has some influence on the nonlinear solver, even though there is no clear trend. The second column represents the average number of linear iterations per nonlinear iteration. As one would expect, a higher number of BGS coupling iterations reduces the number of linear iterations necessary to solve the problem. The next three columns give the average setup time, the average solve time and the average overall time for the linear solver per nonlinear iteration. The AMG(BGS) variants have a clear advantage in the setup costs, but the iteration costs are higher. That is, fully coupled ILU(0) that includes Navier-Stokes and Maxwell coupling solver provides some convergence benefits, but at the cost of a very large setup time. The last three columns show the absolute setup, solve and overall solver time for finishing the simulation.

Island coalescence.
The island coalescence problem is a prototype problem used to study magnetic reconnection. Here the initial equilibrium is described by two 3D current tubes (islands in the 2D plane) embedded in a Harris current sheet (as in the HMKH problem of Section 5.2) in a [−1, 1] × [−1, 1] × [−1, 1] domain [24,9]. The initial condition for the island coalescence problem consists of zero fluid velocities (u 0 = 0), zero Lagrange multipliers (ψ = 0) and a Fadeev magnetic equilibrium [24,9] that defines the magnetic field B and the fluid pressure p. More details of this setup can be found in [50]. The structure of this equilibrium is presented in Figure 12a with an iso-surface of p and iso-lines of B at z = 0. The combined magnetic field of the two islands produces Lorentz forces that pull the islands together. The dynamics of island coalescence changes as a function of resistivity. For larger resistivities, the xand o-points monotonically approach each other. For low resistivities, fluid-plasma pressure builds up as the islands approach and a sloshing or bouncing of the o-point position is encountered that leads to lower reconnection rates (for more details on the physics see e.g. [5]). Figure 12 shows different stages of the reconnection event. Clearly evident is the formation of the x-point in the intersecting planes between the islands (see images at t = 4), the development of thin current sheets at that same x-point location (and the corresponding 3D surface), and the movement of the center of the tubes (island o-points) towards the x-point [5,39].
In this study, we have taken ρ = 1, µ = η = 10 −3 , µ 0 = 1 and, using the spacing of the o-points, we have L = 1, resulting in Re = Re m = 10 3 . As in [39] these choices imply that the resistivity η = 1/S, where S is the Lundquist number and is defined as S = µ 0 Lu A /η, where u A is the Alfven velocity. We preform transient simulations of the problem with timestep sizes of ∆t ∈ {0.05, 0.025, 0.0125}. The mesh sizes used were 32 × 32 × 32, 64 × 64 × 64, and 128 × 128 × 128 and were run on 8, 64, and 512 processors respectively (Intel Broadwell E5-2695 (2.1 GHz), 2, 16, and 128 cluster nodes with 128 GB RAM each, Intel Omni-Path high speed interconnect). This provides simulations with a CFL ranging from 1.6 to 12.8. We compare the number  of iterations and timings of the non-restarted GMRES solver, using a relative linear solve tolerance of ε = 10 −3 , when combined with different preconditioning strategies, including the fully-coupled or non-blocked AMG method (FC-AMG) as reference and the AMG(BGS) variants. Figure 13 shows the solver performance for different preconditioning strategies over a time sequence for the island coalesce problem with CFL number of 3.2. As with the HMKH example, we see that while the AMG(BGS) variants require more iterations than the FC-AMG reference, the cost per iteration is low enough that the Fig. 12: Structure of the current tubes in 3D island coalescence problem with S = 2 × 10 4 for the initial condition, and four times in the evolution of the problem t = 2, 3, 4. The 3D current tubes have bent in the z-direction and form current sheets. 5.4. Mixed finite elements. Next we illustrate a formulation for which the hydrodynamics and electromagnetics systems are discretized by differing FE spaces. In this example we consider Q2/Q2 VMS for the hydrodynamics (saddle point and convective stabilization) and Q1/Q1 VMS for the induction (electromagnetics) systems (saddle point and convective stabilization). In this problem the difference in order-of-accuracy is motivated by the desire to minimize the overall computational time while still maintaining higher accuracy for the MHD simulation in appropriate applications. For example, when a liquid metal is the conducting fluid in an MHD generator, the flow Reynolds number can be significantly higher than the corresponding magnetic Reynolds number due to the very high magnetic diffusivity of liquid metals. In general the low magnetic Reynolds number is indicative of diffusive dominated transport for the magnetics in the liquid metal. Thus, a mixed discretization with a lower-order approximation for the induction subsystem may be appropriate. Other cases that employ disparate discretizations for hydrodynamics and magnetics would be various forms of structure preserving methods where, for example, nodal FE are employed for flow variables and face or even edge FE are used for the magnetic field (see e.g. [43,6,48,45]).
Here, we again consider the MHD Generator problem from Section 5.1 with the modest intention of demonstrating the ability of our proposed methods to handle disparate discretizations, which in this case are Q2/Q2 VMS for the hydrodynamics   (c) 128 × 128 × 128 mesh, ∆t = 0.0125s, 512 processors and Q1/Q1 VMS for the induction (electromagnetics) systems. The system is difficult to approach through standard fully coupled AMG methods due to the mixed FE spaces with DoFs that are no longer co-located. The blocked approach outlined in this manuscript allows for the separate construction of aggregates for the hydrodynamics block and the electromagnetics block. The monolithic multigrid hierarchy then naturally provides us with coupling between the blocks on all levels of the hierarchy.  The study is carried out on the same set of physical, geometrical and solver parameters as in Section 5.1, with varying viscosities µ ∈ {0.006, 0.01, 0.04}. The relaxation method is a blocked Gauss-Seidel with a damping parameter of 0.6. For the sub-block solves, a single iteration of Additive Schwarz with an overlap of 1 was used to generate approximate sub-block solutions. The current implementation we are using lacks parallel load rebalancing, which is problematic on higher core counts. To circumvent this issue, the maximum number of AMG levels was capped at 4 levels, as further coarsening of the 2048 processor case requires rebalancing. The coarsest level problem, is still relatively large in the 2048 processor case (16,384 rows after 3 levels of refinement, or 8 DoFs per processor), so the coarse level solve is handled with an iteration of the smoother instead of a direct solve. We explore the weak scaling of the method in Table 7, showing iterations and timings for various preconditioner configurations. For comparison, we also consider an additive Schwarz domain decomposition method (DD-Schwarz) with overlap 1 for the entire 2 × 2 block system as a preconditioner. The only other possible AMG option without using the multiphysics framework considers the entire system as a scalar PDE. This non-blocked AMG approach performs so poorly that results are not given. The use of Blocked AMG provides a significant reduction in setup time over the monolithic additive Schwarz, as the block off-diagonal terms are no longer considered in the factorization. While the number of linear iterations does degrade as the problem size increases, the ability to apply Blocked AMG to this mixed FE space problem provides a significant linear solve time speed up compared to the use of Additive Schwarz on the entire 2×2 block system. Additional work is needed to better understand smoother and aggregation choices in the mixed FE case. For example, the increase in iterations as problem size increases for high viscosity indicates potential inefficiencies with our handling of the Q2 hydrodynamics problem, as the low viscosity case has a more consistent iteration count as the problem size increases.
Future work will consider correlated coarsening algorithms where the Q1/Q1 ag-gregates influence the Q2/Q2 aggregation scheme. In this example there is some partial overlap in the location of DoFs on the mesh. Some mesh nodes have 8 DoFs, corresponding to four DoFs for hydrodynamics and four DoFs for electromagnetics, Others only have 4 DoFs, corresponding only to the hydrodynamics. A natural extension is to force the aggregation scheme to preserve this partial co-location aspect of the hydrodynamics and electromagnetics Dofs. In future work, we plan to incorporate some ability to partially share aggregation information between the two AMG invocations. In this case, one might share the aggregate root (or central) vertices generated during the AMG invocation for electomagnetics. These root vertices could then be used to construct an initial set of aggregates for the hydrodynamics. As there are more hydrodynamic unknowns, many hydrodynamic unknowns might remain unaggregated and so further aggregation would be needed to complete the set of aggregates for the hydrodynamics.

Conclusion.
A new framework for developing multiphysics multigrid preconditioners is developed and demonstrated on a number of MHD problems. The key idea is to develop the multigrid components in a block fashion that mirrors the blocks in a multiphysics system. Our approach has been to develop block smoothers and apply them to a multigrid hierarchy constructed using block restriction/prolongation operators. In many cases, the blocked multiphysics multigrid hierarchy allows for faster solution times than a non-blocked approach. For mixed spatial discretizations, the multiphysics framework provides the only genuine avenue to leverage pre-existing multigrid software to produce a monolithic multigrid preconditioner. Here, the AMG engine is invoked multiple times for different sub-blocks and the resulting individual grid transfers are combined into one composite operator that can be employed in a monolithic AMG fashion. Ultimately, the run time benefits are much greater because no other multilevel preconditioning option is available. While this paper has focused on specific examples and MHD, the goal of the framework is to be able to easily construct, adapt, and tailor different monolithic multigrid preconditioners to various PDE systems.