A constrainedtransport magnetohydrodynamics algorithm with nearspectral resolution.
Abstract
Numerical simulations including magnetic fields have become important in many fields of astrophysics. Evolution of magnetic fields by the constrained transport algorithm preserves magnetic divergence to machine precision, and thus represents one preferred method for the inclusion of magnetic fields in simulations. We show that constrained transport can be implemented with volumecentered fields and hyperresistivity on a highorder finite difference stencil. Additionally, the finitedifference coefficients can be tuned to enhance highwavenumber resolution. Similar techniques can be used for the interpolations required for dealiasing corrections at high wavenumber. Together, these measures yield an algorithm with a wavenumber resolution that approaches the theoretical maximum achieved by spectral algorithms. Because this algorithm uses finite differences instead of fast Fourier transforms, it runs faster and isn’t restricted to periodic boundary conditions. Also, since the finite differences are spatially local, this algorithm is easily scalable to thousands of processors. We demonstrate that, for lowMachnumber turbulence, the results agree well with a highorder, nonconstrainedtransport scheme with Poisson divergence cleaning.
1 Introduction
Many astrophysical flows involve dynamically significant magnetic fields, such as molecular clouds, accretion disks, the galactic dynamo, jets, galaxy clusters, stellar dynamos and coronae, the solar wind and the interstellar medium. These problems tend to be threedimensional, multiscale and turbulent, so there is an ongoing interest in developing highresolution and efficient magnetohydrodynamics (MHD) algorithms for them. In this paper, we outline an extension of the constrained transport algorithm (Evans & Hawley 1988) to the combination of higher spatial order and zonecentered grids, and with resolutionenhanced tuned derivatives. We then describe how these measures fit together to yield an algorithm that closely approaches the theoretical maximum wavenumber resolution of spectral algorithms.
1.1 Constrained transport
The induction equation for a magnetic field and a velocity field in ideal MHD is
(1) 
Analytically, this equation conserves magnetic divergence: . However, this may or may not be the case for a finitedifference treatment of this equation. Tóth (2000) reviews the methods taken by various algorithms to treat the divergence in MHD simulations. A spectral code explicitly projects the Fourier components so that . For a finite difference code, the magnetic field can be evolved by a constrained transport scheme that preserves the magnetic divergence to machine precision (Evans & Hawley 1988). Alternatively, if the discretization doesn’t conserve magnetic divergence, the divergence can be removed with measures such as periodic use of a Poisson solver (Brackbill & Barnes 1980), adding a divergence diffusion term to the magnetic evolution, or following an artificial and independently evolving divergence field (Dedner et al. 2002) to propagate divergence away from where it is produced and then dissipate it. The Powell scheme (Powell et al. 1999) adds a source term to advect divergence rather than let it grow in place. A finite difference code can also employ a vector potential such that , in which case the magnetic divergence is automatically zero. This requires the use of a higherorder advection algorithm to ensure accurate secondderivatives, as is done in the Pencil code (Brandenburg & Dobler 2002).
We denote any finite difference scheme for MHD that explicitly conserves the magnetic divergence to machine precision as constrained transport (CT), and any scheme that does not as unconstrained transport (UT). Several variations of CT are possible. If the electric field is differenced as a curl: then the magnetic divergence is preserved to machine precision for most grid types (see Appendix). Evans & Hawley (1988) introduced CT for staggered grids and Tóth (2000) showed that it works for centered grids as well (see the Appendix for explanation of centered and staggered grids). Londrillo & Del Zanna (2000) further showed that highorder CT is possible on staggered grids with a radiustwo stencil. In this paper, we show that volumecentered CT is possible on arbitrarily large stencil sizes, with hyperresistivity, and that the resolution of this algorithm at moderately high order approaches the theoretical maximum exhibited by a spectral code.
2 Algorithms
Our algorithm is based on a constrained transport scheme, plus measures to enhance the resolution and maintain stability. High wavenumber resolution is achieved by a combination of highorder and tuned finite differences plus hyperdiffusivity, and stability is achieved by RungeKutta timestepping and hyperdiffusivity.
2.1 Timestepping
A highorder timestepping scheme for the evolution equations is essential for the stability of most algorithms. The time update for a variable is , where represents some estimate of . One example is the secondorder RungeKutta scheme, which estimates and then identifies . Another class of algorithms maintains the conservation of mass and momentum by computing fluxes through zone boundaries. A variety of techniques exist for timeextrapolating the fluxes at , such as piecewise parabolic advection (Colella and Woodward 1984), Total Variation Diminishing (Harten 1983), Riemann solvers (Toro 1999), the method of characteristics (Stone & Norman 1992a, Hawley & Stone 1995), and many more. MHD poses a challenge to time extrapolation because there are seven or eight wavemode characteristics, depending on the technique used for treating magnetic divergence. In particular, the wellknown methodofcharacteristics algorithm interpolates along the Alfvén characteristic while neglecting the fast and slow mode characteristics. For our simulations, we use RungeKutta for timeextrapolation because it doesn’t invoke any diffusive spatial interpolations (§ 2.2), and because it automatically captures all three MHD wavemode types. In our demonstration implementation, we use secondorder RungeKutta while the Pencil code (Brandenburg & Dobler 2002) uses thirdorder, although either order has proven successful.
2.2 Diffusivity
A common class of algorithms is based on momentum fluxes that are timeextrapolated with upwind spatial interpolations. The errors from the interpolations required for these flux transport algorithms produce an intrinsic diffusivity that can stabilize the evolution, even in the absence of any explicit diffusive terms. The nature and magnitude of the diffusivity has been characterized in Zhong (1998) and Dobler et al (2006).
RungeKutta timestepping, on the other hand, has no spatial interpolations, and thus no intrinsic diffusivity. One then generally needs an explicit stabilizing diffusivity. One has various options for the form of this diffusivity, with Laplacian or hyperLaplacian typically chosen. These diffusivities have the benefit that their magnitude is easily characterized, and the diffusive coefficient can be tuned to have the minimum value necessary to preserve stability. Consider:
(2) 
Let the Fourier components be where is the wavenumber. They evolve according to
(3) 
The term is the Laplacian viscosity and the others are higherorder hyperdiffusivities. Specifically, and They differ in that the higher the order, the more selective they are in diffusing highk structure while preserving lowk structure. Equivalently stated, highorder hyperdiffusivities erase smallscale structure without affecting the larger scales. For many physical applications, such as the turbulent magnetic dynamo, the large scale structure is unaffected by the nature of the smallscale diffusivity (Haugen & Brandenburg 2004), and so in these instances, the use of hyperdiffusivity instead of Laplacian diffusivity enhances the wavenumber resolution.
The operator is spherically symmetric in Fourier space while the operator is not. This affects the maximumpossible timestep because in order to be advectively stable, the Courant condition implies that the product must be less than a given value, and so the highk corners of the 3D Fourier cube are the most vulnerable to advective instability. In these corners, the term delivers more diffusion than the term, but with a cost of twice as many finite difference operations.
Hyperdiffusivity acts together with highorder finite differences to enhance the resolution of a simulation. Highorder finite differences allows structure to be finite differenced with less error, and hyperdiffusivity allows this structure to evolve with less dissipation than it would with Laplacian diffusivity.
Hyperdiffusivity can benefit the timestep as well as the resolution. Suppose one chooses the Laplacian term to provide the dominant fraction of the diffusivity. One typically then evaluates by trial and error the minimum value of and the maximum timestep that one can get away with to maintain stability. Once these values are chosen, the addition of a small measure of hyperdiffusivity, small enough so as not to contribute significantly to the total diffusivity, but large enough to affect the highestk structure, increases the maximum stable timestep by about 50%. This technique works because the highestwavenumber structure is the most vulnerable to instability, and so the timestep depends most critically on the value of the diffusivity for these wavenumbers.
One could additionally note that “maximum stable timestep” is not sharply defined. For instance, a timestep at the cusp of stability might be unstable, but only after perhaps more than 10 timesteps. Also, the maximum stable timestep depends on the maximum value of the velocity in the simulation, which is changing in time.
2.3 Magnetohydrodynamic equations
The equations of incompressible MHD with diffusivity and hyperdiffusivity are
(4) 
(5) 
(6) 
V  Velocity  B  Magnetic field 


J  Current ()  

Laplacian Viscosity  Laplacian resistivity  

Hyperviscosity for  Hyperresistivity for  

Hyperviscosity for  Hyperresistivity for  

Hyperviscosity for  Hyperresistivity for  

Hyperviscosity for  Hyperresistivity for  

Hyperviscosity for  Hyperresistivity for  

Divergence viscosity  Divergence resistivity 
.
The variables are defined in Table (1). The magnetic equation (5) is written as a curl so that the contrained transport scheme can preserve the magnetic divergence (Appendix). The remaining term proportional to outside of the curl has the effect of diffusing away any magnetic divergence present, and does nothing otherwise (§ 2.4). Similarly, the term proportional to term can be used to help quell kinetic divergence in incompressible simulations, although for the simulations in this paper we use a spectral projection to accomplish this.
We clarify the meaning of equation (5), by assuming that
(7) 
The terms in equation (5) are seen to have the role of diffusivities in equation (7). In the simulations in § 3, we use equation (5) for constrained transport and equation (7) for unconstrained transport.
Time centering of the staggeredgrid and centeredgrid constrained transport equations encounter a similar circumstance as with the kinetic equation. In the staggeredgrid constrained transport configuration employed by Hawley and Stone (1995), the electric fields are spatially interpolated from the nearest 8 velocity and magnetic field vectors, and they are also timeinterpolated to the next half step with an Alfvén wave Method of Characteristics scheme (Hawley and Stone 1995). The implicit diffusivity inherent in these interpolations stabilizes the magnetic field evolution, even in the absence of explicit diffusivity. For RungeKutta timestepping on a centered grid, there are no spatial interpolations in constructing the electric field, and no accompanying intrinsic diffusivity. Some form of diffusivity is then generally required to maintain stability.
2.4 Divergence
In a constrained transport simulation, the magnetic divergence remains zero to machine precision. For unconstrained transport, we set the magnetic divergence to zero with a Poisson projection in Fourier space once every four timesteps. In tests, we found that the evolution is virtually identical whether the divergence is removed once every timestep or once every four timesteps (Maron 2004), and in both cases the fractional magnetic divergence remains below one percent. For both constrained and unconstrained transport simulations of incompressible flow that we describe below, the kinetic divergence is removed once every timestep in Fourier space, although for quasiincompressible flow, it suffices to do this only once every few timesteps (Maron 2004).
The divergence diffusivity term is helpful for reducing the effect of magnetic divergence in the timesteps between Poisson projections. It diffuses any magnetic divergence present, but otherwise does nothing. If each Fourier component has the form , then the term evolves as and . The divergence diffusivity is also helpful for a constrained transport simulation. Without it, CT preserves the magnetic divergence but does nothing to remove it. With the divergence diffusivity, any divergence present is diffused away.
2.5 Initialization
The constrained transport algorithm evolves the magnetic field in such a way as to conserve magnetic divergence. If the initial conditions have zero divergence, the divergence remains zero indefinitely. Even if any monopoles do grow slowly, they can be removed at negligible computational expense with a Fourier projection, say, once every thousand timesteps.Care must be taken with the initial conditions, however, because the divergence depends on the method for approximately evaluating derivatives. One may use, for example, a spectral or a finite difference divergence operator. For constrained transport to work, the derivatives used for making the initial conditions divergenceless must be the same as those used in the simulation. Even an analytic function with vanishing divergence may not have vanishing numerical divergence.
To initialize the magnetic field, we apply the following procedure. In three dimensions, let the wavenumbers for the magnetic field Fourier components be , the length of each side of the simulation volume be , and the number of grid zones on each side be . If the magnetic divergence is defined spectrally, then the constraint on the magnetic field is . For a finite difference derivative the constraint is slightly different. If we use the highorder finite difference from § 2.6 (eq. 11), we can define an adjusted dimensionless wavenumber by
(8) 
and then the finite difference condition for zero divergence is
(9) 
The initial conditions for a constrained transport simulation should satisfy this condition, which is easily implemented in Fourier space.
2.6 Highwavenumber finite differences
Highorder spatial derivatives can enhance the wavenumber resolution of a simulation. To quantify this, define a function on a periodic grid , with an integer. Then construct a finite difference derivative at using a radius stencil. The familiar result for the gradient on a radiusone stencil is , which is obtained from fitting a polynomial of order to at . For an order 4 polynomial on a radiustwo stencil,
(10) 
For a stencil of order ,
(11) 
where . The coefficients for a radiusthree stencil are .
Consider the finitedifference error at for a Fourier mode . (Cosine modes can be ignored because they don’t contribute to the derivative at .) We scale the wavenumber to grid units so that corresponds to the maximum (Nyquist) wavenumber expressible on the grid. The finite difference formula (eq. 11) gives
(12) 
whereas the correct value is . The spectral procedure of taking the derivative by transforming to Fourier space and back gives the correct value up to , but at a cost of floating point operations per grid point per transform, where is the number of grid points, whereas the finite difference derivative (equation 11) on a radius stencil costs floating point operations. Figure (1) exhibits the accuracy of finite difference derivatives of different orders as a function of wavenumber. The wavenumber resolution increases with polynomial order.
The polynomial fit can be tuned by adjusting the coefficients to enhance highwavenumber accuracy substantially at the expense of a negligible loss of accuracy at low wavenumbers (Maron 2004). Table 2 gives the coefficients for finite difference operators at various orders, along with the maximum wavenumber for which they are 0.5% accurate. In Figure 1 we also show the wavenumber accuracy for a radiusthree stencil using tuned coefficients, designed to have a relative precision of % for to . We will see in § 3.1.2 that a simulation based on this tuned radiusthree scheme performs better than a radiusfour simulation with conventional coefficients. We note that power spectra do not reveal a difference between turbulence simulations with tuned and untuned coefficients, because errors in the derivative manifest as an advective dispersion rather than as a diffusivity. One has instead to examine the fields in real space rather than in Fourier space.
Operation  

(P1) 
0  0.50000  0.12  
(P2) 
0  0.66667  0.08333  0.24  
(P3) 
0  0.75000  0.15000  0.01667  0.34  
(P4) 
0  0.80000  0.20000  0.03810  0.00357  0.40  
(P5) 
0  0.83333  0.23810  0.05952  0.00992  0.00079  0.44  
(P8) 
0  0.88889  0.31111  0.11313  0.03535  0.00870  0.00155  0.00018  0.00001  0.56 
(T3) 
0  0.81796  0.21324  0.03683  0.50  
(T8) 
0  0.95951  0.42312  0.22746  0.12450  0.06461  0.03004  0.01156  0.00273  0.76 
(T8) 
0  0.96685  0.43635  0.24410  0.14168  0.07976  0.04148  0.01884  0.00536  0.80 
(P3) 
2.72222  1.5  0.15  0.01111  
(P3) 
9.33333  6.5  2.  0.16667  
(P3) 
20.  15.  6.  1.  
Interpolation (P3) 
0  0.58594  0.09766  0.01172  0.35  
Interpolation (T3) 
0  0.60103  0.12312  0.02214  0.50  
Interpolation (T8) 
0  0.63099  0.19580  0.10149  0.05770  0.03248  0.01709  0.00793  0.00233  0.80 

.
2.7 Dealiasing
Nonlinear terms in the MHD equations such as suffer an aliasing error for highwavenumber structure. We first illustrate this with a 1D example before treating the 3D case (Canuto et al. 1987). Let and be discrete functions on a onedimensional grid with grid points: , and with accompanying Fourier modes ranging from The Fourier expansion of is
(13) 
and similarly for . Let and each be composed of a single Fourier mode with . The product is then composed of a single Fourier mode with , but in this discrete representation, this is equivalent to . This remapping of high wavenumber modes to low wavenumber modes is known as aliasing error. Note that is within the Nyquist limit of while the product mode is not. While and might be resolvable on the grid, their product need not be. This can be fixed by truncating before and after the product all modes outside to ensure that no modes in the product will exceed the Nyquist limit. A spectral code can make this truncation because the fields are in Fourier space, while a finite difference code cannot. As a practical matter it suffices to set the diffusivities high enough so that negligible structure exists outside the aliasing limit.
Alternatively, the aliasing problem can be remedied with a staggeredgrid correction (Canuto et al. 1987), and then the truncation isn’t necessary. Additionally, this allows us to simulate structure beyond We first exhibit this procedure in 1D and then extend it to 3D below. To implement the staggeredgrid correction, first construct the usual product . Then interpolate the centered grids and to the staggered grids and , multiply to produce and transform back to the centered grid . Finally, combine the centered and staggered results as . This result is free from aliasing error for all Fourier modes, and so no Fourierspace truncation is necessary. In Fourier space, the fields on the staggered grid can be computed exactly by applying a phase shift to each Fourier mode. A finite difference scheme can accomplish this with the highorder interpolation discussed in § 2.8.
In three dimensions the aliasing error can be eliminated by truncating the Fourier modes outside , where and . The gridshift correction is more complicated. Seven grid shifts are needed to completely correct the error, but it is more efficient to settle for a more limited correction involving just one shift (Canuto et al. 1987). Interpolate and to the staggered grid , multiply, return to the centered grid , and as in the onedimensional case, average the result with the product carried out on the original centered grid. This yields an aliasfree result if accompanied by a Fourier truncation of all modes outside , which is 94 % of the Nyquist limit. For a finite difference code, the diffusivity can be set to minimize the energy outside the truncation zone. When the threedimensional, staggeredgrid, aliasing correction is combined with highorder, volumecentered, constrained transport, magnetic divergence is still preserved to machine precision. This scheme uses the highorder finite differences and interpolations discussed in § 2.6 and § 2.8.
2.8 Interpolation
Highorder interpolation allows implementation of the staggered aliasing correction discussed in § 2.7. It can also be used for doubling the grid size. In this case, the new points in the doubled threedimensional grid can be generated with a set of interpolations along and diagonal to the three grid axes. Doubling the grid is also useful for interpolating to points that are not exactly halfway between grid points. For example, the grid can be doubled with the highorder interpolation discussed above, and then a simpler algorithm can be used to further interpolate to a point anywhere within the refined grid. We used this technique to trace the path of magnetic fieldlines in our studies of cosmic ray diffusion (Maron 2003), and also for doubling the resolution of turbulent dynamo simulations in Maron et al. (2004). Lastly, we have found that highorder tuned derivatives and interpolations are useful for postsimulation analysis of timeslice data without the need for computationally expensive fast Fourier transforms.
A highorder interpolation to a point halfway between two grid points can be accomplished with
(14) 
The error in the interpolation as a function of wavenumber can be calculated by applying equation (14) to a cosine mode centered on : . The interpolated value is
(15) 
whereas the correct value is . Table 2 gives coefficients for an interpolation based on a radiusthree stencil from a polynomial fit to , in the row labeled “Interpolation (P3).” Tuning the coefficients can improve the wavenumber resolution of the interpolation in a similar manner as was done for the derivative. Tuned coefficients for radius 3 and 8 stencils are given in Table (2) in the rows labeled T3 and T8, with the column labeled giving the maximum wavenumber for which the interpolation is accurate.
2.9 Operation count
A finite difference derivative on a radius stencil costs multiplies and adds. Since add and multiply units tend to come in pairs on most modern machines, this is effectively floating point operations per grid cell. A code with finite difference convolutions per timestep involves floating point operations.
The computation cost per timestep scales as the number of derivatives computed per grid element per time update. For constrained transport, one needs to calculate three derivatives for , three for , nine for , six for and six for . Diffusion terms such as and aren’t included in the tally because they need only be applied every few timesteps (Maron 2004). The same is true for the Fourierspace Poisson projection to make the fields divergenceless. Maron (2004) found that the results were effectively identical whether these terms were applied once every timestep or once every four timesteps.
In Table 3 we list the number of derivatives computed per timestep for various classes of algorithms. Pencil is a vector potential code (Brandenburg & Dobler 2002) that calculates directly from to take advantage of the cache memory. This involves 15 distinct 2level derivatives computed from a 2D stencil. For a radiusS finite difference, a 2D derivative involves adds, compared to adds for a conventional 1D derivative. We don’t consider multiples because they are always less numerous than adds, and because CPUs tend to feature add and multiply units in pairs. As a result, a 2D derivative corresponds to times the computational effort of a conventional 1D derivative. Thus, the calculation of with corresponds effectively to about conventional derivatives. Alternatively, a vector potential algorithm can instead calculates from and then from in separate stages, at a possible cost of extra memory latency. In table 3, we denote a vector potential algorithm that calculates directly from as “1stage,” and a vector potential algorithm that calculates from and then from as “2stage.” The constrained transport algorithm evolves and according to equations (4) and (5), and the unconstrained transport algorithm evolves and according to equations (4) and (7). The spectral algorithm computes 15 Fourier transforms instead of derivatives and thus doesn’t directly compare with the finite difference codes.
A vector potential code such as Pencil can be adapted to the algorithm described here by substituting equation (5) for the vector potential equation, and by changing the finite difference coefficients from the polynomialbased values to the tuned values. This allows one to compute boundary conditions using the magnetic field instead of the vector potential. However, there is an extra round of interprocessor communication associated with calculating from and then applying the curl.
Algorithm  Terms  Derivatives 

Constrained transport 
, , , ,  27 
Unconstrained transport 
, , , ,  30 
2stage vector potential 
, , , ,  27 
1stage vector potential (Pencil) 
, , , ,  120 
Spectral 
15 Fourier transforms  N/A 

A convenient unit for execution speed is Kilogrid Elements per CPU GigaFLOPS per Second (KEGS). A constrained transport code with a radiusthree stencil and two RungeKutta steps per timestep would run ideally at a speed of 1850 KEGS. The actual speed is slower of course because of overhead and because the finite difference convolutions don’t run at the peak floating point speed. In benchmarks, we have observed that the Pencil code, which has 3 RungeKutta stages, runs at a speed of around 100 KEGS for serial operation and
50 KEGS for massively parallel operation. For comparison, the highlyoptimized spectralMHD code Tulku (Maron & Goldreich 2001), has a speed of around 80 KEGS in serial and 40 KEGS in parallel (Maron 2004).
3 Simulations
3.1 The turbulent nonhelical dynamo
Our first test model is the turbulent, nonhelical, MHD dynamo (Batchelor 1950), which is the magnetic analogue of the Kolmogorov cascade for hydrodynamic turbulence. The Kolmogorov cascade is the longterm, steady state of hydrodynamic turbulence that is forced at the large scale and dissipated by viscosity at the small scale, and it has an energy spectrum of . The nonhelical MHD dynamo has the same setup but includes a spatially homogeneous magnetic field with unit magnitude and zero mean. Maron et al. (2004) found that the steady state of this system has a kinetic energy spectrum of , where the energy is dominantly in largescale eddys, and a magnetic spectrum with the energy predominantly in the smallestscale (resistivescale) magnetic structures. The kinetic and magnetic spectra are shown in figure 2.
Models
We ran a set of simulations of forced homogeneous MHD turbulence (Table 4) with and without constrained transport and at various spatial orders to test the effectiveness of highorder constrained transport. Each simulation has the same timestep, viscosity and resistivity. The grid used in all cases has zones covering a periodic unit cube. The forcing power, density, initial RMS velocity and magnetic field are all unity, and they remain so in the longterm steady state.
The initial state is the same for all simulations. It was taken from a simulation of forced magnetized turbulence (§ 3.1) in the longterm steady state. The magnetic divergence is zero in Fourier space: . The initial state was modified slightly for the constrained transport simulations to make it magnetically divergenceless according to the finitedifference derivative (§ 2.5). This amounts to a very small adjustment in the fields, much less than the difference between the fields after a crossing time of evolution.
Constrained transport needs some form of diffusivity, either Laplacian or hyperdiffusivity, to maintain stability. With diffusivity, it is stable indefinitely, whereas it is unstable without it. To establish this, we evolved the MHD equations with constrained transport for ten crossing times stably. The diffusivities used in these models are Laplacian viscosity , hyperviscosity and hyperresistivity , which we empirically find are the minimum values that maintain stability for on a grid. Expressed in dimensionless form: and these are and
Tag  Method  Stencil  Finite difference  norm  

radius  technique  
CT1  Constrained transport  1  Polynomial  0.12  0.561 
CT2  Constrained transport  2  Polynomial  0.24  0.193 
CT3  Constrained transport  3  Polynomial  0.34  0.083 
CT4  Constrained transport  4  Polynomial  0.40  0.042 
CT5  Constrained transport  5  Polynomial  0.44  0.023 
CT8  Constrained transport  8  Polynomial  0.56  0 
CT3t1  Constrained transport  3  Tuned  0.50  0.034 
CT3t2  Constrained transport  3  Tuned  0.52  0.037 
UT8  Unconstrained transport  8  Polynomial  0.56  0.267 

The forcing is the same as used by Maron et al. (2004). A random forcing field is added to the velocity every timestep. The spectrum of the forcing field is , truncated 2.5 lattice units from the origin in Fourier space, and the Fourier components have random phases. The forcing power, simulation volume and density are unity, which yields RMS velocity and magnetic fields of order unity (Maron 2004).
Results
The diffusivities are given in § 3. We plot the kinetic and magnetic spectra in figure 2(a). The spectra are very similar for constrained and unconstrained transport simulations, and also for different orders. However the spectra alone do not distinguish between simulations of different orders because an error in the derivative manifests itself as an advective dispersion rather than as a diffusivity (§ 2.6). One instead has to examine the fields in real space.
In Figure 2(b), we compare the magnetic fields at t=0.4 crossing times. For the comparison, we examine the difference between the fields integrated over space by computing the norm between simulations and
(16) 
Stone et. al. (1992b) argue that this kind of comparison is more meaningful than merely plotting the overlay of both fields. The constrained transport simulation with polynomialbased finite differences on a radiuseight stencil (CT8 in Table 4) serves as the basis of comparison. We compare the constrained transport simulations to an unconstrained transport simulation on a radiuseight tuned finitedifference stencil (UT8). We use UT8 as a standin for the spectral algorithm because of its high wavenumber resolution.
The spectral algorithm delivers the highestattainable resolution because spectral derivatives are exact for all wavenumbers. With this, a 3D spectral simulation without an aliasing gridshift correction can resolve structure up to k=2/3, and with a gridshift correction it can resolve up to k=.94 (Canuto 1987). The spectral algorithm can also set the magnetic divergence to zero in Fourier space at negligible cost. Unconstrained transport does not explicitly conserve magnetic divergence and so in model UT8 the divergence is cleaned with a Fourier projection every timestep. We also tried applying the correction every fourth timestep and with virtually identical results. The radiuseight stencil of UT8 yields derivatives that are accurate up to k=.56.
The norms given in Table (2) show how the simulations progressively approach the CT8 result as the stencil size increases. The match is poor for CT1 and better for CT3. We also note that the radiusthree simulation with tuned derivatives (CT3t) performs better than the radiusfour simulation with polynomialbased derivatives, establishing the effectiveness of tuned derivatives. This can also be qualitatively seen in Figure (2), where we see that the fields for CT8 and UT8 are closely aligned (Fig. 2), and that they also closely resemble those for CT3.
We attribute the remaining differences between CT8 and UT8 to the fact that the magnetic divergence is removed spectrally in UT8, while it is handled by constrained transport in CT8. Collectively, the highorder constrained and unconstrained transport simulations (CT3, CT8 & UT8) more closely resemble each other than they do the loworder constrained transport simulations (CT1 & CT2). We conclude that CT3 is already a good approximation to the spectral algorithm.
3.2 Comparison of the constrained transport and vector potential techniques
We adapted the vector potential code Pencil to run in CT mode and used it to compare the vector potential and CT techniques. We ran an Alfvén wave on a grid with zero viscosity and resistivity. (Figure 3). After ten crossing times, both the vector potential and CT techniques yield wave profiles that agree with each other to within 1 percent. The shape of the profiles are also wellmatched with the initial conditions, with a phase error of 10 percent.
We also used both techniques to run a turbulent dynamo simulation We started with an initially weak magnetic field in the form of a Beltrami wave and applied helical forcing until it grew to a steady state. The box size is , the density is unity, the forcing power is equal to the viscosity is equal to and the resistivity is equal to The RMS magnetic field strength is plotted in Figure (4). After 30 crossing times, the values for for the CT and vector potential techniques agree to 1 percent (Figure 4).
3.3 Oblique Alfvén wave test
We ran an Alfvén wave test where the propagation axis is oblique to the grid axes, with the initial conditions in Gardiner & Stone (2005):
(17) 
(18) 
The simulation volume is a unit cube, modeled on a grid of size . The velocity field is quasiincompressible, with the divergence removed spectrally every 4 timesteps. The kinetic and magnetic diffusivities are all set to zero for this linear problem. We ran two simulations: one with thirdorder polynomial finite differences and another with thirdorder tuned finite differences from Table (2). After the wave has traveled 16 times around the periodic box, the waveform remains almost indistinguishable from the initial conditions, with the tuned finite differences yielding a more precise result than the polynomial finite differences (figure 5).
4 Summary
We have developed a new version of the constrained transport algorithm that uses volumecentered fields and hyperresistivity on a highorder finite difference stencil, with tuned finite difference coefficients to enhance highwavenumber resolution. Highorder interpolation allows implementation of staggered dealiasing. Together, these measures yield a wavenumber resolution that approaches the ideal value achieved by the spectral algorithm.
Volume centered fields are desirable because then , and all reside at the same grid location, allowing to be constructed directly from the cross product of and without interpolation. For staggered fields, and reside at the zone faces and on the edges, and so constructing involves spatial interpolation, which reduces wavenumber resolution.
Highorder stencils and tuned finite difference coefficients both enhance the wavenumber resolution of finite differences. For a radiusthree stencil with tuned coefficients, derivatives can be computed to a relative precision of up to a Nyquistscaled wavenumber of . Without tuning, this would be for a radiusthree stencil. A radiusone stencil derivative such as is used in Zeus (Stone & Norman 1992a) is only accurate up to . The spectral derivative is precise up to , although in practice it is limited to because of aliasing. Aliasing limits a finitedifference code to unless the finitedifference grid shift aliasing correction is used (§ 2.7).
Hyperresistivity is desirable because it is more effective than Laplacian resistivity in diffusing highwavenumber modes while at the same time preserving lowwavenumber modes. The fact that hyperresistivity can be written as a curl allows its inclusion into CT. If Laplacian diffusivity were used instead, too much highwavenumber structure would be diffused for the highorder or tuned derivatives to matter.
The resolution of the algorithm described here approaches that of a spectral code, but because it uses finite differences, it runs faster than a spectral code and isn’t restricted to periodic boundary conditions. Also, since the finite differences are local, it is easily scalable to thousands of processors. The spectral algorithm is more difficult to scale to large numbers of processors because it involves alltoall communications between processors. A finite difference code only passes information between processors whose subgrids are adjacent in physical space. Lastly, because the code works with the magnetic field rather than the vector potential, boundary conditions are often easier to implement.
Appendix A Appendix
Constrained transport expresses the magnetic induction equation as a pure curl plus a divergence diffusivity: where is defined in equation (5). The term serves to diffuse away magnetic divergence, and the finite differences are arranged so that Thus, the curl term does not contribute to the evolution of the magnetic divergence, and if the initial conditions are divergencefree, the magnetic divergence remains zero throughout the evolution.
To see how constrained transport works, denote the vector field by = where are integers specifying the locations of grid cell centers. There are two basic grid types: “centered” and “staggered” (figure 6). For a centered grid, scalar and vector quantities are located at cell centers. For a staggered grid, scalar quantities are located at cell centers and vector quantities at cell faces. For instance, we would index the components of as
The finite divergence divergence of the curl of is which consists of terms such as One can straightforwardly see that this is zero for finite differences of the form eq. 11 for both centered and staggered grids. Thus, constrained transport can be coordinated with highorder and tuned finite differences, as well as with hyperresistivity.
For a staggered grid, the “” vectors are located at cell edges, whereas the and vectors from which they are constructed are found at cell faces. A staggered grid CT scheme therefore involves spatial interpolation, one example being the method of characteristics scheme for timeinterpolating Alfvén waves. We use volumecentered fields and RungeKutta timestepping because, among other reasons, no interpolation is required.
Footnotes
 affiliation: badger@amnh.org
 affiliation: mordecai@amnh.org
 affiliation: Also Department of Astronomy, University of Virginia, Charlottesville, VA; joishi@amnh.org
References
 Batchelor, G. K. 1950, Proc. Roy. Soc. Lon., A201, 405
 Brackbill, J. & Barnes, D. 1980, J. Comput. Phys., 35, 426
 Brandenburg, A., & Dobler, W. 2002, Comput. Phys. Comm., 147, 471
 Canuto, C., Hussaini, M. Y., Quarteroni, A., Zang, T. A., Jr. 1987, Spectral Methods in Fluid Dynamics (Heidelberg: SpringerVerlag)
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174
 Dedner, A., Kemm, F., Kroner, D., Munz, C.D., Schnitzer, T. & Wesenberg, M. 2002, J. Comput. Phys., 175, 645
 Dobler, W., Stix, M. & Brandenburg, A. 2006, ApJ 638, 336
 Evans, C.R. & Hawley, J.F. 1988, ApJ, 332, 659
 Gardiner, T. & Stone, J., J. Comput. Phys., 205, 509
 Harten, A. 1983, J. Comput. Phys., 49, 357
 Haugen N. E. L., & Brandenburg, A. 2004, PRE 70, 026405
 Hawley, J. F., & Stone, J. M. 1995, Comp. Phys. Comm., 89, 127
 Londrillo, P. & Del Zanna, L. 2000, ApJ, 530, 508
 Maron, J., 2004, J. Comput. Phys., submitted (astroph/0402606)
 Maron, J., Chandran, B. & Blackman, E. 2004, PRL, 92, 045001
 Maron, J. & Goldreich, P. 2001, ApJ, 554, 1175
 Powell, K.G., Roe, P.L., Linde, T.J., Gombosi, T.J. & De Zeeuw, D.L. 1999, J. Comput. Phys., 154, 284
 Stone, J. M. & Norman, M. L. 1992a, ApJS, 80, 791
 Stone, J., Hawley, J.F., Evans, C.R. & Norman, M. 1992b, ApJ, 388, 415
 Tóth, G. 2000, J. Comput. Phys., 161, 605
 Toro, E. F. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (Heidelberg: SpringerVerlag)
 Zhong, X. 1998, J. Comp. Phys. 144, 662