# A comparison of spectral element and finite difference methods using statically refined non

A comparison of spectral element and ?nite di?erence methods using statically re?ned nonconforming grids for the MHD island coalescence instability problem

arXiv:0711.3868v1 [physics.comp-ph] 25 Nov 2007

C. S. Ng1 , D. Rosenberg2 , K. Germaschewski1 , A. Pouquet2 , and A. Bhattacharjee1

ABSTRACT

A recently developed spectral-element adaptive re?nement incompressible magnetohydrodynamic (MHD) code [Rosenberg, Fournier, Fischer, Pouquet, J. Comp. Phys. 215, 59-80 (2006)] is applied to simulate the problem of MHD island coalescence instability (MICI) in two dimensions. MICI is a fundamental MHD process that can produce sharp current layers and subsequent reconnection and heating in a high-Lundquist number plasma such as the solar corona [Ng and Bhattacharjee, Phys. Plasmas, 5, 4028 (1998)]. Due to the formation of thin current layers, it is highly desirable to use adaptively or statically re?ned grids to resolve them, and to maintain accuracy at the same time. The output of the spectral-element static adaptive re?nement simulations are compared with simulations using a ?nite di?erence method on the same re?nement grids, and both methods are compared to pseudo-spectral simulations with uniform grids as baselines. It is shown that with the statically re?ned grids roughly scaling linearly with e?ective resolution, spectral element runs can maintain accuracy signi?cantly higher than that of the ?nite di?erence runs, in some cases achieving close to full spectral accuracy. 1. Introduction spectral methods are often not optimized. Therefore, it is of great interest to develop numerical schemes that can combine high accuracy and high spatial resolution. Adaptive spectral element methods have the potential to do just that, providing spectral-like accuracy that can be applied e?ciently to resolve isolated structures. In this paper, we concentrate on comparing the accuracy of simulation results from a spectral element based AMR code (SE)(Rosenberg et al. 2006; Rosenberg et al. 2006) and a ?nite di?erence based AMR code (FD) (Germaschewski et al. 2006; Bhattacharjee et al. 2005) on an astrophysical problem that requires high spatial resolution as well as high accuracy, the so–called MHD island coalescence instability (MICI) problem (Ng & Bhattacharjee 1998). In order to make meaningful comparisons, we let each code run on essentially the same non-uniform grid (i.e. with the total degrees of freedom in the problem ?xed) that is re?ned a priori in regions of the grid where the current sheets will form in the MICI. A separate pseudo-spectral code (PS), running on uni-

In many hydrodynamic or magnetohydrodynamic (MHD) applications in astrophysics or space physics, it is essential that a numerical simulation resolve the development of sharp spatial structures accurately. While pseudo-spectral methods generally can maintain high accuracy, they are mainly applied on more regular geometry and require more uniform grids, which can make it di?cult to reach high resolution in order to resolve sharp isolated structures especially in ?ows dominated by such structures. Static or adaptive mesh re?nement (AMR) methods can put more grid points in and around isolated structures in order to resolve them, but may not achieve similar high-accuracy if low order spatial discretizations are used (Rosenberg et al. 2006). They are particularly useful in bounded ?ows, where pseudo–

1 Space Science Center, University of New Hampshire, 39 College Road, Durham, NH 03824, USA 2 TNT/IMAGe, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307-3000, USA

1

form grids, is also used to provide a baseline for the comparisons. There is little in the literature about attempts to compare FD and high order methods. Often, it is simply accepted that the ?xed lower order FD methods will be less accurate, but that, due to their simplicity and e?ciency, h–type grid re?nement can always be carried out to improve the accuracy for any real problem. We do not attempt in this work to compare performance metrics of the methods, preferring instead to focus on their ability to produce accurate solutions. To an extent, then, it is clear that h-re?nement will improve solution accuracy for many problems; however, previous work (Rosenberg et al. 2006) suggests that local high order may be required in certain instances. We consider this issue of h-re?nement in FD brie?y in this work. Previous comparative work in one dimension (Basdevant et al. 1986) demonstrates that spectral methods, including SE, are more accurate than FD methods due primarily to dispersion problems in the latter. Nevertheless, this work concludes that the SE method is not well suited to the calculation of thin internal layers (structures), particularly when their location is unknown. This assertion is made because, while accuracy is found to be good, Basdevant et al. 1986 ?nd that polynomial orders are required to be inordinately high even in one dimension. This conclusion is refuted in later work (Mavriplis 1994), which found that the SE method is indeed well suited to this type of problem provided adaptivity is used (see also (Rosenberg et al. 2006)). More recent work (St-Cyr et al. 2007) also demonstrates that the SE method exhibits smaller errors than a loworder ?nite volume scheme at comparable resolutions for nearly all of a series of tests of a shallow water model on a cubed–sphere grid that is either adaptively or statically re?ned. However, to the best of our knowledge, ours is the ?rst work that performs a quantitative comparison of FD with SE methods in divergence–free nonlinear ?ows in two space dimensions with nonconforming (statically) adaptive grids. We continue our introductory remarks on the MICI in section 2, where we present details of the problem, and provide a motivation for our work. In the process, we discuss important properties of the MICI as well as some aspects of MHD ?ows

that will serve our later discussion. We provide details about the simulation set up as well as some diagnostic measures used in the comparisons in section 3. In section 4 we present the numerical methods used in su?cient detail to elucidate the results. Our results are presented in section 5, where we o?er several types of accuracy comparisons relevant to MHD ?ows, and the MICI, in particular. We conclude in section 6 with a summary of our ?ndings, and some discussion about the relative advantages of the methods. 2. MHD Island Coalescence Instability

In this section, we provide a brief background for the MICI problem that motivates us to perform our comparative study in this particular simulation con?guration. It is well known that a substantial part of our universe is composed of systems of plasmas, ionized gases, and conducting ?uids. Magnetic ?elds, both ?uctuating and large scale, play an important role in the physics of these systems. In fact, large-scale magnetic ?elds have long been observed to exist in the solar corona (stellar coronae, see e.g. (Parker 1979) and references therein), interstellar space (within a galaxy, see e.g., (Forman et al. 1985; Grimes et al. 2005)), galaxy clusters (see e.g., (Kellogg et al. 1973; Sarazin 1986)). We consider the representation of magnetohydrodynamics (MHD), as a starting point of discussion. For simplicity, we consider the incompressible MHD equations with constant mass density ρ0 , ?t u + u · u=? p+j×b+ν × (u × b) + η

2 2

u,

(1) (2) (3)

?t b =

b

· u = 0,

·b=0

where u and b are the velocity and magnetic ?eld √ (in Alfv?n velocity units, b = B/ ?0 ρ0 with B e the induction and ?0 the permeability); j = × b is the current density; p is the pressure divided by the mass density; and now the normalized viscosity ν is basically the inverse of Rv , and resistivity η is basically the inverse of S. In a dimensionless form, in which all physical quantities are measured by their typical values, the MHD equations generally have dissipation terms involving higher spatial derivatives of the 2

?eld quantities with strength characterized by the inverse of dimensionless parameters: Rv = V L/ν is the Reynolds number (with V a typical ?ow speed, L a typical length scale), S = VA L/η is the Lundquist number (with VA a typical Alfv?n e speed). We can immediately see that for most astronomical length scales L, both Rv and S are very large numbers such that the dissipation terms in the MHD equations can be thought usually to be ignored, (i. e., the ideal MHD equations), except possibly in regions where there exist steep spatial gradients. In ideal MHD, in which ν = η = 0 in Eqs. (1)(3), several quadratic quantities, namely energy, magnetic helicity and cross helicity are conserved exactly in three dimensions. Moreover, a magnetic ?eld line is carried by the ?ow velocity so that the topology of the magnetic ?eld con?guration cannot be changed in ideal evolution. Thus, many important physical processes, such as magnetic ?eld line reconnection, local ?uid heating, and particle acceleration due to parallel electric ?elds, are disallowed. In nonideal MHD, such processes can occur within boundary layers, which are regions of high spatial gradients in the current density and vorticity, described by singular perturbation theory. The tendency for formation of current and vortex singularities in the ideal equations, if and when it occurs, is a phenomenon of great interest because the sites of singularity formation are precisely the sites where astrophysically signi?cant physical processes such as heating and particle acceleration can occur. E. N. Parker has argued for over three decades that current sheets, or tangential discontinuities of the magnetic ?eld, do generally exist in a magnetic equilibrium (Parker 1972; Parker 1979; Parker 1994). If Parker is correct, then the formation of current sheets, and the subsequent ?uid heating and particle acceleration from the release of magnetic energy due to the rapid dissipation of these sheets can have very signi?cant astronomical consequences. One observable consequence is the production of X-rays. The solar corona was the earliest astronomical object that was observed to be emitting X-rays. Based on these observations, the temperature of the solar corona is estimated to be of the order of a million degrees with peak radiation at wavelengths about 30 A (i. e., in the soft X-ray range). While there may be many

di?erent heating mechanisms involved in heating the corona, it is argued that magnetic ?elds must play an important role among these processes. See (Ng & Bhattacharjee 2008) for our recent discussion on solar coronal heating theory and more references, e.g., (Klimchuk 2006). In Parker’s model, a solar coronal loop is treated as a straight ideal plasma column, bounded by two perfectly conducting end-plates representing the photosphere. The footpoints of the magnetic ?eld in the photosphere are frozen (linetied). Initially, there is a uniform magnetic ?eld along the z direction. To simplify the consideration, we may keep the footpoints of the magnetic ?eld on one of the plates (z = 0 ) ?xed, while the footpoints on the other plate (z = L) are subjected to slow, random motions that deform the initially uniform magnetic ?eld. The footpoint motions are assumed to take place on a time scale much longer than the characteristic time for Alfv?n wave propagation between z = 0 e and z = L, so that the plasma can be assumed to be in static equilibrium nearly everywhere, if such equilibrium exists, during such random evolution. For a given equilibrium, a footpoint mapping can be de?ned by following ?eld lines from one plate to the other. Since the plasma is assumed to obey the ideal MHD equations, the magnetic ?eld lines are frozen in the plasma and cannot be broken during the twisting process. Therefore, the footpoint mapping must be continuous for smooth footpoint motion. (Parker 1972) claimed that if a sequence of random footpoint motions renders the mapping su?ciently complicated, there will be no smooth equilibrium for the plasma to relax to, and tangential discontinuities of the magnetic ?eld (or current sheets) must develop. Parker’s claim has stimulated considerable debate (van Ballegooijen 1985; Antiochos 1987; Zweibel & Li 1987; Longcope & Strauss 1994; Cowley et al. 1997) that continues to this day. For review and extensive references, see (Low 1990; Browning 1991; Parker 1994). The ?rst signi?cant objections to Parker’s claim of non-equilibrium was raised by van Ballegooijen (1985), who argued that smooth equilibria must always exist as long as the footpoint motion is smooth (or continuous). van Ballegooijen (1985) developed his argument on a reduced form of the MHD equations (referred to hereafter as the RMHD equations) originally 3

derived by Strauss (1976) for a low β plasma. These equations, which also provide the basis for many later developments in this work, are: ?t ? + [φ, ?] = ?z J + [A, J] + ν ?t A + [φ, A] = ?z φ + η

2 ⊥ A, 2 ⊥ ?,

(4) (5)

where A is the ?ux function so that the magnetic ?eld is expressed as b = ? + ⊥ A × ?, φ z z is the stream function so that the velocity ?eld is expressed as u = ⊥ φ × ?, ? = ? 2 φ is z ⊥ the z-component of the vorticity, J = ? 2 A ⊥ is the z-component of the current density, and [φ, A] = ?y φ?x A ? ?y A?x φ. An ideal magnetostatic equilibrium solution of Eqs. (4) and (5) is obtained by setting all explicitly time-dependent terms, as well as φ and η to zero. We then obtain ?z J + [A, J] = 0 , (6)

which can also be written as b · J = 0. This implies that the current density J must be constant along a given magnetic ?eld-line in an ideal static equilibrium. Based on this set of equations, Longcope and Strauss (1994) argued that even when the magnetic equilibrium is unstable, e.g. due to MICI, it will only relax to a second equilibrium (assuming that it exists) with very thin current layers with thickness less than about 10?7 of the large scale. However, another point of view was raised by Ng and Bhattacharjee (1998), who argued based on a mathematical theorem on the RMHD system that for a given ?xed footpoint mapping between z = 0 and z = L, there exists only one smooth equilibrium. This means that an unstable equilibrium will relax ideally to a ?nal state with current sheets (tangential discontinuities). This scenario has very di?erent implications than those predicted by the Longcope and Strauss 1994, since energy dissipation and other energetic e?ects can be much stronger for the case with current sheets, than the case with smooth but thin current layers. Therefore, it is very important to determine which of these two scenarios should actually occur. However, due to the fact that the current layers predicted by (Longcope & Strauss 1994) are very thin, it is beyond our computational ability if the full simulation volume is to resolve to the same small scale. To provide a resolution to this problem with current computer architectures, one may 4

need to apply AMR techniques that put more grid points in the regions where distinct structures appear. At the same time, to ensure that any such numerical study is actually representative of the true dynamic solution, one needs to make sure that the numerical scheme used can maintain a reasonably high accuracy. Hence, accuracy becomes a particularly important factor in the choice of the numerical method for the Parker problem of MICI. Because ?nite di?erence–based schemes are usually of low order truncation, accuracy typically decreases as higher order derivatives are taken. Also because of low order, these methods can be di?usive (and dispersive). For the RMHD equations (4) and (5), one needs to use the spatial derivatives of J, or third order spatial derivatives of A. A couple of questions arise in simulations based on such schemes: Will the lower accuracy in calculating these higher spatial derivatives change drastically the dynamical properties of the problem? Will the numerical di?usion preclude the formation of a current sheet, and lead instead to a smooth residual current layer? For example, in the Parker problem, one will need a simulation that is accurate enough so that we can show con?dently whether there is a formation of a true current sheet as predicted by Ng & Bhattacharjee (1998), or if the current layers actually tend to ?xed (but small) thickness (as small as 10?7 of the large scale) as predicted by Longcope & Strauss (1994). In order to increase the reliability of the simulations, it is thus of great interest to look for schemes that can maintain higher accuracy, and can make use of AMR techniques to resolve sharp features at the same time. Spectral element methods have the potential to ful?ll such requirements. This is what motivates us to perform the present comparative study. 3. Simulation set up and diagnostics

While the main problem of interest is the three dimensional (3D) MICI problem represented by Eqs. (4 )-(6), for simplicity we will instead simulate the two dimensional (2D) version of the problem for the purpose of this comparative study. Examining this problem in 2D should not a?ect greatly di?erences in accuracy among di?erent codes. The 2D MICI problem can be viewed as

the limit of the 3D problem with L → ∞. In this limit, it has been shown that current singularities must form for the ideal equations (η = 0) (Longcope & Strauss 1993). This is bene?cial for our present study since we know that there must be speci?c structures, and we know where they will appear. We can then compare how well these sharp structures are resolved in di?erent schemes. Therefore, we will simulate the set of equations without the z-dependence, or ?t ? + [φ, ?] = [A, J] + ν ?t A + [φ, A] = η

2 ⊥ A, 2 ⊥ ?,

(7) (8)

where the de?nitions for the variables are the same as in Eqs. (4 )-(6). Periodic boundary conditions are used in both the x? and y? directions. To be speci?c, the simulation domain is set to be 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1. The initial equilibrium is chosen to be A(x, y, t = 0) = A0 sin(2πx) sin(2πy) , (9)

The grids are generated separately for each Rv we consider, and shown in the following sections. The equivalent resolution Req for the FD and SE simulations is de?ned as that which would be achieved if the most ?nely resolved subdomain covered the entire domain. If E0 is the linear number of elements then for SE, equivalent linear resolution, Req , is computed from the initial E0 × E0 grid, the number of re?nement levels, , and the expansion polynomial degree, N , such that Req = [2 N E0 ]. For all SE runs, E0 = 8, N = 8, and varies with Rv . For the FD method, in order to make more direct comparisons, we use the same grids as used in SE, with a 8 × 8 uniform grid within each element when 8th order polynomials are used in SE. For the PS method, a uniform grid of Req × Req collocation points is used. Table 1 contains a list of the viscosity (resistivity) used and the corresponding Req . In the absence of external forcing, viscosity and magnetic resistivity, the 2D MHD equations (7) and (8) have three ideal invariants: the energy 1 2 u + b2 = EK + EM , (11) 2 composed of the kinetic and magnetic energy, the L2 norm of the magnetic potential E= M = A2 , and the cross helicity H = u·b . (13) (12)

with a small initial ?ow, which in terms of the stream function is φ(x, y, t = 0) = φ0 [cos(2πx) ? cos(2πy)] , (10)

where A0 = 0.4 and φ0 = 0.002. φ0 is chosen small enough so that there is a clear linear phase in the growth of the MICI. Because SE evolves a di?erent form for the equations (section 4.1), these initial conditions are converted into conditions on b by using b = × A ?, and on u by using the relation z u = × φ ?. z Fig. 1 shows the contour plots of A and φ at t = 0 as given by Eqs. (9) and (10). The velocity deduced from this stream function, has the initial tendency to push the two islands at the lower right and upper left toward each other and towards the center. Eventually, these two islands will merge with each other; hence, the nomenclature island coalescence. The forms in (9) and (10) should preserve additional symmetries with respect to the lines x = y and x = 1 ? y: A(x, y) = A(?x, ?y) = A(y, x), φ(x, y) = φ(?x, ?y) = ?φ(y, x) . Speci?cally, φ = 0 on the lines x = y and x = 1 ? y. These symmetries are not incorporated into the simulation schemes; however, they can provide information about how well a numerical scheme preserves them. 5

For all functions φ in these de?nitions, we have φ ≡ φ dx2 . In 2D assuming bi-periodicity, it is easy to show using Eqs. (7) – (8) that dE = ?ν ?2 ? η J 2 , dt (14)

Table 1: Parameters used in the simulations described in the following sections; Req is the linear grid resolution, and ν (η) is the viscosity (resistivity). Req 128 256 512 ν=η 2 × 10?3 1 × 10?3 3 × 10?4

dM = ?η b2 /2 , dt and dH = ?(η + ν) J? . dt

(15)

(16)

into the d-dimensional variational form of equations (17)-(18) on a domain D, and using appropriate quadrature rules, we arrive at a set of semi– discrete equations written in terms of spectral element operators: M dZ± j dt = ?MC Z± + DT p± j j ?ν± LZ± ? ν LZj j Dj Z± j = 0, (19) (20)

Equations (14)-(16) must hold for any 2D periodic MHD ?ow, and they serve as critical diagnostics for any numerical treatment of the incompressible MHD equations. In addition to the conservation laws, there are several other quantities of interest in diagnosing numerical solutions to MICI as de?ned later when we present our numerical results. 4. Numerical methods

We provide in this section the numerical methods we use for carrying out the simulations on the MICI problem. These methods have been described in detail in other sources, but we provide enough explanation so that the simulation setup, results, and discussion will be comprehensible, and the paper reasonably self–contained. 4.1. Spectral element method

This method evolves Eqs. (1)-(3) in time, as part of the Geophysics/Astrophysics Spectral Element Adaptive Re?nement (GASpAR) code, and has been described in detail in (Rosenberg et al. 2006). Here we present aspects of that description that relate to the results discussed in this work. Equations (1)-(3) are solved in Els¨sser form a (Els¨sser 1950): a ?t Z± + Z · Z± + p ? ν ±

2

Z± ? ν

2

· Z± = 0 ,

Z =0 (17) (18)

where Z± = u ± b and ν ± = 1 (ν ± η). Thus, 2 we solve for u and b in terms of Z± . All components of Z± and those of the velocity and magnetic ?eld reside on Gauss–Lobatto–Legendre (GL) nodes, while pressure resides on Gauss– Legendre (G) nodes. These choices follow from the ?nite dimensional subspaces to which these quantities are restricted: Z± (and u and b) are expanded in terms of N th order GL polynomials, while p is expanded in terms of N ? 2th order G polynomials. Substituting these expansions

for the j th component. The variables Z± represent values of the Z± collocated at the GL node points, and p± are values of the pressures at the G node points. Note that, even though the continuous equations contain only a single pressure, Eq. (19) contains one for each Els¨sser veca tor because the constraints (20) are enforced separately on them. The operators M, L, and C, are the well–known mass matrix, weak Laplacian and advection operators, respectively (e.g. , (Rosenberg et al. 2006), and references therein), and Dj represent the Stokes operators, in which the GL basis function and its derivative operator are interpolated to the G node points, and multiplied by the G quadrature weights. All d–dimensional operators are computed as tensor products of their component 1D operators. Note that because di?erent expansion bases are used for the vector quantities as are used for pressure, a staggered grid is implied. Hence, this method is referred to as the PN ? PN ?2 formulation. It was chosen to prevent spurious pressure modes (Maday et al. 1992; Fischer 1997). Note also the e?ect of the Dj , which act on so–called v-grid (vector) quantities: they take a derivative that itself resides on the G nodes; hence, the discrete divergence, Eq. (20) resides on the p-grid. The e?ect of the transposed Stokes operator, DT , on the other j hand, is to compute a derivative–residing on the v–grid–of a p-grid quantity. The code has an adaptive mesh capability that we utilize minimally in this work. The connectivity and adaptivity algorithms were described in (Rosenberg et al. 2006). No dynamic adaptivity is used here (see section 3). Instead, the nonconforming grid is constructed initially by turning o? the re?nement criteria, and selecting the elements we want to re?ne explicitly. The nonconforming grid is then used in a static con?guration through-

6

out the simulation. Nonconforming in this context means that there are at most two child elements adjacent to a coarser neighbor; an element is re?ned by dividing it isotropically into 2 × 2 child elements, each of which contains the same polynomial order as the parent. 4.1.1. Time stepping with a constraint

In our simulations we must resolve all temporal (and spatial) scales, so we use an explicit Runge–Kuttta (RK) method for integrating Eqs. (19)-(20) in time. The speci?c RK algorithm is that presented in (Canuto et al. 1988, p. 109), known to be valid to second order in ?t (Brachet et al. 2007), which we use for all computations. At each stage, we can write (recall eq. (19)): Z± j = Z±,n ? j 1 ?t M?1 (MC Z± ? DT p± j j k ± +ν± LZj + ν LZj ). (21)

are used. In this case, we must ensure that all quantities in the subspace represented by the GL grid remain continuous across element interfaces. The speci?c way in which this is carried out in the code is described in (Rosenberg et al. 2006), and entails exchanging interface data. Using the Boolean scatter matrix, Ac , the interpolation matrix from global to local degrees of freedom, Φ, and the masking matrix (that enforces homogeneous boundary conditions), Π, that were presented there (see also Fischer et al. 2002 ), it is found that the local Stokes operators in equations (19)-(20) must be replaced with Dj → DL,j ΦAc Π, where DL,j = diagk (Dk ), and the Dk are the maj j trices from above, and k indexes the elements in the domain. It is su?cient for our purposes to refer to the local form of the Stokes operators, and to simply observe that communication occurs when they are applied. 4.2. Finite di?erence method

We require that each RK stage obey eq. (20) in its discrete form, so multiplying eq. (21) by Dj , summing, and setting the term Dj Z± = 0 leads j to the following pseudo-Poisson equation for the pressures, p± :

± Dj M?1 DT p± = Dj gj , j

(22)

where the quantity

± gj =

1 ?t M?1 (MC Z± +ν± LZ± +ν LZj )?Z±,n j j j k

is the remaining inhomogeneous contribution. Rosenberg et al. (2007) describe how M?1 is computed. Equation (22) is solved using a preconditioned conjugate gradient method (PCG); for all computations reported here, we use a block Jacobi preconditioner computed using a fast diagonalization method. Thus, at each time-step, two RK stages are computed, and each stage solves eq. (22) twice, once for Z+ , and once for Z? leading to 4 pseudoPoisson solves at each time step. 4.1.2. Communication

While equations (19)-(20) are correct for a single element, strictly speaking, they are not complete when multiple subdomains (elements) 7

The Magnetic Reconnection Code (MRC) is a suite of codes (Germaschewski et al. 2006; Bhattacharjee et al. 2005) which integrate various reduced and extended ?uid models of plasma ?ows. In this paper, the 2D AMR version of the code integrating the equations of RMHD has been used. The MRC employs a hierarchical quadtree based approach in block-structured adaptive mesh re?nement. At each re?nement time, every (square or rectangular) box is checked as to whether the data in that box requires additional resolution. The re?nement criterion needs to be selected appropriately for the given problem, from simple evaluation of maxima or gradients to a Richardson extrapolation based approach. If the local resolution in a box is considered insu?cient, it is subdivided into 2 × 2 smaller boxes, which have physically twice the resolution but the number of grid points in each box remains the same it was on the coarse parent box. As opposed to the alternate approach of allowing for arbitrary rectangular patches of re?nement, this method results in simpler data management, easier optimization (each box is always the same size of n × m grid points, important for cache con-

siderations, etc) and more e?cient load balancing, even though it introduces some ine?ciency in that some additional areas are re?ned where the higher resolution is not required. For load-balancing, we use a space-?lling Hilbert Peano curve which connects all the boxes at various levels of re?nement. The boxes along this one-dimensional curve are then evenly distributed to the available processors. Since each patch has the same number of grid points, the computational load is evenly spread, and as the Hilbert-Peano curve has the property that in a certain averaged sense, patches which are close in the two-dimensional domain are also close on the one-dimensional curve, spatially close regions are clustered onto the same processor; that is, it maintains data locality. Since communication is only necessary between spatial neighbors to exchange boundary data, most of the needed data is available on the local processor, expensive MPI communication is only required for boundary data transfer across the boundaries between clusters on di?erent processors, which is e?ectively minimized. An additional di?culty in using AMR to integrate the reduced models – not encountered in the purely hyperbolic systems – is the need to solve elliptic equations, e.g. solving for the stream function φ from the vorticity ?. To discretize this sub-problem, we rewrite the Laplacian 2 as the divergence of a gradient · and apply a conventional ?nite volume method for a conservation law. Combined with correcting the ?uxes at ?ne-coarse boundaries, this provides an exactly numerically conservative discretized expression. To e?ciently solve this elliptic problem we use a variation of the fast adaptive composite (FAC) method (Germaschewski et al. 2004; McCormick 1989), which, due to the multilevel character, provides an iterative solver with very fast convergence. Alternatively, the discretized Laplacian on the AMR hierarchy can be calculated explicitly as a sparse matrix, and PETSc’s (Balay et al. 2004) rich supply of solvers are available to solve Poisson’s equation. In particular, SuperLU (Li & Demmel 2003) in many cases proves to be very fast, since the expensive LU decomposition only needs to be done once and then can be applied for many time-steps until the AMR hierarchy 8

of boxes changes. This is the method used in the simulations for this study. Our AMR code can also be run in fully implicit mode using Crank-Nicholson time-stepping, however no preconditioner has been developed yet, so we are using a simple unpreconditioned NewtonKrylov solver which does not achieve optimal performance, but is still faster than a Runge-Kutta explicit scheme. Most of the runs presented here are using the implicit time integrator. Some runs have also been cross-checked by using an adaptive time-step Runge-Kutta explicit method, which is included with PETSc’s. Spatial discretization is using a regular second order central di?erences, which is equivalent to a conservative ?nite volume method, with no upwinding employed. For the purpose of the present study, the dynamic re?nement capability of the AMR code is turned o? so that we can use the same grids that were also used in SE. This was done because it is di?cult to set re?nement methods and criteria to be the same in both the FD and SE codes. Using static re?nement provides us with nonconforming grids over which we can exert complete control of the number of degrees of freedom. Comparison of re?nement criteria in the SE and FD codes and the e?ect on the resulting dynamics is left for future work. 4.3. Pseudo-spectral method

The PS code is based on fast Fourier transform (FFT) on a 2D bi-periodic domain. It is de-aliased by the standard 2/3 rule. The nonlinear term is calculated in the physical space on a uniform grid of collocation points. A second order predictor-corrector method is used for time integration. The code is parallelized using a parallel version of the FFT. For this study, results from SE and FD methods for the Req case are compared with those from the PS code with Req × Req collocation points. The results from the PS code are themselves checked with PS runs with higher resolutions, up to 2048 × 2048, which con?rm that the results with the original resolution are already well resolved.

5.

Numerical results

Fig. 1.— Contour plots of (a) the initial ?ux function A, and (b) the initial stream function φ. Positive (including zero) contours are solid and negative contours are broken. Contours levels are from ?0.4 to 0.4 with an increment of 0.025 in (a), and ?0.004 to 0.004 with an increment of 0.00025 in (b).

Fig. 2.— Contour plots of (a) A, and (b) φ at t = 1 for the Req = 256 case. Positive (including zero) contours are solid and negative contours are broken. Contours levels are from -0.4 to 0.4 with an increment of 0.025 in (a), and -0.1 to 0.1 with an increment of 0.00625 in (b). Note the sheet formation (a) and the ?uid acceleration (b).

In this section, we compare simulation solutions by SE and FD methods, using solutions by PS as references. Note that each code (SE and FD) has been veri?ed separately on a suite or problems including analytical solutions, but the purpose of the present work is to explore the fully developed nonlinear regime of the MICI problem, when sharp current and vorticity layers appear. The main dynamics of MICI is due to the attractive force between two islands with the same sign of current. The initial small perturbation in the ?ow velocity breaks the unstable equilibrium. After a short transient, a linear phase appears when kinetic energy increases exponentially so that the ?ow velocity pushes the two islands together further. Eventually a sharp current layer appears between the two islands when the dynamics enters the nonlinear phase at time t ? 1, in units of the Alfv?n time. Magnetic reconnection e then proceeds faster, producing a strong out?ow velocity, as well as sharp vorticity layers. Fig. 2 shows contour plots of A and φ at t = 1 for the Req = 256 case. We see that the two islands with positive A are pushed towards each other with some A contours (magnetic ?eld lines) already reconnected, as compared to Fig. 1. The stream function φ shows strong out?ows as indicated by the concentration of stream lines. The grids used in our simulations for the three cases are shown in Fig. 3 to 5. The re?nements with each grid are imposed so as to resolve structures produced by the above dynamics, for the di?erent dissipation levels indicated in Table 1. In Fig. 3, we also plot the color contours of the current density J at t = 1.3 produced by the SE run for the Req = 128 case, with red representing the positive end and blue representing the negative end of J values (as in the following). Only one level of re?nement is used that covers the region containing stronger J. Note that within each square, polynomials of order N = 8 are used for SE, while an 8 × 8 uniform grid is used for FD (the same as in the other two cases). In Fig. 4, the color contours of the vorticity ? at t = 1.3 produced by the SE run are plotted along with the grid for the Req = 256 case. Note that there are now two levels of re?nement. In Fig. 5, the color contours of J at t = 0.93 produced by the SE run

9

Fig. 4.— The grid used in the Req = 256 case. The background is the color contour plot of the vorticity ? at t = 1.3 produced by the SE run, with red representing the positive end and blue representing the negative end of ? values. Two levels of re?nement are used here.

Fig. 3.— The statically re?ned grid used in the Req = 128 case. The background is the color contour plot of the current density, J, at t = 1.3 produced by the SE run, with red representing the positive end and blue representing the negative end of J values. Within each square, polynomials of the order of N = 8 are used for SE, while an 8×8 uniform grid is used for FD. Only one level of re?nement is needed at this (low) Reynolds number.

Fig. 5.— The grid used in the Req = 512 case. The background is the color contour plot of the current density J at t = 0.93 produced by the SE run, with red representing the positive end and blue representing the negative end of J values. Now, three levels of re?nement are used.

10

are plotted along with the grid for the Req = 512 case, with one further level of re?nement. At the same time, we found that the largest squares have to be re?ned for this case due to a much smaller dissipation in this case. Before presenting our numerical results, let us look at the degrees of freedom (DOF) of these grids at the di?erent Req levels, as shown in Fig. 6 in red asterisks; they follow roughly a linear scaling proportional to Req . Also plotted are DOF of the PS runs using uniform grids, in blue diamonds, which follow a scaling of R2 as they should. Since eq as Req increases, the di?erence between the two scalings can be very large, using adaptive grid re?nement has the potential to provide considerable savings in memory and/or CPU usage, if the linear scaling of the DOF in these adaptive grids continues to hold for even larger Req , and thus for large Rv and S. 5.1. Accuracy of the conservation laws Fig. 6.— DOF of the grids shown in Fig. 3 to 5 at the di?erent Req levels, in red asterisks, and those of the PS runs using uniform grids, in blue diamonds. A dashed line showing a R2 dependence eq in the PS runs, and a dotted line showing a linear dependence of Req of the statically–re?ned runs are also plotted.

We start our comparison by looking at how each code preserves conservation laws as shown in Eqs. (14) and (15). We do not include the conservation law for cross helicity, Eq. (16), in our comparison, since H = 0 exactly based on our initial conditions in Eqs. (9) and (10), and thus dH/dt = 0; all three codes preserve H = 0 and dH/dt = 0 well during the duration of the runs for all three cases. This is more due to how well the structure of each code preserves symmetries, rather than due to numerical accuracy. Also, since both the left hand side (LHS) and right hand side (RHS) of Eq. (16) are close to zero, taking the difference between the two to see fractional changes will not yield meaningful results. In Figs. 7 to 9, the fractional di?erence between the LHS and RHS of (a) Eq. (14), and (b) Eq. (15) are plotted as functions of time, for the three values of Req . The fractional di?erence (or error) ? is de?ned as | LHS ? RHS | / | RHS |, with the time derivative in the LHS calculated by taking ?nite di?erence of the output in time, thus providing an overestimate of the error. In these ?gures, black curves are for PS runs, red curves are for SE runs and blue curves are for FD runs. For the Req = 512 case, we present results up to t ? 1, since both SE and FD experience larger error after this time, probably due to the fact

Fig. 7.— Fractional error in energy conservation ˙ law in (a) ?E, and for the L2 norm of the mag˙ netic potential, ?M in (b) for the Req = 128 case as functions of time. Black curves are for the PS run, red curves are for the SE run, and blue curves are for the FD run (same in the other ?gures below) .

11

Fig. 8.— Fractional error in conservation laws in ˙ ˙ (a) ?E, and (b) ?M for the Req = 256 case as functions of time. See Fig. 7 for de?nitions.

Fig. 9.— Fractional error in conservation laws in ˙ ˙ (a) ?E, and (b) ?M for the Req = 512 case as functions of time. See Fig. 7 for de?nitions.

that we are using a ?xed adaptive grid and could not follow the development of small scales closely enough. To achieve more accurate results, both codes would have to employ dynamic adaptive re?nement. However, it is di?cult to make certain the two codes re?ne and coarsen in the same way in order to make meaningful comparisons; hence, this investigation is left for future work. For all three levels of Req , we see that PS results (black) preserve conservation laws the best, as expected. This is why we may use them as baselines ˙ for comparisons. The error level in ?M as shown ˙ ≡ d A2 /dt) is around or in (b) panels (with M ˙ slightly over 10?6 , while the error ?E as shown in (a) panels is around or slightly below 10?5 , since energy involves one more spatial derivative than A and thus is less accurate in this computation. For SE runs (red), the accuracy turns out to be quite good, given that they are running on stati˙ cally re?ned grids. The error ?E is more or less one order of magnitude above that of PS runs for all three cases, but still is in a low level of about ˙ 10?4 . For the error ?M , the di?erence between SE and PS becomes larger (about two orders of magnitude), at a level of about 10?3 . The reason ˙ ˙ ?M is greater than ?E for SE runs is due to the fact that b and u are the primary variables in the computations, and A is a ?eld quantity derived from b, by solving the equation 2 A = ?J. This process introduces error in A and thus the error ˙ ˙ ?M is greater than ?E. For FD runs (blue), the code does capture quantitatively the main dynamical evolution of the MICI problem, although with a lower accuracy. ˙ For ?M , the error level is up to above 10?2 , and ˙ the error ?E can go above that, even reaching as ˙ high as 10?1 or more. The error in ?M is at a ˙ level slightly lower than that of ?E in this version of the FD code, similar to the trend of PS, since it uses A as primary variable. From these comparisons of the conservation laws, we see that the accuracy of SE is higher as expected for spectral–based methods; it sometimes approaches that of PS runs using uniform grids. This con?rms that SE can deliver near spectral accuracy, the main advantage of employing such a scheme. Using the same adaptive grid, FD runs show lower accuracy. It is conceivable that the accuracy of the FD scheme can be improved somewhat by algorithm modi?cations. However, it 12

is unlikely that it can be improved to the level of the SE runs, when constrained to the same grids. 5.2. Accuracy with respect to higher order derivatives

In the above comparison of conservation laws, we see that the accuracy level can change with respect to ?eld quantities involving a di?erent order of spatial derivatives of A. We now look into this issue further by comparing integrated ?eld quantities. Instead of presenting comparisons for all three levels of Req , we will only present ?gures based on the Req = 128 case, and just mention that similar conclusions can be drawn for the other two cases. In Figs. 10 to 13, we show in the (a) panels the time series of the quantities A2 , E = 1 b2 + u2 , 2 J 2 , and ( J)2 , for the PS run (black), SE run (red), and FD run (blue), plotted in this order for the Req = 128 case. These four plots involve ?eld quantities of increasing order of spatial derivatives of the magnetic potential A. When plotting this way, it is not easy to see the di?erence between the curves when they are close to each other. In fact, red curves almost totally cover black curves for all cases, and blue curves almost totally cover the other two in A2 and E. So in the (b) panels, we plot the fractional di?erence ? between these values and those from a “converged” PS run using a much higher resolution with a uniform grid of 2048 × 2048. Again, black is for the PS run (fractional di?erence between the Req = 128 PS run and the 2048 PS run), red is for the SE run, and blue is for the FD run. The fractional di?erence (or error) ? is de?ned as, e.g., ? A2 ≡| A2 ? A2 2048 | / | A2 2048 |. We note from the black curves that the PS run with a 128 × 128 grid is indeed a highly converged run, in the sense that the fractional error with respect to the 2048 × 2048 is very small. The error ? A2 is at a level of about or slightly over 10?8 . This increases to around 5 × 10?7 for ?E, around 10?6 to 10?5 for ? J 2 , and with a maximum value slightly below 10?3 for ? ( J)2 . The trend of obtaining less accurate results for quantities involving higher derivatives is expected. Still, for all quantities that are important to the numerical integration, i.e. up to J, the accuracy of the 128 × 128 run is high enough for the whole duration of the simulation. Again, this shows that the 13

Fig. 10.— (a) Time series of A2 for the PS (black), SE (red), and FD (blue) runs plotted in this order for the Req = 128 case. All three runs overlap. (b) Fractional error in ? A2 , as compared with a PS run with 2048 × 2048 resolution.

Fig. 11.— (a) Time series of total energy E for the PS (black), SE (red), and FD (blue) runs plotted in this order for the Req = 128 case. Again, all three runs are nearly coincident. (b) Fractional error in ?E, as compared with a PS run with 2048 × 2048 resolution.

Fig. 12.— (a) Time series of J 2 for the PS (black), SE (red), and FD (blue) runs plotted in this order for the Req = 128 case. Di?erences in the dissipation η J 2 are observed. (b) Fractional error in ? J 2 , as compared with a PS run with 2048 × 2048 resolution.

PS run can be used as a baseline for comparison. For the SE run, the error level is higher than that of the PS run. This is qualitatively similar to the comparison of conservation laws, except now it is higher by a somewhat larger amount, sometimes about two orders of magnitude. Speci?cally, ? A2 is at a level of about or a little over 10?5 , ?E is up to about 5 × 10?4 , ? J 2 is about or a little over 10?3 , and ? ( J)2 is about or a little over 10?2 . This level is still acceptably low, but is somewhat higher than those in the error of conservation laws as shown above, especially for ? J 2 and ? ( J)2 . However, we need to keep in mind that the errors in the conservation laws are indicating how well a code simulates the equation self-consistently at each moment, but the errors in comparing with results from a converged solution are accumulated over time, and can thus be larger than the former. For the FD run, the error level is again higher than that of the SE run, sometimes by an order of magnitude or more. The error level of ? A2 is about 10?3 or slightly above, while ?E is at a level or slightly over 10?2 . These two are still reasonably small and so we do not observe appreciable di?erences in their respective (a)-panel plots. However, ? J 2 and ? ( J)2 become large enough to be observable in the time series plots themselves. Quantitatively, ? J 2 is at a level of 10?1 or a slightly below, while ? ( J)2 is at a level of 10?1 or somewhat above. Although such high error levels do not seem to alter the overall dynamics of the solutions qualitatively, they are at a level high enough to be of some concern. 5.3. Accuracy of current layer width

that of the SE run. Again, this accuracy level is qualitatively reasonable. However, if we need to investigate the more di?cult problem of the Parker’s model (3D MICI) and need to determine whether l → 0 or not in the η → 0 limit, a 10% error could lead to signi?cant uncertainty. 5.4. Quantitative summary

One quantity of great interest in the MICI problem is the width of current layers, which can be de?ned as l≡ J 2 / ( J)2

1/2

.

(23)

In Fig. 14, we show the current layer width l de?ned in Eq. (23) for the PS (black), SE (red), and FD (blue) runs, plotted in this order for the Req = 128 case in panel (a), and the error ?l, as compared with the 2048 × 2048 PS run in (b). We see that the error level ?l is similar to that of ? ( J)2 , at about 10?1 or slightly over for the FD run, about an order of magnitude higher than 14

Let us summarize the comparisons from Fig. 7– 13 by plotting the maximum fractional error over the plotted time period for all cases. In Fig. 15, we show in (a) the maximum fractional error over ˙ the plotted period in Fig. 7 to Fig. 9 of ?E (black triangles for PS, red asterisks for SE, blue dia˙ monds for FD) and ?M (black plus signs for PS, red squares for SE, blue crosses for FD) for the three Req cases. In (b), we show the maximum fractional error over the plotted period in Fig. 10 to Fig. 13 of ? A2 (n = 0), ?E (n = 1), ? J 2 (n = 2), and ? ( J)2 (n = 3) versus n, the order of spatial derivative with respect to A (black triangles for PS, red asterisks for SE, blue diamonds for FD), for the Req = 128 case. From these two plots, we see clearly the separation between the three groups (black for PS, red for SE, and blue for FD). While PS has the lowest level of error as expected (due to a fully spectral treatment and to the fact that uniform grids are used), SE achieves an error level somewhat in between PS and FD. In some cases, the error in SE approaches that of PS, indicating near spectral accuracy. So far we have looked at comparisons of spatially integrated quantities. It is conceivable that the accuracy of ?eld quantities at particular spatial points may not follow similar trends. Here we present one example that compares values of the maximum current density, Jmax , over the entire box. In Fig. 16, we show the maximum current density value over the periodic box, Jmax , for the PS (black), SE (red), and FD (blue) runs, plotted in this order for the Req = 128 case in panel (a), and the error ?Jmax , as compared with the 2048 × 2048 PS run in (b). We see in this comparison that even the PS run has an error level at around 10?2 , much higher than the errors of other quantities we have shown so far with this method. This is because the grid used in the 2048 × 2048 run is much ?ner than what is used in the Req = 128 PS run. So, the value of Jmax from a much higher resolution run

Fig. 13.— (a) Time series of ( J)2 for the PS (black), SE (red), and FD (blue) runs plotted in this order for the Req = 128 case. PS and SE data still overlap, but FD the FD result is noticeably di?erent. (b) Fractional error in ? ( J)2 , as compared with a PS run with 2048 × 2048 resolution.

Fig. 15.— (a) Maximum fractional error over the ˙ plotted period in Fig. 7 to Fig. 9 of ?E (black triangles for PS, red asterisks for SE, blue dia˙ monds for FD) and ?M (black plus signs for PS, red squares for SE, blue crosses for FD) for the three Req cases, when compared to a PS run at a grid resolution of 2048x × 2048. (b) Maximum fractional error over the plotted period in Fig. 10 to Fig. 13 of ? A2 (n = 0), ?E (n = 1), ? J 2 (n = 2), and ? ( J)2 (n = 3) versus n, the order of spatial derivative with respect to A (black triangles for PS, red asterisks for SE, blue diamonds for FD), for the Req = 128 case.

Fig. 14.— (a) Time series of the current layer width l for the PS (black), SE (red), and FD (blue) runs plotted in this order for the Req = 128 case. (b) Fractional error in ?l, as compared with a PS run with 2048 × 2048 resolution.

Fig. 16.— (a) Time series of Jmax for the PS (black), SE (red), and FD (blue) runs plotted in this order for the Req = 128 case. (b) Fractional error in ?Jmax , as compared with a PS run with 2048 × 2048 resolution.

15

cannot be located at the collocation points of the lower resolution run, and thus it reaches a slightly higher value. With this consideration, the fact that both the PS and SE runs return ?Jmax to about 1 % accuracy is actually quite good. At the same time, ?Jmax of the FD run is about an order of magnitude higher, of about the same level as ? J 2 . 5.5. Comparison with FD at higher resolution

We must point out that in the above comparison, we require the FD runs to use the same grid as the SE runs, with similar DOF. This is done in order to obtain constrained comparisons that are easier to interpret. However, this will, of course, make the FD runs intrinsically less accurate, as we have so far seen, because FD schemes have a lower order truncation than spectral schemes. In reality, one can compensate for this and increase the accuracy of FD runs by using higher resolution or a greater number of DOF. This will of course require more computational resources. However, with FD schemes usually being more e?cient (e.g. , scalable), this may certainly be considered a reasonable way to obtain more accurate solutions. Here we study brie?y this possibility by running the FD code using higher resolutions. In Fig. 17, the blue curve is the fractional er˙ ror ?E of the FD run as shown in Fig. 7 (a). The ˙ green curve is ?E of a FD run using the same grid as shown in Fig. 3 but with a 16 × 16 uniform grid within each square instead of 8 × 8 as used in the ˙ blue curve. The purple curve is ?E of a FD run using a uniform grid of 256 × 256 resolution. We see that the green curve is mostly at a level below the blue curve, when the linear resolution is doubled. The decrease is sometimes quite high, but at other times just marginal. For the FD run using a 256×256 uniform grid, which has a cell size equals to the smallest cell size of run represented by the green curve, the error is substantially (about an order of magnitude) lower than that of the blue ˙ curve. However, when compared with ?E of the SE run as shown in the red curve in Fig. 7 (a), this is still about one to two orders of magnitude higher. Of course, one can continue increasing the DOF of the FD run to try to reach the accuracy level of the SE run. Following the argu16

Fig. 17.— The blue curve is the fractional error ˙ ?E of the FD run as shown in Fig. 7 (a). The ˙ green curve is ?E of a FD run using the same grid as shown in Fig. 3 but with a 16 × 16 uniform grid within each square instead of 8 × 8 as used in ˙ the blue curve. The purple curve is ?E of a FD run using a uniform grid of 256 × 256 resolution.

ment in (Rosenberg et al. 2006), the linear error bound scaling for the SE case can be written as (Deville et al. 2002, p. 273)

SE

? hmin(p,s) p?s ,

where p is the polynomial expansion order, h ? 1/E is the uniform element length scale, and s is the smoothness of the exact solution. For the moment we neglect prefactors in the error term. We assume s = p so that derivatives can be computed up to the order of the method, suggesting a reasonably smooth function. The error for the FD method is β ?β , FD ? h ? N where N is the total linear (equivalent) number of grid cells, and β is the nominal truncation order (typically ≈ 2). Equating the logarithmic errors yields log N = p log(pE) + prefactor terms. β

This scaling relation provides the number of FD cells required to achieve roughly the error of the SE method at order p. We cannot actually use this scaling directly for comparisons between the FD and SE methods, however, because the prefactor terms are unknown, and may depend critically on the ?ows. Still, we can illustrate the scaling for the case of the SE run in Fig. 7 (a). We have that E = 16 (the equivalently-resolved number of elements) and p = 8, implying that log N ≈ 16/β, so with β = 2, we ?nd that achieving comparable accuracy in the FD method could lead to the requirement for a catastrophically large number of cells, except, perhaps, in the case where the ?lling factor of the ?ne grid is extremely small. This point will require further study. 6. Discussion and Conclusion

which can change e?ciency greatly. Instead, we concentrate on comparing accuracy obtained by the two schemes, mainly to answer the question of whether using SE scheme has any advantage that deserves putting more e?ort into further development. We believe the results of this paper have shown that indeed the SE method has substantial advantages and needs to be investigated further. We have shown that using statically re?ned grids with DOF of roughly linear scaling, SE method can produce results with high accuracy that can sometimes even approach the spectral accuracy of the PS method. This can be potentially very helpful in the investigation of the 3D MICI problem in the ideal limit, which will require adaptive grids that can resolve distinct features. The higher accuracy of the SE method can potentially provide a reliable de?nitive answer to the important question of whether a true current sheet forms in the Parker’s problem of solar and stellar coronal heating. Our main conclusion that SE methods can produce simulations with accuracy somewhat in between PS and FD methods is not surprising in itself. However, our results yield important quantitative information in the context of statically re?ned grids in problems with distinct spatial structures. The accuracy of FD runs presented here is mainly for comparison purposes, since we have imposed the restriction of using the same static grids with the same DOF. This is by no means a suggestion that FD method cannot be used in the investigation of the MICI, or similar problems. Indeed we have also studied brie?y that accuracy in FD method can also be increased by using more DOF. However, we have not studied the trade o? of doing this with respect to the increase of the usage of computational resources. This may be an important topic for future investigation. We acknowledge helpful discussions with Amik St.Cyr at NCAR. Computer time was provided by NCAR and UNH. This research is supported in part by a NSF grant AST-0434322. The NSF grant CMG-0327888 at NCAR also supported this work in part and is gratefully acknowledged.

It is worth noting that at the present stage, the spectral element method described in this paper is more costly in terms of computational time than the FD and PS methods with which we compare our results, the latter being optimal for periodic boundary conditions. We have not tried to compare CPU usage of the di?erent codes due to the fact that the SE and FD codes are running on di?erent computational platforms; furthermore, both codes are still under continuous development, 17

REFERENCES Antiochos, S. K. 1987, ApJ, 312, 886 Balay, S., Buschelman, K., Eijkhout, V., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C., Smith, B. F., & Zhang, H. 2004, PETSc Users manual, ANL-95/11 - Revision 2.1.5, Argonne National Laboratory Basdevant, C, Deville, M., Haldenwang, P., Lacroix, J. M., Ouazzani, J., Peyret, R., Orlandi, P., & Patera, A. 1986, Comp. Fluids, 14, 23 Bhattacharjee, A., Germaschewski, K., & Ng, C. S, 2005, Phys. Plasmas, 12, 042305 Brachet, M., Mininni, P. D., Rosenberg, D., & Pouquet, A. 2007, In preparation Browning, P. K. 1991, Plasma Phys. Controlled Fusion, 33, 539 Canuto, C., Hussaini, M. Y., Quarteroni, A., Zang, T. A. 1988, Spectral methods in ?uid dynamics, (New York, Springer–Verlag) Cowley, S. C., Longcope, D. W., & Sudan, R. N. 1997, Phys. Reports, 283, 227 Deville, M. O., Fischer, P. F., & Mund, E. H. 2002, High-Order Methods for Incompressible Fluid Flow, (Cambridge, Cambridge University Press) Els¨sser, W. M. 1950, Phys. Rev., 79, 183 a Fischer, P. F., Kruse, G. W., & Lottis, F. 2002, J. Sci. Comput., 17, 81 Fischer, P. F. 1997, J. Comp. Phys., 133, 84 Forman, W., Jones, C., Tucker, W. 1985, ApJ, 293, 102 Germaschewski, K., Bhattacharjee, A., Grauer, R. , Keyes, D., & Smith, B. 2004, in Lecture Notes in Computational Sciences and Engineering, Vol. 41, Adaptive Mesh Re?nement—Theory and Applications, ed. T. Plewa, T. Linde, & V. G. Weirs, (Springer, New York), 115 Germaschewski, K., Bhattacharjee, A., &, Ng, C. S. 2006, in ASP Conference Ser. 359, Numerical Modeling of Space Plasma Flows, ed. N. B. Pogorelov and G. P. Zank, 151 18

Grimes, J. P., Heckman, T., D. Strickland, D., & Ptak, A. 2005, ApJ, 628, 187 Kellogg, E., Murray, S., Giacconi, R., Tananbaum, T., & Gursky, H. 1973, ApJ, 185, L13 Klimchuk, J. A. 2006, Sol. Phys., 234, 41 Li, X. S. & Demmel, J. W. 2003, ACM Trans. Mathematical Software, 29, 110 Longcope, D. W., & Strauss, H. R. 1993, Phys. Fluids B, 5, 2858 Longcope, D. W., & Strauss, H. R. 1994, ApJ, 437, 851 Low, B. C 1990, ARA&A, 28, 491 Maday, Y., Patera, A. T., & R?nquist, E. M. 1992, The PN -PN ?2 method for the approximation of the Stokes problem, Pulications du Laboratoire d’Analyse Num?rique, Universit? Pierre et e e Marie Curie Mavriplis, C. 1994, Comput. Meth. Appl. Mech. Engrg., 116, 77 McCormick, S. F. 1989, Multilevel Adaptive Methods for Partial Di?erential Equations, Frontiers in Applied Mathematics, Vol. 6, (SIAM, Philadelphia, PA) Ng. C. S., & Bhattacharjee, A. 1998, Phys. Plasma, 5, 4028 Ng, C. S. & A. Bhattacharjee, A. 2008 ApJ, in press. Parker, E. N. 1972, ApJ, 174, 499 Parker, E. N. 1979, Cosmical Magnetic Fields, (Clarendon Press, Oxford) Parker, E. N. 1994, Spontaneous Current Sheets in Magnetic Fields: with Applications to Stellar X-Rays, (Oxford University Press, New York) Rosenberg, D., Pouquet, A., & Mininni, P. D. 2007, New J. Phys., 9, 304 Rosenberg, D., Fournier, A., Fischer, P., Pouquet, A. 2006, J. Comp Phys., 215, 59 Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1

St-Cyr, A., Jablonowski, C., Dennis, J. M., Tufo, H. M., Thomas, S. J. 2007, Mon. Weather Rev., in press Strauss , H. R. 1976, Phys. Fluids, 19, 134 van Ballegooijen, A. A. 1985, ApJ, 298, 421 Zweibel, E. G., & Li, H.-S. 1987, ApJ, 312, 423

A This 2-column preprint was prepared with the AAS L TEX macros v5.2.

19