0 Implicit Numerical Methods for Magnetohydrodynamics

A fluid description of the plasma is obtained by taking velocity moments of the kinetic equations (Vlasov or Fokker-Planck equations) for electrons and ions and employing certain closure assumptions. A hierarchy of MHD models can be derived. Generally, if the time scales of interest are larger than the electron-ion collision time scales, then one may model the plasma as a single fluid. Furthermore, the fluid description of a plasma is valid when the length scales under investigation are larger than the Debye length; and the frequencies are smaller than the cyclotron frequency. The Debye length argument can also be cast in terms of a frequency: namely the plasma frequency. In addition, it is a standard assumption that the speeds involved are much smaller than the speed of light. The oft-used term “resistive MHD” is a single-fluid model of a plasma in which a single velocity and pressure describe both the electrons and ions. The resistive MHD model of a magnetized plasma does not include finite Larmor radius (FLR) effects, and is based on the simplifying limit in which the particle collision length is small compared with the macroscopic length scales.


Introduction
A fluid description of the plasma is obtained by taking velocity moments of the kinetic equations (Vlasov or Fokker-Planck equations) for electrons and ions and employing certain closure assumptions.A hierarchy of MHD models can be derived.Generally, if the time scales of interest are larger than the electron-ion collision time scales, then one may model the plasma as a single fluid.Furthermore, the fluid description of a plasma is valid when the length scales under investigation are larger than the Debye length; and the frequencies are smaller than the cyclotron frequency.The Debye length argument can also be cast in terms of a frequency: namely the plasma frequency.In addition, it is a standard assumption that the speeds involved are much smaller than the speed of light.The oft-used term "resistive MHD" is a single-fluid model of a plasma in which a single velocity and pressure describe both the electrons and ions.The resistive MHD model of a magnetized plasma does not include finite Larmor radius (FLR) effects, and is based on the simplifying limit in which the particle collision length is small compared with the macroscopic length scales.

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 (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 κ,a n dt h e 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,w h i c hi st h e ratio of momentum to thermal diffusivity.The non-dimensionalization is carried out using a characteristic length scale, L, and the Alfvén speed U 0 = B 0 / √ μ 0 ρ 0 ,whereB 0 , ρ 0 ,andμ 0 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 and the stress-strain tensor relation 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).

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

61
Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.com 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. (3) This can be rewritten as We then subtract a term from either side of the above equation as where a 0 is a constant coefficient chosen mainly from stability considerations.Furthermore, a 0 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 0 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 0 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: 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

62
Topics in Magnetohydrodynamics www.intechopen.com The above equation is implicit and is solved iteratively.Let U n+1,k denote the k-th iteration of the solution at the n + 1-th time level.Rewrite the above equation as where ∂U ∂t A truncated Taylor series expansion yields: where . 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.Ri,j ≡ R(U i,j , U i+ 1 2 ,j , ). Substituting all back into equation (8) gives The above equation is linear and iterated until ∆U k 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.

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.

General approach
The entire ideal MHD (or resistive MHD and beyond) can be written as a nonlinear function as follows: where U n+1 is the vector of unknowns at time step n + 1.For example, if we use a θ-scheme, one can write the nonlinear function as: 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 8MN.T h ea b o v e nonlinear systems can be solved using an inexact Newton-Krylov solver.Apply the standard Newton's method to the above nonlinear system gives where is the Jacobian; and δU k ≡ U n+1,k+1 − U n+1,k ,a n dk is the iteration index in the Newton method.For the two dimensional system the Jacobian matrix is 8MN × 8MN 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 where r 0 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 64 Topics in Magnetohydrodynamics www.intechopen.commost 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.
It the approximate solution U n+1,k is "close" to the true solution U * of the nonlinear system, the convergence is quadratic, i.e., where C is a constant independent of U n+1,k 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, If we impose the condition that lim k→∞ η k = 0 then convergence is super-linear, and if η k is constant for all k 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 η k are recommended which take into account how well the nonlinear system is converging.
where γ 1 = 0.9 and γ 2 = 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,

65
Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.com 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: 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).

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 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 The Krylov iteration is then used to solve one of where X = −P −1 f is computed prior to the Krylov solve or V = P −1 W is computed after the Krylov solve.Scalable convergence of the method then depends on the spectrum of the preconditioned operator (JP −1 or P −1 J), as opposed to the original Jacobian operator J.
Hence, an optimal preconditioning strategy will satisfy the two competing criteria: 1. P ≈ J, to help cluster the spectrum of the preconditioned operator.
2. Application of P −1 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 increasing number of processors in a parallel simulation.

66
Topics in Magnetohydrodynamics www.intechopen.com 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: • 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 and storing the Jacobian matrix.• "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 themselves and in this sense this type of preconditioning is "physics-based".

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.

A model illustration
Consider the following hyperbolic system Differencing with backward Euler we get Substitute the second equation into the first to obtain: 67 Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.comwhich 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.
Stiff off-diagonal blocks L and U are now shifted to the diagonal via the Schur complement

Application to resistive MHD
We begin by examining the linearized resistive MHD equations.These are written as which illustrates the couplings between the various unknowns.In NK the Jacobian has the following coupling 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 where

68
Topics in Magnetohydrodynamics www.intechopen.com 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: where P S = D v − LM −1 U is the Schur complement.The exact Jacobian inverse require M −1 and P −1 S .The following predictor-corrector algorithm is proposed.
Multigrid is impractical for P S because of the M −1 factor and hence some simplifications are desirable.For the velocity update and the corrector part in the above equations, we can treat where P SI = D v − ∆tLU and is block-diagonally dominant.Multigrid is employed to compute the inverse of P SI and M.

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): 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 n solves the nonlinear residual equation g(U)= 0. q n determines the method's order of accuracy and at q n = {1, 2} the method is stable for any ∆t n , with stability decreases as q n increases.α n,i and β n,i are fixed parameters for a given method order q n .In this approach ∆t n , q n are adaptively chosen at each time step to balance solution accuracy, solver convergence, and temporal stability (Hindmarsh (2000)).

69
Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.com Alternatively one may also use a θ-scheme 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 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 fi+ 1 2 ,j,k is referred to as the numerical flux through the face {i + 1 2 , j, k} and is computed as a linear combination of the fluxes at cell centers as Reynolds et al. give the options for several spatial difference schemes.For a second-order central difference implementation, m = 0, n = 1a n da 0 = a 1 = 1 2 ; for a fourth-order central difference approximation, m = 1, n = 2, and a −1 = a 2 = −1 12 , a 0 = a 1 = 7 12 ;a n d for tuned second-order central differences, a −1 = a 2 = −0.197, a 0 = a 1 = 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 where B α is the α-component of the magnetic field, and the terms Bα 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.

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,wesolve the related system (JP −1 )(P δU)=−g, i.e., the right preconditioning approach is adopted.
Since MHD stiffness results from fast hyperbolic and diffusive effects, we set 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.
P h : 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.
This nonlinear problem has Jacobian with, e.g., J F (U)= ∂ ∂U F(U).Using the notation (•) to denote the location at which the action of the linear operator takes place, e.g.
Omitting the explicit dependence on U from the notation, and introducing nonsingular matrices L F , L G and L H , re-write the Jacobian system (46) as 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, (47)

71
Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.com Denote these components as P = P F P G P H P local .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 J F , J G and J H .H en c e,Pu = b may be solved via the steps (i) P F χ = b, (ii) P G w = χ, (iii) P H v = w,and(iv)P local 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 P F , P G and P H from (47) of the form, e.g.(x-direction) To this point L F , L G ,a n dL H are still unspecified.These are n × n matrices (n = 7o r8f o r 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, where R F ≡ L −1 F are the right eigenvectors (n × n column matrix), and λ k are the eigenvalues of J F .Through pre-multiplication of (48) by L F ,gives Defining the vector of characteristic variables w = L F χ, decouple the equations as , where w k denotes the k-th element of the characteristic vector w,andβ = L F b.
Spatial discretization of each of the characteristic variables w k 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 k j .F o re x a m p l et h e tridiagonal version due to a O(∆x 2 ) finite-difference discretization is Reynolds et al. use a second order centered finite-volume approximation, with resulting systems for each w k that are tridiagonal.Moreover, the above approach results not only in tridiagonal systems for each characteristic variable w k , 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 characteristic variable w k .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 ΛF = diag(λ 1 ,...,λ q ,0,...,0).Omission of the (n − q) slowest waves in this fashion amounts to a further approximation of the preconditioner to the original discretized PDE system.Writing PF as the x-directional preconditioner based on q waves, consider χ − χ p ,whereχ solves P F χ = b and χ solves PF χ = b,i.e.
where ĴF = R F ΛF L F .Left-multiplying by L F and proceeding as before, to obtain Since the eigenvector matrices L F and R F may be renormalized as desired, and the eigenvalues are ordered so that λ i ≥ λ j ,f o ri < j, the dominant error from preconditioning only the q fastest waves is approximately |γλ q+1 /∆x| 1 −|γλ q+1 /∆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 P local u = v, Note that for spatially homogeneous Jacobians, ∂ x (R F Λ 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 J F , J G and J H .In keeping with the previous discretization approaches, approximate this system as, e.g.
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.

73
Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.comP d : Diffusive MHD Preconditioner: P d solves the remaining diffusive effects within the implicit system, Setting P d to be the Jacobian of this operator, and then its structure may be exploited for efficient and accurate solution.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: where W = y ρ , y ρv , y B ,0 T .
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.

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.
They report that auto-differentiation tools can lead to an improvement in the accuracy of the computed Jacobian compared with a finite difference approach.

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 u ← S(A, u, b), the restriction operator R k+1 k ,whichmaps residuals and current solutions from the fine grid space k to the coarse grid space k + 1( t h e rows of R k+1 k are the discrete representation, on the fine grid, of the coarse grid functions), and the prolongation operator P k k+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 μ 1 and μ 2 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 75 Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.com 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 ( Ã 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 76 Topics in Magnetohydrodynamics www.intechopen.com /* restriction of residual to coarse grid */ w k+1 ← R k+1 k (u k ) /* restriction of residual to coarse grid */ u k+1 ← MGF(A k+1 , w k+1 , r k+1 + A k+1 w k+1 ) /* recursive multigrid application */ u k+1 ← u k+1 − w k+1 /* convert solution to an increment */ /* accurate solve of coarsest grid */ return u k 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.

Linear wave propagation
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 O(ǫ 2 ).Linear waves may be initialized in 2D using the following procedure.First choose a background quiescent equilibrium state as Ũ =( ρ, ρu, b, e) T ,w h e r eρ = 1, u = 0, b =(cos α cos θ,sinα sin θ,0) T .Here, θ = tan −1 k y k x , in which the ratio k y k x gives the direction of wave propagation and α is the orientation of the constant magnetic field.We project these equilibrium conserved quantities to characteristic variables via W = L Ũ,w h e r eL is the left eigenvector matrix of the linearized MHD system.The k−th linear wave is setup by perturbing the k−th characteristic, w k = w k + ǫ cos πk x x + πk y y .The initial condition is then set as U(x, y,0)=R( Ũ)W,w h er eR is the right eigenvector matrix.Periodic boundary conditions should be implemented in both the x-a n dy-directions.This procedure can be easily extended to three dimensions.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 256 2 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  for several tests, indicating that the wave-structure based preconditioner is optimal for such cases.
Here we reproduce results from Reynolds et al. (2010) of a slow magnetosonic wave of amplitude ǫ = 10 −5 in a periodic domain chosen as [0, 2] × [0, 2].T h e w a v e is propagated until a final time of 10 units.The equilibrium state chosen in Ũ = (1, 0, cos α cos θ,sinα sin θ,0,0.1)T , α = −44.5 o .Two different propagation directions are: θ = 0, 45 o , i.e., the wave propagates aligned with the x− axis, and along the diagonal.Results for the wave propagation are shown in Figure 3.The total number of Krylov iterations is plotted for different time step sizes and spatial discretizations (horizontal axis).For the linear wave propagating aligned with the mesh, the preconditioner is nearly exact, and hence the Krylov iterations remain nearly constant as the mesh is refined, as compared with the non-preconditioned tests that increase rapidly.For the oblique propagation case, the directional splitting does not appear to significantly affect the preconditioner accuracy, again resulting in nearly constant Krylov iterations with mesh refinement.

Magnetic reconnection in 2D
Magnetic reconnection (MR) refers to the breaking and reconnecting of oppositely-directed magnetic field lines in a plasma.In this process, magnetic field energy is converted to plasma kinetic and thermal energy.A test which has gained a lot of popularity in testing MHD codes is the so-called GEM reconnection challenge problem (Birn & et al. (2001a)).The initial conditions consist of a perturbed Harris sheet configuration as described in Birn & et al. (2001a).Reynolds et al. (2006) computed the GEM reconnection challenge problem with a Newton-Krylov method without preconditioning and reproduced the expected Sweet-Parker scaling for the reconnection rate for Lundquist numbers ranging from S = 200 − 10 4 .Furthermore, for a mesh resolution of 512 × 256 their implicit method (without preconditioning) achieved a speedup of about 5.6 compared with an explicit method.
The GEM reconnection problem was also chosen for extensive testing by the nonlinear multigrid method developed by Adams et al. (2010).In fact, this work also extended the GEM problem by including a guide field in the third direction, thereby increasing the stiffness induced by more than a factor of five.Adams et al. also reported on the scalability of 78 Topics in Magnetohydrodynamics www.intechopen.comtheir approach by demonstrating good weak scaling up to 32K processors on a CRAY XT-5 supercomputer.a few residual calculations (work units) per time step, with the largest time step that can accurately resolve the dynamics of the problem.In this study, the solver is fixed at one iteration of FAS F-cycle with two defect corrected V(1,1) cycles at each level, as described in Section §4, and with a work complexity of about 18 work units per time step.There are three applications of the fine grid operator in residual calculations and defect correction in FAS multigrid, and three fine grid work units in the smoothers and residual calculations in each of the two V(1,1) cycle, plus lower-order work in restriction/prolongation and FAS terms.This results in about ten work units on the fine grid.Each successive grid is four times smaller (in 2D), and F-cycles process the second grid twice, the third grid three times, and so on, resulting in the equivalent of about eight additional work units for a total of 18 work units (there are actually fewer total work units in 3D because the coarse grids are relatively smaller).
The smoother is nonlinear Gauss-Seidel with one iteration per grid point and red-black (or checkerboard) ordering.Even though Adams et al. (2010) use defect correction, they demonstrate a second-order rate of convergence on several important diagnostic quantities: these are the kinetic energy, reconnection rate and reconnected magnetic flux as shown in Figure 5,

Ideal Kelvin-Helmholtz instability
This test is generally a hard test for implicit solvers because the growth rate of the instability is high and the dynamics becomes nonlinear very quickly stressing all aspects of an implicit solver.On a 256 2 mesh, Chacón (2008b) reports a speedup in excess of three for the implicit solver compared with an explicit one, and nearly ten Krylov iterations per time step (and nearly five Newton steps per time step) for a time step which was 156 times larger than that for an explicit method.Reynolds et al. (2010) also performed the ideal Kelvin-Helmholtz instability (KHI) test with their wave preconditioner NK solver, and reported a speed up over a factor of three compared with simulations without using a preconditioner for a 256 2 mesh.They did not report comparisons with an explicit time stepping solver.In their 2D simulations, the number of Krylov iterations per time step ranged from 6-13 and the number of Newton steps ranged from 1-3 per time step.We hasten to add that this is not meant to be a comparison to 3, with the associated preconditioned Krylov iteration counts in the range of 6-13 per time step.Solver results for these tests are shown in Figure 7.For all time step sizes and all spatial discretizations the preconditioner results in significantly fewer linear iterations, with the disparity growing as the mesh is refined.

Other examples
There are a variety of other test cases reported in the literature ranging from ideal to resistive MHD.Chacón (Chacón (2004;2008a;b)) reports on 2D tearing instability test cases and 81 Implicit Numerical Methods for Magnetohydrodynamics www.intechopen.comdemonstrates a speedup ranging from 8 − 15 for a 128 2 mesh Chacón (2008a).Another example of a good verification test case in 3D is that of 3D island coalescence (Chacón (2008a)).Reynolds et al. (2006) reported on a 3D ideal MHD problem which models pellet fueling in tokamaks.

Conclusion
In this chapter, we discussed the need for implicit algorithms for resistive magnetohydrodynamics.
We highlighted two broad classes of nonlinear methods: Newton-Krylov and nonlinear multigrid.We illustrated two Newton-Krylov approaches for MHD which are essentially very similar in the overall approach, but differed in the preconditioning strategies for expediting the iterative solution steps in the Krylov linear solver stage of the overall method.One preconditioning strategy is based on a "parabolization" approach while the other utilizes the local wave structure of the underlying hyperbolic waves in the MHD PDEs.The literature on the use of nonlinear multigrid for MHD is essentially sparse and therein we focused on a defect-correction approach coupled with a point-wise Gauss-Seidel smoother utilizing a first order upwind approach.Both approaches are valid and have their place, but it is clear that the nonlinear multigrid approach for MHD is still relatively new and could be further developed.

Future challenges
In this chapter, we have focused exclusively on methods for single fluid resistive MHD.Future challenges will lie in the area of implicit methods for more complicated extended MHD models with FLR effects, several of which exhibit dispersive wave phenomena such as Whistler, Kinetic Alfvén waves, and gyroviscous waves.These dispersive high frequency waves essentially make the stable explicit time step proportional to the square of the mesh spacing, i.e., ∆t ∝∆x 2 ; and hence the benefit from implicit methods is much more than those for single fluid MHD.Some progress in using Newton Krylov approaches for Hall-MHD has been reported by Chacón & Knoll (2003).However more work is required for general geometry, and inclusion of all dispersive wave families.Research in the area of nonlinear multigrid is essentially unexplored for extended MHD.Another interesting challenge in developing implicit methods for MHD is the combination of JFNK or FAS methods with adaptive mesh refinement (AMR).Some progress towards JFNK with AMR has been reported by Philip et al. (2008) on reduced incompressible MHD in 2D.Combining implicit methods with AMR will help mitigate not only the temporal stiffness issues but also help effectively resolve the range of spatial scales in MHD.
U) ≡∇•F(U) −∇•F d (U), and the solution vector U ≡ U(x, t) is U = {ρ, ρu, B, e} T and the hyperbolic flux vector F(U) and the diffusive fluxes F d (U) are given by To solve P d y = b for y =[y ρ , y ρv , y B , y e ] T : 1. Update y ρ = b ρ 2. Solve (I − γD ρv ) y ρv = b ρv for y ρv 3. Solve (I − γD B ) y B = b B for y B 4. Update be = b e + γ L ρ y ρ + L ρv y ρv + L B y B 5. Solve (I − γD e ) y e = be for y e .

Figure 4 (Fig. 4 .
Figure 4 (reproduced from Adams et al. (2010)) from shows a time sequence of current density J z field during reconnection.The goal is to develop solvers with a complexity equivalent to