A simple unsplit Godunov method for multidimensional MHD James M. Stone * , Thomas Gardiner 1 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, United States article info Article history: Received 3 April 2008 Received in revised form 20 June 2008 Accepted 24 June 2008 Available online 3 July 2008 Communicated by G.F. Gilmore PACS: 95.30.Qd 95.75.Pq Keywords: Hydrodynamics MHD Methods: numerical abstract We describe a numerical algorithm based on Godunov methods for integrating the equations of com- pressible magnetohydrodynamics (MHD) in multidimensions. It combines a simple, dimensionally- unsplit integration method with the constrained transport (CT) discretization of the induction equation to enforce the divergence-free constraint. We present the results of a series of fully three-dimensional tests which indicate the method is second-order accurate for smooth solutions in all MHD wave families, and captures shocks, contact and rotational discontinuities well. However, it is also more diffusive than other more complex unsplit integrators combined with CT. Thus, the primary advantage of the method is its simplicity. It does not require a characteristic tracing step to construct interface values for the Rie- mann solver, it is straightforward to extend with additional physics, and it is suitable for use with nested and adaptive meshes. The method is implemented as one of two dimensionally unsplit MHD integrators in the Athena code, which is freely available for download from the web. Ó 2008 Elsevier B.V. All rights reserved. 1. Introduction Numerical methods for compressible magnetohydrodynamics (MHD) have become indispensable for the study of a wide variety of problems in astrophysics, such as turbulence and angular momentum transport in accretion disks (Balbus, 2003), the pro- duction and collimation of jets and winds (de Gouveia Dal Pino, 2005), and the properties of supersonic turbulence in the ISM (McKee and Ostriker, 2007). Although a wide variety of numerical algorithms for integrating the equations of MHD are possible, there has been a substantial amount of effort in recent years to extend Godunov methods to MHD (Dai and Woodward, 1998; Ryu et al., 1998; Falle et al., 1998; Powell et al., 1999; Balsara and Spicer, 1999; Tóth, 2000,; Dedner et al., 2002; Pen et al., 2003; Londrillo and Del Zanna, 2004; Ziegler, 2004; Crockett et al., 2005; Fromang et al., 2006; Mignone et al., 2007; Cunningham et al., 2007), pri- marily because such methods are very good for shock capturing, and flows involving contact and/or rotational discontinuities. In addition, single-step Eulerian methods that solve the equations of motion in the conservative form are better suited for use with static or adaptive mesh refinement (SMR or AMR, respectively). In a recent series of papers (Gardiner and Stone, 2005, 2008, here- after GS05 and GS08, respectively), we have described the extension of the directionally unsplit corner transport upwind (CTU) integra- tor of Colella (1990) to MHD using the constrained transport (CT) method of Evans and Hawley (1988) to enforce the divergence-free constraint on the magnetic field. We refer to this combination as the CTU + CT algorithm. Several novel ingredients of the MHD CTU + CT algorithm were shown to be important. In particular, spatial recon- struction schemes that contain a directionally-split time advance (such as the piecewise parabolic method, PPM, of Colella and Wood- ward, 1984, hereafter CW), must include additional source terms in multidimensional MHD compared to the one-dimensional (1D) case. Moreover, the transverse flux gradients used in CTU to correct the interface states in multidimensions also require similar terms. The nature and form of these terms is documented in detail in Sec- tions 3.1 and 4.1.2 of GS05 for two dimensions (2D), and Sections 3.1 and 4.2 in GS08 for three dimensions (3D). In addition, it was shown that simple arithmetic averaging is not sufficient to compute the edge-centered electric fields required by CT from the face-cen- tered fluxes returned by the Riemann solver. Instead, simple recon- struction algorithms were developed and tested (GS05). More recently, we have presented a detailed description of the implementation of the CTU + CT algorithm in Athena, a new code for astrophysical MHD (Stone et al., in press, hereafter S08). The results of a comprehensive hydrodynamic and MHD test suite were used to demonstrate the fidelity of the CTU + CT method in 1D, 2D, and 3D. During the course of developing the CTU + CT algorithm in Athe- na, we have experimented with other dimensionally unsplit integra- tors for MHD. We have found that a simple predictor–corrector scheme (see the appendix in Falle, 1991) similar to the MUSCL–Han- cock scheme described by van Leer (2006; see also Toro 1999) can form the basis of a robust Godunov method for MHD when com- bined with CT. We refer to this combination as the VL + CT algorithm. 1384-1076/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.newast.2008.06.003 * Corresponding author. Fax: +1 609 258 8226. E-mail address: jmstone@princeton.edu (J.M. Stone). 1 Present address: 3915 Rayado Pl NW, Albuquerque, NM 87114, United States. New Astronomy 14 (2009) 139–148 Contents lists available at ScienceDirect New Astronomy journal homepage: www.elsevier.com/locate/newast