A000 FVM

!!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!!
Eidgenössische Ecole polytechnique fédérale de Zurich

!!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!! Technische Hochschule Politecnico federale di Zurigo
!!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!! Zürich Swiss Federal Institute of Technology Zurich
Approximate Riemann solvers and stable

high-order finite volume schemes for
multi-dimensional ideal MHD
F.G. Fuchs∗ , A.D. McMurry∗ , S. Mishra, N.H. Risebro∗ and K. Waagan†
Research Report No. 2009-37

November 2009
Seminar für Angewandte Mathematik

Eidgenössische Technische Hochschule
CH-8092 Zürich
Switzerland
∗
Center of Mathematics for Applications, University of Oslo, 1053 Blindern, Norway
†
Center for Scientific Computation and Mathematical Modeling, The University of
Maryland, CSCAMM, College Park, MD 20742-3289, USA
APPROXIMATE RIEMANN SOLVERS AND STABLE HIGH-ORDER FINITE VOLUME
SCHEMES FOR MULTI-DIMENSIONAL IDEAL MHD
F. G. FUCHS, A. D. MCMURRY, S. MISHRA, N. H. RISEBRO, AND K. WAAGAN
Abstract. We design stable and high-order accurate finite volume schemes for the ideal MHD equations in
multi-dimensions. We obtain excellent numerical stability due to some new elements in the algorithm. The
schemes are based on three- and five-wave approximate Riemann solvers of the HLL-type, with the novelty that
we allow a varying normal magnetic field. This is achieved by considering the semi-conservative Godunov-Powell
form of the MHD equations. We show that it is important to discretize the Godunov-Powell source term in the
right way, and that the HLL-type solvers naturally provide a stable upwind discretization. Second-order versions
of the ENO- and WENO-type reconstructions are proposed, together with precise modifications necessary to
preserve positive pressure and density. Extending the discrete source term to second order while maintaining
stability requires non-standard techniques, which we present. The first- and second-order schemes are tested
on a suite of numerical experiments demonstrating impressive numerical resolution as well as stability, even on
very fine meshes.
1. Introduction
Many interesting problems in astrophysics, solar physics and engineering involve macroscopic plasma models
and are usually described by the equations of ideal magneto-hydrodynamics (MHD). In these models, the
variables of interest are the mass density of the plasma ρ, the velocity field u = (u1 , u2 , u3 )T , the magnetic
field B = (B1 , B2 , B3 )T , the pressure p and the total energy E. The unknowns obey the following conservation
(balance) laws (see [32] for details),
(i.) Conservation of mass: Mass of a plasma is conserved, resulting in
ρt + div(ρu) = 0.
(ii.) Faraday’s law: The magnetic flux across a surface S bounded by a curve δS is given by Faraday’s law
! !
d
− B · dS = E · dl.
dt S δS
By using Stokes theorem and the fact that the electric field in a co-moving frame is zero and assuming
zero resistivity, Faraday’s law leads to
(1.1) Bt + curl(B × u) = −u(divB).
The above equation is termed the magnetic induction equation and can also be written in the divergence
form:
Bt + div(u ⊗ B − B ⊗ u) = −u(divB).
(iii.) Conservation of Momentum: In differential form, the conservation of momentum is
(ρu)t + div(ρu ⊗ u + pI) = J × B,
where J denotes the current density and I the 3 × 3 identity matrix. The Lorentz force exerted by the
magnetic field is given by J × B. Under the assumptions of ideal MHD, Ampere’s law expresses the
current density as
J = curl(B).
Date: November 16, 2009.

Key words and phrases. Conservation laws, MHD, divergence constraint, upwinded source terms, high-order schemes.
1
2 F. G. FUCHS, A. D. MCMURRY, S. MISHRA, N. H. RISEBRO, AND K. WAAGAN
Using standard vector identities results in the following semi-conservative form,

" #
1 2
(ρu)t + div ρu ⊗ u + (p + B )I − B ⊗ B = −B(divB).
2
(iv.) Conservation of energy: Defining the hydrodynamic energy of an ideal gas as
p 1
E hd = + ρu2 ,
γ−1 2
and using the conservation of this energy results in
Ethd + div((E hd + p)u) = J · (B × u).
The right hand side represents the change in energy due to the magnetic field. By using standard vector
identities and Ampere’s law, we obtain
∂B
J · (B × u) = (B · − (u · B)(divB) − div(B · B)u − (u · B)B).
∂t
Defining the total energy of the plasma as E = E hd + 12 B2 , energy conservation takes the form
1
Et + div((E + p + B2 )u − (u · B)B) = −(u · B)(divB).
2
Combining all the above balance laws leads to the semi-conservative form of the ideal MHD equations:
ρt + div(ρu) = 0,
1 2
(ρu)t + div(ρu ⊗ u + (p + |B| )I − B ⊗ B) = −B(divB),
(1.2) 2
Bt + div(u ⊗ B − B ⊗ u) = −u(divB),
1
Et + div((E + p + |B|2 )u − (u · B)B) = −(u · B)(divB).
2
The semi-conservative form (1.2) is also called the Godunov-Powell form, and the source on the right-hand side
of (1.2) is termed the Godunov-Powell source term.
Magnetic monopoles have not been observed in nature (although their existence has been hypothesized in a
number of quantum regimes by both the unified field theory as well as string theory). Hence, it is common to
assume that the magnetic field is solenoidal, i.e., it satisfies the divergence constraint,
(1.3) div(B) ≡ 0.
Under this constraint, the source terms in (1.2) are zero and the constraint is explicitly added to the equations
to obtain the conservative form of the ideal MHD equations:
ρt + div(ρu) = 0,
1 2
(ρu)t + div(ρu ⊗ u + (p + |B| )I − B ⊗ B) = 0,
2
(1.4) Bt + div(u ⊗ B − B ⊗ u) = 0,
1
Et + div((E + p + |B|2 )u − (u · B)B) = 0,
2
div(B) = 0.
Taking divergence on both sides of the magnetic induction equation (1.1) yields,
(1.5) (divB)t + div(u(divB)) = 0.
This means that any solenoidal initial magnetic field remains divergence free. Hence, for smooth solutions, the
semi-conservative form (1.2) is equivalent to the standard form (1.4), provided that the initial magnetic field is
divergence free.
Despite their formal equivalence, the two forms, (1.4) and (1.2), differ in some respects. We believe that
it is more natural to consider and discretize the semi-conservative form (1.2), since the derivation from first
principles gives (1.2). Furthermore, (1.2) is Galilean invariant, whereas the standard form (1.4) is not.
APPROXIMATE RIEMANN SOLVERS FOR MHD 3
From a mathematical perspective, the semi-conservative form (1.2) was shown to be symmetrized by the
physical entropy in [19], leading to stability estimates [6]. The standard form (1.4) is not symmetrizable and it
is not clear if it is possible to obtain any well-posedness results for (1.4) whereas estimates on the entropy may
pave the way for obtaining well-posedness results for the symmetrizable form (1.2) ,[6].
The standard form (1.4) (semi-conservative form (1.2)) is a system of conservation (balance) laws. Eigenvalue
analysis, see [33, 6], shows that the system is hyperbolic but not strictly hyperbolic. Solutions are quite
complicated and contain interesting discontinuities like shock waves, contact discontinuities, compound and
intermediate shocks. Even for “simple” initial value problems, such as the Riemann problem, we do not have
existence or well-posedness results. This means that numerical simulations are the main tools to study solutions
of these equations.
Finite volume schemes are among the most efficient and widely used numerical methods for the numerical
solution of conservation (balance) laws. In these methods (see [24]), the computational domain is divided into
cells, and cell-averages of the conserved quantities are evolved by integrating the balance law over the cell.
The update requires numerical fluxes, defined in terms of exact or approximate solutions of Riemann problems
(along the normal direction) at each cell interface. Higher-order spatial accuracy is recovered by employing non-
oscillatory piecewise polynomial reconstructions like second-order TVD [43], higher-order ENO [22] and WENO
[36] methods. Higher-order temporal accuracy results from using stability-preserving Runge-Kutta methods
[21].
Finite volume schemes for the MHD equations have been employed with a fair amount of success. The
Riemann problem is too complicated to solve exactly [40], and approximate Riemann solvers are employed to
define numerical fluxes. Linearized solvers of the Roe-type [33, 12] have been devised but often give negative
pressures and densities, see [18] for examples. An alternative is to use approximate Riemann-solvers of the HLL-
type [23]. These solvers approximate the wave-structure of the Riemann problem (consisting up to eight waves),
by a smaller number of waves. Three-wave HLL solvers have been designed in [26, 20, 8] and five-wave solvers
in [29, 8]. Some of these solvers [20, 29, 8] are proved to preserve positive pressures and densities. They also
typically ensure that the second law of thermodynamics is not violated, referred to as entropy stability, which
is not easily achieved with linearized solvers. Numerical results showing the robustness of these HLL-solvers,
particularly in one space dimension, have been presented, see [29, 8]. We remark that Roe-type Riemann solvers
are compatible with the algorithmic framework presented here, and would conceivably benefit from it.
However, the extension of one-dimensional numerical fluxes to multi-dimensional MHD in standard form
(1.4) is not straightforward. The divergence constraint (1.3) in one-dimension implies that the normal magnetic
field must be a constant in space. HLL-type solvers like the ones described in [20, 29] use this information in
their definitions of speeds and states. For multi-dimensional MHD, the magnetic field in each normal direction
is no longer constant. Consequently, it is not trivial to extend the HLL-solvers in this case. One possible
solution consists in using an average of the normal magnetic field across each interface in the expressions. This
somewhat arbitrary choice may destroy the stability properties of the solvers.
Another highly non-trivial aspect in several dimensions is the treatment of the divergence constraint divB = 0.
Standard finite volume schemes will generate divergence errors, and these can induce instabilities, see [42]. A
popular method to remove divergence is the projection method of [10]. In this method the non-solenoidal
approximate magnetic field is corrected at each time step by using an elliptic solver. The method leads to
solenoidal fields, but is computationally expensive on account of the elliptic solver. Other stability issues with
this method have been discussed in [34, 42].
A popular divergence cleaning method consists in staggering the discretizations of the velocity and magnetic
field, leading to methods with a (discrete) divergence free magnetic field. Variants of this method have been
proposed in [16, 13, 2, 27, 34, 35, 37] and references therein. Staggering of the variables leads to complications
when parallelizing the method, and in designing variable grid methods. Unstaggered variants of this approach
have been proposed in [38, 39, 28]. A thorough comparison of divergence cleaning methods was performed in
[42]. It must be noted that controlling a particular discrete form of the divergence does not generally lead to
any control on other discrete forms as the solutions of the ideal MHD equations are not necessarily smooth.
Therefore, in the presence of shocks, the values of a particular discrete form do not necessarily imply anything
about the quality of the simulation results. Examples illustrating this phenomenon are presented in [17].
Another potential problem with both the staggering and projection techniques lies in their numerical stability.
The nonlinear stability analysis techniques (i.e. entropy stability, positivity) that have been quite successful for
finite volume schemes do not apply. In a recent paper [18], we provided examples where the projection method
was quite stable on coarse meshes, but exhibited instabilities when the mesh was refined. We have also observed
similar behavior for some staggered mesh methods.
A third alternative for divergence cleaning was proposed in [31] and consists of discretizing the semi-
conservative Godunov-Powell form (1.2). In [31, 32], a linearized solver is used to define numerical fluxes
and a simple centered discretization is proposed for the Godunov-Powell source term in (1.2). Note that (1.5)
suggests the initial divergence errors will be transported out of the domain by the flow. In a recent paper [17],
examples were constructed showing that a centered discretization of the Godunov-Powell source term can lead
to numerical instabilities, even for the simple case of the linear magnetic induction equation (1.1). Hence, the
Godunov-Powell source term has to be suitably “upwinded” to obtain stable discretizations. A related method
is the Generalized Lagrange multiplier method of [14].
The above discussion brings out the relative strengths and weaknesses of the three available methods for
divergence cleaning. All the three methods might lead to problems with numerical stability. It is mostly
manifested on very fine meshes as the numerical diffusion is not sufficient to stabilize the schemes. Since the
projection method can be computationally expensive and staggering based methods are quite complicated to
implement on parallel computers (and with adaptive mesh refinement), we focus on a stable approximation of the
semi-conservative Godunov-Powell form (1.2). The fact that this form is Galilean invariant and symmetrizable
makes it particularly desirable to discretize. We are also motivated by related approaches in [8, 44] and [17],
which contain stability estimates not available for staggered schemes, as well as strong numerical evidence of
robustness. Note that a straightforward central discretization of the source term may be unstable, so we will
upwind the source term in a suitable manner.
Our aim in this paper is to design a stable and high-order accurate finite volume scheme for the semi-
conservative form (1.2) of ideal MHD in multi-dimensions. Our numerical algorithm consists of the following
ingredients
• We derive HLL-type three-wave and five-wave approximate Riemann solvers for the semi-conservative
form (1.2). The fluxes are defined in terms of approximate solutions to Riemann problems (in the
normal direction) at each interface. The main difference between existing solvers and our approach lies
in the fact that we allow the normal magnetic field to vary across the interface. The resulting solvers
extend the three-wave HLLC solver of [20] and the highly popular HLLD five-wave solver of [29] to
non-constant normal magnetic fields and hence trivially to multi-dimensions.
• We discretize the Godunov-Powell source term in (1.2) by using the states and speeds of the HLL solvers
to calculate the source term in each direction. This is simply and naturally achieved by taking the usual
cell averages. Thus, the source term is automatically upwinded.
• Second-order spatial accuracy is obtained by designing suitable ENO and WENO-type reconstructions.
Standard ENO-WENO reconstructions have been used to obtain high-order accurate schemes for the
MHD equations, see [41, 5, 4] and references therein. However, these reconstructions have to be modified
to ensure that the resulting schemes preserve positive pressures and densities. We rely on the results of
[7, 30] and of a recent paper [44] to design these modifications.
• The upwind discrete Godunov-Powell source term is extended to second order with techniques from [1]
and [44].
• Second-order temporal accuracy is obtained by using Runge-Kutta methods.
The above ingredients are combined to obtain a second-order finite volume scheme for MHD equations (1.2)
in multi-dimensions. The resulting schemes are very simple to implement. Although we are unable to provide
rigorous stability proofs, we validate the resulting schemes on a wide variety of benchmark numerical experi-
ments. The numerical results on a sequence of meshes (including very fine meshes) demonstrate that both the
first- and the second-order versions of our schemes are numerically stable.
As stated earlier, stability is the key to numerically resolve details on very fine meshes. We would like to
mention that highly resolved solutions are not widely reported in the literature, leaving the stability at those
resolutions open for questioning.
Related discretizations of a modified (partial) form of (1.2) based on operator splitting and on relaxation
were proposed in [18] and [9] respectively. The solvers of [9] were proved to be positivity preserving and entropy
stable. Second-order positivity preserving extensions of this approach have been presented and tested in a recent
paper [44].
The rest of this paper is organized as follows: the numerical schemes are presented in Section 2 and the
numerical experiments are reported in Section 3.
2. Numerical Schemes
For notational simplicity, we focus on the semi-conservative form of the MHD equations (1.2) in two space
dimensions,
(2.1) Wt + f (W)x + g(W)y = s1 (W, Wx ) + s2 (W, Wy ),
where
T
W = (ρ, ρu1 , ρu2 , ρu3 , B1 , B2 , B3 , E) ,
is the vector of conserved variables, and the fluxes are given by
   
ρ ρ
 B2   ρu1 u2 − B1 B2 
 ρu21 + π1 − 21   
   B 2

 ρu1 u2 − B1 B2   ρu22 + π2 − 22 
   
 ρu1 u3 − B1 B3   ρu1 u3 − B1 B3 
f (W) =   , g(W) =  ,
 0   u2 B1 − u1 B2 
   
 u1 B2 − u2 B1   0 
   
 u1 B3 − u3 B1   u2 B3 − u3 B2 
2 2
B B
(E + π1 )u1 − u1 21 − B1 (u2 B2 + u3 B3 ) (E + π2 )u2 − u2 22 − B2 (u1 B1 + u3 B3 )
where we have defined
B22 + B32 B 2 + B32
(2.2) π1 = p + , π2 = p + 1 .
2 2
Similarly, the Godunov-Powell source terms in (1.2) can be written explicitly as
   
0, 0,
 B2   −B1 (B2 )y 
 −( 21 )x   
   B22 
 −B2 (B1 )x   −( 2 )y 
   
 −B 3 (B1 )x   −B 3 (B 2 )y 
s1 (W, Wx ) =   , s2 (W, Wy ) =  
 −u1 (B1 )x   −u1 (B2 )y 
   
 −u2 (B1 )x   −u2 (B2 )y 
   
 −u3 (B1 )x   −u3 (B2 )y 
B12 B22
−u1 ( 2 )x − (u2 B2 + u3 B3 )(B1 )x −u2 ( 2 )y − (u1 B1 + u3 B3 )(B2 )y
B2 B2
Note that we have used the chain rule B1 (B1 )x = ( 21 )x and B2 (B2 )y = ( 22 )y . While true for smooth solutions,
this formula may no longer hold when the magnetic field has discontinuities. However, using this formulation
proved to be robust in practice.
We approximate (2.1) in a domain x = (x, y) ∈ [Xl , Xr ] × [Yb , Yt ]. For simplicity, the domain is discretized
by a uniform grid in both directions with mesh sizes ∆x and ∆y, respectively. We set xi = Xl + i∆x and
yj = Yb + j∆y, and Ii,j = [xi−1/2 , xi+1/2 ) × [yj−1/2 , yj+1/2 ). The cell average of the unknown vector W at time
tn in the cell Ii,j is denoted Wi,j
n
.
A standard finite volume scheme (first-order in both space and time) is obtained by integrating the balance
law (2.1) over the cell Ii,j and the time interval [tn , tn+1 ] with tn+1 = tn + ∆tn , where the time-step ∆tn is
determined by a suitable CFL condition. The resulting fully-discrete form of the scheme is
∆tn n ∆tn n
(2.3) n+1
Wi,j = Wi,j
n
− (Fi+1/2,j − Fni−1/2,j ) − (Gi,j+1/2 − Gni,j−1/2 ) + ∆tn (S1i,j + S2i,j ).
∆x ∆y
The flux Fni+1/2,j = F(Wi,j

n
, Wi+1,j
n
) and the source S1i,j are hence determined from a solution to the following
Riemann problem,
*
WL x < 0,
(2.4) Wt + f (W)x = s (W, Wx ), W(x, 0) =
1
WR x > 0.
Similarly, the flux G and the approximation of the Godunov-Powell source term S2 are given in terms of a
Riemann problem in the y-direction. The Riemann problem (2.4) has an intricate solution involving up to eight
waves. Therefore the solution to (2.4) is in practice replaced by a so called approximate Riemann solver (see
e.g. [24]). The purpose of this section is to derive two such solvers. They will be given as functions of x/t
(as the exact solution of (2.4) is self-similar), and we must ensure local conservation in order to end up with a
scheme of the semi-conservative form (2.3).
2.1. HLL three-wave solver. To derive this solver we follow the approach of [20]. The approximate solution
and fluxes for (2.4) are given by
 

 WL if xt ≤ sL , 
 FL if xt ≤ sL ,

W∗ if s < x < s , 
F∗ if s < x < s ,
L M L M
(2.5) WHLL3 = L t FHLL3 (WL , WR ) = L t


 W ∗
R if sM < x
t < sR , 

 F ∗
R if sM < x
t < sR,
 
WR if sR ≤ t ,x
FR if sR ≤ t .x
Note that we do not enforce F = f (W). The reason is that we allow π1 of (2.2) to be a free variable. It play a
B 2 +B 2
role similar to the relaxation pressure in [8]. For consistency we set π1 = p + 2 2 3 in FL and FR . The outer
wave speeds sL and sR model the fast magneto-sonic waves and are defined as in [20, 15], i.e.,
(2.6) sL = min {u1L − cf L , ū1 − c̄f } , sR = max {u1R + cf R , ū1 + c̄f } ,
where ū1 and c̄f are the normal velocity and the fast wave speed of the Jacobian matrix A((WL + WR )/2)
respectively. This choice is important for numerical stability and accuracy.
In order to describe the solver, we need to determine the speed of the middle wave sM and the intermediate
states WL ∗
, WR∗
. We follow [20] in letting the middle wave model a material contact discontinuity. Hence, the
velocity field and the tangential magnetic fields are assumed to be constant across the middle wave. This allows
defining u∗ = u∗L = u∗R , B2∗ = B2L ∗
= B2R∗
and B3∗ = B3L∗
= B3R
∗
. Furthermore, the difference in our solver
and the three-wave solver of [20] lies in the fact that we consider a non-constant normal magnetic field B1 . The
normal magnetic field B1 only jumps across the middle wave (modeling the linear degenerate ”divergence wave”
implied by (1.5)) and is constant across the outer waves.
We will impose local conservation across each wave to determine the various states. Local conservation across
the outermost waves means that
(2.7) ∗
sL WL − F∗L = sL WL − FL , sR WR − FR = sR WR
∗
− F∗R .
Conservation across the middle wave sM involves taking the source term s1 in (2.4) into account. The conser-
vation relation is given by,
(2.8) ∗
sM WR ∗
− sM WL = F∗R − F∗L + s1,∗
where
 
0
 2 2 
 − (B1R ) −(B
2
1L )

 
 −B2∗ (B1R − B1L ) 
(2.9) s1,∗ = ,
 −B3∗ (B1R − B1L ) 
 
 −u∗ (B1R − B1L ) 
2 2
∗ (B1R ) −(B1L )
−u1 2 − (u∗2 B2∗ + u∗3 B3∗ )(B1R − B1L )
amounts to integrating the source s1 in (2.4) across the middle wave (as described in the next section). The
above expression follows from the assumption that B1 jumps only across the middle wave while the velocity
field and tangential components of the magnetic field remain constant. The use of the source term in the above
conservation relations is the key difference in our approach and the one used in [20].
For any middle speed sM , a straightforward application of the conservation relations (2.7) determines unique
values of intermediate densities given by
u1θ − sθ
(2.10) ρ∗θ = ρθ , θ ∈ {L, R}
sM − sθ
Using conservation across all the three waves (adding (2.7) and (2.8)) results in the global conservation relation,
(2.11) FR − FL = sR WR − sL WL + (sM − sR )WR
∗
+ (sL − sM )WL
∗
+ s1,∗ .
We can use the intermediate density states (2.10) and global conservation (2.11) to obtain
π1R − π1L + ρR u1R (u1R − sR ) − ρL u1L (u1L − sL )
sM = u∗1 = .
ρR (u1R − sR ) − ρL (u1L − sL )
Similarly, one uses local conservation (2.7) across the two outer waves to obtain the intermediate “relaxed”
pressures,
(2.12) ∗
π1θ = π1θ + ρθ (u1θ − sθ )(u1θ − sM ),
for θ ∈ {L, R}. Note that conservation across the middle wave automatically implies that π1L
∗
= π1R
∗
, and that
(2.12) confirms this assertion. The next step is to determine the tangential velocity and magnetic field. Using
global conservation across the wave fan (2.11), we obtain that the intermediate values u∗σ and Bσ∗ satisfy the
following two linear equations,
αu∗σ − βBσ∗ = cσ , −βu∗σ − ζBσ∗ = dσ , σ ∈ {2, 3},
where
cσ = ρR uσR (u1R − sR ) − ρL uσL (u1L − sL ) − (B1R BσR − B1L BσL ),
(2.13) dσ = BσR (sR − u1R ) − BσL (sL − u1L ) − (B1L uσL − B1R uσR ),
α = ρR (u1R − sR ) − ρL (u1L − sL ), ζ = sR − sL , β = B1R − B1L .
Solving the linear system (2.13), the intermediate tangential components of velocity and magnetic field are
obtained as,
ζcσ − βdσ −αdσ − βcσ
(2.14) u∗σ = , Bσ∗ = .
αζ + β 2 αζ + β 2
Remark 2.1. In general, the denominator; αζ + β 2 , in (2.14) can become small, leading to a degeneracy in the
states. A simple calculation shows that αζ + β 2 '= 0 if (ρR cR
f + ρL cf )(sR − sL ) > (B1R − B1L ) . This condition
L 2
can be ensured by “widening” the wave fan slightly by modifying the fast wave speeds in (2.6). The resulting
conditions are
1 1
(2.15) sR ≥ u1R + (max((u1L − u1R ), 0)) + / cf R , sL ≤ u1L − (max((u1L − u1R ), 0)) − / cf L ,
2 2
where 0
γpθ B 2
B 2
+ B 2 γpθ + |B|2θ 2 γpθ B1θ2
c2f θ =
/ + 1θ (1 + +) + 2θ 3θ
+ ( ) −4 , θ ∈ {L, R},
ρθ ρθ ρθ ρθ ρ2θ
for some small positive +. Using the conditions (2.15) to widen the wave fan ensures that the denominator
αζ + β 2 is never zero and the states are well defined.
Finally, the intermediate total energy states are determined by local conservation relations (2.7) to be
2
B1θ
Eθ (u1θ − sθ ) + π1θ u1θ − π1θ
∗
sM + 2 (u1θ
− sM ) + B1θ (B2θ u2θ + B3θ u3θ − B2θ
∗ ∗
u3θ )
∗ ∗
u2θ − B3θ
Eθ∗ = ,
sM − sθ
for θ ∈ {L, R}. Hence, all the intermediate states are determined explicitly. The intermediate fluxes are obtained
in terms of the intermediate states by local conservation (2.7),
F∗L = FL + sL (WL
∗
− WL ), F∗R = FR + sR (WR
∗
− WR ).
Combining the above expressions for the states and the fluxes, we write down our explicit flux formula for the
three-wave solver as


 Fi,j , if sL,i+1/2,j > 0,

F∗ , if sL,i+1/2,j ≤ 0 ∧ sM,i+1/2,j ≥ 0,
(2.16) FH
i+1/2,j =  ∗
3 i,j

 Fi+1,j , if sM,i+1/2,j < 0 ∧ sR,i+1/2,j ≥ 0,

Fi+1,j , if sR,i+1/2,j < 0.
Note that this may be discontinuous at sM,i+1/2,j = 0 according to (2.8). Hence our choice of FH3 in that case
is merely a convention. It is the proper addition of the source term which ensures that the scheme is continuous.
Remark 2.2. If we assume that the normal magnetic field B1 is constant i.e., B1L = B1R , then the three-wave
solver defined above reduces to the three-wave solver derived in [20]. Hence, our three-wave solver extends the
standard three-wave solver of [20] for the case of a non-constant normal magnetic field.
2.2. Discretization of the Godunov-Powell source term. In this section we explain why we must have the
jump condition (2.8) across the middle wave, and specify the discrete source S1,n
i,j in (2.3). The discrete source
must be consistent with the Godunov-Powell source term in x-direction s1 (W, Wx ). It will be determined
from our solution of the Riemann problem (2.4) along the x-direction at the cell interfaces (xi+1/2 , yj ). The
HLL three-wave approximate Riemann solver of the previous section provide us with the assumptions we need:
The normal magnetic field jumps only across the contact-discontinuity modeled by the middle wave, while the
velocity field and the tangential components of the magnetic field are constant across the middle wave.
Let T be a quantity that is constant across the middle wave and takes the value T ∗ , then
(2.17) (T Bx1 )(x, t) = T ∗ (B1R − B1L )δ(x + tu∗1 ),
where δ denotes the Dirac delta function. If we assume that |u∗1 | ∆tn ≤ ∆x, integrating T Bx1 over (−∆x, 0)
yields
! 0
(2.18) T Bx1 dx = T ∗ (B1R − B1L )1{u∗ <0} ,
1
−∆x
where 1A denotes the characteristic function of the set A. Integration over (0, ∆x) leads to
! ∆x
T Bx1 dx = T ∗ (B1R − B1L )1{u∗ >0} ,
1
0
under the same restriction, |u∗1 | ∆tn
≤ ∆x. Similarly, by again using the assumption that B1 jumps only across
the contact and T remains constant across it, we obtain that
! 0
B2 B 2 − B1L
2
T ( 1 )x dx = T ∗ 1R 1{u∗ <0} ,
−∆x 2 2 1
and
! ∆x
B12 B 2 − B1L
2
(2.19) T( )x dx = T ∗ 1R 1{u∗ >0} .
0 2 2 1
Hence, we can derive (2.8) from (2.18)-(2.19) by observing that we must have
! ∆x
s1,∗ = s1 (WH3 , WxH3 )dx.
−∆x
The final scheme is defined by evolving the piecewise constant function Wi,j according to the approximate
Riemann solver (2.5), and then taking the cell average of the conserved quantities. Hence, the scheme is
determined by (2.7)-(2.8), yielding (2.16), and
(2.20) S1,n 1,∗ 1,∗
i,j = si−1/2,j 1{(sM,i−1/2,j ≥0)} + si+1/2,j 1{(sM,i+1/2,j <0)} ,
where s1,∗
i±1/2,j is defined in (2.9). For the case that sM,i+1/2,j = 0 our choice here was dictated by our choice in
(2.16). Integration along the y-direction is taken care of by the midpoint rule.
We emphasize that the discrete Godunov-Powell source term in each cell consists of contributions from
Riemann solutions at the bordering interfaces and depends on the sign of the middle wave at each interface.
Thus, the Godunov-Powell source term is suitably upwinded. Note that assuming the normal magnetic field B1
to be constant for the whole domain leads to the source term being zero. This approach is novel and is very
different from the usual centered discretization of the Godunov-Powell source term ([32] and other references
therein). A related upwind discretization of a partial form of the source term s1 was presented in [9].
Remark 2.3. As stated earlier, when the normal magnetic field is constant i.e, B1L = B1R , the three-wave
solver presented here reduces to the HLLC three-wave solver of [20] which is shown to preserve positive pressure
and density. When the normal magnetic field is no longer constant, we are unable to provide a rigorous proof that
the three-wave solver preserves positive pressure and density. However, extensive numerical testing illustrates
that upwinding the source term leads to a scheme that preserves positive density and pressure.
2.3. HLL five-wave solver. The three-wave solver of the previous section does not model Alfvén waves
precisely, and instead diffuse these waves more than necessary. Alfvén waves can be approximated better by
extending the three-wave solver to an HLL type five-wave solver. In addition to the three waves with wave
speeds sL , sR and sM , we add two new waves with speeds s∗L and s∗R respectively with the requirement that
sL ≤ s∗L ≤ sM ≤ s∗R ≤ sR . Hence, solution of the Riemann problem at each interface is approximated by
four intermediate states WL ∗
, WL∗∗
, WR
∗∗
and WR∗
. A five-wave solver for ideal MHD equations with constant
normal magnetic fields was developed in [29] and we will extend this solver to the case of non-constant normal
magnetic fields below.
The outer-wave speeds sL and sR are determined by (2.6) (using the correction (2.15)) as in the three-wave
solver. We assume that the normal velocity is constant across the three inner waves i.e.,
(2.21) sM = u∗1L = u∗∗
1L = u1R = u1R ,
∗∗ ∗
as the inner waves model a contact discontinuity and Alfvén waves, and the normal velocity remains constant
across all three of them. Similarly, the waves with speeds s∗L and s∗R model Alfvén waves. Hence, as in [29], the
density and the “relaxed” pressures are constant across them leading to
θ = ρθ ,
ρ∗∗ = π1θ
∗ ∗∗ ∗
π1θ , θ ∈ {L, R}.
Furthermore, the wave with speed sM models a contact discontinuity and the tangential components of the
velocity and the magnetic field remain constant across it leading to
(2.22) σL = uσR = uσ ,
u∗∗ ∗∗ ∗∗ ∗∗
BσL = BσR
∗∗
= Bσ∗∗ ,
for σ ∈ {2, 3}. The normal magnetic field should only jump across the middle wave,i.e,
(2.23) ∗
B1θ = B1θ
∗∗
= B1θ , θ ∈ {L, R}.
As a result of (2.22) and (2.23), the Godunov-Powell source term takes the same form as in (2.9) with u∗σ , Bσ∗
being replaced by u∗∗
σ and Bσ for σ ∈ {2, 3}. We denote this source term by s
∗∗ 1,∗∗
.
Using local conservation (2.7) across the outermost waves and (2.21), we obtain unique values of ρ∗θ of the
form,
u1θ − sθ
(2.24) ρ∗θ = ρθ , θ ∈ {L, R}.
sM − sθ
Conservation across the entire wave fan leads to the following relation,
(2.25) FR − FL = sR WR − sL WL + (s∗R − sR )WR
∗
+ (sM − s∗R )W∗∗ + (s∗L − sM )WL
∗∗
+ (sL − s∗L )WL
∗
+ s1,∗∗ .
We can use the intermediate density states (2.24) and global conservation (2.25) to obtain the following expres-
sion for the middle speed,
π1R − π1L + ρR u1R (u1R − sR ) − ρL u1L (u1L − sL )
u∗1,L = u∗1,R = sM = .
ρR (u1R − sR ) − ρL (u1L − sL )
Note that this is the same expression as the middle speed in the three-wave solver. Similarly, one uses local
conservation (2.7) across the two outer waves to obtain the intermediate “relaxed” pressures,
∗
π1θ = π1θ + ρθ (u1θ − sθ )(u1θ − sM ),
for θ ∈ {L, R}. Note that π1L

∗
= π1R
∗
. Following [29], we can use conservation across the outer-most waves (2.7)
and obtain a system of two linear equations for u∗σθ and Bσθ
∗
for σ = {2, 3}. These equations can be explicitly
solved to obtain,
B1θ Bσθ (u1θ − sM ) ρθ (u1θ − sθ )2 − (B1θ )2
(2.26) u∗σθ = uσθ + , ∗
Bσθ = Bσθ .
ρθ (u1θ − sθ )(sM − sθ ) − (B1θ )2 ρθ (u1θ − sθ )(sM − sθ ) − (B1θ )2
Remark 2.4. The denominator in the above states (2.26) can become zero. This will typically occur in the
√
degenerate case that the Alfvén speed |B1 |/ ρ approaches the fast speed cf . Therefore, it is natural to switch
to the three-wave solver of the previous section when
1 1
1ρθ (u1θ − sθ )(sM − sL ) − (B1θ )2 1 < +ρθ c2f θ
for some small + > 0.
The intermediate energy states can be determined by local conservation (2.7) as
2
B1θ
Eθ (u1θ − sθ ) + π1θ u1θ − π1θ
∗
sM + 2 (u1θ
− sM ) + B1θ (B2θ u2θ + B3θ u3θ − B2θ
∗ ∗ ∗ ∗
u2θ − B3θ u3θ )
Eθ∗ = ,
sM − sθ
for θ ∈ {L, R}. The local conservation relations across the new waves imply that
|B1L | |B1R |
s∗L = sM − 2 ∗ , s∗R = sM + 2 ∗ .
ρL ρR
Hence the Alfvén wave speeds are accurately represented. Note that we use essentially the same Alfvén speeds as
in [29], but account for the variation of the normal magnetic field across the contact discontinuity. Furthermore,
conservation across the new waves reduces to
2 2
(2.27) Bσ∗∗ − BσL
∗
= sign (B1L ) ρ∗L (u∗∗σ − uσL ),
∗
Bσ∗∗ − BσR∗
= −sign (B1L ) ρ∗R (u∗∗σ − uσR ).
∗
These relations are identical to the exact Alfvén wave jump conditions. They imply that
2 2
sign (B1L ) ρ∗L u∗σL + sign (B1R ) ρ∗R u∗σR + BσR
∗
− BσL∗
uσ =
∗∗
2 ∗ 2 ∗
sign (B1L ) ρL + sign (B1R ) ρR
(2.28) 2 ∗ ∗ 2 2
sign (B1L ) ρR BσL + sign (B1R ) ρ∗L BσR ∗
+ ρ∗L ρ∗R (u∗σR − u∗σL )
Bσ =
∗∗
2 2
sign (B1L ) ρ∗R + sign (B1R ) ρ∗L
Remark 2.5. Observe that if
(2.29) sign (B1R ) '= sign (B1L ) and ρ∗R = ρ∗L ,
the relations (2.27) can not be solved, and the formulas (2.28) break down. In this case we relax to the HLL
three-wave solver of the previous section. This should be seen in light of the discontinuity in the jump conditions
(2.27) in the non strictly hyperbolic case B1 = 0. Furthermore, if B1L = 0 and B1R '= 0 or vice versa, we get
meaningful formulas, but we need to check that the jump conditions across the middle wave hold. They become
∗
SM WL − F∗L = SM WR
∗∗
R +s
− F∗∗ 1,∗∗
,
which is easily verified. When B1L = B1R = 0, it simply means that s∗L = sM = s∗R , hence we do not need to
calculate the ∗∗-states. The jump conditions are again easily verified, if we assume that the source s1,∗∗ is zero
(which is the limiting value away from (2.29)). Hence, the numerical fluxes and sources can be calculated in
the same manner as for non-zero B1 .
Finally, the remaining energy states are given by conservation across the Alfvén waves resulting in,
B1θ (B2θ
∗ ∗
u2θ + B3θ
∗ ∗
2 − B3 u3 )
u3θ − B2∗∗ u∗∗ ∗∗ ∗∗
Eθ∗∗ = Eθ∗ +
s∗θ − sM
for θ ∈ {L, R}. This completes a description of the states of the five-wave solver. The corresponding fluxes
can be determined by local conservation and the numerical flux is obtained similar to the formula (2.16). We
remark that whenever B1L = B1R , the above solver reduces to the five-wave solver derived in [29].
For the discretization of the corresponding Godunov-Powell source term in this case, we use exactly the same
arguments as in the case of the three-wave solver (as B1 jumps only across the middle wave where the velocity
and the tangential magnetic fields are constant) to obtain that
S1,n 1,∗∗ 1,∗∗

i,j = si−1/2,j 1{sM,i−1/2,j ≥0} + si+1/2,j 1{sM,i+1/2,j <0} ,
where s1,∗∗
i±1/2,j is defined as in (2.9) with the ∗∗ replacing the ∗ states.
2.4. Fluxes and sources in the y-direction. We have completed a description of the numerical flux F and
discretized source S1 in the x-direction. In-order to complete the scheme (2.3), we need to specify the numerical
flux in the y-direction G and the corresponding Godunov-Powell source term S2 . This is straightforward as the
form of equations in each direction is similar.
The numerical flux G is defined in terms of both a three-wave solver and a five-wave solver. The three-wave
solver is analogous to the states and fluxes obtained in (2.5) with normal velocity u2 and normal magnetic field
B2 . Similarly, the discretized source term S2 is similar to S1 defined in (2.20). The only change is to replace
the normal velocity and magnetic fields to u2 and B2 respectively. The five-wave solver is analogously defined.
The specification of G, S2 completes the description of the scheme (2.3).
2.5. Second-order accurate schemes. The finite volume scheme (2.3) is first-order accurate in both space
and time. For practical applications, we need higher order of accuracy. We will design a finite volume scheme
based on (2.3) that is second-order accurate in both space and time. The semi-discrete form of this scheme is
given by,
d 1 1 /2 ,
/1 + S
(2.30) Wi,j = F i,j = − (Fi+1/2,j − Fi−1/2,j ) − (Gi,j+1/2 − Gi,j−1/2 ) + S
dt ∆x ∆y i,j i,j
where Wi,j (t) is the cell-average of the unknown at time t. We will define the numerical fluxes F, G and the
sources S/1, S
/ 2 below.
It is standard [24] to replace the piecewise constant approximation Wi,j with non-oscillatory piecewise
linear reconstructions in-order to obtain second-order spatial accuracy. There are a variety of reconstructions
including the popular TVD-MUSCL limiters [43]. We will use second-order ENO reconstruction [22] and WENO
reconstruction [36] as these procedures can be easily extended to obtain even higher-order schemes.
2.5.1. ENO Reconstruction: Given the cell averages Wi,j , we reconstruct in the primitive variables Vi,j =
{ρi,j , ui,j , Bi,j , pi,j }. Define the ENO-differences in each direction as
* *
Vi+1,j − Vi,j if Γxi,j ≤ 1, Vi,j+1 − Vi,j if Γyi,j ≤ 1,
D Vi,j =
x
D Vi,j =
y
Vi,j − Vi−1,j otherwise. Vi,j − Vi,j−1 otherwise.
where
|ψ(Vi+1,j ) − ψ(Vi,j )| |ψ(Vi,j+1 ) − ψ(Vi,j )|
Γxi,j = , Γyi,j = ,
|ψ(Vi,j ) − ψ(Vi−1,j )| |ψ(Vi,j ) − ψ(Vi,j−1 )|
and ψ for some function ψ called the global smoothness indicator. We use ψ(V) = ρ + B2 .
Note that for piecewise linear reconstruction, the ENO procedure reduces to providing a limiter for the slopes
in each direction. The reconstructed piecewise linear function is each cell Ii,j is denoted by
1 x 1 y
Vi,j (x, y) = Vi,j + D Vi,j (x − xi ) + D Vi,j (y − yj ).
∆x ∆y
The reconstructed conservative variables can be easily obtained by transforming the reconstructed primitive
variables.
2.5.2. WENO procedure: As an alternative to the above reconstruction, consider the following cell-gradients
3 4 5 6
D̄x Vi,j = ωi,j
x
(Vi+1,j − Vi,j ) + 1 − ωi,j
x
(Vi,j − Vi−1,j ) ,
3 4 6
y y 5
D̄y Vi,j = ωi,j (Vi,j+1 − Vi,j ) + 1 − ωi,j (Vi,j − Vi,j−1 ) ,
where the weights are given by,
a0i,j 1 2
x
ωi,j = , a0i,j = x,0 , a1i,j = x,1 ,
a0i,j + a1i,j 3(+ + βi,j ) 3(+ + βi,j )
y b0i,j 1 2
ωi,j = , b0i,j = y,0 , b1i,j = y,1 ,
b0i,j + b1i,j 3(+ + βi,j ) 3(+ + βi,j )
where + is a small positive number, and the parameters are given by
x,0
βi,j = (ψ(Vi+1,j ) − ψ(Vi,j ))2 , β x,1 = (ψ(Vi,j ) − ψ(Vi−1,j ))2 ,
y,0
βi,j = (ψ(Vi,j+1 ) − ψ(Vi,j ))2 , β y,1 = (ψ(Vi,j ) − ψ(Vi,j−1 ))2 ,
and the indicator function ψ is defined above. The corresponding linear reconstruction is given by,
1 x 1 y
(2.31) Vi,j (x, y) = Vi,j + D Vi,j (x − xi ) + D Vi,j (y − yj ).
∆x ∆y
Note that the choice of weights implies that the WENO approximation (2.31) is third-order accurate for smooth
solutions (at least in one-space dimension).
ENO and WENO reconstructions have been used extensively to compute approximate solutions of the ideal
MHD equations (1.4), see [41, 5] and other references therein. However, both the ENO and WENO recon-
structions suffer from a common problem: the reconstructed densities and pressures may not be positive. Since
the positivity of density and pressure is absolutely essential for obtaining any physically meaningful results,
the reconstructions have to be modified further. A simple modification consists of further limiting the slope in
either direction. Let Dx ρ and Dx p be the ENO-gradients of density and pressure in the x-direction, we follow
[30] and introduce the following Perthame-Shu clipping,
(2.32) Dx,c ρi,j = max {−ωρi,j , min {ω, Dx ρi,j }} , Dx,c pi,j = max {−ωpi,j , min {ω, Dx pi,j }} ,
where ω < 2 is a positive parameter. Simple calculations show that the above modification ensures that the
reconstructed pressure and density remain positive. We choose ω = 1.9 in our simulations. Similarly in the
y-direction, we have
Dy,c ρi,j = max {−ωρi,j , min {ω, Dy ρi,j }} , Dy,c pi,j = max {−ωpi,j , min {ω, Dy pi,j }} .
However, this only guarantees positivity of the reconstructed variables, not for the updated ones.
A further modification of the gradients was suggested in a recent paper [44] to ensure positivity preserving
updated states. Denote
(Dx,c Vi,j , Dy,c Vi,j ) = ({Dx,c ρi,j , Dx ui,j , Dx Bi,j , Dx,c pi,j } , {Dy,c ρi,j , Dy ui,j , Dy Bi,j , Dy,c pi,j }) ,
and
" " ##
1 2 2 1 1 2 2
=
Lxi,j ρi,j |D ui,j | + |D Bi,j | + (min {0, ρi,j (ui,j · D ui,j )})
x x x
+ (Dx,c ρi,j ) |Dx ui,j | ,
8 2 2ρi,j
" " ##
y 1 2 2 1 1 2 2
Li,j = ρi,j |D ui,j | + |D Bi,j | + (min {0, ρi,j (ui,j · D ui,j )})
y y y
+ (Dy,c ρi,j ) |Dy ui,j | ,
8 2 2ρi,j
pi,j
Ri,j = .
γ−1
Then, we further modify the gradients by
0 0
7 x Vi,j = Dx,c Vi,j R i,j 7 y Vi,j = Dy,c Vi,j Ri,j
(2.33) D , D .
max{Lxi,j , Ri,j } max{Lyi,j , Ri,j }
Thus D7 x, D
7 y are the gradients that we use for the final reconstruction. In [44] it was shown that this results in
a positive scheme provided that the underlying first-order scheme is positive, and a particular discretisation of
the source term is used.
The reconstructed states are now given by
(2.34) V 7 x Vi,j (x − xi ) + 1 D
7 i,j (x, y) = Vi,j + 1 D 7 y Vi,j (y − yj ).
∆x ∆y
A similar procedure can be used for the WENO-gradients to obtain a modified WENO reconstruction.
8 i,j (x, y). De-
The reconstructed primitive variables correspond to the reconstructed conservative function W
fine the point-values,
Wi,j
E 8 i,j (xi+1/2 , yj ),
=W Wi,j
W 8 i,j (xi−1/2 , yj ),
=W
Wi,j
N 8 i,j (xi , yj+1/2 ),
=W Wi,j
S 8 i,j (xi , yj+1/2 ).
=W
We can use the above defined values to define the second-order numerical fluxes as
4 E 5 4 N 5
Fi+1/2,j = F Wi,j , Wi+1,j
W
, Gi,j+1/2 = G Wi,j , Wi,j+1
S
,
where F,and G are given by either the three-wave solver or the five-wave solver of the previous section. Similarly,
the second-order source terms can be calculated as
S1i,j = s1,∗ 1,∗
i−1/2,j 1{sM,i−1/2,j ≥0} + si+1/2,j 1{sM,i+1/2,j <0} ,
where s1,∗
i+1/2,j is defined as in (2.9), but with the values Wi,j , Wi+1,j replaced by Wi,j , Wi+1,j . The source
E W
S2i,j in the y-direction is defined analogously. Observe that for smooth solutions, the discretized source S1i,j
vanishes to truncation order with (B1E )i,j − (B1W )i+1,j . Hence, we need to add an extra term for second-order
consistency. However, this term should vanish when S1i,j becomes significant at jumps (see e.g., [1] for an
analogous situation). We suggest the following simple modification,
 
0
/ 1 = S1 +  Bi,j  1 D
 
S 7x 1
i,j i,j ui,j · Bi,j  ∆x Bi,j .
ui,j
The term S / 2 in the y-direction is analogously defined. This way of discretizing the source was found to be very
i,j
stable in [44], while a similar but more complicated form also gave a provably positive scheme. Note that S / 1,2 ,
i,j
are consistent second-order discretizations of the Godunov-Powell source terms s1,2 . Hence, we have completed
a description of the (formally) second-order accurate in space semi-discrete scheme (2.30).
2.6. Time Stepping: The standard scheme for a first order approximation in time is the forward-Euler time
stepping, formally written
n+1
Wi,j = Wi,j
n
+ ∆tn F ni,j
where F ni is defined in (2.30). For overall second-order schemes, we use the second-order strong-stability
preserving Runge-Kutta (SSP) time stepping [21],
∗
Wi,j = Wi,j
n
+ ∆tn F ni,j ,
∗∗
Wi,j = Wi,j
∗
+ ∆tn F ∗i,j ,
1
n+1
Wi,j = (Wi,jn
+ Wi,j
∗∗
).
2
The time step is determined by a standard CFL condition. This completes our description of the finite volume
schemes for (2.1).
3. Numerical Experiments
We will validate the first- and second-order finite volume schemes on a series of numerical experiments in
both one- and two-space dimensions. We test a total of six schemes:
H3 First order with the HLL three-wave solver,
H5 first order with the HLL five-wave solver,
H3 E second order with HLL three-wave solver and ENO reconstruction,
H3 W second order with HLL three-wave solver and WENO reconstruction,
H5 E second order with HLL five-wave solver and ENO reconstruction,
H5 W second order with HLL five-wave solver and WENO reconstruction.
All the second order schemes use the positivity preserving modifications (2.34). The first order schemes are
evolved with a CFL number of 0.45 (which is theoretically sound due to excluding wave interactions in the
cells), and the second order schemes use a CFL number of 0.9. In all our computations we use γ = 5/3.
Regarding the measurement of errors, if we have a reference solution available, then we define the relative
error as
,α − αref ,
100 × ,
,αref ,
where α is (a component of) the numerical approximation and αref is (the same component of) the reference
solution, and , · , is some (usually L1 ) norm.
3.0.1. Brio-Wu shock tube. This is a standard one-dimensional numerical test case for ideal MHD [11]. The
initial conditions are given by
*
(1.0, +1.0, 0, 0, 0, 0.7, 0, 1.0), if x < 0.5,
(ρ, ρu, p, B) =
(0.3, −0.3, 0, 0, 0, 0.7, 0, 0.2), otherwise.
The computational domain is (x, t) ∈ [0, 1.5] × [0, 0.5] with Neumann boundary conditions. Note that the
normal magnetic field is constant. Therefore, the H3 scheme and H5 scheme reduce to the solvers presented
in [20] and [29] respectively. However, the four higher-order schemes are different from those presented in the
literature. We present the computed total energy and the magnetic field B2 with the schemes at time t = 0.5
in figure 1. The reference solution in this case is calculated with the H5 W second-order scheme on a mesh with
Figure 1. Results for the Brio-Wu shock tube with 200 grid points at t = 0.5. Reference solu-
tion is the H5 W scheme with 3200 grid points. Left: Energy Right: Magnetic field component
B2
3200 points. The solution is quite complicated containing shock waves, contact discontinuities and rarefaction
waves. As expected, the H5 scheme is more accurate than the H3 schemes. The second-order schemes are
clearly more accurate than the first-order schemes. The differences are also illustrated in Table 1, showing the
relative percentage errors in the L1 norm of the total energy on a sequence of meshes. The table confirms the
observations obtained from the figure. Both the formal first-order accurate schemes have an average convergence
rate close to 0.7 . The second-order schemes are more accurate than the first-order schemes with considerably
smaller errors. The H5 E and H5 W are slightly more accurate than the corresponding H3 E and H3 W scheme.
The second-order schemes have an average rate of convergence around 1. This is to be expected as the solution
contains discontinuities, and the order of accuracy deteriorates near these.
M H3 rate H3 E rate H3 W rate H5 rate H5 E rate H5 W rate

50 2.2e-00 1.3e-00 1.3e-00 1.7e-00 1.2e-00 1.2e-00
100 1.5e-00 0.6 7.9e-01 0.7 7.2e-01 0.9 1.2e-00 0.5 7.1e-01 0.8 6.5e-01 0.9
200 9.8e-01 0.6 5.0e-01 0.7 4.6e-01 0.6 7.4e-01 0.7 4.4e-01 0.7 4.1e-01 0.7
400 5.3e-01 0.9 1.9e-01 1.4 1.7e-01 1.4 3.7e-01 1.0 1.6e-01 1.5 1.4e-01 1.6
800 3.4e-01 0.6 1.1e-01 0.8 9.7e-02 0.8 2.2e-01 0.8 9.6e-02 0.7 8.3e-02 0.8
1600 1.9e-01 0.8 5.4e-02 1.0 4.4e-02 1.1 1.2e-01 0.9 4.8e-02 1.0 3.9e-02 1.1
Table 1. Relative percentage errors for the total energy in L1 at time t = 0.5 and the order
of convergence for the Brio-Wu shock tube for various mesh sizes M using the H5 W scheme
with 3200 grid points as a reference solution.
3.0.2. Super-fast expansion. A crucial stability criteria is the ability of a numerical scheme to preserve positive
density and pressure. This numerical experiment [8, 18] is configured to highlight positivity aspects of schemes.
We consider the one-dimensional form of the MHD equations (1.4) in the computational domain [0, 1] with
initial data, *
(1.0, −3.1, 0, 0, 0, 0.5, 0, , 1.0), if x < 0.5,
(ρ, ρu, p, B) =
(1.0, 3.1, 0, 0, 0, 0.5, 0, 1.0), otherwise.
The exact solution involves a super-fast expansion with one left moving rarefaction and one right moving
rarefaction wave separated by a region of very low density and pressure. Linearized solvers like the Roe solver
crash almost instantly on account of negative pressure [18].
Since the normal magnetic field is constant, the first-order H3 and H5 schemes reduce to the standard HLLC
solver of [20] and HLLD solver of [29], respectively. These solvers are shown to be positivity preserving. The
second-order H3 E, H3 W , H5 E and H5 W solvers are designed to preserve positive density and pressure. As we
are unable to provide a rigorous proof of this fact, we test all the four second order schemes for the super-fast
expansion and find them to be positivity preserving for all tested resolutions. Since there are no differences in
results between the three- and five-wave solvers in this test case, we present the density, computed with H3 ,
H3 E and H3 W at time t = 0.1 for 200 mesh points in figure 2. The results show that both the first and second
order solvers preserve positivity. The second-order solvers are slightly more accurate than the first-order H3
solver. Similar results were obtained for finer resolutions.
We remark that the positive preserving modifications (2.33) are absolutely essential to obtain stability. If we
just use the Perthame-Shu clipping (2.32), the second-order schemes crash for very small times on account of
negative pressures. This test further reinforces the case for modifying gradients in second-order reconstruction,
as presented in [44].
3.0.3. Godunov-Powell Magnetic Advection. This test case is constructed to demonstrate the performance of
the schemes in one-dimension with a non-constant normal magnetic field. Consider the one-dimensional semi-
conservative form of the MHD equations (1.2) with the following initial data,
(ρ, u1 , u2 , u3 , B1 , B2 , B3 , p) = (1, 1, 0, 0, 1 − c sin(2πx), 0.5, 0, 0, 0.5),
where c is a constant. A straightforward calculation shows that the exact solution is given by,
(ρ, u1 , u2 , u3 , B1 , B2 , B3 , p) = (1, 1, 0, 0, 1 − c sin(2π(x − t)), 0.5, 0, 0, 0.5),
Figure 2. The pressure for super-fast expansion test with H3 ,H3 E and H3 W solvers for 200
mesh points.
Thus, the solution consists of advection of the normal magnetic field due to the presence of the Godunov-Powell
source term. Note that putting c = 0 gives a constant solution. The aim is to take non-zero values of c and test
the schemes. It turns out that any standard scheme like the three-wave solver of [20] or the five-wave solver of
[29] crashes almost immediately if c > 0.001 due to unphysical state values (B1 was evaluated as a local average
in the numerical fluxes). This is not unexpected since these schemes are based on the form (1.4) which needs
the normal magnetic field to be a constant in space. Hence, this test case is vital in showing the effectiveness
of our new schemes based on the semi-conservative form.
We consider this problem on a domain [0, 1] and c = 1 with periodic boundary conditions. Figure 3 shows the
approximate normal magnetic field B1 at time t = 1. The H3 and H5 schemes approximate the solution like any
standard first-order scheme for linear advection. The second-order schemes resolve the solution much better,
and there is very little visual difference between the exact solution and the second-order schemes at this mesh
resolution. Further study of the approximation is reported in the Table 2, which shows the relative percentage
Figure 3. B1 distribution for Powell magnetic advection test with 200 grid points.
errors in L1 for the magnetic field B1 on a sequence of meshes. There is no difference between the H3 and H5
solvers. Both the first-order schemes have the expected rate of convergence 1. Among the second-order schemes,
the WENO based H3 W and H5 W are more accurate than their ENO based counterparts. In fact, the observed
order of convergence for the ENO schemes is 2, but the WENO schemes approach an order of convergence close
to 3. This is on account of the design of a WENO scheme where third-order is achieved for piecewise linear
reconstruction. Note that the positivity-preserving modifications have not reduced the orders of accuracy of
either the ENO or the WENO schemes. Thus, this experiment demonstrates that our schemes approximate
the semi-conservative form (1.2) in one-dimension quite well, particularly for problems with perturbations from
constant normal magnetic fields (standard schemes fail in these cases, even for very small perturbations of
constant normal magnetic fields). It also illustrates that the semi-conservative Godunov-Powell form (1.2) is
stable with respect to perturbations in the constant normal magnetic field whereas the standard form (1.4) is
not even well-defined.

50 9.4e-00 3.2e-00 2.2e-00 9.4e-00 3.2e-00 2.2e-00
100 5.1e-00 0.9 9.9e-01 1.7 6.0e-01 1.9 5.2e-00 0.9 9.9e-01 1.7 6.0e-01 1.9
200 2.8e-00 0.9 2.9e-01 1.8 1.5e-01 2.0 2.8e-00 0.9 2.9e-01 1.8 1.5e-01 2.0
400 1.5e-00 0.9 8.7e-02 1.7 3.5e-02 2.1 1.5e-00 0.9 8.7e-02 1.7 3.5e-02 2.1
800 8.1e-01 0.9 2.4e-02 1.9 7.3e-03 2.3 8.1e-01 0.9 2.4e-02 1.9 7.3e-03 2.3
1600 4.1e-01 1.0 6.4e-03 1.9 1.2e-03 2.6 4.1e-01 1.0 6.4e-03 1.9 1.2e-03 2.6
3200 2.1e-01 1.0 1.7e-03 1.9 1.5e-04 3.0 2.1e-01 1.0 1.7e-03 1.9 1.5e-04 3.0
Table 2. Relative percentage errors for B1 at time t = 1 and the order for the Powell magnetic
advection test for various mesh sizes M .
3.0.4. Rotor Problem: We start considering two-dimensional numerical experiments with this standard example
(introduced in [2], considered in [42] among others). The computational domain is (x, t) ∈ [0, 1]2 × [0, 0.295]
with Neumann boundary conditions. The initial data are given by

10.0
 if r < 0.1,
ρ = 1 + 9f (r) if 0.1 ≤ r < 0.115,


1.0 otherwise,
with r(x) = |x − (0.5, 0.5)| and
23 − 200r
f (r) = .
3
The other variables are initially

(−(10y − 5)ρ, (10x − 5)ρ) if r < 0.1,
4 1 5 
ρu , ρu = (−(10y − 5)f (r)ρ, (10x − 5)f (r)ρ)
2
if 0.1 ≤ r < 0.115,


(0.0, 0.0) otherwise,
4 3 1 2 3 5 4 √ 5
ρu , B , B , B , p = 0.0, 2.5/ π, 0.0, 0.0, 0.5 .
This describes a dense rotating region surrounded by static plasma with a uniform magnetic field. The pressure
drops to very low values in the center. The main difficulty in the numerical solution of this problem is the low
pressure, particularly on fine meshes. As stated in the introduction, most results presented in the literature
show the approximation obtained on relatively coarse meshes. On coarse meshes, the numerical dissipation
is large and provides some stability. However, computations on fine meshes lead to crashes due to negative
pressures (see [18] for illustrations of this). We compute with all the six schemes and show the computed
pressures at time t = 0.295 and on a mesh with 200 × 200 mesh points in Figure 4. The figure shows that both
the first-order schemes provide a stable but diffusive approximation and the H5 scheme is more accurate than
the H3 scheme. The second-order schemes are much more accurate and capture the shocks and smooth regions
sharply at this resolution. For a more elaborate quantitative study of this problem, we compute solutions on a
very fine 4000 × 4000 mesh and find that all the six schemes are stable and approximate the solution very well.
We show the computed pressure with the H5 W scheme on a 4000 × 4000 mesh in Figure 5 and use it as the
reference solution.
We tabulate the relative percentage errors in L1 for the pressure with respect to this reference solution and
present them in table 3. The table shows that the first-order H5 scheme is (about thirty percent) more accurate
than the first-order H3 scheme. However, the second-order schemes are much more accurate than the first-order
(a) H3 (b) H5
(c) H3 E (d) H5 E
(e) H3 W (f) H5 W
Figure 4. Pressure for the rotor problem on a 200 × 200 mesh at time t = 0.295 all scaled to
the extrema of the pressure obtained for a 200 × 200 mesh.
Figure 5. This figure shows the pressure for the rotor problem on a 4000 × 4000 mesh at time
t = 0.295 using the H5 W scheme.
schemes. Sometimes, the gain in accuracy is an order of magnitude by using a second-order scheme. There is a
gain in accuracy using the H5 solver together with a second-order scheme. Similarly, the WENO-based schemes
H3 W and H5 W are more accurate than their ENO-based counterparts. The observed rate of convergence for
the first-order schemes is around 0.5 and that of the second-order schemes is better than 1 (except at the lowest
resolutions).

50 10.0e-00 5.3e-00 5.1e-00 9.2e-00 4.6e-00 4.4e-00
100 7.7e-00 0.4 3.1e-00 0.8 2.9e-00 0.8 6.6e-00 0.5 2.5e-00 0.9 2.3e-00 0.9
200 5.3e-00 0.5 1.6e-00 1.0 1.4e-00 1.1 4.5e-00 0.6 1.3e-00 0.9 1.1e-00 1.1
400 3.5e-00 0.6 7.6e-01 1.1 6.3e-01 1.2 2.8e-00 0.7 6.0e-01 1.1 5.1e-01 1.1
800 2.2e-00 0.7 3.5e-01 1.1 2.8e-01 1.2 1.7e-00 0.7 2.7e-01 1.2 2.0e-01 1.4
Table 3. Relative percentage errors in L1 for the pressure at time t = 0.295 for the rotor
problem for various mesh sizes M using the H5 W scheme on a 1600 × 1600 mesh as a reference
solution.
Another key issue in numerical simulations of multi-dimensional MHD is the behavior of divB. Note that
the initial magnetic field is divergence free, hence the solutions of (2.1) are expected to remain divergence free.
However, the schemes we use do not preserve any discrete divergence. We consider the standard second order
discrete divergence,
(B1 )i+1,j − (B1 )i−1,j (B2 )i,j+1 − (B2 )i,j−1
(3.1) div(B)i,j = +
2∆x 2∆y
and present the L1 norm of the above discrete divergence in Table 4. As expected, all the six schemes produce
a nonzero discrete divergence. However, the values are quite small, and seem to be decreasing with increasing
mesh size, albeit quite slowly. Note that the divergence values are higher with the second-order schemes than
with the first-order schemes. This is not unexpected as the second-order schemes resolve the shocks within fewer
mesh points and hence generate a larger discrete divergence. These data primarily point towards the difficulty
of numerically evaluating derivatives at under-resolved flow features. The key point is that divergence errors
are not effecting the stability of the schemes as all the schemes are stable even at the finest mesh resolution of
4000×4000 mesh. This is unusual, and we have not come across other papers presenting solutions on comparable
mesh sizes. As an example, the reference solutions in [42] were computed on a 400 × 400 mesh. We would like to
point out that the positivity preserving modifications (2.32) and (2.33) are absolutely essential for stability of
the second-order schemes. The computations crash almost immediately if these modifications are not employed.
Similarly the upwinding of the source term in (1.2) is absolutely essential as the central discretization of this
term leads to negative pressures almost immediately, even for the first-order schemes for fine mesh resolutions.
M H3 H3 E H3 W H5 H5 E H5 W
50 2.2e-01 2.5e-01 2.1e-01 2.2e-01 2.4e-01 2.2e-01
100 1.7e-01 1.9e-01 1.8e-01 1.7e-01 2.0e-01 1.9e-01
200 1.4e-01 1.8e-01 1.6e-01 1.5e-01 1.9e-01 1.6e-01
400 1.3e-01 1.8e-01 1.5e-01 1.2e-01 1.6e-01 1.4e-01
800 1.1e-01 1.7e-01 1.5e-01 9.6e-02 1.5e-01 1.2e-01
1600 1.0e-01 1.6e-01 1.4e-01 8.2e-02 1.4e-01 1.2e-01
Table 4. The L1 -norm of the discrete divergence at time t = 0.295 for the rotor problem for
various mesh sizes M .
3.0.5. Orszag-Tang vortex. This commonly used benchmark test [42] has initial conditions given by
4 5
(3.2) (ρ, ρu, B, p) = γ 2 , −γ 2 sin(πy), γ 2 sin(πx), 0, − sin(πy), sin(2πx), 0, γ .
The computational domain is (x, t) ∈ [0, 2]2 × [0, 1] with periodic boundary conditions.
Even though the initial data are smooth, the solution develops shocks near the diagonals and a current
sheet in the center of the domain. The solution also has interesting smooth features. We compute with all the
six schemes and present the computed pressure at the final time on a 200 × 200 mesh in Figure 6. Both the
first-order H3 and H5 schemes are stable but dissipative. The shocks are smeared and the central vortex is
not resolved. The H5 scheme is better at approximating the solution than the H3 scheme. The second-order
schemes resolve the solution far better. The resolution of the shocks with the second-order schemes is very
impressive. The smooth regions are also resolved quite accurately. We compute the a reference solution using
the H5 W scheme on a 4000 × 4000 mesh. Figure 7 shows the pressure of the reference solution. We observe that
the H5 W scheme at this very fine mesh is stable and resolves the shocks as well as the central current sheet very
well. Table 5 shows relative percentage L1 errors in pressure on a sequence of meshes. As shown in the table,
the first-order H5 scheme has lower errors than the first-order H3 scheme. Similarly, the second-order schemes
significantly outperform (by an order of magnitude) the first-order schemes. In particular, the WENO-based
schemes have lower errors than the ENO-based ones. Note that the second-order H5 W scheme is the most
accurate with respect to the errors and has the best rate of convergence of about 1.3. The results are consistent
with those obtained for the rotor problem.
M H3 rate H3 E rate H3 W E rate H5 rate H5 E rate H5 W rate

50 3.6e+01 1.5e+01 1.4e+01 2.9e+01 1.2e+01 1.1e+01
100 2.9e+01 0.3 8.8e-00 0.8 8.0e-00 0.8 2.3e+01 0.3 7.1e-00 0.8 6.2e-00 0.8
200 2.2e+01 0.4 5.0e-00 0.8 4.4e-00 0.9 1.6e+01 0.5 4.0e-00 0.8 3.4e-00 0.9
400 1.5e+01 0.6 2.6e-00 0.9 2.4e-00 0.9 1.1e+01 0.5 2.0e-00 1.0 1.6e-00 1.1
800 1.0e+01 0.6 1.5e-00 0.8 1.3e-00 0.9 7.2e-00 0.6 8.3e-01 1.3 5.9e-01 1.4
Table 5. Percentage relative errors in L1 for pressure at time t = π for the Orszag-Tang vortex
for various mesh sizes M using the H5 W scheme on a 1600 × 1600 mesh as a reference solution.
As in the rotor problem, we tabulate the L1 discrete divergence (exact divergence is zero as we have divergence
free initial data) in table 6. The divergence values are low again and decrease very slowly (or remain constant)
as the mesh is refined. Again, the discrete divergence values don’t seem to affect the scheme’s performance.
(a) H3 (b) H5
(c) H3 E (d) H5 E
(e) H3 W (f) H5 W
Figure 6. Pressure for the Orszag-Tang vortex on a 200 × 200 mesh at time t = 1 scaled to
the extrema of the pressure in the reference solution.
Figure 7. This figure shows the computed pressure for the Orszag-Tang vortex using the H5 W
scheme on a 4000 × 4000 mesh at time t = π.
M H3 H3 E H3 W H5 H5 E H5 W
50 7.7e-02 2.5e-01 2.2e-01 1.5e-01 1.9e-01 1.7e-01
100 9.2e-02 2.6e-01 2.4e-01 1.6e-01 1.6e-01 1.5e-01
200 9.9e-02 2.7e-01 2.4e-01 1.4e-01 1.7e-01 1.5e-01
400 9.3e-02 2.0e-01 2.0e-01 1.3e-01 1.6e-01 1.3e-01
800 8.7e-02 1.8e-01 1.5e-01 1.1e-01 1.2e-01 1.0e-01
1600 8.4e-02 1.5e-01 1.2e-01 1.0e-01 1.0e-01 9.2e-02
Table 6. Divergence in L1 at time t = π for the Orszag-Tang vortex for various mesh sizes M .
3.0.6. Cloud-Shock Interaction. This is a benchmark test describing the interaction of a dense region (cloud) at
rest with a moving shock. The computational domain is (x, t) ∈ [0, 1]2 × [0, 0.06] with artificial Neumann type
boundary conditions. The initial conditions consist of a shock moving to the right initially located at x = 0.05,
and a circular cloud of density ρ = 10 and radius r = 0.15 centered at x = (0.25, 0.5).

3.86859 if x < 0.05,

ρ = 10.0 if |x − (0.25, 0.5)| < 0.15,


1.0 otherwise,
*
(11.2536, 0, 0) if x < 0.05,
u=
(1.0, 0, 0) otherwise,
*
167.345 if x < 0.05,
p=
1.0 otherwise,
*
(0, 2.18261820, −2.18261820) if x < 0.05,
B=
(0, 0.56418958, 0.56418958) otherwise.
The cloud is initially in hydrostatic equilibrium with the surrounding fluid. The bubble stays stationary whereas
the shock travels towards it, hits it and starts interacting with it. This interaction generates a bow shock in
the front, tail shocks in the rear and we expect the creation of interesting turbulent-like structures where the
cloud interacts with the shock. Again, we compute the solutions with all the six schemes and present numerical
results of the total energy on a very fine mesh 1600 × 1600 mesh with all the six schemes in Figure 8. The
figure illustrates the stability of all the six schemes at this fine mesh resolution. The first-order schemes are
a bit dissipative and do not approximate the turbulent like structures that are visible in the second order
computations. The second order schemes are much more accurate with good resolution of the shocks as well
as the turbulent-like structures in the bubble. The differences between the first-order and the second-order
schemes are considerable while there are some minor differences between the ENO and WENO approximations.
In this test case we found the discrete divergence to be larger than in the other test cases, with L1 norms in
the range 0.1–4.0. Furthermore we did not observe a strong convergence to zero as the mesh was refined. This
might be due to the strong shocks present in this test case. However, there is no sign that this affects either
the stability or the resolution of the schemes.
3.0.7. Isothermal blast wave. This test case is taken from [3]. It easily exhibits spurious behavior if the equations
are not properly discretised. The equation of state is isothermal in this case, i.e, p = ρ, the computational domain
is [0, 1] × [0, 1] and the initial data are
5
{u, B1 , B2 , B3 }t=0 = {0, √ , 0, 0}
π
and *
100 if |x − 0.5|2 + |y − 0.5|2 ≤ (0.05)2 ,
ρt=0 =
0. otherwise
The solution takes the form of a blast wave spreading out from the high density cloud. We test our upwind
versions of the Godunov-Powell source term to ascertain whether upwinding the source term imparts numerical
stability in this test case. For the sake of clarity, we present the results only with the H3 and H3 W schemes.
Similar results were obtained with the H5 schemes and with the ENO reconstruction.
In order to illustrate the necessity of the upwinded Godunov-Powell source term, we also considered the
central source discretisation of ([32]) , and schemes with the Godunov-Powell source set to zero. This lead to
the following schemes:
N P H3 First order HLL three-wave solver with Godunov-Powell source term in (2.3) set to zero.
N P H3 W Second-order WENO version of N P H3 .
CP H3 First order HLL three-wave solver, central discretization of the Godunov-Powell source term.
CP H3 W Second-order WENO version of CP H3 .
Note that both N P H3 and CP H3 use the same numerical flux as H3 and the only difference is in the source
term.
We display the computed density (on a logarithmic scale) at time t = 0.09 with the above four schemes and
H3 , H3 W on a uniform mesh of 400 × 400 mesh points in figure 9. The results presented in figure 9 reveal that
the first-order version of the HLL three-wave solver works quite well, even when the Godunov-Powell source
term is set to zero. This might be on account of the numerical flux, which incorporates some information about
the discontinuities in the normal magnetic field. However, the second-order version of the scheme with zero
Godunov-Powell source terms led to instabilities as shown in the top right of figure 9. The central discretization
of the Godunov-Powell source term in (2.3) was unstable even at first-order. The instabilities are manifested
as oscillations in the middle left of figure 9. The second-order version of this scheme is even more unstable and
features large amplitude oscillations as shown in the middle right of figure 9 (Note that scales in each of the sub-
figures are different). These oscillations continue to grow as the mesh is refined leading to a crash of the central
Powell source based scheme on a finer 800 × 800 mesh. On the other hand, both the H3 and H3 W schemes
show no evidence of the formation of any instabilities (check the bottom row of figure 9) at this mesh resolution
or at even finer mesh resolutions. Furthermore, the second-order scheme approximates the shock waves quite
sharply. These results indicate that not only should the Godunov-Powell source term be introduced, it should
also be discretized in a careful upwind manner in-order to obtain numerical stability. The schemes presented in
this paper are tailor made for these needs.
(a) H3 (b) H5
(c) H3 E (d) H5 E
(e) H3 W (f) H5 W
Figure 8. Energy distribution for the cloud-shock interaction on a 1600 × 1600 mesh at time t = 0.06.
(a) N P H3 (b) N P H3 W
(c) CP H3 (d) CP H3 W
(e) H3 (f) H3 W
Figure 9. Density (on a logarithmic scale) for the Isothermal blast wave test case at time
t = 0.09 on a 400 × 400 mesh. Left Column: First-order Right Column: Second-order WENO.
Top row: Zero Godunov-Powell source. Middle row: Central Godunov-Powell source. Bottom
row: Upwind Godunov-Powell source.
4. Conclusion
We describe highly robust finite volume schemes for the MHD equations based on the semi-conservative
Godunov-Powell form (1.2). This form is symmetrizable as well as Galilean-invariant, enabling stability estimates
([6]). Furthermore, the form of the semi-conservative version of the equations in several dimensions is similar to
that in one-space dimension, enabling us to design stable one-dimensional schemes that are easily generalized
to multidimensional data. The following new numerical elements are found necessary to obtain the desired
stability:
• We design HLL-type three-wave and five-wave approximate Riemann solvers that can handle non-
constant normal magnetic fields. They can be thought of as extensions of the highly popular three-wave
and five-wave solvers of [20] and [29] respectively.
• These HLL-solvers naturally impose an upwind discretisation of the Godunov-Powell source term.
• The second-order ENO and WENO-type reconstructions we use are modified to preserve positive pres-
sures and densities.
• We provide a second order extension of the upwind discrete source term.
We test the schemes on a variety of numerical experiments in both one- and two-space dimensions. The
results obtained by the schemes (particularly the second-order versions) were impressive, both with respect to
stability and accuracy. In particular, the schemes were able to compute the solutions for the advection of varying
normal magnetic field in one-dimension and the expected orders of convergence were obtained. A super-fast
expansion test illustrated the positivity preservation abilities of the schemes.
The benchmark two-dimensional numerical experiments showed that the schemes were stable, even on very
fine meshes. It is well known that computing two-dimensional MHD on very fine meshes may lead to stability
problems and the schemes of this paper were able to handle these fine mesh resolutions. The accuracy of the
schemes was very satisfactory on all tests. Among the schemes, the HLL five-wave solver was better in terms
of resolution (even at second-order) than the HLL three-wave solver. Similarly the WENO-based schemes were
more accurate than the ENO-based schemes. The highly resolved solutions (obtained on very fine meshes) can
serve as benchmark reference solutions for future computations. Discrete values of divB were low in general and
reduced as the mesh was refined. Non-zero divergence values did not affect either the stability or the accuracy
of the resulting solutions.
Based on the numerical evidence, we conclude that using the semi-conservative form of the MHD equations,
with very careful discretizations of the fluxes and the Godunov-Powell source term together with proper high-
order reconstructions is a very appealing strategy for designing robust schemes for MHD equations. The
numerical tests also indicated that all these ingredients- the introduction of the Godunov-Powell source term,
its careful upwind discretization and suitable positivity preserving reconstructions are all necessary for the design
of stable schemes. These schemes are easier to comprehend and cheaper computationally than the projection
method. Similarly, they are easier to parallelize and to use in conjunction with adaptive mesh refinement than
staggering based methods.
These robust schemes are trivial to extend to three dimensions and we aim to employ them for more realistic
astrophysical simulations in a forthcoming paper. Also, designing higher than second-order schemes is a work
in progress.
References
[1] E. Audusse, F. Bouchut, M. O. Bristeau, R. Klien and B. Perthame. A fast and stable well-balanced scheme with hydrostatic
reconstruction for shallow water flows. SIAM. Jl. Sci. Comp, 25 (6), 2004, 2050 - 2065.
[2] D. S. Balsara and D. Spicer. A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields
in magnetohydrodynamic simulations. J. Comp. Phys., 149(2):270-292, 1999.
[3] D. S. Balsara. Total variation diminishing algorithm for adiabatic and isothermal magnetohydrodynamics. Ap. J. Suppl., 116,
119-, 1998.
[4] D. S. Balsara. Divergence-free reconstructions of magnetic fields and WENO schemes for magnetohydrodynamics. J. Comput.
Phys., 228 (14), 5040-5056, 2009.
[5] D. S. Balsara, T. Rumpf, M. Dumbser and C. D. Munz. Efficient high accuracy ADER-WENO schemes for hydrodynamics
and divergence-free MHD. J. Comput. Phys., 228 (7), 2480-2516, 2009.
[6] T. J. Barth. Numerical methods for gas dynamics systems. In An introduction to recent developments in theory and numerics
for conservation laws, D. Kröner, M. Ohlberger, C. Rohde (eds), Springer, 1999.
[7] C. Berthon Why the MUSCL-Hancock scheme is L1 -stable Numer. Math., 104, 27-46, 2006.
[8] F. Bouchut, C. Klingenberg and K. Waagan. A multi-wave HLL approximate Riemann solver for ideal MHD based on relaxation
I- theoretical framework. Numer. Math., 108 (1), 7-42, 2007.
[9] F. Bouchut, C. Klingenberg and K. Waagan. A multi-wave HLL approximate Riemann solver for ideal MHD based on relaxation
II- Numerical experiments. Preprint, 2008.
[10] J. U. Brackbill and D. C. Barnes. The effect of nonzero divB on the numerical solution of the magnetohydrodynamic equations.
J. Comp. Phys., 35:426-430, 1980.
[11] M. Brio and C. C. Wu. An upwind differencing scheme for the equations of ideal MHD. J. Comp. Phys, 75 (2), 1988, 400 -
422.
[12] P. Cargo and G. Gallice. Roe matrices for ideal MHD and systematic construction of Roe matrices for systems of conservation
laws. J. Comp. Phys, 136 (2), 1997, 446 - 466.
[13] W. Dai and P. R. Woodward. A simple finite difference scheme for multi-dimensional magnetohydrodynamic equations. J.
Comp. Phys., 142(2):331-369, 1998.
[14] A. Dedner, F. Kemm, D. Kröner, C. D. Munz, T. Schnitzer and M. Wesenberg. Hyperbolic divergence cleaning for the MHD
equations. J. Comp. Phys., 175, 645-673, 2002.
[15] B. Einfeldt. On the Godunov type methods for gas dynamics. SIAM. J. Num. Anal., 25 (2), 1988, 294 - 318.
[16] C. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flow: a constrained transport method. Astrophys. J., 332:659,
1998.
[17] F. Fuchs, K. H. Karlsen, S. Mishra and N.H. Risebro. Stable upwind schemes for the Magnetic Induction equation. Preprint,
Math. Model. Num. Anal, to appear.
[18] F. Fuchs, S. Mishra and N. H. Risebro. Splitting based finite volume schemes for the ideal MHD equations. J. Comp. Phys.,
228 (3), 2009, 641-660.
[19] S. K. Godunov. The symmetric form of magnetohydrodynamics equation. Num. Meth. Mech. Cont. Media, 1:26-34, 1972.
[20] K. F. Gurski. An HLLC-type approximate Riemann solver for ideal Magneto-hydro dynamics. SIAM. J. Sci. Comp., 25(6),
2004, 2165-2187.
[21] S. Gottlieb, C. W. Shu and E. Tadmor. High order time discretizations with strong stability property. SIAM. Review, 43, 2001,
89 - 112.
[22] A. Harten, B. Engquist, S. Osher and S. R. Chakravarty. Uniformly high order accurate essentially non-oscillatory schemes.
J. Comput. Phys., 1987, 231-303.
[23] A. Harten, P. D. Lax and B. Van Leer. On upstream differencing and Godunov type schemes for hyperbolic conservation laws.
SIAM Rev., 25 (1), 1983, 35-61.
[24] R. J. LeVeque. Finite volume methods for hyperbolic problems. Cambridge university press, Cambridge, 2002.
[25] R. J. LeVeque. Wave propagation algorithms for multi-dimensional hyperbolic systems, J. Comp. Phys., 131, 327-353, 1997.
[26] T. J. Linde. A three adaptive multi fluid MHD model for the heliosphere. Ph.D thesis, University of Michigan, Ann-Arbor,
1998.
[27] P. Londrillo and L. del Zanna On the divergence-free condition in Godunov-type schemes for ideal magnetohydrodynamics:
the upwind constrained transport method. J. Comp. Phys. 195 (1), 2004, 17 - 48
[28] S. Mishra and E. Tadmor, Constraint preserving schemes using potential-based fluxes. III. Genuinely multi-dimensional central
schemes for MHD equations. Preprint, 2009.
[29] T. Miyoshi and K. Kusano. A multi-state HLL approximate Riemann solver for ideal magneto hydro dynamics. J. Comp.
Phys. 208 (1), 2005, 315 - 344.
[30] B. Perthame and C. W. Shu. On positivity preserving finite volume schemes for Euler equations. Numer. Math., 73 (1), 119 -
130, 1996.
[31] K. G. Powell. An approximate Riemann solver for magneto-hydro dynamics (that works in more than one space dimension).
Technical report, 94 -24, ICASE, Langley, VA, 1994.
[32] K. G. Powell, P. L. Roe. T. J. Linde, T . I. Gombosi and D . L. De zeeuw, A solution adaptive upwind scheme for ideal MHD.
J. Comp. Phys, 154(2), 284 - 309, 1999
[33] P. L. Roe and D. S. Balsara. Notes on the eigensystem of magnetohydrodynamics. SIAM. J. Appl. Math., 56 (1), 1996, 57 -
67.
[34] J. Rossmanith. A wave propagation method with constrained transport for shallow water and ideal magnetohydrodynamics.
Ph.D thesis, University of Washington, Seattle, 2002.
[35] D. S. Ryu, F. Miniati, T. W. Jones and A. Frank. A divergence free upwind code for multidimensional magnetohydrodynamic
flows. Astrophys. J., 509(1):244-255, 1998.
[36] C. W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory schemes - II, J. Comput. Phys., 83, 1989, 32 -
78.
[37] J. M. Stone, T.A. Gardiner, A. Thomas, P. Teuben, J. F. Hawley and J. B. Simon. Athena: A New Code for Astrophysical
MHD. The Astrophysical Journal Supplement Series 178, (1), 2008, 137-177
[38] M. Torrilhon and M. Fey. Constraint-preserving upwind methods for multidimensional advection equations. SIAM. J. Num.
Anal., 42(4):1694-1728, 2004.
[39] M. Torrilhon. Locally divergence preserving upwind finite volume schemes for magnetohyrodynamic equations. SIAM. J. Sci.
Comp., 26(4):1166-1191, 2005.
[40] M. Torrilhon. Uniqueness conditions for Riemann problems of ideal Magneto-hydrodynamics. J. Plasma. Phys., 69 (3), 253-276,
2003.
[41] M. Torrilhon and D.S.Balsara. High-order WENO schemes: investigations on non-uniform convergence for MHD Riemann
problems. J. Comput. Phys., 201 (2), 586-600, 2004.
[42] G. Toth. The divB = 0 constraint in shock capturing magnetohydrodynamics codes. J. Comp. Phys.,161:605-652, 2000.
[43] B. Van Leer. Towards the ultimate conservative difference scheme. V. J. Comp. Phys., 32, 101-136, 1979.
[44] K. Waagan. A positive MUSCL-Hancock scheme for ideal MHD. J. Comp. Phys,. To appear.
(Franz Georg Fuchs) Centre of Mathematics for Applications (CMA)

University of Oslo
P.O. Box 1053, Blindern
N–0316 Oslo, Norway
E-mail address: franzf@math.uio.no
URL: http://folk.uio.no/franzf/
(Andrew D. McMurry) Centre of Mathematics for Applications (CMA)

University of Oslo
E-mail address: a.d.mcmurry@ifi.uio.no
(Siddhartha Mishra) Seminar for Applied Mathematics (SAM),

D-Math, ETH Zürich,
HG G. 57.2, Rämistrasse, 101,
Zürich-8092, Switzerland.
E-mail address: smishra@sam.math.ethz.ch
(Nils Henrik Risebro) Centre of Mathematics for Applications (CMA)

University of Oslo
E-mail address: nilshr@math.uio.no
(Knut Waagan) Center for Scientific Computation and Mathematical Modeling (CSCAMM)
The University of Maryland
CSCAMM
4146 CSIC Building #406
Paint Branch Drive
College Park MD 20742-3289
USA
E-mail address: knutwa@ucar.edu
Research Reports
No. Authors/Title
09-37 F.G. Fuchs, A.D. McMurray, S. Mishra, N.H. Risebrom, K. Waagan

Approximate Riemann solvers and stable high-order finite volume
schemes for multi-dimensional ideal MHD
09-36 Ph. LeFloch, S. Mishra
Kinetic functions in magnetohydrodynamics with resistivity and hall
effects
09-35 U.S. Fjordholm, S. Mishra
Vorticity preserving finite volume schemes for the shallow water equations
09-34 S. Mishra, E. Tadmor
Potential based constraint preserving genuinely multi-dimensional
schemes for systems of conservation laws
Constraint preserving schemes using potential-based fluxes.
III. Genuinely multi-dimensional central schemes for for MHD equations
II. Genuinely multi-dimensional central schemes for systems of conserva-
tion laws
I. Multidimensional transport equations
09-30 D. Braess, S. Sauter, C. Schwab
On the justification of plate models
09-29 D. Schötzau, C. Schwab, T. Wihler
hp-dGFEM for second-order elliptic problems in polyhedra.
II: Exponential convergence
09-28 D. Schötzau, C. Schwab, T. Wihler
hp-dGFEM for second-order elliptic problems in polyhedra.
I: Stability and quasioptimality on geometric meshes
09-27 A. Moiola, R. Hiptmair, I. Perugia
Approximation by plane waves
09-26 M. Karow, E. Kokiopoulou, D. Kressner
On the computation of structured singular values and pseudospectra
09-25 M. Durán, M. Guarini, C.F. Jerez-Hanckes
Hybrid FEM/BEM modeling of finite-sized photonic crystals for semi-
conductor laser beams

A000 FVM

Uploaded by

Copyright:

Available Formats

A000 FVM

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

A000 FVM

Uploaded by

Copyright:

Available Formats

!!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!! !!!!!!!!!!!

Eidgenössische Ecole polytechnique fédérale de Zurich

Approximate Riemann solvers and stable

F.G. Fuchs∗ , A.D. McMurry∗ , S. Mishra, N.H. Risebro∗ and K. Waagan†

Research Report No. 2009-37

Seminar für Angewandte Mathematik

F. G. FUCHS, A. D. MCMURRY, S. MISHRA, N. H. RISEBRO, AND K. WAAGAN

Date: November 16, 2009.

Using standard vector identities results in the following semi-conservative form,

The flux Fni+1/2,j = F(Wi,j

for θ ∈ {L, R}. Note that π1L

S1,n 1,∗∗ 1,∗∗

M H3 rate H3 E rate H3 W rate H5 rate H5 E rate H5 W rate

M H3 rate H3 E rate H3 W rate H5 rate H5 E rate H5 W rate

M H3 rate H3 E rate H3 W rate H5 rate H5 E rate H5 W rate

M H3 rate H3 E rate H3 W E rate H5 rate H5 E rate H5 W rate

(Franz Georg Fuchs) Centre of Mathematics for Applications (CMA)

(Andrew D. McMurry) Centre of Mathematics for Applications (CMA)

(Siddhartha Mishra) Seminar for Applied Mathematics (SAM),

(Nils Henrik Risebro) Centre of Mathematics for Applications (CMA)

09-37 F.G. Fuchs, A.D. McMurray, S. Mishra, N.H. Risebrom, K. Waagan

You might also like