Table Of Content

An Unsplit Godunov Method for Ideal MHD via Constrained Transport Thomas A. Gardinera & James M. Stonea,b aDepartment of Astrophysical Sciences 5 Princeton University 0 0 Princeton, NJ 08544 2 bProgram in Applied and Computational Mathematics n Princeton University a J Princeton, NJ 08544 5 2 1 Abstract v 7 We describe a single step, second-order accurate Godunov scheme for ideal MHD 5 5 based on combining the piecewise parabolic method (PPM) for performing spatial 1 reconstruction, the corner transport upwind (CTU) method of Colella for multidi- 0 mensional integration, and the constrained transport(CT) algorithm for preserving 5 0 the divergence-free constraint on the magnetic field. We adopt the most compact / h form of CT, which requires the field be represented by area-averages at cell faces. p We demonstrate that the fluxes of the area-averaged field used by CT can be made - o consistentwiththefluxesofthevolume-averaged fieldreturnedbyaRiemannsolver r if they obey certain simple relationships. We use these relationships to derive new t s algorithms for constructing the CT fluxes at grid cell corners which reduce exactly a : to the equivalent one-dimensional solver for plane-parallel, grid-aligned flow. We v i show that the PPM reconstruction algorithm must include multidimensional terms X for MHD, and we describe a number of important extensions that must be made r a to CTU in order for it to be used for MHD with CT. We present the results of a variety of test problems to demonstrate the method is accurate and robust. Key words: PACS: 1 Introduction In recent years a variety of numerical algorithms for multidimensional magnetohydrodynamics (MHD) based on Godunov’s method have been developed [3,9,10,16,23,25,28]. There are two important extensions to the basic hydro- dynamical algorithm that are required for MHD. The first is an extension Preprint submitted to Elsevier Science 2 February 2008 of the Riemann solver used to compute the fluxes of each conserved quan- tity to MHD; the second is a method by which the divergence free constraint B = 0 is imposed upon the numerically evolved magnetic field. Of the two, ∇· the latter has emerged as the more difficult to achieve. The fact that it is important to ensure the numerically evolved field satisfies the divergence free constraint was first noted by Brackbill & Barnes [4]. They pointed out that the Lorentz force is not orthogonal to B if B = 0, and ∇· 6 that this could lead to incorrect dynamics in a well defined test problem. More recently Tóth [32] has shown that in some circumstances it is possible to get the wrong jump conditions across MHD shocks if the constraint is not satisfied (this is also evident in the method of [16]). Currently, there are three methods by which the divergence free constraint is applied in Godunov schemes. The first is to use a Hodge projection to clean the magnetic field of any divergence after each time step (e.g. [2,9,32]). The second is to extend the system of conservation laws with an evolutionary equation for the divergence designed to minimize the accumulation of error in anyonelocation.Examplesofthisapproachincludetheeight-wavescheme[26] and the GLM-MHD scheme [12]. Finally, the third is to design the difference equations for the magnetic field to explicitly conserve magnetic flux, and so preserve the divergence free constraint. The latter method termed constrained transport (CT) by [15] has proved successful in other MHD algorithms [6,31] and is the method adopted here. The most compact CT difference formulae are built upon area-averaged magnetic field components located at the faces of a grid cell, rather than volume- averaged field components located at grid cell centers (CT algorithms built upon cell-centered fields have been developed in [32], however they require averaging over a stencil which is larger than that used to compute the fluxes). The need for a staggered grid is often thought of as a disadvantage of CT. In fact,however, itreflectsoneofthemost attractivepropertiesofCT:thefunda- mental conserved property of the magnetic field in MHD is the magnetic flux (which is an area- rather than volume-average) and by design, CT conserves the magnetic flux (and therefore it preserves the divergence free character of the magnetic field) in an integral sense, over the smallest discretization scale, the grid cell size. We would like to combine CT with a finite volume shock capturing method. It may at first seem inconsistent to build a numerical algorithm which mixes a finite volume approach (which conserves integrals of volume-averaged values) with CT (which conserves integrals of area-averaged magnetic fluxes at grid cell faces). In this paper we show that, provided the fluxes of the volume- and area-averaged fields obey certain simple relationships, the finite-volume and CT approaches can be made consistent (see also [23]). Most importantly, the 2 relationships we derive allow us to develop algorithms for constructing fluxes of the face-centered field (located at grid cell corners in 2D) from the fluxes of the volume-averaged field computed by a Riemann solver (located at grid cellfaces).Wedemonstratethatmultidimensionalalgorithmsdevelopedinthis way reduce exactly to the equivalent one-dimensional solver for plane-parallel, grid-aligned flow. This is one difference between the methods developed here and previous implementations of CT in Godunov schemes [3,11,28]. We combine the CT algorithms developed in this paper with the piecewise parabolic method (PPM) [8] using Roe’s linearization as the MHD Riemann solver [5]. An essential ingredient of PPM is a spatial reconstruction step to compute time-advanced estimates of the conserved variables at grid faces. In this paper we show that for MHD, this reconstruction step must include multidimensional terms in the induction equation (used to reconstruct the transverse components of the field). We argue that dimensionally split MHD algorithms cannot preserve the divergence free constraint between each one- dimensional update, therefore we adopt the unsplit corner transport upwind (CTU) algorithm of Colella [7] to develop a multidimensional algorithm. How- ever, there are a number of important extensions that must be made to CTU to make it suitable for MHD with CT. These include using a CT update for the magnetic field during the predict step, and inclusion of multidimensional terms along with the transverse flux gradients used to predict multidimensional fluxes. We describe these extensions in detail. Theresulting two-dimensionalMHDPPMalgorithmusesasinglestepupdate, is second order accurate, and is fully conservative. Hence it is ideally suited for use on a statically or adaptively refined mesh. Moreover, the two-dimensional algorithm reduces exactly to the base one-dimensional algorithm for planar, grid-aligned flows. Although this paper will only describe the combination of CT with CTU, there is no reason why the CT algorithm described here could not be combined with other unsplit methods. To simplify the discussion we confine ourselves to a two-dimensional algorithm in this paper. The paper is organized as follows. In 2 we review the finite volume and § CT methods in order to demonstrate the relationships between the area- and volume-averaged magnetic fields and their fluxes. In 3.3 we use these re- § lationships to derive algorithms for constructing the fluxes of face-centered area-averaged fields needed by CT from the fluxes of cell-centered volume- averaged fields returned by a Riemann solver. We present a single step first- order Godunov method and use it to test different algorithms for constructing the fluxes used in CT. This method forms the first half of the second-order CTU + CT integration algorithm developed in section 4. In 5 we present § § a variety of tests which demonstrate the linear and nonlinear behavior of the scheme. Finally, in 6 we conclude. § 3 2 Constrained Transport in Finite Volume Schemes The equations of ideal magnetohydrodynamics (MHD) can be written in conservative form as ∂ρ + (ρv)=0 (1) ∂t ∇· ∂ρv ∗ + (ρvv BB)+ P =0 (2) ∂t ∇· − ∇ ∂B + (vB Bv)=0 (3) ∂t ∇· − ∂E ∗ + ((E +P )v B(B v))=0 (4) ∂t ∇· − · where ρ is the mass density, ρv is the momentum density, B is the magnetic field, and E is the total energy density. The total pressure P∗ P +(B B)/2 ≡ · where P is the gas pressure, and the total energy density E is related to the internal energy density ǫ via E ǫ+ρ(v v)/2+(B B)/2 . (5) ≡ · · Throughout this paper we will assume an ideal gas equation of state for which P = (γ 1)ǫ, where γ is the ratio of specific heats. Unless otherwise stated, − we take γ = 5/3. None of the main results described in this paper depend directly upon the equation of state. Note also that we have chosen a system of units in which the magnetic permeability µ = 1. In addition to the evolutionary conservation laws equations (1) through (4), themagneticfieldmust alsoobeythedivergence freeconstraint,i.e. B = 0. ∇· It is of paramount importance that the numerically evolved field satisfy this constraint at all times, otherwise, for example, the system of equations for the conservative variables is inconsistent with the same system written in terms of the primitive variables, i.e. (ρ, v, B, P). The algorithm described in this paper is built upon finite volume (FV) methods, in which the conserved variables are averaged over grid cell volumes. On the other hand, the CT method is built upon area-averaging of the magnetic field, leading to difference equations for the magnetic flux throughthe surfaces of grid cells. In the following subsections, we briefly review the FV and CT schemes in order to arrive at consistent relationships between the volume- and area-averaged magnetic field and their associated numerical fluxes. 4 2.1 Finite Volume Method Consider a regular, two dimensional cartesian grid with grid cell (i,j) centered at (x ,y ) and of size (δx,δy). The conservative system of equations for ideal i j MHD can be written in vector form as ∂q + f = 0 (6) ∂t ∇· where ρ   ρv x     ρv   y      ρvz  q   (7)   ≡  B   x       B   y       Bz       E        is the vector of conserved variables and ρv ρv x y     ρv2 +P∗ B2 ρv v B B  x − x   y x − y x   ρv v B B   ρv2 +P∗ B2   x y − x y   y − y           ρvxvz BxBz   ρvyvz ByBz  f  −  ,f  −  ,(8) x   y   ≡  0  ≡  v B B v     y x y x     −       v B B v   0   x y − x y             vxBz Bxvz   vyBz Byvz   −   −      (E +P∗)v B (B v) (E +P∗)v B (B v)  x x   y y   − ·   − ·      are the flux vectors. Integrating over the volume of grid cell (i,j) and the time interval δt = tn+1 tn and applying Gauss’s theorem we obtain − δt δt qn+1 = qn + Fn+1/2 Fn+1/2 + Fn+1/2 Fn+1/2 (9) i,j i,j δx x,i−1/2,j − x,i+1/2,j δy y,i,j−1/2 − y,i,j+1/2 (cid:16) (cid:17) (cid:16) (cid:17) 5 the integral form of the evolution equation. The conserved quantities yi+δy/2xi+δx/2 1 qn q(x,y,tn)dxdy , (10) i,j ≡ δxδy yi−Zδy/2xi−Zδx/2 are averaged over the grid cell volume and the fluxes tn+1yi+δy/2 1 n+1/2 F f (x δx/2,y,t)dydt (11) x,i±1/2,j ≡ δyδt x i ± tZn yi−Zδy/2 tn+1xi+δx/2 1 n+1/2 F f (x,y δy/2,t)dxdt (12) y,i,j±1/2 ≡ δxδt y i ± tZn xi−Zδx/2 are averaged over the surface area of a grid cell face and the time interval δt. Typically one approximates the flux integrals in equations (11) and (12) to some order of accuracy, while maintaining strict conservation by evolving the conserved quantities throughequation(9).Notewhen written inthis form,the components of the flux vectors are non-zero for the transverse components of the magnetic field only, meaning that directionally split updates of the volume averaged field based on these fluxes will not generally satisfy the divergence free constraint between directional sweeps. This suggests directionally split algorithms are inappropriate for MHD. 2.2 Constrained Transport Method In the CT method, the integral form of the induction equation is based on area rather then volume averages. Starting from the differential form of the induction equation, ∂B + = 0 (13) ∂t ∇×E where the electric field = v B in ideal MHD, one may integrate over the E − × bounding surface of a grid cell and use Stoke’s Theorem to obtain δt Bn+1 = Bn + n+1/2 n+1/2 (14) x,i±1/2,j x,i±1/2,j δy Ez,i±1/2,j−1/2 −Ez,i±1/2,j+1/2 (cid:16) (cid:17) δt Bn+1 = Bn n+1/2 n+1/2 (15) y,i,j±1/2 y,i,j±1/2 − δx Ez,i−1/2,j±1/2 −Ez,i+1/2,j±1/2 (cid:16) (cid:17) 6 as the integral form of the evolution equation. The magnetic field components yi+δy/2 1 Bn B (x δx/2,y,tn)dy (16) x,i±1/2,j ≡ δy x i ± yi−Zδy/2 xi+δx/2 1 Bn B (x,y δy/2,tn)dx (17) y,i,j±1/2 ≡ δx y i ± xi−Zδx/2 are averaged over the grid cell bounding faces and tn+1 1 n+1/2 (x δx/2,y δy/2,t)dt (18) Ez,i±1/2,j±1/2 ≡ δt Ez i ± i ± tZn is averaged over the time interval δt. Note the fundamental representation of the magnetic field is an area-average at cell faces. Although CT-like difference formulae are possible based on volume averaged fields at cell centers [32], they preserve a discretization of the divergence on a different (larger) stencil than used to compute the fluxes. Just as in finite volume methods, one typically approximates the electric field (flux)integralinequation(18)tosomeorderofaccuracyandappliesequations (14)and(15)to evolve themagnetic fieldcomponents in time.Nevertheless, in a manner exactly analogous to the finite volume method, conservation of magnetic flux is strictly enforced, implying that the net magnetic charge interior to a grid cell vanishes at time tn+1 if it did so at time tn. As such, the preser- vation of B = 0 (in an integral sense) in the CT method is as fundamental ∇· as, e.g., the conservation of mass in a finite volume method. Moreover, in the CT approach, the magnetic field components are averaged over the smallest dimensional volume necessary so as to transform the differential equation into its integral form. In this way one maintains the maximal “point-wise” infor- mation possible, thereby minimizing the dissipation inherent in the averaging process. 2.3 Consistency of the CT & FV Methods To build a numerical scheme based on the CT method for the magnetic flux, and a FV method for the remaining conserved quantities, it is very important that the surface and volume averaged magnetic field components (and their fluxes) be coupled in a consistent manner. One common approach [11,28,3] 7 (which we also follow) is to define the volume-averaged magnetic field components at cell centers as equal to the average of the area-averaged values at cell faces, i.e. 1 Bn = Bn +Bn (19) x,i,j 2 x,i−1/2,j x,i+1/2,j (cid:16) (cid:17) 1 Bn = Bn +Bn (20) y,i,j 2 y,i,j−1/2 y,i,j+1/2 (cid:16) (cid:17) which is sufficient for second order accuracy. However, as shown below this choice implies a specific relationship between the numerical fluxes for the induction equation as integrated in the CT and FV formulations. The expressions which describe the coupling of the fluxes in the CT & FV methods are a direct result of requiring that equations (19) and (20) also hold at time tn+1. Subtracting equation (19) from an equivalent expression at time tn+1, and substituting equations (9) and (14) for the time differences of the volume- and area-averaged magnetic field respectively we find 1 n+1/2 n+1/2 n+1/2 n+1/2 eˆ F F = (21) Bx · y,i,j−1/2 − y,i,j+1/2 2 Ez,i−1/2,j−1/2 −Ez,i−1/2,j+1/2 (cid:16) (cid:17) (cid:16) (cid:17) 1 n+1/2 n+1/2 + . (22) 2 Ez,i+1/2,j−1/2 −Ez,i+1/2,j+1/2 (cid:16) (cid:17) where eˆ is a unit vector for the B component of the flux vector. It follows Bx x that 1 n+1/2 n+1/2 n+1/2 eˆ F = + (23) Bx · y,i,j±1/2 2 Ez,i−1/2,j±1/2 Ez,i+1/2,j±1/2 (cid:16) (cid:17) whichisconsistentwiththeobservationthatthefluxaveragesusedinthefinite volume method (equation 12) are spatial averages of the electric fields used in the constrained transport method (equation 18). Repeating this analysis for B we find y 1 n+1/2 n+1/2 n+1/2 eˆ F = − + (24) By · x,i±1/2,j 2 Ez,i±1/2,j−1/2 Ez,i±1/2,j+1/2 (cid:16) (cid:17) where eˆ is a unit vector for the B component of the flux vector. Function- By y ally, equations (23) and (24) imply that one must replace the Godunov fluxes for the volume averaged x- and y-components of the magnetic field with the average of the corner centered , regardless of the details of the CT algorithm z E used to compute the latter. Thus, equations (23) and (24) can be thought of as a corrector step that makes the predicted Godunov fluxes given by the Riemann solver consistent with the CT fluxes. Clearly, the CT algorithm used 8 to compute the corner centered will directly impact the accuracy and sta- z E bility of the underlying Godunov scheme; in 3.2 we address the problem of § constructing CT algorithms and the properties they should possess. 3 First Order CT Godunov Scheme Inthissectionwewillconstructandtestasinglestep,twodimensional,firstor- der integration algorithm based upon the piecewise parabolic method (PPM). The simplicity of this algorithm allows us to develop and test two of the most important elements of this paper: the calculation of the interface states in MHD and the systematic construction of CT algorithms. Moreover, the re- sulting algorithm (apart from the update step) is essentially the first half of the CTU + CT integration algorithm described in 4. § In section 3.1 we describe the calculation of the interface states in the PPM algorithm for the system of ideal MHD. This interface state calculation in- volves a characteristic evolution of a dimensionally split system. In this step we will find the appearance of truly multidimensional terms in the induction equation which are proportional to ∂B /∂x and ∂B /∂y. When these terms x y are included in the dimensionally split, linearized system of equations which are used to perform the characteristic evolution, they take the appearance of “source terms”. We present two gedanken experiments which show that these multidimensional terms are essential to accurately predicting the time evolution in the interface states. In 3.2 we address the question of the consistency of the CT algorithm with § the underlying finite volume integration algorithm. Repeating the arguments which lead to equations (23) and (24), we present an example of a CT algorithm which has insufficient dissipation and fails to reduce to the underlying integration algorithm for plane-parallel, grid-aligned flows. We proceed to describe a systematic approach to constructing a CT algorithm which can be applied with any approximate Riemann solver and reduces exactly to the underlying finite volume integration algorithm for plane-parallel, grid-aligned flows. In 3.3 we present a numerical study comparing three, surprisingly sim- § ple, CT algorithms each of which differs only in its dissipation for truly multidimensional problems. 3.1 Calculating the Interface States The algorithms presented in this paper are built upon the the piecewise parabolic method (PPM). For a thorough discussion of PPM or its linear 9 variant PLM and its implementation we refer the reader to the excellent de- scriptions in [8,29,24]. Roughly speaking the PPM algorithm can be broken down into threesteps: spatialreconstruction, characteristic evolution, andflux evaluation. The purpose of these first two steps is to calculate a one-sided esti- mate of the time averaged state at the left- or right-hand sides of a particular grid cell interface. With these interface states in hand, the interface flux may be calculated via either an exact or approximate Riemann solver. The calculation of the interface states in PPM is performed in primitive variables, and is a one-dimensional algorithm. However, in the two dimensional (x,y) system of equations for ideal MHD there appear terms proportional to ∂B /∂x and ∂B /∂y which are not present in the truly one-dimensional sys- x y tem. In primitive variables, these terms only appear in the induction equation, which in component form is ∂B ∂ x + (v B B v ) = 0 (25) y x y x ∂t ∂y − ∂B ∂ y + (v B B v ) = 0 (26) x y x y ∂t ∂x − ∂B ∂ ∂ z + (v B B v )+ (v B B v ) = 0 . (27) x z x z y z y z ∂t ∂x − ∂y − Using the magnetic charge constraint ( B = 0) these terms can be elimi- ∇· nated from equation (27), giving ∂B ∂ ∂v ∂ ∂v z z z + (v B ) B + (v B ) B = 0 . (28) x z x y z y ∂t ∂x − ∂x ∂y − ∂y However, no such simplification can be made to equations (25) and (26). It is natural to ask just how important these terms are and what role they play in the evolution of the magnetic field. A few gedanken experiments quickly show that they are absolutely essential, and at times are the dominant term in the equation. Perhaps the most trivial example of a case in which these terms play an important role is for stationary solutions. As a concrete example, consider a circularly polarized Alfvén wave oriented at some oblique angle to the grid. One may always choose a reference frame in which the Alfvén wave is stationary. Because the wave is oriented oblique to the grid, ∂B /∂x and ∂B /∂y x y are non-zero throughout the domain (except at extrema). In addition, for a standing Alfvén wave, the velocities are of the order of the Alfvén speed. In this case, the term v (∂B /∂y) in equation (25) is not a small term and in x y fact, it must exactly balance the remaining term ∂/∂y(v B ) B (∂v /∂y) y x y x − 10