A Appendix

# A constrained-transport magnetohydrodynamics algorithm with near-spectral 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 volume-centered fields and hyperresistivity on a high-order finite difference stencil. Additionally, the finite-difference coefficients can be tuned to enhance high-wavenumber 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 low-Mach-number turbulence, the results agree well with a high-order, non-constrained-transport 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 three-dimensional, multiscale and turbulent, so there is an ongoing interest in developing high-resolution 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 zone-centered grids, and with resolution-enhanced 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

 ∂tB=\boldmath∇×(V×B). (1)

Analytically, this equation conserves magnetic divergence: . However, this may or may not be the case for a finite-difference 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 higher-order advection algorithm to ensure accurate second-derivatives, 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 high-order CT is possible on staggered grids with a radius-two stencil. In this paper, we show that volume-centered 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.

In § 2, we discuss the specifics of the algorithm, and in § 3, we describe test simulations that demonstrate the capabilities of CT with high-order spatial derivatives.

## 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 high-order and tuned finite differences plus hyperdiffusivity, and stability is achieved by Runge-Kutta timestepping and hyperdiffusivity.

### 2.1 Timestepping

A high-order 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 second-order Runge-Kutta 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 time-extrapolating 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 well-known method-of-characteristics algorithm interpolates along the Alfvén characteristic while neglecting the fast and slow mode characteristics. For our simulations, we use Runge-Kutta for time-extrapolation 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 second-order Runge-Kutta while the Pencil code (Brandenburg & Dobler 2002) uses third-order, although either order has proven successful.

### 2.2 Diffusivity

A common class of algorithms is based on momentum fluxes that are time-extrapolated 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).

Runge-Kutta 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 hyper-Laplacian 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:

 ∂tV=ν2\boldmath∇2V−ν4\boldmath∇4V+ν6\boldmath∇6V−ν[4]\boldmath∇[4]V. (2)

Let the Fourier components be where is the wavenumber. They evolve according to

 ∂t^V(k)=−ν2k2^V(k)−ν4k4^V(k)−ν6k6^V(k).−ν[4](k4x+k4y+k4z)^V(k). (3)

The term is the Laplacian viscosity and the others are higher-order hyperdiffusivities. Specifically, and They differ in that the higher the order, the more selective they are in diffusing high-k structure while preserving low-k structure. Equivalently stated, high-order hyperdiffusivities erase small-scale 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 small-scale 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 maximum-possible 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 high-k 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 high-order finite differences to enhance the resolution of a simulation. High-order 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 highest-k structure, increases the maximum stable timestep by about 50%. This technique works because the highest-wavenumber 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

 ∂tV=−V⋅\boldmath∇V+B⋅\boldmath∇B+ν2% \boldmath∇2V−ν4% \boldmath∇4V−ν[4]% \boldmath∇[4]V+ν6% \boldmath∇6V+ν[6]% \boldmath∇[6]V+νD% \boldmath∇\boldmath∇⋅V (4)
 ∂tB=\boldmath∇×[V×B−η2J+η4\boldmath∇2J−η6% \boldmath∇4J−η2,[4]% \boldmath∇[4]J]+ηD\boldmath∇\boldmath∇⋅B (5)
 J=\boldmath∇×B (6)

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

 ∂tB=B⋅\boldmath∇V−V⋅\boldmath∇B−B\boldmath∇⋅V
 +η2\boldmath∇2B−η4\boldmath∇4B+η6\boldmath∇6B+η2,[4]\boldmath∇2\boldmath∇[4]B+ηD\boldmath∇\boldmath∇⋅B (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 staggered-grid and centered-grid constrained transport equations encounter a similar circumstance as with the kinetic equation. In the staggered-grid 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 time-interpolated 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 Runge-Kutta 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 quasi-incompressible 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 high-order finite difference from § 2.6 (eq. 11), we can define an adjusted dimensionless wavenumber by

 k∗i=S∑j=−Smjsin(kiLij/Ni), (8)

and then the finite difference condition for zero divergence is

 ^B(k)⋅k∗=0. (9)

The initial conditions for a constrained transport simulation should satisfy this condition, which is easily implemented in Fourier space.

### 2.6 High-wavenumber finite differences

High-order 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 radius-one stencil is , which is obtained from fitting a polynomial of order to at . For an order 4 polynomial on a radius-two stencil,

 f′(0)∼112f−2−23f−1+23f1−112f2. (10)

For a stencil of order ,

 f′(0)∼S∑j=−Smjfj (11)

where . The coefficients for a radius-three stencil are .

Consider the finite-difference 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

 f′k(0)∼2S∑j=1mjsin(πjk) (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 high-wavenumber 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 radius-three 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 radius-three scheme performs better than a radius-four 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.

### 2.7 Dealiasing

Nonlinear terms in the MHD equations such as suffer an aliasing error for high-wavenumber 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 one-dimensional grid with grid points: , and with accompanying Fourier modes ranging from The Fourier expansion of is

 Aj=N/2−1∑s=−N/2^Asexp(2πijs/N) (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 staggered-grid 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 staggered-grid 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 Fourier-space 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 high-order interpolation discussed in § 2.8.

In three dimensions the aliasing error can be eliminated by truncating the Fourier modes outside , where and . The grid-shift 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 one-dimensional case, average the result with the product carried out on the original centered grid. This yields an alias-free 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 three-dimensional, staggered-grid, aliasing correction is combined with high-order, volume-centered, constrained transport, magnetic divergence is still preserved to machine precision. This scheme uses the high-order finite differences and interpolations discussed in § 2.6 and § 2.8.

### 2.8 Interpolation

High-order 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 three-dimensional 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 high-order 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 high-order tuned derivatives and interpolations are useful for post-simulation analysis of timeslice data without the need for computationally expensive fast Fourier transforms.

A high-order interpolation to a point halfway between two grid points can be accomplished with

 fj+1/2=S∑i=1mi(fj+i+fj+1−i). (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

 f∗j+1/2=S∑i=1mjcos(πik) (15)

whereas the correct value is . Table 2 gives coefficients for an interpolation based on a radius-three 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 Fourier-space 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 2-level derivatives computed from a 2D stencil. For a radius-S 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 “1-stage,” and a vector potential algorithm that calculates from and then from as “2-stage.” 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 polynomial-based 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 inter-processor communication associated with calculating from and then applying the curl.

A convenient unit for execution speed is Kilo-grid Elements per CPU GigaFLOPS per Second (KEGS). A constrained transport code with a radius-three stencil and two Runge-Kutta 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 Runge-Kutta stages, runs at a speed of around 100 KEGS for serial operation and

50 KEGS for massively parallel operation. For comparison, the highly-optimized spectral-MHD 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 long-term, 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 large-scale eddys, and a magnetic spectrum with the energy predominantly in the smallest-scale (resistive-scale) 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 high-order 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 long-term 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 long-term 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 finite-difference 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

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

 A2ij=∫[By(i)−By(j)]2d(Vol)∫By(i)2dVol+∫By(j)2d(Vol). (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 polynomial-based finite differences on a radius-eight stencil (CT8 in Table 4) serves as the basis of comparison. We compare the constrained transport simulations to an unconstrained transport simulation on a radius-eight tuned finite-difference stencil (UT8). We use UT8 as a stand-in for the spectral algorithm because of its high wavenumber resolution.

The spectral algorithm delivers the highest-attainable resolution because spectral derivatives are exact for all wavenumbers. With this, a 3D spectral simulation without an aliasing grid-shift correction can resolve structure up to k=2/3, and with a grid-shift 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 radius-eight 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 radius-three simulation with tuned derivatives (CT3t) performs better than the radius-four simulation with polynomial-based 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 high-order constrained and unconstrained transport simulations (CT3, CT8 & UT8) more closely resemble each other than they do the low-order 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 well-matched 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):

 V={0,0.1sin(2πx+4πy),0.1cos(2πx+4πy)} (17)
 B={1,0.1sin(2πx+4πy),0.1cos(2πx+4πy)} (18)

The simulation volume is a unit cube, modeled on a grid of size . The velocity field is quasi-incompressible, 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 third-order polynomial finite differences and another with third-order 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 volume-centered fields and hyperresistivity on a high-order finite difference stencil, with tuned finite difference coefficients to enhance high-wavenumber resolution. High-order 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.

High-order stencils and tuned finite difference coefficients both enhance the wavenumber resolution of finite differences. For a radius-three stencil with tuned coefficients, derivatives can be computed to a relative precision of up to a Nyquist-scaled wavenumber of . Without tuning, this would be for a radius-three stencil. A radius-one 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 finite-difference code to unless the finite-difference grid shift aliasing correction is used (§ 2.7).

Hyperresistivity is desirable because it is more effective than Laplacian resistivity in diffusing high-wavenumber modes while at the same time preserving low-wavenumber modes. The fact that hyperresistivity can be written as a curl allows its inclusion into CT. If Laplacian diffusivity were used instead, too much high-wavenumber structure would be diffused for the high-order 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 all-to-all 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.

We received support for this work from NSF Career grant AST99-85392, NSF grants AST03-07854 and AST06-12724, and NASA grant NAG5-10103. We acknowledge stimulating discussions with E. Blackman, A. Brandenburg, B. Chandran, and J. Stone, and we also acknowledge the referee, Wolfgang Dobler, for thorough comments that improved the paper.”

## 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 divergence-free, 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 high-order 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 time-interpolating Alfvén waves. We use volume-centered fields and Runge-Kutta timestepping because, among other reasons, no interpolation is required.

### Footnotes

2. affiliation: mordecai@amnh.org
3. affiliation: Also Department of Astronomy, University of Virginia, Charlottesville, VA; joishi@amnh.org

### References

1. Batchelor, G. K. 1950, Proc. Roy. Soc. Lon., A201, 405
2. Brackbill, J. & Barnes, D. 1980, J. Comput. Phys., 35, 426
3. Brandenburg, A., & Dobler, W. 2002, Comput. Phys. Comm., 147, 471
4. Canuto, C., Hussaini, M. Y., Quarteroni, A., Zang, T. A., Jr. 1987, Spectral Methods in Fluid Dynamics (Heidelberg: Springer-Verlag)
5. Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174
6. Dedner, A., Kemm, F., Kroner, D., Munz, C.-D., Schnitzer, T. & Wesenberg, M. 2002, J. Comput. Phys., 175, 645
7. Dobler, W., Stix, M. & Brandenburg, A. 2006, ApJ 638, 336
8. Evans, C.R. & Hawley, J.F. 1988, ApJ, 332, 659
9. Gardiner, T. & Stone, J., J. Comput. Phys., 205, 509
10. Harten, A. 1983, J. Comput. Phys., 49, 357
11. Haugen N. E. L., & Brandenburg, A. 2004, PRE 70, 026405
12. Hawley, J. F., & Stone, J. M. 1995, Comp. Phys. Comm., 89, 127
13. Londrillo, P. & Del Zanna, L. 2000, ApJ, 530, 508
14. Maron, J., 2004, J. Comput. Phys., submitted (astro-ph/0402606)
15. Maron, J., Chandran, B. & Blackman, E. 2004, PRL, 92, 045001
16. Maron, J. & Goldreich, P. 2001, ApJ, 554, 1175
17. Powell, K.G., Roe, P.L., Linde, T.J., Gombosi, T.J. & De Zeeuw, D.L. 1999, J. Comput. Phys., 154, 284
18. Stone, J. M. & Norman, M. L. 1992a, ApJS, 80, 791
19. Stone, J., Hawley, J.F., Evans, C.R. & Norman, M. 1992b, ApJ, 388, 415
20. Tóth, G. 2000, J. Comput. Phys., 161, 605
21. Toro, E. F. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (Heidelberg: Springer-Verlag)
22. Zhong, X. 1998, J. Comp. Phys. 144, 662
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters