Application of ADER Scheme in MHD Simulation ∗

The Arbitrary accuracy Derivatives Riemann problem method (ADER) scheme is a new high order numerical scheme based on the concept of finite volume integration, and it is very easy to be extended up to any order of space and time accuracy by using a Taylor time expansion at the cell interface position. So far the approach has been applied successfully to flow mechanics problems. Our objective here is to carry out the extension of multidimensional ADER schemes to multidimensional MHD systems of conservation laws by calculating several MHD problems in one and two dimensions: (i) Brio-Wu shock tube problem, (ii) Dai-Woodward shock tube problem, (iii) Orszag-Tang MHD vortex problem. The numerical results prove that the ADER scheme possesses the ability to solve MHD problem, remains high order accuracy both in space and time, keeps precise in capturing the shock. Meanwhile, the compared tests show that the ADER scheme can restrain the oscillation and obtain the high order non-oscillatory result.


Introduction
The classical Riemann problem is the Cauchy problem for a system of hyperbolic conservation laws, with initial condition consisting of two constant states separated by a discontinuity.The self-similar solution of this Riemann problem was first used by Godunov [1−2] to construct his first-order upwind numerical flux.
Since then, methods to solve the classical Riemann problem have been studied, for example, in Ref. [3].
One of the developed methods is to extend the local Riemann problem to the Generalized Riemann problem at cell interface position.The so-call Generalized Riemann problem, in which the initial condition consists of two polynomials of the first degree (vectors) separated by a discontinuity at the interface, was first constructed by Ben-Artzi et al. [4] in one dimension.
According to Toro's concept of the Generalized Riemann problem [5] , the initial condition consists of two kth order polynomial functions, which are denoted as the corresponding Generalized Riemann problem by GRP k .For example, GRP 0 means that all first and higher-order spatial derivatives of the initial condition for the GRP away from the origin vanish identically, which corresponds to the classical piece-wise constant data Riemann problem.Similarly, GRP 1 means that all second and higher-order spatial derivatives of the initial condition for the GRP away from the origin vanish identically, which corresponds to the piecelinear data Riemann problem.The most general case is that the initial conditions are two arbitrary but infinitely differentiable functions, which are denoted by GRP ∞ .
A major simplification to the GRP methodology comes with the Modified GRP (MGRP) scheme, proposed by Toro et al. [6] In this scheme the GRP is solved by two conventional Riemann problems, namely one non-linear problem for the leading term for state variables and one linear problem for gradients of state variables.The Arbitrary accuracy Derivatives Riemann problem method (ADER) approach can be regarded as a further development of the MGRP scheme in that it breaks the barrier of second-order accuracy and allows the construction of arbitrarily high-order accurate schemes, both in time and space.This approach was put forward firstly by Toro and his collaborators [7] , where they achieved 10th order of accuracy in both space and time.The ADER scheme has been applied successfully to the hydrodynamic systems.Titarev et al. [3,5,8,15] applied this scheme to the scalar advection-reaction-diffusion equations.Toro et al. [9] used TVD fluxes for the highorder ADER schemes.Toro et al. [10−14] used this scheme for linear and non-linear hyperbolic conservation laws from one-to three-dimensional hydrodynamic problem.Motivated by the ADER's successful application to hydrodynamic problem, we try to implement the scheme for MHD simulation.We apply the scheme to several MHD problems which are calculated in one and two dimensions: (i) Brio-Wu shock tube problem, (ii) Dai-Woodward shock tube problem, (iii) Orszag-Tang MHD vortex problem.The numerical results show that the ADER scheme possesses the ability to solve MHD supported by prob-lem, and remains high order accuracy both in space and time.
1 Solution Method

Numerical Scheme in One Dimension
Consider the one dimensional ideal MHD equations of conservation form: where is the vector of unknown conservative variables and F (U ) is the flux vector, where ρ is the density, is the total energy, and p is the pressure.Integrating Eq. ( 1) over a space-time control volume in x-t space we obtain the following one-step finite-volume scheme: Here U n i is the cell average of the solution at time level t n , and can be given by the initial value.F i+1/2 is the time average of the physical flux at cell interface x i+1/2 .They can be given respectively: For convenience, we replace t by τ = t−t n .Then we can change the numerical flux as: Then we can obtain the value of U (x i+1/2 , τ) as follows: To evaluate the leading term, we write a Taylor expansion of the interface state in time where U (x i+1/2 , 0 + ) = lim t→0 + U (x i+1/2 , t) and K is the order of the approximation.
Eq. ( 5) consists of two parts: the leading term U (x i+1/2 , 0 + ) and the higher order terms ∂ k ∂t k U (x i+1/2 , 0 + ), which need to be solved respectively.The leading term accounts for the firstinstant interaction of the boundary extrapolated values U L (x i+1/2 ) and U R (x i+1/2 ), and can be obtained by the conventional Riemann problem: A key ingredient here is the availability of approximate Riemann solver to provide this first term in the expansion.In this paper, we use the HLL's Riemann solver, which was put forward by Harten, Lax and Van Leer [16] .
The higher order terms are evaluated in two steps.First we express all time derivatives via spatial derivatives by means of the Cauchy-Kowalewski procedure [3] .For system Eq.( 1) the procedure yields the following expressions: and so on.From Eq. ( 7), we can obtain all the high order time derivatives required by Eq. ( 5) if the spatial derivatives are known.Thus, our next step is to calculate the spatial derivatives.The evolution equations for the spatial derivatives can be derived by differentiating the governing Eq. ( 1).In general, the obtained evolution equation for each spatial derivative is in non-conservative form and contains a nonlinear source term H depending on spatial derivatives: where the coefficient matrix A(U ) is precisely the Jacobian matrix ∂F ∂U .For the Taylor expansion Eq. ( 5) we need the solution of Eq. ( 8) Note that the coefficient matrix A i+1/2 is the same for all derivatives and has to be evaluated only once.
Finally, having found all spatial derivatives we form the Taylor expansion Eq. ( 5).Then there are two options now existing to evaluate the numerical flux.The first option is the state-expansion ADER [15] , in which the approximate state Eq. ( 5) is inserted into the definition of the numerical flux Eq. ( 3) and then an appropriate Kth-order accurate quadrature is used for time integration: Here γ l and ω l are properly scaled nodes and weights of the rule and K α is the number of nodes.
The second option to evaluate the numerical flux is the flux-expansion ADER [5,12] , in which we seek Taylor expansion of the physical flux at x i+1/2 : From Eq. (3) and Eq. ( 11) the numerical flux is now given by The leading term F (x i+1/2 , 0 + ) accounts for the first interaction of left and right boundary extrapolated values and is computed as a certain monotone flux of the conventional Riemann problem Eq. ( 6) for the leading term of the state expansion Eq. ( 5).The remaining higher order time derivatives of the flux in Eq. ( 12) are expressed via time derivatives of the intercell state ∂ k ∂t k U (x i+1/2 , 0 + ), which are known from Eq. ( 5).
In our calculations, we choose the second option Eq. (12).And then, the solution is advanced in time by updating the cell averages according to the onestep formula Eq. (2).

Numerical Scheme in Two Dimensions
The process is the same as one-dimensional situation.Consider the following two dimensional nonlinear system of conservation laws: Integration of Eq. ( 13) over a space-time control vol- produces the following one-step finite-volume scheme: Here U n i,j is the cell average of the solution at time level: And F i+1/2,j , G i,j+1/2 are the space-time averages of the physical fluxes at the cell interfaces: While describing the procedure to evaluate the numerical flux in two dimensions we concentrate on F i+1/2,j ; the expression for G i,j+1/2 is obtained in an entirely analogous manner.
The evaluation of the ADER numerical flux consists of the following steps.First we discretize the spatial integrals over the cell faces in Eq. ( 16) using a tensor product of a suitable Gaussian numerical quadrature.The expression for the numerical flux in the x coordinate direction then reads Where y α is the integration points over the cell face [y j−1/2 , y j+1/2 ] and K α is the weights computed by the Gauss numerical quadrature.Normally, we use the two-point Gaussian quadrature for third and fourth order schemes and a higher-order Gaussian quadrature for fifth and higher order schemes.
The following steps are provided only for xdirection computation and the same procedure follows for y-direction computation.
We write Taylor series expansion in time where τ = t − t n .The leading term U (x i+1/2 , y α , 0 + ) is the Godunov state of the conventional Riemann problem To evaluate higher-order terms we first express all time derivatives by spatial derivatives by means of the Cauchy-Kowalewski procedure.We can derive homogeneous evolution equations and the initial conditions for each spatial derivative The spatial derivatives at (x − x i+1/2 )/τ = 0 are then the Godunov states of the following linearized Riemann problem with piece-wise constant initial data: After solving Eq. ( 20) for 1 k + l K, we form the Taylor expansion Eq. ( 18) for the interface state at the Gaussian integration point (x i+1/2 , y α ).The flux of the state-expansion ADER scheme is obtained by inserting the approximate state Eq. ( 18) into formula Eq. ( 17) and using an appropriate Kth-order accurate quadrature for time integration: For the flux expansion ADER schemes we write Taylor time expansion of the physical flux at each point From Eq. ( 16) and Eq. ( 22) the numerical flux is given by Entirely analogous to the one-dimensional case, the leading term F (x i+1/2 , y α , 0 + ) is computed from Eq. ( 19) using a monotone upwind flux.The remaining higher order time derivatives of the flux in Eq. (22)   are expressed via time derivatives of the intercell state U (x i+1/2 , y α , 0 + ) which are given by the Taylor expansion Eq. ( 18).The solution is advanced in time by updating the cell averages according to the one-step formula Eq. ( 14).
In the ADER scheme solving method, the key process is to solve the Riemann problem.There are many methods to solve the Riemann problem, such as the MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) scheme [17] , PPM (piecewise parabolic method) scheme [18] , and HLL scheme [13] .
The HLL scheme used in the present paper to solve the Riemann problem can be described as follows: The interface value and the numeral flux have the form as: where where vxi and cfxi are the Roe average of velocity and fast magnetic sound velocity.ρi = Here H is the enthalpy, has the form as: where z , then the Roe average are: The fast magnetic sound velocity where and then cfxi = c fxi (ρ i , vαi , Bβi , Hi ).

Brio-Wu Shock Tube Problem
Proposed by Brio and Wu [19−20] , this MHD problem has been a benchmark test for solving the one dimensional ideal MHD equations.Similar to Sod's shocktube problem [19] , the initial condition is composed of two distinct constant states: In addition, B x = 0.75, γ = 2.0 is a constant value in the whole domain [0, 1].According to Brio and Wu et al. [19] , numerical solutions with a uniform mesh of 800 grid points and the CFL=0.8 are obtained.Figure 1  And the CFL number used for these three calculations is the same.It is clear to see that, mesh refinement can improve numerical accuracy.

Dai-Woodward Shock Tube Problem
Dai and Woodward's shock tube problem has been studied by Dai and Woodward [21] and Bouchut et al. [20] The initial data [20] are given as:  By comparison with that of Bouchut et al.'s results [20] , again the MHD solution with ADER scheme is quite good.Under the same CFL condition, different resolutions have been examined with 1000, 2000, 3000 grid points respectively.Figure 2(e) shows comparison of density profiles.The accuracy is nearly the same with different resolution.So in this case, the mesh refinement does not improve numerical accuracy.

Orszag-Tang MHD Vortex Problem
As our preliminary effort to apply the ADER method for multidimensional MHD problem, the Orszag-Tang MHD vortex problem [22] is adopted here.This compressible flow proposed by Orszag and Tang contains significant features of MHD turbulence, and has been taken as a standard numerical test [22−26] .The flow involves complex evolution due to interactions between shock waves and the vortices.
Here we used γ = 5/3.The problem was calculated with a CFL number of 0.2.The computational domain is [0, 2π] × [0, 2π] with periodic boundary condition in both x and y directions.According to Zhang et al. [25] , we use a uniform mesh of 200×200 grid points.
We have successfully run our calculation from t = 0.5 to t = 3, which is the final time in most of previously published results.Figure 3(a), 3(b) show the density and pressure contours at t = 0.5 using the same contour level as those in Jiang and Wu [26] .Clearly, our results are similar to Jiang and Wu's results.
For the sake of comparison, we show the density contours at t = 0.5, t = 1, and t = 3 using the same contour level, shown in Figure 3 To compare the mesh refinement, we calculated the density of 200 × 200, 300 × 300 and 400 × 400 grid points, see Figure 3(c) (e) (f).Evidently, these results are similar to those obtained formerly [22−26] .
It is clear that the effect of mesh refinement is not very palpable.In MHD simulation, how to control the numerical error of magnetic field divergence is a challenge [27] .In order to see how this numerical error evolves, Figure 3(g) expresses the evolution of the is the number of mesh nodes in the computational domain) [27] .This figure tells that the numerical error of magnetic field divergence is tolerable.

Comparison Between HLL and ADER
In order to see the difference of ADER scheme and HLL scheme, we calculate the following cases: using second-order HLL method only and ADER method introduced in Section 2, respectively.We test Brio-Wu shock tube problem through these two methods.Figure 4 number is small, the two methods show no differences.
But when the CFL number is large, the HLL scheme turns up some oscillation at the slow shock, while the ADER scheme does a good job.In conclusion, the HLL scheme has higher oscillatory result, when CFL number is large, and the ADER scheme can avoid this oscillatory and capture the shock for large CFL number.

Conclusion
In this paper, we have presented the ADER scheme for ideal MHD equations in one and two spatial dimension.Contrast to modern Godunov-type, the ADER method is one-step, fully discrete scheme for which the reconstruction procedure is carried out only once per time step, and can be implemented in the In conclusion, the ADER scheme can be applied in the MHD numerical simulation, and keeps very high order accuracy both in space and time.
According to these conclusions, future applications of the ADER scheme to the MHD numerical simulations will be a promising work, such as its applications to magnetic reconnection, solar wind simulation and solar disturbances in solar-terrestrial space [23,27−28] .At the same time, improving the ADER scheme's order by using different reconstruction methods, such as Essentially Non-Oscillatory (ENO), WENO or CWENO scheme [23,29−34] , is also an emphasis in future work.
/2 at the time τ = 0 + .Therefore, the source term comes into effect for τ > 0 only, and can be neglected.Additionally, we linearize the equation around the leading term U (x i+1/2 , 0 + ) of the time expansion Eq. (5) and replace the piece-wise polynomial initial data by left and right boundary extrapolated values of spatial derivatives at x i+1/2 .The described simplifications result in the following linear conventional Riemann problem for the spatial derivatives U (k) : (a)∼(d) show numerical solution of density, velocity v x , magnetic field B y and pressure at 400 time steps.From this figure, we can see that the wave family consists of right-moving and left-moving waves.The left-moving waves are a fast rarefaction wave, denoted by FR in the figure, and a slow compound wave, denoted by SM.The right-moving waves included a contact discontinuity, denoted by C, a Slow Shock (SS), and a Fast Rarefaction wave, FR.Our result shows sharp resolution of all shock waves.

Figure 1 (
Figure 1(e) shows the result of grid refinement

Fig. 1
Fig. 1 MHD solution by ADER scheme for Brio-Wu shock tube problem.(a) density profile, (b) velocity vx profile, (c) magnetic field By profile, (d) pressure profile, (e) comparison between different resolutions

Fig. 2
Fig. 2 MHD solution by ADER scheme for Dai-Woodward shock tube problem.(a) density profile, (b) velocity vy profile, (c) magnetic field By profile, (d) pressure profile, (e) comparison between different resolution (a) (c) (d) respectively.Figure 3(a) (c) (d) show very complex shock structure.At t = 1, an intermediate shock is formed at the shock front in the region of π < x < 1.5π and 0 < y < 0.75π.
(a) (b) are the density solved by HLL scheme and ADER scheme with CFL= 0.3 and CFL= 0.6 at t = 0.1.Through the two figures we can see that, when the Courant-Friedrichs-Lewy (CFL)

Fig. 3 Fig. 4
Fig. 3 MHD solution by the ADER method for Orszag-Tang MHD vortex problem.(a) Density contours of 200 ×200 grid points at t = 0.5.(b) Pressure contours of 200×200 grid points at t = 0.5.(c) Density contours of 200×200 grid points at t = 1.0.(d) Density contours of 200×200 grid points at t = 3.0.(e) Density contours of 300×300 grid points at t = 1.0.(f) Density contours of 400×400 grid points at t = 1.0.(g) Evolution of the error for ∇ • B in the three cases of the mesh refinement