Adaptive mesh refinement with spectral accuracy for magnetohydrodynamics in two space dimensions

We examine the effect of accuracy of high-order spectral element methods, with or without adaptive mesh refinement (AMR), in the context of a classical configuration of magnetic reconnection in two space dimensions, the so-called Orszag–Tang (OT) vortex made up of a magnetic X-point centred on a stagnation point of the velocity. A recently developed spectral-element adaptive refinement incompressible magnetohydrodynamic (MHD) code is applied to simulate this problem. The MHD solver is explicit, and uses the Elsässer formulation on high-order elements. It automatically takes advantage of the adaptive grid mechanics that have been described elsewhere in the fluid context (Rosenberg et al 2006 J. Comput. Phys. 215 59–80); the code allows both statically refined and dynamically refined grids. Tests of the algorithm using analytic solutions are described, and comparisons of the OT solutions with pseudo-spectral computations are performed. We demonstrate for moderate Reynolds numbers that the algorithms using both static and refined grids reproduce the pseudo-spectral solutions quite well. We show that low-order truncation—even with a comparable number of global degrees of freedom—fails to correctly model some strong (sup-norm) quantities in this problem, even though it satisfies adequately the weak (integrated) balance diagnostics.


Introduction
In geophysical and astrophysical flows, the Reynolds numbers are large, and nonlinear terms appearing in the primitive equations lead to strong mode coupling and multiple scale interactions. Moreover, magnetic fields are often in quasi-equipartition with the velocity, as is the case in the Solar wind (with in fact a slight excess of magnetic energy [1,2]), or in the interstellar medium. However, direct numerical simulations (DNS) using regular grids cannot, even to this day, deal with such large Reynolds numbers R v . Doubling the grid resolution (and thus multiplying the Reynolds number by roughly a factor two) comes at a cost of increase in needed computer time by a factor of sixteen in three-dimensions (3D), assuming the temporal scheme is explicit; even when taking Moore's law into account, such an increase in R v can only be achieved roughly every six years. Thus, one is led to resort to more sophisticated techniques, such as, for example, turbulence modelling (see e.g. [3]). However, in the case of coupling to a magnetic field, and using the magnetohydrodynamic (MHD) approximation valid for the description of the largescale dynamics, few such techniques have been developed and tested (see however [4,7]).Another venue is to develop adaptive mesh refinement (AMR) methods. In this context, we examine in this paper the accuracy of an AMR code using spectral elements by comparing its output to exact solutions in simplified cases and to computations using a pseudo-spectral code at the same Reynolds numbers on a classical configuration of magnetic reconnection in 2D geometry, the Orszag-Tang (OT) vortex [8].
We set-up the equations in the next section 2, and in 3, we describe the numerical method for MHD developed within the context of the adaptive spectral element method presented in [9] for the 2D Burgers equation. This section also presents test results for the method. In section 4, we apply the method to the OT configuration, and compare it with pseudo-spectral results. We consider effects of low-order versus high-order local approximations in section 5, and section 6 is the conclusion, in which we summarize the results, and offer some observations about the performance of the method and some directions for future work.

Equations, code and simulations
For an incompressible fluid with constant mass density ρ 0 , the MHD equations read: where u and b are the velocity and magnetic field (in Alfvén velocity units, b = B/ √ µ 0 ρ 0 with B the induction and µ 0 the permeability); P is the pressure divided by the mass density, and ν and η are the kinematic viscosity and the magnetic resistivity (note that, in this paper, we adopt the following notation in order to be consistent with [9] and references therein: we use vector signs for continuous vector fields, boldface for the collocated (discretized) quantities, and sans serif bold for discrete operators). The mode with the largest wavevector in the Fourier transform of u at initial time is k 0 . We also define the viscous dissipation wavenumber as k ν = ( /ν 3 ) 1/4 , where ∼ U 3 0 /L 0 is the energy dissipation rate, with U 0 the rms velocity and L 0 the integral length scale (see below). The Kolmogorov scale at which dissipation sets in is defined as l D = 2π/k ν ; the expression for k ν is based on a kinetic energy spectrum E V (k) ∼ ε 2/3 k −5/3 . A large separation between the two scales (k −1 0 k −1 ν ) is required for the flow to reach a turbulent state with significant nonlinear interactions.
In the absence of external forcing, viscosity and magnetic resistivity, the MHD equations in two space dimensions (2D) conserve the total energy: together with the cross helicity H c = 1 2 u · b dxdy , and the L 2 norm of the magnetic potential M a = 1 2 a 2 dxdy , with b = ∇ × a. The Reynolds number is defined as R v = U 0 L 0 /ν, where the integral length scale of the flow is defined as The large scale turnover time can then be defined as τ NL = U 0 /L 0 . We can also introduce the Taylor based Reynolds number R λ = Uλ/ν, where the Taylor length scale λ is given by Length scales built with b and its energy spectrum E M (k) can also be defined in a similar fashion; for example, the magnetic Reynolds number is R m = U 0 L 0 /η.

Algorithm description for MHD
For this work, we use a spectral element method [10], encapsulated within the Geophysical /Astrophysical Spectral-elementAdaptive Refinement (GASpAR) code.Aspects of this code-in particular those regarding the dynamic grid refinement-have been described elsewhere [9], where results were presented for the multi-dimensional Burgers (advection-diffusion) equation.
Here, we extend that development in several distinct ways in order to solve (1)-(3), specifically, by adding the pressure term and the Lorentz force in the momentum equation and by taking into account the magnetic induction equation, and the divergence-free conditions on velocity and magnetic field. We solve equations (1)-(3) in Elsässer form [11]: where Thus, we solve for u and b in terms of Z ± . Equations (7) and (8) can be recast into an equivalent d-dimensional variational form on domain D by defining the following function spaces: where w = u, b. Let Z ± and P and their test functions, ζ ± and q be restricted to finitedimensional subspaces of these spaces: The basis for the velocity expansion in P N is the set of Lagrange interpolating polynomials on the Gauss-Lobatto-Legendre (GL) quadrature nodes, and the basis for the pressure is the set of Lagrange interpolants on the Gauss-Legendre (G) quadrature nodes. Then, the equations (7) and (8) can be written in weak form as [12]: where C ± := Z ± · ∇ is the continuous advection operator, ·, · GL , represents the inner product using quadrature on the GL nodes, and ·, · G indicates the inner product using quadrature on the G nodes, as in [9] appendix A. Thus, we use a staggered grid, where the quantities ( u, b, Z ± ), on the GL nodes are continuous, while those on the G nodes (P) are not. This so-called P N − P N−2 formulation was chosen to prevent spurious pressure modes [12,13].
In the spectral element method, functions in U N and Y N−2 are represented as expansions in terms of tensor products of basis functions within each subdomain, or element, E k , the non-overlapping union of which comprises the domain [10]: D = K k=1 E k . By expanding Z ± and P in terms of their basis functions, inserting these expansions into (15), and using the appropriate quadrature rules, we arrive at a set of semi-discrete equations in terms of spectral element operators: for the jth components of momentum. In this equation, M, L and C are the well-known mass matrix, weak Laplacian and advection operators, respectively (e.g. [9], and references therein).
The variables Z ± represent values of the Z ± collocated at the GL node points, while P ± are values of the pressures at the G node points. Similarly, we denote the collocated test functions in Y N−2 by q. The Stokes operators, D j , arise from the quadrature rule in (15), in which the GL basis function and its derivative must be interpolated to the G node points, and multiplied by the G quadrature weights. In 2D, . For regular rectangular elements, of lengths L k j , are, respectively, the weighted 1D interpolation matrix mapping GL points to the G points, and the weighted 1D GL derivative matrix interpolated to the G points [12]. In these definitions σ i are the weights corresponding to the G node points. Just like with the mass and Laplacian operators, the Stokes operators are written as tensor products of 1D operators. The above equations (16) and (17), are correct for a single element, but they are not complete when we have multiple subdomains. In this case, we must ensure that all quantities in U N remain continuous across element interfaces. The manner in which this is done for advection-diffusion on (non)conforming elements was described in [9]. Using the Boolean scatter matrix, A c , the interpolation matrix from global to local degrees of freedom (d.o.f.), , and the masking matrix (that enforces homogeneous boundary conditions), , that were presented there, we find that we must replace the local Stokes operators above with where D L,j = diag k (D k j ), and the D k j are the matrices from above. In much of what follows, we will continue to use the local form of the Stokes operators, and simply recognize that the multiple subdomain form can be imposed.
Note in (16) the presence of different pressures for Z ± . As we show below, we will maintain the divergence constraints by solving for the pressures for both Z ± . While analytically these pressures are the same, they serve as Lagrange multipliers [14] for their respective fields, Z ± . Given that each field has its own constraint that is solved independently of the other, numerically they will in general be different.

Time stepping
Because we want to resolve all time scales (as well as spatial scales), we choose an explicit method for integrating (16) and (17) in time. The method we use is the mth-order Runge-Kutta (RK), and since the right-hand sides of the equations have no explicit dependence on time, we can use the following formulation [15, p 109] in order to solve dU/dt = F(U).
Considering one component of the Elsässer variables and following the above RK algorithm, we write for each iteration (recall (16)) We insist that each RK stage obeys (17) in its discrete form, so multiplying (19) by D j , summing, and setting the term D j Z ± j = 0 leads to the following pseudo-Poisson equation for the pressures, P ± : where the remaining inhomogeneity The pseudo-Laplacian operator on the left-hand side of (20) also arises from a second-order implicit time descretization of the spatially discretized equations (16) and (17) as shown in [12,16]. In this formulation, even if Z ± j is not divergence-free, the partial update after this RK stage will be. The inverse mass operator, M −1 , must be computed for a given grid configuration. For conforming elements, this matrix is trivially inverted since it is diagonal. For nonconforming discretizations, the mass matrix is lumped in order to recover a diagonal matrix that can be inverted straightforwardly [16]. Equation (20) is solved using a preconditioned conjugate gradient (PCG) method [17,18] (see [9] for modifications required for nonconforming elements). For preconditioning, we use a block Jacobi preconditioner computed using either a fast diagonalization method [19], or a direct inversion.
Thus, at each time step, m RK stages are computed, and each stage solves (20) twice, once for Z + , and once for Z − leading to 2m pseudo-Poisson solutions at each time step. The MHD solver is encapsulated within a class as are all solvers in the code.
Currently the solver considers the Elsässer variables, Z ± , to be auxiliary variables, so a transformation is done from the native u, b, and back. However, it is easy to create a new constructor for the solver object, so that Z ± are the primary quantities of interest. In all results presented in this work, we choose m = 2 for the RK order.

Tests with exact solutions
To test the algorithm described above, we first consider two test problems that have analytic solutions. The first case tests the algorithm without a magnetic field ( b ≡ 0)-i.e. solves the Navier-Stokes equations-and the second case tests the full MHD algorithm. Throughout the remainder of the paper, the symbol p will refer to the polynomial order in each coordinate direction within an element.
For Navier-Stokes, we use a steady state solution of a Kovasznay [20] flow field behind a periodic array of cylinders (see also [19]): where We initialize the grid with the solution, and march to a steady state (typically to t ≈ 2), where we compare the numerical solution with the analytic solution. For the following test, the grid we use is conforming, and non-adaptive. The computational grid has lower left and upper right vertices (−0.5, −0.5), (1, 1.5), and is comprised of 2 × 4 elements, with higher resolution near the array on the inflow axis. The time step is fixed at 1 × 10 −3 , while R v = 1/ν = 40. Dirichlet boundary conditions are set from the analytic solution. The left panel of figure 1, which shows streamlines = cst over the entire grid, where u = ∇ × ẑ, illustrates the global grid and its elemental decomposition. In the right panel of this figure, we plot the L 2 norm of the error as a function of polynomial degree, p, in order to show convergence. Clearly, the solution is converging spectrally as desired.
The next test looks at a steady Hartmann flow, consisting of a flow of viscous conducting fluid between two parallel plates, with separation 2a, and with a magnetic field, B 0 , applied in a direction perpendicular to the plates. The solution is [21]: where the Hartmann number, and δ represents a boundary layer thickness over which the velocity compared with that in the central region decreases at the top and bottom plates.
Again, we initialize a grid of K = 4 × 2-elements-whose lower left and upper right global grid vertices are (0, 0), and (4, 2)-with the analytic solution, at fixed p, and we march to a steadystate, comparing the numerical and analytic solutions. For these runs, we use a Courant-limited time step. Dirichlet boundary conditions are again set from the analytic solution. To drive the flow, we add a constant x-momentum forcing, γ x , to the right-hand side of (1) (or, equivalently, to (7)) such that For all our tests, we set H a = 4, R v ≡ u 0 a ν = 40, u 0 = 1, and B 0 = 1; the choice of H a > 1 suggests a flow where the magnetic field is dynamically significant. In figure 2 we present our convergence results. As expected, again we see spectral convergence.

The OT vortex
The OT vortex ( [8]) is a geometrically simple initial configuration with a magnetic X-point centred at a stagnation point of the velocity. It is best expressed in terms of the stream function (with u = ∇ × ẑ) and the magnetic potential. It reads = 2α[cos(2πx) + cos(2πy)], a z = 2 cos(2πx) + cos(2πy).
The total energy E T = E V + E M , for α = 1, is equally divided initially between its kinetic and magnetic components E V and E M , both equal to 2, and the initial correlation coefficientρ c is equal to 41%, whereρ c = 2 u· b u 2 +b 2 . This configuration is known to develop several current sheets in a time of order unity in these units; such current sheets further destabilize in time through a tearing mode instability embedded within a turbulent flow and develop numerous complex small-scale structures (see e.g. [22]).
In this section, we apply the algorithm described in section 3 to the OT problem. For adaptive runs, a spectral estimator refinement criterion is used [9,23,24] that estimates the solution error. If, in a given element, this error is greater than a specified tolerance, ε est , the element is tagged for refinement. The nominal resolution of the adaptive runs is given by an equivalent resolution, which is computed by where N 0 is the initial number of elements in either direction (the base grid), and max is the maximum refinement level. All computations are performed on a periodic grid of dimension [0, 1] 2 . We compare the spectral element solutions with those obtained from a well-characterized pseudo-spectral code that has been used to produce numerous results cited in the literature [25,26].
The total dissipation in the flow is defined as where ω = ∇ × u is the vorticity and j = ∇ × b is the current density. It is a global quantity, as is the total energy E T , and is characteristic of the dynamical evolution of the flow as a whole: as small-scale gradients in both the velocity and the magnetic field develop through nonlinear interactions, current and vorticity sheets form and dissipation sets in at a time of order unity. The temporal evolution of E T and D T are shown in figure 3 for the GASpAR MHD run and the pseudospectral code for a fiducial run with R v = 10 3 π; the initial grid is K = 8 × 8 with max = 3 and p = 8 (implying N eq = 512), and ε est = 1 × 10 −5 . The adaptive code reproduces the temporal characteristic times (as well as secondary maxima in D T at a lower Reynolds number, shown in figure 10 in the context of a fixed grid, see below); it also reproduces the amplitude of the global phenomenon of reconnection of magnetic field lines and ensuing dissipation of energy. The total number of d.o.f. normalized by the number of modes in the pseudo-spectral run used in this comparison during the temporal evolution of the flow is also shown in figure 3 in black. The d.o.f. quickly increases because the initial grid was coarse (K = 8 × 8 elements, at p = 8). After this, the d.o.f. grows regularly in anticipation of the peak in total dissipation at a slightly later time; note that it is about one-half that of the pseudo-spectral computation. Other runs at lower Reynolds numbers indicate that, beyond this peak, the number of d.o.f. is roughly constant until about t = 1.8, where the grid appears to anticipate a secondary peak in D T that begins at about t = 2.5, and corresponds to renewed reconnection of current layers [22] (not shown).
As the relative number of d.o.f. increases toward unity, the AMR becomes less efficient. In this run, the grid refines on the native solver variables u and b, and by the time the nonlinear regime begins, the small scales dominate the flow requiring finer elements to remain on the grid to resolve them. The number of d.o.f. also depends on the criterion of refinement and on ε est ; note that here ε est is set tightly so that the grid is more likely to refine. It will be left to future work to investigate this and other refinement criteria more systematically, in order to determine what range of parameters, and, indeed, which criteria, are the most robust and efficient for adaptive modelling of flows such as OT. We have freedom in the code also to define new variables to which to apply our refinement criteria, providing yet another avenue for investigation (see also the discussion below around figure 6). As an illustration, we also show in figure 3 the number of d.o.f. for a spectral criterion of refinement, based this time on the vorticity and the current. In this case, the refinement of the grid starts at a later time, once the strong gradients have developed, and thus the run is roughly twice as fast; however, once we approach the saturation of the growth of small scale production, the number of d.o.f. for the two refinement criteria become comparable.
Conservation of energy (and that of the other quadratic invariants) is an essential feature in the detailed dynamical evolution of a turbulent flow; it is at the foundation of the concept of energy (or invariant) cascade and leads, assuming a constant energy dissipation rate within the cascade, to the celebrated Kolmogorov energy spectrum E V (k) ∼ k −5/3 . Note that in MHD, the power law followed by the energy spectrum in the inertial range is less clear, and other spectra can be postulated a priori on the basis of Alfvén wave propagation, including an anisotropic component of the spectrum linked with the bi-dimensionalization of the flow in the presence of a strong uniform magnetic field. Such power laws are barely observable, due to a variety of reasons. In the laboratory or in observations of geophysical flows, the instrumentation has cutoff frequencies, and in numerical simulations, the resolution is barely sufficient to be able to resolve the ranges necessary to follow the evolution of the flow, namely the energy-containing range The black curves that represent the spectra computed from the pseudo-spectral method, and the red curves obtained with GASpAR cannot be distinguished except at the very highest wavenumbers. The upper curves are the magnetic energy, while the lower are the kinetic energy. The line E(k) ∝ k −5/3 corresponding to a Kolmogorov spectrum is given as a reference.
around k ∼ k 0 , the inertial range where loss-less transfer of energy occurs, and the dissipation range. For example, for fluids in 3D, it was shown in [27] that the Kolmogorov energy spectrum, on resolutions of up to 1024 3 points, occurs on a small range of wavenumbers (slightly more than two octaves) and is followed by a shallower power law, named a bottleneck effect.
The energy spectra close to the maximum of enstrophy, at t ≈ 0.85, are given in figure 4; the spectra for the AMR run are computed on an irregular grid using the algorithm derived in [28]. We see that the agreement is quite good between the adaptive spectral element run and the pseudo-spectral run. It is interesting to note that, while the pseudo-spectral case uses dealiazing, no explicit dealiazing is required in the GASpAR run.
A further test of the code is to verify that the nonlinear terms of the primitive equations do preserve the invariants, as spectral methods are expected to be conservative. This can be done by computing the time derivative of the invariant (say, using an algorithm of order 2), and comparing it to its theoretical value. In the case of the total energy, for example, the dissipation D T given in equation (21) is the theoretical value. In figure 5 we show the degree to which the energy is preserved and observe again an excellent agreement.
A snapshot of the current density shown in figure 6 near the peak in total dissipation indicates that, as expected, the grid refines in the region of the current sheets (and of vorticity sheets, which are known to be almost co-located). However, when compared with runs from other authors [29] using finite difference methods and configuration space refinement on the gradients of the basic variables, there appears in our run to be more refinement than is initially expected outside of the  sheets. In the case of the spectral refinement criterion applied to u and b, which are native to the solver, the criterion appears to capture the variation in the curvature of u and b near the current sheets to an extent that may not be required since, with the same ε est , half the d.o.f. are needed at early times with the spectral criterion based on j and ω (see figure 3 and subsequent text). , and on j and ω (blue). The green curve is for a run that also uses the spectral error criterion on j and ω, but has p = 16. The black (lower) curve is for the pseudo-spectral run.
We now examine the temporal evolution of the current maximum J max in this fiducial run, as shown in figure 7. As is well known, the development of singularities in such flows may be diagnosed by a 1/[t * − t] behaviour [30], with t * being the time of singularity. As such, the temporal evolution of extrema is an important feature of turbulent flows; from a physical point of view, such extreme events are also the locus of anomalous diffusion and concentration of (say, chemical) tracers so that a faithful reproduction of such events may be required. The linear phase corresponding to an exponential growth of the current instability is well reproduced but as we enter the nonlinear phase, studied in the literature in the context of the nonlinear growth of the tearing mode instability [31], errors appear when using the AMR codes, errors that are comparable for the two refinement criteria just described (see figure 7). As already shown in 3D [32], the nonlinear phase is highly nonlocal, which indicates that many scales are interacting in the development of current instabilities. This numerical discrepancy (not observable in the L 2 norm) might be remedied by tightening the refinement criteria. However, the accuracy of the code, as measured primarily by the polynomial order in each element for a fixed number of elements and without refining the grid, is a parameter that must also be examined, and this is done in the next section.

High versus low-order
We now consider the behaviour of the OT solutions when the number of global d.o.f. is kept (roughly) constant, while the polynomial truncation (degree) varies. In the first series of runs, we have ν = η = 0.025 (R v = 80π), and these runs are compared with a pseudo-spectral run with a resolution of 128 2 grid points. Here, as stated before, we use a static conforming mesh Table 1. Parameters used in the simulations described in section 5; p is the polynomial order in each direction for each element, N eq is the linear grid resolution, and E is the number of elements in each coordinate direction such that the total number of elements is K = E × E. for each run so that the results will not be affected by refinement and coarsening criteria, since we want to focus on the effect the order of the method has on the results. In table 1, we present the relevant run parameters. Figure 8 shows the profile of J max as a function of time for each of the N eq ≈ 128 runs. It is seen immediately that the low-order truncation does not yield accurate values for J max , particularly as the peak in total dissipation is reached near t = 1.1. Clearly, J max converges to the correct solution as p increases (compare the red curve for p = 5 and the black curve for the pseudo-spectral run). The maximum of vorticity behaves in much the same way as J max .
At this stage, it is desirable to compare the accuracy for a pseudo-spectral code with N points per linear dimension, ps , to the error, se , of a spectral element code with E elements per linear dimension, each with polynomials of order p. Omitting prefactors with slower variations in the truncation orders, we have for the pseudo-spectral error [15, p 400] that ps ∼ x N ∼ 1/N N , and for the spectral element error bound [19, p 273], that where h ∼ 1/E is the uniform element length, and s is the smoothness of the exact solution. It is clear that in practice s p, since the derivatives for s > p cannot be computed. We thus choose s = p so that the function is sufficiently smooth to allow spectral convergence.
Equating the logarithmic errors immediately shows the relationship between N, E and p, namely: The above scaling argument allows for choosing a range of values of polynomial order under simple assumptions. Let us say the Reynolds number is doubled from a well-adjusted run; roughly speaking, we need to double the resolution of the pseudo-spectral code from N to 2N grid points per linear dimension. Let us also assume that we double the number of elements in the spectral element code. Then, under the reasonable assumption of large N, we find that the polynomial order needs to be doubled as well. Though empirical, this criterion does indicate that the polynomial order needs to increase with the Reynolds number. It should be noted that, whereas the pseudo-spectral code, in the preceding example of figure 8, uses a truncation of order 128, the equivalent spectral element code uses polynomials of order 5, substantially smaller, in order to achieve comparable accuracy.  (see table 1). For each case, ν = η = 0.025. The black curve is the pseudo-spectral case; next to it, the red curve is the p = 5 run, followed by dark blue (p = 4) and green (p = 3). Note that the p = 5 and pseudo-spectral results are nearly coincident.
As a specific example, an examination of figure 8 indicates that, with E = 26, setting p = 5 leads to a satisfactory computation of J max . Let us now double the resolution of the pseudo-spectral code to N = 256 and take E = 52. As can be seen in figure 9, it is indeed the case that p 8 gives an accurate representation of the dynamics of the flow at the enhanced Reynolds number. According to the scaling relationship, if we were to retain p = 5, we would need E ≈ 10 000 in order to reach the same level of accuracy at the higher R v . Applying these scaling relations to the fiducial adaptive run in section 4, by bootstrapping from N eq = 256 with p = 8 to N eq = 512, suggests that we set p ∼ 16 in order to see comparable accuracy to the pseudo-spectral results. Indeed, as we see in figure 7, J max is better behaved with the higher-order elements, confirming our scaling relation in the case of adaption as well.
For completeness, we show in figure 10 the corresponding L 2 norms for the total energy (top) and the total generalized enstrophy ω 2 + j 2 for the N eq ≈ 256 cases; the different runs cannot be discerned.

Discussion and conclusion
We have presented an explicit spectral element method for solving the equations of incompressible MHDs. The method is developed within the context of an existing spectral element code (GASpAR) that provides many of the spectral element operators required of the MHD algorithm, and that also offers an adaptive, nonconforming mesh algorithm. The new operators that arise in the explicit MHD treatment, have been defined. Here, we have described tests that compare the numerical results with analytic solutions and establish that the method achieves spectral convergence in the case of conforming elements (for preliminary tests on a static reconnection problem, see [33]). We have then applied the method to a challenging problem in the literature, the so-called OT flow, which allows for magnetic reconnection and the development of current sheets. The OT runs are compared with well-tested pseudo-spectral solutions as a baseline, and found to agree well. We find that the quadratic L 2 diagnostic quantities are insensitive to variations in polynomial degree, while the sup-norm quantities, such as the maximum current density, are not accurate at low-order truncations. Because such sup-norm quantities are the foundation for criteria of the development (or not) of singularities in Navier-Stokes and MHD flows [30,34], it may be of some importance to be able to solve for them accurately.
It is worth pointing out that the spectral element method described in this paper is more costly in terms of computational time than the pseudo-spectral method with which we compare our results, the latter being optimal for periodic boundary conditions. This is despite the fact that the spectral element method requires only nearest-neighbour communication, while the pseudospectral method requires all-to-all exchange of data in the fast Fourier transform algorithm. This performance issue can be traced directly to the solutions of the pseudo-Poisson equation (20) that are required in order to maintain the divergence constraints. The pseudo-Laplacian operator is known to be ill-conditioned, primarily because of a 1D null space. Indeed, for the N eq = 512 2 equivalent runs, we see typical iteration counts for each RK stage of ∼250, which further increase-albeit reasonably slowly-in the OT runs as the enstrophy maximum is approached; we see no significant reduction in the PCG iteration count when we attempt to remove this null-space, but more work is needed in this area. (It should be noted that the scaling of the spectral element code on multiprocessors is very good, but we refer here to a single processor performance.) We remark that other formulations can be used to solve the MHD equations with the spectral-element method. In particular, in 2D, the vorticity-stream function-vector potential formulation produces manifestly solenoidal fields. However, a Poisson-type equation has still to be solved to obtain the stream function from the vorticity. Other formulations that can be readily extended to 3D (such as the Elsässer variables used here, or the velocity-magnetic field formulation) also yield Poisson-type equations for incompressible flows. Although we expect the ill-conditioning of these spectral element Laplacian operators to be similar to the pseudo-Laplacian described here, different formulations will be explored in the future. Furthermore, when dealing with incompressible MHD fluids at low magnetic Prandtl number, as encountered in the laboratory and in the liquid core of the earth, there is a need for accurate simulations of the generation of magnetic fields in turbulent flows with complex boundaries, in conjunction with several ongoing experiments (see e.g., [35]- [38]). Spectral element codes, which easily encompass a variety of boundary conditions and geometries may be useful in this context.
There are a number of ways to speed up the spectral element solutions. One is to implement a more sophisticated preconditioner, and we are making progress in this regard that will be reported elsewhere. Another is to relax the degree to which the divergence constraints are maintained, by increasing the PCG convergence tolerance to a less stringent value. Still another, as alluded to in section 4, is to find refinement criteria that optimize the number of d.o.f. for a desired quadrature or truncation error. We have seen that the adaptive mesh code can provide a substantial savings in work because the number of d.o.f., as shown in several instances (e.g. [9]), can be reduced. The code is presently being extended to include compressibility in view of the many MHD applications we have in mind in the astrophysical context (solar wind, magnetosphere, solar convection zone and corona), in which case the performance problems associated with this divergence-free constraint should be alleviated.