Application of Central Upwind Scheme for Solving Special Relativistic Hydrodynamic Equations

The accurate modeling of various features in high energy astrophysical scenarios requires the solution of the Einstein equations together with those of special relativistic hydrodynamics (SRHD). Such models are more complicated than the non-relativistic ones due to the nonlinear relations between the conserved and state variables. A high-resolution shock-capturing central upwind scheme is implemented to solve the given set of equations. The proposed technique uses the precise information of local propagation speeds to avoid the excessive numerical diffusion. The second order accuracy of the scheme is obtained with the use of MUSCL-type initial reconstruction and Runge-Kutta time stepping method. After a discussion of the equations solved and of the techniques employed, a series of one and two-dimensional test problems are carried out. To validate the method and assess its accuracy, the staggered central and the kinetic flux-vector splitting schemes are also applied to the same model. The scheme is robust and efficient. Its results are comparable to those obtained from the sophisticated algorithms, even in the case of highly relativistic two-dimensional test problems.


Introduction
The special relativistic hydrodynamical models can be used to simulate many high energy phenomena in astrophysics, including accretion flows, gamma-ray bursts, and jet flows [1,2]. Moreover, free-electron laser technology, high energy particles beams and heavy-ion collisions can be modeled by using special relativistic hydrodynamics (SRHD) [3]. Since last decade, several upwind high-resolution shock-capturing (HRSC) schemes were being applied to solve relativistic hydrodynamical models. The schemes have high order accuracy in smooth regions of the simulated flow and resolve sharp discontinuous profiles in the shock regions [1,[4][5][6][7][8][9]. A part from the kinetic scheme which is based on kinetic equation and equilibrium distribution function [10][11][12][13][14], all these techniques are based on the macroscopic continuum description. Moreover, Qamar and Yousaf have implemented the space-time CESE method [15,16] and discontinuous Galerkin finite element method [17] for solving these RHD equations.
Central schemes could serve as a common numerical technique to solve several scientific and engineering problems due to avoiding the specific eigenstructure of the problem [18]. The schemes have been successfully applied to solve problems in computational fluid dynamics, astrophysics, metrology, semiconductors, shallow flow and multi-component flows [13,19]. The first-order LaxFriedrichs scheme is the the back bone of such schemes. The central Nessyahu-Tadmor (NT) scheme [20] is a Riemann-solver-free high resolution staggered central scheme which do not need any information about the eignstructure of the problem. Thus, the method can be implemented as a black box solver to any system of conservation laws. However, this family of central schemes suffers from excessive numerical viscosity when a sufficiently small time step is imposed, e.g., due to the presence of degenerate diffusion terms. To overcome this deficiency, Kurganov and Tadmor [18,21] improved the NT scheme by using the correct information of local propagation speeds and obtained the semi-discrete central upwind scheme. Similar to the staggered NT scheme, it enjoys the benefits of high resolution, simplicity and robustness. However, the central upwind scheme reduces large amount of numerical dissipation present in the NT central schemes. This scheme has already been applied to different problems namely, two-layer shallow water equations [22] and Hamilton Jacobi equations [23].
In this article, the central upwind scheme is implemented for solving the SRHD equations in one and two space dimensions. For validation and comparison, the numerical results of central(NT) and KFVS schemes are also presented [20,24]. Several numerical cases studies are carried out. It was found that central upwind scheme gives comparable solutions to the KFVS scheme and is superior over the central NT scheme.
The present paper is organized as follows. Section 2 gives a brief description of SRHD equations. The one-dimensional central upwind scheme is briefly presented in Section 3. Afterwards, the scheme is extended to two-space dimensions. In Section 4, one and two-dimensional numerical test problems are given. Conclusions and remarks are offered in Section 5.

Description of Special Relativistic Hydrodynamics
In this section, we give a general technical description of the SRHD equation as they will be used in the analytic description of different problems and development of our numerical code.
The SRHD equation are obtained from the local conservation laws of the stress-energy tensor T μν and the mass flux vector N μ @N m @x m ¼ 0 and @T mn @x m ¼ 0; ð1Þ as @ @x m denotes for the covariant derivative associating with the four space-time coordinates x μ = (t, x i ). Throughout Greek indices (e.g., μ, ν) represent the space-time components while Latin indices run from 1 to 3 and units in which c (the speed of light) = 1 are used. The components of N μ and T μν on a coordinate basis are where, u μ representing the 4-velocity of fluid and ρ the rest-mass density, p is the pressure in a locally inertial reference frame and g μν = diag(1, −1, −1, −1) is the metric tensor in Minkowski Space. The system composed by the continuity equation and the equation of motion must be supplemented by an equation of state relating the pressure, e.g.
where, Γ is adiabatic index with G ¼ 5 3 in the case of mildly relativistic.
In Minkowski Spacetime and cartesian coordinates, the local conservation laws describing the motion of a relativistic fluid can be cast as a first-order flux-conservative system of the form @W @t þ @FðWÞ @x þ @GðWÞ @y ¼ 0 ; ð4Þ by introducing the conservative variables W = (W 1 , W 2 , W 3 , W 4 ) T and the fluxes F and G are given by where, v i is the 3-velocity and γ is the Lorentz factor satisfying g ¼ u 0 ¼ 1 ffiffiffiffiffiffiffiffiffi The 3-velocity components are obtained from the spatial components of the 4-velocity as v i ¼ u i u 0 (normalized as u μ u μ = 1) and the conservative variables are related with physical or primitive variables, w = (ρ, v 1 , v 2 , p) T .

Recovery of physical variables
We need to convert conservative variables to physical variables. Here, W m , m = 1, 2, 3, 4, will denote the conservative variables while are the physical variables. From Eqs (5)-(7), it is concluded that each W m and each flux functions are the functions of the physical variables. To show that each physical variable is an implicit function of the conservative variables W m , let and U ¼ ðe þ pÞg 2 and G 1 ¼ G=ðG À 1Þ : Then with the help of Eqs (3) and (5), we get and By eliminating the parameter p from Eqs (10) and (11), one has Moreover, by using Eqs (5), (8), (9) and (12), we have and ð14Þ Thus ½UðD; E; gÞ 2 ð1 À g À2 Þ À ðQÞ For the solution of nonlinear Eq (15), the Newton routine was employed to find its root γ. In our numerical computations, four to six iterations were required to obtain a root of tolerance 10 −6 . According to Eqs (14) and (15), the parameter γ is an implicit function of D, E and (Q) 2 in W m , m = 1, 2, 3, 4. Moreover, by using direct results from Eqs (5), (8), (9) and (11), i.e., It is seen that the physical variables are implicit functions of the conservative variables.
2 Numerical Schemes 2.1 One-dimensional central upwind scheme In this section, the semi discrete central-upwind scheme is derived [18]. The above SRHD model in one space dimension can be expressed as Before applying the scheme, it is necessary to discretize the computational domain. Let N denote the number of discretization points and ðx iÀ 1 2 Þ i2f1;ÁÁÁ;Nþ1g are partitions of the given interval [0, x max ]. For each i = 1, 2, Á Á Á, N, Δx is a constant width of each mesh interval, x i denote the cell centers, and x iAE 1 2 refer as cell boundaries. We assign, Moreover, By integrating Eq (17) over the interval The numerical fluxes are defined as Here, W + and W − are the point values of the piecewise linear reconstructionW ¼ ðr;ru;Ẽ Ã Þ for W, namely The numerical derivatives W x i are at least first-order approximations of W x (x i , t) and are computed using a nonlinear limiter that would ensure the non-oscillatory nature of reconstruction (23). A possible computation of these slopes is given by family of discrete derivatives parameterized with 1 θ 2, for example Where, Δ denotes central differencing and MM is the min-mod nonlinear limiter. Further, the local one sided speed at x iþ 1 2 is given as To obtain second order accuracy in time, we use a second order TVD Runge-Kutta scheme to solve Eq (21). Denoting the right-hand side of Eq (21) as L(W), a second order TVD Runge-Kutta scheme updates W through the following two stages where, W n is a solution at previous time step t n and W n+1 is updated solution at next time step t n+1 . Moreover, Δt represents the time step.

One-dimensional central scheme
We briefly give the high-resolution non-oscillatory central schemes of Nessyahu and Tadmor [20] for the solution of current SRHD model. These are the predictor-corrector type methods consist of two steps. In the predictor step, we use the non-oscillatory piecewise-linear reconstructions of the cell averages to predict the midpoint values. In the corrector step, staggered averaging along with the predicted mid-values are employed to get the updated cell averaged solution. The scheme can be summarized as Corrector : W nþ1 where, ξ = Δt/Δx. Moreover, 1 Dx F x ðW i Þ stands for an approximate numerical derivatives of the flux The fluxes F x (W i ) are computed by the same manner as for W x .

Two-dimensional central upwind scheme
Consider the two-dimensional SRHD equation (c.f. Eq (4)) Let N x and N y denote large integers in the x and y-directions, respectively. We assume a Cartesian mesh with a rectangular domain [x 0 , x max ] × [y 0 , y max ] which is covered by cells C ij The representative coordinates of the population in cell C ij are represented by (x i , y j ). Here, and The cell averaged values of conserved variables W i,j (t) at any time t is and the piecewise linear interplant is where, χ i,j is the characteristic function for the corresponding cell ðx iÀ 1 2 ; x iþ 1 2 Þ Â ðy jÀ 1 2 ; y jþ 1 2 Þ, (W x ) i,j and (W y ) i,j are the approximations of the x-and y-derivatives of W at the cell centers (x i , y j ). The generalized MM limiter is used for the computation of these partial derivatives to avoid oscillations where, 1 θ 2 and MM is defined above. After integrating the two-dimensional SRHD Eq (31) over the control volume C ij , the two-dimensional extension of the scheme can be expressed as Here, where, the intermediate values are defined as

Two-dimensional central scheme
This scheme was proposed by Jaing and Tadmor [25]. The scheme has a two-step predictorcorrector type. We start with the cell averages, W n i;j , the first-order predictor step is used to get midpoint values, W nþ 1 2 i;j , followed by the second-order corrector step to compute the new cell averages W nþ1 i;j . Similar to the 1D case, no exact (approximate) Riemann solver is required. The non-oscillatory property of the scheme depends upon the reconstructed discrete slopes, W x , W y , F x (w), and G y (w). At each time step, the grid is staggered to skip the flux calculation at the cell interfaces. The scheme can be written as below.
In the predictor step one has to calculate the mid-point values where, ξ = Δt/Δx and η = Δt/Δy. This step is followed by a corrector step to get the updated values at the next time step This completes the derivation of numerical schemes.

Numerical Tests
Here, some test problems are presented to validate our numerical scheme.

One-dimensional test problems
In this section, one-dimensional numerical test problems are considered. The results of proposed scheme are validated against those obtained from the central and KFVS schemes [20,24]. For the considered problems, the computational costs of all schemes were of the order of few seconds. Before going to the shock tube problems, we analyze the accuracy of schemes for smooth initial data in Problem 1.
The domain is taken to be 0 x 1, and the final simulation time t is 0.5. If h = Δx then L 1norm is defined as Here, α denotes its order, W is the exact solution and W h is the approximate solution. We de-  1. Moreover, all schemes have second order convergence rates. Shock tube problems. Here, the central upwind scheme is applied to simulate the discontinuous profiles of shock-tube problems. Given a numerical 1D pipe of unit length having N grid points, (ρ, v, p) and (ρ, v, p) being the two constant states on the left (0 < x 0.5) and on the right (0.5 < x 1), respectively. These states are with respect to a diaphragm. This diaphragm is placed initially at the middle of the pipe and then pull out. Shocks, contact discontinuities and rarefaction waves are the typical patterns which can be seen in the subsequent evolution. In the relativistic setup, these attributes does not change qualitatively, since the configuration of the characteristics is the same but density jumps could not be limited by a function of the adiabatic index and rarefaction waves do not produce straight profiles due to the nonlinear Lorentz transformation [4]. To validate our proposed central upwind scheme, the results are compared with those of high resolution central(NT) and KFVS methods [20,24] as well as with the exact Riemann solver [8,26].  The spatial domain is taken as [0, 1] with 500 grid points and t = 0.35 is the final time of simulation. The flow pattern is similar to that observed in Problem 2 but is more extreme now. Relativistic effects reduce the post-shock state to a thin dense shell which has width about 1% of the grid size at t = 0.35. The fluid velocity is v = 0.96c. The jump in the density is 10.6 for the exact Riemann solution while 3.9 in the case of the central upwind solution. A similar patterns were observed in [4] by using third order HHL and LLF solvers. The results are given in Fig 3. Here, the results of the central upwind and KFVS schemes seems better than the central(NT) scheme and are equivalent to those available in [4]. Problem 4: Perturbed relativistic shock tube flow. This problem was considered by [4,14]. The conditions at time t = 0 are specified as (ρ, v, p) L = (5, 0, 50) for 0 x 0.5 and (ρ, v, p) R = (ρ R , 0, 5) for 0.5 x 1.0. Here, a perturbed density field of a sinusoidal wave, ρ R = 2

Two-Dimensional Tests
Here, the 2D test problems are presented. The 2D numerical simulations are more complicated as compared to the 1D case. Here, both components of velocity are interpolated separately and, thus, the condition (v) 2 < 1 may not satisfy in the case of ultra-relativistic due to the effect of numerical errors.
Problem 5: EOC in two space dimensions. In this problem, we compare the EOCs of suggested numerical methods in two space dimensions. The initial data are taken as ; with s ¼ 0:13 and m ¼ 0 The domain is chosen to be 0 x 1, 0 y 1 and the total simulation time t is 0.5. Table 2 gives L 1 −errors and EOCs for the central upwind, the central (NT) and the KFVS methods. It can be observed that the central upwind scheme produces less errors in the solution as compared to the other schemes. Moreover, all schemes have second order convergence rates. Problem 6: Two-dimensional shock tube problem. Here, a 2D shock-tube problem is investigated [4]. A square domain of unit side length is taken which is divided into 4 quadrants.   The periphery of quadrants specify two 1D shocks and two contact discontinuities. These are symmetric with respect to the main diagonal. It is to be noted that we have not considered an exact 1D shock across the N and E interfaces, and this may be identified by investigating the emerged discontinuities in Fig 5, focusing toward the NE direction with their complete Riemann fan. In the remanning domain, the structure emerges with curved shock fronts and a Central Upwind Scheme for Special Relativistic Hydrodynamic Equations complex pattern in the SW quadrant, reminiscent of an oblique jet with converging internal shock fronts and a bow shock. The lines in the SW corner with respect to the bow shock are actually due to spurious waves produced at the S and W interfaces by the numerical diffusion term appearing in the energy equation, right and left states have jump in kinetic energy and cannot be eliminated even with a Roe-type solver [4]. The results are presented in Fig 5. Our results are identical to those available in [4]. Finally, Fig 6 gives a comparison of density and pressure obtained from the central upwind and the central(NT) schemes. It can be observed that the central upwind and KFVS schemes give more resolved solutions as compared to the central (NT) scheme. Problem 7: Diffracting shock waves and the forward facing step. In this problem, the shock diffraction at a corner is considered. This is a relativistic analog of the classical experimental results presented by van Dyke [27]. Here, the dimensions of the square region are less than those presented in [27]. The initial conditions are given in Fig 7. The results are presented in Fig 8. There is a shock wave placed at x = 0 and a rolled-up vortex is produced at the corner. A grid of 300 × 300 points is taken and the total time for simulation is 1.6. Finally, Figs 9 and 10 give a comparison of the central upwind, KFVS and central(NT) schemes. It can be observed that the central upwind and KFVS schemes gives more resolved solutions as compared to the central (NT) scheme. Moreover, the density plots of central (NT) scheme are oscillatory.
Problem 8: Cylindrical Explosion Problem. Here, a two-dimensional problem is examined. We consider a squared domain [0, 1] × [0, 1] which is divided into 200 × 200 mesh points. The initial data are constant in two regions separated by a circle of radius 0.2 with center (1/2, 1/2). Inside the circle ρ is 10.0 and p is 10, while outside the circle ρ is 1 and p is 0.01. The velocities are zero everywhere. In the solution, a circular shock wave propagates outwards from the origin, followed by circular contact discontinuity traveling in the same direction, and a circular rarefaction wave propagating towards the origin. Fig 11 shows the results of the central upwind scheme at final time t = 0.2. In Fig 12, the results of the central upwind scheme are compared with the central(NT) and KFVS schemes results. It can be observed that the results of all schemes agreeing well. However, the central upwind scheme has resolved solution than other two schemes.
Problem 9: Explosion in a box. This test problem corresponds to a 2D Riemann problem within a square box having unit side length. The walls are reflecting. The initial velocities are

Conclusions
In this article, the central upwind scheme was applied to solve special relativistic Euler equations in one and two space dimensions. The suggested method has capability to accurately capture the discontinues profiles of relativistic fluid and avoids numerical diffusion and dispersion in the solution. It was found that suggested scheme is second order accurate for smooth initial data in one and two space dimensions. The numerical errors in the solutions are less for the Central Upwind Scheme for Special Relativistic Hydrodynamic Equations central upwind scheme and are larger for the central (NT) scheme. In most of the considered problems, the numerical results of KFVS and central upwind schemes are comparable and central (NT) scheme produced diffusive results. However, in some problems, the results of central upwind scheme were found to be superior than the other two schemes. It was found that the computational costs of all schemes were of the order of few seconds in the one-dimensional case and of the order few minutes in the two-dimensional case.