**5. Numerical test cases**

### **5.1 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 <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**, *<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 kx* , in which the ratio *ky kx* 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* = *LU*˜ , where *L* is the left eigenvector matrix of the linearized MHD system. The *k*−th linear wave is setup by perturbing the *<sup>k</sup>*−th characteristic, *<sup>w</sup><sup>k</sup>* <sup>=</sup> *<sup>w</sup><sup>k</sup>* <sup>+</sup> *�* cos *πkxx* + *πkyy* . The initial condition is 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 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 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

their approach by demonstrating good weak scaling up to 32K processors on a CRAY XT-5

Implicit Numerical Methods for Magnetohydrodynamics 79

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

(a) (b) (c)

Fig. 4. Time sequence of current density,*Jz* during reconnection at time (a)t=0, (b)t=15 and (c)t=60. Parameters for this test are in Adams et al. (2010). Figure obtained from authors of

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

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 2562 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 2562 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

supercomputer.

Reference (Adams et al. (2010)).

**5.3 Ideal Kelvin-Helmholtz instability**

Figure 5,

Fig. 3. Krylov iterations for the linear wave tests: *x*-directional (left) and oblique (right). Figure obtained from authors of Reference (Reynolds et al. (2010)).

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 *�* <sup>=</sup> <sup>10</sup>−<sup>5</sup> in a periodic domain chosen as [0, 2] <sup>×</sup> [0, 2]. The wave is propagated until a final time of 10 units. The equilibrium state chosen in *U*˜ = (1, 0, cos *<sup>α</sup>* cos *<sup>θ</sup>*, sin *<sup>α</sup>* sin *<sup>θ</sup>*, 0, 0.1)*T*, *<sup>α</sup>* <sup>=</sup> <sup>−</sup>44.5*o*. Two different propagation directions are: *<sup>θ</sup>* <sup>=</sup> 0, 45*o*, i.e., the wave propagates aligned with the *<sup>x</sup>*<sup>−</sup> 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.

#### **5.2 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 *<sup>S</sup>* <sup>=</sup> <sup>200</sup> <sup>−</sup> 104. Furthermore, for a mesh resolution of 512 <sup>×</sup> 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 20 Will-be-set-by-IN-TECH

Fig. 3. Krylov iterations for the linear wave tests: *x*-directional (left) and oblique (right).

for several tests, indicating that the wave-structure based preconditioner is optimal for such

Here we reproduce results from Reynolds et al. (2010) of a slow magnetosonic wave of amplitude *�* <sup>=</sup> <sup>10</sup>−<sup>5</sup> in a periodic domain chosen as [0, 2] <sup>×</sup> [0, 2]. The wave is propagated until a final time of 10 units. The equilibrium state chosen in *U*˜ = (1, 0, cos *<sup>α</sup>* cos *<sup>θ</sup>*, sin *<sup>α</sup>* sin *<sup>θ</sup>*, 0, 0.1)*T*, *<sup>α</sup>* <sup>=</sup> <sup>−</sup>44.5*o*. Two different propagation directions are: *<sup>θ</sup>* <sup>=</sup> 0, 45*o*, i.e., the wave propagates aligned with the *<sup>x</sup>*<sup>−</sup> 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

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 *<sup>S</sup>* <sup>=</sup> <sup>200</sup> <sup>−</sup> 104. Furthermore, for a mesh resolution of 512 <sup>×</sup> 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

Figure obtained from authors of Reference (Reynolds et al. (2010)).

resulting in nearly constant Krylov iterations with mesh refinement.

500

1000

1500

Krylov iterations

2000

2500

No Prec, Δt=5e−2 FW Prec, Δt=5e−2 No Prec, Δt=1e−1 FW Prec, Δt=1e−1 No Prec, Δt=2e−1 FW Prec, Δt=2e−1

64^2 128^2 256^2 <sup>0</sup>

Total Krylov Iteration (oblique propagation)

mesh size

64^2 128^2 256^2 <sup>0</sup>

Total Krylov Iteration (x−directional propagation)

mesh size

**5.2 Magnetic reconnection in 2D**

500

1000

Krylov iterations

cases.

1500

No Prec, Δt=5e−2 FW Prec, Δt=5e−2 No Prec, Δt=1e−1 FW Prec, Δt=1e−1 No Prec, Δt=2e−1 FW Prec, Δt=2e−1 their approach by demonstrating good weak scaling up to 32K processors on a CRAY XT-5 supercomputer.

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

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,

#### **5.3 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 2562 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 2562 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

Fig. 6. Snapshots of *Bx* (left) and *Bz* (right) in the 2D Kelvin–Helmholtz test at *t* = 2. Figure

Implicit Numerical Methods for Magnetohydrodynamics 81

Total Krylov Iterations (2D Kelvin−Helmholtz)

64^2 128^2 256^2

mesh size

Fig. 7. Krylov iterations for the 2D Kelvin–Helmholtz tests. Figure obtained from authors of

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 used, the preconditioner results in significantly fewer linear iterations, with

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

obtained from authors of Reference (Reynolds et al. (2010)).

No Prec, Δt=2.5e−3 FW Prec, Δt=2.5e−3 No Prec, Δt=5e−3 FW Prec, Δt=5e−3 No Prec, Δt=1e−2 FW Prec, Δt=1e−2

0

**5.4 Other examples**

Reference (Reynolds et al. (2010)).

the disparity growing as the mesh is refined.

0.5

1

1.5

Krylov iterations

2

2.5

<sup>3</sup> x 104

Fig. 5. Order of spatial accuracy for simulating magnetic reconnection usinga nonlinear multigrid approach. Error in peak kinetic energy and kinetic energy, reconnection flux rate and reconnection rate, high viscosity cases, *Bz* = 0 (left) and low viscosity *Bz* = 5 (right). Figure obtained from authors of Reference (Adams et al. (2010)).

between KHI simulations by Chacón and Reynolds et al. because their respective setup and solver tolerances were not necessarily identical. However, these test results are an indication of the type of speed up and code performance one may expect with strongly nonlinear MHD cases with a Newton-Krylov approach. Reynolds et al. (2010) also performed simulation tests on the 3D version of the ideal KHI.

Here we reproduce results from Reynolds et al. (2010) for the 2D KHI test. We set the computational domain to −5 4 , 5 4 × −1 2 , 1 2 × −5 4 , 5 4 , with periodic boundary conditions in the *x*- and *z*-directions, and homogeneous Neumann boundary conditions in the *y*-direction. We initialize the constant fields *ρ* = 1, *b* = (0.1, 0, 10)*T*, *p* = 0.25, and *uy* = *uz* = *By* = 0. We then set *ux* = <sup>1</sup> <sup>2</sup> tanh(100*y*) + <sup>1</sup> <sup>10</sup> cos(0.8*πx*) + <sup>1</sup> <sup>10</sup> sin(3*πy*) + <sup>1</sup> <sup>10</sup> cos(0.8*πz*). This problem employs the resistive MHD equations, with resistivity, viscosity, and heat conduction coefficients set to 10−4, and all runs are taken to a final time of *Tf* = 2. As previous results on this problem suggest that the instability growth rate is independent of the size of the resistivity, such small parameters are natural since the instability is predominantly driven by nonlinear (hyperbolic) effects (Jones, Gaalaas, Ryu & Frank (1997); Knoll & Brackbill (2002)). Moreover, for these parameters *Tf* = 2 is well within the nonlinear evolution regime for this problem. Snapshots of the *x* and *z* components of the (initially homogeneous) magnetic field at *t* = 2 are shown in Figure 6 for a 2562 mesh simulation computed with a time step of Δ*t* = 0.0025. Throughout this simulation, the number of nonlinear iterations ranged from 1 22 Will-be-set-by-IN-TECH

=5.0, low viscosity, Convergence, Δ t=.1

GEM Bz

Peak K.E. K.E. at t=28.0

quadratic

 −5 4 , 5 4 × −1 2 , 1 2 × −5 4 , 5 4 

Reconnection flux at t=28.0 Reconnection rate at t=28.0

Figure obtained from authors of Reference (Adams et al. (2010)).

<sup>2</sup> tanh(100*y*) + <sup>1</sup>

<sup>256</sup> <sup>512</sup> 1,024 2,048 4,096 8,192 16,384 10−7

N grid (X direction)

between KHI simulations by Chacón and Reynolds et al. because their respective setup and solver tolerances were not necessarily identical. However, these test results are an indication of the type of speed up and code performance one may expect with strongly nonlinear MHD cases with a Newton-Krylov approach. Reynolds et al. (2010) also performed simulation tests

Here we reproduce results from Reynolds et al. (2010) for the 2D KHI test. We set the

the *x*- and *z*-directions, and homogeneous Neumann boundary conditions in the *y*-direction. We initialize the constant fields *ρ* = 1, *b* = (0.1, 0, 10)*T*, *p* = 0.25, and *uy* = *uz* = *By* =

problem employs the resistive MHD equations, with resistivity, viscosity, and heat conduction coefficients set to 10−4, and all runs are taken to a final time of *Tf* = 2. As previous results on this problem suggest that the instability growth rate is independent of the size of the resistivity, such small parameters are natural since the instability is predominantly driven by nonlinear (hyperbolic) effects (Jones, Gaalaas, Ryu & Frank (1997); Knoll & Brackbill (2002)). Moreover, for these parameters *Tf* = 2 is well within the nonlinear evolution regime for this problem. Snapshots of the *x* and *z* components of the (initially homogeneous) magnetic field at *t* = 2 are shown in Figure 6 for a 2562 mesh simulation computed with a time step of Δ*t* = 0.0025. Throughout this simulation, the number of nonlinear iterations ranged from 1

<sup>10</sup> cos(0.8*πx*) + <sup>1</sup>

, with periodic boundary conditions in

<sup>10</sup> cos(0.8*πz*). This

<sup>10</sup> sin(3*πy*) + <sup>1</sup>

Fig. 5. Order of spatial accuracy for simulating magnetic reconnection usinga nonlinear multigrid approach. Error in peak kinetic energy and kinetic energy, reconnection flux rate and reconnection rate, high viscosity cases, *Bz* = 0 (left) and low viscosity *Bz* = 5 (right).

10−6

on the 3D version of the ideal KHI.

computational domain to

0. We then set *ux* = <sup>1</sup>

10−5

10−4

Error

10−3

10−2

Fig. 6. Snapshots of *Bx* (left) and *Bz* (right) in the 2D Kelvin–Helmholtz test at *t* = 2. Figure obtained from authors of Reference (Reynolds et al. (2010)).

Fig. 7. Krylov iterations for the 2D Kelvin–Helmholtz tests. Figure obtained from authors of Reference (Reynolds et al. (2010)).

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 used, the preconditioner results in significantly fewer linear iterations, with the disparity growing as the mesh is refined.

#### **5.4 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

Böhmer, K., Gross, W., Schmitt, B. & Schwarz, R. (1984). Defect corrections and Hartree–Fock

Implicit Numerical Methods for Magnetohydrodynamics 83

Brandt, A. & Dinar, N. (1979). Multigrid solutions to elliptic flow problems, *in* S. Parter

Chacón, L. (2004). A non-staggered, conservative, divB=0, finite volume scheme for

Chacón, L. (2008a). An optimal, parallel, fully implicit Newton-Krylov solver

Chacón, L. (2008b). Scalable parallel implicit solvers for 3D magnetohydrodynamics, *Journal*

Chacón, L. & Knoll, D. A. (2003). A 2D high-*β* Hall MHD implicit nonlinear solver, *J. Comput.*

Dick, E. (1991). Second order formulation of a multigrid method for steady Euler equations

Eisenstat, S. C. & Walker, H. F. (1996). Choosing the forcing terms in an inexact Newton

Harned, D. S. & Kerner, W. (1985). Semi-implicit method for three-dimensional compressible

Harned, D. S. & Schnack, D. D. (1986). Semi-implicit method for long time scale

Harten, A. (1983). High-resolution schemes for hyperbolic conservation-laws, *J. Comput. Phys.*

Hemker, P. W. (1986). Defect correction and higher order schemes for the multigrid solution

Hill, D. J. & Pullin, D. I. (2004). Hybrid tuned center-difference-WENO method for large eddy simulations in the presence of strong shocks, *J. Comput. Phys.* 194: 435–450. Hindmarsh, A. (2000). The PVODE and IDA algorithms, *Technical Report UCRL-ID-141558*,

Hindmarsh, A., Brown, P., Grant, K., Lee, S., Serban, R., Shumaker, D. & Woodward, C. (2005).

Jones, O. S., Shumlak, U. & Eberhardt, D. S. (1997). An implicit scheme for nonideal

Jones, T. W., Gaalaas, J. B., Ryu, D. & Frank, A. (1997). The MHD Kelvin Helmholtz instability II. the roles of weak and oblique fields in planar flows, *Ap.J.* 482: 230–244.

magnetohydrodynamics, *J. Comput. Phys.* 130(2): 231–242.

magnetohydrodynamic computations in three dimensions, *J. Comput. Phys.* 65: 57–70.

of the steady Euler equations, *in* W. Hackbusch & U. Trottenberg (eds), *Multigrid*

SUNDIALS, SUite of Nonlinear and DIfferential/ALgebraic equation Solvers, *ACM*

Greenbaum, A. (1997). *Iterative Methods for Solving Linear Systems*, SIAM, Philadelphia. Greenbaum, A., Pták, V. & Strakous, Z. (1996). Any nonincreasing convergence curve is

through defect correction, *J. Comput. Appl. Math.* 35: 159–168.

possible for gmres, *SIAM J. Matrix Anal. Appl.* 17(3): 465–469.

magnetohydrodynamic simulation, *J. Comput. Phys.* 60: 62–75.

*Applications*, Computing Suppl. 5, Springer–Verlag, Vienna, pp. 193–209. Brandt, A. (1977). Multi-level adaptive solutions to boundary value problems, *Math. Comput.*

31: 333–390.

pp. 53–147.

*Physics Comm.* 163: 143–171.

*of Physics: Conference Series* 125: 012041.

method, *SIAM J. Sci. Comput.* 17(1): 16–32.

*Methods II*, Springer–Verlag, Berlin, pp. 149–165.

*Trans. Math. Software* 31: 363–396.

15: 056103–056103–12.

*Phys.* 188: 573–592.

49(3): 357–393.

LLNL.

method, *in* K. Böhmer & H. J. Stetter (eds), *Defect Correction Methods: Theory and*

(ed.), *Numerical Methods for Partial Differential Equations*, Academic Press, New York,

3D implicit extended magnetohydrodynamics in curvilinear geometries, *Computer*

for three-dimensional visco-resistive magnetohydrodynamics, *Phys. Plasmas*

demonstrates a speedup ranging from 8 <sup>−</sup> 15 for a 1282 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.
