**1.1 Scope of this chapter**

The scientific literature has numerous instances of methods and techniques to solve the MHD system of equations. To limit the scope of this chapter, we focus our discussion to single fluid resistive and ideal MHD. Although single fluid resistive (or ideal) MHD is in a sense the simplest fluid model for a plasma, these equations constitute a system of nonlinear partial differential equations, and hence pose many interesting challenges for numerical methods and simulations. In particular, there is a vast amount of literature devoted to numerical methods and simulations of resistive MHD wherein the time stepping method is explicit or semi-implicit. For example, in simulating MHD flows with shocks, shock-capturing methods from hydrodynamics have been tailored to MHD and have been very successfully used (see for example Reference Samtaney et al. (2005)). Such aforementioned shock-capturing methods almost exclusively employ explicit time stepping. This is entirely sensible given that the flow speeds are of the same order as, or exceed the fast wave speeds. In several physical situations, the diffusive time scales are much larger than the advective time scale. In these cases, the Lundquist number is large (*S* >> 1) and the diffusion terms are usually much smaller than the hyperbolic or wave-dominated terms in the equations. Usually the diffusion terms become important in thin boundary layers or thin current sheets within the physical domain. We are interested in computing such flows but with the additional constraint that

where *R*(*U*) ≡∇· *F*(*U*) −∇· *Fd*(*U*), and the solution vector *U* ≡ *U*(**x**, *t*) is

*<sup>F</sup>*(*U*) =

*Fd*(*U*) =

and the stress-strain tensor relation

**2. Brief survey of implicit methods for MHD**

∇ · *B* = 0.

**2.1 Early approaches**

*ρu* , *ρuu* +

*uB* − *Bu* ,

0 , *Re*−1**ø**¯¯ , *S*−<sup>1</sup>

*Re*−1**ø**¯¯ · *<sup>u</sup>* <sup>+</sup> *<sup>γ</sup>*

characteristic length scale, *L*, and the Alfvén speed *U*<sup>0</sup> = *B*0/

 *p* + 1 2 *B* · *B* ¯¯ **I** − *BB* ,

*e* + *p* +

*γ* − 1

and the hyperbolic flux vector *F*(*U*) and the diffusive fluxes *Fd*(*U*) are given by

1 2 *B* · *B* 

*κ Re Pr* <sup>∇</sup>*<sup>T</sup>* <sup>+</sup> *<sup>η</sup>*

*<sup>η</sup>*∇*<sup>B</sup>* <sup>−</sup> *<sup>η</sup>*(∇*B*)*<sup>T</sup>*

In the above equations *ρ* is the density, *u* is the velocity, *B* is the magnetic field, *p* and *T* are the pressure and temperature respectively, and *e* is the total energy per unit volume of the plasma. The plasma properties are the resistivity *η*, the thermal conductivity *κ*, and the viscosity *μ*, which have been normalized, respectively, by a reference resistivity *ηR*, a reference conductivity *κR*, and a reference viscosity *μR*. The ratio of specific heats is denoted by *γ* and generally fixed at 5/3 in most MHD simulations. The non-dimensional parameters in the above equations are the Reynolds number, defined as *Re* ≡ *ρ*0*U*0*L*/*μR*, the Lundquist number, defined as *S* ≡ *μ*0*U*0*L*/*ηR*, and the Prandtl number, denoted by *Pr*, which is the ratio of momentum to thermal diffusivity. The non-dimensionalization is carried out using a

the characteristic strength of the magnetic field, a reference density, and the permeability of

1 2 *B* · *B*,

 − 2 3 *<sup>μ</sup>*∇ · *<sup>u</sup>*¯¯ **I**.

free space, respectively. The equations are closed by the following equation of state

*<sup>γ</sup>* <sup>−</sup> <sup>1</sup> <sup>+</sup> *<sup>ρ</sup>* 2 *u* · *u* +

<sup>∇</sup>*<sup>u</sup>* + (∇*u*)*<sup>T</sup>*

Finally, a consequence of Faraday's law is that an initially divergence-free magnetic field must lead to a divergence-free magnetic field for all times, which corresponds to the lack of observations of magnetic monopoles in nature. This solenoidal property is expressed as

In the limit of zero resistivity, conductivity and viscosity, the equations of resistive MHD reduced to those of ideal MHD. These equations are similar to those written above with *κ* = *μ* = *η* = 0. Ideal MHD equations are hyperbolic PDEs (although not strictly hyperbolic).

An implicit treatment of the fast magnetosonic wave coupled with arguments of large length scales dominantly in a certain direction allows one to investigate long-time scale phenomena

*<sup>e</sup>* <sup>=</sup> *<sup>p</sup>*

**ø**¯¯ = *μ*  *<sup>U</sup>* <sup>=</sup> {*ρ*, *<sup>ρ</sup>u*, *<sup>B</sup>*,*e*}*<sup>T</sup>*

Implicit Numerical Methods for Magnetohydrodynamics 61

*u* − *B*(*B* · *u*)

 ,

*S* 1 2 *<sup>T</sup>* ,

<sup>∇</sup>(*<sup>B</sup>* · *<sup>B</sup>*) <sup>−</sup> *<sup>B</sup>*(∇*B*)*<sup>T</sup>*

*<sup>T</sup>* .

√*μ*0*ρ*0, where *<sup>B</sup>*0, *<sup>ρ</sup>*0, and *<sup>μ</sup>*<sup>0</sup> are

(2)

the wave speeds are also much larger than the fluid advective speeds. In such cases, implicit time stepping methods are preferred to overcome the stiffness induced by the fast waves and march the flow simulations forward in time at time steps dictated by accuracy rather than stability constraints.

A brief outline of this chapter is as follows. In Section 2 we will provide a brief survey of implicit numerical methods for MHD. Following this, we will focus on two broad classes of nonlinearly implicit methods: Newton-Krylov (Section 3) and FAS or nonlinear multigrid (Section 4). Instead of writing a survey of a large number of implicit methods, we will present details of implicit methods with explanations of two different Newton-Krylov approaches (differing in the preconditioning approach) and on one nonlinear multigrid implementation for MHD. The chapter will close with a section on simulation test cases, and finally a conclusion section.

#### **1.1.1 Rationale for implicit treatment**

In compressible MHD we encounter the fast magnetosonic, Alfvén, and the slow magnetosonic waves. Typically, plasma confinement devices, such as tokamaks, stellarators, reversed field pinches etc. are characterized by a long scale length in the direction of the magnetic field, and shorter length scale phenomena in the direction perpendicular to the field. For example, in a tokamak, the magnetic field is dominantly along the toroidal direction and consequently the long length scale is mostly along the toroidal direction whereas the short scales are in the radial-poloidal plane. It is known that the Alfvén wave is a transverse wave with fastest propagation along the magnetic field. The fast magnetosonic, i.e., the fast compressional wave, is also anisotropic with the fastest propagation perpendicular to the magnetic field. In explicit methods, the time step is restricted by the familiar CFL condition. Several MHD phenomena are studied for long-time behavior where long-time is of the order of resistive or a combination of resistive-Alfvén time scales. For such investigations, the CFL condition implies an overly restrictive time step which translates to an enormous number of time steps. It is advantageous and desirable to design numerical schemes which allow us to have time steps larger than that imposed by the CFL condition, and yet the computational cost of each time step is only slightly larger than the explicit case. Finally, we remark that as we progress towards petascale computing and beyond to exascale, it is well recognized that breakthroughs and discoveries in science will be well-enabled by massive computations. However, hardware capability alone will not be sufficient and must be complemented by a large increase in efficiency by the development of clever algorithms. Implicit methods may prove beneficial as simulations increase in scale, since explicit methods can succumb to poor parallel weak scaling (Keyes et al. (2006)).

#### **1.2 Resistive MHD equations**

The single-fluid resistive MHD equations couple the equations of hydrodynamics and resistive low-frequency Maxwell's equations, and may be written in conservation form as

$$\frac{\partial \mathcal{U}}{\partial t} + \underbrace{\nabla \cdot F(\mathcal{U})}\_{\text{hyperbolic terms}} = \underbrace{\nabla \cdot F\_d(\mathcal{U})}\_{\text{diffusive terms}},$$
 
$$\frac{\partial \mathcal{U}}{\partial t} + R(\mathcal{U}) = 0 \tag{1}$$

2 Will-be-set-by-IN-TECH

the wave speeds are also much larger than the fluid advective speeds. In such cases, implicit time stepping methods are preferred to overcome the stiffness induced by the fast waves and march the flow simulations forward in time at time steps dictated by accuracy rather than

A brief outline of this chapter is as follows. In Section 2 we will provide a brief survey of implicit numerical methods for MHD. Following this, we will focus on two broad classes of nonlinearly implicit methods: Newton-Krylov (Section 3) and FAS or nonlinear multigrid (Section 4). Instead of writing a survey of a large number of implicit methods, we will present details of implicit methods with explanations of two different Newton-Krylov approaches (differing in the preconditioning approach) and on one nonlinear multigrid implementation for MHD. The chapter will close with a section on simulation test cases, and finally a

In compressible MHD we encounter the fast magnetosonic, Alfvén, and the slow magnetosonic waves. Typically, plasma confinement devices, such as tokamaks, stellarators, reversed field pinches etc. are characterized by a long scale length in the direction of the magnetic field, and shorter length scale phenomena in the direction perpendicular to the field. For example, in a tokamak, the magnetic field is dominantly along the toroidal direction and consequently the long length scale is mostly along the toroidal direction whereas the short scales are in the radial-poloidal plane. It is known that the Alfvén wave is a transverse wave with fastest propagation along the magnetic field. The fast magnetosonic, i.e., the fast compressional wave, is also anisotropic with the fastest propagation perpendicular to the magnetic field. In explicit methods, the time step is restricted by the familiar CFL condition. Several MHD phenomena are studied for long-time behavior where long-time is of the order of resistive or a combination of resistive-Alfvén time scales. For such investigations, the CFL condition implies an overly restrictive time step which translates to an enormous number of time steps. It is advantageous and desirable to design numerical schemes which allow us to have time steps larger than that imposed by the CFL condition, and yet the computational cost of each time step is only slightly larger than the explicit case. Finally, we remark that as we progress towards petascale computing and beyond to exascale, it is well recognized that breakthroughs and discoveries in science will be well-enabled by massive computations. However, hardware capability alone will not be sufficient and must be complemented by a large increase in efficiency by the development of clever algorithms. Implicit methods may prove beneficial as simulations increase in scale, since explicit methods can succumb to poor

The single-fluid resistive MHD equations couple the equations of hydrodynamics and resistive low-frequency Maxwell's equations, and may be written in conservation form as

> <sup>=</sup> ∇ · *Fd*(*U*) *di f f usive terms*

,

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>R</sup>*(*U*) = <sup>0</sup> (1)

stability constraints.

conclusion section.

**1.1.1 Rationale for implicit treatment**

parallel weak scaling (Keyes et al. (2006)).

*∂U*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ · *<sup>F</sup>*(*U*) *hyperbolic terms*

*∂U*

**1.2 Resistive MHD equations**

where *R*(*U*) ≡∇· *F*(*U*) −∇· *Fd*(*U*), and the solution vector *U* ≡ *U*(**x**, *t*) is

$$\boldsymbol{U} = \left\{ \rho, \rho \mathbf{u}, \mathbf{B}, \boldsymbol{e} \right\}^T$$

and the hyperbolic flux vector *F*(*U*) and the diffusive fluxes *Fd*(*U*) are given by

$$\begin{split} F(\mathcal{U}) &= \left\{ \rho u \,, \rho u \boldsymbol{u} + \left( p + \frac{1}{2} \mathbf{B} \cdot \mathbf{B} \right) \bar{\mathbf{I}} - \mathbf{B} \boldsymbol{B} \right. \\ & \left. u \boldsymbol{B} - \mathbf{B} u \,, \left( \boldsymbol{\varepsilon} + p + \frac{1}{2} \mathbf{B} \cdot \mathbf{B} \right) \boldsymbol{u} - \mathbf{B} (\mathbf{B} \cdot \boldsymbol{u}) \right\}^{\mathrm{T}}, \\ F\_{d}(\mathcal{U}) &= \left\{ 0 \,, \operatorname{Re}^{-1} \bar{\mathbf{o}} \; \operatorname{S}^{-1} \left( \eta \nabla \mathbf{B} - \eta (\nabla \mathbf{B})^{\mathrm{T}} \right) \right\}, \\ & \qquad \operatorname{Re}^{-1} \bar{\mathbf{o}} \cdot \boldsymbol{u} + \frac{\gamma}{\gamma - 1} \frac{\kappa}{\operatorname{Re} \operatorname{Pr}} \nabla T + \frac{\eta}{\mathcal{S}} \left( \frac{1}{2} \nabla (\mathbf{B} \cdot \mathbf{B}) - \mathbf{B} (\nabla \mathbf{B})^{\mathrm{T}} \right) \right\}^{\mathrm{T}}. \end{split} \tag{2}$$

In the above equations *ρ* is the density, *u* is the velocity, *B* is the magnetic field, *p* and *T* are the pressure and temperature respectively, and *e* is the total energy per unit volume of the plasma. The plasma properties are the resistivity *η*, the thermal conductivity *κ*, and the viscosity *μ*, which have been normalized, respectively, by a reference resistivity *ηR*, a reference conductivity *κR*, and a reference viscosity *μR*. The ratio of specific heats is denoted by *γ* and generally fixed at 5/3 in most MHD simulations. The non-dimensional parameters in the above equations are the Reynolds number, defined as *Re* ≡ *ρ*0*U*0*L*/*μR*, the Lundquist number, defined as *S* ≡ *μ*0*U*0*L*/*ηR*, and the Prandtl number, denoted by *Pr*, which is the ratio of momentum to thermal diffusivity. The non-dimensionalization is carried out using a characteristic length scale, *<sup>L</sup>*, and the Alfvén speed *<sup>U</sup>*<sup>0</sup> = *<sup>B</sup>*0/√*μ*0*ρ*0, where *<sup>B</sup>*0, *<sup>ρ</sup>*0, and *<sup>μ</sup>*<sup>0</sup> are the characteristic strength of the magnetic field, a reference density, and the permeability of free space, respectively. The equations are closed by the following equation of state

$$e = \frac{p}{\gamma - 1} + \frac{\rho}{2}\mu \cdot \mu + \frac{1}{2}B \cdot B\_A$$

and the stress-strain tensor relation

$$\bar{\mathfrak{o}} = \mu \left( \nabla \mathfrak{u} + \left( \nabla \mathfrak{u} \right)^{T} \right) - \frac{2}{3} \mu \nabla \cdot \mathfrak{u} \bar{\mathbf{I}}.$$

Finally, a consequence of Faraday's law is that an initially divergence-free magnetic field must lead to a divergence-free magnetic field for all times, which corresponds to the lack of observations of magnetic monopoles in nature. This solenoidal property is expressed as ∇ · *B* = 0.

In the limit of zero resistivity, conductivity and viscosity, the equations of resistive MHD reduced to those of ideal MHD. These equations are similar to those written above with *κ* = *μ* = *η* = 0. Ideal MHD equations are hyperbolic PDEs (although not strictly hyperbolic).

#### **2. Brief survey of implicit methods for MHD**

#### **2.1 Early approaches**

An implicit treatment of the fast magnetosonic wave coupled with arguments of large length scales dominantly in a certain direction allows one to investigate long-time scale phenomena

The above equation is implicit and is solved iteratively. Let *Un*+1,*<sup>k</sup>* denote the k-th iteration of

Implicit Numerical Methods for Magnetohydrodynamics 63

(3*Un*+1,*k*+<sup>1</sup>

*k i*,*j* + *∂ <sup>∂</sup><sup>U</sup> ∂t k i*,*j ∂U*

*∂R<sup>k</sup> i*,*j <sup>∂</sup><sup>U</sup>* <sup>Δ</sup>*U<sup>k</sup> i*,*j*

*<sup>i</sup>*,*<sup>j</sup>* +

fluxes with respect to the solution vector is difficult to evaluate for second order upwind schemes. Hence, at this stage an approximation is made, i.e., such terms are evaluated with a first order scheme. The first order hyperbolic flux divergence, denoted as *R*¯ is at point (*i*, *j*)

system above is a large banded matrix and will be generally expensive to invert. Instead Jones et al. recommend the use of further approximations and using a lower-upper Gauss-Seidel (LU-SGS) technique. The hyperbolic fluxes *R* are evaluated with the Harten's approximate Riemann solver (Harten (1983)), applied with the framework of the eight-wave scheme developed by Powell et al. (1999). One philosophical concern about implicit upwinding methods is as follows. Generally, upwind methods are based on the solution of a Riemann problem at cell faces; such a solution is self-similar in time, i.e. depends only on *x*/*t* for times until the waves from neighboring cell faces Riemann problems start interacting. In traditional explicit upwind methods, this problem is avoided because we are operating within the CFL limit. In an implicit method, the CFL limit is violated and waves from neighboring Riemann solvers will interact. One may adopt the viewpoint that upwind methods are, in a sense, providing dissipation proportional to each wave and decrease the dispersion error which are the bane of central difference schemes. Adopting this viewpoint, one may ignore the

 *Rk <sup>i</sup>*,*<sup>j</sup>* +

<sup>=</sup> <sup>−</sup>*Rn*+1,*k*+<sup>1</sup>

*<sup>i</sup>*,*<sup>j</sup>* <sup>−</sup> <sup>4</sup>*U<sup>n</sup>*

*<sup>i</sup>*,*<sup>j</sup>* <sup>+</sup> *<sup>U</sup>n*−<sup>1</sup>

Δ*U<sup>k</sup>*

*<sup>i</sup>*,*<sup>j</sup>* . The partial derivative of the divergence of the hyperbolic

 *∂U ∂t*

*k i*,*j* *<sup>i</sup>*,*<sup>j</sup>* , (8)

*<sup>i</sup>*,*<sup>j</sup>* ). (9)

*<sup>i</sup>*,*<sup>j</sup>* (10)

(12)

, (11)

2 ,*j* , *Ui*<sup>−</sup> <sup>1</sup> 2 ,*j* , *Ui*,*j*+<sup>1</sup> 2 , *Ui*,*j*<sup>−</sup> <sup>1</sup> 2 ).

*ij* is driven to zero. The matrix in the linear

the solution at the *n* + 1-th time level. Rewrite the above equation as

 *∂U ∂t*

A truncated Taylor series expansion yields:

*<sup>i</sup>*,*<sup>j</sup>* <sup>=</sup> *<sup>U</sup>n*+1,*k*+<sup>1</sup>

*<sup>k</sup>*+<sup>1</sup> *i*,*j*

 *∂U ∂t*

*<sup>i</sup>*,*<sup>j</sup>* <sup>−</sup> *<sup>U</sup>n*+1,*<sup>k</sup>*

 *∂R*ˆ *<sup>k</sup> i*,*j <sup>∂</sup><sup>U</sup>* <sup>+</sup>

The above equation is linear and iterated until Δ*U<sup>k</sup>*

interactions between neighboring Riemann problems.

Substituting all back into equation (8) gives

*<sup>k</sup>*+<sup>1</sup> *i*,*j*

*Rk*+<sup>1</sup> *<sup>i</sup>*,*<sup>j</sup>* <sup>=</sup> *<sup>R</sup><sup>k</sup>*

is coupled to the neighboring four points in 2D. *<sup>R</sup>*ˆ*i*,*<sup>j</sup>* <sup>≡</sup> *<sup>R</sup>*ˆ(*Ui*,*j*, *Ui*<sup>+</sup> <sup>1</sup>

3*I* 2Δ*t*  Δ*U<sup>k</sup> <sup>i</sup>*,*<sup>j</sup>* = −

= *∂U ∂t*

where

where Δ*U<sup>k</sup>*

 *∂U ∂t*

> <sup>≡</sup> <sup>1</sup> 2Δ*t*

*<sup>k</sup>*+<sup>1</sup> *i*,*j*

in MHD in a computationally efficient manner. The approach discussed here was developed by Harned & Kerner (1985). We begin our discussion with a model problem which exposes the philosophy behind the implicit treatment of the fastest waves in MHD. Consider the following hyperbolic system of equations.

$$\begin{aligned} \frac{\partial u}{\partial t} &= -a \frac{\partial v}{\partial x}, \\ \frac{\partial v}{\partial t} &= -a \frac{\partial u}{\partial x}. \end{aligned} \tag{3}$$

This can be rewritten as

$$\frac{\partial^2 u}{\partial t^2} = a^2 \frac{\partial^2 u}{\partial x^2}. \tag{4}$$

We then subtract a term from either side of the above equation as

$$a\frac{\partial^2 u}{\partial t^2} - a\_0^2 \frac{\partial^2 u}{\partial x^2} = a^2 \frac{\partial^2 u}{\partial x^2} - a\_0^2 \frac{\partial^2 u}{\partial x^2},\tag{5}$$

where *a*<sup>0</sup> is a constant coefficient chosen mainly from stability considerations. Furthermore, *a*<sup>0</sup> is something which mimics the behavior of *a*, perhaps in some limit. The underlying idea of the semi-implicit methods discussed here is this: the term containing *a*<sup>0</sup> on the left hand side of eqn. (4) is treated implicitly, while the same term on the right hand side of eqn. (4) is treated explicitly. Moreover, the cost of solving the linear system stemming from the implicit treatment of the term containing *a*<sup>0</sup> should be small relative to the total cost of evolving the entire system. Harned and Kerner generalized the above approach to the MHD system in a slab geometry, with implicit treatment of the fast compressional wave. Furthermore, their method was applicable to a case where the scale lengths in the z-direction are much longer than those in the x-y plane. The fastest time scale is then due to the fast compressional wave in the x-y plane. The method for the implicit treatment of the shear Alfvén wave was proposed by Harned & Schnack (1986). The procedure is somewhat similar to the one adopted for the fast compressional wave except the linear term which is subtracted on the velocity evolution equation has a different form. The implicit treatment of the shear Alfvén wave is, in general, more problematic and required certain ad-hoc heuristics to be employed for stability (See Reference Harned & Schnack (1986) for details).

The above approaches may be classified as linearly implicit. An example of a nonlinearly implicit method is the work of Jones, Shumlak & Eberhardt (1997) on an upwind implicit method for resistive MHD. Their method applied to ideal MHD equations may be written as:

$$\frac{\partial \mathcal{U}}{\partial t} = -R(\mathcal{U}) = -\left(\frac{\partial F}{\partial x} + \frac{\partial G}{\partial y}\right) \tag{6}$$

where *R*(*U*) is the divergence of the hyperbolic fluxes (here in this 2D discussion, *F* and *G* denote the fluxes in the *x*− and *y*− directions, respectively). The above equation (or rather system of equations) is discretized in time as

$$\frac{1}{2\Delta t}(3\mathcal{U}\_{i,j}^{n+1} - 4\mathcal{U}\_{i,j}^{n} + \mathcal{U}\_{i,j}^{n-1}) = -\mathcal{R}\_{i,j}^{n+1} \tag{7}$$

The above equation is implicit and is solved iteratively. Let *Un*+1,*<sup>k</sup>* denote the k-th iteration of the solution at the *n* + 1-th time level. Rewrite the above equation as

$$\left(\frac{\partial \mathcal{U}}{\partial t}\right)\_{i,j}^{k+1} = -\mathcal{R}\_{i,j}^{n+1,k+1}\,'\,\tag{8}$$

where

(3)

(6)

*<sup>i</sup>*,*<sup>j</sup>* (7)

4 Will-be-set-by-IN-TECH

in MHD in a computationally efficient manner. The approach discussed here was developed by Harned & Kerner (1985). We begin our discussion with a model problem which exposes the philosophy behind the implicit treatment of the fastest waves in MHD. Consider the

> *∂v ∂x* ,

*∂u ∂x* .

*<sup>∂</sup>x*<sup>2</sup> <sup>−</sup> *<sup>a</sup>*<sup>2</sup> 0 *∂*2*u*

*<sup>∂</sup>x*<sup>2</sup> . (4)

*<sup>∂</sup>x*<sup>2</sup> , (5)

*∂u <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

*∂v <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

*∂*2*u*

We then subtract a term from either side of the above equation as

*∂*2*u <sup>∂</sup>t*<sup>2</sup> <sup>−</sup> *<sup>a</sup>*<sup>2</sup> 0 *∂*2*u*

Reference Harned & Schnack (1986) for details).

system of equations) is discretized in time as

*∂U*

1 2Δ*t*

(3*Un*+<sup>1</sup>

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*R*(*U*) = <sup>−</sup>

*<sup>i</sup>*,*<sup>j</sup>* <sup>−</sup> <sup>4</sup>*U<sup>n</sup>*

as:

*<sup>∂</sup>t*<sup>2</sup> <sup>=</sup> *<sup>a</sup>*<sup>2</sup> *<sup>∂</sup>*2*<sup>u</sup>*

*<sup>∂</sup>x*<sup>2</sup> <sup>=</sup> *<sup>a</sup>*<sup>2</sup> *<sup>∂</sup>*2*<sup>u</sup>*

where *a*<sup>0</sup> is a constant coefficient chosen mainly from stability considerations. Furthermore, *a*<sup>0</sup> is something which mimics the behavior of *a*, perhaps in some limit. The underlying idea of the semi-implicit methods discussed here is this: the term containing *a*<sup>0</sup> on the left hand side of eqn. (4) is treated implicitly, while the same term on the right hand side of eqn. (4) is treated explicitly. Moreover, the cost of solving the linear system stemming from the implicit treatment of the term containing *a*<sup>0</sup> should be small relative to the total cost of evolving the entire system. Harned and Kerner generalized the above approach to the MHD system in a slab geometry, with implicit treatment of the fast compressional wave. Furthermore, their method was applicable to a case where the scale lengths in the z-direction are much longer than those in the x-y plane. The fastest time scale is then due to the fast compressional wave in the x-y plane. The method for the implicit treatment of the shear Alfvén wave was proposed by Harned & Schnack (1986). The procedure is somewhat similar to the one adopted for the fast compressional wave except the linear term which is subtracted on the velocity evolution equation has a different form. The implicit treatment of the shear Alfvén wave is, in general, more problematic and required certain ad-hoc heuristics to be employed for stability (See

The above approaches may be classified as linearly implicit. An example of a nonlinearly implicit method is the work of Jones, Shumlak & Eberhardt (1997) on an upwind implicit method for resistive MHD. Their method applied to ideal MHD equations may be written

where *R*(*U*) is the divergence of the hyperbolic fluxes (here in this 2D discussion, *F* and *G* denote the fluxes in the *x*− and *y*− directions, respectively). The above equation (or rather

*<sup>i</sup>*,*<sup>j</sup>* <sup>+</sup> *<sup>U</sup>n*−<sup>1</sup>

*<sup>i</sup>*,*<sup>j</sup>* ) = <sup>−</sup>*Rn*+<sup>1</sup>

 *∂F ∂x* + *∂G ∂y* 

following hyperbolic system of equations.

This can be rewritten as

$$\left(\frac{\partial \mathcal{U}}{\partial t}\right)\_{i,j}^{k+1} \equiv \frac{1}{2\Delta t} (3\mathcal{U}\_{i,j}^{n+1,k+1} - 4\mathcal{U}\_{i,j}^{n} + \mathcal{U}\_{i,j}^{n-1}).\tag{9}$$

A truncated Taylor series expansion yields:

$$
\left(\frac{\partial \mathcal{U}}{\partial t}\right)\_{i,j}^{k+1} = \left(\frac{\partial \mathcal{U}}{\partial t}\right)\_{i,j}^k + \frac{\partial \left(\frac{\partial \mathcal{U}}{\partial t}\right)\_{i,j}^k}{\partial \mathcal{U}} \Delta \mathcal{U}\_{i,j}^k \tag{10}
$$

$$R\_{i,j}^{k+1} = R\_{i,j}^k + \frac{\partial R\_{i,j}^k}{\partial \mathcal{U}} \Delta \mathcal{U}\_{i,j\prime}^k \tag{11}$$

where Δ*U<sup>k</sup> <sup>i</sup>*,*<sup>j</sup>* <sup>=</sup> *<sup>U</sup>n*+1,*k*+<sup>1</sup> *<sup>i</sup>*,*<sup>j</sup>* <sup>−</sup> *<sup>U</sup>n*+1,*<sup>k</sup> <sup>i</sup>*,*<sup>j</sup>* . The partial derivative of the divergence of the hyperbolic fluxes with respect to the solution vector is difficult to evaluate for second order upwind schemes. Hence, at this stage an approximation is made, i.e., such terms are evaluated with a first order scheme. The first order hyperbolic flux divergence, denoted as *R*¯ is at point (*i*, *j*) is coupled to the neighboring four points in 2D. *<sup>R</sup>*ˆ*i*,*<sup>j</sup>* <sup>≡</sup> *<sup>R</sup>*ˆ(*Ui*,*j*, *Ui*<sup>+</sup> <sup>1</sup> 2 ,*j* , *Ui*<sup>−</sup> <sup>1</sup> 2 ,*j* , *Ui*,*j*+<sup>1</sup> 2 , *Ui*,*j*<sup>−</sup> <sup>1</sup> 2 ). Substituting all back into equation (8) gives

$$
\left[\frac{\partial \hat{R}^k\_{i,j}}{\partial \boldsymbol{U}} + \frac{3\boldsymbol{I}}{2\Delta t}\right] \Delta \boldsymbol{U}^k\_{i,j} = -\left[\boldsymbol{R}^k\_{i,j} + \left(\frac{\partial \boldsymbol{U}}{\partial t}\right)^k\_{i,j}\right] \tag{12}
$$

The above equation is linear and iterated until Δ*U<sup>k</sup> ij* is driven to zero. The matrix in the linear system above is a large banded matrix and will be generally expensive to invert. Instead Jones et al. recommend the use of further approximations and using a lower-upper Gauss-Seidel (LU-SGS) technique. The hyperbolic fluxes *R* are evaluated with the Harten's approximate Riemann solver (Harten (1983)), applied with the framework of the eight-wave scheme developed by Powell et al. (1999). One philosophical concern about implicit upwinding methods is as follows. Generally, upwind methods are based on the solution of a Riemann problem at cell faces; such a solution is self-similar in time, i.e. depends only on *x*/*t* for times until the waves from neighboring cell faces Riemann problems start interacting. In traditional explicit upwind methods, this problem is avoided because we are operating within the CFL limit. In an implicit method, the CFL limit is violated and waves from neighboring Riemann solvers will interact. One may adopt the viewpoint that upwind methods are, in a sense, providing dissipation proportional to each wave and decrease the dispersion error which are the bane of central difference schemes. Adopting this viewpoint, one may ignore the interactions between neighboring Riemann problems.

most commonly used Krylov methods are GMRES (Generalized Minimum Residual) and BiCGStab (Bi conjugate gradient stabilized) which can both handle non-symmetric linear systems. GMRES is very robust but generally is heavy on memory usage, while BiCGStab has a lower memory requirement, it is less robust given that the residual is not guaranteed to

Implicit Numerical Methods for Magnetohydrodynamics 65

It the approximate solution *Un*+1,*<sup>k</sup>* is "close" to the true solution *U*<sup>∗</sup> of the nonlinear system,

where *C* is a constant independent of *Un*+1,*<sup>k</sup>* and *U*∗. This result assumes that the linear system is solved exactly. If the linear systems are solved inexactly as in the Newton-Krylov method, then *Itol*, the linear system tolerance, has to be carefully chosen. In inexact NK, *Itol* = *<sup>η</sup>k*||F*k*||. Quadratic convergence is retained if *<sup>η</sup><sup>k</sup>* <sup>=</sup> *<sup>C</sup>*||F*k*||. If we impose the condition that *limk*→∞*η<sup>k</sup>* <sup>=</sup> 0 then convergence is super-linear, and if *<sup>η</sup><sup>k</sup>* is constant for all *<sup>k</sup>* then convergence is linear. Since Newton's method may be viewed as a linear model of the original nonlinear system, the model is a better approximation as the solution is approached. When "far" from the solution, it is not essential to solve the linear system to machine-zero convergence. The following choices for *η<sup>k</sup>* are recommended which take into account how well the nonlinear


2

where *γ*<sup>1</sup> = 0.9 and *γ*<sup>2</sup> = 2 as recommended by Eisenstat & Walker (1996). The first of these choices is how well the linear model agreed with the nonlinear system at the prior step, while

In examining the Krylov methods, we notice that these require only matrix-vector products. Thus it is never necessary to store the entire Jacobian matrix. Hence the term "Jacobian-Free Newton-Krylov" (abbreviated JFNK) is frequently encountered in the literature. Furthermore,


 


, (19)

1. Begin by guessing the solution *Un*<sup>+</sup>1,0. Typically the initial guess is *Un*+1,0 = *Un*.

*<sup>J</sup>*(*U<sup>k</sup>* )*δU<sup>k</sup>* <sup>=</sup> −F(*Un*+1,*k*) so that ||*J*(*Uk*)*δU<sup>k</sup>* <sup>+</sup> <sup>F</sup>(*Un*+1,*k*)|| <sup>&</sup>lt; *Itol*.

decrease monotonically.

Steps in a typical NK solver are the following:

(a) Using a Krylov method, approximately solve for *δUk*,

3. Update the Newton iterate, *Un*+1,*k*+<sup>1</sup> = *Un*+1,*<sup>k</sup>* + *λδU<sup>k</sup>*

*η<sup>k</sup>* =

*η<sup>k</sup>* = *γ*<sup>1</sup>

 

the second uses a measure of the rate of convergence of the linear system.


2. For each Newton iteration *k* = 1, 2, ···

Each Krylov iteration requires: a. One matrix-vector multiply with J

4. Test for convergence ||F(*Un*+1,*k*+1)|| <sup>&</sup>lt; *f tol*.

b. One preconditioner solve

the convergence is quadratic, i.e.,

system is converging.

#### **2.2 Modern approaches**

We, somewhat arbitrarily, classify *fully implicit* numerical methods (in distinction from semi-implicit or linearly implicit) as "modern". The main feature which distinguishes these approaches from semi-implicit or linear implicit is the ability to allow for very large time steps. The modern era of fully nonlinearly implicitly solvers was ushered in the since the early-mid nineties. Broadly the fully implicit methods can be classified as: (a) Newton-Krylov and (b) nonlinear multigrid (also known as FAS, i.e., full approximation scheme). An early example of an implicit Newton-Krylov-Schwarz method applied to aerodynamics was by Keyes (1995). Several papers subsequently appeared in the mid-late nineties and in early part of this century in fluid dynamics. Newton-Krylov methods found applicability in MHD in the early 2000s. In the subsequent sections, we will elaborate on both the Newton-Krylov and nonlinear multigrid as applied to MHD.

### **3. Newton-Krylov (NK) methods for MHD**

#### **3.1 General approach**

The entire ideal MHD (or resistive MHD and beyond) can be written as a nonlinear function as follows:

$$\mathcal{F}(\mathcal{U}^{n+1}) = 0,\tag{13}$$

where *Un*+<sup>1</sup> is the vector of unknowns at time step *n* + 1. For example, if we use a *θ*-scheme, one can write the nonlinear function as:

$$\mathcal{F}(\mathcal{U}^{n+1}) = \mathcal{U}^{n+1} - \mathcal{U}^n + \theta \Delta t R(\mathcal{U}^{n+1}) + (1 - \theta)\Delta t R(\mathcal{U}^n) = 0,\tag{14}$$

where *R*(*U*) is the divergence of the fluxes (see equation 1). For compressible MHD, on a two dimensional *N* × *M* mesh, the total number of unknowns would then be 8*MN*. The above nonlinear systems can be solved using an inexact Newton–Krylov solver. Apply the standard Newton's method to the above nonlinear system gives

$$
\delta \mathcal{U}^k = -\left[ \left( \frac{\partial \mathcal{F}}{\partial \mathcal{U}} \right)^{n+1,k} \right]^{-1} \mathcal{F}\_\prime \tag{15}
$$

where *<sup>J</sup>*(*Un*+1,*k*) <sup>≡</sup> *<sup>∂</sup>*<sup>F</sup> *∂U n*+1,*<sup>k</sup>* is the Jacobian; and *<sup>δ</sup>U<sup>k</sup>* <sup>≡</sup> *<sup>U</sup>n*+1,*k*+<sup>1</sup> <sup>−</sup> *<sup>U</sup>n*+1,*k*, and *<sup>k</sup>* is the iteration index in the Newton method. For the two dimensional system the Jacobian matrix is 8*MN* × 8*MN* which, although sparse, is still impractical to invert directly.

In NK methods, the linear system at each Newton step is solved by a Krylov method. In Krylov methods, an approximation to the solution of the linear system *JδU* = −F is obtained by iteratively building a Krylov subspace of dimension *m* defined by

$$\mathcal{K}(r\_{0\prime}f) = \text{span}\{r\_{0\prime}f r\_{0\prime}f^2 r\_{0\prime} \cdots \prime, f^{m-1} r\_0\},\tag{16}$$

where *r*<sup>0</sup> is the initial residual of the linear system. The Krylov method can be either: one in which the solution in the subspace minimizes the linear system residual, or two in which the residual is orthogonal to the Krylov subspace. Within Newton-Krylov methods the two 6 Will-be-set-by-IN-TECH

We, somewhat arbitrarily, classify *fully implicit* numerical methods (in distinction from semi-implicit or linearly implicit) as "modern". The main feature which distinguishes these approaches from semi-implicit or linear implicit is the ability to allow for very large time steps. The modern era of fully nonlinearly implicitly solvers was ushered in the since the early-mid nineties. Broadly the fully implicit methods can be classified as: (a) Newton-Krylov and (b) nonlinear multigrid (also known as FAS, i.e., full approximation scheme). An early example of an implicit Newton-Krylov-Schwarz method applied to aerodynamics was by Keyes (1995). Several papers subsequently appeared in the mid-late nineties and in early part of this century in fluid dynamics. Newton-Krylov methods found applicability in MHD in the early 2000s. In the subsequent sections, we will elaborate on both the Newton-Krylov

The entire ideal MHD (or resistive MHD and beyond) can be written as a nonlinear function

where *Un*+<sup>1</sup> is the vector of unknowns at time step *n* + 1. For example, if we use a *θ*-scheme,

where *R*(*U*) is the divergence of the fluxes (see equation 1). For compressible MHD, on a two dimensional *N* × *M* mesh, the total number of unknowns would then be 8*MN*. The above nonlinear systems can be solved using an inexact Newton–Krylov solver. Apply the standard

*<sup>n</sup>*+1,*<sup>k</sup>*

 *<sup>∂</sup>*<sup>F</sup> *∂U*

iteration index in the Newton method. For the two dimensional system the Jacobian matrix is

In NK methods, the linear system at each Newton step is solved by a Krylov method. In Krylov methods, an approximation to the solution of the linear system *JδU* = −F is obtained

where *r*<sup>0</sup> is the initial residual of the linear system. The Krylov method can be either: one in which the solution in the subspace minimizes the linear system residual, or two in which the residual is orthogonal to the Krylov subspace. Within Newton-Krylov methods the two

<sup>F</sup>(*Un*+1) = *<sup>U</sup>n*+<sup>1</sup> <sup>−</sup> *<sup>U</sup><sup>n</sup>* <sup>+</sup> *<sup>θ</sup>*Δ*tR*(*Un*+1)+(<sup>1</sup> <sup>−</sup> *<sup>θ</sup>*)Δ*tR*(*Un*) = 0, (14)

−<sup>1</sup>

<sup>2</sup>*r*0, ··· , *<sup>J</sup>*

is the Jacobian; and *<sup>δ</sup>U<sup>k</sup>* <sup>≡</sup> *<sup>U</sup>n*+1,*k*+<sup>1</sup> <sup>−</sup> *<sup>U</sup>n*+1,*k*, and *<sup>k</sup>* is the

<sup>F</sup>(*Un*+1) = 0, (13)

F, (15)

*<sup>m</sup>*−1*r*0}, (16)

**2.2 Modern approaches**

**3.1 General approach**

where *<sup>J</sup>*(*Un*+1,*k*) <sup>≡</sup>

as follows:

and nonlinear multigrid as applied to MHD.

**3. Newton-Krylov (NK) methods for MHD**

one can write the nonlinear function as:

Newton's method to the above nonlinear system gives

*n*+1,*<sup>k</sup>*

 *<sup>∂</sup>*<sup>F</sup> *∂U* *<sup>δ</sup>U<sup>k</sup>* <sup>=</sup> <sup>−</sup>

8*MN* × 8*MN* which, although sparse, is still impractical to invert directly.

K(*r*0, *J*) = *span*{*r*0, *Jr*0, *J*

by iteratively building a Krylov subspace of dimension *m* defined by

most commonly used Krylov methods are GMRES (Generalized Minimum Residual) and BiCGStab (Bi conjugate gradient stabilized) which can both handle non-symmetric linear systems. GMRES is very robust but generally is heavy on memory usage, while BiCGStab has a lower memory requirement, it is less robust given that the residual is not guaranteed to decrease monotonically.

Steps in a typical NK solver are the following:

	- (a) Using a Krylov method, approximately solve for *δUk*, *<sup>J</sup>*(*U<sup>k</sup>* )*δU<sup>k</sup>* <sup>=</sup> −F(*Un*+1,*k*) so that ||*J*(*Uk*)*δU<sup>k</sup>* <sup>+</sup> <sup>F</sup>(*Un*+1,*k*)|| <sup>&</sup>lt; *Itol*. Each Krylov iteration requires: a. One matrix-vector multiply with J b. One preconditioner solve

It the approximate solution *Un*+1,*<sup>k</sup>* is "close" to the true solution *U*<sup>∗</sup> of the nonlinear system, the convergence is quadratic, i.e.,

$$||\mathcal{U}^{n+1,k+1} - \mathcal{U}^\*|| \le \mathbb{C}||\mathcal{U}^{n+1,k} - \mathcal{U}^\*||^2,\tag{17}$$

where *C* is a constant independent of *Un*+1,*<sup>k</sup>* and *U*∗. This result assumes that the linear system is solved exactly. If the linear systems are solved inexactly as in the Newton-Krylov method, then *Itol*, the linear system tolerance, has to be carefully chosen. In inexact NK, *Itol* = *<sup>η</sup>k*||F*k*||. Quadratic convergence is retained if *<sup>η</sup><sup>k</sup>* <sup>=</sup> *<sup>C</sup>*||F*k*||. If we impose the condition that *limk*→∞*η<sup>k</sup>* <sup>=</sup> 0 then convergence is super-linear, and if *<sup>η</sup><sup>k</sup>* is constant for all *<sup>k</sup>* then convergence is linear. Since Newton's method may be viewed as a linear model of the original nonlinear system, the model is a better approximation as the solution is approached. When "far" from the solution, it is not essential to solve the linear system to machine-zero convergence. The following choices for *η<sup>k</sup>* are recommended which take into account how well the nonlinear system is converging.

$$\eta^k = \frac{\left| ||\mathcal{F}^k - ||{f^{k-1} \delta \mathcal{U}^{k-1} + \mathcal{F}^{k-1}}|| \right|}{||{\mathcal{F}^{k-1} ||} }\tag{18}$$

$$\eta^k = \gamma\_1 \left( \frac{||\mathcal{F}^k||}{||\mathcal{F}^{k-1}||} \right)\_2^\gamma \tag{19}$$

where *γ*<sup>1</sup> = 0.9 and *γ*<sup>2</sup> = 2 as recommended by Eisenstat & Walker (1996). The first of these choices is how well the linear model agreed with the nonlinear system at the prior step, while the second uses a measure of the rate of convergence of the linear system.

In examining the Krylov methods, we notice that these require only matrix-vector products. Thus it is never necessary to store the entire Jacobian matrix. Hence the term "Jacobian-Free Newton-Krylov" (abbreviated JFNK) is frequently encountered in the literature. Furthermore,

We note that the approximations used in the preconditioner should have no effect on the overall accuracy of the nonlinear system. It can be shown that JFNK method applied with right preconditioning preserves the conservation properties of the equations written in conservation form (Chacón (2004)) regardless of the nonlinear convergence tolerances. However, one cannot prove this for left preconditioning unless the solution is converged to

Implicit Numerical Methods for Magnetohydrodynamics 67

• Algebraic preconditioners: The nature of such preconditioners is of the "black-box" type. These represent a close representation of the Jacobian and are obtained using relatively inexpensive algebraic techniques such as stationary iterative techniques, incomplete LU decomposition, multigrid techniques etc. These preconditioners typically require forming

• "Physics-based" preconditioners: These preconditioners are derived from other techniques such Picard iteration, or by semi-implicit techniques. They do not require forming and storing the entire Jacobian matrix and can be harnessed for Jacobian-free implementations. The form of the preconditioners here generally tend to exploit the structure of the PDEs

In this section, we essentially reproduce the work by Chacón (2008a), wherein a JFNK approach for resistive MHD with physics based preconditioners has been developed. The approach, given below, essentially relies on the trick of "parabolization" and using a Schur complement approach. Parabolization refers to the technique by which a hyperbolic system

themselves and in this sense this type of preconditioning is "physics-based".

is converted to a parabolic one which is then amenable to multigrid techniques.

*∂u <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

*∂v <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

*un*+<sup>1</sup> = *u<sup>n</sup>* + *a*

*vn*+<sup>1</sup> = *v<sup>n</sup>* + *a*

<sup>2</sup> *∂*<sup>2</sup> *∂x*<sup>2</sup>  *∂v ∂x* ,

*∂u ∂x*

 *∂v ∂x*

 *∂u ∂x*

*un*+<sup>1</sup> = *u<sup>n</sup>* + *a*Δ*t*

*<sup>n</sup>*+<sup>1</sup> ,

*<sup>n</sup>*+<sup>1</sup>

 *∂v ∂x* *<sup>n</sup>*

. (22)

. (23)

, (24)

machine precision (Chacón (2004)).

and storing the Jacobian matrix.

**3.3 JFNK method for resistive MHD I**

Consider the following hyperbolic system

Differencing with backward Euler we get

Substitute the second equation into the first to obtain:

*<sup>I</sup>* <sup>−</sup> *<sup>a</sup>*2Δ*<sup>t</sup>*

**3.3.1 A model illustration**

Preconditioners can be divided into two broad classes:

for complicated nonlinear systems such as those arising in MHD, the Jacobian entries are not even known analytically. Instead one can conveniently evaluate the Jacobian vector product using first order finite differences as follows:

$$J(\mathcal{U}^k)\delta \mathcal{U}^k \approx \frac{\mathcal{F}(\mathcal{U}^{n+1,k} + \sigma \delta \mathcal{U}^k) - \mathcal{F}(\mathcal{U}^{n+1,k})}{\sigma},\tag{20}$$

where *σ* is typically used as the square root of machine zero. The above expression assumes that F is sufficiently differentiable, a property which is easily violated in upwind methods with its myriad switches, and limited-reconstruction methods. The beauty of the Newton-Krylov method as outlined above is that it only relies on the evaluation of the nonlinear function F. For a detailed review of the field of JFNK see the review paper by Knoll & Keyes (2004).

#### **3.2 Preconditioners**

Since all operations in the Newton-Krylov context require only linear complexity operations, the key component required for scalability of fully implicit simulations using this technology is an optimal preconditioning strategy for the inner Krylov linear solver (Kelley (1995); Knoll & Keyes (2004)). In Newton-Krylov algorithms, at each Newton iteration a Krylov iterative method is used to solve Jacobian systems of the form

$$J(\mathcal{U})V = -\mathcal{F}(\mathcal{U}), \qquad J(\mathcal{U}) \equiv I + \gamma \frac{\partial}{\partial \mathcal{U}} (\mathcal{R}(\mathcal{U})), \qquad \gamma = \theta \Delta t. \tag{21}$$

The number of iterations required for convergence of a Krylov method depends on the eigenstructure of *J*, where systems with clustered eigenvalues typically result in faster convergence than those with evenly distributed eigenvalues (Greenbaum (1997); Greenbaum et al. (1996); Trefethen & Bau (1997)). Unfortunately, for a fixed Δ*t*, as the spatial resolution is refined the distribution of these eigenvalues spreads, resulting in increased numbers of Krylov iterations and hence non-scalability of the overall solution algorithm. The role of a preconditioning operator *P* is to transform the original Jacobian system (21) to either

$$J P^{-1} P V = -f \text{ (right prec.)}, \quad \text{or} \quad \quad P^{-1} J V = -P^{-1} f \text{ (left prec.)}.$$

The Krylov iteration is then used to solve one of

$$(JP^{-1})W = -f,\qquad \text{or} \qquad (P^{-1}f)V = X\_{\prime\prime}$$

where *<sup>X</sup>* <sup>=</sup> <sup>−</sup>*P*−<sup>1</sup> *<sup>f</sup>* is computed prior to the Krylov solve or *<sup>V</sup>* <sup>=</sup> *<sup>P</sup>*−1*<sup>W</sup>* is computed after the Krylov solve. Scalable convergence of the method then depends on the spectrum of the preconditioned operator (*JP*−<sup>1</sup> or *P*−<sup>1</sup> *J*), as opposed to the original Jacobian operator *J*. Hence, an optimal preconditioning strategy will satisfy the two competing criteria:


8 Will-be-set-by-IN-TECH

for complicated nonlinear systems such as those arising in MHD, the Jacobian entries are not even known analytically. Instead one can conveniently evaluate the Jacobian vector product

*<sup>J</sup>*(*Uk*)*δU<sup>k</sup>* <sup>≈</sup> <sup>F</sup>(*Un*+1,*<sup>k</sup>* <sup>+</sup> *σδUk*) − F(*Un*+1,*k*)

where *σ* is typically used as the square root of machine zero. The above expression assumes that F is sufficiently differentiable, a property which is easily violated in upwind methods with its myriad switches, and limited-reconstruction methods. The beauty of the Newton-Krylov method as outlined above is that it only relies on the evaluation of the nonlinear function F. For a detailed review of the field of JFNK see the review paper by

Since all operations in the Newton-Krylov context require only linear complexity operations, the key component required for scalability of fully implicit simulations using this technology is an optimal preconditioning strategy for the inner Krylov linear solver (Kelley (1995); Knoll & Keyes (2004)). In Newton-Krylov algorithms, at each Newton iteration a Krylov

The number of iterations required for convergence of a Krylov method depends on the eigenstructure of *J*, where systems with clustered eigenvalues typically result in faster convergence than those with evenly distributed eigenvalues (Greenbaum (1997); Greenbaum et al. (1996); Trefethen & Bau (1997)). Unfortunately, for a fixed Δ*t*, as the spatial resolution is refined the distribution of these eigenvalues spreads, resulting in increased numbers of Krylov iterations and hence non-scalability of the overall solution algorithm. The role of a preconditioning operator *P* is to transform the original Jacobian system (21) to either

*JP*−1*PV* <sup>=</sup> <sup>−</sup> *<sup>f</sup>* (right prec.), or *<sup>P</sup>*−<sup>1</sup> *JV* <sup>=</sup> <sup>−</sup>*P*−<sup>1</sup> *<sup>f</sup>* (left prec.).

(*JP*−1)*<sup>W</sup>* <sup>=</sup> <sup>−</sup> *<sup>f</sup>* , or (*P*−<sup>1</sup> *<sup>J</sup>*)*<sup>V</sup>* <sup>=</sup> *<sup>X</sup>*,

where *<sup>X</sup>* <sup>=</sup> <sup>−</sup>*P*−<sup>1</sup> *<sup>f</sup>* is computed prior to the Krylov solve or *<sup>V</sup>* <sup>=</sup> *<sup>P</sup>*−1*<sup>W</sup>* is computed after the Krylov solve. Scalable convergence of the method then depends on the spectrum of the preconditioned operator (*JP*−<sup>1</sup> or *P*−<sup>1</sup> *J*), as opposed to the original Jacobian operator *J*.

2. Application of *P*−<sup>1</sup> should be much more efficient than solution to the original system, optimally with linear complexity as the problem is refined and with no dependence on an

Hence, an optimal preconditioning strategy will satisfy the two competing criteria:

1. *P* ≈ *J*, to help cluster the spectrum of the preconditioned operator.

increasing number of processors in a parallel simulation.

*∂ ∂U R*(*U*) 

iterative method is used to solve Jacobian systems of the form

The Krylov iteration is then used to solve one of

*J*(*U*)*V* = −F(*U*), *J*(*U*) ≡ *I* + *γ*

*<sup>σ</sup>* , (20)

, *γ* = *θ*Δ*t*. (21)

using first order finite differences as follows:

Knoll & Keyes (2004).

**3.2 Preconditioners**

We note that the approximations used in the preconditioner should have no effect on the overall accuracy of the nonlinear system. It can be shown that JFNK method applied with right preconditioning preserves the conservation properties of the equations written in conservation form (Chacón (2004)) regardless of the nonlinear convergence tolerances. However, one cannot prove this for left preconditioning unless the solution is converged to machine precision (Chacón (2004)).

Preconditioners can be divided into two broad classes:


#### **3.3 JFNK method for resistive MHD I**

In this section, we essentially reproduce the work by Chacón (2008a), wherein a JFNK approach for resistive MHD with physics based preconditioners has been developed. The approach, given below, essentially relies on the trick of "parabolization" and using a Schur complement approach. Parabolization refers to the technique by which a hyperbolic system is converted to a parabolic one which is then amenable to multigrid techniques.

#### **3.3.1 A model illustration**

Consider the following hyperbolic system

$$\begin{aligned} \frac{\partial u}{\partial t} &= a \frac{\partial v}{\partial x}, \\ \frac{\partial v}{\partial t} &= a \frac{\partial u}{\partial x}. \end{aligned} \tag{22}$$

Differencing with backward Euler we get

$$u^{n+1} = u^n + a \left(\frac{\partial v}{\partial x}\right)^{n+1},$$

$$v^{n+1} = v^n + a \left(\frac{\partial u}{\partial x}\right)^{n+1}.\tag{23}$$

Substitute the second equation into the first to obtain:

$$\left(I - a^2 \Delta t^2 \frac{\partial^2}{\partial \mathbf{x}^2}\right) u^{n+1} = u^n + a \Delta t \left(\frac{\partial v}{\partial \mathbf{x}}\right)^n \tag{24}$$

The matrix *M* above is relatively easy to invert and is amenable to multigrid. The Schur

Implicit Numerical Methods for Magnetohydrodynamics 69

where *PS* <sup>=</sup> *Dv* <sup>−</sup> *LM*−1*<sup>U</sup>* is the Schur complement. The exact Jacobian inverse require *<sup>M</sup>*−<sup>1</sup>

Multigrid is impractical for *PS* because of the *M*−<sup>1</sup> factor and hence some simplifications are desirable. For the velocity update and the corrector part in the above equations, we can treat

where *PSI* = *Dv* − Δ*tLU* and is block-diagonally dominant. Multigrid is employed to compute

In this section, we discuss yet another NK approach to resistive MHD. This section is essentially based on the work by Reynolds et al. (2006) in which they have developed a fully implicit Jacobian-Free NK method for compressible MHD. The main difference between this section and the previous one is in the preconditioning strategy employed during the Krylov

The resistive MHD equations are rewritten in a form which allows a method-of-lines

*qn* ∑ *i*=1 

where *R*(**U**) is defined using the divergence of the fluxes (both hyperbolic and diffusion terms) as in equation (1). The time-evolved state **U***<sup>n</sup>* solves the nonlinear residual equation *g*(**U**) = 0. *qn* determines the method's order of accuracy and at *qn* = {1, 2} the method is stable for any Δ*tn*, with stability decreases as *qn* increases. *αn*,*<sup>i</sup>* and *βn*,*<sup>i</sup>* are fixed parameters for a given method order *qn*. In this approach Δ*tn*, *qn* are adaptively chosen at each time step to balance solution accuracy, solver convergence, and temporal stability (Hindmarsh (2000)).

approach. Reynolds et al. use a BDF method (up to fifth order accurate):

*<sup>g</sup>*(**U***n*) <sup>≡</sup> **<sup>U</sup>***<sup>n</sup>* <sup>−</sup> <sup>Δ</sup>*tnβn*,0*R*(**U***n*) <sup>−</sup>

*<sup>δ</sup>v*<sup>∗</sup> <sup>=</sup> <sup>−</sup>*P*−<sup>1</sup>

*M*−1*L* 0

0 *P*−<sup>1</sup> *S*

*<sup>δ</sup>u*<sup>∗</sup> <sup>=</sup> <sup>−</sup>*M*−1F*<sup>u</sup>* (*Predictor*) (34)

*<sup>δ</sup><sup>u</sup>* <sup>=</sup> *<sup>δ</sup>u*<sup>∗</sup> <sup>−</sup> *<sup>M</sup>*−1*Uδ<sup>v</sup>* (*Corrector*). (36)

 *<sup>I</sup>* <sup>−</sup>*M*−1*<sup>U</sup>* 0 *I*

*<sup>S</sup>* [F*<sup>v</sup>* − *Lδv*∗] (*Velocity update*) (35)

*<sup>δ</sup>u*<sup>∗</sup> <sup>=</sup> <sup>−</sup>*M*−1F*<sup>u</sup>* (37)

*SI* [F*<sup>v</sup>* − *Lδv*∗] (38)

*δu* = *δu*<sup>∗</sup> − Δ*tUδv*, (39)

*αn*,*i***U***n*−*<sup>i</sup>* + Δ*tnβn*,*iR*(**U***n*−*<sup>i</sup>*

) 

, (40)

. (33)

complement analysis of the above 2 × 2 system is given below:

*<sup>δ</sup>v*<sup>∗</sup> <sup>=</sup> <sup>−</sup>*P*−<sup>1</sup>

*I* 0

<sup>−</sup>*LM*−<sup>1</sup> *<sup>I</sup>*

*<sup>S</sup>* . The following predictor-corrector algorithm is proposed.

<sup>−</sup><sup>1</sup> =

 *M U L Dv*

and *P*−<sup>1</sup>

*<sup>M</sup>*−<sup>1</sup> <sup>≈</sup> <sup>Δ</sup>*t*. This gives

the inverse of *PSI* and *M*.

step.

**3.4 NK method for resistive MHD II**

which is is much better conditioned because the parabolic operator is diagonally dominant. Multigrid techniques usually perform well on elliptic and parabolic operators do poorly on hyperbolic operators which are diagonal submissive.

We now turn to parabolization by the Schur complement approach.

$$
\begin{bmatrix} D\_1 & U \\ L & D\_2 \end{bmatrix} = \begin{bmatrix} I \ \mathcal{U} D\_2^{-1} \\ 0 & I \end{bmatrix} \begin{bmatrix} D\_1 - \mathcal{U} D\_2^{-1} L & 0 \\ 0 & D\_2 \end{bmatrix} \begin{bmatrix} I & 0 \\ D\_2^{-1} L & I \end{bmatrix} . \tag{25}
$$

Stiff off-diagonal blocks *L* and *U* are now shifted to the diagonal via the Schur complement *<sup>D</sup>*<sup>1</sup> <sup>−</sup> *UD*−<sup>1</sup> <sup>2</sup> *<sup>L</sup>*. Applied to the model system above, *<sup>D</sup>*<sup>1</sup> <sup>−</sup> *UD*−<sup>1</sup> <sup>2</sup> *L* = � *<sup>I</sup>* <sup>−</sup> *<sup>a</sup>*2Δ*<sup>t</sup>* 2 *∂*<sup>2</sup> *∂x*<sup>2</sup> � .

#### **3.3.2 Application to resistive MHD**

We begin by examining the linearized resistive MHD equations. These are written as

$$
\delta\rho = L\_{\theta}(\delta\rho, \delta\sigma) \tag{26}
$$

$$
\delta p = L\_p(\delta p, \delta \sigma) \tag{27}
$$

$$
\delta \mathbf{B} = L\_{\mathbf{B}}(\delta \mathbf{B}, \delta \mathbf{v}) \tag{28}
$$

$$
\delta\mathfrak{v} = L\_{\upsilon}(\delta\mathfrak{v}, \delta\mathbf{B}, \delta\mathfrak{p}, \delta\mathfrak{p}),
\tag{29}
$$

which illustrates the couplings between the various unknowns. In NK the Jacobian has the following coupling

$$J\delta\mathcal{U} = \begin{bmatrix} D\_{\rho} & 0 & 0 & \mathcal{U}\_{pv} \\ 0 & D\_p & 0 & \mathcal{U}\_{pv} \\ 0 & 0 & D\_B & \mathcal{U}\_{Bv} \\ L\_{\rho v} & L\_{pv} & L\_{Bv} & D\_v \end{bmatrix} \begin{bmatrix} \delta\rho \\ \delta p \\ \delta \mathbf{B} \\ \delta \mathbf{B} \\ \delta \mathbf{v} \end{bmatrix} . \tag{30}$$

which shows that the momentum equations are intimately coupled with other equations but that the density is only coupled with the velocity nonlinearly and so on. The diagonal blocks are of the "advection-diffusion" type and clearly amenable to multigrid and easily inverted. The off-diagonal terms denoted by *L* and *U* contain the hyperbolic couplings. The above Jacobian is rewritten as

$$J\delta U = \begin{bmatrix} M & U \\ L & D\_v \end{bmatrix} \begin{pmatrix} \delta \mathfrak{u} \\ \delta \mathfrak{v} \end{pmatrix},\tag{31}$$

where

$$
\delta\mathfrak{u} = \begin{pmatrix} \delta\rho \\ \delta p \\ \delta \mathbf{B} \end{pmatrix}, \quad M = \begin{pmatrix} D\_{\rho} & 0 & 0 \\ 0 & D\_{p} & 0 \\ 0 & 0 & D\_{B} \end{pmatrix}. \tag{32}
$$

10 Will-be-set-by-IN-TECH

which is is much better conditioned because the parabolic operator is diagonally dominant. Multigrid techniques usually perform well on elliptic and parabolic operators do poorly on

� �*D*<sup>1</sup> <sup>−</sup> *UD*−<sup>1</sup>

Stiff off-diagonal blocks *L* and *U* are now shifted to the diagonal via the Schur complement

which illustrates the couplings between the various unknowns. In NK the Jacobian has the

*D<sup>ρ</sup>* 0 0 *Uρ<sup>v</sup>* 0 *Dp* 0 *Upv* 0 0 *DB UBv Lρ<sup>v</sup> Lpv LBv Dv*

which shows that the momentum equations are intimately coupled with other equations but that the density is only coupled with the velocity nonlinearly and so on. The diagonal blocks are of the "advection-diffusion" type and clearly amenable to multigrid and easily inverted. The off-diagonal terms denoted by *L* and *U* contain the hyperbolic couplings. The above

> � � *δu δv*

⎛

⎜⎜⎝

�

⎞

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

*δρ δp δB δv*

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

We begin by examining the linearized resistive MHD equations. These are written as

<sup>2</sup> *L* 0 0 *D*<sup>2</sup>

� � *I* 0 *D*−<sup>1</sup> <sup>2</sup> *L I*

�

<sup>2</sup> *L* =

*δρ* = *Lρ*(*δρ*, *δv*) (26)

*δp* = *Lp*(*δp*, *δv*) (27)

*δB* = *LB*(*δB*, *δv*) (28)

*δv* = *Lv*(*δv*, *δB*, *δρ*, *δp*), (29)

�

*<sup>I</sup>* <sup>−</sup> *<sup>a</sup>*2Δ*<sup>t</sup>*

2 *∂*<sup>2</sup> *∂x*<sup>2</sup> � .

, (30)

, (31)

⎟⎟⎠ . (32)

. (25)

hyperbolic operators which are diagonal submissive.

� = �

� *D*<sup>1</sup> *U L D*<sup>2</sup>

**3.3.2 Application to resistive MHD**

*<sup>D</sup>*<sup>1</sup> <sup>−</sup> *UD*−<sup>1</sup>

following coupling

Jacobian is rewritten as

where

We now turn to parabolization by the Schur complement approach.

*I UD*−<sup>1</sup> 2 0 *I*

<sup>2</sup> *<sup>L</sup>*. Applied to the model system above, *<sup>D</sup>*<sup>1</sup> <sup>−</sup> *UD*−<sup>1</sup>

*JδU* =

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

*JδU* =

*δu* =

⎛

*δρ δp δB* ⎞

⎜⎜⎝

� *M U L Dv*

⎟⎟⎠ , *<sup>M</sup>* <sup>=</sup>

The matrix *M* above is relatively easy to invert and is amenable to multigrid. The Schur complement analysis of the above 2 × 2 system is given below:

$$
\begin{bmatrix} M \ U \\\\ L \ D\_{\overline{v}} \end{bmatrix}^{-1} = \begin{bmatrix} I & 0 \\\\ -LM^{-1} & I \end{bmatrix} \begin{bmatrix} M^{-1}L & 0 \\\\ 0 & P\_{\overline{S}}^{-1} \end{bmatrix} \begin{bmatrix} I - M^{-1}U \\\\ 0 & I \end{bmatrix} . \tag{33}
$$

where *PS* <sup>=</sup> *Dv* <sup>−</sup> *LM*−1*<sup>U</sup>* is the Schur complement. The exact Jacobian inverse require *<sup>M</sup>*−<sup>1</sup> and *P*−<sup>1</sup> *<sup>S</sup>* . The following predictor-corrector algorithm is proposed.

$$
\delta u^\* = -M^{-1} \mathcal{F}\_u \quad \text{(Predictor)}\tag{34}
$$

$$
\delta\boldsymbol{\sigma}\* = -P\_{\mathcal{S}}^{-1}[\mathcal{F}\_{\mathcal{V}} - L\delta\boldsymbol{\sigma}^\*] \quad \text{(Velocity } \upmu\text{)}\tag{35}
$$

$$
\delta\mathfrak{u} = \delta\mathfrak{u}^\* - M^{-1}\mathcal{U}\delta\mathfrak{v} \quad \text{(Corrector)}.\tag{36}
$$

Multigrid is impractical for *PS* because of the *M*−<sup>1</sup> factor and hence some simplifications are desirable. For the velocity update and the corrector part in the above equations, we can treat *<sup>M</sup>*−<sup>1</sup> <sup>≈</sup> <sup>Δ</sup>*t*. This gives

$$
\delta \mu^\* = -M^{-1} \mathcal{F}\_{\mu} \tag{37}
$$

$$
\delta\mathfrak{v}\* = -P\_{SI}^{-1}[\mathcal{F}\_{\mathcal{v}} - L\delta\mathfrak{v}^\*] \tag{38}
$$

$$
\delta\mathfrak{u} = \delta\mathfrak{u}^\* - \Delta t \mathsf{U} \delta\mathfrak{v}, \tag{39}
$$

where *PSI* = *Dv* − Δ*tLU* and is block-diagonally dominant. Multigrid is employed to compute the inverse of *PSI* and *M*.

#### **3.4 NK method for resistive MHD II**

In this section, we discuss yet another NK approach to resistive MHD. This section is essentially based on the work by Reynolds et al. (2006) in which they have developed a fully implicit Jacobian-Free NK method for compressible MHD. The main difference between this section and the previous one is in the preconditioning strategy employed during the Krylov step.

The resistive MHD equations are rewritten in a form which allows a method-of-lines approach. Reynolds et al. use a BDF method (up to fifth order accurate):

$$\log(\mathbf{U}^{n}) \equiv \mathbf{U}^{n} - \Delta t\_{n} \beta\_{n,0} \mathbb{R}(\mathbf{U}^{n}) - \sum\_{i=1}^{q\_{n}} \left[ \alpha\_{n,i} \mathbf{U}^{n-i} + \Delta t\_{n} \beta\_{n,i} \mathbb{R}(\mathbf{U}^{n-i}) \right],\tag{40}$$

where *R*(**U**) is defined using the divergence of the fluxes (both hyperbolic and diffusion terms) as in equation (1). The time-evolved state **U***<sup>n</sup>* solves the nonlinear residual equation *g*(**U**) = 0. *qn* determines the method's order of accuracy and at *qn* = {1, 2} the method is stable for any Δ*tn*, with stability decreases as *qn* increases. *αn*,*<sup>i</sup>* and *βn*,*<sup>i</sup>* are fixed parameters for a given method order *qn*. In this approach Δ*tn*, *qn* are adaptively chosen at each time step to balance solution accuracy, solver convergence, and temporal stability (Hindmarsh (2000)).

Since MHD stiffness results from fast hyperbolic and diffusive effects, we set

*<sup>h</sup> <sup>P</sup>*−<sup>1</sup>

*<sup>d</sup>* <sup>=</sup> *<sup>J</sup>*(**U**)−<sup>1</sup> <sup>+</sup> <sup>O</sup>(Δ*<sup>t</sup>*

This operator-splitting approach, widely used as a stand-alone solver, is used to accelerate

Implicit Numerical Methods for Magnetohydrodynamics 71

*Ph***: Ideal MHD Preconditioner:** The ideal MHD preconditioner discussed here essentially exploits the local wave structure of the underlying hyperbolic portion of the PDEs. Hence this approach may be dubbed a "wave-structure"-based preconditioner. For linear multistep time integration approaches, it is convenient to first rewrite the nonlinear problem (40) in the form

where the terms *F*(*U*), *G*(*U*) and *H*(*U*) denote the *x*, *y* and *z* directional hyperbolic fluxes, and the term *g* incorporates previous time-level information into the discretized problem.

[*I* + *γ∂x*(*JF*(*U*)(·))] *V* = *V* + *γ∂x*(*JF*(*U*)*V*).

Omitting the explicit dependence on *U* from the notation, and introducing nonsingular

 *JF L*−<sup>1</sup> *F LF*(·)

 *JGL*−<sup>1</sup> *G LG*(·)

 *JH L*−<sup>1</sup> *H LH*(·) 

 *JGL*−<sup>1</sup> *G* 

*<sup>I</sup>* <sup>+</sup> *<sup>γ</sup>JGL*−<sup>1</sup>

 *JGL*<sup>−</sup><sup>1</sup> *G* 

The preconditioning scheme in this approach is based on the assumption that the majority of the stiffness found in the Jacobian is a result of a small number of very fast hyperbolic waves. To develop an approach for separately treating only these fast waves, consider the preconditioning matrix, *P*, constructed using a directional and operator-based splitting of J,

*∂xF*(*U*) + *∂yG*(*U*) + *∂zH*(*U*)

*<sup>∂</sup>x*(*JF*(*U*)(·)) + *<sup>∂</sup>y*(*JG*(*U*)(·)) + *<sup>∂</sup>z*(*JH*(*U*)(·))

*<sup>∂</sup><sup>U</sup> F*(*U*). Using the notation (·) to denote the location at which the action

 + *∂<sup>z</sup> JH L*−<sup>1</sup>

*<sup>G</sup> <sup>∂</sup><sup>y</sup>* (*LG*(·)) <sup>+</sup> *JH <sup>L</sup>*−<sup>1</sup>

*<sup>G</sup> ∂<sup>y</sup>* (*LG*(·))

*LG* + *γ∂<sup>z</sup>*

 *JH L*−<sup>1</sup> *H LH* 

*LG*(·) + *∂<sup>z</sup>*

2).

+ *g* = 0, (45)

*<sup>H</sup> LH*(·)

*<sup>H</sup> ∂<sup>z</sup>* (*LH*(·))

 *JH L*−<sup>1</sup> *H LH*(·) .

*<sup>I</sup>* <sup>+</sup> *<sup>γ</sup>JH <sup>L</sup>*−<sup>1</sup>

*<sup>H</sup> ∂<sup>z</sup>* (*LH*(·))

(47)

, (46)

*P*−<sup>1</sup> = *P*−<sup>1</sup>

convergence of the more stable and accurate implicit NK approach.

*f*(*U*) = *U* + *γ*

matrices *LF*, *LG* and *LH*, re-write the Jacobian system (46) as

 + *∂<sup>y</sup> JGL*−<sup>1</sup> *<sup>G</sup> LG*(·)

*<sup>F</sup> ∂<sup>x</sup>* (*LF*(·)) + *∂<sup>x</sup>*

*<sup>G</sup> ∂<sup>y</sup>* (*LG*(·)) + *∂<sup>y</sup>*

*<sup>H</sup> ∂<sup>z</sup>* (*LH*(·)) + *∂<sup>z</sup>*

*<sup>F</sup> <sup>∂</sup><sup>x</sup>* (*LF*(·)) <sup>+</sup> *JGL*−<sup>1</sup>

*LF* + *γ∂<sup>y</sup>*

*LF*(·) + *∂<sup>y</sup>*

This nonlinear problem has Jacobian

of the linear operator takes place, e.g.

*J* = *I* + *γ*

= *I* + *γ*

= *I* + *γ*

*P* = 

> *I* + *γ∂<sup>x</sup>*

 *∂x JF L*−<sup>1</sup> *<sup>F</sup> LF*(·)

 *JF L*−<sup>1</sup>

+ *JGL*−<sup>1</sup>

+*JHL*−<sup>1</sup>

 *JF L*−<sup>1</sup>

+*∂<sup>x</sup> JF L*−<sup>1</sup> *F* 

*I* + *γJF L*−<sup>1</sup>

= *J* + *O*(*γ*2).

 *JF L*−<sup>1</sup> *F* 

*<sup>F</sup> ∂<sup>x</sup>* (*LF*(·))

with, e.g., *JF*(*U*) = *<sup>∂</sup>*

*J*(*U*) = *I* + *γ*

Alternatively one may also use a *θ*-scheme

$$\log(\mathcal{U}^{\boldsymbol{\eta}}) = \mathcal{U}^{\boldsymbol{\eta}} - \mathcal{U}^{\boldsymbol{\eta}-1} + \Delta t \left( \theta \mathcal{R}(\mathcal{U}^{\boldsymbol{\eta}}) + (1 - \theta) \mathcal{R}(\mathcal{U}^{\boldsymbol{\eta}-1}) \right), \tag{41}$$

where *θ* = 0.5 corresponds to a Crank-Nicholson approach. The inexact Jacobian-Free NK approach is adopted to solve the nonlinear function *g*(**U***n*). The divergence of the fluxes in (1) is discretized using the following finite difference form

$$
\left(\frac{\partial f}{\partial \mathbf{x}}\right)\_{i,j,k} = \frac{\vec{f}\_{i+\frac{1}{2},j,k} - \vec{f}\_{i-\frac{1}{2},j,k}}{\Delta \mathbf{x}},\tag{42}
$$

where *f* may represent either the hyperbolic or the parabolic fluxes, and Δ*x* is the mesh spacing in the x-direction (assumed uniform). The quantity ˜ *f i*+ <sup>1</sup> <sup>2</sup> ,*j*,*<sup>k</sup>* is referred to as the numerical flux through the face {*<sup>i</sup>* <sup>+</sup> <sup>1</sup> <sup>2</sup> , *j*, *k*} and is computed as a linear combination of the fluxes at cell centers as

$$\tilde{f}\_{i+\frac{1}{2},j,k} = \sum\_{\nu=-m}^{n} a\_{\nu} f\_{i+\nu,j,k}.\tag{43}$$

Reynolds et al. give the options for several spatial difference schemes. For a second-order central difference implementation, *m* = 0, *n* = 1 and *a*<sup>0</sup> = *a*<sup>1</sup> = <sup>1</sup> <sup>2</sup> ; for a fourth-order central difference approximation, *<sup>m</sup>* <sup>=</sup> 1, *<sup>n</sup>* <sup>=</sup> 2, and *<sup>a</sup>*−<sup>1</sup> <sup>=</sup> *<sup>a</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup><sup>1</sup> <sup>12</sup> , *<sup>a</sup>*<sup>0</sup> <sup>=</sup> *<sup>a</sup>*<sup>1</sup> <sup>=</sup> <sup>7</sup> <sup>12</sup> ; and for tuned second-order central differences, *<sup>a</sup>*−<sup>1</sup> = *<sup>a</sup>*<sup>2</sup> = −0.197, *<sup>a</sup>*<sup>0</sup> = *<sup>a</sup>*<sup>1</sup> = 0.697 (Hill & Pullin (2004)). These central difference approximations are free of dissipation errors, except perhaps near domain boundaries. They do, however, suffer from dispersion errors. Consequently, physical phenomena that are not well resolved can suffer from ringing. The dispersion errors can be minimized by using schemes such as the tuned-second order scheme, mentioned above, which has lower dispersion error than the central difference schemes. The numerical approximation to the divergence ∇ · *B* is written as

$$\begin{split} \nabla \cdot \mathbf{B} &= \frac{\tilde{\mathcal{B}}\_{i+\frac{1}{2},j,k}^{x} - \tilde{\mathcal{B}}\_{i-\frac{1}{2},j,k}^{x}}{\Delta x} + \frac{\tilde{\mathcal{B}}\_{i,j+\frac{1}{2},k}^{y} - \tilde{\mathcal{B}}\_{i,j-\frac{1}{2},k}^{y}}{\Delta y} + \frac{\tilde{\mathcal{B}}\_{i,j,k+\frac{1}{2}}^{z} - \tilde{\mathcal{B}}\_{i,j,k-\frac{1}{2}}^{z}}{\Delta z} \\ &+ \mathcal{O}(\Delta x^{p}) + \mathcal{O}(\Delta y^{p}) + \mathcal{O}(\Delta z^{p}) \end{split} \tag{44}$$

where *B<sup>α</sup>* is the *α*-component of the magnetic field, and the terms *B*˜ *<sup>α</sup>* are evaluated as shown in equation (43), and *p* is the order of the spatial derivatives. If the numerical approximation of ∇ · *B* is ensured to be zero at *t* = 0 then it can be easily shown that the numerical fluxes, as computed above, ensure the solenoidal property of the magnetic field in the discrete sense is automatically satisfied. This conservation property of preserving the solenoidal nature of the magnetic field in an implicit method is generally very desirable.

#### **3.4.1 Preconditioner formulation**

The preconditioner strategy, overall, uses an operator split approach to separate the wave-dominated portion from the diffusion portion. Instead of solving *J δ***U** = −*g*, we solve the related system (*JP*−1)(*<sup>P</sup> <sup>δ</sup>***U**) = <sup>−</sup>*g*, i.e., the right preconditioning approach is adopted. 12 Will-be-set-by-IN-TECH

where *θ* = 0.5 corresponds to a Crank-Nicholson approach. The inexact Jacobian-Free NK approach is adopted to solve the nonlinear function *g*(**U***n*). The divergence of the fluxes in (1)

where *f* may represent either the hyperbolic or the parabolic fluxes, and Δ*x* is the mesh

*n* ∑*ν*=−*m*

Reynolds et al. give the options for several spatial difference schemes. For a second-order

for tuned second-order central differences, *<sup>a</sup>*−<sup>1</sup> = *<sup>a</sup>*<sup>2</sup> = −0.197, *<sup>a</sup>*<sup>0</sup> = *<sup>a</sup>*<sup>1</sup> = 0.697 (Hill & Pullin (2004)). These central difference approximations are free of dissipation errors, except perhaps near domain boundaries. They do, however, suffer from dispersion errors. Consequently, physical phenomena that are not well resolved can suffer from ringing. The dispersion errors can be minimized by using schemes such as the tuned-second order scheme, mentioned above, which has lower dispersion error than the central difference schemes. The numerical

<sup>2</sup> ,*<sup>k</sup>* <sup>−</sup> *<sup>B</sup>*˜ *<sup>y</sup>*

Δ*y*

where *B<sup>α</sup>* is the *α*-component of the magnetic field, and the terms *B*˜ *<sup>α</sup>* are evaluated as shown in equation (43), and *p* is the order of the spatial derivatives. If the numerical approximation of ∇ · *B* is ensured to be zero at *t* = 0 then it can be easily shown that the numerical fluxes, as computed above, ensure the solenoidal property of the magnetic field in the discrete sense is automatically satisfied. This conservation property of preserving the solenoidal nature of the

The preconditioner strategy, overall, uses an operator split approach to separate the wave-dominated portion from the diffusion portion. Instead of solving *J δ***U** = −*g*, we solve the related system (*JP*−1)(*<sup>P</sup> <sup>δ</sup>***U**) = <sup>−</sup>*g*, i.e., the right preconditioning approach is adopted.

*<sup>i</sup>*,*j*−<sup>1</sup> 2 ,*k*

<sup>+</sup>O(Δ*xp*) + <sup>O</sup>(Δ*yp*) + <sup>O</sup>(Δ*zp*) (44)

+ *B*˜ *z i*,*j*,*k*+<sup>1</sup> 2 <sup>−</sup> *<sup>B</sup>*˜ *<sup>z</sup> <sup>i</sup>*,*j*,*k*−<sup>1</sup> 2

*<sup>θ</sup>R*(*Un*)+(<sup>1</sup> <sup>−</sup> *<sup>θ</sup>*)*R*(*Un*−1)

<sup>Δ</sup>*<sup>x</sup>* , (42)

*a<sup>ν</sup> fi*<sup>+</sup>*ν*,*j*,*k*. (43)

Δ*z*

*f i*+ <sup>1</sup>

<sup>2</sup> , *j*, *k*} and is computed as a linear combination of the

, (41)

<sup>2</sup> ,*j*,*<sup>k</sup>* is referred to as the

<sup>2</sup> ; for a fourth-order

<sup>12</sup> ; and

<sup>12</sup> , *<sup>a</sup>*<sup>0</sup> <sup>=</sup> *<sup>a</sup>*<sup>1</sup> <sup>=</sup> <sup>7</sup>

Alternatively one may also use a *θ*-scheme

numerical flux through the face {*<sup>i</sup>* <sup>+</sup> <sup>1</sup>

approximation to the divergence ∇ · *B* is written as

<sup>2</sup> ,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>B</sup>*˜ *<sup>x</sup>*

Δ*x*

magnetic field in an implicit method is generally very desirable.

*<sup>i</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> ,*j*,*k* + *B*˜ *y i*,*j*+ <sup>1</sup>

*B*˜ *x i*+ <sup>1</sup>

∇ · *B* =

**3.4.1 Preconditioner formulation**

fluxes at cell centers as

*<sup>g</sup>*(*Un*) = *<sup>U</sup><sup>n</sup>* <sup>−</sup> *<sup>U</sup>n*−<sup>1</sup> <sup>+</sup> <sup>Δ</sup>*<sup>t</sup>*

*∂ f ∂x* 

spacing in the x-direction (assumed uniform). The quantity ˜

˜ *f i*+ <sup>1</sup> <sup>2</sup> ,*j*,*<sup>k</sup>* <sup>=</sup>

central difference implementation, *m* = 0, *n* = 1 and *a*<sup>0</sup> = *a*<sup>1</sup> = <sup>1</sup>

central difference approximation, *<sup>m</sup>* <sup>=</sup> 1, *<sup>n</sup>* <sup>=</sup> 2, and *<sup>a</sup>*−<sup>1</sup> <sup>=</sup> *<sup>a</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup><sup>1</sup>

*i*,*j*,*k* = ˜ *f i*+<sup>1</sup> <sup>2</sup> ,*j*,*<sup>k</sup>* <sup>−</sup> ˜ *f <sup>i</sup>*−<sup>1</sup> <sup>2</sup> ,*j*,*k*

is discretized using the following finite difference form

Since MHD stiffness results from fast hyperbolic and diffusive effects, we set

$$P^{-1} = P\_h^{-1} P\_d^{-1} = f(\mathbf{U})^{-1} + \mathcal{O}(\Delta t^2).$$

This operator-splitting approach, widely used as a stand-alone solver, is used to accelerate convergence of the more stable and accurate implicit NK approach.

*Ph***: Ideal MHD Preconditioner:** The ideal MHD preconditioner discussed here essentially exploits the local wave structure of the underlying hyperbolic portion of the PDEs. Hence this approach may be dubbed a "wave-structure"-based preconditioner. For linear multistep time integration approaches, it is convenient to first rewrite the nonlinear problem (40) in the form

$$f(\mathcal{U}) = \mathcal{U} + \gamma \left[\partial\_{\mathcal{X}} F(\mathcal{U}) + \partial\_{\mathcal{Y}} G(\mathcal{U}) + \partial\_{\mathcal{Z}} H(\mathcal{U})\right] + \mathcal{g} = 0,\tag{45}$$

where the terms *F*(*U*), *G*(*U*) and *H*(*U*) denote the *x*, *y* and *z* directional hyperbolic fluxes, and the term *g* incorporates previous time-level information into the discretized problem. This nonlinear problem has Jacobian

$$J(\mathcal{U}) = I + \gamma \left[ \partial\_x \left( f\_F(\mathcal{U})(\cdot) \right) + \partial\_y \left( f\_G(\mathcal{U})(\cdot) \right) + \partial\_z \left( f\_H(\mathcal{U})(\cdot) \right) \right],\tag{46}$$

with, e.g., *JF*(*U*) = *<sup>∂</sup> <sup>∂</sup><sup>U</sup> F*(*U*). Using the notation (·) to denote the location at which the action of the linear operator takes place, e.g.

$$\left[I + \gamma \partial\_{\mathfrak{X}} (J\_F(\mathcal{U})(\cdot))\right] V = V + \gamma \partial\_{\mathfrak{X}} (J\_F(\mathcal{U})V).$$

Omitting the explicit dependence on *U* from the notation, and introducing nonsingular matrices *LF*, *LG* and *LH*, re-write the Jacobian system (46) as

*J* = *I* + *γ ∂x JF L*−<sup>1</sup> *<sup>F</sup> LF*(·) + *∂<sup>y</sup> JGL*−<sup>1</sup> *<sup>G</sup> LG*(·) + *∂<sup>z</sup> JH L*−<sup>1</sup> *<sup>H</sup> LH*(·) = *I* + *γ JF L*−<sup>1</sup> *<sup>F</sup> ∂<sup>x</sup>* (*LF*(·)) + *∂<sup>x</sup> JF L*−<sup>1</sup> *F LF*(·) + *JGL*−<sup>1</sup> *<sup>G</sup> ∂<sup>y</sup>* (*LG*(·)) + *∂<sup>y</sup> JGL*−<sup>1</sup> *G LG*(·) +*JHL*−<sup>1</sup> *<sup>H</sup> ∂<sup>z</sup>* (*LH*(·)) + *∂<sup>z</sup> JH L*−<sup>1</sup> *H LH*(·) = *I* + *γ JF L*−<sup>1</sup> *<sup>F</sup> <sup>∂</sup><sup>x</sup>* (*LF*(·)) <sup>+</sup> *JGL*−<sup>1</sup> *<sup>G</sup> <sup>∂</sup><sup>y</sup>* (*LG*(·)) <sup>+</sup> *JH <sup>L</sup>*−<sup>1</sup> *<sup>H</sup> ∂<sup>z</sup>* (*LH*(·)) +*∂<sup>x</sup> JF L*−<sup>1</sup> *F LF*(·) + *∂<sup>y</sup> JGL*−<sup>1</sup> *G LG*(·) + *∂<sup>z</sup> JH L*−<sup>1</sup> *H LH*(·) .

The preconditioning scheme in this approach is based on the assumption that the majority of the stiffness found in the Jacobian is a result of a small number of very fast hyperbolic waves. To develop an approach for separately treating only these fast waves, consider the preconditioning matrix, *P*, constructed using a directional and operator-based splitting of J,

$$\begin{split} P &= \left[ I + \gamma I\_F L\_F^{-1} \eth\_X \left( L\_F(\cdot) \right) \right] \left[ I + \gamma I\_G L\_G^{-1} \eth\_Y \left( L\_G(\cdot) \right) \right] \left[ I + \gamma I\_H L\_H^{-1} \eth\_Z \left( L\_H(\cdot) \right) \right] \\ & \left[ I + \gamma \eth\_X \left( I\_F L\_F^{-1} \right) L\_F + \gamma \eth\_Y \left( I\_G L\_G^{-1} \right) L\_G + \gamma \eth\_Z \left( I\_H L\_H^{-1} \right) L\_H \right] \\ &= I + O(\gamma^2). \end{split} \tag{47}$$

characteristic variable *wk*. In cases where the wave speeds can be estimated, a pre-defined cutoff to the number of waves included in the preconditioner can be set. This reduction allows for significant savings in preconditioner computation. For those waves that are not preconditioned, approximate them as having wave speed equal to zero, i.e. solving with the approximation <sup>Λ</sup><sup>ˆ</sup> *<sup>F</sup>* <sup>=</sup> *diag*(*λ*1,..., *<sup>λ</sup>q*, 0, . . . , 0). Omission of the (*<sup>n</sup>* <sup>−</sup> *<sup>q</sup>*) slowest waves in this fashion amounts to a further approximation of the preconditioner to the original discretized

Implicit Numerical Methods for Magnetohydrodynamics 73

*χ* + *γJFRF∂x*(*LFχ*) = *b*, *χ*ˆ + *γ* ˆ*JFRF∂x*(*LFχ*ˆ) = *b*,

*w* + *γ*Λ*F∂xw* = *LFb*, *w*ˆ + *γ*Λˆ *<sup>F</sup>∂xw*ˆ = *LFb*,

*w*ˆ *<sup>k</sup>* = (*LFb*)*<sup>k</sup>*

Since the eigenvector matrices *LF* and *RF* may be renormalized as desired, and the eigenvalues are ordered so that *λ<sup>i</sup>* ≥ *λj*, for *i* < *j*, the dominant error from preconditioning only the *q*

> <sup>|</sup>*γλq*<sup>+</sup>1/Δ*x*<sup>|</sup> <sup>1</sup> − |*γλq*<sup>+</sup>1/Δ*x*<sup>|</sup>

Hence omission of waves with small eigenvalues compared to the dynamical time scale

The remaining component of the split preconditioner (47) comprises the local system *Plocalu* =

*I* + *γ∂<sup>x</sup>* (*JFRF*) *LF* + *γ∂<sup>y</sup>* (*JGRG*) *LG* + *γ∂<sup>z</sup>* (*JHRH*) *LH*

*I* + *γ∂<sup>x</sup>* (*RF*Λ*F*) *LF* + *γ∂<sup>y</sup>* (*RG*Λ*G*) *LG* + *γ∂<sup>z</sup>* (*RH*Λ*H*) *LH*

Note that for spatially homogeneous Jacobians, *∂x*(*RF*Λ*F*) = 0 (similarly for *y* and *z*), so this system reduces to *u* = *v*. Hence this component may optionally be included to correct for spatial inhomogeneity in *JF*, *JG* and *JH*. In keeping with the previous discretization

These solves are spatially decoupled (with respect to *u*), resulting in a block-diagonal matrix

(i.e. *γλ* � 1) will not significantly affect preconditioner accuracy.

**Local Non-Constant Coefficient Correction Solve:**

.

*RF*,*i*+1Λ*F*,*i*+<sup>1</sup> − *RF*,*i*−1Λ*F*,*i*−<sup>1</sup>

where ˆ*JF* = *RF*Λˆ *<sup>F</sup> LF*. Left-multiplying by *LF* and proceeding as before, to obtain

*w<sup>k</sup>* + *γλk∂xw<sup>k</sup>* = (*LFb*)*<sup>k</sup>*

*w*ˆ *<sup>k</sup>* + *γλk∂xw*ˆ *<sup>k</sup>* = (*LFb*)*<sup>k</sup>*

*<sup>F</sup>* as the *x*-directional preconditioner based on *q* waves, consider

, *k* = 1, . . . , *n*

, *k* = 1, . . . , *q*

, *k* = *q* + 1, . . . , *n*.

 *u* = *v*

 *LF*,*i*.

 *u* = *v*.

*<sup>F</sup>χ*ˆ = *b*, i.e.

PDE system. Writing *P*ˆ

�*<sup>χ</sup>* <sup>−</sup> *<sup>χ</sup>*ˆ�*p*, where *<sup>χ</sup>* solves *PF<sup>χ</sup>* <sup>=</sup> *<sup>b</sup>* and *<sup>χ</sup>*<sup>ˆ</sup> solves *<sup>P</sup>*<sup>ˆ</sup>

⇔

fastest waves is approximately

approaches, approximate this system as, e.g.

*γ∂<sup>x</sup>* (*RF*Λ*F*) *LF* <sup>≈</sup> *<sup>γ</sup>*

2Δ*x* 

whose solution requires only *n* × *n* dense linear solves at each spatial location.

⇔

*v*,

Denote these components as *P* = *PFPGPHPlocal*. Through constructing the operator *P* as a product in this manner, the preconditioner solve consists of 3 simpler, 1-dimensional implicit advection problems, along with one additional correction for spatial variations in the directional Jacobians *JF*, *JG* and *JH*. Hence, *Pu* = *b* may be solved via the steps (i) *PF χ* = *b*, (ii) *PG w* = *χ*, (iii) *PH v* = *w*, and (iv) *Plocal u* = *v*. Note that the splitting (47) is not unique, and that in fact these operations can be applied in any order. The technique for efficient solution of each of the above systems is presented in the ensuing paragraphs.

#### **Directional Preconditioner Solves:**

First consider solution of the three preconditioning systems *PF*, *PG* and *PH* from (47) of the form, e.g. (*x*-direction)

$$P\_{\mathcal{F}}\chi = b \quad \Leftrightarrow \quad \chi + \gamma \, J\_{\mathcal{F}} L\_{\mathcal{F}}^{-1} \partial\_{\mathcal{X}} \left( L\_{\mathcal{F}} \chi \right) = b. \tag{48}$$

To this point *LF*, *LG*, and *LH* are still unspecified. These are *n* × *n* matrices (*n* = 7 or 8 for compressible MHD depending upon whether the seven- or eight-wave formulation is used) whose rows are the left eigenvectors of the respective Jacobians, giving the identities,

$$L\_{\sf F} \sf I\_{\sf F} = \Lambda\_{\sf F} L\_{\sf F} \qquad \Lambda\_{\sf F} = \mathsf{diag}(\lambda^1, \dots, \lambda^n)\_{\sf r} \qquad \sf I\_{\sf F} \mathsf{R}\_{\sf F} = \mathsf{R}\_{\sf F} \Lambda\_{\sf F}$$

where *RF* <sup>≡</sup> *<sup>L</sup>*−<sup>1</sup> *<sup>F</sup>* are the right eigenvectors (*<sup>n</sup>* <sup>×</sup> *<sup>n</sup>* column matrix), and *<sup>λ</sup><sup>k</sup>* are the eigenvalues of *JF*. Through pre-multiplication of (48) by *LF*, gives

$$L\_F \chi + \gamma \, L\_F \rfloor\_F R\_F \partial\_{\chi} \left( L\_F \chi \right) = L\_F b \qquad \Leftrightarrow \qquad L\_F \chi + \gamma \, \Lambda\_F \partial\_{\chi} \left( L\_F \chi \right) = L\_F b.$$

Defining the vector of characteristic variables *w* = *LFχ*, decouple the equations as ,

$$\{\!\!\!\/:w + \gamma\Lambda\_F\partial\_\mathbf{x}w = L\_Fb\} \qquad \Leftrightarrow \qquad w^k + \gamma\lambda^k \partial\_\mathbf{x} w^k = \mathfrak{G}^k, \; k = 1, \dots, n\_k$$

where *w<sup>k</sup>* denotes the *k*-th element of the characteristic vector *w*, and *β* = *LFb*.

Spatial discretization of each of the characteristic variables *w<sup>k</sup>* in the same manner as the original PDE (1), results in a tightly-banded linear system of equations (tridiagonal, pentadiagonal, etc., depending on the method), to solve for the values *w<sup>k</sup> <sup>j</sup>* . For example the tridiagonal version due to a *O*(Δ*x*2) finite-difference discretization is

$$w\_{\dot{j}}^{k} + \frac{\gamma \lambda\_{\dot{j}}^{k}}{2\Delta x} \left(w\_{\dot{j}+1}^{k} - w\_{\dot{j}-1}^{k}\right) = \beta\_{\dot{j}}^{k}.\tag{49}$$

Reynolds et al. use a second order centered finite-volume approximation, with resulting systems for each *w<sup>k</sup>* that are tridiagonal. Moreover, the above approach results not only in tridiagonal systems for each characteristic variable *wk*, but the systems are in fact *block tridiagonal*, where each block corresponds to only one spatial {*x*, *y*, *z*} row that is decoupled from all other rows through the domain in the same direction. Thus solution of these linear systems can be very efficient, as the computations on each row may be performed independently of one another.

Furthermore, since the initial assumption was that the stiffness of the overall system resulted from a few very fast waves, one may not construct and solve the above systems for each 14 Will-be-set-by-IN-TECH

Denote these components as *P* = *PFPGPHPlocal*. Through constructing the operator *P* as a product in this manner, the preconditioner solve consists of 3 simpler, 1-dimensional implicit advection problems, along with one additional correction for spatial variations in the directional Jacobians *JF*, *JG* and *JH*. Hence, *Pu* = *b* may be solved via the steps (i) *PF χ* = *b*, (ii) *PG w* = *χ*, (iii) *PH v* = *w*, and (iv) *Plocal u* = *v*. Note that the splitting (47) is not unique, and that in fact these operations can be applied in any order. The technique for efficient solution

First consider solution of the three preconditioning systems *PF*, *PG* and *PH* from (47) of the

To this point *LF*, *LG*, and *LH* are still unspecified. These are *n* × *n* matrices (*n* = 7 or 8 for compressible MHD depending upon whether the seven- or eight-wave formulation is used)

*LF JF* = Λ*<sup>F</sup> LF*, Λ*<sup>F</sup>* = diag(*λ*1,..., *λn*), *JFRF* = *RF*Λ*F*,

*LFχ* + *γ LF JFRF∂<sup>x</sup>* (*LFχ*) = *LFb* ⇔ *LFχ* + *γ* Λ*F∂<sup>x</sup>* (*LFχ*) = *LFb*.

*<sup>w</sup>* <sup>+</sup> *<sup>γ</sup>*Λ*F∂xw* <sup>=</sup> *LFb* <sup>⇔</sup> *<sup>w</sup><sup>k</sup>* <sup>+</sup> *γλk∂xw<sup>k</sup>* <sup>=</sup> *<sup>β</sup>k*, *<sup>k</sup>* <sup>=</sup> 1, . . . , *<sup>n</sup>*,

Spatial discretization of each of the characteristic variables *w<sup>k</sup>* in the same manner as the original PDE (1), results in a tightly-banded linear system of equations (tridiagonal,

Reynolds et al. use a second order centered finite-volume approximation, with resulting systems for each *w<sup>k</sup>* that are tridiagonal. Moreover, the above approach results not only in tridiagonal systems for each characteristic variable *wk*, but the systems are in fact *block tridiagonal*, where each block corresponds to only one spatial {*x*, *y*, *z*} row that is decoupled from all other rows through the domain in the same direction. Thus solution of these linear systems can be very efficient, as the computations on each row may be performed

Furthermore, since the initial assumption was that the stiffness of the overall system resulted from a few very fast waves, one may not construct and solve the above systems for each

*<sup>F</sup>* are the right eigenvectors (*<sup>n</sup>* <sup>×</sup> *<sup>n</sup>* column matrix), and *<sup>λ</sup><sup>k</sup>* are the eigenvalues

*j*−1 = *β<sup>k</sup>*

*<sup>F</sup> ∂<sup>x</sup>* (*LFχ*) = *b*. (48)

*<sup>j</sup>* . For example the

*<sup>j</sup>* . (49)

*PF<sup>χ</sup>* <sup>=</sup> *<sup>b</sup>* <sup>⇔</sup> *<sup>χ</sup>* <sup>+</sup> *<sup>γ</sup> JF <sup>L</sup>*−<sup>1</sup>

whose rows are the left eigenvectors of the respective Jacobians, giving the identities,

Defining the vector of characteristic variables *w* = *LFχ*, decouple the equations as ,

where *w<sup>k</sup>* denotes the *k*-th element of the characteristic vector *w*, and *β* = *LFb*.

pentadiagonal, etc., depending on the method), to solve for the values *w<sup>k</sup>*

*γλ<sup>k</sup> j* 2Δ*x wk <sup>j</sup>*+<sup>1</sup> <sup>−</sup> *<sup>w</sup><sup>k</sup>*

tridiagonal version due to a *O*(Δ*x*2) finite-difference discretization is

*wk <sup>j</sup>* +

of each of the above systems is presented in the ensuing paragraphs.

of *JF*. Through pre-multiplication of (48) by *LF*, gives

**Directional Preconditioner Solves:**

form, e.g. (*x*-direction)

where *RF* <sup>≡</sup> *<sup>L</sup>*−<sup>1</sup>

independently of one another.

characteristic variable *wk*. In cases where the wave speeds can be estimated, a pre-defined cutoff to the number of waves included in the preconditioner can be set. This reduction allows for significant savings in preconditioner computation. For those waves that are not preconditioned, approximate them as having wave speed equal to zero, i.e. solving with the approximation <sup>Λ</sup><sup>ˆ</sup> *<sup>F</sup>* <sup>=</sup> *diag*(*λ*1,..., *<sup>λ</sup>q*, 0, . . . , 0). Omission of the (*<sup>n</sup>* <sup>−</sup> *<sup>q</sup>*) slowest waves in this fashion amounts to a further approximation of the preconditioner to the original discretized PDE system. Writing *P*ˆ *<sup>F</sup>* as the *x*-directional preconditioner based on *q* waves, consider �*<sup>χ</sup>* <sup>−</sup> *<sup>χ</sup>*ˆ�*p*, where *<sup>χ</sup>* solves *PF<sup>χ</sup>* <sup>=</sup> *<sup>b</sup>* and *<sup>χ</sup>*<sup>ˆ</sup> solves *<sup>P</sup>*<sup>ˆ</sup> *<sup>F</sup>χ*ˆ = *b*, i.e.

$$(\chi + \gamma I\_F R\_F \partial\_\chi (L\_F \chi) = b, \qquad \hat{\chi} + \gamma \hat{f}\_F R\_F \partial\_\chi (L\_F \hat{\chi}) = b\_\nu$$

where ˆ*JF* = *RF*Λˆ *<sup>F</sup> LF*. Left-multiplying by *LF* and proceeding as before, to obtain

$$w + \gamma \Lambda\_F \partial\_x w = L\_F b\_\prime \qquad \Updot{w} + \gamma \hat{\Lambda}\_F \partial\_x \grave{w} = L\_F b\_\prime$$

$$\Leftrightarrow$$

$$w^k + \gamma \lambda^k \partial\_x w^k = \left(L\_F b\right)^k \quad k = 1, \dots, n$$

$$\not\ni\_\prime^k + \gamma \lambda^k \partial\_x \grave{w}^k = \left(L\_F b\right)^k \quad k = 1, \dots, q$$

$$\not\ni\_\prime^k = \left(L\_F b\right)^k \quad k = q+1, \dots, n.$$

Since the eigenvector matrices *LF* and *RF* may be renormalized as desired, and the eigenvalues are ordered so that *λ<sup>i</sup>* ≥ *λj*, for *i* < *j*, the dominant error from preconditioning only the *q* fastest waves is approximately

$$\frac{|\gamma\lambda^{q+1}/\Delta x|}{1-|\gamma\lambda^{q+1}/\Delta x|}.$$

Hence omission of waves with small eigenvalues compared to the dynamical time scale (i.e. *γλ* � 1) will not significantly affect preconditioner accuracy.

#### **Local Non-Constant Coefficient Correction Solve:**

The remaining component of the split preconditioner (47) comprises the local system *Plocalu* = *v*,

$$\begin{aligned} \iff & \begin{bmatrix} I + \gamma \partial\_{\mathbf{x}} \left( f\_{F} \mathcal{R}\_{F} \right) L\_{F} + \gamma \partial\_{\mathcal{Y}} \left( f\_{G} \mathcal{R}\_{G} \right) L\_{G} + \gamma \partial\_{\mathcal{Z}} \left( f\_{H} \mathcal{R}\_{H} \right) L\_{H} \end{bmatrix} \boldsymbol{u} = \boldsymbol{v} \\ \iff & \begin{bmatrix} I + \gamma \partial\_{\mathbf{x}} \left( \mathcal{R}\_{F} \Lambda\_{F} \right) L\_{F} + \gamma \partial\_{\mathcal{Y}} \left( \mathcal{R}\_{G} \Lambda\_{G} \right) L\_{G} + \gamma \partial\_{\mathcal{Z}} \left( \mathcal{R}\_{H} \Lambda\_{H} \right) L\_{H} \end{bmatrix} \boldsymbol{u} = \boldsymbol{v}. \end{aligned}$$

Note that for spatially homogeneous Jacobians, *∂x*(*RF*Λ*F*) = 0 (similarly for *y* and *z*), so this system reduces to *u* = *v*. Hence this component may optionally be included to correct for spatial inhomogeneity in *JF*, *JG* and *JH*. In keeping with the previous discretization approaches, approximate this system as, e.g.

$$\gamma \partial\_{\mathcal{X}} \left( \mathcal{R}\_F \Lambda\_F \right) L\_F \approx \frac{\gamma}{2 \Delta x} \left( \mathcal{R}\_{F,i+1} \Lambda\_{F,i+1} - \mathcal{R}\_{F,i-1} \Lambda\_{F,i-1} \right) L\_{F,i}.$$

These solves are spatially decoupled (with respect to *u*), resulting in a block-diagonal matrix whose solution requires only *n* × *n* dense linear solves at each spatial location.

They report that auto-differentiation tools can lead to an improvement in the accuracy of the

Implicit Numerical Methods for Magnetohydrodynamics 75

The literature on using nonlinear multigrid for MHD is quite sparse. Here we focus on the recent work by Adams et al. (2010). Multigrid methods are motivated by the observation that a low resolution discretization of an operator can capture modes or components of the error that are expensive to compute directly on a highly resolved discretization. More generally, any poorly locally-determined solution component has the potential to be resolved with coarser representation. This process can be applied recursively with a series of coarse grids, thereby requiring that each grid resolve only the components of the error that it can solve efficiently. These coarse grids have fewer grid points, typically about a factor of two in each dimension, such that the total amount of work in multigrid iterations can be expressed as a geometric sum that converges to a small factor of the work on the finest mesh. These concepts can be applied to problems with particles/atoms or pixels as well as the traditional grid or cell variables considered here. Multigrid provides a basic framework within which particular multigrid

Geometric multigrid is useful because it has the potential to be very efficient especially if the geometric domains of interest are simple enough that explicit coarse grids can be practically constructed even if, for instance, unstructured grids are used. Geometric multigrid not only provides a powerful basis on which to build a specific solution algorithm, but also allows for the straightforward use of nonlinear multigrid, or full approximation scheme (FAS) multigrid (Brandt (1977)) and matrix-free implementations. Given that the MHD problems are nonlinear, FAS multigrid is very efficient in that it solves the nonlinear set of equation directly and obviates the need of an outer (Newton) iteration. This is a critical component in attaining textbook efficiency on nonlinear problems. Figure 1 shows the standard multigrid

residuals and current solutions from the fine grid space *k* to the coarse grid space *k* + 1 (the

the fine grid. Common notation for this multigrid V-cycle is V(*μ*1,*μ*2), where *μ*<sup>1</sup> and *μ*<sup>2</sup> are the number of pre- and post-smoothing steps, respectively. The canonical model problem is the Laplacian, for which point-wise Gauss-Seidel smoothers combined with linear interpolation for the restriction and prolongation operators generate method that reduce error by about an order of magnitude per V(1,1) cycle. This is theoretically optimal in that this rate of residual reduction is independent of mesh size and the amount of work in a V-cycle is given by a geometric sum that converges to about five work units. This so-called textbook efficiency has been observed, if not proven, for multigrid methods in a wide variety of applications (see

A concept used to determine if a point-wise smoothing method exists is *h-ellipticity*. Brandt & Dinar (1979) first introduced h-ellipticity and it is described in Trottenberg et al. (2000). H-ellipticity is the minimum Fourier symbol of the high half (in at least one dimension) of the spectrum of a discrete operator divided by the maximum Fourier symbol

*<sup>k</sup>* are the discrete representation, on the fine grid, of the coarse grid functions),

*<sup>k</sup>*+1, which maps the current solution from the coarse grid to

*<sup>k</sup>* , which maps

FAS V-cycle and uses the smoother *<sup>u</sup>* <sup>←</sup> *<sup>S</sup>*(*A*, *<sup>u</sup>*, *<sup>b</sup>*), the restriction operator *<sup>R</sup>k*+<sup>1</sup>

computed Jacobian compared with a finite difference approach.

**4. Nonlinear multigrid method for MHD**

methods can be developed for particular problems.

Trottenberg et al. (2000) and references therein for details).

rows of *Rk*+<sup>1</sup>

and the prolongation operator *P<sup>k</sup>*

*Pd***: Diffusive MHD Preconditioner:** *Pd* solves the remaining diffusive effects within the implicit system,

$$
\partial\_t \mathbf{U} - \nabla \cdot \mathbf{F}\_{\upsilon} = 0.
$$

Setting *Pd* to be the Jacobian of this operator,

$$\begin{aligned} P\_d &= J\_{\mathcal{v}}(\mathbf{U}) = I - \overline{\gamma} \frac{\partial}{\partial \mathbf{U}} \left( \nabla \cdot \mathbf{F}\_{\mathcal{v}} \right) \\ &= \begin{bmatrix} I & 0 & 0 & 0 \\ 0 & I - \overline{\gamma} D\_{\rho \mathbf{v}} & 0 & 0 \\ 0 & 0 & I - \overline{\gamma} D\_{\mathbf{B}} & 0 \\ -\overline{\gamma} L\_{\rho} & -\overline{\gamma} L\_{\rho \mathbf{v}} & -\overline{\gamma} L\_{\mathbf{B}} & I - \overline{\gamma} D\_{\mathcal{C}} \end{bmatrix} \end{aligned}$$

and then its structure may be exploited for efficient and accurate solution. To solve *Pd y* = *b* for *y* = [*yρ*, *yρ***v**, *y***B**, *ye*] *T*:


.

5. Solve (*<sup>I</sup>* <sup>−</sup> *<sup>γ</sup>De*) *ye* <sup>=</sup> ˜ *be* for *ye*.

Due to their diffusive nature, steps 2, 3 and 5 are solved using a system-based geometric multigrid solver. Step 4 may be approximated through one finite-difference, instead of constructing and multiplying by the individual sub-matrices:

$$L\_{\rho}y\_{\rho} + L\_{\rho \mathbf{v}}y\_{\rho \mathbf{v}} + L\_{\mathbf{B}}y\_{\mathbf{B}} = \frac{1}{\sigma} \left[ \nabla \cdot \mathbf{F}\_{\mathbf{v}}(\mathcal{U} + \sigma \mathcal{W}) - \nabla \cdot \mathbf{F}\_{\mathbf{v}}(\mathbf{U}) \right]\_{\mathbf{t}} + \mathcal{O}(\sigma)\_{\mathbf{t}}$$

where *W* = � *<sup>y</sup>ρ*, *<sup>y</sup>ρ***v**, *<sup>y</sup>***B**, 0�*<sup>T</sup>*

As far as implementation details are concerned, Reynolds et al. employ the SUNDIALS software library (Hindmarsh et al. (2005)). Within SUNDIALS, extensive use is made of the CVODE ordinary differential equations integration package, as well as KINSOL for nonlinear solution of algebraic systems.

#### **3.5 Next steps**

Once the preconditioner is in place, several heuristic ideas may be applied to further decrease computational time. Some of these ideas discussed in Reynolds et al. (2010) are: freezing the Jacobian for a few time steps, freezing the computations of the eigen-values and eigen-vector, preconditioning only the fastest waves, eliminating the local solve etc. Depending on the physical problem under investigation, these heuristic ideas can reduce the computational time significantly. Reynolds et al. (2011) have generalized their approach to tokamak geometry wherein the poloidal plane is discretized using a curvilinear mesh. The form of the equations solved are similar to the ones discussed by Samtaney et al. (2007). The complexity of generating the Jacobian for the Newton-Krylov method makes it an attractive candidate for automatic differentiation tools. This is, in fact, employed by Reynolds & Samtaney (2012) and Reynolds et al. (2011) for implicit solution of the resistive MHD in the tokamak geometry. 16 Will-be-set-by-IN-TECH

*Pd***: Diffusive MHD Preconditioner:** *Pd* solves the remaining diffusive effects within the

*∂t***U** −∇· **F***<sup>v</sup>* = 0.

*<sup>∂</sup>***<sup>U</sup>** (∇ · **F***v*)

and then its structure may be exploited for efficient and accurate solution. To solve *Pd y* = *b*

�

*<sup>σ</sup>* [∇ · **F***v*(*U* + *σW*) −∇· **F***v*(**U**)]*<sup>e</sup>* + *O*(*σ*),

Due to their diffusive nature, steps 2, 3 and 5 are solved using a system-based geometric multigrid solver. Step 4 may be approximated through one finite-difference, instead of

As far as implementation details are concerned, Reynolds et al. employ the SUNDIALS software library (Hindmarsh et al. (2005)). Within SUNDIALS, extensive use is made of the CVODE ordinary differential equations integration package, as well as KINSOL for nonlinear

Once the preconditioner is in place, several heuristic ideas may be applied to further decrease computational time. Some of these ideas discussed in Reynolds et al. (2010) are: freezing the Jacobian for a few time steps, freezing the computations of the eigen-values and eigen-vector, preconditioning only the fastest waves, eliminating the local solve etc. Depending on the physical problem under investigation, these heuristic ideas can reduce the computational time significantly. Reynolds et al. (2011) have generalized their approach to tokamak geometry wherein the poloidal plane is discretized using a curvilinear mesh. The form of the equations solved are similar to the ones discussed by Samtaney et al. (2007). The complexity of generating the Jacobian for the Newton-Krylov method makes it an attractive candidate for automatic differentiation tools. This is, in fact, employed by Reynolds & Samtaney (2012) and Reynolds et al. (2011) for implicit solution of the resistive MHD in the tokamak geometry.

*I* 0 00 0 *I* − *γDρ***<sup>v</sup>** 0 0 0 0 *I* − *γD***<sup>B</sup>** 0 −*γL<sup>ρ</sup>* −*γLρ***<sup>v</sup>** −*γL***<sup>B</sup>** *I* − *γDe*

⎤ ⎥ ⎥ ⎦

implicit system,

for *y* = [*yρ*, *yρ***v**, *y***B**, *ye*]

5. Solve (*<sup>I</sup>* <sup>−</sup> *<sup>γ</sup>De*) *ye* <sup>=</sup> ˜

solution of algebraic systems.

1. Update *y<sup>ρ</sup>* = *b<sup>ρ</sup>*

4. Update ˜

where *W* = �

**3.5 Next steps**

Setting *Pd* to be the Jacobian of this operator,

*T*:

2. Solve (*I* − *γDρ***v**) *yρ***<sup>v</sup>** = *bρ***<sup>v</sup>** for *yρ***<sup>v</sup>** 3. Solve (*I* − *γD***B**) *y***<sup>B</sup>** = *b***<sup>B</sup>** for *y***<sup>B</sup>**

*be* = *be* + *γ* �

*Pd* <sup>=</sup> *Jv*(**U**) = *<sup>I</sup>* <sup>−</sup> *<sup>γ</sup> <sup>∂</sup>*

=

⎡ ⎢ ⎢ ⎣

*L<sup>ρ</sup> y<sup>ρ</sup>* + *Lρ***<sup>v</sup>** *yρ***<sup>v</sup>** + *L***<sup>B</sup>** *y***<sup>B</sup>**

*be* for *ye*.

constructing and multiplying by the individual sub-matrices:

*L<sup>ρ</sup> y<sup>ρ</sup>* + *Lρ***<sup>v</sup>** *yρ***<sup>v</sup>** + *L***<sup>B</sup>** *y***<sup>B</sup>** = <sup>1</sup>

*<sup>y</sup>ρ*, *<sup>y</sup>ρ***v**, *<sup>y</sup>***B**, 0�*T*.

They report that auto-differentiation tools can lead to an improvement in the accuracy of the computed Jacobian compared with a finite difference approach.

### **4. Nonlinear multigrid method for MHD**

The literature on using nonlinear multigrid for MHD is quite sparse. Here we focus on the recent work by Adams et al. (2010). Multigrid methods are motivated by the observation that a low resolution discretization of an operator can capture modes or components of the error that are expensive to compute directly on a highly resolved discretization. More generally, any poorly locally-determined solution component has the potential to be resolved with coarser representation. This process can be applied recursively with a series of coarse grids, thereby requiring that each grid resolve only the components of the error that it can solve efficiently. These coarse grids have fewer grid points, typically about a factor of two in each dimension, such that the total amount of work in multigrid iterations can be expressed as a geometric sum that converges to a small factor of the work on the finest mesh. These concepts can be applied to problems with particles/atoms or pixels as well as the traditional grid or cell variables considered here. Multigrid provides a basic framework within which particular multigrid methods can be developed for particular problems.

Geometric multigrid is useful because it has the potential to be very efficient especially if the geometric domains of interest are simple enough that explicit coarse grids can be practically constructed even if, for instance, unstructured grids are used. Geometric multigrid not only provides a powerful basis on which to build a specific solution algorithm, but also allows for the straightforward use of nonlinear multigrid, or full approximation scheme (FAS) multigrid (Brandt (1977)) and matrix-free implementations. Given that the MHD problems are nonlinear, FAS multigrid is very efficient in that it solves the nonlinear set of equation directly and obviates the need of an outer (Newton) iteration. This is a critical component in attaining textbook efficiency on nonlinear problems. Figure 1 shows the standard multigrid FAS V-cycle and uses the smoother *<sup>u</sup>* <sup>←</sup> *<sup>S</sup>*(*A*, *<sup>u</sup>*, *<sup>b</sup>*), the restriction operator *<sup>R</sup>k*+<sup>1</sup> *<sup>k</sup>* , which maps residuals and current solutions from the fine grid space *k* to the coarse grid space *k* + 1 (the rows of *Rk*+<sup>1</sup> *<sup>k</sup>* are the discrete representation, on the fine grid, of the coarse grid functions), and the prolongation operator *P<sup>k</sup> <sup>k</sup>*+1, which maps the current solution from the coarse grid to the fine grid. Common notation for this multigrid V-cycle is V(*μ*1,*μ*2), where *μ*<sup>1</sup> and *μ*<sup>2</sup> are the number of pre- and post-smoothing steps, respectively. The canonical model problem is the Laplacian, for which point-wise Gauss-Seidel smoothers combined with linear interpolation for the restriction and prolongation operators generate method that reduce error by about an order of magnitude per V(1,1) cycle. This is theoretically optimal in that this rate of residual reduction is independent of mesh size and the amount of work in a V-cycle is given by a geometric sum that converges to about five work units. This so-called textbook efficiency has been observed, if not proven, for multigrid methods in a wide variety of applications (see Trottenberg et al. (2000) and references therein for details).

A concept used to determine if a point-wise smoothing method exists is *h-ellipticity*. Brandt & Dinar (1979) first introduced h-ellipticity and it is described in Trottenberg et al. (2000). H-ellipticity is the minimum Fourier symbol of the high half (in at least one dimension) of the spectrum of a discrete operator divided by the maximum Fourier symbol

**function** *u* ← *MGF*(*Ak*, *uk*, *fk*) **if coarse grid** *k* + 1 **exists** *rk* ← *fk* − *Akuk rk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup>

*wk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup>

*uk* <sup>←</sup> *uk* <sup>+</sup> *<sup>P</sup><sup>k</sup>*

*uk* <sup>←</sup> *MGV*(*A*˜

*uk* <sup>←</sup> *<sup>A</sup>*−<sup>1</sup>

**else**

**return** *uk*

Newton method.

**5. Numerical test cases**

**5.1 Linear wave propagation**

*<sup>b</sup>* = (cos *<sup>α</sup>* cos *<sup>θ</sup>*, sin *<sup>α</sup>* sin *<sup>θ</sup>*, 0)*T*. Here, *<sup>θ</sup>* = tan−<sup>1</sup> *ky*

perturbing the *<sup>k</sup>*−th characteristic, *<sup>w</sup><sup>k</sup>* <sup>=</sup> *<sup>w</sup><sup>k</sup>* <sup>+</sup> *�* cos

easily extended to three dimensions.

*<sup>k</sup>* (*rk*) /\* restriction of residual to coarse grid \*/

Implicit Numerical Methods for Magnetohydrodynamics 77

*uk*<sup>+</sup><sup>1</sup> ← *uk*<sup>+</sup><sup>1</sup> − *wk*<sup>+</sup><sup>1</sup> /\* convert solution to an increment \*/

*<sup>k</sup> fk* /\* accurate solve of coarsest grid \*/

less than the incremental error on the model problem (Trottenberg et al. (2000)). Multigrid can thus achieve discretization error with a work complexity of a few residual calculations. An additional advantage of the FAS multigrid algorithm is that it is an effective global nonlinear solver in that it does not suffer from the problem of limited radius of convergence of a standard

Linear wave propagation refers to the initialization of low amplitude magnetosonic or Alfvén waves and computing their evolution using the nonlinear equations. If the amplitude is small (O(*�*)), these waves will propagate linearly with nonlinear effects essentially being <sup>O</sup>(*�*2). Linear waves may be initialized in 2D using the following procedure. First choose a background quiescent equilibrium state as *U*˜ = (*ρ*, *ρu*, *b*,*e*)*T*, where *ρ* = 1, *u* = **0**,

*kx*

wave propagation and *α* is the orientation of the constant magnetic field. We project these equilibrium conserved quantities to characteristic variables via *W* = *LU*˜ , where *L* is the left eigenvector matrix of the linearized MHD system. The *k*−th linear wave is setup by

then set as *U*(*x*, *y*, 0) = *R*(*U*˜ )*W*, where *R* is the right eigenvector matrix. Periodic boundary conditions should be implemented in both the *x*- and *y*-directions. This procedure can be

Chacón (2004) tested the evolution of a magnetosonic wave to verify that the method had low dissipation using a Newton-Krylov approach but without any preconditioning. Reynolds et al. (2006) also tested the evolution of a slow magnetosonic wave propagating 45 deg to the mesh, with a Newton-Krylov solver without preconditioning. Numerical tests at 2562 mesh resolution in 2D indicated that even without preconditioning, the implicit NK method yielded over a factor of ten decrease in CPU time. Reynolds et al. (2010) reported a further benefit of over a factor of five decrease in CPU time when the wave-structure based preconditioner was employed for the linear wave propagation test. Furthermore, for linear waves aligned with the mesh, the preconditioned solves converged in one Krylov iteration

, in which the ratio *ky*

*πkxx* + *πkyy*

*kx* gives the direction of

. The initial condition is

Fig. 2. Nonlinear FAS multigrid *F-cycle* algorithm with defect correction

*<sup>k</sup>* (*uk*) /\* restriction of residual to coarse grid \*/ *uk*<sup>+</sup><sup>1</sup> ← *MGF*(*Ak*+1, *wk*+1,*rk*<sup>+</sup><sup>1</sup> + *Ak*<sup>+</sup>1*wk*<sup>+</sup>1) /\* recursive multigrid application \*/

*<sup>k</sup>*+1(*uk*<sup>+</sup>1) /\* prolongation of coarse grid correction \*/

*<sup>k</sup>*, *uk*, *fk* <sup>−</sup> *Akuk* <sup>+</sup> *<sup>A</sup>*˜ *kuk*) /\* low-order V-cycle, defect correction \*/

**function** *u* ← *MGV*(*Ak*, *uk*, *fk*) **if coarse grid** *k* + 1 **exists** *uk* ← *S*(*Ak*, *uk*, *fk*) /\* pre-smoothing \*/ *rk* ← *fk* − *Akuk rk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup> *<sup>k</sup>* (*rk*) /\* restriction of residual to coarse grid \*/ *wk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup> *<sup>k</sup>* (*uk*) /\* restriction of current solution to coarse grid \*/ *uk*<sup>+</sup><sup>1</sup> ← *MGV*(*Ak*+1, *wk*+1,*rk*<sup>+</sup><sup>1</sup> + *Ak*<sup>+</sup>1*wk*<sup>+</sup>1) /\* recursive multigrid application \*/ *uk*<sup>+</sup><sup>1</sup> ← *uk*<sup>+</sup><sup>1</sup> − *wk*<sup>+</sup><sup>1</sup> /\* convert solution to an increment \*/ *uk* <sup>←</sup> *uk* <sup>+</sup> *<sup>P</sup><sup>k</sup> <sup>k</sup>*+1(*uk*<sup>+</sup>1) /\* prolongation of coarse grid correction \*/ *uk* ← *S*(*Ak*, *uk*, *fk*) /\* restriction of residual to coarse grid \*/ **else** *uk* <sup>←</sup> *<sup>A</sup>*−<sup>1</sup> *<sup>k</sup> fk* /\* post-smoothing \*/ **return** *uk*

Fig. 1. Nonlinear FAS multigrid *V-cycle* algorithm

of the operator. An h-ellipicity bounded well above zero is a necessary and sufficient condition for the existence of a point-wise smoother for an operator with a symmetric stencil (Trottenberg et al. (2000)). An important result of h-ellipticity is that effective point-wise smoothers (eg, Gauss-Seidel and distributive Gauss-Seidel) can be constructed for upwind discretizations of hyperbolic systems with no restriction on the time step, whereas point-wise Gauss-Seidel is unstable for a central difference scheme for a large time step. Adams et al. observed textbook multigrid efficiency with standard multigrid methods (e.g., point-wise Gauss-Seidel smoothers) using a first-order upwinding method for ideal and resistive MHD. First-order accuracy is, however, generally not sufficient for many applications, and second-order schemes are required for efficiency. These stable low-order smoothers have been used extensively with a higher-order operator via a defect correction scheme, which is identical to preconditioning, but is more amenable to a nonlinear solve (Atlas & Burrage (1994); Böhmer et al. (1984); Dick (1991); Hemker (1986); Koren (1991)).

An additional requirement of an optimal solver is to be able to reduce the algebraic error to the order of the discretization (or truncation) error for steady-state problems. For transient problems the solver needs to reduce the algebraic error to below the *incremental* error – that is, the product of the truncation error of the time integration scheme and the spacial truncation error. Reducing algebraic error far below that of the incremental error is computationally wasteful, though potentially useful for debugging. There is generally no need to spend resources to reduce the algebraic error far below the incremental error. This observation leads to our definition of an optimal solver as one that can reduce the error to less than the incremental error with a few work units per time step. This is an ambitious goal in that it requires both scalability and small constants in the actual computational costs. In fact, this results in a solver in which the rate of reduction in the residual actually increases as the mesh is refined, because the truncation error decreases. This goal can be achieved by using a multigrid V-cycle within what is called an F-cycle iteration (Trottenberg et al. (2000)). Figure 2 shows the standard nonlinear multigrid F-cycle with defect correction to accommodate the nonlinear V-cycle with a lower-order operator (*A*˜ is the first-order upwinding operator) for which our point-wise Gauss-Seidel smoother is stable. The complexity of an F-cycle is asymptotically similar to a V-cycle, and it can be proven to result in a solution with algebraic error that is 18 Will-be-set-by-IN-TECH

*<sup>k</sup>* (*rk*) /\* restriction of residual to coarse grid \*/

*uk*<sup>+</sup><sup>1</sup> ← *MGV*(*Ak*+1, *wk*+1,*rk*<sup>+</sup><sup>1</sup> + *Ak*<sup>+</sup>1*wk*<sup>+</sup>1) /\* recursive multigrid application \*/

*<sup>k</sup>*+1(*uk*<sup>+</sup>1) /\* prolongation of coarse grid correction \*/

*uk*<sup>+</sup><sup>1</sup> ← *uk*<sup>+</sup><sup>1</sup> − *wk*<sup>+</sup><sup>1</sup> /\* convert solution to an increment \*/

*<sup>k</sup> fk* /\* post-smoothing \*/

(1994); Böhmer et al. (1984); Dick (1991); Hemker (1986); Koren (1991)).

*uk* ← *S*(*Ak*, *uk*, *fk*) /\* restriction of residual to coarse grid \*/

of the operator. An h-ellipicity bounded well above zero is a necessary and sufficient condition for the existence of a point-wise smoother for an operator with a symmetric stencil (Trottenberg et al. (2000)). An important result of h-ellipticity is that effective point-wise smoothers (eg, Gauss-Seidel and distributive Gauss-Seidel) can be constructed for upwind discretizations of hyperbolic systems with no restriction on the time step, whereas point-wise Gauss-Seidel is unstable for a central difference scheme for a large time step. Adams et al. observed textbook multigrid efficiency with standard multigrid methods (e.g., point-wise Gauss-Seidel smoothers) using a first-order upwinding method for ideal and resistive MHD. First-order accuracy is, however, generally not sufficient for many applications, and second-order schemes are required for efficiency. These stable low-order smoothers have been used extensively with a higher-order operator via a defect correction scheme, which is identical to preconditioning, but is more amenable to a nonlinear solve (Atlas & Burrage

An additional requirement of an optimal solver is to be able to reduce the algebraic error to the order of the discretization (or truncation) error for steady-state problems. For transient problems the solver needs to reduce the algebraic error to below the *incremental* error – that is, the product of the truncation error of the time integration scheme and the spacial truncation error. Reducing algebraic error far below that of the incremental error is computationally wasteful, though potentially useful for debugging. There is generally no need to spend resources to reduce the algebraic error far below the incremental error. This observation leads to our definition of an optimal solver as one that can reduce the error to less than the incremental error with a few work units per time step. This is an ambitious goal in that it requires both scalability and small constants in the actual computational costs. In fact, this results in a solver in which the rate of reduction in the residual actually increases as the mesh is refined, because the truncation error decreases. This goal can be achieved by using a multigrid V-cycle within what is called an F-cycle iteration (Trottenberg et al. (2000)). Figure 2 shows the standard nonlinear multigrid F-cycle with defect correction to accommodate the nonlinear V-cycle with a lower-order operator (*A*˜ is the first-order upwinding operator) for which our point-wise Gauss-Seidel smoother is stable. The complexity of an F-cycle is asymptotically similar to a V-cycle, and it can be proven to result in a solution with algebraic error that is

*<sup>k</sup>* (*uk*) /\* restriction of current solution to coarse grid \*/

*uk* ← *S*(*Ak*, *uk*, *fk*) /\* pre-smoothing \*/

**function** *u* ← *MGV*(*Ak*, *uk*, *fk*) **if coarse grid** *k* + 1 **exists**

> *rk* ← *fk* − *Akuk rk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup>

*wk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup>

*uk* <sup>←</sup> *uk* <sup>+</sup> *<sup>P</sup><sup>k</sup>*

*uk* <sup>←</sup> *<sup>A</sup>*−<sup>1</sup>

Fig. 1. Nonlinear FAS multigrid *V-cycle* algorithm

**else**

**return** *uk*

**function** *u* ← *MGF*(*Ak*, *uk*, *fk*) **if coarse grid** *k* + 1 **exists** *rk* ← *fk* − *Akuk rk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup> *<sup>k</sup>* (*rk*) /\* restriction of residual to coarse grid \*/ *wk*<sup>+</sup><sup>1</sup> <sup>←</sup> *<sup>R</sup>k*+<sup>1</sup> *<sup>k</sup>* (*uk*) /\* restriction of residual to coarse grid \*/ *uk*<sup>+</sup><sup>1</sup> ← *MGF*(*Ak*+1, *wk*+1,*rk*<sup>+</sup><sup>1</sup> + *Ak*<sup>+</sup>1*wk*<sup>+</sup>1) /\* recursive multigrid application \*/ *uk*<sup>+</sup><sup>1</sup> ← *uk*<sup>+</sup><sup>1</sup> − *wk*<sup>+</sup><sup>1</sup> /\* convert solution to an increment \*/ *uk* <sup>←</sup> *uk* <sup>+</sup> *<sup>P</sup><sup>k</sup> <sup>k</sup>*+1(*uk*<sup>+</sup>1) /\* prolongation of coarse grid correction \*/ *uk* <sup>←</sup> *MGV*(*A*˜ *<sup>k</sup>*, *uk*, *fk* <sup>−</sup> *Akuk* <sup>+</sup> *<sup>A</sup>*˜ *kuk*) /\* low-order V-cycle, defect correction \*/ **else** *uk* <sup>←</sup> *<sup>A</sup>*−<sup>1</sup> *<sup>k</sup> fk* /\* accurate solve of coarsest grid \*/ **return** *uk*

Fig. 2. Nonlinear FAS multigrid *F-cycle* algorithm with defect correction

less than the incremental error on the model problem (Trottenberg et al. (2000)). Multigrid can thus achieve discretization error with a work complexity of a few residual calculations. An additional advantage of the FAS multigrid algorithm is that it is an effective global nonlinear solver in that it does not suffer from the problem of limited radius of convergence of a standard Newton method.
