Splitting Based Scheme for Three-dimensional MHD with Dual Time Stepping ∗

A new hybrid numerical scheme of combining an E-CUSP (Energy-Convective Upwind and Split Pressure) method for the fluid part and the Constrained Transport (CT) for the magnetic induction part is proposed. In order to avoid the occurrence of negative pressure in the reconstructed profiles and its updated value, a positivity preserving method is provided. Furthermore, the MHD equations are solved at each physical time step by advancing in pseudo time. The use of dual time stepping is beneficial in the computation since the use of dual time stepping allows the physical time step not to be limited by the corresponding values in the smallest cell and to be selected based on the numerical accuracy criterion. This newly established hybrid scheme combined with positivity preserving method and dual time technique has demonstrated the accurateness and robustness through numerical experiments of benchmark problems such as the 2D Orszag-Tang vortex problem and the 3D shock-cloud interaction problem.


Introduction
In recent years, the Convective Upwind and Split Pressure (CUSP) family schemes have achieved great success in CFD simulations.The characteristic of CUSP schemes is that it simultaneously consider the convective upwind characteristics and avoid the complex matrix dissipation such as that of the Roe's flux difference splitting scheme.The CUSP family can be basically categorized into two types: the H-CUSP and E-CUSP [1−2] .The H-CUSP scheme has total enthalpy in the energy equation in the convective vector, so the scheme is not fully consistent with the disturbance propagation that may affect the stability and robust of the H-CUSP scheme.While, the E-CUSP schemes use total energy in the convective terms, which result in splitting the eigenvalues of the Jacobian to convection (velocity) and waves (speed of sound).
In reality, the fluid part of the MHD leads to visiting the neighboring zones that have connectivity with a zone of interest, we are able to identify minimum and maximum values for density and pressure that should bound the reconstructed profile within the zone of interest [5] .Then, weighted mean of the reconstructed conserved variable and the conserved variable is used to correct the reconstructed variable.The self-adjusting positivity preserving method is easy and inexpensive to implement and the result show that this method can efficiently enhance the robustness.
In the three-dimensional simulation, time steps should be dictated by numerical stability, so they are much smaller than required by accuracy considerations, which will increase computer time for conditionally stable time-marching schemes.To improve computational efficiency, schemes which can use larger possible time-step sizes permitted by accuracy considerations should be taken into account.
Here, we use the dual time step scheme to update the flow and magnetic conserved variables.The use of dual time stepping allows the physical time not to be limited by the corresponding values in the smallest cell and to be selected based on the numerical accuracy criterion [8] .It is found to be efficient in terms of memory required and computing effort per time step.This is due to the fact that no matrix manipulation is required by the scheme [9] .Dual time step, which do not modify the original transient evolution of the governing equation, adds to the governing equation a pseudo-time derivative.It uses the pseudo-time steady-state solution to approach the physical-time solution.The advantage of dual time step is that it can speed up the computational efficiency and extend the range of plasma value.

Numerical Method for MHD Equations
The ideal MHD equation for inviscid flow can be expressed as: where ρ is the flow density, V = (u, v, w) is the flow velocity, B is the magnetic field, p is the pressure, ρe is the total energy, and p t = 1 2 B 2 + p.In this paper, following the idea of Fuchs [3] we split the MHD equations into a fluid part Eq. ( 1)∼( 3) and a magnetic induction part Eq. ( 4), and employ the E-CUSP for the fluid part and CT scheme for the magnetic induction part in order to get a scheme for full MHD equations [3] .The fluid part leads to an extended Euler system with magnetic forces as source terms.This set of equations is approximated by E-CUSP [4] .The magnetic part is modeled by the magnetic induction equations which are solved using the constrained transport method [6−7] , while preserving the divergence-free constraint.These two sets of schemes can be combined either component by component, or by using an operator splitting procedure to produce a full scheme for the MHD equations.This is, piecing them together, the resulting update full scheme for MHD equations reads . In the present paper, we implement each item of our hybrid scheme by using E-CUSP for the fluid part, CT for the magnetic induction part and dual time integration strategy.
We can rewrite the governing Eq. ( 1)∼( 3) in the Cartesian coordinates with compact form as with each term having its usual meaning.The semi-discretized conservative three-dimensional MHD equation reads In an E-CUSP scheme for hydrodynamics, the flux are divided into two parts: the convective terms and the acoustic waves propagating in each direction at subsonic regime.The basic idea has been extended by Shen et al. [4] to the MHD because the MHD wavelike structure is analogous to that of hydrodynamics.
Following the E-CUSP scheme of Zha et al. [10] , we decompose the flux F to convective and generalized wave fluxes: where The numerical flux at the interface is evaluated according to that given by Shen et al [4] .
For numerical simulation of the MHD equation, how to preserve the divergence-free condition of the magnetic filed ∇ • B = 0 is a curial issue.There are several approaches to deal with this problem [11−12,14−15] , the constraint transport is one of the approaches.In this paper, the constraint transport scheme suggested by Ziegler [6−7] is used.
As usual, in order to achieve second order spatial accuracy, flux reconstruction technique should be employed.Here, the Monotone Upstream-Centered (MUSCL) reconstruction scheme is used to reconstruct the flow variable.Following Ziegler [6−7] , the reconstruction of magnetic field can be achieved by using values at the centers of the cell surfaces.
Positivity of density and pressure may be lost in the unknowns' updating step at next time step somewhere within a zone or during the reconstruction procedure in obtaining the reconstructed values at the cell interface.In order to avoid the negativeness of density and pressure in the reconstruction and updating procedure, we take the following methods.First, to preserve positivity for constructed variables, we 2015, 35 (1)   follow the self-adjusting positivity-preserving method proposed by Balsara [5] .His method examines the local flow to identify regions with strong shocks.The permitted range of densities and pressures is also obtained at each zone by examining neighboring zones.
The range is expanded if the solution is free of strong shocks in order to accommodate higher order nonoscillatory reconstructions.The density and pressure are then brought into the permitted range.For details, refer to Balsara [5] .Second, to insure the positivity of the updated variables for each new time level, if the pressure from the energy is negative, we turn to solve the entropy equation Furthermore, If the pressure is still negative, we use the method of Keppens et al. [16] by suggesting the positivity fixing strategy as follows: (1) identify all cells that represent physical states surrounding a faulty cell; (2) convert those cells to primitive variables; (3) For all but the magnetic field components, replace the faulty cell values by the average of surrounding physical state cells.Finally, after all the above methods, if there are still some points which have negative pressure, we use the values at the previous time step to replace those for this time step in these points.
To improve the computational efficiency, we use a special implicit scheme which is matric-free implicit dual time-stepping scheme to advance time marching.The scheme added to the governing equations a pseudo-time derivative that emulates the original physical-time derivative.In reality, the implicit dual time step scheme [18] can use a lager time step than the explicit scheme, even dozens of time step that of explicit schemes.Thus the scheme can speed up the simulation.
In the implicit physical time steeping, the second order backward difference in time formula is used to discretize ∂/∂t in Eq. ( 6): where the superscript n + 1 is the time level at (n + 1)δt, and all the variables with superscript n + 1 are evaluated at this time level.Then, we add the derivative of a fictitious pseudo-time τ to the above Eq.( 8): where, n is the physical time level, m is the pseudotime level (the number of subiterations), δτ is the pseudotime step size, and δt is the physical time step.
δt will be set to be some integer times the explicit time step Δt constrained by the usual Courant-Friedrichs- Here, c fx , c fy , c fz are defined as usual to be the fast magnetosonic speeds in x, y, z directions, respectively.
After the residual terms in the righthand of the above equation are linearized at the m + 1 pseudo-time level with respect to the previous pseudotime level m, we arrive at an unfactored implicit form as follows. .
The matrix B and C are like A. Similarly, the above dual time procedure applies to the CT scheme for the magnetic induction part.
These equations are marched in pseudo time until the derivation of U i,j,k with respect to τ converges to zero.In simulation, the equations are iterated in pseudo time so that U n+1,m+1 approaches the physi- This method described above is named the dual time hybrid scheme.
In the dual time stepping procedure, a physicaltime accurate solution is generated upon convergence towards pseudo-time steady-state per physical-time step [17] .In practice, it is not necessary for the the derivation of U i,j,k with respect to τ to approach zero.As pointed out by Zhao and Li [17] , the convergence criterion can be prescribed by requiring that (CFL sub is the CFL number, V i,j,k represent the cell volume, λ A , λ B , λ C represent the maximum eigenvalue in x, y, z direction, respectively) is used such that the matrix involved in the solver is diagonal dominant and the Gauss-Seidel line iteration can converge.
Since backward differentiation in time integration involves three time steps, namely n − 1, n and n + 1, a startup procedure is required.Hence, for the first iteration, the explicit second-order Runge-Kutta time stepping that involves two time steps, n and n+1, is used to integrate Eq. ( 4) and Eq. ( 6) with variable reconstruction given above (hereafter called the explicit hybrid scheme).

Explicit Hybrid Scheme
The interaction of a strong shock wave with a density clump in a magnetic environment is a problem of astrophysical relevance [6] .The computational domain is a Cartesian box given by (x, y, z) ∈ [−0.5, 0.5] 3 with the initial condition from Ziegler [6] .There is a discontinuity at x = 0.1 with left and right states: used for all the variables.
The simulation are carried out on 100×100×100 grid and stopped at t = 0.06 after a violent collision between the shock and clump has taken place.The resulting density distribution and the magnitude of B x in the x-y plane runs in Fig. 1.In our simulation, negative pressure has been seen in the reconstructed profiles, but after the remedy with the selfadjusting positivity-preserving method, the negative pressure in the reconstructed profiles has gone.Our results are in good agreement with the computations in Ziegler [6] .

2D: Orszag-Tang Vortex with Dual Time Hybrid Scheme
Because the Orszag-Tang vortex problem has interactions of multiple shock waves generated as the vortex evolves, it is considered as one of the standard models to validate a MHD numerical method [4,19] .
The Orszag-Tang vortex problem is given by the giving initial conditions on a square domain From Fig. 2 and 3, it is evident that at t = 3, a fast shock front is formed in the region of 1.25π < x < 1.5π and 0 < y < 0.75π, a slow shock front is formed in the region of 0.5π < x < π and 0.5π < y < 0.75π.
We found that our results are almost identical to those of Ref. [4, 20-24], although different grids are used.
Fig. 4 shows the influence of different CFL conditions on the pressure distributions along the line of y = 1.0 with the different mesh grids: 200 × 200 (Fig. 4a) and 400 × 400 (Fig. 4b).From Fig. 2 and 3, we can not discover any difference of the pressure and B x with different CFL.But from Fig. 4, we can find the dual time stepping scheme will increase the maximum and the minimum value of the vortex problem slightly without changing the general structure of the problem.In total, the scheme works well.
From our experiments, it is found that the convergence criterion for the inner iteration loop is significant in simulating unsteady MHD flows.In the inner iteration, using a given convergence criterion is suggested; if using a maximum iteration step, to ensure the time accuracy and at the same time improve the computational efficiency, we suggest that the maximum iterations are limited in different areas.Near large gradient regions, increasing the inner maximum

Conclusions
Here is proposed a second order accuracy hybrid numerical scheme of combing E-CUSP scheme for the fluid part in MHD equations and the constrained transport method for the magnetic induction part.
This newly established scheme can avoid the complex eigenstructure of the Jacobian matrices.Backward  more, in a pseudo-time, the acceleration techniques, such as multigrid [25] , time-derivation preconditioning [26] , can be used without changing the properties of the physical-time.These considerations, combined with the present method's application to the numerical simulation of solar corona [11][12][13]15] , will be left for our future work.

Fig. 1 Fig. 2
Fig. 1 Pressure and Bx contours in the x-y plane of shock-cloud interaction problem with the explicit hybrid scheme with grid mesh 100 × 100 × 100 at t = 0.06 scheme can extend the limitation of the plasma beta in the strong magnetic field region, such as solar coronal magnetic field MHD reconstruction, because of the admittance of enlarged CFL condition.What's

Fig. 3 Fig. 4
Fig. 3 Pressure contour (a) and Bx contours (b) of Orszag-Tang MHD problem with different CFL valueswith grid mesh 400 × 400 at t = 3.0.The CFL number of the first row is 0.3, the second row is 6.0, the third row is 9.0, and the fourth row is 18.0