Computer Physics Communications

Computer Physics Communications ELSEV1ER

Computer PhysicsCommunications89 (1995) 127-148

MOCCT: A numerical technique for astrophysical MHD John F. Hawley a, James M. Stoneb a Virginia Institute for Theoretical Astronomy, Department of Astronomy, University of Virginia, PO Box 3818, Charlottesville, VA 22903, USA b Department of Astronomy, University of Muryland, College Park, MD 20742, USA

Abstract

Magnetic fields are important components of many astrophysical systems. These systems are studied with numerical simulations of the equations of magnetohydrodynamics (MHD). We briefly review the MHD assumptions, some of the numerical approaches currently in use, and some of the problems to which MHD codes are being applied. We describe in detail the Method of Characteristics-Constrained Transport scheme. This is a time-explicit, operator split, nonconservative technique that uses information propagated along Alfvrn characteristics to compute electromagnetic forces at intermediate time levels for the induction equation and for the Lorentz force. The advantages and difficulties with this scheme are investigated with a series of test problems.

1. Introduction

Magnetic fields are ubiquitous in astrophysical systems. They are present in planets, in the sun and other stars, and they are observed in interstellar and intergalactic space. The magnetic field strengths measured in many of these sites often have magnetic energy densities as large or larger than other energy densities in the system (such as the thermal or kinetic energy densities). Hence the study of magnetized plasmas is of great importance to the field of astrophysical fluid dynamics. Despite their importance, many of the most fundamental questions regarding the origin, effect, and evolution of magnetic fields in astrophysical systems remain unanswered. For example, the basic question of how strong cosmical magnetic fields originate is still uncertain. It is thought that strong astrophysical magnetic fields are generated via amplification of an initially small seed field through a dynamo process. This process is associated with bulk fluid motions that occur in such diverse places as the core of the Earth, the convection zone of the Sun and other stars, and within the interstellar medium of the galaxy itself. Despite the wide range of physical conditions that result in dynamo action, we have only the most basic ideas as to how the dynamo mechanism actually operates. In addition to the dynamo studies there are many other areas of active research within astrophysics that involve detailed modeling of magnetic processes. One dramatic example of a phenomenon that is thought to require magnetic fields is the astrophysical jet. Jets are seen in both stellar and galactic systems, on scales ranging from parsecs to megaparsecs. They are generally interpreted as beams of high energy plasma shot out from a central power source, or "engine". On the stellar scale this engine is often a protostar; for galactic 0010-4655/95/$09.50 (~) 1995 Elsevier Science B.V. All rights reserved SSD10010-4655(95)00190-1

128

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

systems it is believed to be a massive black hole. Since the galactic jets are observed in synchrotron radio emission, there is direct evidence for the presence of magnetic fields. There is also good evidence that jets are collimated by strong magnetic fields, and that magnetic processes are involved in launching them as well. Naturally the list of astrophysical problems in which magnetic fields are thought to play an important role goes on, ranging from the dynamics of the sun's atmosphere and the interaction of the solar wind with planetary magnetospheres, to the dynamics and evolution of interstellar clouds, star-forming regions, neutron stars, and active galactic nuclei. Each of these is a complex system, exhibiting both intricate spatial structure and time variability. To study these phenomena we must solve time-dependent partial differential dynamical equations in several spatial dimensions. Although many potential astrophysical investigations will require physics besides that governing the magnetic fields (such as radiative transport, fully and partially ionized plasma physics, nuclear reactions, molecular physics, etc.) in this paper we will restrict our attention to the equations of magnetohydrodynamics (MHD). As we discuss more fully in Section 2, MHD considers the macroscopic interaction between electromagnetic fields and a conducting plasma subject to a variety of assumptions that are nonetheless valid in many astrophysical systems. The perfect MHD limit results from setting the resistivity of the plasma to zero. Astrophysical plasmas are often well approximated by this limit, unlike many terrestrial applications which use the imperfect MHD equations that retain a finite resistivity. Magnetic field strengths in astrophysical settings range from extremely weak, where the magnetic field is carried along by a relatively unperturbed fluid flow (passive MHD), to the field-dominated limit where hydrodynamic forces can be neglected. The strength of the field is often parameterized by the plasma fl value, defined as the ratio of the gas to the magnetic pressure, fl = /°gas/(B2/87r). In a typical astrophysical context fl ,-~ 1 and neither the magnetic field nor the compressible gas dynamics can be neglected. This stands in contrast to magnetically-confined fusion research which focuses on the fl << l limit. Although MHD represents a dramatic simplification to the general problem of the dynamics of a plasma, it is still a great deal more complicated than one-fluid compressible hydrodynamics so often used in astrophysical studies. However, even one-fluid hydrodynamical studies are challenging for astrophysical researchers. Timedependent multidimensional hydrodynamics is sufficiently complicated that it is tractable only through the use of numerical codes run on supercomputers. Moreover, typical astrophysical problems involve a disparate range of time and length scales making it difficult to carry out large-scale global simulations. This remains true for MHD which is further complicated by new wave families, discontinuities, and other physical effects, all of which increase the difficulty of achieving an accurate solution. As an example, Maxwell's equations consist of two evolution equations and two constraint equations. Although the analytic evolution equations automatically satisfy the constraints, this need not be true for numerical implementations of the equations, i.e., numerical errors can result in the generation of unphysical magnetic monopoles. Another difficulty is that in nonrelativistic problems the propagation speed of MHD waves can become much larger than the fluid velocities or the sound speed (for a low fl plasma); this can bring the Courant limit of a time-explicit numerical computation down to the point where further evolution is impossible. The behavior of MHD flows changes dramatically from the passive field to the strong field limit, and any attempt to develop a numerical scheme that can go easily from one limit to the other will face many daunting obstacles. Many years of effort have gone into the search for good numerical algorithms for compressible hydrodynamics to bring the art to its present relatively mature form. Considerably less time has been spent on developing the more complicated MHD algorithms, so that while substantial progress has been made, it remains an area of active algorithmic development and testing. In this paper we review some numerical techniques for solving the MHD equations, and describe in detail the particular scheme that we are currently using. It is important to begin (in Section 2) with a description of the applications for which this code is written since this determines the range in fl over which the techniques must work, as well as the performance specifications required to answer the physically-motivated questions. In order to highlight more clearly the reasoning behind our particular algorithm, we present in Section 2 a general discussion of issues in testing and evaluating MHD methods. We also briefly review some alternate

J.E Hawley, J.M. Swne / Computer Physics Communications 89 (1995) 127-148

129

numerical approaches to the MHD equations that have been used in astrophysical research. Having outlined some astrophysical applications for which one would like to have an MHD code, we then proceed in Section 3 to describe in detail the scheme that we are currently using. However, as we (and others) continue to test and improve MHD algorithms, what we present here will not be the last word on numerical MHD. In Section 4 wc sum up with a discussion of future applications and directions in astrophysical MHD.

2. Background 2.1. Astrophysical M H D There are many different algorithms available to solve the equations of hydrodynamics; we expect a similar diversity for MHD. It is unlikely that any one numerical technique will prove adequate for the diverse range of applications required for astrophysical studies. One's goal should not be a search for the best scheme, but for the most appropriate scheme for a specific problem. In this section we briefly describe some of the ideas behind approaches currently in use or being studied for use in astrophysical problems. This review will be somewhat subjective since our impressions of the fortes and foibles of these techniques were obtained without having extensive personal experience with most of them. The literature lacks detailed information about astrophysical MHD codes, in part because this area of research is relatively new. This paper is an effort to provide such information about our MHD techniques. Hence our intention in this review has less to do with rating or comparing techniques than it does with providing the background that motivated our algorithm's development. Standard magnetobydrodynamics is based upon several assumptions. The first is that the system is nonrelativistic and all velocities are much less than the speed of light c. Among other things this means that the displacement currents in Maxwell's equations can be neglected. While there are important astrophysical systems where this condition will be violated, such as accretion disks near black holes (where orbital speeds are near c), it is acceptable in many, more prosaic contexts. The next approximation is that the local net electric charge is very small, i.e., there is local charge neutrality. We further assume that the plasma is completely ionized and can be described with a one-fluid treatment. This approximation is certainly adequate for many astrophysical situations, although, again, there are important exceptions, such as cold molecular clouds where a two-fluid ion-neutral plasma is required, And finally, an MHD treatment of the plasma is only appropriate when a ttuid description is valid, that is, when the length scales of interest are sufficiently large compared to the plasma length scales, e.g., particle gyroradii. With these standard assumptions one arrives at the equations of MHD. In conservative form they are (in Gaussian units) c),p + X7 • ( p v ) = 0,

(2.1)

i~t(S) + • .

(2.2)

3te + V .

( S v + I ( P + B2/87r) - BB/4rr) = 0, Pv( ½pv2 + ( y _7 l

atB = ~7 × E.

+

x E =0,

(2.3) (2.4)

The symbols have their usual meaning, with p the fluid density, P the fluid pressure, v the velocity, e the total energy, B the magnetic field, S =-- pv the momentum, o- the conductivity, and I is the unit vector. We have assumed an equation of state P = pc(y-

1),

(2.5)

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

130

where • is the specific internal energy. We also assume that the electric field E is specified from Ohm's law

c E = - v x B + J/o-.

(2.6)

The current J is given by 4 7 r j = V x B.

(2.7)

c

In addition the magnetic field is constrained by the no-monopole condition V . B = 0.

(2.8)

It is important to note that V' • B is a constraint, not an evolution equation. Analytically this constraint is preserved for all time, but this cannot be generally guaranteed in numerical evolution. The equations that we evolve with our finite difference code are modified from those above, first by assuming the ideal MHD limit and dropping the terms that depend upon the conductivity. Second, we do not use the conservative form for the momentum conservation equation (2.2) and the energy conservation equation (2.3). Instead, we construct a nonconservative momentum equation that represents the magnetic force as a source term,

at(S)+(V.Sv)=[~----~(VxB×B)-XTP] 1 ( B .V)B-V = ~--~

( P + 8~-) .

(2.9,

The advantage of (2.9) is that terms in the magnetic force that are proportional to V • B have been removed analytically. The energy equation is written in terms of the internal energy pc,

O,(pe) + V . ( p • v ) = - ( P + Q ) ( V . v),

(2.10)

where Q is an artificial bulk viscosity included to account for nonadiabatic shock heating. One can choose to evolve either the conservative or nonconservative form of the equations and there are advantages and disadvantages to both approaches. The most obvious advantage of the conservative form is that it is possible to construct straightforward finite difference operators that will conserve total mass, energy and momentum. One potential difficulty with the total energy equation is that the internal energy (and hence the pressure) must be obtained from the total energy by subtracting off the contributions from kinetic and magnetic energies. Thus, conservation alone doesn't guarantee entropy conservation (for an adiabatic flow), or that entropy will only increase in a given fluid element (for dissipative flow). In addition, if a flow is highly supersonic, the truncation error in total energy may exceed the internal energy, leading to obvious errors such as negative pressures. Brackbill and Barnes (1980) point out a potential problem with the conservative form of the momentum equation. Since an arbitrary numerical implementation of the induction equation (2.4) does not necessarily guarantee that the constraint ~7. B = 0 will be satisfied, any resulting nonzero divergence of B leads to unphysical accelerations due to monopoles in the direction of the magnetic field if the conservative form of the momentum Eq. (2.2) is used. In this form, the force contains terms proportional to V • B which, while analytically zero, may be numerically nonnegligible. Such forces can ruin a simulation even if they appear only at the level of the truncation error. Brackbill and Barnes (1980) demonstrate that this problem is avoided if an equation of the form of (2.9) is used and those force terms are removed analytically prior to finite-differencing the momentum equation.

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

131

It is nevertheless useful to consider how one solves systems of equations in conservative form with a finite-difference technique. Such systems can be formulated as

cgt

+ - -

cgxi

-0,

(2.11)

where U represents a column matrix of the fundamental variables, F ( U ) a matrix of fluxes constructed out of those variables, and i goes over all spatial dimensions. In a finite difference formulation U is discretized in both space and time, located at time levels labeled by n, and spatial locations labeled by integers i, j, and k. The fluxes are at half-integral space and time locations so that (2.11) can be written as a centered second-order difference approximations. The difficulty, of course, is in obtaining accurate flux values at the appropriate timeand space-centered locations from the available data, i.e., U n, the fundamental variables at timestep n. This requires an accurate technique to interpolate in space and extrapolate in time. There are many numerical approaches that have been applied to equations of the form of (2.11). For MHD one of the most popular is the two-step Lax-Wendroff technique, or related variations (see Roache, 1972). For example, Shibata and Uchida (1985) have studied the transients associated with the spin-up of a vertical magnetic field threaded through an accretion disk using a 2.5 dimensional (2.5D) Lax-Wendroff MHD code. (The term "2.5 dimensional" refers to a code that fully evolves two spatial dimensions and retains the vector components of the third, symmetric direction. An example is the retention of angular momentum and toroidal field in an axisymmetric computation in spherical or cylindrical coordinates.) One of the most extensive applications of this approach to three-dimensional (3D) MHD has been by groups investigating the global structure of the Earth's magnetosphere. For examples see Kageyama, Watanabe, and Sato (1992), Watanabe and Sato (1990), Fedder and Lyon (1987), and Ogino et al. (1986). These authors discuss some of the difficulties associated with their numerical approaches. The two-step Lax-Wendroff technique is used by Ogino et al (1986), while Watanabe and Sato (1990) and Kageyama, Watanabe and Sato (1992) use a Runge-Kutta time evolution. However, as discussed by Kageyama et al. (1992) both the Lax-Wendroff and their Runge-Kutta technique exhibit considerable phase error, manifested primarily as high-frequency ringing near discontinuities. We discuss the issue of phase error below in Section 3.3. This ringing can be smoothed through extensive use of numerical viscosities (Ogino et al., 1986), but, of course, at the cost of greater diffusion error. For example, Watanabe and Sato (1990) carry out Alfv6n wave test problems that show one of the disadvantages of Lax-Wendroff schemes that use an explicit numerical viscosity: large numerical damping on short wavelengths. These are well-known drawbacks of these schemes. Another approach that makes use of the conservative form of the equations is the Godunov technique (Godunov, 1959). In the Godunov technique the time-advanced fluxes F ( U ) are obtained through the solution of nonlinear characteristic equations. For compressible hydrodynamics, this approach has met with considerable success. The best example is the Piecewise Parabolic Method (PPM) code of Colella and Woodward (1984) which uses third-order piecewise representations of the fundamental variables and a nonlinear Riemann-solver to obtain fluxes F that are second-order in time. The technique is extended to more than one dimension through directional splitting (Strang, 1968; Gottlieb, 1972). In directional splitting a multidimensional finite-difference operator Lx~.. is constructed out of symmetrized one dimensional operators, i.e., Lx~,/_7.= 7(L~L~,I + L:.Lx)U. So long as the one-dimensional operators Lx and Ly are each second-order and time-symmetric, the procedure will be second-order. Intuitively, one can see that the successive application of ID operators will account for cross terms (such as c)x0:,) whereas constructing one operator out of separate, entirely ID components will not. Given the success of the hydrodynamic Godunov approach, it is natural to try to extend the technique to MHD. Brio and Wu (1988), Zachary and Colella (1992), Zachary, Malagoli, and Colella (1994), and Dai and Woodward (1994a,b) have made such an attempt. The algorithm is still under development, but initial results certainly appear promising. The difficulty is that the MHD Riemann problem is considerably more complex than the hydrodynamic version. There are additional wave families, each with its own left- and right-going characteristic. These new families are the fast and slow magnetosonic waves and the Alfv6n waves, all of which

132

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

can lead to new types of discontinuities and rarefactions. Many of these wave families become degenerate in the high /3 (low field) limit, further adding to the complexity of developing a general numerical Riemann solver. When these techniques are fully developed one can anticipate that they will have all the advantages seen in hydrodynamic versions, namely highly accurate shock and wave propagation. The disadvantages of this approach include their complexity, and the difficulties encountered if one wishes to extend the method to include additional physics in a self-consistent fashion. One of the more ambitious attempts to apply numerical techniques to an astrophysical MHD problem has recently been carried out by Nordlund et al. (1992) in their study of magnetized convection in the sun. Their code was 3D, and included both magnetic fields and radiative transport. Their approach was straightforward finite-differencing of the nonconservative equations, using a vector potential A, defined by B = V x A as the fundamental variable rather than the magnetic field itself. As discussed below, this is one way of ensuring that the solenoidal constraint is preserved in a numerical simulation. Numerical derivatives were obtained using splines over several grid points. Explicit dissipation and resistivity were added to eliminate dispersion error and to maintain smooth solutions. This work is an example of how a complicated problem, involving many physical effects, can be studied with a direct finite difference scheme. It could not have been done with a Godunov-type scheme since the appropriate Riemann-solver simply does not now exist. By comparing the hydrodynamic portion of this code with the results from a PPM simulation on a turbulence test problem, however, Stein (private communication) has found that while the PPM scheme was less diffusive, particularly at smaller scales, the large-scale features of the two simulations were quite similar. There are several alternative numerical schemes that have also been investigated. One such is the Method of Characteristics (MOC). An MOC test code was constructed by Lou, Rosner and Ulmschneider (1987). While instructive, such an approach has had limited application in hydrodynamics, mainly because it is incapable of obtaining solutions that contain discontinuities, an almost certain state of affairs in astrophysical simulations. A more radical departure is offered by abandoning finite-difference methods altogether. Philips and Monaghan (1985) have attempted to incorporate MHD into the smooth-particle hydrodynamics (SPH) technique. SPH has seen considerable use for hydrodynamic astrophysical problems, although its extension to the MHD system may strain the limits of the ability of discrete particles to model a continuum. One advantage of the particle approaches is that they are Galilean invariant, unlike more traditional finite difference techniques where the truncation error is a function of the velocity.

2.2. Accretion disks Our efforts to develop a multidimensional MHD code are motivated by our interest in the structure and evolution of accretion disk systems. An accretion disk consists of gravitationally bound gas in orbit around a compact star or black hole. The gas orbits that central star in nearly circular, Keplerian, coplanar orbits forming a disk that extends radially outward. These disks are generally thin in the sense that their vertical thickness H(R) is much less than the radial distance R from the star. This in turn requires that the sound speed of the gas, c,., be much less than the orbital velocity, RS2, where .(2 is the angular velocity. Note that since the orbital speeds are supersonic, shocks are a distinct possibility. Now the gas in the disk would simply orbit indefinitely were it not for torques that remove angular momentum from a given fluid element. This allows that fluid element to drop to a lower orbit. As the gas spirals in, it releases gravitational binding energy that is thermalized and radiated through the disk surface. The accretion disk is believed to play a fundamental role in such diverse systems as protostars, mass-transferring binary star systems, and the central engines of active galactic nuclei and quasars. The basic properties of accretion disks are reviewed by Pringle ( 1981 ). While there are many possible applications of a numerical MHD code to accretion disk theory, we are particularly interested in the nature of the torque that decouples the angular momentum from a fluid element and allows it to accrete. Ordinary molecular viscosity is an example of a hydrodynamic torque; in accretion disks, however, the molecular viscosity is infinitesimal. Magnetic fields are a promising alternative. We can

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

133

classify the influence of a magnetic field as an external or an internal torque. An example of an external torque is a rotating magnetized wind coming off of the disk and carrying away angular momentum. In order for this wind to be a significant sink of angular momentum, the magnetic fields must be moderately strong, on order of equipartition or stronger. An internal torque, on the other hand, results in the radial transfer of angular momentum out through the disk. The effective viscosity is proportional to the magnetic energy density, so this mechanism also seems to require near equipartition strength magnetic fields (o be effective. If the fields are weak they would have little influence. However, this conclusion is premature. It turns out that accretion disks are linearly unstable to an angular momentum transfer process in the presence of a weak magnetic field (Balbus and Hawley, 1991). The addition of a rotational velocity to the MHD equations has added a new wave family, an incompressible magnetorotational wave. What is remarkable is that there is a branch of the dispersion relation that is unstable for certain wavenumbers. The instability results from the ability of the magnetic tension force to transfer angular momentum from fluid elements in low orbits with high angular velocity to fluid elements in higher orbits with lower angular velocities. In orbital mechanics a body that is accelerated in the direction of its orbit will gain angular momentum and move to a higher orbit which results in a decrease in its angular velocity (Kepler's law). A magnetic field, by attempting to enforce corotation between such fluid elements (Ferraro's law) will decelerate the lower fluid element and accelerate the higher fluid element, transferring angular momentum outward in the process. Paradoxically, it is the weak magnetic field that is unstable, because the centrifugal forces resulting from the displacement must dominate over the ability of the magnetic field to hold the fluid elements together. If the field is strong the magnetic tension provides a restoring force that leads to rotationallymodified Alfv6n waves. From a linear analysis (Balbus and Hawley, 1991) we find that perturbations are unstable for wavenumbers k when ds22 (k - va) 2 < - dln"--R'

(2.12)

where/2 A is the Alfv6n velocity. Since the magnetic field strength enters only in combination with the wavenumbet, there will always be an unstable wavelength no matter how weak the field. (Strictly speaking this conclusion holds only down to wavelengths for which resistive effects become significant. In accretion disks, however, the resistive lengthscale and the corresponding field strengths are minute.) Thus whether or not magnetic tension lorces are important is not determined by the value of the integrated magnetic energy in comparison to the pressure or kinetic energies, i.e. /3. For any magnetic field strength there are length scales for which the magnetic tension forces are significant. In a Keplerian disk, with ,(2 ~ R -3/2, the disk will be unstable for all wavenumbers k such that ( k . VA) / ~ < "g/3.

(2.13)

For weak fields, i.e. VA << C~ << R/2 and/3 >> 1, there will be unstable modes with wavelength smaller than the disk scale height, H ,-~ cs/£2. If the magnetic field becomes too strong, the shortest unstable wavelength would not "fit" within one disk scale height, leading to a strongly magnetized disk. The existence of this magnetorotational instability indicates that MHD is essential for the structure and evolution of an accretion disk. Further, since the magnetic field enters the dispersion relation only as the dot product with the wavenumber, the instability is independent of the direction of the magnetic field. In practical terms this means that there are both axisymmetric and nonaxisymmetric forms of the instability, and their study requires a fully three-dimensional code. If we wish to simulate a realistic accretion disk we need an ideal MHD code that can handle dynamically important magnetic fields on both small and large length scales. Such field strengths range from /3 ~ I (equipartition) to ,8 >> l (weak). It is also of obvious importance to include the compressive, centrifugal, gravitational, and buoyant forces that are already known to be important in an accretion disk. Hence we require a code that can do compressible magnetohydrodynamics, with gravity and noninertial forces, in 3D. The code we

134

J.F. Hawley, J.M. Stone~ComputerPhysics Communications 89 (1995) 127-148

describe below was designed with these issues in mind, and has been used successfully for several investigations of MHD accretion disk physics.

3. The MOC-CT scheme 3.1. Overview

As the background discussion in Section 2 should make clear, there are many techniques for solving the MHD equations currently in use and in development. The scheme that we are using and that is described here is the Method of Characteristic Constrained Transport (MOCCT) scheme. We are presently using MOCCT for a variety of applications including, as described above, accretion disk studies. The MOCCT scheme is also undergoing continuing development. The MOCCT code incorporates the following general framework: ( 1 ) The equations are solved in the nonconservative form using operator-split, time-explicit finite differencing on a staggered grid, in which scalar variables are located at zone centers and vectors are located at zone boundaries. (2) The hydrodynamic terms are evolved following the already established procedures of codes such as ZEUS3D (Stone and Norman, 1992a; see also Norman and Winkler, 1986). We wish to maintain unchanged as much of the hydrodynamic algorithm as possible, and incorporate the new MHD terms in a manner that is consistent with the hydrodynamics. (3) We use the Constrained Transport (CT) formalism described by Evans and Hawley (1988; hereafter EH) which treats the magnetic fields as fluxes through zone faces at staggered grid locations. These fluxes are evolved using electromotive forces (EMFs) computed along the boundaries of the zone face. Through this procedure the V- B = 0 constraint is preserved to machine-roundoff precision in the numerical evolution. (4) The EMFs used in the induction equation are computed from time-advanced velocity and magnetic field estimates that are obtained with upwinded AlfvEn characteristic equations. (5) The magnetic pressure force terms of MHD are handled in a manner analogous to the hydrodynamic pressure terms, while the noncompressive transverse Lorentz forces are computed with the Alfv6n characteristic equations. We begin with a brief review of the hydrodynamic techniques. A detailed and specific discussion can be found in the cited literature, in particular Stone and Norman (1992a). The essential components to this approach are staggered grids and operator splitting. On a staggered grid the variables are located at different points in space, and possibly in time as well. (The leapfrog scheme provides one example of a grid that is staggered in both space and time.) The aim is to center naturally the derivatives that enter into evolution equations. For example, the evolution of the momentum is determined in part by pressure gradients, VP. By positioning the momentum between the spatial locations of the pressure one naturally obtains a centered derivative. Operator splitting is the technique in which the equations are broken down into successive steps, each of which operates on the output of the previous step. The full set of operations constitute the complete timestep. (See Bowers and Wilson, 1991, for a discussion of the advantages and drawbacks of operator splitting.) We divide the hydrodynamic equations into a Lagrangian and an Eulerian update. The Lagrangian phase consists of those terms that affect a specific fluid element, namely the acceleration terms, e.g., pressure and gravitational accelerations, and terms such as compressive heating and radiative losses. Since we solve the nonconservative internal energy equation we also include an artificial viscosity Q to account for shock heating. The Eulerian advection update moves the fluid elements through the fixed finite-difference grid according to the updated velocity field. For example, to evolve one-dimensional sound waves in a fixed grid one would update the momentum at spatial location i + 1/2 using the pressure gradient ~'P = (P/+I - P i ) / A x , solve for the updated velocity, and then evolve the pressure using this updated velocity in the compression term ~7 - v. This centers these derivatives in space and time resulting in a scheme that is in principle second order. All variables would then be advected through the grid according to the current velocity, thus completing the timestep. Operator splitting has a distinct advantage in that one can easily incorporate more physics into a scheme by adding additional terms to the Lagrangian update.

J.E Hawle),; J.M. Stone~Computer Physics Communications 89 (1995) 127-148

135

We note that there are additional differences of approach within this same basic framework, and these might vary from code to code or application to application. For example, one must choose either momentum S = pv or velocity v as the fundamental variable of evolution. In spherical or cylindrical coordinates one may also choose between angular momentum l, velocity v,b, or angular velocity ,(2, and that choice can have specific and important consequences (Norman, Wilson and Barton, 1980). In our experience, the best choice varies from problem to problem. The Eulerian phase of the timestep requires the computation of advection terms. We take as our model the continuity equation (2.1). The advection problem has been described in great detail in the cited literature to which the reader should refer for additional information. Here we review certain pertinent facts that will be used in developing the MHD algorithm. The finite difference form of the ID advection equation is At F, ,,+1/2 ,+1/21 p~i,+l = pn _ ~ b~ pVx)i+l/2 _ ( p v x ) i _ l / 2 j .

(3.1)

where Pi" represents the zone-averaged density in zone i at timestep n, and I'Pbx)i+l/2" , ,,+1/2 is a density flux through the i + 1/2 zone face at time level n + 1/2. The fluxes must be obtained from the values PI'; on a n+l/2

staggered mesh the velocity is already appropriately located in space and time so that v.~ i+1/2 is known. We make the approximation that the fluid is simply advected along at this centered velocity over the timestep. n+l/2 We can therefore predict the value Pi+l/2 by interpolating back along the advective characteristic a distance n~ I/2

l,~ ,+1/2 At and computing the average value of the density within this domain of dependence. This procedure depends upon establishing some functional representation of the density p ( x ) consistent with the zone average values that we have. The simplest scheme is to assume that the density is constant throughout the zone. The interpolation always yields the upwind value of the density, i.e., Pi for L'x > 0 and Pi+l for t'x "~ 0. This is the first order donor-cell scheme; donor-cell is noted for its great stability and large numerical diffusion. A signiticant improvement can be obtained by computing and making use of the slope A p / A x in the zone, i.e., a second-order representation. Lax-Wendroff is an example of one such second order scheme. The drawback of such higher order schemes is that they typically have substantial dispersion error that results in unacceptable high frequency oscillations in the advected function. This can be dealt with through the addition of an explicit numerical dissipation that filters out the short wavelength components. However, this runs contrary to one's aim in implementing a second-order scheme which is to reduce the numerical diffusion. The monotonic schemes such as those described by van Leer (1977; see also Sweby, 1984) provide a much better alternative to the "Lax-Wendroff plus viscosity" approach. In these schemes, a nonlinear term is added to the interpolation that damps precisely those 2Ax wavelength components that produce high frequency ringing. Ill van Leer's scheme the variables are represented by piecewise-linear functions that are allowed to be discontinuous at zone interfaces. For example, from x = Xi_l/2 to X = Xi+l/2, the density is given by the function p ( x ) = Pi + ( A p / A x ) i ( x

- xi),

(3.2)

where (Ap/Ax)i is a zone-centered slope. The nonlinearity enters through the monotoniciO' condition which requires that the value of the function at the edge of the zone, say p(Xi+l/2), not exceed the zone average value in the next zone, Pi+l. This condition prevents the development of new, numerically-generated local minima or maxima in the density. An example of one such monotonic slope is given by the harmonic average ol' two face-centered slopes,

(Ap/Ax)i =

(Pi+l--Pi)(Pi--Pi-1) 1 L(P,--~ --)5 +(-~--~//-~) '

2 [ ~x

(3.3)

At local maxima or minima, where (Pi+l - - P i ) ( P i - Pi-l) < O, the slope is set to zero. van Leer (1977) gives several suitable monotonic slopes which can be used instead of (3.3) as desired.

136

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

Once one has the density distribution (3.2), the flux term for the zone face at i + 1/2 in Eq. (3.1) will be the time-averaged value of the density that flows through that zone face in time At. This is obtained by integrating the density function from Xi+l/2 to xi+l/2 - VxAt. For the piecewise-linear function this is n+l/2 Pi+ 1/2 = Pi + I ( A x

VxAt)(Ap/Ax)i,

for vx > O, and n+l/2 = Pi+l -- I ( A x + Vxat) ( A p / A x ) i, Pi+l/2

(3.4)

for Vx < 0. One computes these values for all zone faces, and then carries out the evolution step (3.1) to complete the update of the continuity equation. For purposes of comparison, consider PPM. The PPM scheme represents variable distributions through the zones by third-order parabolas (Colella and Woodward, 1984). As discussed in Section 2.1, PPM is an example of a Godunov-type scheme in which all variables are zone-centered zone averages. To obtain the flux at the interface in the Eulerian version of PPM one interpolates back from the zone face not only along the advective characteristic v (as discussed above) but also along the characteristics v + c and v - c, where c is the sound speed. The averaged fluid values within these characteristic domains of dependence are the time-averaged fluid values that will interact at the zone face during the time interval At. These averaged values are used as input states to a nonlinear Riemann solver that resolves the input states into unique time- and space-centered values consistent with the equations of motion. Thus all the physics in a given problem must be incorporated directly into the solver. Conservative fluxes for all the equations are constructed from these values, and the update of Eq. (2.11 ) completes the timestep. Thus the Godunov-type scheme relies on a quite different approach from an operator-split finite differencing scheme such as we have outlined here. For a comparison of the performance of artificial viscosity schemes with conservative Riemann-solving schemes see Woodward and Colella (1984). As mentioned in Section 2.1, the Godunov approach could be employed for the full MHD equations and several groups are working on such an algorithm. However, the basic MHD Riemann solver is still under development, and, in any case, we wish to study systems that have more physics and physical wave families than these first MHD Riemann solvers are designed to handle. Hence, instead of following that track, we will try to make use of some of the ideas of a Godunov-type scheme to obtain a more accurate staggered-grid scheme. 3.2. Constrained transport

We begin the discussion of the MHD algorithm with the equation for the evolution of the magnetic field, i.e., the induction equation (2.4). For simplicity we will assume infinite conductivity, although nothing in the general CT approach depends upon that assumption. In this limit the induction equation is c~,B+ V × ( v × B) = 0 .

(3.5)

One of the challenging aspects of MHD is to evolve this equation subject to the solenoidal constraint on the magnetic field V . B = 0. In a finite difference calculation there are several ways one might choose to define the numerical divergence of a vector, only one of which can be exactly equal to zero. This can make it slightly ambiguous as to what one means by "divergence-less field" Brackbill and Barnes (1980) suggest that what is most important is that the magnetic field produces no nonphysical forces due to diverging magnetic fields. They show how this can be guaranteed by removing all terms that depend on the divergence of B from the force computation directly. However, the divergence free condition also implies a conservation law, i.e., the magnetic flux through a surface enclosing a volume should always be zero, and for accurate field evolution one would like to maintain this property regardless of how one computes the forces. This implies the need for a numerical

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127--148

137

conservation formulation of the induction equation that is explicitly multidimensional. This issue remains the subject of investigation, and here we summarize several of the options that have been considered. ( 1) Just evolve the magnetic field directly using the best available technique without regard to the constraint equation. The constraint will be violated because of truncation error, but this should vanish as Ax ~ 0. There are two concerns with this approach. First, one should check that the truncation error does not accumulate; it should remain bounded throughout the calculation. One should also determine the order of the error in V • B and how that order compares with that of the code as a whole. Second, because magnetic monopoles are completely unphysical, even if the divergence of the field is proportional to the truncation error its effect may be significant. The nonconservative form of the equations has the advantage of eliminating terms that depend on V - B prior to differencing. (2) One can evolve the vector potential, defined by B = V x A, instead of the magnetic field. This is advantageous primarily in two dimensions since it reduces the evolution of two magnetic field terms to the evolution of one scalar. The primary disadvantage is that the Lorentz force depends on second derivatives of the vector potential. This can be a real problem if the numerical scheme is based upon terms that are only continuous in first derivatives. Anomalous current reversals are quite easy to generate in such circumstances. (3) Another option is to fix up the magnetic field after an evolutionary step so as to force satisfaction of the constraint. One might, for example, evolve only two field components and use the constraint to get the third. This could be a problem if one is continuously dumping accumulated truncation error into one of the magnetic field components. Another approach would be to add an artificial "flux-cleaning" term using a scalar potential 05 defined by V205 + V - B = 0; the new field B = B - V'05 will be divergenceless (Brackbill and Barnes, 1980; Ramshaw, 1983). This approach was used by Zachary et al. (1993) in their 2D Godunov MHD scheme along with nonconservative momentum equations. They note, however, that there was no apparent difference between tests run with flux cleaning and tests without, suggesting that if the nonconservative momentum equations are used, monopoles may not be too severe a problem. (4) Another technique has been suggested by van Putten (1993). This entails replacing Faraday's equations plus constraints with a more general form for which the constrained MHD equations are a subset. An interesting facet of van Putten's equations is that they are fully general relativistic and may constitute the appropriate approach for the black hole electrodynamics problem. He has shown that for certain types of numerical schemes the constraint can be satisfied numerically throughout a computation (van Putten, 1994). (5) Finally, one can enforce the constraint on the numerical evolution equations through the use of staggered grid locations for the field and constraint-preserving finite difference equations. This guarantees that a unique finite-difference definition of the divergence will be zero throughout the evolution. This is the approach we have adopted. The details are given in EH who refer to this general method of numerical constraint preservation as "Constrained Transport" (CT). We now describe this approach in more detail. The CT approach starts by rewriting the induction equation into integral form using Stokes's theorem

c?tcl)s= f (v x B) dl,

(3.6)

8S

where q0s is the magnetic flux through surface S bounded by OS. If we consider a 3D finite difference zone aligned with the directions (x, y, z) labeled by indices (i, j, k), with edges of length Ax, Ay, and Az, then the magnetic flux piercing a given face, e.g. that located at Xi+l/2, is f

qbi+t/2 = [ Bxdydz. 2xy~z

(3.7)

J.F Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

138

|,j,k+l/2

/ ~y

Ex

(Bx)i÷ t/2,j,k

(x,,yj,~)~y)i,j+ I/2,k

AZ

~z

X Y

AY

"--...L_--- t f

AX

Fig. I. Locations of variables in the finite-difference zone. Scalars such as pressure and density are located at the zone center, labeled by integer values (i, j, k). The magnetic fields are regarded as area-averaged fluxes through zone faces at half-integer locations. Velocities are also located at zone faces. Electromotive forces (EMFs) are located along zone edges. The time evolution of the magnetic flux through a given zone face is obtained by a right-handed summation of the four edge EMFs.

Similar definitions can be made for the remaining zone faces. This is illustrated in Fig. 1. The scalar quantities are considered to be volume averages centered in the zone at the point labeled i,j,k. Vector quantities are located at the faces, with both velocity and magnetic components located at the same face, e.g., Vx and Bx are both located at zone x-boundary. The constraint simply requires that the total flux emerging through the 6 faces of the zone vanishes, ¢Pi+1/2 -- C1)i--1/2 ~- ¢Pj+l/2 -- % - 1 / 2 "1- ¢~k+1/2 -- ~ k - l / 2 = 0.

(3.8)

Substituting the definition (3.7) into (3.8) and dividing by the cell volume AxAyAz produces the finite difference definition of the constraint • • B = 0. By regarding the magnetic field as a flux through a zone face we recast the field's evolution entirely in terms of a line integral around that face, i.e, along the zone boundaries. The total rate of change of the flux through the surface equals minus the total electromotive force (EMF), -~" = p × B, around the contour 3S. An EMF at a given edge, say the one running in the z-direction, is defined, z(~+]/2) P

~zi+l/2,j+l/2

=

- -

/

(vxBy -- vyBx) dz.

(3.9)

z(k--l/2)

The extension to the other directions is obvious. The finite difference EMFs surrounding the surface at i - 1/2 are Cy i+1/2,k÷1/2' Ey i+U2,k-U2' Ez i+1/2,j+1/2 and Ez i+1/2,j-1/2. From these definitions we can write the finitedifference version of (3.6), taking the line integral around the surface in the right-hand sense, as ~.~_( ¢15n+1

At "ri+J/z - 4'~+1/2) = + ~ ' i+l/2'k+l/zA Y -- Ev i+l/2.k_l/2A y -- ~z i+1/2,j+1/2 A z "]'-~z i+1/2,j--1/2 A z,

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

139

for B~, and for By and Bz we have I t @r~+l

n

A~ ~ .]+1/2 -- ¢~j+1/2 ) = -- gX j+l/2,k+l/2 m x ÷ g.~ j+l/2,k-I/2 A x

÷ ~Z i+1/2,j+1/2 AZ --~z i-1/2,j+1/2 mg, L/,q~n+l At ~k+l/2

_ ~ibn k-~ 1/2) = -- gY i+l/2,k+l/2AY ÷ gY i-1/2,k+l/2Ay

÷ £_~,j+l/2,k+l/2Ax -- £.,: j_I/Z,k+I/zAX.

(3.10)

Division of both sides of these finite-difference equations by the appropriate area, e.g. A y A z for q~i+~/2, yields the CT evolution equation for the magnetic fields. If one substitutes the finite-difference evolution equation (3.10) into the constraint equation (3.8), one sees that numerical conservation of flux follows naturally since each EMF edge contribution appears twice in the sum, once with positive sign and once with negative. If the initial field satisfies the constraint, so will the field generated by (3.10). EH were concerned with obtaining such a conservative formalism for general coordinates with an eye to general relativistic calculations. This is achieved by evolving not the magnetic field directly, but instead the magnetic vector density/3 i = Bi/x/~, where x/Y is the determinant of the three metric (e.g., r 2 sin 0 in spherical coordinates). By subsuming the coordinate information in the definition of the vector density we reduce the induction equation in curvilinear coordinates back to the simple Cartesian form described here. This approach will work so long as the zone volume is well approximated by v / ~ A x A y A z , where Ax, Ay, and AZ are general curvilinear coordinate elements. Near a coordinate singularity, e.g. near r = 0 in spherical coordinates, additional care is required if the truncation error associated with the geometry is to be kept to a minimum (see, e.g., M6nchmeyer and Mtiller, 1989). Note that the addition of a finite resistivity can be accomplished by computing the current J along the zone edges from the curl of the magnetic field (already naturally centered), and subtracting r/J from the appropriate EMFs (where "r/ is the magnetic diffusivity). 3.3. ID M O C - C T Up to this point we have not specified how to obtain the EMFs. As long as one stays within the CT framework, any technique could be used to compute EMFs and the constraint would still be preserved. Obviously there will be good and bad procedures for computing EMFs as determined by numerical stability and consistency requirements, as well as the size of the truncation error. To begin our search for a suitable EMF we consider a 1.5 dimensional (1.5D) system in (x, y) where Oy = 0 and assume a uniform Bx field. For this simple system, OtBx = 0. The induction equation for B,. can be written OlBy = BxOxCy - 0x(cxBy).

(3.11 )

The second term is just an advection term corresponding to the physical transport of the transverse field by the velocity Vx. The first term, however, is different. It represents the shearing of the Bx field by a gradient in C'y. It is this first term which, when coupled to the transverse Lorentz force, leads to the propagation of Alfvdn waves at the Alfvdn velocity Ca = Bx/x/4zrp. Since the EMF is located at an advanced time level, we must obtain estimates for v and B at the n + ½ time level. If the induction equation consisted of nothing more than the advection of magnetic field by fluid motion we would simply replace the density in (3.4) with By and vx to obtain average values for the magnetic field and transverse velocity within the domain defined by the characteristic speed Vx. This was the approach used by EH lbr their passive field advection test problems, and in the limit Vx >> Ca this is entirely adequate. However, as we shall see, when evolving systems in which the magnetic field is dynamically important we need to take into account the Alfv6n waves produced by the transverse Lorentz force.

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

140

The Lorentz force is ([ x7 × B ] × B )/4zr. This in turn can be subdivided into the magnetic pressure and the transverse force terms,

( [~7 x B ] x B ) / 4 ~

=

- ~ ( B 2 / 8 ~ ) + (B. "~B)/4~r.

(3.12)

It is the second term on the right, the transverse tension force, that is responsible for Alfv6n waves. Ignoring for the moment the pressure terms, the momentum equation,

8tUy = 4~pSxBy - 8x(VxVy),

(3.13)

consists of the tension force and an advection term. The advection term is familiar, but the transverse force term is why MHD has fundamentally different properties from hydrodynamics. This implies that caution is in order when we try to extrapolate from past experience with numerical hydrodynamics. To develop a scheme we need to start with simple applications, and the 1.5D system makes an admirable testing ground for MHD algorithms. Stone et al. (1992) describe a number of such 1.5D tests, including a rotating solar wind problem, the Brio and Wu (1988) MHD shock tube, and pure Alfv6n wave propagation problems. It is the latter which are most enlightening in developing MHD algorithms. In a full MHD problem such as the shock tube or the solar wind evolution, hydrodynamic and magnetic pressures are more important than the magnetic tension force, and it is possible to solve those test problems to apparently acceptable accuracy without ever rigorously testing the Alfv6n terms. In the limit of the pure Alfv6n wave problem, however, these potentially troublesome terms from the MHD equations are isolated for study. Indeed, our first attempts to create a finite difference MHD scheme, based on a straightforward extension of the original CT scheme of EH, had difficulty with the pure Alfv6n wave test problem (Stone and Norman, 1992b), although it handled the shock tube and MHD wind problems reasonably well. C. Evans (private communication) showed that it was possible to understand these initial difficulties from a v o n Neumann analysis of the Alfv6n wave system. Consider this limit in which the MHD equations reduce to a simple wave equation of the form

3,(b = 3xW, 3tw = 8xq~,

(3.14)

assuming unit Alfv6n velocity. One obvious approach to differencing this system is with operator-splitting and a staggered mesh so that (fy~:r+ 1

j

n

-~b~ = o ' ( w

n+2l-

n+2I-

l - w . 1), J+~ J-~

(3.15)

and W

1 -- W

.i+ ~

J+

1 ----"O ' ( t ~ J + l

-- t~;),

(3.16)

where o- = A t / A x is the Courant number. This is the staggered leap-frog scheme and it is stable, with amplification factor equal to unity, so long as the Courant number is less than one. However, the standard von Neumann analysis yields a numerical dispersion relation o~(k), where w and k are the frequency and wavenumber in which the group velocity d w / d k goes to zero as the wavelength approaches 2Ax. This means that an Aifv6n wave computed with the leapfrog scheme will be very dispersive, leaving high frequency components in its wake as it travels through the grid. Things are worse if the terms & and w are located at the same point in space, as they are in the centering described above for the MOCCT scheme. The scheme is still stable but the dispersion relation has a maximum at the wavenumber kmax/2, where kmax = 7T/Ax, and then turns over. Thus, for the shortest wavelength on the grid, the group velocity is minus the physical propagation speed. In the absence of sufficient dissipation

J.E Hawley, J.M. Stone~ComputerPhysics Communications 89 (1995) 127-148

141

to eliminate the short wavelength components, significant ringing and numerical noise will quickly develop throughout a calculational grid. Waves propagating off the grid through boundaries pose a great problem since truncation errors at boundaries often generate noise, and this dispersion relation guarantees that this noise will propagate back onto the grid. This explains the early difficulties with Alfv6n waves, such as the test shown in Fig. 18 of Stone and Norman (1992b). It is worth noting that these undesirable numerical effects were present in other test problems where pressure terms were more important than Alfv6n waves, but at sufficiently subtle levels that they were not easily noticed (although they could be seen clearly by plotting the error rather than the variables themselves). This demonstrates the importance of continual testing with many different problems in different physical limits. Dispersion errors such as these can be dealt with through the use of an aggressive numerical diffusion term. Many of the finite-difference codes used in previous work depend on just such diffusion terms. This may not be the best strategy, however, if one wishes to evolve wave phenomena over a fairly wide dynamic range. We would like to work toward a scheme in which just enough damping is applied to prevent these dispersion errors, but not so much as to compromise the solution. This was, of course, the idea behind such algorithms as van Leer's monotonicity condition for hydrodynamics. Clearly we must take greater care when formulating a scheme to compute both the EMFs and the Lorentz force. Let's return to the EMF problem and consider what will happen over time At at the zone edge located at xi+l/2 where the EMF is computed. First the magnetic field will be swept past by any velocity field that is present (advection). In the usual way the time average over At of the advected field will be the average field within the domain of dependence determined by Vxat (as in [3.4]). However, simple advection is not the only way that the magnetic field evolves. There will also be left- and right-traveling Alfv6n waves if there is a transverse field to support such waves. These waves will have domains of dependence determined by wave speeds (u + VA)At and (v - UA)At. Using (3.4) with magnetic field substituted for density and the left- and right-traveling AlfvEn wave speeds for the velocity, we can compute the average magnetic field values over these domains for the time interval At. These are designated By and B +, respectively. Similar averages can be obtained for Vy, i.e., Vy and u +. Now the solution to (3.11) and (3.13) is given by the Alfv6n characteristic equation Dt;.~, i Dt

1 DBy _ O, 4 x / - ~ Dt

(3.17)

where

D Dt

c9 Bx - cgt + (Ux 5: ~

)a_a_ ax

At xi+l/2 where we wish to compute the EMF, this equation must be satisfied for both the left- and right-going Alfv6n wave over the timestep At. There are two equations for the unknowns By and v~,, the time-advanced magnetic field and velocity. Along the + characteristic B v. and Uy evolve from B + and Oy+to B~,* and v~,,*"along the - characteristic B~7n and o~7n also evolve to By and Vy. Using Eq. (3.17) we can write

e ,• =

[B:+ +

a; +

Vy) ] ,

(3.18a)

and V•y = 7,[+Vy +

Vy + ( B y+ - B y )

.

(3.18b)

Thus the time-averaged values of the magnetic field and velocity carried to the zone face by the fluid flow and the Alfv6n waves, along with the characteristic equations, give us a prediction for By and Vy at an advanced time level that will, by design, account for the Alfv6n waves.

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

142

Here we describe in detail the steps in the MOCCT algorithm for the 1.5D problem in (x,y), with y as the symmetry direction. For convenience we have incorporated the factor of ~ into the definition of the magnetic field in all following equations. The density p, By field and Uy velocity are all located at zone centers designated by the integer xi locations. (In 2D By and Uy will be located at the half-integer Yj+1/2 zone faces.) The velocity Vx and magnetic field Bx are located at half-integer zone boundaries xi+v2. Note that in the Cartesian 1.5D problem Bx is a constant in space and time. The induction equation for By is relatively simple, involving only one EMF. To compute that EMF £i+1/2 at location xi+l/2, first compute the characteristic wave speeds

C%1/2 = Ux i+1/2

--

Bx i+l/2/v@i+l/2,

Ci;1/2 = Ux i+1/2

@

Bx i+l/2/V/-'Pi+l/2 ,

where

(3.19)

Pi+l/2 is a suitable average to the location Xi+l/2, e.g.,

Pi+l/2 =- ½(Pi+I

+

(3.20)

Pi).

In multidimensions both Bx and Vx would need to be averaged in the transverse (in this case y) direction as well. Next compute the time-averaged values of the y field and velocity within the domain of dependence for these two wave speeds using the piecewise-continuous van Leer functions. For example, consider B +,

B,+=Byi+½(Ax-c+At)(~)

i

,

if c + > 0, and B ,+ = B r i + j "

½( Ax + c+At) ( ABy )

(3.21)

\ AX //i+1 '

if c + < 0. The t e r m ( A B y / A x ) i is the monotonic van Leer slope centered on xi as defined, for example, in (3.3). Similar time-averages are obtained for By and v~. Once all four input averages are computed, we can solve the characteristic equation for By and Cy at zone face i + 1/2,

8:* = ~' (SCNv~l v,+-vTl + le,+ + 8 7 1 ) ,

(3.22a)

v~*,= ½ (SGN [B~+, - Byl/V'7 + IVy+ + Vy]),

(3.22b)

and

where SGN is the sign of Bx. Note that this algorithm reduces to simple advection in the limit that VA << Vx. Finally, construct the EMFs from the starred values, --Ei+I/2

= Vx i+l/2By - VyBx

i+1/2"

(3.23)

The magnetic field is then updated via the CT step, By n+l = By n ..]_At (Ei+I/2 -- Ei-1/2) lax.

(3.24)

Because both the induction equation and the Lorentz force are implicated in the production of undesirable dispersion error, it is important to carry out a similar procedure for the transverse Lorentz force. We again require magnetic fields at a time advanced level, now for the momentum difference equation,

Syin+l

= S n

yi + A t B x ( B y i*+ l / 2 - B y i -*u 2 ) / A x "

(3.25)

J.F. Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

143

Here the B~ terms are those responsible for Alfv6n waves and again we use the Alfvdn characteristics to improve our estimate. Note that in problems more complicated than the simple 1.5D system we must use an appropriate average of Bx after it is updated in the induction equation routine. The By values for (3.25) are obtained as described above using updated magnetic fields and characteristic speeds of ±Bx/x/-P. Note that the characteristic speed no longer includes the vx velocity because we regard this update as a Lagrangian step; the advection of momentum and magnetic field have already been performed for this timestep. The magnetic pressure term poses no special problem. In the 1.5 dimensional case the By field provides an acceleration term for the x momentum Sx, S n~l i~1/2 --= S.~ ,,i+1/2 - A t ~ ( O y i + l

-Oyi)/Ax,

(3.26)

where ~ is an appropriately centered spatial average of By at location i + 1/2. In higher dimensions the difference in B). will also have to be averaged in the y-direction. The pressure acceleration is a Lagrangian step and is carried out after momentum advection and the evolution of the magnetic field by the induction equation. This is required for proper centering in time of the terms in the wave equation. The 1.5D version of MOCCT has been tested extensively and many of these tests are described in Stone et al. (1992). The scheme reproduces all the features of the Brio and Wu (1988) magnetic shock tube. Comparing our results with those of the MHD Godunov-type codes we find that MOCCT has larger errors in the post-shock values; this is consistent with the performance of an artificial viscosity code versus a Godunov code on ordinary hydrodynamic shock tubes. The Alfvdn wave tests are quite satisfactory. Alfvdn waves are propagated on the grid in both stationary and moving backgrounds with, if anything, less dissipation than seen in the passive advection of a pulse with the van Leer algorithm. Stone et al. (1992) compute a rotating, magnetic solar wind problem and a magnetic-braking problem and compare with analytic solutions. For even moderately resolved grids the relative error is on order 0.1%, and Stone et al. (1992) demonstrate convergence toward the analytic solution with increasing resolution. 3.4. Multidimensional M O C C T

In principle, the extension of the 1.5D MOCCT scheme to muitidimensions is straightlorward. In either 2D or 3D the MOCCT step is done in planar slices. For example, the computation of the z-component of the EMF involves only the x and y components of B and v. This is illustrated in Fig. 2. Consider a 3D system in Cartesian coordinates (x, y, z). A given EMF, say £z i+l/2,j+l/2,k, is computed in the (x, y) plane located at zk and centered at the zone corner xi+v2, Yj+J/2. First, as described above, one computes the B~ and v~ values using u~ and the Alfvdn speed averaged to the zone corner. Then one repeats the procedure, computing B~ and ;.,. •* from the t,,. velocity and Alfvdn speed similarly averaged to the zone corner. One can then construct an EMF, - £ : i~J/2,j+l/2.k = v~B,*, - v3*Bx.

(3.27a)

Similarly one defines --'~x i,.i+l/2,k+l/2 = urO,. - vzB~.,

(3.27b)

* * - ~ y i+l/2,j,k+l/2 = Uz*Ox* - VxBz,

(3.27c)

and

each computed in the appropriate 2D plane. These can now be inserted into the CT evolution equation (3.10). The extension to multidimensions of the 1D scheme for the magnetic force is also obvious. For example, the acceleration of the z momentum Sz by the transverse magnetic forces is given by

J.E Hawley,J.M. Stone~ComputerPhysics Communications89 (1995) 127-148

144

(~,) (e,)t,j+./...

AY

(~)tn+l/sj+l/~k)

pl P

(B,)o.,/~j.~)

] (x,, ~, z,)

f I I f I

f~x

r

Fig. 2. A two-dimensional (x, y) slice showing the location of variables used in the computationof the z componentof the EME Sz

i,j,k+l/2 n+l

Sz i,j,k+l/2 -~- At [ -~v(Bz. n

=

i,j+l - Bz ,,j*.)/Ay + Bx(Bz i+l,j *

- -

*

Bz i,j)/Ax ] , *

- -

(3.28)

where Bx and By are four-point averages to the spatial location of Sz of the magnetic field after an update by constrained transport. The starred magnetic field values are each computed with a 1D characteristic computation as described above. Similarly, the magnetic pressure acceleration is

Szn+,

,

i,j,k+l/2 = Sz i,j,k+l/2 -- A t

[-~xfABx)+-~y(ABy)] ~x Az } \ mz }J '

(3.29)

where the differences and averages are all computed using the updated magnetic field values. The multidimensional MOCCT scheme has undergone testing on a variety of problems, including wave advection tests (both pressure and Aifv6n) and Brio and Wu (1988) shock problems in all directions. However, these tests are just variations on the 1D tests. True multidimensional tests are more difficult to come by. The solar coronal ejection problem in Stone et al. (1992) is one example. The simulations described in Hawley and Balbus ( 1991 ) provide another test by comparing numerical growth rates for the MHD instability in rotationally supported systems with those obtained analytically from perturbation theory. The MOCCT scheme does well with all of these these problems. We have found, however, a test problem on which the scheme as described above has difficulties. This is the evolution of a current sheet. A current sheet is present whenever the magnetic field changes discontinuously. The existence of current sheets is a problem in general for finite difference MHD, as is any discontinuity, since they cannot be represented accurately in finite-sized zones. In the absence of resistivity a current sheet should not diffuse away, nor should field lines reconnect across it. However, since a finite difference code cannot resolve structures on scales smaller than a grid zone, numerical reconnection will occur. It is worth noting that the finite grid size sets an effective reconnection length scale even in the ideal MHD limit. Because magnetic flux is viewed as a finite integral over a surface, one loses information each time oppositely-directed fields are carried into the same grid zone. This is anomalous numerical reconnection and the magnetic energy thus lost does not go into heat via Joule heating. Some of this magnetic energy is carried off in waves, however, so it need not be entirely lost to the system. The corresponding numerical magnetic Reynolds number is roughly equal to the number of grid zones used to resolve a given spatial dimension,

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

145

although the effect is not strictly a magnetic resistivity; stationary fields will not diffuse away, for example. This numerical Reynolds number does, however, set a limit on how small a physical resistivity can be modeled with a finite-difference code. The diffusion length scale set by the physical resistivity must be greater than the size of the grid zones, i.e., it must be resolved. Given the relatively crude resolution currently possible in 3D astrophysical simulations, and the very high magnetic Reynolds numbers associated with most astrophysical plasmas, there is no need to include a physical resistivity term. The current sheet test problem is designed to see what an algorithm will do with a perturbed current sheet. Although we do not have an analytic solution for this test problem, it is nevertheless useful to test the response of the code to the presence of a current sheet. A typical current sheet test is run in two dimensions using a periodic, square grid upon which there is a uniform magnetic field that discontinuously reverses direction at some point. This current sheet is then perturbed by adding a small, spatially-varying transverse velocity. In one test we used a 64 x 64 zone grid over an x- and y-domain from 0 to 1. The background pressure and density are constants ( P = 10 -6, p = 1). From 0 < x < 0.5 we set By equal to a large positive value (we used lv/]-0P), and from 0.5 < x < 1 we set By negative. Because of the periodic boundary conditions there are two current sheets. These initial conditions are perturbed with a periodically varying sub-Alfv6nic transverse velocity, Vx (= 10-4By sin(27ry)). We then observe how the current sheet evolves. One expects it to oscillate with an Alfv6n frequency (so long as the amplitudes are very small). Any numerical reconnection along the current sheet will show up as loss of magnetic energy and/or the generation of waves at the reconnection site. We found that numerical reconnection occurred in amounts that depend upon such factors as the relative strength of the magnetic field and the grid resolution. However, some of the reconnection events were accompanied by the generation of a strong local maxima in the EMF at one grid corner location which, in turn, created a strong field loop around that corner. Sometimes this strong field loop lingers and sometimes it is promptly destroyed by reconnection, but even then it produces anomalous sound waves that add to the numerical noise at the current sheet. These strong field loops can also significantly reduce the timestep, a serious problem for multidimensional simulations since one then spends inordinate time evolving an anomalously strong field loop. These difficulties can be avoided by adding a large enough resistivity so that the resistive length scale is larger than a grid zone. The current sheet is then no longer a discontinuity. If this were done, however, we would no longer be computing the ideal MHD equations. In fact the goal for astrophysical MHD is to reduce the resistivity as much as possible. Hence we again seek a solution to this difficulty that does not involve increasing the numerical diffusivity. A problem similar to that seen in the current sheet test problem was observed by D. Clarke (private communication) in tests of 3D MOCCT with strongly shearing MHD jet flows in an extreme limit: fl ~ 10 l°. In Clarke's case, this was an especially severe problem, since these field generation events generated anomalous strong field loops in an otherwise utterly weak field domain. Clarke was able to identify the source of the difficulty as stemming from two things: (1) the presence of discontinuities at the zone boundaries due to the piecewise-linear monotonic representations of the functions, and (2) a potential inconsistency between the velocity and magnetic field used to estimate the Alfv6n wave domain of dependence, (B,F), and the time-advanced values, (B*, c,*). Let's consider these two points in more detail. First, the piecewise linear functions representing the magnetic and velocity fields are often discontinuous at the zone boundaries. Physically such discontinuities will generate an Alfv6n wave which will, in turn, be included in the time-advanced field and velocity estimates. The strength of that wave, however, will depend on the magnitude of the transverse magnetic field; a negligible field will produce a negligible effect. Now recall that when estimating the time-advanced magnetic field, e.g. B,*., we used an Alfv6n characteristic cone determined by the transverse magnetic field strength and velocity, obtained using averaged values, v-~-Bx/v/-fi. Similarly, B--~yand ~ are used when computing B x and vx. However, each starred field component, e.g., B~*, is a solution of the characteristic problem only with the transverse field estimate, Bx, that was used in establishing the characteristic cone. This Bx need not be consistent with the transverse starred

146

J.F. Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

field estimate B x. In most cases it is, but near current sheets or near strong shear layers they may be quite different. To make this more concrete, consider what happens near the current sheet in the problem described above. There is a strong discontinuity in By. If Bx/v/-fi > vx then the plus and minus characteristics domains straddle the discontinuity and there will be a large difference between B~+ and By due to the discontinuity. By the characteristic Eq. (3.18) a significant Vy will be generated. This is the correct solution: a rotational discontinuity in the field will produce an Alfv6n wave; however, the amplitude of that wave will depend on Bx. If, however, in the computation of B x through the sweep in the y-direction one obtains a dramatically different value of Bx from that used to compute the Alfv6n speed in the x-direction, then the two values B x and v~*, will not be consistent with one another and an EMF constructed out of their product will be inaccurate. This is fundamentally a multidimensional problem; in simple 1.5D only one Alfv6n wave direction exists and there is only one characteristic equation to solve. The root of this difficulty lies with the use of a ID directionallysplit procedure to create components of an inherently multidimensional operation, i.e., evolving the induction equation subject to the divergence-free field constraint. Clarke addresses this difficulty with an algorithm he calls Consistent MOCCT. The characteristic problem is done implicitly in two directions at once so that the starred values are the same as the values used to compute the domain of dependence in the time-averaging. (In the present notation this means that B x and B~ are the same.) Thus the two 1D MOC procedures are replaced with one 2D implicit MOC calculation. From tests and simulations, Clarke (private communication) has demonstrated that this procedure solves the anomalous EMF problem. It has the drawback usually associated with any implicit scheme of much greater complexity and computational expense. We have developed an alternative procedure that remains time-explicit, but corrects the difficulty by using our I D sweeps consistently. To demonstrate this approach, consider the computation of Ez from x- and ycomponents of the velocity and magnetic field. To obtain B**. and Vy we need characteristic cones projected in the x direction back from the zone corner at velocities ~ + B--~/v/-fi and b T - ~ - ~ / v ~ " To estimate these velocities we assume that a good first guess of Fxx and Bx at the zone corner can be obtained by using velocities that have already been advected in the Eulerian transport step and assuming passive advection of the magnetic field by that velocity. (This is in contrast to the simple spatial averages used originally.) Essentially these estimates are those used in the passive advection formulation of CT (EH). The advected values are then used to formulate the cones for the characteristic problem in the y-direction. The characteristic problem is carried out as usual, producing B~*. and v~*..The process is repeated for v~ and B~ using advection-estimates ~ and BT'. Finally we form the EMF by combining starred estimates only with the values with which they are consistent. The two 1D EMFs are averaged to obtain a final EMF, ( Vx By + v-TB ~*

v,, Bx - v~B~ ).

(3.30)

The remaining EMF components are similarly dealt with and the constrained transport proceeds as usual. In the torsional MHD force routine the starred values were already coupled directly with the magnetic field value used to compute the Alfv6n speed, hence no changes are required there. We can now summarize the full timestep in the multidimensional MOCCT scheme: ( 1 ) Eulerian transport of momentum, energy and density, (2) constrained transport of the magnetic fields, (3) transverse Lorentz force using updated fields, (4) acceleration by gas and magnetic pressure, followed by the velocity update, and (5) PdV work and viscous heating. As usual, the timestep is limited by the Courant condition, which now must include the fast magnetosonic speed as well as the sound speed. Undoubtably there are still improvements that could be made within this framework, but in tests this algorithm performs very well. The current sheet test problem is significantly improved by using the EMFs as defined by (3.30). There are no anomalous current spikes or fields anywhere along the current sheet. There are few other apparent differences from the original formulation for the EMF, (3.27). The only change observed in MHD

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

147

shock test problem is that the post-shock oscillations were slightly smaller in the original form of the scheme. The only proof of the procedure, however, will be further applications and tests.

4. Applications and directions In this paper we have described the current implementation of the MOCCT finite difference algorithm for multidimensional MHD. We believe that the 3D MOCCT scheme is a significant improvement over MHD schemes based on the donor cell or Lax-Wendroff techniques. We have been able to simulate systems both with weak fields B >> 1 and relatively strong fields/3 < 1, as well as models where the magnetic field is found in only part of the system, i.e. hydrodynamic regions and MHD regions in the same problem. To date the MOCCT scheme has been employed in a number of astrophysical applications including MHD jets, and shocks hitting magnetized interstellar clouds. In our own work, MOCCT has been used for studies of magnetized accretion disks. The behavior of the MHD accretion disk instability has been investigated through 2.5D axisymmetric simulations (Hawley and Balbus, 1991; 1992). The results were compared directly with analytic linear perturbations calculations. Hawley and Balbus (1992) found excellent agreement between numerical and analytic growth rates. Fully 3D simulations of the instability have also been carried out (Hawley, Gammie and Balbus, 1995). The resolution that can be used for such simulations is necessarily crude, but there is still good agreement between analytic and numerical growth rates and wave frequencies. Stone and Norman (1994) have simulated an accretion disk in 2.5D with axisymmetry. The disk is initially in Keplerian equilibrium and threaded by a vertical magnetic field. In the weak field limit (/3 > 1) the disk is subject to the MHD instability of Balbus and Hawley (1991). The transfer of angular momentum within the disk quickly leads to infall. In the strong field limit (/3 < 1) the predominate means of angular momentum transfer is up and out of the disk through torsional Alfv6n waves. This too leads to a rapid disk infall. There are several future directions for MHD code development that will be important for astrophysical simulations. One of the advantages of the MOCCT scheme is that it is relatively straightforward to add additional physics. For example, simple radiative cooling could be added after the P d V work step. Another example would be finite resistivity. To date we have used only the ideal MHD equations, i.e., infinite conductivity, in our astrophysical applications. We have not included any explicit resistivity, nor is any needed in order for the code to perform well. Of course, there is no difficulty in adding finite resistivity to this formulation and we have run such cases as test problems. Finite resistivity may be appropriate when considering the evolution of protostellar disks. An important astrophysical extension of the present algorithm would be to two-fluid systems consisting of both ions and neutrals. Such a system is appropriate to describe dense regions of cold molecular clouds where the ionization fraction is small. Two fluid systems of ions and neutrals have several additional wave families whose properties will have to be accounted for properly. Another challenging task is to develop techniques that will work both for MHD regions and electro-magnetic regions. An example application is to the study of accretion disks around black holes. Magnetic fields can, in principle, tap directly into the rotational energy of a black hole. The process is primarily electrodynamic, although an accretion disk is required to anchor the field around the hole. Such a model has been invoked as a possible scenario for radio jet production. For this we would need a fully relativistic algorithm that included displacement currents. Such a code is still beyond the current state of the art, but we are hopeful that the efforts of the many people now investigating astrophysical MHD algorithms will point the way toward this future.

148

J.E Hawley, J.M. Stone~Computer Physics Communications 89 (1995) 127-148

Acknowledgements We w o u l d like to a c k n o w l e d g e D a v i d Clarke, Charles Evans, and M i k e N o r m a n for contributions to this a l g o r i t h m d e v e l o p m e n t effort, and Charles G a m m i e for a careful reading o f the manuscript. This w o r k is partially supported by N S F grant P H Y - 9 0 1 8 2 5 1 and N A S A grants N A G W - 1 5 1 0 and N A G W - 2 3 7 6 . C o m p u t a t i o n s were carried out on the Cray YMP, Cray 2, and C M 2 systems o f the National Center for S u p e r c o m p u t i n g Applications.

References S.A. Balbus and J.E Hawley, 1991, Astrophys. J. 376, 214. R.L. Bowers and J.R. Wilson, 1991, Numerical Modeling in Applied Physics and Astrophysics (Jones and Bartlett, Boston). J.U. Brackbill and D.C. Barnes, 1980, J. Comput. Phys 35, 426. M. Brio and C.C. Wu, 1988, J. Comput. Phys. 75, 400. P. Colella and P.R. Woodward, 1984, J. Comput. Phys. 54, 174. W. Dai and P.R. Woodward, 1994a, J. Comput. Phys. 111,354. W. Dai and P.R. Woodward, 1994b, J. Comput. Phys. 115, 485. C.R. Evans and J.E Hawley, 1988, Astrophys. J. 332, 659 (EH). J.A. Fedder and J.G. Lyon, 1987, Geophys. Res. Lett. 14, 880. S.K. Godunov, 1959, Matematicheskii Sbornik 47, 271. D. Gottlieb, 1972, J. Numer. Anal. 9, 650. J.E Hawley and S.A. Balbus, 1991, Astrophys. J. 376, 223. J.E Hawley and S.A. Balbus, 1992, Astrophys. J. 400, 595. J.F. Hawley, C.F. Gammie and S.A. Balbus, 1995, Astrophys. J. 440, 742. A. Kageyama, K. Watanabe and T. Sato, 1992, J. Geophys. Res. 97, 3929-3943. Y.Q. Lou, R. Rosner and P. Ulmschneider, 1987, Astrophys. J. 315, 349. R. Mtinchmeyer and E. Mtiller, 1989, Astron. Astrophys. 317, 351. A. Nordlund, A. Brandenburg, R.L. Jennings, M. Rieutord, J. Ruokolainen, R.E Stein and I. Tuominen, 1992, Astrophys. J. 392, 647. M.L. Norman, J.R. Wilson and R.T. Barton, 1980, Astrophys. J. 239, 968. M.L. Norman and K-H.A. Winkler, 1986, in: Astrophysical Radiation Hydrodynamics, K-H. Winkler and M. Norman, eds (Reidel, Dordrecht) p. 187. T. Ogino, R.J. Walker, M. Ashour-Abdalla and J.M. Dawson, 1986, J. Geophys. Res. 91, 10029-10045. G.J. Philips and J.J. Monaghan, 1985, Mon. Not. R. Astron. Soc. 216, 883. J. Pringle, 1981, Ann. Rev. Astron. Astrophys. 19, 137. J.D. Ramshaw, 1983, J. Comput. Phys. 52, 592. P.J. Roache, 1972, Computational Fluid Dynamics (Hermosa, Albuquerque). K. Shibata, Y. Uchida, 1985, Publ. Astron. Soc. Jpn 37, 31-46. J. Stone, J.F. Hawley, C.R. Evans and M. Norman, 1992, Astrophys. J. 388, 415. J. Stone and M. Norman, 1992a, Astrophys. J. Supp. 80, 753. J. Stone and M. Norman, 1992b, Astrophys. J. Supp. 80, 791. J. Stone and M. Norman, 1994, Astrophys. J. 433, 746. G. Strang, 1968, SIAM J. Numer. Anal. 5, 506. P.K. Sweby, 1984, SIAM J. Numer. Anal. 21,995. B. van Leer, 1977, J. Comput. Phys. 23, 276. M. van Putten, 1993, J. Comput. Phys. 105, 339. M. van Putten, 1994, SIAM J. Numer. Anal., in press. K. Watanabe, and T. Sato, 1990, J. Geophys. Res. 95, 75-88. P. Woodward and P. Colella, 1984, J. Comput. Phys. 54, 115. A.L. Zachary and P. Colella, 1992, J. Comput. Phys. 99, 341. A.L. Zachary, A. Malagoli and P. Colella, 1994, SIAM J. Sci. Comput. 15, 263.