Minimum Principle on Specific Entropy and High-Order Accurate Invariant Region Preserving Numerical Methods for Relativistic Hydrodynamics
MMinimum Principle on Specific Entropy and High-Order Accurate InvariantRegion Preserving Numerical Methods for Relativistic Hydrodynamics
Kailiang Wu * February 9, 2021
Abstract
This paper first explores Tadmor’s minimum entropy principle for the special relativistichydrodynamics (RHD) equations and incorporates this principle into the design of robust high-order discontinuous Galerkin (DG) and finite volume schemes for RHD on general meshes.The proposed schemes are rigorously proven to preserve numerical solutions in a global in-variant region constituted by all the known intrinsic constraints: minimum entropy principle,the subluminal constraint on fluid velocity, the positivity of pressure, and the positivity of rest-mass density. Relativistic effects lead to some essential difficulties in the present study, whichare not encountered in the non-relativistic case. Most notably, in the RHD case the specific en-tropy is a highly nonlinear implicit function of the conservative variables, and, moreover, thereis also no explicit formula of the flux in terms of the conservative variables. In order to over-come the resulting challenges, we first propose a novel equivalent form of the invariant region,by skillfully introducing two auxiliary variables. As a notable feature, all the constraints in thenovel form are explicit and linear with respect to the conservative variables. This provides ahighly effective approach to theoretically analyze the invariant-region-preserving (IRP) prop-erty of numerical schemes for RHD, without any assumption on the IRP property of the exactRiemann solver . Based on this, we prove the convexity of the invariant region and establishthe generalized Lax–Friedrichs splitting properties via technical estimates, lying the founda-tion for our rigorous IRP analyses. It is rigorously shown that the first-order Lax–Friedrichstype scheme for the RHD equations satisfies a local minimum entropy principle and is IRPunder a CFL condition. Provably IRP high-order accurate DG and finite volume methods arethen developed for the RHD with the help of a simple scaling limiter, which is designed byfollowing the bound-preserving type limiters in the literature. Several numerical examplesdemonstrate the effectiveness and robustness of the proposed schemes.
Keywords: minimum entropy principle, invariant region, bound-preserving, relativistic hy-drodynamics, discontinuous Galerkin, finite volume, high-order accuracy * Department of Mathematics, Southern University of Science and Technology, Shenzhen 518055, P.R. China.( [email protected] ). a r X i v : . [ m a t h . NA ] F e b Introduction
In the study of the fluid dynamics, when the fluid flow moves close to the speed of light or whensufficiently strong gravitational fields are involved, the special or general relativistic effect has tobe taken into account accordingly. Relativistic hydrodynamics (RHD) plays an important role in awide range of applications, such as astrophysics and high energy physics, and has been applied toinvestigate astrophysical scenarios from stellar to galactic scales.In the framework of special relativity, the motion of ideal relativistic fluid is governed by theconservation of mass density D , momentum m , and energy E . In the laboratory frame, the d –dimensional special RHD equations can be written into the following nonlinear hyperbolic system ∂ U ∂ t + d ∑ i = F i ( U ) ∂ x i = . (1)with the conservative vector U and the flux F i defined by U = ( D , mmm (cid:62) , E ) (cid:62) = ( ρ W , ρ hW vvv (cid:62) , ρ hW − p ) (cid:62) , (2) F i = ( Dv i , v i mmm (cid:62) + p e (cid:62) i , m i ) (cid:62) = ( ρ W v i , ρ hW v i vvv (cid:62) + p e (cid:62) i , ρ hW v i ) (cid:62) , (3)where and hereafter the geometrized unit system is employed so that the speed of light c = ρ denotes the rest-mass density, p is the thermal pressure, the column vector vvv = ( v , . . . , v d ) (cid:62) represents the velocity field of the fluid, W = / (cid:112) − (cid:107) vvv (cid:107) denotes the Lorentzfactor with (cid:107) · (cid:107) denoting the vector 2-norm, h = + e + p ρ stands for the specific enthalpy with e being the specific internal energy, and the row vector e i denotes the i -th column of the identitymatrix of size d . To close the system (1), an equation of state (EOS) is needed. We focus on theideal EOS: p = ( Γ − ) ρ e . (4)Here the constant Γ ∈ ( , ] denotes the adiabatic index, for which the restriction Γ ≤ U and the flux F i are explicitly expressedin terms of the primitive quantities V : = ( ρ , vvv (cid:62) , p ) (cid:62) . However, unlike the non-relativistic case,for RHD there are no explicit formulas for either the flux F i or the primitive vector V in terms of U . Therefore, in order to update the flux F i ( U ) in the computation, one has to first perform theinverse transformation of (2) from the conservative vector U to the primitive vector V . Given aconservative vector U = ( D , mmm (cid:62) , E ) (cid:62) , we can compute the values of the corresponding primitivequantities { p ( U ) , vvv ( U ) , ρ ( U ) } as follows: first solve a nonlinear algebraic equation (cid:107) mmm (cid:107) E + p + D (cid:115) − (cid:107) mmm (cid:107) ( E + p ) + p Γ − − E = , p ∈ [ , + ∞ ) , (5)2y utilizing certain root-finding algorithm to get the pressure p ( U ) ; then calculate the velocity andrest-mass density by vvv ( U ) = mmmE + p ( U ) , ρ ( U ) = D (cid:113) − (cid:107) vvv ( U ) (cid:107) . (6)Both the physical significance and the hyperbolicity of (1) require that the following constraints ρ > , p > , (cid:107) vvv (cid:107) < c = , (7)always hold. In other words, the conservative vector U must stay in the admissible state set G : = (cid:110) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : ρ ( U ) > , p ( U ) > , (cid:107) vvv ( U ) (cid:107) < (cid:111) , (8)where the functions ρ ( U ) , p ( U ) , and vvv ( U ) are highly nonlinear and have no explicit formulas, asdefined above. It was observed in [27] and rigorously proven in [38, Lemma 2.1] that the set G isconvex and is exactly equivalent to the following set G : = (cid:26) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : D > , E > (cid:113) D + (cid:107) mmm (cid:107) (cid:27) . (9)Moreover, if U ∈ G , then the nonlinear equation (5) has a unique positive solution [38].Due to the nonlinear hyperbolic nature of the equations (1), discontinuous solutions can de-velop from even smooth initial data, and weak solutions must therefore be considered. As well-known, weak solutions are not uniquely defined in general; the following inequality, known asthe entropy condition, is usually imposed as an admissibility criterion to select the “physicallyrelevant” solution among all weak solutions: ∂ E ∂ t + d ∑ i = ∂ F i ∂ x i ≤ , (10)which is interpreted in the sense of distribution. Here E ( U ) is a strictly convex function of U and called an entropy function, and F i ( U ) is the associated entropy fluxes such that the relation ∂ E ∂ U ∂ F i ∂ U = ∂ F i ∂ U holds. Entropy solutions are defined to be weak solutions which in addition satisfy(10) for all entropy pairs ( E , F i ) . For the (non-relativistic) gas dynamics equations, Tadmor [31]proved that entropy solutions satisfy a local minimum principle on the specific entropy S ( xxx , t ) = log (cid:0) p ρ − Γ (cid:1) : S ( xxx , t + τ ) ≥ min { S ( yyy , t ) : (cid:107) yyy − xxx (cid:107) ≤ τ v max } , (11)where τ > v max denotes the maximal wave speed. This implies the spatial minimum ofthe specific entropy, min xxx S ( xxx , t ) , is a nondecreasing function of time t , and S ( xxx , t ) ≥ min xxx S ( xxx , ) .(Entropy principles were also shown by Guermond and Popov [9] with viscous regularizationof the non-relativistic Euler equations.) In this paper, we will explore such a minimum entropy3rinciple in the RHD case and prove that it also holds for entropy solutions of (1). The entropyprinciple, along with the intrinsic physical constraints in (8), imply a global invariant region forthe solution of the RHD equations (1) with initial data U ( xxx ) , that is, Ω S : = (cid:110) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : ρ ( U ) > , p ( U ) > , (cid:107) vvv ( U ) (cid:107) < , S ( U ) ≥ S (cid:111) , (12)where S : = ess inf xxx S ( U ( xxx )) .It is natural and interesting to explore robust numerical RHD schemes, whose solutions alwaysstay in the invariant region Ω S , i.e., satisfy the minimum entropy principle at the discrete leveland also preserve the intrinsic physical constraints (7). Note that, to obtain a well-defined spe-cific entropy for the numerical solution, it is necessary to first guarantee the positivity of pressureand rest-mass density. The subluminal constraint on the fluid velocity is also crucial for the rel-ativistic causality, because its violation would yield imaginary Lorentz factor. In fact, violatingany of the constraints (7) will cause numerical instability and the break down of the computa-tion. Therefore, the preservation of the minimum entropy principle should be considered togetherwith the constraints (7). Recent years have witnessed some advances in developing high-ordernumerical schemes, which provably preserve the constraints (7), for the special RHD [38, 28, 41]and the general RHD [32]. Those works were motivated by [44, 45, 42, 16] on designing bound-preserving high-order schemes for scalar conservation laws and the non-relativistic Euler equa-tions. More recently, bound-preserving numerical schemes were also developed for the specialrelativistic magnetohydrodynamics (RMHD) in one and multiple dimensions [40, 37], as exten-sion of the positivity-preserving MHD schemes [33, 34, 35]. In addition, a flux-limiting approachwhich preserves the positivity of the rest-mass density was designed in [29]. A reconstructiontechnique was proposed in [1] to enforce the subluminal constraint on the fluid velocity. A flux-limiting entropy-viscosity approach was developed in [7] for RHD, based on a measure of theentropy generated by the solution. Systematic review of numerical RHD schemes is beyond thescope of the present paper; we refer interested readers to the review articles [24, 4, 25] and a lim-ited list of some recent works [15, 48, 39, 3, 2, 36, 26] as well as references therein. Yet, up tonow, there is still no work that studied the minimum entropy principle for the RHD equations (1)at either the PDE level or the numerical level, and high-order schemes which provably preservethe invariant region (12) have not yet been developed for RHD.For the non-relativistic counterparts such as the compressible Euler system, the minimum en-tropy principle and invariant-region-preserving (IRP) numerical schemes have been well studied inthe literature. Tadmor [31] discovered (11) and proved, for the compressible Euler equations of thegas dynamics, that first-order approximations such as the Godunov and Lax-Friedrichs schemessatisfy a minimum entropy principle. Using a slope reconstruction with limiter, Khobalatte andPerthame [21] developed second-order kinetic schemes that preserve a discrete minimum principlefor the specific entropy. It was also observed in [21] that enforcing the discrete minimum entropy4rinciple could help to damp numerical oscillations near the discontinuities. Zhang and Shu [46]proposed a framework of enforcing the minimum entropy principle for high-order accurate finitevolume and discontinuous Galerkin (DG) schemes, by extending their positivity-preserving high-order schemes [45, 47], for the (non-relativistic) gas dynamics equations. The resulting high-orderschemes in [46] were proven to preserve a discrete minimum entropy principle and the positiv-ity of density and pressure, under a condition accessible by a simple bound-preserving limiterwithout destroying the high-order accuracy. Lv and Ihme [23] proposed an entropy-bounded DGscheme for the Euler equations on arbitrary meshes. Guermond, Popov, and their collaborators(cf. [10, 11, 8, 12, 13]) developed the IRP approximations in the context of continuous finite el-ements with convex limiting for solving general hyperbolic systems including the compressibleEuler equations. Jiang and Liu proposed new IRP limiters for the DG schemes to the isentropicEuler equations [20], the compressible Euler equations [17], and general multi-dimensional hyper-bolic conservation laws [18]. Gouasmi et al. [6] proved a minimum entropy principle on entropysolutions to the compressible multicomponent Euler equations at the smooth and discrete levels.The aim of this paper is twofold. The first is to show that the minimum entropy principle (11),which was originally demonstrated by Tadmor [31] for the (non-relativistic) gas dynamics, is alsovalid for the RHD equations (1) with the ideal EOS (4). A key point in the present study is toprove a condition on smooth function H ( S ) such that the entropy function E ( U ) = − D H ( S ) is strictly convex. The second goal is to develop high-order accurate IRP DG and finite volumemethods which provably preserve the numerical solutions in the invariant region Ω S , i.e., preservea discrete minimum entropy principle and the intrinsic physical constraints (7). In fact, achievingthese two goals is nontrivial. Due to the nonlinearity and the implicit form of the function S ( U ) , itis not easy to study the convexity of the entropy function E ( U ) in the RHD case; see Proposition2.1. Also, analytically judging whether an arbitrarily given state U belongs to Ω S is already a dif-ficult task; it is more challenging to design and analyze numerical schemes that provably preservethe solutions in Ω S . We will address the difficulties via a novel equivalent form of the invariantregion; see Theorem 3.1. As a notable feature, all the constraints in the novel form are explicit and linear with respect to the conservative variables. This provides a highly effective approach totheoretically analyze the IRP property of RHD schemes. Based on this, we will prove the con-vexity of the invariant region (Section 3.2) and establish the generalized Lax–Friedrichs splittingproperties via highly technical estimates (Section 3.3), which lie the foundation for analyzing ourIRP schemes in Sections 4–5. The high-order accurate IRP schemes are constructed with the aidof a simple scaling limiter, which is designed by following the bound-preserving type limiters andframeworks in the literature [45, 46, 28, 18, 41].It is worth noting that the proposed IRP analysis approach has some essential and significantdifferences from those in the literature (cf. [46, 10, 8, 19]). For example, some standard analyseswere often based on the IRP property of the exact Riemann solver for the studied equations, while5ur IRP analyses do not require any assumption on the IRP property of the exact (or any approx-imate) Riemann solver. This makes our analysis approach potentially extensible to some othercomplicated physical systems for which the exact Riemann solver is not easily available.This paper is organized as follows. We study the minimum entropy principle at the PDE levelin Section 2. After establishing some auxiliary theories for IRP analysis in Section 3, we presentthe IRP schemes in Sections 4–5 for one- and multi-dimensional RHD equations, respectively.Numerical examples are provided in Section 6 and will confirm that incorporating the minimumentropy principle into a scheme could be helpful for damping some undesirable numerical oscilla-tions, as observed in e.g. [21, 46, 19] for some other systems. Section 7 concludes the paper.
We first explore the space of admissible entropy functions. Let H ( S ) be a function of the specificentropy S . As shown in [2] for d =
1, for any smooth function H ( S ) , the smooth solutions of theRHD equations (1) satisfy ∂∂ t (cid:16) − ρ W H ( S ) (cid:17) + d ∑ i = ∂∂ x i (cid:16) − ρ W v i H ( S ) (cid:17) = . (13)This implies that ( E , F i ) = (cid:0) − D H ( S ) , − Dv i H ( S ) (cid:1) is an entropy–entropy flux pair, if E ( U ) = − D H ( S ) is a strictly convex function of the conservative variables U ∈ G .It is well-known that, for a special choice H ( S ) = S or H ( S ) = S Γ − , the corresponding E ( U ) is a valid entropy function for the RHD; see, for example, [7, 2, 3, 36]. However, it is unclear,for a general H ( S ) , what is the condition on H ( S ) such that the corresponding E ( U ) is strictlyconvex. This has not been addressed for the RHD case in the literature and is now explored in thefollowing proposition. Proposition 2.1.
For a smooth function H ( S ) , the corresponding E ( U ) = − D H ( S ) is a strictlyconvex entropy function if and only if H (cid:48) ( S ) > , H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) > . (14) Proof.
We study the convexity of E ( U ) by investigating the positive definiteness of the associatedHessian matrix E uu : = (cid:18) ∂ E ∂ u i ∂ u j (cid:19) ≤ i , j ≤ d + , It is certainly reasonable to assume that the exact Riemann solver preserves the invariant domain. In fact, this isprovenly true for a number of other systems, but has not yet been proven for the RHD equations (1). Rigorous analysison the IRP property of the exact Riemann solver, including the preservation of the constraints (7), is highly nontrivial. u i denotes the i th component of U . A straightforward computation gives ∂ E ∂ u i ∂ u j = − H (cid:48) ( S ) (cid:18) ∂ D ∂ u i ∂ S ∂ u j + ∂ D ∂ u j ∂ S ∂ u i (cid:19) − D H (cid:48)(cid:48) ( S ) ∂ S ∂ u i ∂ S ∂ u j − D H (cid:48) ( S ) ∂ S ∂ u i ∂ u j , which implies that E uu = − H (cid:48) ( S ) (cid:16) e S (cid:62) u + S u e (cid:62) + DS uu (cid:17) − D H (cid:48)(cid:48) ( S ) S u S (cid:62) u = − H (cid:48) ( S ) A + D Γ (cid:0) H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) (cid:1) S u S (cid:62) u , (15)where e = ( , (cid:62) d + ) (cid:62) , d + denotes the zero vector of length d +
1, and A : = e S (cid:62) u + S u e (cid:62) + DS uu + Γ DS u S (cid:62) u . Since S cannot be explicitly formulated in terms of U , direct derivation of S u and S uu is difficult.Let us consider the primitive variables V = ( ρ , vvv (cid:62) , p ) (cid:62) . Note that both S and U can be explicitlyformulated in terms of V , then it is easy to derive ∂ S ∂ V = (cid:16) − Γ / ρ , (cid:62) d , / p (cid:17) , ∂ U ∂ V = W ρ W vvv (cid:62) W vvv ρ hW I d + ρ hW vvvvvv (cid:62) Γ W Γ − vvvW ρ hW vvv (cid:62) Γ W Γ − − , where I d denotes the identity matrix of size d . The inverse of the matrix ∂ U ∂ V gives ∂ V ∂ U = ρ h ( − c s (cid:107) vvv (cid:107) ) ρ h ( − ( Γ − ) (cid:107) vvv (cid:107) ) W − − ρ ( + ( Γ − ) (cid:107) vvv (cid:107) ) vvv (cid:62) ρ Γ (cid:107) vvv (cid:107) ( Γ − ) W − vvv A Γ (cid:0) (cid:107) vvv (cid:107) − (cid:1) vvv − ( Γ p + ( Γ − ) ρ ) W − − ( Γ p + ( Γ − ) ρ ) vvv (cid:62) Γ p (cid:0) + (cid:107) vvv (cid:107) (cid:1) + ( Γ − ) ρ with c s = (cid:113) Γ p ρ h denoting the acoustic wave speed in the RHD case (note that 0 < c s < A : = (cid:0) − (cid:107) vvv (cid:107) (cid:1) (cid:104)(cid:0) − c s (cid:107) vvv (cid:107) (cid:1) I d + (cid:0) Γ − + c s (cid:1) vvvvvv (cid:62) (cid:105) . It follows that S (cid:62) u = ∂ S ∂ U = ∂ S ∂ V ∂ V ∂ U = Γ − p (cid:16) − hW − , − vvv (cid:62) , (cid:17) . The derivative of S (cid:62) u with respect to V gives S uv = Γ ρ (cid:112) − (cid:107) vvv (cid:107) Γ − p hW vvv (cid:62) Γ − p (cid:112) − (cid:107) vvv (cid:107) d − Γ − p I d Γ − p vvv (cid:62) d − Γ − p . A = e S (cid:62) u + S u e (cid:62) + DS uv ∂ V ∂ U + Γ DS u S (cid:62) u = − Γ ph ( h − )( − c s (cid:107) vvv (cid:107) ) a a vvv (cid:62) a a vvv A a vvva a vvv (cid:62) a with a : = h ( Γ − ) W − > , a : = ( h − )( Γ − ) , a : = − ( Γ − )( h + ( h − ) (cid:107) vvv (cid:107) ) , A : = ( h − )( − c s (cid:107) vvv (cid:107) ) W I d + W (cid:18) ( h − )( − c s (cid:107) vvv (cid:107) ) + h ( Γ − )( h − ) (cid:19) vvvvvv (cid:62) , a : = W (cid:0) h ( − c s (cid:107) vvv (cid:107) ) − Γ ( h − ) (cid:1) , a : = W (cid:0) ( h − )( Γ − ) (cid:107) vvv (cid:107) + ( Γ − ) h (cid:1) . Let us define the invertible matrix P = (cid:62) d − a a vvv I d d − a a (cid:62) d . A straightforward computation gives P A P (cid:62) = − Γ ph ( h − )( − c s (cid:107) vvv (cid:107) ) (cid:32) a (cid:62) d + d + a A (cid:33) (16)with a : = ( h − ) W (cid:0) − c s (cid:107) vvv (cid:107) (cid:1) >
0, and A : = (cid:32) ( − (cid:107) vvv (cid:107) ) I d + vvvvvv (cid:62) − vvv − vvv (cid:62) (cid:107) vvv (cid:107) (cid:33) . (17)Note that P S u = Γ − p (cid:16) − hW − , ( h − ) vvv (cid:62) , ( − h )( + (cid:107) vvv (cid:107) ) (cid:17) (cid:62) = : Γ − p bbb . (18)Combining equations (15), (16) and (18) gives P E uu P (cid:62) = a H (cid:48) ( S ) A + a (cid:0) H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) (cid:1) bbb bbb (cid:62) (19)with a : = Γ − ph ( h − )( − c s (cid:107) vvv (cid:107) ) > a : = D ( Γ − ) p Γ >
0, and A : = (cid:32) a (cid:62) d + d + a A (cid:33) . A defined in (17). The matrix ( − (cid:107) vvv (cid:107) ) I d + vvvvvv (cid:62) is symmetric andits eigenvalues consist of 1 and 1 − (cid:107) vvv (cid:107) , which are all positive, implying that ( − (cid:107) vvv (cid:107) ) I d + vvvvvv (cid:62) is positive definite. Furthermore, a straightforward calculation shows det ( A ) =
0. Therefore, A is positive semi-definite, and rank ( A ) = d . Since a > a >
0, it follows that A is positivesemi-definite, and rank ( A ) = d +
1. Hence, there exists a rank- ( d + ) matrix A ∈ R ( d + ) × ( d + ) such that A (cid:62) A = A . (20)Because E uu and P E uu P (cid:62) are congruent, it suffices to prove that the matrix P E uu P (cid:62) is posi-tive definite if and only if H ( S ) satisfies the condition (14). (i) . First prove the condition (14) is sufficient for the positive definiteness of P E uu P (cid:62) . Be-cause A and bbb bbb (cid:62) are both positive semi-definite, by (19) we know that P E uu P (cid:62) is positivesemi-definite under the condition (14). It means zzz (cid:62) P E uu P (cid:62) zzz ≥ , ∀ zzz ∈ R d + . (21)Hence, it suffices to show zzz = when zzz (cid:62) P E uu P (cid:62) zzz =
0. Using (19) and (20), we have zzz (cid:62) P E uu P (cid:62) zzz = a H (cid:48) ( S ) (cid:107) A zzz (cid:107) + a (cid:0) H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) (cid:1) (cid:12)(cid:12)(cid:12) bbb (cid:62) zzz (cid:12)(cid:12)(cid:12) = , which implies A zzz = d + and bbb (cid:62) zzz =
0. Then A zzz = A (cid:62) A zzz = . Let zzz = : ( z ( ) , zzz ( ) , z ( ) ) (cid:62) with zzz ( ) ∈ R d . From A zzz = we can deduce that a z ( ) = a A ( zzz ( ) , z ( ) ) (cid:62) = . It further yields z ( ) = ( − (cid:107) vvv (cid:107) ) zzz ( ) + vvvvvv (cid:62) zzz ( ) − z ( ) vvv = d (22) − vvv (cid:62) zzz ( ) + (cid:107) vvv (cid:107) z ( ) = . (23)Combining z ( ) = bbb (cid:62) zzz = ( h − ) vvv (cid:62) zzz ( ) + ( − h )( + (cid:107) vvv (cid:107) ) z ( ) = , which, together with (23), imply vvv (cid:62) zzz ( ) = z ( ) =
0. Substituting it into (22) gives zzz ( ) = d .Therefore, we have zzz = when zzz (cid:62) P E uu P (cid:62) zzz =
0. This along with (21) yield that P E uu P (cid:62) ispositive definite under the condition (14). This completes the proof of sufficiency. (ii) . Then prove the condition (14) is necessary for the positive definiteness of P E uu P (cid:62) .Assume that P E uu P (cid:62) is positive definite, then zzz (cid:62) P E uu P (cid:62) zzz = a H (cid:48) ( S ) zzz (cid:62) A zzz + a (cid:0) H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) (cid:1) (cid:12)(cid:12)(cid:12) bbb (cid:62) zzz (cid:12)(cid:12)(cid:12) > , ∀ zzz ∈ R d + \ { } . (24)Note that the matrix A does not have full rank. There exist two vectors zzz , zzz ∈ R d + \ { } suchthat A zzz = and bbb (cid:62) zzz =
0, respectively. It follows from (24) that0 < zzz (cid:62) P E uu P (cid:62) zzz = a H (cid:48) ( S ) zzz (cid:62) A zzz = a H (cid:48) ( S ) (cid:107) A zzz (cid:107) , < zzz (cid:62) P E uu P (cid:62) zzz = a (cid:0) H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) (cid:1) (cid:12)(cid:12)(cid:12) bbb (cid:62) zzz (cid:12)(cid:12)(cid:12) , H (cid:48) ( S ) > H (cid:48) ( S ) − Γ H (cid:48)(cid:48) ( S ) >
0, respectively. This completes the proof ofnecessity. (cid:4)
Remark 2.1.
Proposition 2.1 and equation (13) imply that there exists a family of (generalized)entropy pairs ( E , F i ) associated with the d-dimensional ( ≤ d ≤ ) RHD equations (1) , E ( U ) = − D H ( S ) , F i ( U ) = − Dv i H ( S ) , i = , . . . , d , (25) generated by the smooth functions H ( S ) satisfying (14) . Our found condition (14) is consistentwith the one derived by Harten [14, Section 2] for the 2D non-relativistic Euler equations. How-ever, the analysis in the RHD case does not directly follow from [14] and is more difficult. Due tothe complicated structures of the matrix E uu , some standard approaches for investigating its posi-tive definiteness, e.g., checking the positivity of its leading principal minors, can be intractable inour RHD case. We are now in a position to verify that Tadmor’s minimum entropy principle (11) also hold for theRHD system (1). We consider the convex entropy E ( U ) = − D H ( S ) established in Section 2.1for all smooth functions H ( S ) satisfying (14).Assume that U ( xxx , t ) is an entropy solution of the RHD equations (1). According to [30, The-orem 4.1] and following [31, Lemma 3.1], we have, for all smooth functions H ( S ) satisfying(14), (cid:90) (cid:107) xxx − xxx (cid:107)≤ R D ( xxx , t + τ ) H ( S ( xxx , t + τ )) d xxx ≥ (cid:90) (cid:107) xxx − xxx (cid:107)≤ R + τ v max D ( xxx , t ) H ( S ( xxx , t )) d xxx , ∀ R > , (26)where v max denotes the maximal wave speed in the domain; we can take v max as the speed of light c =
1, a simple upper bound of all wave speeds in the RHD case. Note that the density involved in(26) is D = ρ W , instead of the rest-mass density ρ .Consider a special function H ( S ) [31] defined by H ( S ) : = min { S − S , } , S = min { S ( yyy , t ) : (cid:107) yyy − xxx (cid:107) ≤ R + τ v max } . As observed in [31], the function H ( S ) , although not smooth, can be written as the limit ofa sequence of smooth functions satisfying (14); see also [6, Section 3.1] for a detailed review.Therefore, by passing to the limit, the inequality (26) holds for H = H , which gives (cid:90) (cid:107) xxx − xxx (cid:107)≤ R D ( xxx , t + τ ) min { S ( xxx , t + τ ) − S , } d xxx ≥ (cid:90) (cid:107) xxx − xxx (cid:107)≤ R + τ v max D ( xxx , t ) H ( S ( xxx , t )) d xxx = . Because D ( xxx , t + τ ) >
0, we obtain min { S ( xxx , t + τ ) − S , } = (cid:107) xxx − xxx (cid:107) ≤ R . It leads to S ( xxx , t + τ ) ≥ S = min { S ( yyy , t ) : (cid:107) yyy − xxx (cid:107) ≤ R + τ v max } , ∀(cid:107) xxx − xxx (cid:107) ≤ R , (27)10hich yields the local minimum entropy principle (11). In particular, it implies that the spatialminimum of the specific entropy, min xxx S ( xxx , t ) , is a nondecreasing function of time t , yielding S ( xxx , t ) ≥ min xxx S ( xxx , ) , ∀ t ≥ . (28)In the above derivation it is implicitly assumed that U ( xxx , t ) always satisfies the physical con-straints (7). The entropy principle (28) and the constraints (7) constitute the global invariant region Ω S , defined in (12), for entropy solutions of the RHD equations (1). In order to analyze the local minimum entropy principle of numerical schemes, we introduce a(more general) “local” invariant region for an arbitrarily given σ : Ω σ : = (cid:110) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : ρ ( U ) > , p ( U ) > , (cid:107) vvv ( U ) (cid:107) < , S ( U ) ≥ σ (cid:111) . (29)The special choice σ = S = ess inf xxx S ( xxx , ) corresponds to the global invariant region Ω S in (12).It is evident that the following “monotonicity” holds for Ω σ . Lemma 3.1 ( Monotonic decreasing). If σ ≥ σ , then Ω σ ⊆ Ω σ . Thanks to G = G proved in [38, Lemma 2.1], we immediately have Lemma 3.2 ( First equivalent form).
The invariant region set Ω σ is equivalent to Ω ( ) σ = (cid:26) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : D > , E > (cid:113) D + (cid:107) mmm (cid:107) , S ( U ) ≥ σ (cid:27) . (30)Note that the specific entropy S = log ( p ρ − Γ ) is a nonlinear function of ( ρ , p ) , and, as men-tioned in Section 1, the functions ρ ( U ) and p ( U ) are already highly nonlinear and without explicitformulas. The combination of these nonlinear functions leads to S ( U ) , which is certainly a highlynonlinear function and also cannot be explicitly formulated in terms of U . This causes it difficultto study the minimum entropy principle at the numerical level and explore IRP schemes for RHD.In order to overcome the challenges, several important properties of the invariant region Ω σ willbe derived in this section. To address the difficulties arising from the nonlinearity of S ( U ) , we discover the following novelequivalent form of Ω σ . 11 heorem 3.1 ( Second equivalent form).
The invariant region set Ω σ is equivalent to Ω ( ) σ = (cid:110) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : D > , ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ , ∀ vvv ∗ ∈ B ( ) , ∀ ρ ∗ ∈ R + (cid:111) , (31) where B ( ) : = { xxx ∈ R d : | xxx | < } denotes the open unit ball centered at in R d , and ϕ σ ( U ; vvv ∗ , ρ ∗ ) : = E − mmm · vvv ∗ − D (cid:113) − (cid:107) vvv ∗ (cid:107) + e σ (cid:18) ρ Γ ∗ − ΓΓ − D ρ Γ − ∗ (cid:113) − (cid:107) vvv ∗ (cid:107) (cid:19) . (32)Before the proof of Theorem 3.1, we mention a crucial feature of the above equivalent form Ω ( ) σ . Note that all the “nonlinear” constraints in Ω ( ) σ or Ω σ are equivalently transformed into a linear constraint ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ Ω ( ) σ are explicit and linear with respect to U , although two (additional) auxiliary variables vvv ∗ and ρ ∗ are introducedhere. Such linearity makes Ω ( ) σ very useful for analytically verifying the IRP property of RHDschemes. This becomes a key to our IRP analysis, which is significantly different from the standardbound-preserving and IRP analysis techniques in e.g. [45, 46, 18].We first give two lemmas, which will be used in the proof of Theorem 3.1. Lemma 3.3.
For any η > − , it holds ( η Γ + ) Γ ≥ ( η + ) , where the (constant) adiabatic index Γ ∈ ( , ] .Proof. Consider the function f ( x ) = ( η x + ) x with x ∈ ( , ] . Note that η x + > f ( x ) > x satisfies x f ( x ) f (cid:48) ( x ) = − log ( η x + ) + η x η x + = log (cid:18) + − η x η x + (cid:19) + η x η x + ≤ , where we have used the elementary inequality log ( + z ) ≤ z for z > −
1. It follows that f (cid:48) ( x ) ≤ x ∈ ( , ] . This implies f ( Γ ) = ( η Γ + ) Γ = f ( ) − (cid:82) Γ f (cid:48) ( x ) d x ≥ f ( ) = ( η + ) . (cid:4) Lemma 3.4.
For any vvv , vvv ∗ ∈ B ( ) , it holds Γ ( − vvv · vvv ∗ ) − (cid:107) vvv (cid:107) − Γ + ≥ (cid:32) (cid:112) − (cid:107) vvv (cid:107) (cid:112) − (cid:107) vvv ∗ (cid:107) (cid:33) − Γ . (33) Proof.
Note that η : = − vvv · vvv ∗ − (cid:107) vvv (cid:107) − ≥ − (cid:107) vvv (cid:107) − (cid:107) vvv (cid:107) − = + (cid:107) vvv (cid:107) − > − . With the aid of Lemma 3.3, we get ( η Γ + ) Γ ≥ √ η + , which yields (cid:18) Γ ( − vvv · vvv ∗ ) − (cid:107) vvv (cid:107) − Γ + (cid:19) Γ ≥ (cid:115) (cid:18) − vvv · vvv ∗ − (cid:107) vvv (cid:107) − (cid:19) + ≥ (cid:112) − (cid:107) vvv ∗ (cid:107) (cid:112) − (cid:107) vvv (cid:107) . Then, by raising both sides to the power of Γ , we obtain (33). (cid:4)
12e are now ready to prove Theorem 3.1.
Proof of Theorem 3.1 . The proof of Ω σ = Ω ( ) σ is split into two parts — showing that Ω ( ) σ ⊆ Ω σ and that Ω σ ⊆ Ω ( ) σ . (i). Prove that U ∈ Ω ( ) σ ⇒ U ∈ Ω σ . When U = ( D , mmm (cid:62) , E ) (cid:62) ∈ Ω ( ) σ , by definition wehave D > ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ , ∀ vvv ∗ ∈ B ( ) , ∀ ρ ∗ ∈ R + . (34)If we take a special ( vvv ∗ , ρ ∗ ) = (cid:16) mmm √ D + (cid:107) mmm (cid:107) , D √ D + (cid:107) mmm (cid:107) (cid:17) , which satisfy vvv ∗ ∈ B ( ) and ρ ∗ >
0, andthen we obtain 0 ≤ ϕ σ (cid:32) U ; mmm (cid:112) D + (cid:107) mmm (cid:107) , D (cid:112) D + (cid:107) mmm (cid:107) (cid:33) = E − (cid:113) D + (cid:107) mmm (cid:107) − e σ Γ − (cid:32) D (cid:112) D + (cid:107) mmm (cid:107) (cid:33) Γ < E − (cid:113) D + (cid:107) mmm (cid:107) , which implies the second constraint in G . This, along with D >
0, yield U ∈ G = G . Therefore,the corresponding primitive quantities of U satisfy ρ ( U ) > , p ( U ) > , (cid:107) vvv ( U ) (cid:107) < . (35)Taking another special ( vvv ∗ , ρ ∗ ) = ( vvv ( U ) , ρ ( U )) , in (34) gives0 ≤ ϕ σ ( U ; vvv ( U ) , ρ ( U ))= E − mmm · vvv − D (cid:113) − (cid:107) vvv (cid:107) + e σ (cid:18) ρ Γ − ΓΓ − D ρ Γ − (cid:113) − (cid:107) vvv (cid:107) (cid:19) = Γ − ( p − e σ ρ Γ ) , which, together with Γ >
1, imply p ≥ e σ ρ Γ . It follows that S ( U ) = log ( p ρ − Γ ) ≥ σ . Combiningit with (35), we obtain U ∈ Ω σ . (ii). Prove that U ∈ Ω σ ⇒ U ∈ Ω ( ) σ . When U = ( D , mmm (cid:62) , E ) (cid:62) ∈ Ω σ , the correspondingprimitive quantities satisfy ρ > , (cid:107) vvv (cid:107) < , p ≥ e σ ρ Γ . (36)This immediately gives D = ρ W = ρ (cid:0) − (cid:107) vvv (cid:107) (cid:1) − > . It remains to prove ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ vvv ∗ ∈ B ( ) and any ρ ∗ >
0. Let us rewrite ϕ σ ( U ; vvv ∗ , ρ ∗ ) as ϕ σ ( U ; vvv ∗ , ρ ∗ ) = Π p + Π , Π : = ΓΓ − (cid:18) − vvv · vvv ∗ − (cid:107) vvv (cid:107) (cid:19) − ≥ (cid:18) − vvv · vvv ∗ − (cid:107) vvv (cid:107) (cid:19) − = ( − (cid:107) vvv (cid:107) ) + ( (cid:107) vvv (cid:107) − vvv · vvv ∗ ) − (cid:107) vvv (cid:107) > , Π : = ρ (cid:32) − vvv · vvv ∗ − (cid:107) vvv (cid:107) − (cid:112) − (cid:107) vvv ∗ (cid:107) (cid:112) − (cid:107) vvv (cid:107) (cid:33) + e σ (cid:32) ρ Γ ∗ − ΓΓ − ρρ Γ − ∗ (cid:112) − (cid:107) vvv ∗ (cid:107) (cid:112) − (cid:107) vvv (cid:107) (cid:33) , where Γ ∈ ( , ] and (cid:107) vvv (cid:107) < Π >
0. Then, using p ≥ e σ ρ Γ gives ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ Π e σ ρ Γ + Π = ρ − (cid:107) vvv (cid:107) (cid:18) − vvv · vvv ∗ − (cid:113) − (cid:107) vvv (cid:107) (cid:113) − (cid:107) vvv ∗ (cid:107) (cid:19) + e σ ρ Γ Π ≥ ρ − (cid:107) vvv (cid:107) (cid:18) − (cid:113) (cid:107) vvv (cid:107) + ( − (cid:107) vvv (cid:107) ) (cid:113) (cid:107) vvv ∗ (cid:107) + ( − (cid:107) vvv ∗ (cid:107) ) (cid:19) + e σ ρ Γ Π = e σ ρ Γ Π , (37)where the Cauchy–Schwarz inequality has been used, and Π : = φ σ (cid:0) ρρ ∗ ; vvv , vvv ∗ (cid:1) with φ σ ( x ; vvv , vvv ∗ ) : = x − Γ − Γ (cid:112) − (cid:107) vvv ∗ (cid:107) ( Γ − ) (cid:112) − (cid:107) vvv (cid:107) x − Γ + + Γ ( − vvv · vvv ∗ )( Γ − ) ( − (cid:107) vvv (cid:107) ) − , x > . It is easy to verify that the function φ σ ( x ; vvv , vvv ∗ ) is strictly decreasing on the interval (cid:18) , √ −(cid:107) vvv (cid:107) √ −(cid:107) vvv ∗ (cid:107) (cid:21) and strictly increasing on (cid:20) √ −(cid:107) vvv (cid:107) √ −(cid:107) vvv ∗ (cid:107) , + ∞ (cid:19) with respect to x . Thefore, we have Π ≥ min x ∈ R + φ σ ( x ; vvv , vvv ∗ ) = φ σ (cid:32) (cid:112) − (cid:107) vvv (cid:107) (cid:112) − (cid:107) vvv ∗ (cid:107) ; vvv , vvv ∗ (cid:33) = Γ − − (cid:32) (cid:112) − (cid:107) vvv (cid:107) (cid:112) − (cid:107) vvv ∗ (cid:107) (cid:33) − Γ + Γ ( − vvv · vvv ∗ ) − (cid:107) vvv (cid:107) − Γ + ≥ , (38)where we have used Γ ∈ ( , ] and the inequality (33) derived in Lemma 3.4. Then, combining(38) with (37), we conclude ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ e σ ρ Γ Π ≥
0. This along with D > U ∈ Ω ( ) σ .The proof is completed. (cid:4) Considering σ → − ∞ in Ω σ and using Theorem 3.1, one can also obtain: Corollary 3.1.
The admissible state set G is equivalent to G : = (cid:26) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : D > , E > mmm · vvv ∗ + D (cid:113) − (cid:107) vvv ∗ (cid:107) , ∀ vvv ∗ ∈ B ( ) (cid:27) . .2 Convexity of invariant region The convexity of an invariant region is a highly desirable property, as it can be used to simplifythe IRP analysis of those numerical schemes that can be reformulated into some suitable convexcombinations; see, for example, [45, 46, 33, 18]. In the RHD case, the convexity of the invariantregion Ω σ is discussed below. Lemma 3.5.
For any fixed σ ∈ R , the invariant region Ω σ is a convex set.Proof. Since Ω σ = Ω ( ) σ , one only needs to show the convexity of Ω ( ) σ . For any U = ( D , mmm (cid:62) , E ) (cid:62) , U = ( D , mmm (cid:62) , E ) (cid:62) ∈ Ω ( ) σ and any λ ∈ [ , ] , we have λ D + ( − λ ) D > vvv ∗ ∈ B ( ) , ρ ∗ ∈ R + , it holds ϕ σ (cid:0) λ U + ( − λ ) U ; vvv ∗ , ρ ∗ (cid:1) = λ ϕ σ ( U ; vvv ∗ , ρ ∗ ) + ( − λ ) ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ , where the linearity of ϕ σ ( U ; vvv ∗ , ρ ∗ ) with respect to U has been used. Hence, we obtain λ U + ( − λ ) U ∈ Ω ( ) σ , and by definition, Ω ( ) σ is a convex set. (cid:4) The following more general conclusion can be shown by using Lemma 3.5 and Lemma 3.1.
Lemma 3.6.
Let σ and σ be two real numbers. For any U ∈ Ω σ and any U ∈ Ω σ , λ U + ( − λ ) U ∈ Ω min { σ , σ } , ∀ λ ∈ [ , ] . Proof.
With the help of Lemma 3.1, we have U i ∈ Ω σ i ⊆ Ω min { σ , σ } , i = ,
2. The proof is thencompleted by using the convexity of Ω min { σ , σ } . (cid:4) Remark 3.1.
An alternative proof of Lemma 3.6 or Lemma 3.5 is based on the first equivalent set Ω ( ) σ . Note that G is a convex set [38]. For any U ∈ Ω ( ) σ ⊂ G and any U ∈ Ω ( ) σ ⊂ G , one has U λ : = λ U + ( − λ ) U ∈ G , ∀ λ ∈ [ , ] . (39) Let E ( U ) = − DS ( U ) . According to Proposition 2.1, E ( U ) is a convex function of U on G . UsingJensen’s inequality gives E ( U λ ) ≤ λ E ( U ) + ( − λ ) E ( U ) for all λ ∈ [ , ] . Let D λ > denotethe first component of U λ . It follows that − D λ S ( U λ ) ≤ − λ D S ( U ) − ( − λ ) D S ( U ) ≤ − λ D σ − ( − λ ) D σ ≤ − D λ min { σ , σ } . This implies S ( U λ ) ≥ min { σ , σ } , which along with (39) yield U λ ∈ Ω min { σ , σ } , ∀ λ ∈ [ , ] . .3 Generalized Lax-Friedrichs splitting properties In the bound-preserving analysis of numerical schemes with the Lax-Friedrichs (LF) flux, thefollowing property (40) is usually expected: U ± F i ( U ) α ∈ Ω σ , ∀ U ∈ Ω σ , ∀ α ≥ α i , (40)where α i denotes a suitable upper bound of the wave speeds in the x i -direction, and, in the RHDcase, it can be taken as the speed of light c = LF splittingproperty . This property is valid for the admissible state set G or G and played an important rolein constructing bound-preserving schemes for the RHD; see [38, 28, 41]. However, unfortunately,the property (40) does not hold in general for the invariant region Ω σ as the entropy principle S ( U ) ≥ σ is included.Since (40) does not hold, we would like to look for some alternative properties that are valid butweaker than (40). By considering the convex combination of some LF splitting terms, we achievethe generalized LF (gLF) splitting properties , whose derivations are highly nontrivial. Built on atechnical inequality (41) constructed in Section 3.3.1, the gLF splitting properties are presented inSection 3.3.2. We first construct an important inequality (41), which will be the key to establishing the gLFsplitting properties.
Theorem 3.2. If U ∈ Ω σ , then for any vvv ∗ = ( v , ∗ , . . . , v d , ∗ ) (cid:62) ∈ B ( ) , any ρ ∗ ∈ R + , and any θ ∈ [ − , ] , it holds ϕ σ (cid:16) U + θ F i ( U ) ; vvv ∗ , ρ ∗ (cid:17) + θ e σ v i , ∗ ρ Γ ∗ ≥ , (41) where i ∈ { , . . . , d } , and the function ϕ σ is defined in (32) .Proof. Due to the relativistic effects, the flux F i ( U ) also cannot be explicitly formulated in termsof U . Therefore, we have to work on the corresponding primitive quantities { ρ , vvv , p } of U , whichsatisfy ρ > (cid:107) vvv (cid:107) < p ≥ e σ ρ Γ because U ∈ Ω σ . We observe that ϕ σ (cid:16) U + θ F i ( U ) ; vvv ∗ , ρ ∗ (cid:17) + θ e σ v i , ∗ ρ Γ ∗ = (cid:98) Π + (cid:98) Π p + e σ (cid:98) Π , with (cid:98) Π : = ρ W ( + θ v i ) (cid:18) − vvv · vvv ∗ − (cid:113) − (cid:107) vvv (cid:107) (cid:113) − (cid:107) vvv ∗ (cid:107) (cid:19) , (cid:98) Π : = ΓΓ − ( + θ v i ) (cid:18) − vvv · vvv ∗ − (cid:107) vvv (cid:107) (cid:19) − ( + θ v i , ∗ ) , (cid:98) Π : = ρ Γ ∗ ( + θ v i , ∗ ) − ΓΓ − ( + θ v i ) ρ W ρ Γ − ∗ (cid:113) − (cid:107) vvv ∗ (cid:107) . vvv · vvv ∗ + (cid:113) − (cid:107) vvv (cid:107) (cid:113) − (cid:107) vvv ∗ (cid:107) ≤ (cid:113) (cid:107) vvv (cid:107) + ( − (cid:107) vvv (cid:107) ) (cid:113) (cid:107) vvv ∗ (cid:107) + ( − (cid:107) vvv ∗ (cid:107) ) = , which implies (cid:98) Π ≥
0. It follows that ϕ σ (cid:16) U + θ F i ( U ) ; vvv ∗ , ρ ∗ (cid:17) + θ e σ v i , ∗ ρ Γ ∗ ≥ (cid:98) Π p + e σ (cid:98) Π . (42)Recalling that Γ ∈ ( , ] , (cid:107) vvv (cid:107) < (cid:107) vvv ∗ (cid:107) <
1, we obtain (cid:98) Π ≥ (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) ( − vvv · vvv ∗ ) − ( + θ v i , ∗ ) = : (cid:101) Π > , (43)where the positivity of (cid:101) Π is deduced by using the Cauchy–Schwarz inequality as follows: (cid:101) Π = (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) − − (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) (cid:34)(cid:18) v i + θ ( − (cid:107) vvv (cid:107) ) ( + θ v i ) (cid:19) v i , ∗ + ∑ j (cid:54) = i v j v j , ∗ (cid:35) ≥ (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) − − (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) (cid:34)(cid:18) v i + θ ( − (cid:107) vvv (cid:107) ) ( + θ v i ) (cid:19) + ∑ j (cid:54) = i v j (cid:35) (cid:107) vvv ∗ (cid:107) = (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) − − (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) (cid:107) vvv ∗ (cid:107) (cid:20) − − (cid:107) vvv (cid:107) + θ v i + θ ( − (cid:107) vvv (cid:107) ) ( + θ v i ) (cid:21) ≥ (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) − − (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) (cid:107) vvv ∗ (cid:107) (cid:34)(cid:18) − − (cid:107) vvv (cid:107) ( + θ v i ) (cid:19) (cid:35) = (cid:20) (cid:18) + θ v i − (cid:107) vvv (cid:107) (cid:19) − (cid:21) ( − (cid:107) vvv ∗ (cid:107) ) ≥ (cid:20) (cid:18) − (cid:107) vvv (cid:107) − (cid:107) vvv (cid:107) (cid:19) − (cid:21) ( − (cid:107) vvv ∗ (cid:107) ) > . (44)Combining (cid:98) Π > p ≥ e σ ρ Γ , we then derive from (42) that ϕ σ (cid:16) U + θ F i ( U ) ; vvv ∗ , ρ ∗ (cid:17) + θ e σ v i , ∗ ρ Γ ∗ ≥ (cid:98) Π e σ ρ Γ + e σ (cid:98) Π = e σ ρ Γ (cid:98) Π , (45)with (cid:98) Π : = (cid:98) Π + ρ − Γ (cid:98) Π = ˆ φ σ (cid:0) ρρ ∗ ; vvv , vvv ∗ (cid:1) , andˆ φ σ ( x ; vvv , vvv ∗ ) : = (cid:98) Π + x − Γ ( + θ v i , ∗ ) − x − Γ + (cid:32) ΓΓ − ( + θ v i ) (cid:112) − (cid:107) vvv ∗ (cid:107) (cid:112) − (cid:107) vvv (cid:107) (cid:33) . The subsequent task is to show that (cid:98) Π is always nonnegative. Define η s : = + θ v i + θ v i , ∗ (cid:18) − vvv · vvv ∗ − (cid:107) vvv (cid:107) (cid:19) − , x s : = ( + θ v i , ∗ ) (cid:112) − (cid:107) vvv (cid:107) ( + θ v i ) (cid:112) − (cid:107) vvv ∗ (cid:107) .
17y studying the derivative of ˆ φ σ with respect to x , we observe that the function ˆ φ σ is strictlydecreasing on the interval ( , x s ] and strictly increasing on [ x s , + ∞ ) . We therefore have (cid:98) Π ≥ min x ∈ R + ˆ φ σ ( x ; vvv , vvv ∗ ) = ˆ φ σ ( x s ; vvv , vvv ∗ )= (cid:98) Π − Γ − x − Γ s ( + θ v i , ∗ ) = + θ v i , ∗ Γ − (cid:16) η s Γ + − x − Γ s (cid:17) . (46)Using the formulation of (cid:101) Π defined in (43), we reformulate η s and observe that η s = (cid:32) (cid:101) Π + θ v i , ∗ + (cid:33) − = (cid:101) Π ( + θ v i , ∗ ) − > − , where the last step follows from the positivity of (cid:101) Π , which has been proven in (44). Thanks toLemma 3.3, we obtain ( η s Γ + ) / Γ ≥ ( η s + ) / , or equivalently, η s Γ + ≥ ( η s + ) Γ , which, along with (46), imply (cid:98) Π ≥ + θ v i , ∗ Γ − (cid:20)(cid:0) η s + (cid:1) Γ − (cid:0) x − s (cid:1) Γ (cid:21) . (47)Next, we would like to show that 2 η s + ≥ x − s . Define a : = θ v i , a : = (cid:113) (cid:107) vvv (cid:107) − θ v i , b : = θ v i , ∗ , and b : = (cid:113) (cid:107) vvv ∗ (cid:107) − θ v i , ∗ . By an elementary inequality ( + a ) ( − b − b ) + ( + b ) ( − a − a ) ≤ ( − a b − a b )( + a )( + b ) , we get x − s = ( + a ) ( − b − b )( + b ) ( − a − a ) ≤ ( − a b − a b )( + a )( + b )( − a − a ) − = (cid:101) η s + (cid:101) η s : = + θ v i + θ v i , ∗ (cid:16) − a b − a b −(cid:107) vvv (cid:107) (cid:17) −
1. With the aid of the Cauchy–Schwarz inequality, we obtain vvv · vvv ∗ = θ v i v i , ∗ + (cid:16) v i (cid:112) − θ (cid:17) (cid:16) v i , ∗ (cid:112) − θ (cid:17) + ∑ j (cid:54) = i v j v j , ∗ ≤ θ v i v i , ∗ + (cid:32)(cid:16) v i (cid:112) − θ (cid:17) + ∑ j (cid:54) = i v j (cid:33) (cid:32)(cid:16) v i , ∗ (cid:112) − θ (cid:17) + ∑ j (cid:54) = i v j , ∗ (cid:33) = a b + a b , which implies η s ≥ (cid:101) η s . It then follows from (48) that x − s ≤ η s +
1, which together with (47) imply (cid:98) Π ≥
0. Then by using (45), we conclude ϕ σ (cid:16) U + θ F i ( U ) ; vvv ∗ , ρ ∗ (cid:17) + θ e σ v i , ∗ ρ Γ ∗ ≥ e σ ρ Γ (cid:98) Π ≥ (cid:4) Subtracting the right-hand side terms from the left-hand side terms leads to − ( a − b − a b + a b ) ≤ emark 3.2. It is worth emphasizing the importance of the last term ( θ e σ v i , ∗ ρ Γ ∗ ) at the left-handside of (41) . This term is very technical, necessary and crucial in deriving the gLF splittingproperties. The value of this term is not always positive or negative. However, without this term,the inequality (41) does not hold. More importantly, this term can be canceled out dexterously inIRP analysis; see the proofs of gLF splitting properties in the theorems in Section 3.3.2. We first present the one-dimensional version.
Theorem 3.3 (
1D gLF splitting). If ˆ U = ( ˆ D , ˆ mmm (cid:62) , ˆ E ) (cid:62) ∈ Ω σ and ˇ U = ( ˇ D , ˇ mmm (cid:62) , ˇ E ) (cid:62) ∈ Ω σ , then forany α ≥ c = and any given i ∈ { , . . . , d } , it holds G i , α ( ˆ U , ˇ U ) : = (cid:18) ˆ U − F i ( ˆ U ) α + ˇ U + F i ( ˇ U ) α (cid:19) ∈ Ω σ . (49) Proof.
The first component of G i , α equals (cid:0) ˆ D (cid:0) − ˆ v i α (cid:1) + ˇ D (cid:0) + ˇ v i α (cid:1)(cid:1) >
0. With the help of Theo-rem 3.2, we obtain, for any vvv ∗ ∈ B ( ) and any ρ ∗ ∈ R + , that2 ϕ σ (cid:0) G i , α ; vvv ∗ , ρ ∗ (cid:1) = ϕ σ (cid:18) U − F i ( U ) α ; vvv ∗ , ρ ∗ (cid:19) − e σ v i , ∗ ρ Γ ∗ α + ϕ σ (cid:18) U + F i ( U ) α ; vvv ∗ , ρ ∗ (cid:19) + e σ v i , ∗ ρ Γ ∗ α ≥ , where we have used the linearity of ϕ σ ( U ; vvv ∗ , ρ ∗ ) with respect to U . Then, using Theorem 3.1concludes that G i , α ∈ Ω ( ) σ = Ω σ . (cid:4) Remark 3.3.
Another approach to show Theorem 3.3 is based on the assumption that the exactRiemann solver preserves the invariant domain, which is reasonable but not yet proven for RHD.Our present analysis approach is very direct, does not rely on any assumption, and would motivatethe further study of IRP schemes for some other complicated physical systems such as the MHDand RMHD equations.
Next, we present the multidimensional gLF splitting property on a general polygonal or poly-hedron cell. For any vector ξξξ = ( ξ , · · · , ξ d ) (cid:62) ∈ R d , we define ξξξ · F ( U ) : = ∑ dk = ξ k F k ( U ) . Theorem 3.4.
For ≤ j ≤ N, let s j > and the unit vector ξξξ ( j ) = ( ξ ( j ) , . . . , ξ ( j ) d ) (cid:62) satisfy N ∑ j = s j ξξξ ( j ) = . (50) Given admissible states U ( i j ) ∈ Ω σ , ≤ i ≤ Q, ≤ j ≤ N, then for any α ≥ c = , it holds U : = N ∑ j = s j N ∑ j = Q ∑ i = s j ω i (cid:18) U ( i j ) − α ξξξ ( j ) · F ( U ( i j ) ) (cid:19) ∈ Ω σ , (51) where the sum of all positive numbers { ω i } Qi = equals one. d -polytope with N edges ( d =
2) or faces ( d = j on the variables in Theorem 3.4 represents the j th edge or face of the polytope, while i stands for the i th (quadrature) point on each edge or face, with ω i denoting the associated quadra-ture weight at that point. Besides, s j and ξξξ ( j ) respectively correspond to the ( d − ) -dimensionalHausdorff measure and the unit outward normal vector of the j th edge or face. One can verify thatthe condition (50) holds naturally. In addition, U ( i j ) stands for the approximate values of U at the i th quadrature point on the j th edge or face. Proof.
Let Q ( j ) ξξξ ∈ R d × d be a rotational matrix associated with the unit vector ξξξ ( j ) and satisfying e (cid:62) Q ( j ) ξξξ = ( ξξξ ( j ) ) (cid:62) , (52)where e = ( , (cid:62) d − ) (cid:62) , and d − denotes the zero vector in R d − . It can be verified that the system(1) satisfies the following rotational invariance property ξξξ ( j ) · F ( U ( i j ) ) = Q − j F ( Q j U ( i j ) ) , (53)where Q j : = diag { , Q ( j ) ξξξ , } . Notice that the matrix Q ( j ) ξξξ is orthogonal. For any fixed j and any vvv ∗ ∈ B ( ) , define (cid:98) vvv ∗ : = Q ( j ) ξξξ vvv ∗ ∈ B ( ) . Utilizing (52) gives (cid:98) v , ∗ = e (cid:62) vvv ∗ = e (cid:62) Q ( j ) ξξξ vvv ∗ = ξξξ ( j ) · vvv ∗ . (54)For U ( i j ) ∈ Ω σ , with the aid of the first equivalent form Ω ( ) σ in (30), one can verify that (cid:98) U ( i j ) : = Q j U ( i j ) ∈ Ω ( ) σ = Ω σ . For any vvv ∗ ∈ B ( ) and any ρ ∗ >
0, we have ϕ σ (cid:16) U ( i j ) − α − ξξξ ( j ) · F ( U ( i j ) ) ; vvv ∗ , ρ ∗ (cid:17) − α − e σ (cid:16) ξξξ ( j ) · vvv ∗ (cid:17) ρ Γ ∗ = ϕ σ (cid:16) Q − j (cid:98) U ( i j ) − α − Q − j F ( (cid:98) U ( i j ) ) ; ( Q ( j ) ξξξ ) − (cid:98) vvv ∗ , ρ ∗ (cid:17) − α − e σ (cid:98) v , ∗ ρ Γ ∗ = ϕ σ (cid:16)(cid:98) U ( i j ) − α − F ( (cid:98) U ( i j ) ) ; (cid:98) vvv ∗ , ρ ∗ (cid:17) − α − e σ (cid:98) v , ∗ ρ Γ ∗ ≥ , (55)where we have used (53)–(54) in the first equality and the orthogonality of Q ( j ) ξξξ in the secondequality, and the inequality follows from Theorem 3.2 for (cid:98) U ( i j ) ∈ Ω σ , (cid:98) vvv ∗ ∈ B ( ) , ρ ∗ ∈ R + and0 < α − ≤
1. It then follows that (cid:32) N ∑ j = s j (cid:33) ϕ σ (cid:0) U ; vvv ∗ , ρ ∗ (cid:1) = N ∑ j = Q ∑ i = s j ω i ϕ σ (cid:16) U ( i j ) − α − ξξξ ( j ) · F ( U ( i j ) ) ; vvv ∗ , ρ ∗ (cid:17) ≥ α − e σ N ∑ j = Q ∑ i = s j ω i (cid:16) ξξξ ( j ) · vvv ∗ (cid:17) ρ Γ ∗ = α − e σ ρ Γ ∗ (cid:32) N ∑ j = s j ξξξ ( j ) (cid:33) · vvv ∗ = , ϕ σ ( U ; vvv ∗ , ρ ∗ ) with respect to U , the inequality (55), and thecondition (50). Therefore, ϕ σ ( U ; vvv ∗ , ρ ∗ ) ≥ vvv ∗ ∈ B ( ) and any ρ ∗ >
0. Note that the firstcomponent of U equals to1 N ∑ j = s j N ∑ j = Q ∑ i = s j ω i D ( i j ) (cid:18) − α ξξξ ( j ) · vvv ( i j ) (cid:19) ≥ N ∑ j = s j N ∑ j = Q ∑ i = s j ω i D ( i j ) (cid:18) − (cid:107) vvv ( i j ) (cid:107) α (cid:19) > . Hence, U ∈ Ω ( ) σ . Thanks to Theorem 3.1, we have U ∈ Ω σ . The proof is completed. (cid:4) As two special cases of Theorem 3.4, the following corollaries show the gLF splitting proper-ties on 2D and 3D Cartesian mesh cells.
Corollary 3.2. If ¯ U i , ˜ U i , ˆ U i , ˇ U i ∈ Ω σ for i = , · · · , Q, then for any α ≥ c = , it holds U : = (cid:16) ∆ x + ∆ y (cid:17) Q ∑ i = ω i (cid:32) G , α ( ¯ U i , ˜ U i ) ∆ x + G , α ( ˆ U i , ˇ U i ) ∆ y (cid:33) ∈ Ω σ , (56) where ∆ x > , ∆ y > , and the sum of the positive numbers { ω i } Qi = equals one. Corollary 3.3. If ¯ U i , ˜ U i , ˆ U i , ˇ U i , ´ U i , ` U i ∈ Ω σ for i = , · · · , Q, then for any α ≥ , it holds U : = (cid:16) ∆ x + ∆ y + ∆ z (cid:17) Q ∑ i = ω i (cid:32) G , α ( ¯ U i , ˜ U i ) ∆ x + G , α ( ˆ U i , ˇ U i ) ∆ y + G , α ( ´ U i , ` U i ) ∆ z (cid:33) ∈ Ω σ , (57) where ∆ x > , ∆ y > , ∆ z > , and the sum of the positive numbers { ω i } Qi = equals one. This section applies the above theories to study the IRP schemes for the RHD system (1) in onespatial dimension. To avoid confusing subscripts, we will use the symbol x to represent the variable x in (1). Let I j = [ x j − , x j + ] and ∪ j I j be a partition of the spatial domain Σ . Denote ∆ x j = x j + − x j − . Assume that the time interval is divided into the mesh { t n } N t n = with t = ∆ t determined by certain CFL condition. Let U nj denote the numerical cell-averagedapproximation of the exact solution U ( x , t ) over I j at t = t n .Let U ( x ) : = U ( x , ) be an initial solution and S : = min x S ( U ( x )) . We are interested innumerical schemes preserving U nj ∈ Ω S for the system (1).Note that the initial cell average U j ∈ Ω S for all j , according to the following lemma.21 emma 4.1. Assume that U ( x ) is an (admissible) initial data of the RHD system (1) on thedomain Σ , i.e., it satisfies U ( x ) ∈ G for all x ∈ Σ . For any I ⊆ Σ , we have U I : = | I | (cid:90) I U ( x ) d x ∈ Ω S , where | I | = (cid:82) I d x > , and S = min x S ( U ( x )) .Proof. Let U I = ( D I , mmm (cid:62) I , E I ) (cid:62) . It is evident that D I >
0. Define q ( U ) : = E − (cid:112) D + (cid:107) mmm (cid:107) , it is easy to show that q ( U ) is a concave function on R d + , and the second constraint in G isequivalent q ( U ) >
0. According to Jensen’s inequality, q ( U I ) ≥ | I | (cid:82) I q ( U ( xxx )) d xxx >
0. Therefore, U I ∈ G = G .Let E ( U ) = − DS ( U ) . According to Proposition 2.1, E ( U ) is a convex function of U on G . Byfollowing [46, Lemma 2.2] and using Jensen’s inequality E (cid:18) | I | (cid:90) I U ( x ) d x (cid:19) ≤ | I | (cid:90) I E ( U ( x )) d x , one obtains S (cid:0) U I (cid:1) ≥ S , which together with U I ∈ G , imply U I ∈ Ω S . (cid:4) Consider the first-order LF type scheme U n + j = U nj − ∆ t ∆ x j (cid:16)(cid:98) F ( U nj , U nj + ) − (cid:98) F ( U nj − , U nj ) (cid:17) , (58)where the numerical flux (cid:98) F ( · , · ) is defined by (cid:98) F ( U − , U + ) = (cid:16) F ( U − ) + F ( U + ) − α ( U + − U − ) (cid:17) , (59)and α denotes the numerical viscosity parameter, which can be taken as α = c =
1, a simple upperbound of all wave speeds in the theory of special relativity.
Lemma 4.2. If U nj − , U nj and U nj + all belong to Ω σ for certain σ , then the state U n + j , computedby the scheme (58) under the CFL condition α ∆ t ≤ ∆ x j , (60) belongs to Ω σ .Proof. We rewrite the scheme (58) in the following form U n + j = ( − λ ) U nj + λ G , α ( U nj + , U nj − ) , G , α is defined in (49), and λ : = α ∆ t / ∆ x j ∈ ( , ] . Thanks to Theorem 3.3, we have G , α ( U nj + , U nj − ) ∈ Ω σ , which leads to U n + j ∈ Ω σ by the convexity of Ω σ (Lemma 3.5). (cid:4) Lemma 4.2 implies a discrete (local) minimum entropy principle of the scheme (58).
Theorem 4.1. If U nj − , U nj and U nj + all belong to G , then the state U n + j , computed by the scheme (58) under the CFL condition (60) , belongs to G and satisfiesS ( U n + j ) ≥ min (cid:8) S ( U nj − ) , S ( U nj ) , S ( U nj + ) (cid:9) . Proof.
By Lemma 4.2 with σ = min j − ≤ k ≤ j + S ( U nk ) , one directly draws the conclusion. (cid:4) The IRP property of the scheme (58) is shown in the following theorem.
Theorem 4.2.
Under the CFL condition (60) , the scheme (58) always preserves U nj ∈ Ω S for allj and n ≥ .Proof. We prove it by the mathematical induction for the time level number n . By Lemma 4.1, wehave U j ∈ Ω S for all j , i.e., the conclusion holds for n =
0. If assuming that U nj ∈ Ω S for all j ,then, by Lemma 4.2 with σ = S , we have U n + j ∈ Ω S for all j . (cid:4) We now study the IRP high-order DG and finite volume schemes for 1D RHD equations (1).For the moment, we use the forward Euler method for time discretization, while high-order timediscretization will be discussed later. We consider the high-order finite volume schemes as wellas the scheme satisfied by the cell averages of a standard DG method, which have the followingunified form U n + j = U nj − ∆ t ∆ x j (cid:16)(cid:98) F ( U − j + , U + j + ) − (cid:98) F ( U − j − , U + j − ) (cid:17) , (61)where (cid:98) F ( · , · ) is taken as the LF flux in (59), and U ± j + = lim ε → ± U h (cid:0) x j + + ε (cid:1) . (62)Here U h ( x ) is a piecewise polynomial vector function of degree k , i.e., U h ∈ V kh : = (cid:110) u = ( u , · · · , u d + ) (cid:62) : u (cid:96) (cid:12)(cid:12) I j ∈ P k , ∀ (cid:96), j (cid:111) , with P k denoting the space of polynomials of degree up to k . Specifically, the function U h ( x ) with (cid:82) I j U h d x = U nj is an approximation to U ( x , t n ) within the cell I j ; it is either reconstructed in thefinite volume methods from { U nj } or directly evolved in the DG methods. The evolution equationsfor the high-order “moments” of U h ( x ) in the DG methods are omitted because here we are onlyconcerned with the IRP property of the schemes.23 .2.1 Theoretical analysis If the polynomial degree k =
0, i.e., U h ( x ) = U nj , ∀ x ∈ I j , then the scheme (61) reduces to thefirst-order scheme (58), which is IRP under the CFL condition (60). When the polynomial degree k ≥
1, the solution U n + j of the high-order scheme (61) does not always belong to Ω S even if U nj ∈ Ω S for all j . In the following theorem, we present a satisfiable condition for achieving theprovably IRP property of the scheme (61) when k ≥ { (cid:98) x ( µ ) j } L µ = be the L -point Gauss–Lobatto quadrature nodes in the interval I j , and the asso-ciated weights are denoted by { (cid:98) ω µ } L µ = with ∑ L µ = (cid:98) ω µ =
1. We require 2 L − ≥ k such that thealgebraic precision of the quadrature is at least k , for example, one can particularly take L = (cid:100) k + (cid:101) . Theorem 4.3.
If the piecewise polynomial vector function U h satisfies U h ( (cid:98) x ( µ ) j ) ∈ Ω S , ∀ µ ∈ { , , · · · , L } , ∀ j , (63) then, under the CFL condition α ∆ t ∆ x j ≤ (cid:98) ω = L ( L − ) , (64) the solution U n + j , computed by the high-order scheme (61) , belongs to Ω S for all j.Proof. The exactness of the L -point Gauss–Lobatto quadrature rule for the polynomials of degree k implies U nj = ∆ x j (cid:90) I j U h ( x ) d x = L ∑ µ = (cid:98) ω µ U h ( (cid:98) x ( µ ) j ) . Noting (cid:98) ω = (cid:98) ω L and (cid:98) x , Lj = x j ∓ , we can then rewrite the scheme (61) into the convex combinationform U n + j = L − ∑ µ = (cid:98) ω µ U h ( (cid:98) x ( µ ) j ) + ( (cid:98) ω − λ ) (cid:18) U + j − + U − j + (cid:19) + λ Ξ − + λ Ξ + , (65)where λ = α ∆ t n / ∆ x ∈ ( , (cid:98) ω ] , and Ξ ± = U ± j + − F ( U ± j + ) α + U ± j − + F ( U ± j − ) α . According to the gLF splitting property in Theorem 3.3 and U ± j + ∈ Ω S by (63), we obtain Ξ ± ∈ Ω S . Using the convexity of Ω S (Lemma 3.5), we conclude U n + j ∈ Ω S from (65). (cid:4) .2.2 Invariant-region-preserving limiter In general, the high-order scheme (61) does not meet the condition (63) automatically. Now, wedesign a simple limiter to effectively enforce the condition (63), without losing high-order accu-racy and conservation. Our limiter is motivated by the existing bound-preserving and physical-constraint-preserving limiters (cf. [45, 46, 38, 28, 41, 32]).Before presenting the limiter, we define q ( U ) : = E − (cid:113) D + (cid:107) mmm (cid:107) , (66)then the second constraint in (30) becomes q ( U ) >
0. As observed in [38], the function q ( U ) is strictly concave. This function will play an important role, similar to the pressure function inthe non-relativistic case, in our limiter for RHD. (Unlike the non-relativistic case, the pressurefunction p ( U ) is not concave in the RHD case [38].)Denote X j : = { (cid:98) x ( µ ) j } L µ = and V kh : = (cid:26) u ∈ V kh : 1 ∆ x j (cid:90) I j u ( x ) d x ∈ Ω S , ∀ j (cid:27) , (cid:101) V kh : = (cid:110) u ∈ V kh : u (cid:12)(cid:12) I j ( x ) ∈ Ω S , ∀ x ∈ X j , ∀ j (cid:111) . For any U h ∈ V kh with U h (cid:12)(cid:12) I j = : U j ( x ) = (cid:0) D j ( x ) , m j ( x ) , E j ( x ) (cid:1) (cid:62) , we define the IRP limiting operator Π h : V kh → (cid:101) V kh by Π h U h (cid:12)(cid:12) I j = (cid:101) U j ( x ) , ∀ j , (67)with the limited polynomial vector function (cid:101) U j ( x ) constructed via the following three steps. Step (i):
First, modify the density to enforce its positivity via (cid:98) D j ( x ) = θ (cid:0) D j ( x ) − D nj (cid:1) + D nj , θ : = min (cid:40) , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) D nj − ε D nj − min x ∈ X j D j ( x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:41) , (68)where ε is a small positive number as the desired lower bound for density, is introduced toavoid the effect of the round-off error, and can be taken as ε = min { − , D nj } . Step (ii):
Then, modify (cid:98) U j ( x ) = (cid:0) (cid:98) D j ( x ) , m j ( x ) , E j ( x ) (cid:1) (cid:62) to enforce the positivity of q ( U ) via ( U j ( x ) = θ (cid:16)(cid:98) U j ( x ) − U nj (cid:17) + U nj , θ : = min (cid:40) , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q ( U j ) − ε q ( U nj ) − min x ∈ X j q ( (cid:98) U j ( x )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:41) , (69)where ε is a small positive number as the desired lower bound for q ( U ) , is introduced toavoid the effect of the round-off error, and can be taken as ε = min { − , q ( U nj ) } .25 tep (iii): Finally, modify ( U j ( x ) to enforce the entropy principle S ( U ) ≥ S via (cid:101) U j ( x ) = θ (cid:16) ( U j ( x ) − U nj (cid:17) + U nj , θ : = min x ∈ X j ˜ θ ( x ) , (70)where, for x ∈ { x ∈ X j : S ( ( U j ( x )) ≥ S } , ˜ θ ( x ) =
1, and, for x ∈ { x ∈ X j : S ( ( U j ( x )) < S } ,˜ θ ( x ) is the (unique) solution to the nonlinear equation S (cid:16) ( − ˜ θ ) U nj + ˜ θ ( U j ( x ) (cid:17) = S , ˜ θ ∈ [ , ) . The limiter Π h is a combination of the bound-preserving limiter (68)–(69) (cf. [28]) and the entropy limiter (70). According to the above definition of the limiter Π h and the Jensen’s inequalityfor the concave function q ( U ) , we immediately obtain the following proposition: Proposition 4.1.
For any U h ∈ V kh , one has Π h U h ∈ (cid:101) V kh . Proposition 4.1 indicates that the limited solution (67) satisfies the condition (63). Note thatsuch type of local scaling limiters keep the conservation (cid:82) I j Π h ( u ) d x = (cid:82) I j u d x , ∀ u ∈ V kh , and donot destroy the high-order accuracy; see [44, 45, 43] for details. Remark 4.1.
The invariant region Ω S can also be reformulated as Ω ( ) S = (cid:110) U = ( D , mmm (cid:62) , E ) (cid:62) ∈ R d + : D > , q ( U ) > , (cid:101) q ( U ) ≥ (cid:111) , where (cid:101) q ( U ) : = D ( S ( U ) − S ) , and we have used D > to reformulate the third constraint S ( U ) ≥ S in (30) . Thanks to Proposition 2.1, the function (cid:101) q ( U ) is strictly concave for U ∈ G . Motivatedby this property and [19], we can also use another (simpler but possibly more restrictive) approachto enforce the entropy principle S ( U ) ≥ S by modifying ( U j ( x ) to (cid:101) U j ( x ) = ˜ θ (cid:16) ( U j ( x ) − U nj (cid:17) + U nj , ˜ θ : = min (cid:40) , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) q ( U nj ) (cid:101) q ( U nj ) − min x ∈ X j (cid:101) q ( ( U j ( x )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:41) . Remark 4.2.
Another way to define an IRP limiting operator is (cid:101) Π h U h (cid:12)(cid:12) I j = θ j ( U j ( x ) − U nj ) + U nj , ∀ j , with θ j = min (cid:26) , D nj − ε D nj − min x ∈ X j D j ( x ) , q ( U nj ) − ε q ( U nj ) − min x ∈ X j q ( U j ( x )) , (cid:101) q ( U nj ) (cid:101) q ( U nj ) − min x ∈ X j (cid:101) q ( U j ( x )) (cid:27) . This simpler def-inition is motivated by the IRP limiter in [19] for the non-relativistic Euler system. It seems,however, unpractical for numerical implementation in our RHD case, because U j ( x ) does notnecessarily belong to G for all x ∈ X j so that (cid:101) q ( U j ( x )) is generally not well-defined in θ j . Π h defined in (67), we modify the high-order scheme (61) into U n + j = U nj − ∆ t ∆ x j (cid:16)(cid:98) F ( (cid:101) U − j + , (cid:101) U + j + ) − (cid:98) F ( (cid:101) U − j − , (cid:101) U + j − ) (cid:17) = : U nj + ∆ t L j ( Π h U h ) , (71)where (cid:101) U ± j + = lim ε → ± Π h U h (cid:0) x j + + ε (cid:1) . (72)Based on Theorem 4.3 and Proposition 4.1, we know that the resulting scheme (71) is IRP underthe CFL condition (64).The scheme (71) is only first-order accurate in time. To achieve high-order accurate IRPscheme in both time and space, one can replace the forward Euler time discretization in (71)with any high-order accurate strong-stability-preserving (SSP) methods [5]. For example, thethird-order accurate SSP Runge-Kutta (SSP-RK) method: U ∗ j = U nj + ∆ t L j ( Π h U h ) , U ∗∗ j = U nj + (cid:16) U ∗ j + ∆ t L j ( Π h U ∗ h ) (cid:17) U n + j = U nj + (cid:16) U ∗∗ j + ∆ t L j ( Π h U ∗∗ h ) (cid:17) , (73)and the third-order accurate SSP multi-step (SSP-MS) method: U n + j = ( U nj + ∆ t L j ( Π h U nh )) + (cid:18) U n − j + ∆ t L j ( Π h U n − h ) (cid:19) . (74)Since SSP methods are formally convex combinations of the forward Euler method, the resultinghigh-order schemes (73) and (74) are still IRP, according to the convexity of Ω S (Lemma 3.5). In this section, we study IRP schemes for the two-dimensional (2D) RHD equations, keeping inmind that the proposed methods and analyses are extensible to the 3D case. Assume that thephysical domain Σ in the 2D space is discretized by a mesh T h . In general, the mesh may beunstructured and consists of polygonal cells. We also partition the time interval into a mesh { t n } N t n = with t = ∆ t determined by certain CFL condition.We will use the capital letter K to denote an arbitrary cell in T h . Let E jK , j = , · · · , N K , denotethe edges of K , and K j be the adjacent cell which shares the edge E jK with K . We denote by ξξξ ( j ) K = (cid:0) ξ ( j ) , K , ξ ( j ) , K (cid:1) the unit normal vector of E jK pointing from K to K j . The notations | K | and | E jK | are used to denote the area of K and the length of E jK , respectively.Let U nK denote the numerical cell-averaged approximation of the exact solution U ( xxx , t ) over K at t = t n . Define U ( xxx ) : = U ( xxx , ) as the initial data and S : = min xxx S ( U ( xxx )) . According toLemma 4.1, the initial cell average U K always belongs to Ω S for all K ∈ T h . We would like toseek numerical schemes maintaining U nK ∈ Ω S for all K ∈ T h and n ≥ .1 First-order scheme Consider the following first-order scheme on the mesh T h for the RHD equations (1): U n + K = U nK − ∆ t | K | N K ∑ j = (cid:12)(cid:12)(cid:12) E jK (cid:12)(cid:12)(cid:12) (cid:98) F (cid:0) U nK , U nK j ; ξξξ ( j ) K (cid:1) , (75)with the numerical flux (cid:98) F taken as the LF flux (cid:98) F (cid:0) U − , U + ; ξξξ (cid:1) = (cid:16) ξξξ · F ( U − ) + ξξξ · F ( U + ) − α ( U + − U − h ) (cid:17) , (76)where the numerical viscosity parameter α is chosen as the speed of light in vacuum c =
1, whichis a simple upper bound of all wave speeds in the theory of special relativity.
Lemma 5.1. If U nK ∈ Ω σ , U nK j ∈ Ω σ , ≤ j ≤ N K , for certain σ , then the state U n + K , computed bythe scheme (75) under the CFL condition α ∆ t | K | N K ∑ j = (cid:12)(cid:12)(cid:12) E jK (cid:12)(cid:12)(cid:12) ≤ , (77) belongs to Ω σ .Proof. Substituting the numerical flux (76) into the scheme (75) and then using the identity N K ∑ j = (cid:12)(cid:12)(cid:12) E jK (cid:12)(cid:12)(cid:12) ξξξ ( j ) K = , we rewrite the scheme (75) in the following form U n + K = ( − λ ) U nK + λ Ξ , where λ : = α ∆ t | K | ∑ N K j = (cid:12)(cid:12)(cid:12) E jK (cid:12)(cid:12)(cid:12) ∈ ( , ] under the condition (77), and Ξ : = ∑ N K j = (cid:12)(cid:12)(cid:12) E jK (cid:12)(cid:12)(cid:12) N K ∑ j = (cid:12)(cid:12)(cid:12) E jK (cid:12)(cid:12)(cid:12) (cid:18) U nK j − α ξξξ ( j ) K · F (cid:0) U nK j (cid:1)(cid:19) . Thanks to the gLF splitting property in Theorem 3.4, we have Ξ ∈ Ω σ , which leads to U n + K ∈ Ω σ by the convexity of Ω σ (Lemma 3.5). (cid:4) Lemma 5.1 implies a discrete (local) minimum entropy principle of the scheme (75).
Theorem 5.1. If U nK ∈ G , U nK j ∈ G , ≤ j ≤ N K , then the state U n + K , computed by the scheme (75) under the CFL condition (77) , belongs to G and satisfiesS ( U n + K ) ≥ min (cid:26) S ( U nK ) , min ≤ j ≤ N K S ( U nK j ) (cid:27) . roof. By Lemma 5.1 with σ = min (cid:8) S ( U nK ) , min ≤ j ≤ N K S ( U nK j ) (cid:9) , we directly draw the conclusion. (cid:4) The IRP property of the scheme (75) is shown in the following theorem.
Theorem 5.2.
Under the CFL condition (77) , the scheme (75) always preserves U nK ∈ Ω S for allK ∈ T h and n ≥ .Proof. The conclusion follows from Lemma 5.1 with σ = S and the principle of mathematicalinduction for the time level number n . (cid:4) This subsection discusses the provably IRP high-order finite volume or DG schemes for the 2DRHD equations (1). We will focus on the first-order forward Euler method for time discretization,and our analysis also works for high-order explicit time discretization using the SSP methods,which are formed by convex combinations of the forward Euler method [5].To achieve ( k + ) th-order accuracy in space, an piecewise polynomial vector function U h ( xxx ) (i.e., U h | K is a polynomial vector of degree k for all K ∈ T h ) is also built, as approximation tothe exact solution U ( xxx , t n ) . It is, either evolved in the DG methods, or reconstructed in the finitevolume methods from the cell averages { U nK : K ∈ T h } . Moreover, the cell average of U h ( xxx ) over K is equal to U nK . A high-order finite volume scheme as well as the scheme satisfied by the cellaverages of a standard DG method can then be written as U n + K = ¯ U nK − ∆ t | K | N K ∑ j = | E jK | (cid:98) F E jK , (78)with (cid:98) F E jK = Q ∑ ν = ω ν (cid:98) F (cid:16) U int ( K ) h ( xxx ( j ν ) K ) , U ext ( K ) h ( xxx ( j ν ) K ) ; ξξξ ( j ) K (cid:17) (79) ≈ | E jK | (cid:90) E jK (cid:98) F (cid:16) U int ( K ) h ( xxx ) , U ext ( K ) h ( xxx ) ; ξξξ ( j ) K (cid:17) d s . Here the superscripts “ext ( K ) ” and “int ( K ) ” indicate that the corresponding limits of U h ( xxx ) at thecell edges are taken from the exterior and interior of K , respectively; the numerical flux (cid:98) F is takenas the LF flux defined in (76); { xxx ( j ν ) K , ω ν } ≤ ν ≤ Q denote the Q -point Gauss quadrature nodes andweights on E jK .In the following theorem, we derive a satisfiable condition for achieving the provably IRPproperty of the scheme (78) when the polynomial degree k ≥
1. Assume that one can exactly decompose the cell average by certain 2D quadrature: U nK = | K | (cid:90) K U h ( xxx ) d xxx = N K ∑ j = Q ∑ ν = ϖ j ν U int ( K ) h ( xxx ( j ν ) K ) + (cid:101) Q ∑ β = (cid:101) ϖ β U int ( K ) h ( (cid:101) xxx ( β ) K ) , (80)29here { (cid:101) xxx ( β ) K } are the (possible) involved quadrature points excluding { xxx ( j ν ) K } in the cell K ; { ϖ j ν } and { (cid:101) ϖ β } are positive weights satisfying ∑ N K j = ∑ Q ν = ϖ j ν + ∑ (cid:101) Q β = (cid:101) ϖ β =
1. Such a quadrature-baseddecomposition was first proposed by Zhang and Shu in [44, 45] on rectangular cells by tensorproducts of Gauss and Gauss–Lobatto quadratures. It can also be designed on triangular cells andmore general polygons, as demonstrated in, e.g., [47, 23]. Define X K : = (cid:110) xxx ( j ν ) K (cid:111) ≤ j ≤ N K , ≤ ν ≤ Q (cid:91) (cid:110)(cid:101) xxx ( β ) K (cid:111) ≤ β ≤ (cid:101) Q . (81) Theorem 5.3.
If the piecewise polynomial vector function U h satisfies U h ( xxx ) ∈ Ω S , ∀ xxx ∈ X K , ∀ K ∈ T h , (82) then, under the CFL condition α ∆ t | E jK || K | ≤ min ≤ ν ≤ Q ϖ j ν ω ν , ≤ j ≤ N K , ∀ K ∈ T h , (83) the solution U n + K , computed by the high-order scheme (78) , belongs to Ω S for all K ∈ T h .Proof. Substituting the decomposition (80) and the numerical flux (79) with (76) into (78), we canrewrite the scheme (78) in the following convex combination form U n + K = N K ∑ j = Q ∑ ν = (cid:32) ϖ j ν − α ∆ t ω ν | E jK || K | (cid:33) U int ( K ) h ( xxx ( j ν ) K )+ (cid:101) Q ∑ β = (cid:101) ϖ β U int ( K ) h ( (cid:101) xxx ( β ) K ) + α ∆ t | K | (cid:32) N ∑ j = | E jK | (cid:33) (cid:16) Ξ int ( K ) + Ξ ext ( K ) (cid:17) , (84)with Ξ int ( K ) : = N K ∑ j = | E jK | N K ∑ j = Q ∑ ν = | E jK | ω ν (cid:18) U int ( K ) h ( xxx ( j ν ) K ) − α ξξξ ( j ) K · F (cid:16) U int ( K ) h ( xxx ( j ν ) K ) (cid:17)(cid:19) , Ξ ext ( K ) : = N K ∑ j = | E jK | N K ∑ j = Q ∑ ν = | E jK | ω ν (cid:18) U ext ( K ) h ( xxx ( j ν ) K ) − α ξξξ ( j ) K · F (cid:16) U ext ( K ) h ( xxx ( j ν ) K ) (cid:17)(cid:19) . Thanks to the gLF splitting property in Theorem 3.4, under the assumption (82) we obtain Ξ int ( K ) ∈ Ω S and Ξ ext ( K ) ∈ Ω S . Using the convexity of Ω S (Lemma 3.5), we conclude U n + K ∈ Ω S fromthe convex combination form (84) under the condition (83). (cid:4) Theorem 5.3 provides a sufficient condition (82) for the high-order scheme (78) to be IRP.The condition (82), which is not satisfied automatically in general, can again be enforced by30 simple IRP limiting operator Π h similar to the 1D case; see Section 4.2.2 with the 1D pointset X j replaced by the 2D point set X K (81) accordingly. With the IRP limiter applied to theapproximation solution (cid:101) U h = Π h U h , the resulting scheme U n + K = ¯ U nK − ∆ t | K | N K ∑ j = Q ∑ ν = | E jK | ω ν (cid:98) F (cid:16)(cid:101) U int ( K ) h ( xxx ( j ν ) K ) , (cid:101) U ext ( K ) h ( xxx ( j ν ) K ) ; ξξξ ( j ) K (cid:17) , is IRP and also high-order accurate in space. As the 1D case, Theorem 5.3 also remains valid if ahigh-order SSP time discretization [5] is used. Remark 5.1.
Assume that the mesh is rectangular with cells { [ x i − / , x i + / ] × [ y (cid:96) − / , y (cid:96) + / ] } and spatial step-sizes ∆ x i = x i + / − x i − / and ∆ y (cid:96) = y (cid:96) + / − y (cid:96) − / in x- and y-directions re-spectively, where ( x , y ) denotes the 2D spatial coordinate variables. Let S xi = { x ( µ ) i } Q µ = and S y (cid:96) = { y ( µ ) (cid:96) } Q µ = denote the Q-point Gauss quadrature nodes in the intervals [ x i − / , x i + / ] and [ y (cid:96) − / , y (cid:96) + / ] respectively. Let (cid:98) S xi = { (cid:98) x ( ν ) i } L ν = and (cid:98) S y (cid:96) = { (cid:98) y ( ν ) (cid:96) } L ν = denote the L-point (L ≥ k + )Gauss–Lobatto quadrature nodes in the intervals [ x i − / , x i + / ] and [ y (cid:96) − / , y (cid:96) + / ] respectively.For the cell K = [ x i − / , x i + / ] × [ y (cid:96) − / , y (cid:96) + / ] , a suitable point set as X K in (81) is given by(cf. [44]) X K = (cid:0)(cid:98) S xi ⊗ S y (cid:96) (cid:1) ∪ (cid:0) S xi ⊗ (cid:98) S y (cid:96) (cid:1) , (85) and the corresponding 2D quadrature [44] satisfying (80) can be constructed as | K | (cid:90) K u ( x ) d x = Q ∑ µ = ∆ x i (cid:98) ω ω µ ∆ x i + ∆ y (cid:96) (cid:16) u (cid:0) x ( µ ) i , y (cid:96) − (cid:1) + u (cid:0) x ( µ ) i , y (cid:96) + (cid:1)(cid:17) + Q ∑ µ = ∆ y (cid:96) (cid:98) ω ω µ ∆ x i + ∆ y (cid:96) (cid:16) u (cid:0) x i − , y ( µ ) (cid:96) (cid:1) + u (cid:0) x i + , y ( µ ) (cid:96) (cid:1)(cid:17) + L − ∑ ν = Q ∑ µ = (cid:98) ω ν ω µ ∆ x i + ∆ y (cid:96) (cid:16) ∆ x i u (cid:0) x ( µ ) i , (cid:98) y ( ν ) (cid:96) (cid:1) + ∆ y (cid:96) u (cid:0)(cid:98) x ( ν ) i , y ( µ ) (cid:96) (cid:1)(cid:17) , ∀ u ∈ P k ( K ) , (86) where { (cid:98) w µ } L µ = are the weights of the L-point Gauss–Lobatto quadrature. If labeling the bottom,right, top and left edges of K as E , E , E and E , respectively, then the identity (86) implies, for ≤ µ ≤ N, that ϖ ( µ ) E j = ∆ x i (cid:98) ω ω µ ∆ x i + ∆ y (cid:96) , j = , ϖ ( µ ) E j = ∆ y (cid:96) (cid:98) ω ω µ ∆ x i + ∆ y (cid:96) , j = , . According to Theorem 5.3,the CFL condition (83) for our positivity-preserving DG schemes on Cartesian meshes becomes α ∆ t (cid:18) ∆ x i + ∆ y (cid:96) (cid:19) ≤ (cid:98) ω = L ( L − ) . (87) In this section, we present numerical tests on several benchmark RHD problems to validate theaccuracy and effectiveness of our IRP DG methods on 1D and 2D uniform Cartesian meshes. The31hird-order SSP-RK method (73) or SSP-MS method (74) will be employed for time discretization.Unless otherwise stated, we use the ideal EOS (4) with Γ = /
3, and set the CFL numbers as 0 . .
15, 0 .
1, respectively, for the second-order ( P -based), third-order ( P -based), fourth-order ( P -based) DG methods with the SSP-RK time discretization; the CFL numbers for the SSP-MS-DGmethods are one-third of those for the SSP-RK-DG methods.For convenience, we refer to the 1D bound-preserving limiter (68)–(69) (cf. [28]) as the BPlimiter . Our IRP limiter (68)–(70) corresponds to a combination of the BP limiter and the entropylimiter (70). The same names/abbreviations will be also used for those corresponding limiters inthe 2D case. We will compare the results with the proposed IRP limiter and those with only theBP limiter.Table 1:
Example 1: Errors at t = . P k -based DG methods( k = , , SSP-RK SSP-MS k N l error order l error order l error order l error order1 10 1.75e-2 – 1.98e-2 – 1.72e-2 – 1.96e-2 –20 3.16e-3 2.47 4.47e-3 2.14 3.13e-3 2.46 4.41e-3 2.1540 8.19e-4 1.95 1.08e-3 2.05 8.19e-4 1.94 1.07e-3 2.0580 1.93e-4 2.08 2.49e-4 2.11 1.92e-4 2.09 2.46e-4 2.12160 4.63e-5 2.06 5.87e-5 2.08 4.61e-5 2.06 5.84e-5 2.08320 1.13e-5 2.04 1.38e-5 2.09 1.12e-5 2.04 1.37e-5 2.092 10 1.24e-3 – 1.47e-3 – 7.76e-4 3.21 9.15e-4 3.3120 1.83e-4 2.76 2.37e-4 2.63 8.40e-5 3.04 9.24e-5 3.0240 2.75e-5 2.74 4.98e-5 2.25 1.02e-5 3.00 1.14e-5 3.0080 4.06e-6 2.76 1.03e-5 2.28 1.27e-6 3.00 1.42e-6 3.00160 5.90e-7 2.78 2.12e-6 2.28 1.59e-7 3.00 1.77e-7 3.00320 9.16e-8 2.69 4.51e-7 2.23 1.99e-8 3.00 2.22e-8 3.003 10 6.12e-5 – 8.32e-5 – 1.92e-5 – 2.24e-5 –20 4.84e-6 3.66 9.94e-6 3.07 1.29e-6 3.90 1.50e-6 3.9040 3.00e-7 4.01 9.89e-7 3.33 7.86e-8 4.04 8.89e-8 4.0780 2.71e-8 3.47 1.26e-7 2.97 4.85e-9 4.02 5.48e-9 4.02160 2.30e-9 3.55 1.52e-8 3.05 3.04e-10 4.00 3.42e-10 4.00320 2.08e-10 3.47 1.84e-9 3.04 1.90e-11 4.00 2.13e-11 4.0132 .1 Example 1: 1D smooth problem To examine the accuracy of our 1D DG methods we first test a smooth problem similar to [28,41, 45, 46]. The initial conditions are ρ ( x , ) = + . (cid:0) π x (cid:1) , v ( x , ) = . p ( x , ) = [ , ] with periodic boundary conditions, so that the exactsolution is ρ ( x , t ) = + . (cid:0) π ( x − . t ) (cid:1) , v ( x , t ) = . p ( x , t ) =
1. In the computations,the domain is partitioned into N uniform cells with N ∈ { , , , , , } . For the P -based DG method, we take (only in this accuracy test) the time step-sizes as ∆ t = . ∆ x and ∆ t = . ∆ x for the third-order SSP-RK and SSP-MS time discretizations respectively, so as tomatch the fourth-order accuracy of spatial discretization.Table 1 lists the numerical errors at t = . P k -based IRP DG methods ( k = , ,
3) at different grid resolutions.As observed in [46, 28], the accuracy degenerates for SSP-RK and k ≥
2, which is due to thelower order accuracy in the RK intermediate stages as explained in [46]. The desired full order ofaccuracy is observed for the SSP-MS time discretization, indicating that the IRP limiter itself doesnot destroy the accuracy for smooth solutions as expected from the analyses in [44, 43, 18].
This example investigates the capability of the 1D IRP DG methods in resolving discontinuoussolutions, by testing two 1D Riemann problems. The computational domain is taken as [ , ] .The initial conditions of the first Riemann problem are ( ρ , v , p )( x , ) = ( . , . , ) , x < . , ( , , ) , x > . . The initial discontinuity will evolve as a left-moving rarefaction wave, a constant discontinuity,and a right-moving shock wave. Figure 1 displays the numerical solutions computed by the fourth-order DG methods with the BP or IRP limiter respectively on a mesh of 320 uniform cells, againstthe exact solution. Note that, in this simulation, we do not use any other non-oscillatory limiterssuch as the TVD/TVB or WENO limiters. We can observe that the numerical results with only theBP limiter exhibit overshoot in the rest-mass density near the contact discontinuity and some smalloscillations. When the IRP limiter is applied (i.e., the entropy limiter (70) is added), the overshootand oscillations in the DG solution are damped. This is consistent with the observation in [21,46, 18] that enforcing the discrete minimum entropy principle could help to oppress numericaloscillations. Figure 2 shows the time evolution of the minimum specific entropy values of the DGsolutions. It is seen that the minimum remains the same for the DG scheme with the IRP limiter,which indicates that the minimum entropy principle is preserved, while the DG scheme with onlythe BP limiter fails to keep the principle. 33 (a) With the BP limiter (b) With the IRP limiter
Figure 1:
Example 2: Solutions (and their close-up) of the first 1D Riemann problem at t = .
4. Thesymbols “ ◦ ” denote the numerical solutions obtained by the fourth-order DG methods with the BP limiter(left) or with the IRP limiter (right) on the mesh of 320 uniform cells, while the solid lines denote the exactsolution. (Here we do not use any other non-oscillatory limiters, e.g. TVD/TVB or WENO limiters.) with the BP limiterwith the IRP limiter -3 with the BP limiterwith the IRP limiter Figure 2:
The first Riemann problem of Example 2: Time evolution of S min ( t ) for the DG solutions withthe IRP limiter (68)–(70) or with the BP limiter (68)–(69). Left: S min ( t ) = min j , µ S ( U h ( (cid:98) x ( µ ) j , t )) ; right: S min ( t ) = min j S ( U j ( t )) . In order to demonstrate the robustness and resolution of the proposed IRP DG methods, wesimulate a ultra-relativistic Riemann problem [38]. The initial conditions are ( ρ , v , p )( x , ) = ( , , ) , x < . , ( , , − ) , x > . , which involve very low pressure and strong initial jump in pressure ( ∆ p : = | p R − p L | / p R ≈ ),so that the simulation of this problem is challenging and the constraint-preserving or BP techniqueshave to be used [38, 28, 41]. The initial discontinuity will evolve as a strong left-moving rarefac-tion wave, a quickly right-moving contact discontinuity, and a quickly right-moving shock wave.More precisely, the speeds of the contact discontinuity and the shock wave are about 0.986956 and0.9963757 respectively, and are very close to the speed of light c =
1. Figure 3 presents the numer-ical solutions computed by the fourth-order DG methods with the BP or IRP limiter respectivelyon a mesh of 400 uniform cells, against the exact solution. Due to the ultra-relativistic effect, aclearly curved profile for the rarefaction fan is yielded (see Figure 3), as opposed to a linear one inthe non-relativistic case. Again, here we do not use any other non-oscillatory limiters such as theTVD/TVB or WENO limiters. We see that both DG methods work very robustly and exhibit sim-ilar high resolution (the numerical solutions are comparable to those obtained by the ninth-orderbound-preserving finite difference WENO methods in [38]). This implies that the use of entropylimiter in the IRP DG method keeps the robustness and does not destroy the high resolution ofthe scheme. However, without the entropy limiter (i.e., with only the BP limiter), the DG schemewould not preserve the minimum entropy principle and thus is not IRP, as confirmed by the plotsin Figure 4. We also remark that if the BP or IRP limiter is not applied, the DG code would break35own due to nonphysical numerical solutions exceeding the set G . To check the accuracy of our 2D IRP DG methods we simulate a smooth problem from [41]. Theinitial data are ( ρ , vvv , p )( x , y , ) = (cid:0) + . ( π ( x + y )) , . / √ , . / √ , − (cid:1) . The computational domain is taken as [ , ] with periodic boundary conditions, so that the exactsolution is ( ρ , vvv , p )( x , y , t ) = (cid:0) + . ( π ( x + y − . √ t )) , . / √ , . / √ , − (cid:1) , which describes the propagation of an RHD sine wave with low density, low pressure, and largevelocity, in the domain [ , ] at an angle 45 ◦ with the x -axis.In the computations, the domain is partitioned into N × N uniform rectangular cells with N ∈ { , , , , } . For the P -based DG method, we take (only in this accuracy test) thetime step-sizes as ∆ t = . (cid:0) ∆ x (cid:1) and ∆ t = . (cid:0) ∆ x (cid:1) for the third-order SSP-RK and SSP-MS timediscretizations respectively, so as to match the fourth-order accuracy of spatial DG discretization.Table 2 lists the numerical errors at t = . P k -based IRP DG methods ( k = , ,
3) at different grid resolutions. Similar tothe 1D case and as also observed in [46, 28], the accuracy slightly degenerates for SSP-RK and k ≥
2. The desired full order of accuracy is observed for the SSP-MS time discretization, confirm-ing that the IRP limiter itself does not destroy the accuracy for smooth solutions as expected.
This example simulates the interaction between a planar shock wave and a light bubble within thedomain [ , ] × [ − , ] . The setup is the same as in [15, 48]. Initially, a left-moving relativisticshock wave is located at x =
265 with the left and right states given by ( ρ , vvv , p )( x , y , ) = ( , , , . ) , x < , ( . , − . , , . ) , x > . A light circular bubble with the radius of 25 is initially centered at ( , ) in front of the initialshock wave. The fluid state within the bubble is given by ( ρ , vvv , p )( x , y , ) = ( . , , , . ) , (cid:113) ( x − ) + y ≤ . (a) With the BP limiter (b) With the IRP limiter Figure 3:
Example 2: Solutions of the second 1D Riemann problem at t = .
45. The symbols “ ◦ ” denotethe numerical solutions obtained by the fourth-order DG methods with the BP limiter (left) or with the IRPlimiter (right) on the mesh of 400 uniform cells, while the solid lines denote the exact solution. (Here wedo not use any other non-oscillatory limiters, e.g. TVD/TVB or WENO limiters.) with the BP limiterwith the IRP limiter with the BP limiterwith the IRP limiter Figure 4:
The second Riemann problem of Example 2: Time evolution of S min ( t ) for the DG solutionswith the IRP limiter (68)–(70) or with the BP limiter (68)–(69). Left: S min ( t ) = min j , µ S ( U h ( (cid:98) x ( µ ) j , t )) ; right: S min ( t ) = min j S ( U j ( t )) . The reflective conditions are specified at both the top and bottom boundaries { y = ± , ≤ x ≤ } , and the inflow (resp. outflow) boundary condition is enforced at the right (resp. left) bound-ary.In this simulation, we also do not use any other non-oscillatory limiters, e.g. TVD/TVB orWENO limiters. Figure 5 shows the numerical results obtained by the fourth-order DG methods,with the BP limiter or the IRP limiter respectively, on a mesh of 650 ×
180 uniform cells. Weobserve serious oscillations developing near the top and bottom boundaries in the DG solutionwith only the BP limiter. However, if the IRP limiter is used, the undesirable oscillations arealmost oppressed, and the discontinuities and some small wave structures including the curling ofthe bubble interface are well resolved. To check the preservation of minimum entropy principle,we plot the time evolution of the minimum specific entropy values of the DG solutions in Figure 6.It shows that the minimum stays the same for the DG scheme with the IRP limiter, which indicatesthat the minimum entropy principle is maintained, while the DG scheme with only the BP limiterdoes not preserve the principle.
In this test, we simulate two Riemann problems of the 2D RHD equations, which have becomebenchmark tests for checking the accuracy and resolution of 2D RHD codes (cf. [15, 48, 38, 41, 2]).38 a) With the BP limiter (b) With the IRP limiter
Figure 5:
Example 4: Schlieren images of the density at t =
90, 180, 270, 360, and 450 (from top tobottom) obtained by the fourth-order DG methods on the mesh of 650 ×
180 cells. (Here we do not use anyother non-oscillatory limiters, e.g. TVD/TVB or WENO limiters.)
90 180 270 360 450-25-20-15-10-505 with the BP limiterwith the IRP limiter with the BP limiterwith the IRP limiter
Figure 6:
Example 4: Time evolution of S min ( t ) for the DG solutions with the IRP limiter or with the BPlimiter. Left: S min ( t ) = min K ∈ T h min xxx ∈ X K S ( U h ( xxx , t )) ; right: S min ( t ) = min K ∈ T h S ( U K ( t )) . -0.6 -0.2 0.2 0.6-1-0.6-0.20.20.61 (a) With the BP limiter (without any othernon-oscillatory limiters) -0.6 -0.2 0.2 0.6-1-0.6-0.20.20.61 (b) With the IRP limiter (without any othernon-oscillatory limiters) Figure 7:
The first Riemann problem of Example 5: The contours of the rest-mass density logarithm at t = . ×
200 cells. (18 equally spaced contour lines from − . − . Example 3: Errors at t = . P k -based DG methods( k = , , SSP-RK SSP-MS k N l error order l error order l error order l error order1 10 3.95e-2 – 4.80e-2 – 4.02e-2 – 4.80e-2 –20 7.62e-3 2.37 1.01e-2 2.24 7.54e-3 2.41 1.00e-2 2.2640 1.65e-3 2.21 2.25e-3 2.17 1.65e-3 2.20 2.24e-3 2.1580 3.85e-4 2.10 5.41e-4 2.05 3.84e-4 2.10 5.40e-4 2.05160 9.49e-5 2.02 1.34e-4 2.01 9.48e-5 2.02 1.34e-4 2.012 10 1.14e-2 – 1.46e-2 – 1.06e-2 – 1.34e-2 –20 3.90e-4 4.88 5.00e-4 4.87 3.80e-4 4.80 4.89e-4 4.7840 4.89e-5 3.00 6.21e-5 3.01 4.11e-5 3.21 5.47e-5 3.1680 6.55e-6 2.90 8.62e-6 2.85 4.90e-6 3.07 6.72e-6 3.03160 7.65e-7 3.10 1.22e-6 2.83 6.08e-7 3.01 8.38e-7 3.003 10 2.51e-4 – 3.75e-4 – 2.47e-4 – 3.70e-4 –20 1.82e-5 3.79 2.65e-5 3.82 1.81e-5 3.77 2.68e-5 3.7940 9.60e-7 4.24 1.43e-6 4.21 9.02e-7 4.33 1.30e-6 4.3780 6.55e-8 3.87 1.04e-6 3.78 5.56e-8 4.02 7.67e-8 4.08160 4.51e-9 3.86 9.14e-9 3.51 3.46e-9 4.01 4.71e-9 4.03The first Riemann problem was proposed in [38], with the initial data given by ( ρ , vvv , p )( x , y , ) = ( . , , , ) (cid:62) , x > , y > , ( . , . , , . ) (cid:62) , x < , y > , ( . , , , . ) (cid:62) , x < , y < , ( . , , . , . ) (cid:62) , x > , y < , where the left and lower initial discontinuities are contact discontinuities, and the upper and rightare shock waves. Because the maximal value of the fluid velocity is very close to the speed of light( c = not use any other non-oscillatory limiters, e.g. TVD/TVBor WENO limiters. We evolve the solution up to t = . ×
200 cells within thedomain [ − , ] . The contours of log ( ρ ) are displayed in Figure 7, obtained by the fourth-orderDG methods, with the BP or IRP limiter respectively. Serious oscillations are observed in the DGsolution with only the BP limiter, while, if the IRP limiter is applied (i.e., the entropy limiter isadded), the oscillations are much reduced. The minimum entropy principle of the DG solutionwith the IRP limiter is validated in Figure 8. As the numerical solution is preserved in the set41 with the BP limiterwith the IRP limiter with the BP limiterwith the IRP limiter Figure 8:
The first Riemann problem of Example 5: S min ( t ) for the DG solutions with the IRP limiter orwith the BP limiter. Left: S min ( t ) = min K ∈ T h min xxx ∈ X K S ( U h ( xxx , t )) ; right: S min ( t ) = min K ∈ T h S ( U K ( t )) . Ω S , the IRP DG scheme exhibits strong robustness in such ultra-relativistic flow simulation; thecomputed flow structures agree well with those reported in [38, 2]. We remark that if the BP orIRP limiter is turned off, the DG code would break down due to nonphysical numerical solutionsviolating the physical constraints (7).The second Riemann problem was first proposed in [15], and its initial condition is given by ( ρ , vvv , p )( x , y , ) = ( . , , , . ) (cid:62) , x > , y > , ( . , . , , ) (cid:62) , x < , y > , ( . , , , ) (cid:62) , x < , y < , ( . , , . , ) (cid:62) , x > , y < . Figures 9(a) and 9(b) display the contours of log ( ρ ) at t = .
8, obtained by the fourth-order DGmethods, with the BP or IRP limiter respectively and without any non-oscillatory limiters, onthe mesh of 200 ×
200 cells within the domain [ − , ] . For comparison and reference, we alsopresent in Figures 9(c–d) the numerical results obtained by the fourth-order DG method with onlythe WENO limiter [48] on two meshes of 200 ×
200 cells and 400 ×
400 cells, respectively. (TheWENO limiter is applied only within some “trouble” cells adaptively identified by the KXRCFindicator [22].) We can see that the result with only the BP limiter is oscillatory, while enforcingthe minimum entropy principle by the IRP limiter helps to damp some of the oscillations. Althoughthe WENO limiter may completely suppress the undesirable oscillations, the resulting numericalsolution in Figure 9(c) is much more dissipative than that computed with the IRP limiter in Figure9(b). Figure 10 shows the time evolution of the minimum specific entropy values of the DGsolutions. One can see that the minimum remains the same for the DG scheme with the IRP limiter.This implies the preservation of minimum entropy principle, which, however, is not ensured by42 (a) With the BP limiter (without any other non-oscillatory limiters); 200 ×
200 cells -0.6 -0.2 0.2 0.6-1-0.6-0.20.20.61 (b) With the IRP limiter (without any othernon-oscillatory limiters); 200 ×
200 cells -0.6 -0.2 0.2 0.6-1-0.6-0.20.20.61 (c) With the WENO limiter; 200 ×
200 cells -0.6 -0.2 0.2 0.6-1-0.6-0.20.20.61 (d) With the WENO limiter; 400 ×
400 cells
Figure 9:
The second Riemann problem of Example 5: The contours of the rest-mass density logarithm at t = . − . − .
426 are displayed.) with the BP limiterwith the IRP limiter with the BP limiterwith the IRP limiter Figure 10:
The second Riemann problem of Example 5: S min ( t ) for the DG solutions with the IRP limiteror with the BP limiter. Left: S min ( t ) = min K ∈ T h min xxx ∈ X K S ( U h ( xxx , t )) ; right: S min ( t ) = min K ∈ T h S ( U K ( t )) . using either only the BP limiter or the WENO limiter.From the above numerical results, we observe that the IRP limiter helps to preserve numericalsolutions in the invariant region Ω S and to damp some numerical oscillations while keeping thehigh resolution. However, as pointed out in [18], the IRP limiter is still a mild limiter and may notcompletely suppress all the oscillations. For problems involving strong shocks, one may use theIRP limiter together with a non-oscillatory limiter (e.g., the WENO limiter), to sufficiently controlall the undesirable oscillations while also preserving the invariant region Ω S . We now further examine the robustness and the IRP property of our DG method by simulating tworelativistic jets. For high-speed jet problems, the internal energy is exceedingly small compared tothe kinetic energy so that negative numerical pressure can be generated easily in the simulations.Since shear waves, strong relativistic shock waves, and ultra-relativistic regions are involved inthis challenging test, we use the IRP or BP limiter together with the WENO limiter to preserve thephysical constraints (7) and suppress all the undesirable oscillations.We first simulate a cold (highly supersonic) jet model from [41]. Initially, a RHD jet (density ρ b = .
1; speed v b = .
99; classical Mach number M b =
50) is injected along the y -direction intothe domain [ − , ] × [ , ] , which is initially filled with a static uniform medium of density ρ a =
1. The fixed inflow condition is enforced at the jet nozzle { ( x , y ) : | x | ≤ . , y = } onthe bottom boundary, while the other boundary conditions are outflow. For this problem, the jetpressure matches that of the ambient medium, the corresponding initial Lorentz factor W ≈ . M r : = M b W / W s ≈ .
37, where W s denotes the Lorentz factorof the acoustic wave speed. The exceedingly high Mach number and jet speed cause the simulation44igure 11: Example 6: The pseudocolor plots of log ( ρ ) for the cold relativistic jet at t = ,
20, and 30(from left to right), respectively. of this problem difficult, so that the BP or constraint-preserving technique is required for high-order schemes to keep the numerical solutions in the invariant region. We set the computationaldomain as [ , ] × [ , ] with the reflective boundary condition at x =
0, and uniformly dividedthe domain into 240 ×
500 cells. The results, simulated by the fourth-order IRP DG method, arepresented in Figure 11. Those plots clearly shows the jet evolution, and the flow structures at thefinal time t =
30 are in good agreement with those computed in [41].We then simulate a pressure-matched hot jet model. Different from [41], the ideal EOS (4)with Γ = is used here. The setups are the same as the above clod jet model, except for thatthe density of the inlet jet becomes ρ b = .
01, and the classical Mach number is set as M b = .
72 and near the minimum Mach number M min = v b / √ Γ − [ , ] × [ , ] and divided into 240 ×
600 uniform cells. Figure 12 gives the numericalresult computed the fourth-order IRP DG method. One can see that the flow structures are clearlyresolved by our method, and the patterns are different from those in the cold jet model.As we have seen, the IRP DG method is very robust in such demanding jet simulations. Tofurther verify the minimum entropy principle, we plot in Figure 13 the time evolution of the min-imum specific entropy values. It is observed that the minimum does not change with time for theDG scheme with the WENO and IRP limiters, which means that the minimum entropy principleis preserved, while using the WENO and BP limiters is unable to maintain the principle.45igure 12:
Example 6: The pseudocolor plots of log ( ρ ) for the hot relativistic jet at t = ,
22, and 33(from left to right), respectively. with WENO and IRP limiterswith WENO and BP limiters (a) for the cold jet with WENO and IRP limiterswith WENO and BP limiters (b) for the hot jet
Figure 13:
Example 6: Evolution of S min ( t ) : = min K ∈ T h min xxx ∈ X K S ( U h ( xxx , t )) for the DG solutions with theIRP limiter or the BP limiter. The WENO limiter is applied right before the IRP or BP limiting procedure. Conclusions
In this work, we showed that the Tadmor’s minimum entropy principle (11) holds for the RHDequations (1) with the ideal EOS (4), and then developed high-order accurate IRP DG and finitevolume schemes for RHD, which provably preserve a discrete minimum entropy principle as wellas the intrinsic physical constraints (7). It was the first time that such a minimum entropy principlewas explored for RHD at the continuous and discrete levels. Due to the relativistic effects, the spe-cific entropy is a highly nonlinear function of the conservative variables and cannot be explicitlyexpressed. This led to some essential difficulties in this work, which were not encountered in thenon-relativistic case. In order to address the difficulties, we first proposed a novel equivalent formof the invariant region. As a notable feature, all the constraints in this novel form became explicitand linear with respect to the conservative variables. This provided a highly effective approach totheoretically analyze the IRP property of numerical RHD schemes. We showed the convexity ofthe invariant region and established the generalized Lax–Friedrichs splitting properties via techni-cal estimates. We rigorously proved that the first-order Lax–Friedrichs type scheme for the RHDequations satisfies a local minimum entropy principle and is IRP under a CFL condition. We thendeveloped and analyzed provably IRP high-order accurate DG and finite volume methods for theRHD. Numerical examples confirmed that enforcing the minimum entropy principle could help todamp some undesirable numerical oscillations and demonstrated the robustness of the proposedhigh-order IRP DG schemes. The verified minimum entropy principle is an important propertythat can be incorporated into improving or designing other RHD schemes. Besides, the proposednovel analysis techniques can be useful for investigating or seeking other IRP schemes for theRHD or other physical systems.
References [1] D. S. B
ALSARA AND
J. K IM , A subluminal relativistic magnetohydrodynamics scheme withADER-WENO predictor and multidimensional Riemann solver-based corrector , Journal ofComputational Physics, 312 (2016), pp. 357–384.[2] D. B
HORIYA AND
H. K
UMAR , Entropy-stable schemes for relativistic hydrodynamics equa-tions , Zeitschrift f¨ur angewandte Mathematik und Physik, 71 (2020), pp. 1–29.[3] J. D
UAN AND
H. T
ANG , High-order accurate entropy stable finite difference schemes forone-and two-dimensional special relativistic hydrodynamics , Advances in Applied Mathe-matics and Mechanics, 12 (2020), pp. 1–29.[4] J. A. F
ONT , Numerical hydrodynamics and magnetohydrodynamics in general relativity ,Living Reviews in Relativity, 11 (2008), p. 7.475] S. G
OTTLIEB , C.-W. S HU , AND
E. T
ADMOR , Strong stability-preserving high-order timediscretization methods , SIAM Review, 43 (2001), pp. 89–112.[6] A. G
OUASMI , K. D
URAISAMY , S. M. M
URMAN , AND
E. T
ADMOR , A minimum entropyprinciple in the compressible multicomponent Euler equations , ESAIM: Mathematical Mod-elling and Numerical Analysis, 54 (2020), pp. 373–389.[7] F. G
UERCILENA , D. R
ADICE , AND
L. R
EZZOLLA , Entropy-limited hydrodynamics: anovel approach to relativistic hydrodynamics , Computational Astrophysics and Cosmology,4 (2017), p. 3.[8] J.-L. G
UERMOND , M. N
AZAROV , B. P
OPOV , AND
I. T
OMAS , Second-order invariant do-main preserving approximation of the Euler equations using convex limiting , SIAM Journalon Scientific Computing, 40 (2018), pp. A3211–A3239.[9] J.-L. G
UERMOND AND
B. P
OPOV , Viscous regularization of the Euler equations and en-tropy principles , SIAM J. Appl. Math., 74 (2014), pp. 284–305.[10] J.-L. G
UERMOND AND
B. P
OPOV , Invariant domains and first-order continuous finite el-ement approximation for hyperbolic systems , SIAM Journal on Numerical Analysis, 54(2016), pp. 2466–2489.[11] J.-L. G
UERMOND AND
B. P
OPOV , Invariant domains and second-order continuous finiteelement approximation for scalar conservation equations , SIAM Journal on Numerical Anal-ysis, 55 (2017), pp. 3120–3146.[12] J.-L. G
UERMOND , B. P
OPOV , L. S
AAVEDRA , AND
Y. Y
ANG , Invariant domains preservingarbitrary Lagrangian Eulerian approximation of hyperbolic systems with continuous finiteelements , SIAM Journal on Scientific Computing, 39 (2017), pp. A385–A414.[13] J.-L. G
UERMOND , B. P
OPOV , AND
I. T
OMAS , Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems , Computer Methods in Ap-plied Mechanics and Engineering, 347 (2019), pp. 143–175.[14] A. H
ARTEN , On the symmetric form of systems of conservation laws with entropy , Journalof Computational Physics, 49 (1983), pp. 151–164.[15] P. H
E AND
H. T
ANG , An adaptive moving mesh method for two-dimensional relativistichydrodynamics , Communications in Computational Physics, 11 (2012), pp. 114–146.[16] X. Y. H U , N. A. A DAMS , AND
C.-W. S HU , Positivity-preserving method for high-order conservative schemes solving compressible Euler equations , Journal of ComputationalPhysics, 242 (2013), pp. 169–180. 4817] Y. J
IANG AND
H. L IU , An invariant-region-preserving (IRP) limiter to DG methods forcompressible Euler equations , in XVI International Conference on Hyperbolic Problems:Theory, Numerics, Applications, Springer, 2016, pp. 71–83.[18] Y. J
IANG AND
H. L IU , Invariant-region-preserving DG methods for multi-dimensional hy-perbolic conservation law systems, with an application to compressible Euler equations ,Journal of Computational Physics, 373 (2018), pp. 385–409.[19] Y. J
IANG AND
H. L IU , Invariant-region-preserving DG methods for multi-dimensional hy-perbolic conservation law systems, with an application to compressible Euler equations ,Journal of Computational Physics, 373 (2018), pp. 385–409.[20] Y. J
IANG AND
H. L IU , An invariant-region-preserving limiter for DG schemes to isentropicEuler equations , Numerical Methods for Partial Differential Equations, 35 (2019), pp. 5–33.[21] B. K
HOBALATTE AND
B. P
ERTHAME , Maximum principle on the entropy and second-orderkinetic schemes , Mathematics of Computation, 62 (1994), pp. 119–131.[22] L. K
RIVODONOVA , J. X IN , J.-F. R EMACLE , N. C
HEVAUGEON , AND
J. E. F
LAHERTY , Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conserva-tion laws , Applied Numerical Mathematics, 48 (2004), pp. 323–338.[23] Y. L
V AND
M. I
HME , Entropy-bounded discontinuous Galerkin scheme for Euler equations ,Journal of Computational Physics, 295 (2015), pp. 715–739.[24] J. M. M
ART ´ I AND
E. M ¨
ULLER , Numerical hydrodynamics in special relativity , Living Re-views in Relativity, 6 (2003), p. 7.[25] J. M. M
ART ´ I AND
E. M ¨
ULLER , Grid-based methods in relativistic hydrodynamics and mag-netohydrodynamics , Living Reviews in Computational Astrophysics, 1 (2015), p. 3.[26] V. M
EWES , Y. Z
LOCHOWER , M. C
AMPANELLI , T. W. B
AUMGARTE , Z. B. E
TIENNE ,F. G. L. A
RMENGOL , AND
F. C
IPOLLETTA , Numerical relativity in spherical coordinates:A new dynamical spacetime and general relativistic MHD evolution framework for the Ein-stein Toolkit , Physical Review D, 101 (2020), p. 104007.[27] A. M
IGNONE AND
G. B
ODO , An HLLC Riemann solver for relativistic flows–I. Hydrody-namics , Monthly Notices of the Royal Astronomical Society, 364 (2005), pp. 126–136.[28] T. Q IN , C.-W. S HU , AND
Y. Y
ANG , Bound-preserving discontinuous Galerkin methods forrelativistic hydrodynamics , Journal of Computational Physics, 315 (2016), pp. 323–347.4929] D. R
ADICE , L. R
EZZOLLA , AND
F. G
ALEAZZI , High-order fully general-relativistic hydro-dynamics: new approaches and tests , Classical and Quantum Gravity, 31 (2014), p. 075012.[30] E. T
ADMOR , Skew-selfadjoint form for systems of conservation laws , J. Math. Anal. Appl.,103 (1984), pp. 428–442.[31] E. T
ADMOR , A minimum entropy principle in the gas dynamics equations , Applied Numeri-cal Mathematics, 2 (1986), pp. 211–219.[32] K. W U , Design of provably physical-constraint-preserving methods for general relativistichydrodynamics , Physical Review D, 95 (2017), 103001.[33] K. W U , Positivity-preserving analysis of numerical schemes for ideal magnetohydrodynam-ics , SIAM Journal on Numerical Analysis, 56 (2018), pp. 2124–2147.[34] K. W
U AND
C.-W. S HU , A provably positive discontinuous Galerkin method for multidi-mensional ideal magnetohydrodynamics , SIAM Journal on Scientific Computing, 40 (2018),pp. B1302–B1329.[35] K. W
U AND
C.-W. S HU , Provably positive high-order schemes for ideal magnetohydrody-namics: analysis on general meshes , Numerische Mathematik, 142 (2019), pp. 995–1047.[36] K. W
U AND
C.-W. S HU , Entropy symmetrization and high-order accurate entropy stablenumerical schemes for relativistic MHD equations , SIAM Journal on Scientific Computing,42 (2020), pp. A2230–A2261.[37] K. W
U AND
C.-W. S HU , Provably physical-constraint-preserving discontinuous Galerkinmethods for multidimensional relativistic MHD equations , arXiv preprint arXiv:2002.03371,(2020).[38] K. W
U AND
H. T
ANG , High-order accurate physical-constraints-preserving finite differenceWENO schemes for special relativistic hydrodynamics , Journal of Computational Physics,298 (2015), pp. 539–564.[39] K. W
U AND
H. T
ANG , A direct Eulerian GRP scheme for spherically symmetric generalrelativistic hydrodynamics , SIAM Journal on Scientific Computing, 38 (2016), pp. B458–B489.[40] K. W
U AND
H. T
ANG , Admissible states and physical-constraints-preserving schemes forrelativistic magnetohydrodynamic equations , Mathematical Models and Methods in AppliedSciences, 27 (2017), pp. 1871–1928. 5041] K. W
U AND
H. T
ANG , Physical-constraint-preserving central discontinuous Galerkin meth-ods for special relativistic hydrodynamics with a general equation of state , The AstrophysicalJournal Supplement Series, 228 (2017), 3.[42] Z. X U , Parametrized maximum principle preserving flux limiters for high order schemessolving hyperbolic conservation laws: one-dimensional scalar problem , Mathematics ofComputation, 83 (2014), pp. 2213–2238.[43] X. Z
HANG , On positivity-preserving high order discontinuous Galerkin schemes for com-pressible Navier-Stokes equations , Journal of Computational Physics, 328 (2017), pp. 301–343.[44] X. Z
HANG AND
C.-W. S HU , On maximum-principle-satisfying high order schemes forscalar conservation laws , Journal of Computational Physics, 229 (2010), pp. 3091–3120.[45] X. Z
HANG AND
C.-W. S HU , On positivity-preserving high order discontinuous Galerkinschemes for compressible Euler equations on rectangular meshes , Journal of ComputationalPhysics, 229 (2010), pp. 8918–8934.[46] X. Z
HANG AND
C.-W. S HU , A minimum entropy principle of high order schemes for gasdynamics equations , Numerische Mathematik, 121 (2012), pp. 545–563.[47] X. Z
HANG , Y. X IA , AND
C.-W. S HU , Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangularmeshes , Journal of Scientific Computing, 50 (2012), pp. 29–62.[48] J. Z
HAO AND
H. T