A Novel High-Order, Entropy Stable, 3D AMR MHD Solver with Guaranteed Positive Pressure
Dominik Derigs, Andrew R. Winters, Gregor J. Gassner, Stefanie Walch
AA Novel High-Order, Entropy Stable, 3D AMR MHD Solver withGuaranteed Positive Pressure
Dominik Derigs a, ∗ , Andrew R. Winters b , Gregor J. Gassner b , Stefanie Walch a a I. Physikalisches Institut, Universit¨at zu K¨oln, Z¨ulpicher Straße 77, 50937 K¨oln b Mathematisches Institut, Universit¨at zu K¨oln, Weyertal 86-90, 50931 K¨oln
Abstract
We describe a high-order numerical magnetohydrodynamics (MHD) solver built upon a novel non-linearentropy stable numerical flux function that supports eight travelling wave solutions. By construction thesolver conserves mass, momentum, and energy and is entropy stable. The method is designed to treat thedivergence-free constraint on the magnetic field in a similar fashion to a hyperbolic divergence cleaningtechnique. The solver described herein is especially well-suited for flows involving strong discontinuities.Furthermore, we present a new formulation to guarantee positivity of the pressure. We present the underlyingtheory and implementation of the new solver into the multi-physics, multi-scale adaptive mesh refinement(AMR) simulation code
FLASH ( http://flash.uchicago.edu ). The accuracy, robustness and computationalefficiency is demonstrated with a number of tests, including comparisons to available MHD implementationsin FLASH . Keywords: magnetohydrodynamics,
FLASH , entropy stable, finite volume schemes, pressure positivity
1. Introduction
Modelling complex non-linear astrophysical phenomena is a central task in the field of astrophysics wherelaboratory experiments are very difficult if not entirely impossible. Examples of interesting phenomenainclude the study of stellar evolution, like the star formation process and supernovae explosions, pre-stellaraccretion discs and many more. Using simulations allows us to study the internals of complex systems thatcannot been seen in experiments and observations.In astrophysics, a flow involving magnetised gas is typically ionized, compressible, and often supersonic.Since the interstellar gas has essentially infinite conductivity [1], we treat the flow by solving the idealmagnetohydrodynamics (MHD) equations. From the hyperbolic nature of the ideal MHD equations, it isknown that discontinuous solutions may develop even from smooth initial data. Obtaining stable numericalresults for the variety of physical flow regimes is extremely challenging, particularly for the natural requirementthat the numerical scheme must be both accurate and robust. In this paper, we present a novel three-dimensional high-order, conservative, quasi-multifluid, entropy stable, eight wave MHD solver developedfor the numerical modelling of MHD flows. It is equally well suited for one, two, or three-dimensionalhydrodynamics (HD) and MHD simulations.The core of the novel MHD solver is the use of entropy stable flux functions developed in [2]. Entropystable algorithms have the benefit that, by construction, the numerical method is nearly isentropic in smoothregions and entropy is guaranteed to be increasing near discontinuities. Thus, the numerics precisely followthe physics of the second law of thermodynamics. Another advantage of entropy stable approximations isthat one can limit the amount of dissipation added to the numerical approximation to guarantee entropy ∗ Corresponding author
Email address: [email protected] (Dominik Derigs)
Preprint submitted to Elsevier May 13, 2016 a r X i v : . [ a s t r o - ph . GA ] M a y tability. The development and investigation of entropy stable algorithms for the ideal MHD equations hasbeen considered by several authors [2, 3, 4, 5].The entropy stable formulation also addresses the issue of divergence cleaning for approximate solutionsof the ideal MHD equations. The proof of entropy stability in [2] required an additional source term thatacts analogously to a hyperbolic divergence cleaning technique [6]. That is, errors introduced into thedivergence-free condition are advected away with the fluid velocity.The scheme handles another major robustness issue in numerical approximations of state-of-the-art high-order MHD solvers – the possible appearance of negative pressures. Negative pressures are a numerical artifactarising due to the problem of finite numerical precision. This phenomenon has been reported frequently inthe literature [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. In current codes, negative pressuresare avoided by adding artificially high amounts of dissipation or by introducing non-conservative low pressurelimits. Negative pressures can arise due the fact that the internal energy is obtained by subtracting thekinetic and magnetic energies from the conserved total energy. In many situations, such as high Mach numberor low plasma β flows ( β ∝ p/ (cid:107) B (cid:107) ), the internal energy can be several orders of magnitude smaller theneither the kinetic or magnetic energies. Thus, discretisation errors in the total energy could be significantenough to result in negative pressures. The inevitable consequence is the failure of the numerical scheme.We describe how the novel solver uses the entropy as an auxiliary equation to eliminate this issue andderive a novel expression for the pressure which completely avoids the subtraction problem. The new pressurepositivity guaranteeing formulation is not tied to any specific numerical flux function. It remains general andit is straightforward to retrofit into any existing HD/MHD schemes if the underlying numerical approximationis constructed in a way that satisfies certain criteria on the entropy (see Sec. 3.6).The new solver achieves high-order accuracy in space and time while remaining attractive from acomputational point of view. The numerical scheme is extended to high-order in space with spatialreconstruction techniques. In particular, we use a third order spatial approximation with the newly developedreconstruction technique of Schmidtmann et al. [24]. High-order accuracy in time is obtained using the familyof strong stability preserving (SSP) Runge-Kutta methods developed by Gottlieb et al. [25].We provide here details of the novel solver as well as its implementation into the multi-scale multi-physicssimulation code FLASH [26, 27].
FLASH is publicly available and has a wide international user base. Theremainder of this paper is organized as follows: Sec. 2 provides the necessary background information todiscuss the novel numerical solver. In Sec. 3 we describe, in detail, the new solver. The most important aspectsof which are the entropy stable numerical fluxes and the new pressure positivity guaranteeing formulation.Sec. 4 presents a variety of numerical results that demonstrate the utility of the new solver. We compare ourresults to already available MHD implementations in
FLASH where applicable. Sec. 5, presents our concludingremarks.
2. Governing Equations and Discretisation
We first provide the necessary background to discuss the novel MHD solver. This includes a briefdescription of the ideal MHD equations, the concept of entropy conservation and stability, and an outline ofthe finite volume scheme used for the spatial discretisation.
The ideal MHD model assumes that a fluid is a good electric conductor and neglects non-ideal effects, e.g. viscosity or resistivity. It is governed by a system of conservation laws ∂∂t (cid:37)(cid:37) u E B + ∇ · (cid:37) u (cid:37) ( u ⊗ u ) + (cid:0) p + (cid:107) B (cid:107) (cid:1) I − B ⊗ Bu (cid:0) E + p + (cid:107) B (cid:107) (cid:1) − B ( u · B ) B ⊗ u − u ⊗ B = 0 , (2.1) ∇ · B = 0 , (2.2)2here (cid:37) , (cid:37) u , and E are the mass, momenta, and total specific energy of the plasma system, p is the thermalpressure, I is the identity matrix, and B is the magnetic field, also referred to as magnetic flux density. Sinceour velocities are non-relativistic, Maxwell’s displacement current may be ignored in the Lorentz force term.We consider the non-dimensional form of the ideal MHD equations. Details concerning physical units can befound in Appendix D.Numerical methods for multidimensional ideal MHD must satisfy some discrete version of the divergence-free condition (2.2). There are several approaches to control the error in ∇ · B and in depth review of manymethods can be found in T´oth [17]. The thermal pressure is related to the conserved quantities through theideal gas law for problems in which relativistic, viscous, and resistive effects can be neglected: p = ( γ − (cid:18) E − (cid:37) (cid:107) u (cid:107) − (cid:107) B (cid:107) (cid:19) (2.3)with the ratio of specific heats γ > ∂∂t ( ∇ · B ) + ∇ · (cid:0) u ( ∇ · B ) (cid:1) = 0 , (2.4)is obtained. From (2.4) we see that the divergence of the magnetic field may be treated as an advectedscalar. Thereby, the robustness and accuracy of a numerical scheme can be significantly improved [28]. Thisimprovement is primarily because the advection of the generated errors prevents the accumulation at fixedlocations. The eigenmode which is advected with the flow in (2.4) is referred to as the divergence wave .We include the Janhunen source term [18] in the ideal MHD equations (2.1) which is proportional to ∇ · B . The use of a source term to control the error in the divergence free condition has known issues, suchas errors can build up at stagnation points as well as in periodic or closed domains. The only mechanismpresent to remove these divergence errors is the numerical dissipation of a scheme, but true hyperbolicdivergence cleaning methods can remove such limitations [6]. However, the Janhunen source term preservesthe conservation of mass, momentum, total energy and allows for the construction of an entropy stableapproximation [2]. We explicitly “clean” magnetic field divergence errors in a post-processing step, as will bedescribed later. The governing equations in conjunction with the Janhunen source term are now a system ofbalance laws ∂∂t (cid:37)(cid:37) u E B + ∇ · (cid:37) u (cid:37) ( u ⊗ u ) + (cid:0) p + (cid:107) B (cid:107) (cid:1) I − B ⊗ Bu (cid:0) E + p + (cid:107) B (cid:107) (cid:1) − B ( u · B ) B ⊗ u − u ⊗ B = − ( ∇ · B ) u . (2.5)Note that the expression “source term” is common in this context, even though the term actually involvesspatial derivatives.To simplify the discussion of the new solver we first consider the modified ideal MHD system (2.5) in onespatial dimension ∂∂t Q + ∂∂x F = Υ , (2.6)where Q = Q ( x, t ) is the vector of conservative variables, F ( Q ) the flux vector, and Υ ( Q ) is the vector3ource term Q = (cid:37)(cid:37)u(cid:37)v(cid:37)wEB B B , F = (cid:37) u(cid:37)u + p + (cid:107) B (cid:107) − B (cid:37) u v − B B (cid:37) u w − B B u (cid:0) E + p + (cid:107) B (cid:107) (cid:1) − B ( u · B )0 u B − v B u B − w B , Υ = − ∂B ∂x uvw . (2.7)In Sec. 3.1 we provide a detailed discussion of the multi-dimensional extension of the solver. This section serves as a brief introduction to entropy and numerical partial differential equations. Athorough review of this topic has been presented by Tadmor [29]. Work specifically related to entropy andthe ideal MHD equations can be found in [2, 4].It is well-known that solutions of balance laws like (2.6) may develop discontinuities in finite time, so weconsider solutions of the balance laws (2.6) in the weak sense. Unfortunately, the weak solution is not unique.Thus, we require an additional admissibility condition on the solution to guarantee that the numericalapproximation will converge to a weak solution that is consistent with the second law of thermodynamics.In the case of ideal MHD a suitable condition can be defined in terms of the physical entropy density, asdefined by Landau [30, p. 315], divided by the constant ( γ −
1) for convenience, i.e. S ( Q ) = (cid:37)sγ − s = ln (cid:0) p(cid:37) − γ (cid:1) , (2.8)where γ is the adiabatic index and s is the entropy per particle. The approximation obeys the second law ofthermodynamics and is based on an entropy condition for two regimes:1. For smooth solutions, one can design numerical methods to be entropy conservative if, discretely,the local changes of entropy are the same as predicted by the continuous entropy conservation law ∂∂t S + ∂∂x F = 0 , (2.9)where we define the corresponding entropy flux F ( Q ) = uS = (cid:37)usγ − . (2.10)2. For discontinuous solutions, the approximation is said to be entropy stable if the entropy alwayspossesses the correct sign and the numerical scheme produces more entropy than an entropy conservativescheme and satisfies the entropy inequality ∂∂t S + ∂∂x F ≥ . (2.11)From the second law of thermodynamics, kinetic as well as magnetic energy can be transformed irreversiblyinto heat (internal energy). If additional dissipation is not included in an entropy conservative method,spurious oscillations will develop near discontinuities as energy is re-distributed at the smallest resolvablescale [31]. A numerical scheme requires a diffusion operator to match such a physical process.For the entropy stable solver discussed in this paper we use the provably entropy stable approximateRiemann solver derived in [2]. 4 .3. Finite Volume Scheme The finite volume method is a discretisation technique for partial differential equations especially usefulfor the approximation of systems of hyperbolic conservations laws. The finite volume method is designed toapproximate conservation laws in their integral form, e.g. , (cid:90) V Q t d x + (cid:90) ∂V F · ˆ n d S = 0 . (2.12)In one spatial dimension we divide the interval, V , into cells V i = (cid:2) x i − / , x i +1 / (cid:3) , (2.13)and the integral equation of a balance law with a source term becomesdd t (cid:90) x i +1 / x i − / Q d x + (cid:2) F ∗ (cid:0) x i +1 / (cid:1) − F ∗ (cid:0) x i − / (cid:1) (cid:3) = (cid:90) x i +1 / x i − / Υ d x. (2.14)A common approximation is to assume a constant solution within the cell [32, p. 436]: (cid:90) x i +1 / x i − / Q d x ≈ (cid:90) x i +1 / x i − / Q i d x = Q i ∆ x i . (2.15)Note that the finite volume solution is typically discontinuous at the boundaries of the cells. To resolve this,we introduce the idea of a “numerical flux”, F ∗ ( Q R , Q L ), often derived from the (approximate) solution of aRiemann problem. The function F ∗ takes the two states of the solution at an element interface and returnsa single flux value. For consistency, we require that F ∗ ( Q , Q ) = F , (2.16)that is, the numerical flux is equivalent to the physical flux if the states on each side of the interface areidentical.Next, we address the discretisation of the source term Υ in (2.14). There is a significant amount offreedom in the source term discretisation. The explicit discretisation of the source term is given in Sec. 3.4.We note that the discrete source term at each left ( i − /
2) and right ( i + 1 /
2) interface will contribute incell i . So, the semi-discrete finite volume method is( Q t ) i + 1∆ x i (cid:2) F ∗ i +1 / − F ∗ i − / (cid:3) = 12 (cid:0) Υ i − / + Υ i +1 / (cid:1) . (2.17)
3. Description of the Novel Entropy Stable MHD Solver
Here we describe the
FLASH implementation of the entropy stable ES solver in three spatial dimensions.In Sec. 3.1 we discuss the extension of the solver to three dimensions using dimensional splitting. Sec. 3.2presents a spatial reconstruction scheme used to achieve a high-order approximation. We describe the explicittime integration technique in Sec. 3.3. The entropy conservative and entropy stable numerical flux functionsare described in Sec. 3.4 and Sec. 3.5, respectively. The new strategy to numerically guarantee the positivityof the pressure is described in Sec. 3.6. The adaptive mesh refinement (AMR) functionality of FLASH andthe new implementation is found in Sec. 3.7. Next, in Sec. 3.8, a brief summary of a quasi-multifluidimplementation is provided. Sec. 3.9 describes how to couple gravity into the entropy stable solver. Thetreatment of the divergence-free condition in higher spatial dimensions is described in Sec. 3.10. Finally,Sec. 3.11 summarizes the MHD update procedure in
FLASH .5 .1. Multi-dimensionality We extend the one-dimensional set of MHD equations (2.6) to two or three spatial dimensions. In thecase of an underlying grid structure that is logically rectangular (like Cartesian grid geometries) a simpleand efficient way of extending the one-dimensional Riemann solver to higher spatial dimensions is to use dimensional splitting . The method of dimensional splitting has become popular in fluid dynamics as it allowsus to apply our knowledge about one-dimensional systems directly to multi-dimensional systems. Using thedimensional splitting method, one-dimensional problems along each coordinate direction are solved in turnto determine the fluxes across the faces of a finite volume cell. It has proven to be an inexpensive way ofextending one-dimensional high-resolution methods to higher dimensions [32, p. 103].We experience that in multi-physics simulations, commonly performed using FLASH , the MHD solveraccounts for less than 10% of the overall CPU time ( e.g. [33]). Thus, an MHD discretisation which allowslarge time steps is beneficial for the overall computational efficiency of the multi-physics framework. It iswell-known that dimensionally split schemes give larger time steps than comparable unsplit schemes wherethe dimensionality directly enters the CFL condition. Although the technique of dimensional splitting reducesthe accuracy of the solver to formally second-order, the overall increase in efficiency is often favourable forpractical applications.If the three-dimensional semi-discrete problem can be written in the form of( Q t ) i + A ( Q ) + B ( Q ) + C ( Q ) = 0 , (3.1)then the total update (3.1) can be split up into an x-sweep ( Q t ) i + A ( Q ) = 0 , (3.2)a y-sweep ( Q t ) i + B ( Q ) = 0 , (3.3)and a z-sweep ( Q t ) i + C ( Q ) = 0 , (3.4)where A ( Q ), B ( Q ), and C ( Q ) are operators for the vector of quantities Q in x , y , and z − directions,respectively. Each of the sweep operators is a compact notation to write the numerical flux and source termcontributions for a given spatial direction. For example, the operator A ( Q ) in three dimensions has the form A ( Q ) = 1∆ x i (cid:16) F ∗ i +1 / ,j,k − F ∗ i − / ,j,k (cid:17) − (cid:0) Υ i − / ,j,k + Υ i +1 / ,j,k (cid:1) . (3.5)Therefore, in each sweep direction, separate solutions of the Riemann problem and source term values arecomputed to update the quantities stored in Q n according to (2.17).To compute the sweeps in y - and z -directions, any direction dependent quantities, i.e. velocity andmagnetic field components, are rotated in order to solve them with the same algorithm that is used for the x -sweep. The finite volume method used by the
FLASH framework approximates the solution with quantities whichare constant within each cell. If one considers these values as point-wise approximations of the solutionlocated at each cell centre, this method computes the numerical interface fluxes at a distance of ∆ x/ not strictly rectangular since cells of different spatial sizes are allowed to coexist on the same grid
6n interface. Rather than using piecewise constant data, we use reconstructed quantities within each cell,( ˜ Q i ) L , R . Reconstruction functions, ( p i ) L , R , allow the computation of the approximated interface quantities( ˜ Q i ) L = Q i −
12 ( p i ) L , and ( ˜ Q i ) R = Q i + 12 ( p i ) R . (3.6)The reconstructed quantities (3.6) are then used to compute high-order accurate numerical fluxes in thefinite volume scheme (2.17), i.e. ˜ F i − / = F ∗ (cid:16)(cid:0) ˜ Q i − (cid:1) R , (cid:0) ˜ Q i (cid:1) L (cid:17) and ˜ F i +1 / = F ∗ (cid:16)(cid:0) ˜ Q i (cid:1) R , (cid:0) ˜ Q i +1 (cid:1) L (cid:17) . (3.7)The resulting high-order accurate semi-discrete approximation, reorganizing (2.17), is of the form( Q t ) i = 1∆ x i (cid:16) ˜ F i − / − ˜ F i +1 / (cid:17) + 12 (cid:0) Υ i − / + Υ i +1 / (cid:1) , (3.8)as illustrated in Figure 1. i − i i + 1cell index:source terms:fluxes: ˜ F i − / ˜ F i +1 / Υ i − / Υ i +1 / ( ˜ Q i − ) L Q i − ( ˜ Q i − ) R ( ˜ Q i ) L Q i ( ˜ Q i ) R ( ˜ Q i + ) L Q i + ( ˜ Q i + ) R Fig. 1.
Graphical representation of the quantities used in (3.7) and (3.8). Reconstructed quantities used for thecomputation of the numerical fluxes are highlighted in blue. The cell-centred quantities are printed in black.
For our reconstruction we use the third order accurate shock capturing limiting procedure for numericalsolutions of hyperbolic conservation laws recently described by Schmidtmann et al. [24]. Their scheme utilizesa local piecewise-parabolic reconstruction away from discontinuities (see Fig. 2) and reads( p i ) L = p ( Q i − , Q i , Q i +1 ) = + 23 Q i − − Q i − Q i +1 = 2 δ i − − δ i + , (3.9)( p i ) R = p ( Q i +1 , Q i , Q i − ) = − Q i − − Q i + 23 Q i +1 = 2 δ i + − δ i − , (3.10)with δ i − = Q i − Q i − and δ i + = Q i +1 − Q i . (3.11)However, such a reconstruction is known to cause oscillations in non-smooth solutions. This can be seen asa direct consequence of Godunov’s Theorem [34]. To avoid oscillations, we use the limiting procedure ofSchmidtmann et al. [24] to switch to a lower-order accurate reconstruction near large gradients, shocks anddiscontinuities. The solution of a system of hyperbolic conservation laws may not be smooth. In such cases inaccuratetime-integration schemes can suffer from poor performance such as an excessively small time step size due tothe presence of spurious oscillations as well as the progressive smearing, clipping or squaring of the numericalapproximation. To alleviate such performance issues, we consider a third order accurate explicit high-order7 − i − i i + 1cell index: ∆ x ̺ x b b b b Interface associated with cell i (˜ ̺ i − ) R (˜ ̺ i ) L Fig. 2.
Principle of our spatial reconstruction. This example shows the parabolic reconstruction of a specific densitypattern. The cell-centred quantities, (cid:37) i − , (cid:37) i − , (cid:37) i , and (cid:37) i +1 , are represented by blue dots. Our scheme uses a localthree-point stencil and is thereby computationally very efficient. strong-stability-preserving (SSP) low-storage Runge-Kutta time-integration scheme [25]. Such schemes arealso referred to in the literature as total variation diminishing (TVD) [35]. However, Gottlieb et al. [36]showed this moniker is misleading as their strong stability property holds in any norm and not only theTVD norm. We complete the discretisation of the reconstructed method (3.8) with the third order SSPRunge-Kutta scheme: Q (cid:48) = Q n + ∆ t · Q t ( Q n ) , (3.12) Q (cid:48)(cid:48) = 34 Q n + 14 (cid:16) Q (cid:48) + ∆ t · Q t ( Q (cid:48) ) (cid:17) , (3.13) Q n +1 = 13 Q n + 23 (cid:16) Q (cid:48)(cid:48) + ∆ t · Q t ( Q (cid:48)(cid:48) ) (cid:17) . (3.14)SSP Runge-Kutta schemes consist of convex combinations of explicit forward Euler integration. Thus, thefamily of methods are guaranteed to be stable under the same time step restriction [25]. We find that thethird order SSP Runge-Kutta time integration enables us to use larger time steps, which is favorable in ourmulti-physics framework.To select a stable time step for a computational run we use the CFL condition∆ t ≤ CFL · min (cid:20) ∆ xλ x max , ∆ yλ y max , ∆ zλ z max (cid:21) , (3.15)where λ d max is the speed of the largest wave at time step n travelling in d = { x, y, z } direction, CFL isthe user-definable CFL coefficient,
CFL ∈ (0 , λ max is known exactly, then the choice CFL = 1 . λ max is usually computed in some approximate way. Thus, a moreconservative choice for the CFL coefficient is typically used in practice ( e.g. CFL = 0 . For the entropy analysis of the ideal MHD equations the divergence-free condition is incorporated into thesystem of conservation laws as a source term [18, 38]. Both the Powell [28] and Janhunen [18] source termstreats the magnetic field as an advected scalar. However, the Janhunen source term remains conservative inthe momentum and total energy equations and restores the positivity of the Riemann problem as well asLorentz invariance [39].The discussion of the entropy conserving numerical flux function of [2] requires the introduction of somenotation. We introduce the jump (cid:74) · (cid:75) , the arithmetic mean ( · ) A as well as the logarithmic mean ( · ) ln of theleft/right states, denoted by ( · ) L and ( · ) R , respectively. These operators are defined as (cid:74) · (cid:75) = ( · ) R − ( · ) L , ( · ) A = ( · ) L + ( · ) R , and ( · ) ln = (cid:74) · (cid:75)(cid:74) ln( · ) (cid:75) . (3.16)8 numerically stable procedure to compute the logarithmic mean is described by Ismail and Roe [40, AppendixB]. For convenience we also introduce z = (cid:114) (cid:37)p , and z = √ (cid:37)p. (3.17) It was shown in [2] that the Janhunen source term can be used to design numerical schemes that guaranteethe discrete conservation of the entropy density for the ideal MHD equations. Guaranteeing this discreteconservation of the entropy density requires a particular discretisation of the Janhunen source term:12 (cid:0) Υ i − / + Υ i +1 / (cid:1) , (3.18)with Υ i − / = − (cid:74) B (cid:75) i − / ( uz ) A ( B ) A (∆ xz B ) A ( vz ) A ( B ) A (∆ xz B ) A ( wz ) A ( B ) A (∆ xz B ) A i − / , and Υ i +1 / = − (cid:74) B (cid:75) i +1 / ( uz ) A ( B ) A (∆ xz B ) A ( vz ) A ( B ) A (∆ xz B ) A ( wz ) A ( B ) A (∆ xz B ) A i +1 / . (3.19) The recently developed provably entropy conserving flux of Winters and Gassner [2] reads: F ∗ , ec = ˆ (cid:37) ˆ u ˆ p + ˆ (cid:37) ˆ u + (cid:16) ˚ B + ˚ B + ˚ B (cid:17) − ˚ B ˆ (cid:37) ˆ u ˆ v − (cid:92) B B ˆ (cid:37) ˆ u ˆ w − (cid:92) B B γγ − ˆ u ˆ p + ˆ (cid:37) ˆ u (cid:0) ˆ u + ˆ v + ˆ w (cid:1) + ˆ u (cid:16) ˆ B + ˆ B (cid:17) − ˆ B (cid:16) ˆ v ˆ B + ˆ w ˆ B (cid:17) u ˆ B − ˆ v ˆ B ˆ u ˆ B − ˆ w ˆ B , (3.20)with the averaged quantities and productsˆ (cid:37) = ( z ) A ( z ) ln , ˆ p = ( z ) A ( z ) A , ˆ p = γ + 12 γ ( z ) ln ( z ) ln + γ − γ ( z ) A ( z ) A , ˆ u = ( z u ) A ( z ) A , ˆ v = ( z v ) A ( z ) A , ˆ w = ( z w ) A ( z ) A , ˆ u = (cid:0) z u (cid:1) A ( z ) A , ˆ v = (cid:0) z v (cid:1) A ( z ) A , ˆ w = (cid:0) z w (cid:1) A ( z ) A , ˆ B = ( B ) A , ˆ B = ( B ) A , ˆ B = ( B ) A , (cid:92) B B = ( B B ) A , ˚ B = (cid:0) B (cid:1) A , ˚ B = (cid:0) B (cid:1) A , ˚ B = (cid:0) B (cid:1) A , (cid:92) B B = ( B B ) A . (3.21)9n the case of smooth solutions, the entropy conserving flux (3.20) conserves the entropy density of thesystem up to the precision of the scheme. In order for the numerical scheme to be applicable for possiblynon-smooth solutions we must extend the purely entropy conserving flux to become an entropy stable flux. Entropy conserving approximations suffer breakdown in the presence of discontinuities, which results inlarge oscillations in post-shock regions. Therefore, we require dissipation to be added to the approximationin an entropy consistent manner to guarantee discrete satisfaction of the entropy inequality (2.11). Thework [2] derived two provably entropy stable approximate Riemann solvers for the ideal MHD equations.In this work we present a new hybrid entropy stable approximation that continuously combines these twoentropy stable fluxes. This introduces explicit non-linearity to permit the calculation of sharp shock frontsand contact discontinuities.
To build an entropy stable approximation we use the entropy conservative approximation (3.20) as abaseline. In particular the work [2] presented two possible dissipation terms that can be added to the entropyconserving scheme:
ES-Roe : a matrix dissipation entropy stabilization. Similar to a Roe type method it selectively appliesdissipation to each of the travelling wave solutions, particularly close to shocks. ES-LLF : a scalar dissipation entropy stabilization. A simple, local Lax-Friedrichs type dissipation mechanism.Due to the simplicity of ES-LLF it cannot distinguish between the various waves present in the MHDflow and can, therefore, lead to a severe smearing of the approximation near discontinuities.Here we outline the construction of the
ES-Roe stabilization. The
ES-LLF stabilization follows almostimmediately. To build the matrix dissipation term we first select the dissipation matrix to be | (cid:98) A | . That isthe absolute value of the flux Jacobian for the ideal MHD 8-wave formulation: (cid:98) A = F Q + P = A + P , (3.22)where A is the flux Jacobian for the homogeneous ideal MHD equations and P is the Powell source term [28]written in matrix form, i.e. P ∂ Q ∂x = B B B u · B u v w ∂∂x (cid:37)(cid:37)u(cid:37)v(cid:37)w(cid:37)eB B B = ∂B ∂x B B B u · B uvw . (3.23)The design of the entropy stable matrix dissipation term requires the specific form of the flux Jacobian(3.22) because it must be possible to relate the eigenvectors of (3.22) to the entropy Jacobian matrix [41].This relationship is referred to as creating entropy scaled eigenvectors, e.g. [41, 42]. To ensure that thisentropy scaling exists, the system of PDEs must be symmetrizable. It is known that the Powell source termrestores the symmetric property to the ideal MHD system [38, 42], whereas the Janhunen source term doesnot restore symmetry. However, both source terms allow to contract the MHD equation to the entropyevolution equation and hence both source terms can be used to construct entropy conserving (or stable)discretisations. We choose the Janhunen source term to construct the entropy conservative discretisation,as this gives us conservation of mass, momentum and energy unlike a method based on the Powell sourceterm. Thus, the consistent symmetric part of the flux is based on the Janhunen source term. As long as the10dditional stabilisation term is guaranteed to dissipate entropy, the scheme is entropy stable. Hence, forthe design of the stabilisation term only, we are considering the flux Jacobian that incorporates the Powellsource term, as this guarantees that the entropy scaled eigenvectors exist for the ideal MHD system, which isnecessary in order to get the Roe type dissipation term.The matrix type stabilization term requires the eigenstructure of the dissipation matrix (3.22) (cid:98) A = (cid:98) RD (cid:98) R − . (3.24)The matrix (cid:98) A supports eight propagating plane-wave solutions: • two fast magnetoacoustic waves ( ± f ), • two slow magnetoacoustic waves ( ± s ), • two Alfv´en waves ( ± a ), • an entropy wave ( E ), • a divergence wave ( D ).It is known that a naively scaled set of right eigenvectors will exhibit several forms of degeneracy that arecarefully described by Roe and Balsara [43]. We follow the same rescaling procedure of Roe and Balsarato improve the numerical behaviour of the fast/slow magnetoacoustic eigenvectors. The matrix of righteigenvectors is (cid:98) R = [ (cid:98) r +f | (cid:98) r +a | (cid:98) r +s | (cid:98) r E | (cid:98) r D | (cid:98) r − s | (cid:98) r − a | (cid:98) r − f ] , (3.25)with the eigenvectors (cid:98) r , and corresponding eigenvalues λ [2, 42, 43]Entropy and Divergence Waves: λ E,D = u (cid:98) r E = uvw (cid:107) u (cid:107) , (cid:98) r D = B , (3.26)Alfv´en Waves: λ ± a = u ± b (cid:98) r ± a = ± (cid:37) β ∓ (cid:37) β ∓ (cid:37) ( β w − β v )0 − (cid:37)β (cid:37)β , (3.27)11agnetoacoustic Waves: λ ± f, ± s = u ± c f,s (cid:98) r ± f = α f (cid:37)α f (cid:37) ( u ± c f ) (cid:37) ( α f v ∓ α s c s β σ ( b )) (cid:37) ( α f w ∓ α s c s β σ ( b ))Ψ ± f α s aβ √ (cid:37)α s aβ √ (cid:37) , (cid:98) r ± s = α s (cid:37)α s (cid:37) ( u ± c s ) (cid:37) ( α s v ± α f c f β σ ( b )) (cid:37) ( α s w ± α f c f β σ ( b ))Ψ ± s − α f aβ √ (cid:37) − α f aβ √ (cid:37) , (3.28)where we introduced several convenience variablesΨ ± s = α s (cid:37) (cid:107) u (cid:107) − aα f (cid:37)b ⊥ + α s (cid:37)a γ − ± α s c s (cid:37)u ± α f c f (cid:37)σ ( b )( vβ + wβ ) , Ψ ± f = α f (cid:37) (cid:107) u (cid:107) aα s (cid:37)b ⊥ + α f (cid:37)a γ − ± α f c f (cid:37)u ∓ α s c s (cid:37)σ ( b )( vβ + wβ ) ,c = b , c , s = 12 (cid:18) ( a + b ) ± (cid:113) ( a + b ) − a b (cid:19) , a = γ p(cid:37) ,b = b + b + b , b ⊥ = b + b , b = B √ (cid:37) , β , , = b , , b ⊥ ,α = a − c c − c , α = c − a c − c , σ ( ω ) = (cid:40) +1 if ω ≥ , − . (3.29)In (3.29), for the wave speed computation c , s , the plus sign corresponds to the fast magnetoacoustic speed, c , and the minus sign corresponds to the slow magnetoacoustic speed, c .The entropy stable dissipation term is built from three components: • Entropy scaled matrix of right eigenvectors: ˚ R = (cid:98) R √ T , where T is the diagonal scaling matrix T = diag (cid:18) (cid:37)γ , p (cid:37) , (cid:37)γ , (cid:37) ( γ − γ , p(cid:37) , (cid:37)γ , p (cid:37) , (cid:37)γ (cid:19) . (3.30)For the complete motivation and details on the entropy scaling of eigenvectors see Barth [42]. • Diagonal matrix of eigenvalues: For
ES-Roe each wave component is weighted with a different eigenvalue,whereas
ES-LLF weights all wave components identically | D ES − Roe | = diag (cid:0) | λ + f | , | λ + a | , | λ + s | , | λ E | , | λ D | , | λ − s | , | λ − a | , | λ − f | (cid:1) , (3.31a) | D ES − LLF | = diag (cid:0) λ max , λ max , λ max , λ max , λ max , λ max , λ max , λ max (cid:1) . (3.31b)The maximum eigenvalue λ max is given by λ max = max (cid:0) | λ + f | , | λ + a | , | λ + s | , | λ E | , | λ D | , | λ − s | , | λ − a | , | λ − f | (cid:1) . (3.32) • Jump in the entropy vector: (cid:74) v (cid:75) The term (cid:74) v (cid:75) is the jump between left and right states of the entropy vector, which is defined as avector field whose components are partial derivatives of the entropy density (2.8) with respect to the12uid quantities Q , v = d S d Q = − (cid:20) γ − sγ − − (cid:37) (cid:107) u (cid:107) p , (cid:37)up , (cid:37)vp , (cid:37)wp , − (cid:37)p , (cid:37)B p , (cid:37)B p , (cid:37)B p (cid:21) (cid:124) , (3.33)with the physical entropy, s , defined in (2.8).The general form of the ES-Roe , and
ES-LLF numerical flux functions is F ES = F ∗ , ec + 12 ˚ R | D | ˚ R T (cid:74) v (cid:75) , (3.34)where F ∗ , ec is the entropy conserving numerical flux (3.20). Note that the only difference between the ES-Roe and
ES-LLF stabilizations is in the selection of the diagonal matrix of eigenvalues D . The matrix ofright eigenvectors and the eigenvalues are discretely computed from the previously defined average quantities(3.21) to ensure consistency, in the presence of vanishing magnetic fields, with the entropy stable Euler solverof Ismail and Roe [40]. Chandrashekar [44] points out that most, if not all, schemes which resolve grid aligned stationary contactdiscontinuities exactly suffer from the carbuncle effect as the profiles around shocks can exhibit spuriousoscillations. This can also be true in our case as our flux function guarantees only the correct sign of theentropy but not necessarily the correct amount of entropy production. However, a flux function must generate enough entropy across a shock to guarantee monotonicity [45]. The usual fix for this problem, i.e. increasingthe amount of induced dissipation, causes poor resolution of features of boundary layers or near shocks.Another possibility is to switch the numerical scheme to a more dissipative one only near shocks and use ahigh resolution Riemann solver in smooth parts of the flow [46]. It is straightforward to implement such anidea in the current context because entropy stable schemes have the freedom to select the eigenvalues thatessentially control the amount of dissipation.The local Lax-Friedrichs type scalar dissipation,
ES-LLF , effectively suppresses the carbuncle phenomenon.However, we want to use the more accurate Roe type matrix dissipation,
ES-Roe , in regions without largepressure jumps to be able to track smooth parts of the solutions with more accuracy. To achieve this goal weconstruct a hybrid entropy stabilization scheme, called
ES-Hybrid , that blends the
ES-Roe and the
ES-LLF scheme continuously. In the hybrid scheme, a new diagonal matrix of eigenvalues is defined as | D ES − Hybrid (Ξ) | = (1 − Ξ) | D ES − Roe | + Ξ | D ES − LLF | , (3.35)with the limitslim Ξ → | D ES − Hybrid (Ξ) | = | D ES − Roe | , lim Ξ → | D ES − Hybrid (Ξ) | = | D ES − LLF | . (3.36)As was done in [44], we define the parameter Ξ ∈ [0 ,
1] using a simple local pressure jump indicatorΞ = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p L − p R p L + p R (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) / . (3.37)From the design of the pressure indicator (3.37), the scheme uses mainly the less dissipative ES-Roe schemefor smooth parts of the flow (but also near e.g. contact discontinuities), while the more dissipative
ES-LLF entropy-stabilization is used near strong shocks. 13 .6. Pressure Positivity Guaranteeing Formulation
We next address the issue that negative pressures may be introduced by a numerical scheme. This hasbeen described in previous publications, e.g. [19, 20, 22]. We present here a general and physically motivatedsolution to the specific numerical issue of negative pressures.In a classical higher order Godunov method, the internal energy and thereby the thermal pressure, p th , isobtained by subtracting the kinetic and magnetic energies from the conserved total energy (2.3). In manysituations, as in high-Mach number or low plasma-beta flows ( β ∝ p/ (cid:107) B (cid:107) ), the internal energy can beseveral orders of magnitude smaller then the kinetic or magnetic energies. Thus, discretisation errors in thetotal energy might be significant enough to result in negative pressures leading to a failure of the numericalscheme. This problem is often addressed by enforcing low pressure limits. However, it is questionable ifthe simulation can then still give a physically meaningful solution. Therefore, it is important to design aconservative pressure positivity guaranteeing scheme that is physically convincing. Ryu et al. [22] state that in regions where the gas is very cold compared to the bulk kinetic energy, theflow cannot be treated using the total energy approach as the errors in calculating the total energy can belarger than the internal energy itself. In order to overcome this difficulty, they solve an entropy conservationequation and extract the pressure directly wherever the internal energy is much less than the kinetic energy, i.e. E th /E kin (cid:28)
1. They present several different criteria used to select whether to compute the pressurefrom the total energy or from their entropy formulation.Balsara and Spicer [19] extend the idea of Ryu et al. to MHD flows and present two strategies to preventnegative pressures. Their “strategy 1” is to use the pressure computed from the entropy density only inthose cells where the thermal pressure could potentially become negative. In all other cells, they use thethermal pressure given by (2.3). Their “strategy 2” uses the pressure computed from the entropy everywhereexcept in regions near strong magnetosonic shocks or a flow configuration that may develop such shocks.They justify the validity of their approach by noting that their work deals with magnetospheric problems.There are no shocks present in a magnetosphere, but there still remains a positivity problem.Li [15] extends the ideas of Balsara and Spicer to an implementation of two new equations: The entropyequation used in [19] and an internal energy equation. Similar to the previously mentioned works, he pointsout that these equations should not be used close to or within regions that contain shocks. His shock detectionscheme sets a floor value for the internal energy in cells near a shock.Balsara [20] presents a general strategy to address problems where the positivity of density and pressureis uncertain. His work addresses the problem that positivity can be lost when using reconstruction schemes.He presents a self-adjusting strategy for enforcing the positivity of the pressure. However, we realize that thepositivity problem of pressure can also be encountered with schemes that are first-order accurate in spaceand as such do not utilize a reconstruction scheme. This is due to the problem of subtracting large numberswith accordingly large discretisation errors as stated above.The work of Ersoy et al. [10] relies on improving the resolution in problematic regions. They utilize adiscrete measure of the entropy in order to stabilize their computation. They use the local entropy production as a mesh refinement criterion on their computational grid.
The current solver is built from an entropic perspective. Thus, at any time in the computation, we cancompute the entropy density as well as the discrete entropy flux. With these tools we determine a value forthe pressure that is guaranteed to be positive. From the computed entropy density for each cell within thecomputational domain, we use (2.8) to derive a new expression for the pressure in the cells: p s = exp (cid:20) γ − (cid:37) S + γ ln( (cid:37) ) (cid:21) . (3.38)From (3.38), it is immediately clear that this “entropy pressure” will always be positive as exp( x ) > , x ∈ R .Hence, our solver fulfils the desired property of being pressure positivity guaranteeing under all circumstances.14e note that our scheme can be used in a similar way as described by Balsara and Spicer [19]. However,our scheme includes a proper treatment of the entropy at shocks. It is applicable in all regions of the flowand not only in sufficiently smooth regions.The current scheme is: Use the “normal” scheme as long as the internal energy is large enough after theupdate with the criterion E int /E tot > smalleint with the user-definable parameter smalleint that defaultsto 0 .
01. If the internal energy is smaller than the criterion, we switch to the entropy pressure formulationwithout violating the conservation of total energy.
It is straightforward for a given semi-discrete finite volume method to compute the entropy update ofthe method. This is because we know how to convert the equations into entropy space. We contract thesemi-discrete equations (possibly including some reconstruction technique) (3.8) with the entropy vector(3.33) to obtain v T ( Q t ) i = v T (cid:18) x i (cid:16) ˜ F i − / − ˜ F i +1 / (cid:17) + 12 (cid:0) Υ i − / + Υ i +1 / (cid:1)(cid:19) . (3.39)From the chain rule and definition of the entropy vector (3.33) we know that v T ( Q t ) i = ( S t ) i . (3.40)So, we have an expression for the time evolution of the entropy density( S t ) i = v T (cid:18) x i (cid:16) ˜ F i − / − ˜ F i +1 / (cid:17) + 12 (cid:0) Υ i − / + Υ i +1 / (cid:1)(cid:19) , (3.41)where, as was shown in [2], the entropy stable fluxes provide a discrete approximation to the spatial derivativeof the entropy flux, i.e. ∂∂x ( uS ) ≈ v T (cid:18) x i (cid:16) ˜ F i − / − ˜ F i +1 / (cid:17) + 12 (cid:0) Υ i − / + Υ i +1 / (cid:1)(cid:19) . (3.42)Thus, we see that (3.41) is a consistent, discrete update for the entropy. This new discrete equation for theentropy density can be added to the MHD system (2.5) and evolved in time with the other fluid quantities.We reiterate that, by construction, the entropy stable approximation will guarantee that entropy is consistentwith the second law of thermodynamics everywhere. Thus, the proposed positive pressure guaranteeingmethod is valid in any region of the flow. The implemented procedure is:1. We update the entropy density S n +1 with the same time integration scheme used to obtain Q n +1 .2. If the updated energies violate the criterion E n +1int /E n +1tot > smalleint , we use (3.38) to get p n +1 s .3. Finally, we recompute the updated internal energy from p n +1 s to make the scheme consistent. FLASH incorporates an adaptive mesh refinement (AMR) strategy using the
PARAMESH library [47], throughwhich the grid is organized in a block structured, oct-tree adaptive grid. The presented entropy stable solveris fully incorporated into
FLASH ’s AMR functionality to optimize computational costs. For completeness webriefly discuss the underlying AMR structure, and parallelization in
FLASH . With AMR, the local spatialresolution can be dynamically controlled. This allows the maximization of the computational efficiency ofthe overall simulation as higher resolution is placed only where it is needed.Parallelization is achieved by dividing the computational domain into several blocks (sub-domains). Ablock contains a number of computational cells (
NXB , NYB , and
NZB in the x , y , and z − direction, respectively).The default block contains NX|Y|ZB = 8. Each block is surrounded by a fixed number of guard cells in15 yz Block
Block boundarySingle cell
Fig. 3.
Two blocks with different levels of refinement(i.e. mesh resolution). A single x sweep is highlightedin red. The guard cells are not shown in this figure. Fig. 4.
An adaptive grid 2D simulation with differentlevels of refinement. The interior cells of one of theblocks are highlighted in yellow. The according guardcells are shown in grey. Guard cells that extend intoblocks having a different grid size are interpolated. each spatial direction, providing the block with information from its neighbouring blocks. The completecomputational domain consists of a number of blocks (most likely with different physical sizes). The three-dimensional structure of the blocks is sketched in Fig. 3, while a simple two-dimensional slice through anadaptive grid is shown in Fig. 4. Three rules apply in the creation of refined blocks:1. A refined block must be one-half of the size of the parent block in each spatial dimension ( e.g. eachrefinement of a block gives 8 additional blocks in three-dimensional computations).2. Refined blocks must fit within the parent block and are not allowed to overlap into other blocks (theyhave to be aligned).3. Blocks sharing a common border are not allowed to differ in more than one level of refinement.Each block contains all information about local and neighbouring cells, making the blocks with the surroundingguard cells self-contained computational domains which allows efficient parallel computation using the MessagePassing Interface (MPI) framework. We configure AMR in such a way that adaptive refinement is allowedafter each two consecutive time steps ( nrefs = 2).
The ability to track the exact composition of a gas is of central importance in astrophysical simulationsas they include detailed chemical networks to treat heating, cooling, as well as molecule formation anddestruction to mimic the behaviour of the interstellar medium (ISM) [33, 48, 49].In order to track the different chemical species in the gas, advection equations of the form ∂X (cid:96) (cid:37)∂t + ∇ · ( X (cid:96) (cid:37) u ) = 0 , (3.43)are solved, where X (cid:96) is the fractional abundance of the (cid:96) th species ( H , H + , H , He , etc.) with the unityconstraint (cid:80) (cid:96) X (cid:96) = 1. For each species the flow of the quantity is calculated by multiplying the fractionalabundances of the species in the cells with the total density fluxes. Our scheme was originally devised for aperfect gas with a constant ratio of specific heats, γ . We generalize our scheme for a multi-species fluid withvariable γ by adopting a mean value of γ at the cell interfaces as suggested by Murawski [50].16e implement the multi-species advection in a similar way as recommended by Plewa & M¨uller [51](known as Consistent Multi-fluid Advection (CMA) method). That is, we ensure that the species fluxes areconsistent during the advection. Note that many existing schemes instead normalize the abundances afterthe advection step. However, as Glover et al. [52] pointed out, this procedure lacks any formal justificationand can lead to large systematic errors in the abundance of the least abundant chemical species.In addition to the multifluid approach using different chemical species, we implement mass tracer fields(also called mass scalars or tracerfields). These are field variables which are advected similar to species massfractions by an equation of the form ∂ψ(cid:37)∂t + ∇ · ( ψ(cid:37) u ) = 0 , (3.44)where ψ is the mass fraction, and ψ(cid:37) is the partial density of the traced mass.Our implementation of the mass tracer fields into the MHD solver allows the use of any number of suchfields. Thus, the mass tracer fields are a flexible tool for tracing different mass quantities according toindividual requirements. For example, a mass tracer field could be used to follow the distribution of metalsin the interstellar gas with virtually no additional computational costs. The inclusion of gravity in the ideal MHD equations (2.1) introduces a force into the right-hand side ofthe momentum equations ∂∂t (cid:37) u + ∇ · (cid:20) (cid:37) ( u ⊗ u ) + (cid:18) p + 12 (cid:107) B (cid:107) (cid:19) I − B ⊗ B (cid:21) = − (cid:37) ∇ φ, (3.45)where the gravitational potential φ satisfies Poisson’s equation ∇ φ = 2 πG(cid:37), (3.46)with the universal gravitational constant G that is an empirical physical constant involved in the calculationsof gravitational forces between two bodies. FLASH provides several algorithms for solving the Poisson equation (3.46). We tested our implementationwith a Barns & Hut tree-based algorithm implemented by R. W¨unsch (
Poisson/BHTree ) [53] and the Fouriertransform-based multigrid algorithm Poisson solver (
Poisson/Multigrid ) [54].
Within the MHD equations (2.1), the divergence free condition of the magnetic field (2.2) is not modelleddirectly. While this constraint is physically fulfilled at any time, we will see that care must be taken to fulfillthis constraint numerically.The extension to higher spatial dimensions, as described in Sec. 3.1, has been performed in a straightforwardmanner by relying on the Cartesian grid structure. In one dimension the divergence-free condition impliesthat the longitudinal component of the magnetic field is constant over time. However, this conclusion doesnot generalize to two and three spatial dimensions.Instead, due to discretisation errors, a non-zero divergence of the magnetic field occurs over time whichinevitably leads to the issue that the conservation of the magnetic flux cannot be maintained. Thesediscretisation errors effectively generate numerical magnetic monopoles that grow exponentially during thecomputation and cause the magnetic field to no longer be solenoidal. From the equations of ideal MHD (2.1)it is clear that these monopoles cause an artificial force parallel to B .In Sec. 2.1, we noted that errors in the divergence-free condition are dealt with by treating the divergenceof the magnetic field as an additional fluid quantity to prevent accumulation of errors when the divergence isnon-zero in the computational domain. The eigenmode which is advected with the flow in (2.4) is referredto as the divergence wave . This procedure might be understood as a form of divergence cleaning for the17agnetic field. However, numerical experiments show that this approach might not be sufficient to maintainadequate divergence-free magnetic fields throughout simulations.Concerning divergence cleaning, there are different techniques available (see e.g. [55]). One particularexample is the elliptic projection, based on the Helmholtz decomposition, originally developed by Chorin[56]. Brackbill and Barnes [57] and Marder [58] developed a projection method in the context of the MHDequations. This method effectively suppresses the growth of unphysical magnetic monopoles locally as shownby Murawski [50] and T´oth [17]. The projection method has successfully been applied by e.g. Zachary etal. [21], Balsara [59], and more recently by Crockett et al. [60]. We implement the projection method fordivergence cleaning as a separate post-processing step and note that our original scheme remains unchanged.The general downside of this scheme may be the high computational costs caused by the projectionapproach. Our implementation is based of the realization that although the divergence problem is of ellipticalcharacter, its influence is only local. Accordingly, we design our implementation of the projection method ina way that is purely local and thereby computationally favourable.To enforce the divergence-free constraint we subtract the portion of the magnetic field that violates ∇ · B = 0. Suppose that the divergence of the magnetic field in the computation is non-zero. An easy fix tothis problem is the addition of a correction field ˜ B such that ∇ · (cid:0) B + ˜ B (cid:1) = 0 . (3.47)To guarantee physical consistency of the magnetic field correction, it is clear that ˜ B must not cause anyadditional current j ∝ ∇ × ˜ B = 0 . (3.48)Hence, we conclude that ˜ B must have the form˜ B = ∇ φ, (3.49)where φ is a scalar potential. Combining (3.47) and (3.49) we obtain∆ φ = ∇ φ = −∇ · B . (3.50)Note that the Laplace operator, ∆, has the physical interpretation for non-equilibrium diffusion as the extentto which a point represents a source or sink of some concentration. The resulting scalar potential can thenbe used to evaluate ˜ B according to (3.49) and, thus, clean the magnetic field B : B (cid:12)(cid:12)(cid:12) ∇· B =0 = B + ˜ B . (3.51)By the projection of the cell-centred magnetic fields onto the space of divergence-free magnetic fields, oneis left with fields at the next time step which are divergence-free to very good approximation. We note thatprojecting the magnetic field in the way described is consistent with the underlying cell-centred scheme.An alternative approach to divergence cleaning that should be mentioned is the constraint transportmethod developed by Evans and Hawley [61] or Balsara and Spicer [62] (reviewed in [17]). In this approach,the divergence-free constraint is satisfied by placing the staggered magnetic field at cell faces instead ofcell centres. On such a grid, the MHD equations can be approximated such that they preserve numericalsolenoidality of the magnetic field by construction. Note that Balsara and Kim [63] found advantages for thestaggered-mesh in their comparison between divergence- cleaning and divergence- free methods for stringenttest cases. However, the staggered grid approach has the downside of being much more expensive in terms ofstorage. In addition, it is not clear if provably stable schemes can be constructed for staggered-meshes [64].The precise implementation of our divergence cleaning approach is described in further detail in AppendixB. 18 .11. MHD Update Procedure On logically Cartesian grid geometries, it is straightforward to solve multi-dimensional problems as setsof one-dimensional problems by using the dimensional split approach. This approach is the principle of thenew ES solver. The MHD equations are solved as one-dimensional problems along each coordinate directionin turn ( x , y , and z − sweeps) in order to determine the fluxes through the finite volume cell surfaces. Block boundariesGuard Cells
Fig. 5.
Principle of the one-dimensional solution update with four guard cells.
Each one-dimensional sweep (like the one highlighted in Fig. 5) works as follows:1. First, the quantities are converted from primitive to conservative form ( e.g. velocity to momentum).2. For y and z − sweeps the solution array is rotated such that we solve this sweep as if it would be an x − sweep. This allow us to use our one-dimensional algorithms without modification.3. For each cell within the array, the reconstructed quantities ( ˜ Q i − ) R and ( ˜ Q i ) L are computed using thespatial reconstruction scheme (see Sec. 3.2) at the current time.4. Then, the entropy stable numerical fluxes as well as the source terms are computed using the algorithmsdescribed in Secs. 3.4 and 3.5.The default behaviour is to use the ES-Hybrid flux due to its flexibility. However, the user can easily changewhich flux function the computation uses with a single switch. Depending on the settings, either the entropyconserving fluxes, F ∗ , ec , the matrix dissipation entropy stable fluxes, F ES − Roe , the scalar dissipation entropystable fluxes, F ES − LLF , or the hybrid entropy stable fluxes, F ES − Hybrid are used.5. After this preparation, the solution array is updated using the time integration scheme described above.6. The updated internal energies are then derived from the updated total energy.We update the total energy as it is a conserved quantity. From the updated total energy, we derive theinternal energy by subtracting the magnetic and kinetic energies as suggested in [65]: E n +1int = E n +1tot − (cid:16) E n +1mag + E n +1kin (cid:17) with E n mag = 12 (cid:13)(cid:13) B n (cid:13)(cid:13) , and E n kin = 12 (cid:37) n (cid:13)(cid:13) u n (cid:13)(cid:13) (3.52)If the computed internal energy fails the criterion E int /E tot > smalleint , then the total energy update isdone with the pressure computed from the entropy density described in Sec. 3.6.7. Finally, the variables are converted to primitive form as other FLASH modules expect primitive variables.8. In higher dimensions, the divergence cleaning procedure, described in Sec. 3.10, is used to diffuse awayerrors in the divergence-free condition as a post-processing step.19 . Numerical Results
We demonstrate the utility, robustness, and accuracy of the new solver by computing the solution toseveral well-known HD and MHD test problems. The version of
FLASH on hand is 4.3 as of 18 th July, 2015.We consider six numerical test cases to test the performance of our new solver and compare to resultsobtained using already available MHD solver implementations for
FLASH . A test that extends the well-knownShu-Osher test to MHD is presented in Sec. 4.1 which is used to test the ES scheme’s artificial dissipation in1D. The propagation of smooth Alfv´en waves is studied in Sec. 4.2. We forgo the presentation of furtherone-dimensional results as we felt multi-dimensional results were more valuable to the present discussion.The application of the entropy stable MHD solver to the shock tube problems of Brio and Wu [66], Ryuand Jones [67], and Torrilhon [68] can be found in Winters and Gassner [2]. In Sec. 4.3 we further explorethe accuracy of the method in multiple spatial dimensions by considering the Orszag-Tang vortex problem.The MHD rotor problem originally proposed by Balsara [62] is investigated in Sec. 4.4. The MHD Rotoris also used in Sec. 4.5 to compare CPU timing and memory consumption of the new ES solver and theother schemes. Sec. 4.6 provides an example of using gravity with the new solver by considering the Jeansinstability. We note that the Jeans instability is a pure HD configuration and demonstrates that the newMHD solver remains applicable to flows with vanishing magnetic fields. Finally, we test the conservation ofthe available MHD schemes using the involving MHD blast wave test discussed in Sec. 4.7. All tests, exceptthe Jeans instability test, are performed using dimensionless units. Each test is run with CFL = 0 . The test proposed by Shu and Osher [69] is commonly used to test a scheme’s ability to resolve small-scalefluid features in the presence of a supersonic shock. A sinusoidal density/entropy perturbation is addeddownstream of a Mach 3 shock wave. The interaction of the shock wave with the perturbations gives riseto complex fluid features as the shock amplifies the initial oscillations. This test is an excellent testbed tomeasure the numerical (artificial) viscosity of a scheme. Additionally, the presence of a supersonic shock isused to demonstrate the robustness and stability of a scheme [70]. We consider a complex MHD version ofthe Shu-Osher problem recently developed by Susanto [71]. We present the initial conditions for this test inTable 1. The left and right boundaries are taken sufficiently far from the initial discontinuity such that theydo not influence the solution. This test has no analytic solution, so we compute a reference solution on ahighly refined grid using the
MHD 8Wave solver for comparison. x ≤ x x > x (cid:37) . . x ) u . v . w p . B B . B { x min , x max } = {− , } Initial shock location x = − t max = 0 . γ = 5 / Table 1
Initial conditions and runtime parameters: MHD Shu-Osher test (1D)
Fig. 6 shows the density at t = 0 . ES-Hybrid solver captures the small-scale flow features much better than the other schemes available in
FLASH . Also, no stability or overshoot problems are visible in the solutions.According to previous investigations, e.g. [73], a scheme is considered “acceptable” for capturing supersonicturbulence if the dynamics can be captured well with 400 cells. However, the entropy stable scheme is alsoable to resolve the dynamics of the flow with a much lower spatial resolution (the result used in Fig. 6 is 208cells in total). 20 − − − . . . . . . . . . . . . . . . . . . . . . ES-Hybrid USM MHD 8Wave Bou5 reference
Fig. 6.
Density of the MHD Shu-Osher problem at t = 0 .
7. These plots be compared to Fig. 1 of [72] or Fig. 3.9 of[71]. We used an adaptive grid resolution of up to 256 cells. The reference solution is computed on a uniform grid of4096 cells.
The smooth Alfv´en wave test [17] is used to compare the accuracy of MHD schemes for smooth flows.The initial circularly polarized Alfv´en wave propagates across a periodic domain. For the 2D test, weincline the smooth Alfv´en wave by an angle of α = 45 ◦ relative to the x -axis. The Alfv´en wave speed is | v A | = B (cid:107) / √ (cid:37) = 1 and thus, the wave is expected to return to its initial state at each time t ∈ N . Thistest is run to a final time t max = 5 . .
6. We introduce additional notation for aperpendicular coordinate x (cid:107) = x · cos( α ) + y · sin( α ) as well as the parallel, B (cid:107) = 1 .
0, and perpendicular, B ⊥ = 0 . πx (cid:107) ), magnetic fields. The field in z -direction is given by B z = 0 . πx (cid:107) ). The initialconditions listed in Table 2 ensure that the magnetic pressure is constant.The ability to propagate Alfv´en waves over long times and distances is crucial for e.g. MHD turbulencesimulations. If the Alfv´en waves are damped strongly because of inherent numerical dissipation in a scheme,the code will fail to capture the resulting turbulence behaviour correctly as MHD turbulence is mainlysustained by Alfv´en waves [74].Density (cid:37) p u B ⊥ · (cid:0) − sin( α ) , cos( α ) , (cid:1) (cid:124) + B z · (cid:0) , , (cid:1) (cid:124) Mag. field B B (cid:107) · (cos( α ) , sin( α ) , (cid:124) + u Domain size { x, y } min = { . , . }{ x, y } max = { / cos( α ) , / sin( α ) } Boundary conditions periodicSimulation end time t max = 5 . γ = 5 / Table 2
Initial conditions and runtime parameters: Smooth Alfv´en wave test (1D, 2D)
In the one dimensional smooth Alfv´en wave test we check the spatial high resolution properties ofour scheme. For sufficiently smooth fields, i.e. in cases where discontinuous features are absent, the usedreconstruction technique is designed to achieve third order accuracy (see Sec. 3.2).21o test the accuracy of our scheme, we run several simulations with varying resolutions and compute the L and L errors for the quantity B ⊥ = B y cos( α ) − B x sin( α ) as described in Appendix E. The obtainederrors are plotted as a function of the number of grid points in logarithmic scale in Fig. 7 and are listed inTable 3. As can be seen, third order accuracy is achieved, already at very low resolutions. N − − − − L / L e rr o r s rd order 2 nd order nd o r d e r EC ES-Hybrid MHD 8Wave USM Bouchut5
Fig. 7. L (solid lines) and L (dashed lines) errors measured with the smooth Alfv´en wave test in 1D. We omit thelines for the ES-Roe and
ES-LLF schemes as they are visually indistinguishable from
ES-Hybrid (cf. Table 3). . . . . . . x − . − . . . . B ⊥ ES-RoeES-LLF ES-HybridEC MHD 8WaveBouchut5 USM exact . . . . . . x .
20 0 .
25 0 . . . . . Fig. 8.
Smooth Alfv´en wave test after five crossing times. For the left plot, we used a fixed grid of 8 cells. For theright plot we use a grid of 16 cells. The exact solution shows the initial configuration at the given resolution.
Fig. 8 shows B ⊥ vs. x ⊥ at time t = 5 for the one dimensional Alfv´en wave test. As we know that thesolution is smooth, we disable the entropy stabilisation term described in Sec. 2.2 and obtain an entropyconserving EC scheme ( ). The EC solution shows very little dissipation. Note that the EC scheme is onlyapplicable to smooth solutions and should not be used for arbitrary flows. We observe that the different ES schemes ( , , and ) resolve the Alfv´en wave with the least dissipation of all tested MHDsolvers (except the entropy conserving scheme) while their results are virtually identical. The MHD 8Wave
USM implementation ( ) [70, 75, 76] are considerablymore diffusive. They show second order convergence. Finally, we note that the
Bouchut5 implementation( ) [12] has the highest measured amount of dissipation for this one dimensional test case.
Fig. 9 shows B ⊥ vs. x ⊥ at time t = 5 for the two dimensional Alfv´en wave test. The EC solution ( )shows again very little dissipation. As before, the ES schemes resolve the Alfv´en wave with the least dissipationof all tested MHD solvers while the ES-Roe scheme ( ) is least dissipative and the
ES-LLF scheme ( )is slightly more diffusive. As expected for smooth problems, the
ES-Hybrid scheme ( ) gives resultsthat are identical to those computed using the
ES-Roe scheme. The
MHD 8Wave implementation ( ) givessimilar results compared to the ES solver but is slightly more diffusive. The unsplit USM implementation( ) shows a higher dissipation compared to the ES or MHD 8Wave implementations and its zero-crossingpoints are clearly shifted at the lower resolution run. We find that the
Bouchut5 implementation ( ) hasthe highest amount of dissipation for this smooth test case. Note that the dissipation of the Alfv´en waves issignificantly reduced in higher dimensions if multidimensional Riemann solvers with sub-structure are used,as was shown by Balsara [74]. . . . . . . . . x ⊥ − . − . . . . B ⊥ ES-RoeES-LLF ES-HybridEC MHD 8WaveBouchut5 USM exact . . . . . . . . x ⊥ .
15 0 . . . . . Fig. 9.
Smooth Alfv´en wave test after five crossing times. These plots be compared to Fig. 8 of [17] or Fig. 4 of [4].For the left plot, we used a fixed grid of 16 ×
16 cells. For the right plot we use a grid of 32 ×
32 cells. The exact solution shows the initial configuration at the given resolution. We use
CFL = 0 . USM solver solution.
In Fig. 10 we plot the evolution of the conserved quantities as well as the individual energies. Looking,for example, at the magnetic energy, E mag , it can be seen that the EC scheme introduces the least amount ofdissipation. It is followed by our entropy stable schemes ES-Roe / ES-Hybrid and
ES-LLF . The
MHD 8Wave and the
USM implementations show higher dissipation while the
Bouchut5 implementation shows the highestamount of dissipation. If one would only look at the internal energy, one might conclude that the
MHD 8Wave solver introduces even less dissipation that the ES schemes. However, one has to be cautious with such aconclusion because both MHD 8Wave as well as
Bouchut5 fail to preserve total energy conservation.We list the computed L and L errors for the quantity B ⊥ in Table 4. They support our conclusionsgiven above, e.g. ES-Roe and
ES-Hybrid give identical solutions for smooth problems. Due to dimensionalsplitting, the obtained results are only second-order accurate in space.23
C ES-Roe ES-Hybrid ES-LLF USM MHD 8Wave Bouchut5 N = L error 7 . · − < . · − = 1 . · − = 1 . · − < . · − ≈ . · − < . · − L error 7 . · − < . · − = 1 . · − = 1 . · − < . · − ≈ . · − < . · − N = L error 1 . · − < . · − = 1 . · − = 1 . · − < . · − < . · − < . · − L error 1 . · − < . · − = 1 . · − = 1 . · − < . · − < . · − < . · − EOC ( L ) 2 .
82 3 .
45 3 .
45 3 .
45 2 .
54 2 .
54 2 . L ) 2 .
78 3 .
41 3 .
41 3 .
41 2 .
52 2 .
52 2 . N = L error 1 . · − < . · − = 1 . · − = 1 . · − < . · − < . · − < . · − L error 1 . · − < . · − = 1 . · − = 1 . · − < . · − < . · − < . · − EOC ( L ) 2 .
95 3 .
22 3 .
22 3 .
22 2 .
68 2 .
65 1 . L ) 2 .
95 3 .
21 3 .
21 3 .
21 2 .
67 2 .
64 1 . N = L error 1 . · − ≈ . · − = 1 . · − = 1 . · − (cid:28) . · − < . · − < . · − L error 1 . · − ≈ . · − = 1 . · − = 1 . · − (cid:28) . · − < . · − < . · − EOC ( L ) 2 .
98 3 .
06 3 .
06 3 .
06 2 .
39 2 .
34 2 . L ) 2 .
99 3 .
06 3 .
06 3 .
06 2 .
39 2 .
34 1 . N = L error 2 . · − ≈ . · − = 2 . · − = 2 . · − (cid:28) . · − < . · − < . · − L error 2 . · − ≈ . · − = 2 . · − = 2 . · − (cid:28) . · − < . · − < . · − EOC ( L ) 2 .
99 3 .
01 3 .
01 3 .
01 2 .
14 2 .
11 2 . L ) 2 .
99 3 .
01 3 .
01 3 .
01 2 .
14 2 .
11 2 . Table 3
Computed errors and experimental order of convergence (EOC) for B after five oscillation of the Alfv´en wave in onedimension ( t = 5 . EC ES-Roe ES-Hybrid ES-LLF MHD 8Wave USM Bouchut5 N = L error 8 . · − < . · − = 1 . · − < . · − < . · − < . · − < . · − L error 9 . · − < . · − = 1 . · − < . · − < . · − < . · − < . · − N = L error 1 . · − < . · − = 2 . · − < . · − < . · − < . · − < . · − L error 2 . · − < . · − = 2 . · − < . · − < . · − < . · − < . · − EOC ( L ) 2 .
19 2 .
50 2 .
51 2 .
55 2 .
75 1 .
83 1 . L ) 2 .
17 2 .
48 2 .
49 2 .
56 2 .
77 1 .
88 1 . N = L error 4 . · − < . · − = 4 . · − < . · − < . · − < . · − < . · − L error 4 . · − < . · − = 5 . · − < . · − < . · − < . · − < . · − EOC ( L ) 2 .
06 2 .
32 2 .
32 2 .
44 2 .
41 1 .
46 1 . L ) 2 .
06 2 .
30 2 .
30 2 .
41 2 .
42 1 .
30 1 . N = L error 1 . · − ≈ . · − = 1 . · − ≈ . · − < . · − < . · − < . · − L error 1 . · − ≈ . · − = 1 . · − ≈ . · − < . · − < . · − < . · − EOC ( L ) 2 .
01 2 .
10 2 .
10 2 .
18 2 .
12 1 .
36 1 . L ) 2 .
01 2 .
10 2 .
10 2 .
16 2 .
09 1 .
36 1 . Table 4
Computed errors and experimental order of convergence (EOC) for B ⊥ after five oscillation of the Alfv´en wave in twodimensions ( t = 5 . . − . − . − . − . − . − . . . E t o t a l × − + 1 . − − − − − − − ( % u ) t o t a l × − ES-RoeES-LLF ES-HybridEC MHD 8WaveBouchut5 USM % t o t a l × − + 2 0 . . . . . . E m ag t . . . . . . . . . . E k i n × − t . . . . . E i n t Fig. 10.
Evolution of the conserved quantities as well as the individual energies in the smooth Alfv´en wave test overfive crossing times (16 ×
16 cells). .3. Orszag-Tang MHD Vortex (2D) The Orszag-Tang vortex problem [77] is a two-dimensional, spatially periodic problem well suited forstudies of MHD turbulence. Thus, it has become a classical test for numerical MHD schemes. It includesdissipation of kinetic and magnetic energy, magnetic reconnection, formation of high-density jets, dynamicalignment and the emergence and manifestation of small-scale structures. The Orszag-Tang MHD vortexproblem starts from non-random, smooth initial data. As the flow evolves it gradually becomes increasinglycomplex, forming intermediate shocks. Thus, this problem demonstrates the transition from initially smoothdata to compressible, supersonic MHD turbulence. The initial data is chosen such that the root mean squarevalues of the velocity and the magnetic fields as well as the initial Mach number are all one. The averageplasma beta is β = .Density (cid:37) . p . /γ Velocity u ( − sin(2 πy ) , sin(2 πx ) , . (cid:124) Mag. field B γ ( − sin(2 πy ) , sin(4 πx ) , . (cid:124) Domain size { x, y } min = { . , . }{ x, y } max = { . , . } Boundary conditions all: periodicAdaptive refinement on density, magnetic fieldSimulation end time t max = 0 . γ = 5 / Table 5
Initial conditions and runtime parameters: Orszag-Tang MHD vortex test
Additionally, we compute the experimental convergence order for the available MHD schemes after 10 % ofthe total runtime. At this time, the solution is already very complex but still smooth. As there is no analyticsolution available, we compare to a high resolution simulation obtained using our entropy-conserving ( EC )scheme on an uniform grid of 2048 × L and L errors as well as the experimentalorder of convergence and list them in Table 6. They coincide with our results presented in the precedingsection. As in Sec. 4.2, we find the results to be at most second-order accurate in space due to dimensionalsplitting. EC ES-Hybrid Bouchut5 USM MHD 8Wave N = L error 6 . · − < . · − < . · − < . · − < . · − L error 8 . · − < . · − < . · − < . · − ≈ . · − N = L error 1 . · − ≈ . · − < . · − ≈ . · − ≈ . · − L error 2 . · − ≈ . · − < . · − ≈ . · − ≈ . · − EOC ( L ) 2 .
04 2 .
10 1 .
80 1 .
95 1 . L ) 2 .
07 2 .
08 1 .
81 1 .
93 1 . N = L error 4 . · − ≈ . · − < . · − < . · − ≈ . · − L error 5 . · − ≈ . · − < . · − < . · − > . · − EOC ( L ) 2 .
00 2 .
03 2 .
02 1 .
91 1 . L ) 1 .
97 2 .
00 2 .
01 1 .
81 1 . Table 6
Computed errors and experimental order of convergence (EOC) for p mag = (cid:107) B (cid:107) before the onset of discontinuitiesin the Orszag-Tang MHD vortex test ( t = 0 . Fig. 11 displays the evolution of the density of a plasma given the initial conditions listed in Table 5.As the solution evolves in time, the initial vortex splits into two vortices. Sharp gradients accumulate andthe vortex pattern becomes increasingly complex due to highly non-linear interactions between multipleintermediate shock waves travelling at different speeds. The result compares very well with results given inthe literature, e.g. [78, 79, 80], as well as with the solution of the Orszag-Tang MHD vortex obtained usingthe
MHD 8Wave , the
Bouchut5 , and the unsplit
USM implementations (shown in Fig. 12).26 . . . . . y t = 0 . t = 0 . t = 0 . . . . . . x . . . . . y t = 0 . . . . . . x t = 0 . . . . . . x t = 0 . . . . . . . . Fig. 11.
Orszag-Tang MHD vortex test: Density plots with superimposed magnetic field directions.
ES-Hybrid scheme with an adaptive grid resolution up to 512 × e.g. Fig. 1 of [78]. The lower right plot ( t = 0 .
5) can be compared to e.g.
Fig. 10of [79] and Fig. 14 of [80]. . . . . . x . . . . . y MHD 8Wave . . . . . x Bouchut5 . . . . . x USM . . . . . . . Fig. 12.
Orszag-Tang MHD vortex test: Density plots with an adaptive grid resolution up to 512 ×
512 at t = 0 . .4. MHD Rotor (2D) The MHD rotor problem [62] describes a rapidly spinning dense cylinder embedded in a magnetized,homogeneous medium at rest. Due to centrifugal forces, the dense cylinder is in non-equilibrium. As therotor spins with the given initial rotating velocity, the initially uniform magnetic field will wind up therotor. The wrapping of the rotor by the magnetic field leads to strong torsional Alfv´en waves launched intothe ambient fluid. Due to the onset and propagation of strong Alfv´en waves, this test is relevant for theunderstanding of star formation. The initial conditions are listed in Table 7. r ≤ r r ∈ ( r , r ) r ≥ r (cid:37) . . . f ( r ) 1 . p . . . B / √ π / √ π / √ πB . . . B . . . u − . y − . f ( r )∆ y . v . x . f ( r )∆ x . w . . . r = (cid:112) ( x − x center ) + ( y − y center ) ,∆ x = ( x − x center ), ∆ y = ( y − y center ),and f ( r ) = r − rr − r Domain size { x, y } min = { . , . }{ x, y } max = { . , . } Inner radius r = 0 . r = 0 . x -center x center = 0 . y -center y center = 0 . t max = 0 . γ = 1 . Table 7
Initial conditions and runtime parameters: MHD rotor test
This test demonstrates that our solver is able to resolve torsional Alfv´en waves, which is particularlyvisible in the plot of the magnetic pressure (right plot in Fig. 13). The Mach number M (see Fig. 13) showsthat the rotor is, up to a certain radial distance, still in uniform rotation. Beyond this radius, the rotorhas exchanged momentum with its environment and decelerated. In the left plot of Fig. 13 we presentthe magnetic field superimposed on the fluid density, (cid:37) . It is clearly seen that the magnetic field basicallymaintains its initial shape outside of the region of influence of the Alfv´en waves. Inside, the magnetic fieldis refracted by the MHD discontinuities. In Fig. 14 we present six snapshots of the evolution of the fluiddensity (as well as the AMR grid) up to the final time t = t max . This is shortly before the torsional Alfv´enwaves leave the computational domain and after the cylinder has rotated almost 180 ◦ . x y % x p tot x M x k B k .
342 6 .
033 9 .
723 0 .
674 1 .
747 2 .
821 0 .
670 2 .
011 3 .
352 0 .
466 1 .
320 2 . Fig. 13.
MHD rotor test: Adaptive grid resolution up to 512 × (cid:37) with overlayedmagnetic field, total pressure p tot = p + p mag , Mach number M with overlayed velocity vectors, and magnetic pressure p mag = (cid:107) B (cid:107) . This plot can directly be compared to Fig. 7 of [2], Fig. 14 of [79], and Fig. 2 of [62]. . . . . y t = 0 . t = 0 . . . . . y t = 0 . t = 0 . .
125 0 .
375 0 .
625 0 . x . . . . y t = 0 .
12 0 .
125 0 .
375 0 .
625 0 . x t = 0 .
15 10 % Fig. 14.
MHD rotor test: Density evolution on a logarithmic scale with superimposed AMR grid. Adaptive gridresolution up to level 8 (up to 1024 × .5. Comparison of Computational Efficiency (2D) We perform a memory and CPU time comparison on a uniform grid. We compare the new ES solverimplementation against the Bouchut 5 wave ( Bouchut5 ) [12], Powell’s 8 wave (
MHD 8Wave ) [28], and theunsplit staggered mesh (
USM ) [75, 76] solver implementations applied to the MHD Rotor problem, describedin the preceding section. For this test we use identical runtime parameters. The AMR grid is fixed to level 5.We present the results of the study in Table 8. Note that the effective computational costs are veryimplementation specific. We see that the ES solver uses slightly less memory than the other schemes. The ES solver needs more computational time per time step since the Runge-Kutta time integration schemeinvolves the full flux computation and spatial reconstruction procedure in each of the intermediate stages.The higher computational costs per time step can be – at least partially – compensated by choosing a larger CFL coefficient. We neglect this benefit here and run all simulations with a fixed
CFL number to give a faircomparison. For
USM we use a second order accurate Roe-type Riemann solver, previously used for thenumerical tests in [76]. Scheme Memory consumption (MB) CPU time (s) ES . Bouchut5 . . MHD 8Wave . . USM . . Table 8
Comparison of computational efficiency. The memory consumption is measured for the whole flash4 process, whilethe computational time corresponds only to the time used by the MHD solver as given by
FLASH ’s code performancesummary.
A particularly simple example of gravitational instability was discovered by Jeans [81]. This phenomenonis of great astrophysical interest in the context of star formation and cosmic structure growth. The configura-tion is a useful test for the coupling of multi-dimensional gravity to hydrodynamics in a computational code.The Jeans instability allows one to study the pressure dominated and gravity dominated limits as well as thenumerical method’s behaviour between the two limits. We start from an infinite homogeneous medium at restand consider a small perturbation in density. We shall suppose that the initial fluctuations in density and pres-sure take place adiabatically, so that p = γ(cid:37) . The initial conditions for this test are summarized in Table 9.We use the direct multigrid fast Fourier transform Poisson solver ( Grid/GridSolvers/Multigrid/fft ) forthe computation of the gravitational source term.Density (cid:37) [g cm − ] (cid:37) · [1 + δ ( x )]Pressure p [dyn cm − ] p · [1 + γδ ( x )]Perturbation δ ( x ) δ · cos ( k · x )Velocity u [cm s − ] Magnetic field B [G] with (cid:37) = 1 . · g cm − , δ = 1 · − ,and p = 1 . · dyn cm − Domain size [cm] { x, y } min = { . , . }{ x, y } max = { . , . } Boundary conditions all: periodicSimulation end time [s] t max = 5 . γ = 5 / k [cm − ] 2 π/ λ with λ = (0 . , , (cid:124) Table 9
Runtime parameters and initial conditions: Jeans Instability test (2D)
We obtain the dispersion relation of a self-gravitating fluid by solving the perturbed wave equations byplanar wave solutions in Fourier space. From the relation, ω = a k − πG(cid:37) , (4.1)30e define the Jeans wavenumber k J = √ πG(cid:37) a ≈ .
75 (4.2)with the gravitational constant, G = 6 . · − cm g − s − , and the adiabatic sound speed a = (cid:112) γ p /(cid:37) ≈ .
29 cm s − where the given numbers correspond to the initial conditions used. The Jeans wavenumbernumber is very important as it defines a scale on which gravitational effect become dominant in astrophysicalsystems. As long as k > k J , the perturbation is stable and oscillates with a real frequency of ω . This isthe case with our chosen initial conditions, as k ≈ > . ≈ k J . However, if k < k J , the perturbationgrows exponentially in time as ω is purely imaginary. An overdense region would become denser and denser,leading to gravitational collapse [82].We compute the oscillation frequency, ω , by measuring the time interval required for the energy toundergo exactly ten oscillations. The analytical expression for the kinetic, internal, and potential energiesare provided in Appendix A. − − − − − − − E i n t +2 . × E k i n . . . . . − . − . − . − . − . − . . E p o t ωt . . . . . ωt ES-Hybrid EC USM MHD 8Wave Bouchut5 exact
Fig. 15.
Jeans Gravitational Instability test: Plot of internal, kinetic, and gravitational energies. (left) The energychanges for all solvers agree well early in the computation. (right) As time progresses we see the dissipation of energyby each solver differs. The ES solver exhibits the least dissipation of the solvers tested and shows the best agreementwith the analytic solution. We use a fixed resolution of 64 ×
64 cells. We give the exact solution in Appendix A.
The resulting kinetic and internal (thermal) energies as functions of ωt are shown in Fig. 15. We performedthe simulations using a uniform resolution of 64 ×
64 grid cells and chose a small CFL coefficient of
CFL = 0 . EC scheme can be used as the solution is smooth. We note that the EC scheme shows essentially no dissipation even at the final time t = 5 .
0. All remaining solvers dissipate energyin some capacity. We see from Fig. 15 that the ES solver is considerably less dissipative than the other testedsolvers. Furthermore, the ES scheme agrees well with the analytic solution, while the other schemes fail tomaintain the exact oscillation period at later times and dissipate much of the energy of the dynamics.31 .7. MHD Blast Wave (2D, 3D) The two-dimensional version of the MHD blast wave problem was studied by [62]. We use an extendedthree-dimensional version to demonstrate the robustness of our scheme in simulations involving regimes withlow thermal pressures and high kinetic as well as magnetic energies in three dimensions. This test problemleads to the onset of strong MHD discontinuities, relevant to astrophysical phenomena where magnetic fieldscan have strong dynamical effects. It describes an initially circular pressure pulse. We choose here a relativemagnitude of 10 for comparison with [62]. The initial conditions used are listed in Table 10. r ≤ r r ∈ ( r , r ) r ≥ r (cid:37) . . . p . . . f ( r ) 0 . B √ π √ π √ π B . . . B . . . u . . . v . . . w . . . r = (cid:112) ( x − x center ) + ( y − y center ) ,and f ( r ) = r − rr − r Domain size { x, y, z } min = {− . , − . , − . }{ x, y, z } max = { . , . , . } Inner radius r = 0 . r = 0 . x center = (0 . , . , . t max = 0 . γ = 1 . Table 10
Initial conditions and runtime parameters: MHD blast wave test
The chosen initial conditions result in a very low plasma- β parameter, β = 2 p/ B ≈ . · − . TheMHD explosion emits fast magnetosonic shock waves propagating with high velocities. The explosion ishighly anisotropic and the displacement of the gas in the transverse y and z − direction is inhibited. Thisleads to the phenomenon that the explosion bubble is strongly distorted according to the initial magneticfield (in the x − direction). Furthermore, we see that the Mach number spans a broad range from 0 up to 42.As can be seen from the results in Fig. 16 showing the density, Mach number, magnetic pressure and theplasma- β at t = 0 .
01, the out-going blast wave shows no grid alignment effect. We also tested the new ES solver with relative pressure magnitudes of 10 and 10 without finding any numerical defects.In Fig. 17, we show zoomed linear density plots of the fast expanding shock front computed with our ES-Hybrid , the
Bouchut5 [12], the unsplit
USM [75, 76], and the
MHD 8Wave [28] solvers. As can be seen, theblast wave front has a elliptical shape as is expected due to the strong influence of the magnetic field. If werun the same simulation without magnetic fields, i.e. the hydrodynamic limit B = , we observe numericaldefects close to the Cartesian grid axes in all simulation runs except the one using the ES solver. Ismail etal. [45] showed that even schemes which have increased dissipation and do not show 1D shock instabilitiescan still suffer from the carbuncle phenomenon in multiple dimensions. We emphasize that there are nonumerical artifacts in the hydrodynamic limit when using the ES-Hybrid solver.Fig. 18 shows the evolution of the conserved quantities as well as the individual energies. We observethat the ES scheme is mass, momentum, and total energy conservative also in this extreme test. The otherMHD solvers fail to maintain the the conservation of total energy by a small amount. Also, the MHD 8Wave solver fails to preserve the conservation of momentum as expected with the addition of the Powell sourceterm. All solvers conserve mass up to machine precision. Fig. 19 shows a 3D variant of this test [70].32 . − . . . . y % M − . − . . . . x − . − . . . . y k B k − . − . . . . x p/ k B k . . . . . . . . − − − − − − Fig. 16.
MHD blast wave test: Adaptive grid resolution up to 512 × (cid:37) , Top right: Machnumber M , Lower left: magnetic pressure p mag = (cid:107) B (cid:107) , Lower right: Plasma- β = p/p mag . The density and magneticpressure plots can be compared to Fig. 13 of [79]. The plots of density and magnetic pressure can directly be comparedto Fig. 4 of [62]. − . − . . . . y ES-Hybrid Bouchut5 MHD 8Wave USM .
20 0 .
25 0 .
30 0 . x − . − . . . . y .
20 0 .
25 0 .
30 0 . x .
20 0 .
25 0 .
30 0 . x .
20 0 .
25 0 .
30 0 . x . . . . . . . . . Fig. 17.
MHD blast wave test: Linear density plots. Top: B = 100 / √ π , Bottom: B = 0 . . . . . . . E t o t a l − ( % u ) t o t a l × − ES-Hybrid USM MHD 8Wave Bouchut5 . . . . . . . % t o t a l E m ag . . . . . t E k i n . . . . . t E i n t .
009 0 . . . . .
009 0 . . . . . .
009 0 . . . . Fig. 18.
MHD blast wave test: Evolution of the total conserved quantities as well as the individual energies in the B x = 100 / √ π simulation run. It can easily be seen that the ES solver preserves the total energy well while the otherschemes fail to preserve the total energy. ig. 19. MHD blast wave test: Adaptive grid resolution up to 128 × ×
128 cells. Three-slice plot. Linear plot ofdensity, (cid:37) , using the same colour range as used in Fig. 16.
5. Conclusion
In this work we describe and test an implementation of a high-order, entropy stable numerical methodfor MHD problems. Entropy stable numerical fluxes serve as the core of the new entropy stable MHD solver.The implementation is implemented as a new module for the multi-physics
FLASH framework [26] offeringadaptive mesh refinement and large-scale, multi-processor simulations.The new scheme is implemented with three entropy stable variants: Roe type dissipation (
ES-Roe ), localLax-Friedrichs type dissipation (
ES-LLF ), and a hybrid dissipation term that uses a pressure switch to use
ES-Roe in smooth regions and the more dissipative
ES-LLF near pressure discontinuities (
ES-Hybrid ). Theintegration in time uses a strong stability preserving (SSP) Runge-Kutta method.The numerical approximation is built from an entropic point of view. Thus, at a given time, it is possibleto compute the current entropy density for all cells in the computational domain. We exploit this additionalknowledge of the entropy and reformulate the computation to guarantee positivity of the pressure whilemaintaining the conservation of the total energy of the system. This reformulation prevents non-physicalnegative pressures which can occur numerically if the internal energy is a small fraction of the total energy.We presented a variety of numerical results for the new entropy stable solver implementation for HDand MHD flows in multiple spatial dimensions. These numerical tests served to demonstrate the flexibilityand robustness of the new solver. The testing included a recently developed shock-tube experiment, smoothAlfv´en wave propagation, the Orszag-Tang MHD vortex, the MHD Rotor, and the strong MHD explosiontest. The coupling of gravity to our new solver has been tested using the Jeans instability test. We comparedthe physical aspects of the numerical results, CPU timing and memory consumption of the new entropystable scheme against the 8-wave, Bouchut 5-wave and unsplit MHD solver implementations already availablefor
FLASH .We found in these comparisons that the newly implemented entropy stable approximation was the mostaccurate in smooth regions of a flow. Also, it was shown that the entropy stable scheme was the only onethat maintains the conservation of total energy during the computation.The new
FLASH
MHD module is freely available upon contact with the corresponding author.35 cknowledgement
We thank the anonymous referee for their useful comments, which helped to improve this article. DominikDerigs and Stefanie Walch acknowledge the support of the Bonn-Cologne Graduate School for Physics andAstronomy (BCGS), which is funded through the Excellence Initiative, as well as the Sonderforschungs-bereich (SFB) 956 on the “Conditions and impact of star formation”. Stefanie Walch thanks the DeutscheForschungsgemeinschaft (DFG) for funding through the SPP 1573 “The physics of the interstellar medium”.This work has been partially performed using the Cologne High Efficiency Operating Platform for Sciences(
CHEOPS ) HPC cluster at the Regionales Rechenzentrum K¨oln (RRZK), University of Cologne, Germany.The software used in this work was developed in part by the DOE NNSA ASC- and DOE Office of ScienceASCR-supported FLASH Center for Computational Science at the University of Chicago. To create some ofthe figures, we have used the free visualization software YT [83]. References [1] D. W. Goldsmith, Thermal Instabilities in Interstellar Gas Heated by Cosmic Rays, Astrophysical Journal 161 (1970)41–54. doi:10.1086/150511 .[2] A. R. Winters, G. J. Gassner, Affordable, Entropy Conserving and Entropy Stable Flux Functions for the Ideal MHDEquations, Journal of Computational Physics 304 (2016) 72–108. doi:10.1016/j.jcp.2015.09.055 .[3] F. Bouchut, C. Klingenberg, K. Waagan, A multiwave approximate Riemann solver for ideal MHD based on relaxation. I:theoretical framework, Numerische Mathematik 108 (1) (2007) 7–42. doi:10.1007/s00211-007-0108-8 .[4] P. Chandrashekar, C. Klingenberg, Entropy stable finite volume scheme for ideal compressible MHD on 2-D cartesianmeshes, submitted to SIAM Journal of Numerical Analysis.[5] J. A. Rossmanith, High-order discontinuous Galerkin finite element methods with globally divergence-free constrainedtransport for ideal MHD, ArXiv e-prints: arXiv:1310.4251.[6] A. Dedner, F. Kemm, D. Kr¨oner, C.-D. Munz, T. Schnitzer, M. Wesenberg, Hyperbolic divergence cleaning for the MHDequations, Journal of Computational Physics 175 (2) (2002) 645–673. doi:10.1006/jcph.2001.6961 .[7] A. J. Christlieb, Y. Liu, Q. Tang, Z. Xu, Positivity-Preserving Finite Difference Weighted ENO Schemes with ConstrainedTransport for Ideal Magnetohydrodynamic Equations, SIAM Journal on Scientific Computing 37 (4) (2015) A1825–A1845. doi:10.1137/140971208 .[8] F. Huazheng, F. Xueshang, Splitting based scheme for three-dimensional MHD with dual time stepping, Chinese Journalof Space Science 35 (1) (2015) 9. doi:10.11728/cjss2015.01.009 .[9] C. M. Xisto, J. C. P´ascoa, P. J. Oliveira, A pressure-based high resolution numerical method for resistive MHD, Journal ofComputational Physics 275 (2014) 323 – 345. doi:10.1016/j.jcp.2014.07.009 .[10] M. Ersoy, F. Golay, L. Yushchenko, Adaptive multiscale scheme based on numerical density of entropy production forconservation laws, Central European Journal of Mathematics 11 (8) (2013) 1392–1415. doi:10.2478/s11533-013-0252-6 .[11] D. S. Spicer, H. Luo, J. E. Dorband, K. M. Olson, P. J. MacNeice, A new 3D, fully parallel, unstructured AMR MHD highorder Godunov code for modeling Sun-Earth connection phenomena, preprint submitted to Journal of Atmospheric andSolar-Terrestrial Physics (4 2013).[12] K. Waagan, C. Federrath, C. Klingenberg, A robust numerical scheme for highly compressible magnetohydrodynamics:Nonlinear stability, implementation and tests, Journal of Computational Physics 230 (9) (2011) 3331 – 3351. doi:10.1016/j.jcp.2011.01.026 .[13] X. Zhang, C.-W. Shu, On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equationson rectangular meshes, Journal of Computational Physics 229 (23) (2010) 8918 – 8934. doi:10.1016/j.jcp.2010.08.016 .[14] V. Wheatley, H. Kumar, P. Huguenot, On the role of Riemann solvers in discontinuous Galerkin methods for magnetohy-drodynamics, Journal of Computational Physics 229 (3) (2010) 660 – 680. doi:10.1016/j.jcp.2009.10.003 .[15] S. Li, A Simple Dual Implementation to Track Pressure Accurately, in: N. V. Pogorelov, E. Audit, G. P. Zank (Eds.),Numerical Modeling of Space Plasma Flows, Vol. 385 of Astronomical Society of the Pacific Conference Series, 2008, pp.273–278.[16] S. Li, H. Li, R. Cen, CosmoMHD: A Cosmological Magnetohydrodynamics Code, The Astrophysical Journal SupplementSeries 174 (2008) 1–12. arXiv:astro-ph/0611863 , doi:10.1086/521302 .[17] G. T´oth, The ∇ · B = 0 Constraint in Shock-Capturing Magnetohydrodynamics Codes, Journal of Computational Physics161 (2) (2000) 605 – 652. doi:10.1006/jcph.2000.6519 .[18] P. Janhunen, A positive conservative method for magnetohydrodynamics based on HLL and Roe methods, Journal ofComputational Physics 160 (2) (2000) 649–661. doi:10.1006/jcph.2000.6479 .[19] D. S. Balsara, D. Spicer, Maintaining pressure positivity in magnetohydrodynamic simulations, Journal of ComputationalPhysics 148 (1) (1999) 133 – 148. doi:10.1006/jcph.1998.6108 .[20] D. S. Balsara, Self-adjusting, positivity preserving high order schemes for hydrodynamics and magnetohydrodynamics,Journal of Computational Physics 231 (22) (2012) 7504–7517. doi:10.1016/j.jcp.2012.01.032 .[21] A. L. Zachary, A. Malagoli, P. Colella, A higher-order Godunov method for multidimensional ideal magnetohydrodynamics,SIAM Journal on Scientific Computing 15 (2) (1994) 263–284. doi:10.1137/0915019 .
22] D. Ryu, J. P. Ostriker, H. Kang, R. Cen, A cosmological hydrodynamic code based on the total variation diminishingscheme, The Astrophysical Journal 414 (1993) 1–19. doi:10.1086/173051 .[23] B. Einfeldt, C.-D. Munz, P. L. Roe, B. Sj¨ogreen, On Godunov-type methods near low densities, Journal of ComputationalPhysics 92 (2) (1991) 273 – 295. doi:10.1016/0021-9991(91)90211-3 .[24] B. Schmidtmann, B. Seibold, M. Torrilhon, Relations between WENO3 and Third-Order Limiting in Finite VolumeMethods, Journal of Scientific Computing, doi:10.1007/s10915-015-0151-z .[25] S. Gottlieb, On high order strong stability preserving Runge-Kutta and multi step time discretizations, Journal of ScientificComputing 25 (1) (2005) 105–128. doi:10.1007/s10915-004-4635-5 .[26] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zingale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Truran, H. Tufo,FLASH: An Adaptive Mesh Hydrodynamics Code for Modeling Astrophysical Thermonuclear Flashes, ApJS 131 (2000)273–334. doi:10.1086/317361 .[27] A. Dubey, L. B. Reid, K. Weide, K. Antypas, M. K. Ganapathy, K. Riley, D. J. Sheeler, A. Siegal, Extensible component-based architecture for FLASH, a massively parallel, multiphysics simulation code, Parallel Computing 35 (10–11) (2009)512 – 522. doi:10.1016/j.parco.2009.08.001 .[28] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, D. L. De Zeeuw, A solution-adaptive upwind scheme for idealmagnetohydrodynamics, Journal of Computational Physics 154 (2) (1999) 284–309. doi:10.1006/jcph.1999.6299 .[29] E. Tadmor, Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependentproblems, Acta Numerica 12 (2003) 451–512. doi:10.1017/s0962492902000156 .[30] L. D. Landau, Fluid Mechanics, Vol. 6, Pergamon, 1959.[31] S. Mishra, Entropy stable high-order schemes for systems of conservation laws, Modern Techniques in the NumericalSolution of Partial Differential Equations.[32] R. J. LeVeque, D. Mihalas, E. A. Dorfi, E. M¨uller, Computational Methods for Astrophysical Fluid Flow: Saas-FeeAdvanced Course 27. Lecture Notes 1997 Swiss Society for Astrophysics and Astronomy, Springer Science & BusinessMedia, 1998. doi:10.1007/3-540-31632-9 .[33] S. Walch, P. Girichidis, T. Naab, A. Gatto, S. C. O. Glover, R. W¨unsch, R. S. Klessen, P. C. Clark, T. Peters, D. Derigs,C. Baczynski, The silcc (simulating the lifecycle of molecular clouds) project – i. chemical evolution of the supernova-drivenism, MNRAS 454 (1) (2015) 238–268. doi:10.1093/mnras/stv1975 .[34] S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics,Matematicheskii Sbornik 89 (3) (1959) 271–306.[35] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Scientific and Statistical Computing 9 (6)(1988) 1073–1084. doi:10.1137/0909073 .[36] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review43 (1) (2001) 89–112. doi:10.1137/S003614450036757X .[37] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 3rd Edition, Springer,2009.[38] S. Godunov, Symmetric form of the equations of magnetohydrodynamics, Numerical Methods for Mechanics of ContinuumMedium 1 (1972) 26–34.[39] P. J. Dellar, A note on magnetic monopoles and the one-dimensional MHD Riemann problem, Journal of ComputationalPhysics 172 (1) (2001) 392–398. doi:10.1006/jcph.2001.6815 .[40] F. Ismail, P. L. Roe, Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks, Journal ofComputational Physics 228 (15) (2009) 5410–5436. doi:10.1016/j.jcp.2009.04.021 .[41] M. L. Merriam, An entropy-based approach to nonlinear stability, NASA Technical Memorandum 101086 (64) (1989)1–154.[42] T. J. Barth, Numerical methods for gasdynamic systems on unstructured meshes, in: D. Kr¨oner, M. Ohlberger, C. Rohde(Eds.), An Introduction to Recent Developments in Theory and Numerics for Conservation Laws, Vol. 5 of Lecture Notes inComputational Science and Engineering, Springer Berlin Heidelberg, 1999, pp. 195–285. doi:10.1007/978-3-642-58535-7_5 .[43] P. L. Roe, D. S. Balsara, Notes on the eigensystem of magnetohydrodynamics, SIAM Journal on Applied Mathematics56 (1) (1996) 57–67. doi:10.1137/S003613999427084X .[44] P. Chandrashekar, Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations, ArXiv e-prints arXiv:1209.4994 .[45] F. Ismail, P. L. Roe, H. Nishikawa, A proposed cure to the carbuncle phenomenon, in: H. Deconinck, E. Dick (Eds.),Computational Fluid Dynamics 2006, Springer Berlin Heidelberg, 2009, pp. 149–154. doi:10.1007/978-3-540-92779-2_21 .[46] J. J. Quirk, A contribution to the great Riemann solver debate, International Journal for Numerical Methods in Fluids18 (6) (1994) 555–574. doi:10.1002/fld.1650180603 .[47] K. Olson, PARAMESH: A parallel, adaptive grid tool, in: A. Deane, A. Ecer, G. Brenner, D. Emerson, J. McDonough,J. Periaux, N. Satofuka, D. Tromeur-Dervout (Eds.), Parallel Computational Fluid Dynamics 2005, Elsevier, Amsterdam,2006, pp. 341 – 348. doi:10.1016/B978-044452206-1/50041-0 .[48] A. Gatto, S. Walch, M.-M. M. Low, T. Naab, P. Girichidis, S. C. O. Glover, R. W¨unsch, R. S. Klessen, P. C. Clark,C. Baczynski, T. Peters, J. P. Ostriker, J. C. Ib´a˜nez-Mej´ıa, S. Haid, Modelling the supernova-driven ism in differentenvironments, MNRAS 449 (1) (2015) 1057–1075. doi:10.1093/mnras/stv324 .[49] S. C. O. Glover, P. C. Clark, Molecular cooling in the diffuse interstellar medium, MNRAS 437 (2014) 9–20. arXiv:1305.7365 , doi:10.1093/mnras/stt1809 .[50] K. Murawski, Analytical and Numerical Methods for Wave Propagation in Fluid Media (Stability, Vibration and Controlof Systems, Series A), World Scientific Pub Co Inc, 2003.
51] T. Plewa, E. M¨uller, The consistent multi-fluid advection method, Astronomy and Astrophysicsp 342 (1999) 179–191. arXiv:astro-ph/9807241 .[52] S. C. O. Glover, C. Federrath, M.-M. Mac Low, R. S. Klessen, Modelling CO formation in the turbulent interstellar medium,Monthly Notices of the Royal Astronomical Society 404 (2010) 2–29. arXiv:0907.4081 , doi:10.1111/j.1365-2966.2009.15718.x .[53] J. Barnes, P. Hut, A hierarchical O(N log N) force-calculation algorithm, Nature 324 (1986) 446–449. doi:10.1038/324446a0 .[54] P. M. Ricker, A Direct Multigrid Poisson Solver for Oct-Tree Adaptive Meshes, The Astrophysical Journal SupplementSeries 176 (2008) 293–300. arXiv:0710.4397 , doi:10.1086/526425 .[55] C. Altmann, Explicit discontinuous galerkin methods for magnetohydrodynamics, Ph.D. thesis, Universit¨at Stuttgart,Holzgartenstr. 16, 70174 Stuttgart (2012).URL http://elib.uni-stuttgart.de/opus/volltexte/2013/7998 [56] A. J. Chorin, The numerical solution of the navier-stokes equations for an incompressible fluid, Bull. Amer. Math. Soc.73 (6) (1967) 928–931.[57] J. U. Brackbill, D. C. Barnes, The effect of nonzero ∇ · b on the numerical solution of the magnetohydrodynamic equations,Journal of Computational Physics 35 (3) (1980) 426 – 430. doi:10.1016/0021-9991(80)90079-0 .[58] B. Marder, A method for incorporating gauss’ law into electromagnetic { PIC } codes, Journal of Computational Physics68 (1) (1987) 48 – 55. doi:10.1016/0021-9991(87)90043-X .[59] D. S. Balsara, Total variation diminishing scheme for adiabatic and isothermal magnetohydrodynamics, The AstrophysicalJournal Supplement Series 116 (1) (1998) 133–153. doi:10.1086/313093 .[60] R. K. Crockett, P. Colella, R. T. Fisher, R. I. Klein, C. F. McKee, An unsplit, cell-centered Godunov method for idealMHD, Journal of Computational Physics 203 (2005) 422–448. doi:10.1016/j.jcp.2004.08.021 .[61] C. R. Evans, J. F. Hawley, Simulation of magnetohydrodynamic flows - A constrained transport method, AstrophysicalJournal 332 (1988) 659–677. doi:10.1086/166684 .[62] D. S. Balsara, D. Spicer, A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fieldsin magnetohydrodynamic simulations, Journal of Computational Physics 149 (2) (1999) 270–292. doi:10.1006/jcph.1998.6153 .[63] D. S. Balsara, J. Kim, A Comparison between Divergence-Cleaning and Staggered-Mesh Formulations for NumericalMagnetohydrodynamics, The Astrophysical Journal 602 (2004) 1079–1090. doi:10.1086/381051 .[64] K. Waagan, A positive MUSCL-Hancock scheme for ideal magnetohydrodynamics, Journal of Computational Physics228 (23) (2009) 8609–8626. doi:10.1016/j.jcp.2009.08.020 .[65] P. Colella, P. R. Woodward, The piecewise parabolic method (PPM) for gas-dynamical simulations, Journal of computationalphysics 54 (1) (1984) 174–201. doi:10.1016/0021-9991(84)90143-8 .[66] M. Brio, C. C. Wu, An upwind differencing scheme for the equations of ideal magnetohydrodynamics, Journal ofComputational Physics 75 (2) (1988) 400–422. doi:10.1016/0021-9991(88)90120-9 .[67] D. Ryu, T. W. Jones, Numerical magnetohydrodynamics in astrophysics: Algorithm and tests for one-dimensional flow,The Astrophysical Journal 442 (1995) 228–258. doi:10.1086/175437 .[68] M. Torrilhon, Uniqueness conditions for Riemann problems of ideal magnetohydrodynamics, Journal of Plasma Physics69 (3) (2003) 253–276. doi:10.1017/s0022377803002186 .[69] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, ii, Journal ofComputational Physics 83 (1) (1989) 32–78. doi:10.1007/978-3-642-60543-7_14 .[70] Flash Center for Computational Science, University of Chicago, FLASH User’s Guide.URL http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p3.pdf [71] A. Susanto, High-Order Finite-Volume Schemes for Magnetohydrodynamics, Ph.D. thesis, University of Waterloo (2014).URL https://uwspace.uwaterloo.ca/handle/10012/8597 [72] R. Balasubramanian, K. Anandhanarayanan, An ideal magnetohydrodynamics solver with artificial compressibility analogydivergence-cleaning, 2014, pp. 1–8.URL [73] V. K. Chakravarthy, K. Arora, D. Chakraborty, A simple hybrid finite volume solver for compressible turbulence,International Journal for Numerical Methods in Fluids 77 (2015) 707–731. doi:10.1002/fld.4000 .[74] D. S. Balsara, Multidimensional Riemann problem with self-similar internal structure. Part I – Application to hyperbolicconservation laws on structured meshes, Journal of Computational Physics 277 (2014) 163–200. doi:10.1016/j.jcp.2014.07.053 .[75] D. Lee, A solution accurate, efficient and stable unsplit staggered mesh scheme for three dimensional magnetohydrodynamics,Journal of Computational Physics 243 (2013) 269–292. doi:10.1016/j.jcp.2013.02.049 .[76] D. Lee, A. E. Deane, An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics, Journal ofComputational Physics 228 (4) (2009) 952–975. doi:10.1016/j.jcp.2008.08.026 .[77] S. A. Orszag, C.-M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, Journal of FluidMechanics 90 (01) (1979) 129–143. doi:10.1017/s002211207900210x .[78] J. Balb´as, E. Tadmor, A Central Differencing Simulation of the Orszag Tang Vortex System, IEEE Transactions on PlasmaScience 33 (2005) 470–471. doi:10.1109/TPS.2005.845282 .[79] P. Londrillo, L. Del Zanna, High-Order Upwind Schemes for Multidimensional Magnetohydrodynamics, The AstrophysicalJournal 530 (2000) 508–524. arXiv:astro-ph/9910086 , doi:10.1086/308344 .[80] W. Dai, P. R. Woodward, A simple finite difference scheme for multidimensional magnetohydrodynamical equations,Journal of Computational Physics 142 (2) (1998) 331 – 369. doi:10.1006/jcph.1998.5944 .[81] J. H. Jeans, The Stability of a Spherical Nebula, Philosophical Transactions of the Royal Society of London A: Mathematical, hysical and Engineering Sciences 199 (312-320) (1902) 1–53. doi:10.1098/rsta.1902.0012 .[82] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover Books on Physics Series, Dover Publications, 1961.[83] M. J. Turk, B. D. Smith, J. S. Oishi, S. Skory, S. W. Skillman, T. Abel, M. L. Norman, yt: A Multi-code AnalysisToolkit for Astrophysical Simulation Data, The Astrophysical Journal Supplement 192 (2011) 9–+. arXiv:1011.3514 , doi:10.1088/0067-0049/192/1/9 . Appendix A. Analytic Solution of the Gravitational Instability Test
We list here the analytic solution for the gravitational instability test as given in [70, 82].Kinetic energy: E kin ( t ) = 18 (cid:37) δ | ω | L k (cid:18) − cos(2 ωt ) (cid:19) (A.1)Internal energy: E int ( t ) − E int (0) = − (cid:37) a δ L (cid:18) − cos(2 ωt ) (cid:19) (A.2)Potential energy: E pot ( t ) = − πG(cid:37) δ L k (cid:18) ωt ) (cid:19) (A.3)where L is the length of the computational domain, and k is the magnitude of the wave vector k . Appendix B. Diffusive Magnetic Field Correction
We present here the equations used for the implementation of the diffusive divergence error methoddescribed in Sec. 3.10. First and second derivatives are approximated by central second-order finite differences.Fig. B.1 illustrates the location of the different indices in the two dimensional case. i − i − i − iii i + 1 i + 1 i + 1 j − j − j − j j jj + 1 j + 1 j + 1 y → x → ∆ x ∆ y Fig. B.1.
Illustration of index locations in 2D ˜ B (cid:48) = ∂ x B + ∂ x ( ∂ y B ) + ∂ x ( ∂ z B ) = B ,i +1 ,j,k − B ,i,j,k + B ,i − ,j,k ∆ x + B ,i +1 ,j +1 ,k − B ,i +1 ,j − ,k y − B ,i − ,j +1 ,k − B ,i − ,j − ,k y x + B ,i +1 ,j,k +1 − B ,i +1 ,j,k − z − B ,i − ,j,k +1 − B ,i − ,j,k − z x (B.1)39 B (cid:48) = ∂ y B + ∂ y ( ∂ x B ) + ∂ y ( ∂ z B ) = B ,i,j +1 ,k − B ,i,j,k + B ,i,j − ,k ∆ y + B ,i +1 ,j +1 ,k − B ,i − ,j +1 ,k x − B ,i +1 ,j − ,k − B ,i − ,j − ,k x y + B ,i,j +1 ,k +1 − B ,i,j +1 ,k − z − B ,i,j − ,k +1 − B ,i,j − ,k − z y (B.2)˜ B (cid:48) = ∂ z B + ∂ z ( ∂ x B ) + ∂ z ( ∂ y B ) = B ,i,j,k +1 − B ,i,j,k + B ,i,j,k − ∆ z + B ,i +1 ,j,k +1 − B ,i − ,j,k +1 x − B ,i +1 ,j,k − − B ,i − ,j,k − x z + B ,i,j +1 ,k +1 − B ,i,j − ,k +1 y − B ,i,j +1 ,k − − B ,i,j − ,k − y z (B.3)˜ B = ∆ x ∆ y ∆ z ∆ x ∆ y + ∆ x ∆ z + ∆ y ∆ z (cid:18) ˜ B (cid:48) , ˜ B (cid:48) , ˜ B (cid:48) (cid:19) (cid:124) (B.4)In 2D, we instead have:˜ B = ∆ x ∆ y ∆ x + ∆ y (cid:18) ˜ B (cid:48) , ˜ B (cid:48) , (cid:19) (cid:124) (B.5)Note that in 2D computations the dark red highlighted parts are zero and can be neglected. Appendix C. Flowchart
In this section we provide flowcharts to detail and outline the workflow of the new ES MHD solver. Wedivide the flowchart description of the algorithm into three parts. The first, shown in Fig. C.4, depicts theglobal solver procedure for a single stage (of the possibly multi-stage) Runge-Kutta method. The secondflowchart in Fig. C.2 displays the workflow for the computation of the numerical fluxes. The third flowchartin Fig. C.3 displays the workflow for the update of the solution array.
StartSpatial reconstructionNumerical flux computationSource term computationCompute maximum ∆ t End StartRotation of solution vectorAdding up of fluxes and source termsComputation of entropy fluxUpdate of conserved quantitiesBackrotation of solution vectorEnd
Fig. C.2.
Flowchart showing the flux computationprocedure of the ES solver. For a full description seeFig. C.4. Fig. C.3.
Flowchart showing the solution update proce-dure of the ES solver. For a full description see Fig. C.4. tartGravity Accelerate all cellsGuard cell filling for all blocks Synchronizationbarrier + MPI com-municationEOS + Conversion to conservative space Prepare current blockFlux computation Calculate the fluxesYes NoMore blocks? Check whether addi-tional blocks have tobe processedAMR flux conservation Correct fluxes atjumps in level ofrefinement1 st stage? Compute ∆ t Compute time step1 st stage + RK? Backup Backup initial solutionduring first stageSolution update Apply the fluxes to thesolution arraySanity check Abort FLASH
Check validityEOS + Conversion to primitive space Postprocess currentblockYes NoMore blocks? Check whether addi-tional blocks have tobe processed ∇ · B cleaning Magnetic field diver-gence treatmentYes NoMore stages? Check whether addi-tional stages have to becomputedEnd R e p e a t e d f o r a ll b l o c k s lb=lb+1 Yes YesOK Failed ̺ ≤ p ≤ R e p e a t e d f o r a ll b l o c k s lb=lb+1 R e p e a t e d f o r a ll s t ag e s n ee d e db y t h e t i m e i n t e g r a t i o n s c h e m e stage=stage+1 Fig. C.4.
Flowchart showing the principle of operation of the ES solver. The steps “Flux computation” and “Solutionupdate” is depicted in separate flowcharts, shown in Fig. C.2 and C.3, respectively. ppendix D. Dimensional MHD equations The dimensional MHD equations are given by ∂∂t (cid:37)(cid:37) u E B + ∇ · (cid:37) u (cid:37) ( u ⊗ u ) + (cid:16) p + (cid:107) B (cid:107) µ (cid:17) I − B ⊗ B µ u (cid:16) E + p + (cid:107) B (cid:107) µ (cid:17) − B ( u · B ) µ B ⊗ u − u ⊗ B = 0 , (D.1)where the thermal pressure is related to the conserved quantities through the ideal gas law: p = ( γ − (cid:18) E − (cid:37) (cid:107) u (cid:107) − (cid:107) B (cid:107) µ (cid:19) . (D.2)The unit system is determined at compilation time through the user-definable parameter HY UNIT in the
ES.h parameter file. In non-dimensional units, (D.1) and (D.2) are identical to (2.1) and (2.3).The chosen unit system leads to different internal conversion factors within our implementation. Theresulting units of the simulation quantities are listed in Table D.1. Depending on the setting, the accordingvacuum permeability constant, µ , is chosen. Note that both the specific energy, E , and the pressure, p , areof the same units.Unit system non-dimensional SI CGS Length (cid:96) t (cid:37) − g cm − Velocities u − cm s − Specific energy E − erg cm − Pressure p − dyn cm − Magnetic field B µ = 1 4 π · − T m J − π G cm erg − Table D.1
Simulation units with different settings of the compilation-time parameter
HY UNIT
Appendix E. Error norms and the experimental order of convergence
The discrete L - and L -errors are defined by (cid:107) ∆ u ( t ) (cid:107) L = 1 N d N (cid:88) i,j =1 (cid:12)(cid:12)(cid:12)(cid:12) u exact i,j − u solution i,j (cid:12)(cid:12)(cid:12)(cid:12) , and (cid:107) ∆ u ( t ) (cid:107) L = (cid:118)(cid:117)(cid:117)(cid:116) N d N (cid:88) i,j =1 (cid:18) u exact i,j − u solution i,j (cid:19) . where N is the number or grid points in each direction, and d is the number of spatial dimensions. Aftercomputing the norms of the errors, we obtain the experimental order of convergence (EOC) usingEOC( i → j ) = log (cid:0) (cid:107) ∆ u i ( t ) (cid:107) (cid:1) − log (cid:0) (cid:107) ∆ u j ( t ) (cid:107) (cid:1) log (cid:0) N i /N j (cid:1) ..