# High-order gas-kinetic scheme for radiation hydrodynamics in equilibrium-diffusion limit

In this paper, a high-order gas-kinetic scheme is developed for the equation of radiation hydrodynamics in equilibrium-diffusion limit which describes the interaction between matter and radiation. To recover RHE, the Bhatnagar-Gross-Krook (BGK) model with modified equilibrium state is considered. In the equilibrium-diffusion limit, the time scales of radiation diffusion and hydrodynamic part are different, and it will make the time step very small for the fully explicit scheme. An implicit-explicit (IMEX) scheme is applied, in which the hydrodynamic part is treated explicitly and the radiation diffusion is treated implicitly. For the hydrodynamics part, a time dependent gas distribution function can be constructed by the integral solution of modified BGK equation, and the time dependent numerical fluxes can be obtained by taking moments of gas distribution function. For the radiation diffusion term, the nonlinear generalized minimal residual (GMRES) method is used. To achieve the temporal accuracy, a two-stage method is developed, which is an extension of two-stage method for hyperbolic conservation law. For the spatial accuracy, the multidimensional weighted essential non-oscillation (WENO) scheme is used for the spatial reconstruction. A variety of numerical tests are provided for the performance of current scheme, including the order of accuracy and robustness.

## Authors

• 2 publications
• 19 publications
• 2 publications
10/12/2020

### A compact high-order gas-kinetic scheme on unstructured mesh for acoustic and shock wave computations

In this paper an even higher-order compact GKS up to sixth order of accu...
02/28/2022

### Three-dimensional discontinuous Galerkin based high-order gas-kinetic scheme and GPU implementation

In this paper, the discontinuous Galerkin based high-order gas-kinetic s...
03/17/2021

### Stability of a numerical scheme for methane transport in hydrate zone under equilibrium and non-equilibrium conditions

In this paper we carry out numerical analysis for a family of simplified...
03/17/2022

### Three-dimensional third-order gas-kinetic scheme on hybrid unstructured meshes for Euler and Navier-Stokes equations

In this paper, a third order gas kinetic scheme is developed on the thre...
02/03/2022

### Anomalous sorption kinetics of self-interacting particles by a spherical trap

In this paper we propose a computational framework for the investigation...
05/13/2020

### Nonlinear Fokker-Planck Acceleration for Forward-Peaked Transport Problems in Slab Geometry

This paper introduces a nonlinear acceleration technique that accelerate...
10/07/2021

### Generic tool for numerical simulation of transformation-diffusion processes in complex volume geometric shapes: application to microbial decomposition of organic matter

This paper presents a generic framework for the numerical simulation of ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The equation of radiation hydrodynamics describes the radiative transport through a fluid with coupled momentum and energy exchange Radiation-1 ; Radiation-2 ; Radiation-3 . Its applications are mainly in high-temperature hydrodynamics, including gaseous stars in astrophysics, supernova explosions, combustion phenomena, reentry vehicles, fusion physics and inertial confinement fusion. The importance of thermal radiation increases as the temperature is raised in the above problems. Such as for the moderate temperature, the role of radiation is primarily one of transporting energy by radiative process. But for the higher temperature, the energy and momentum densities of the radiation field may become comparable to or even dominates the corresponding fluid quantities.

In the last decades, the gas-kinetic scheme (GKS) and based on the Bhatnagar-Gross-Krook (BGK) model BGK-1 ; BGK-2 have been developed systematically for the computations from low speed flows to supersonic ones GKS-Xu1 ; GKS-Xu2 . The gas-kinetic scheme is based on an analytical integral solution of the BGK equation, and gas distribution function at a cell interface provides a multi-scale evolution process from the kinetic particle transport to the hydrodynamic wave propagation. With the two-stage fourth-order method for Lax-Wendroff type flow solvers GRP-high-1 ; GRP-high-2 , the high-order gas-kinetic schemes were developed GKS-high-1 ; GKS-high-2 . The high-order scheme not only reduces the complexity of computation, but also improves the accuracy of the numerical solution. Most importantly, the robustness is as good as the second-order shock capturing scheme. Furthermore, with the discretization of particle velocity space, a unified gas-kinetic scheme (UGKS) has been developed for the flow study in entire Knudsen number regimes from rarefied to continuum ones UGKS-Xu1 ; UGKS-Xu2 ; UGKS-Xu3 . Recently, the UGKS is extended to solve radiative transfer system with both scattering and absorption/emission effects radiative-UGKS-1 ; radiative-UGKS-2 ; radiative-UGKS-3 . The asymptotic preserving (AP) property can be accurately recovered. For the equation of radiation hydrodynamics, a multi-scale scheme is developed, in which GKS is used for the compressible inviscid flow and UGKS is used for the non-equilibrium radiative transfer radiative-UGKS-4 . Due to the possible large variation of fluid opacity in different regions, the transport of photons through the flow system is simulated by the multi-scale scheme.

In this paper, a high-order gas-kinetic scheme is proposed for the equation of radiation hydrodynamic in the equilibrium-diffusion limit. Based on the zeroth-order Chapman-Enskog expansion, the hydrodynamic part of radiation hydrodynamic equation can be obtained from the modified BGK equation with modified equilibrium state. The radiation diffusion term is considered as source term. Since the time scales of radiation diffusion and hydrodynamic part are different and it will make the time step of an explicit scheme very small. Thus, an implicit-explicit (IMEX) scheme is developed for solving the radiation hydrodynamic equation, in which the fluid advection term is treated explicitly and the radiation diffusion is treated implicitly. For the hydrodynamic part, the gas-kinetic solver with the modified equilibrium state is used to solve the compressible flow equations. The nonlinear Newton-GMRES method GMRES-1 ; GMRES-2 is used to deal with the radiation diffusion. To achieve the temporal accuracy, the two-stage third-order temporal discretization is developed, which is an extension of two-stage method was developed for Lax-Wendroff type flow solvers. To achieve the spatial accuracy, the classical weighted essentially non-oscillatory (WENO) WENO-JS ; WENO-Z method reconstruction is used. With the two-stage temporal discretization and WENO reconstruction, a reliable framework was provided for equation of radiation hydrodynamics. Various numerical experiments are carried out to validate the performance of current scheme.

This paper is organized as follows. In Section 2, the equation of radiation hydrodynamics and corresponding BGK model are introduced. The high-order gas-kinetic for RHE is presented in Section 3. Numerical examples are included in Section 4 and the last section is the conclusion.

## 2 Equation of radiation hydrodynamics and BGK model

### 2.1 Equation of radiation hydrodynamics

The equation of radiation hydrodynamics (RHE) describes the motion of flows under a radiation field. It consists of the Euler equations coupling with the radiation momentum and energy sources and the radiation-transport equation, and can be given as follows

 ∂ρ∂t +∇⋅(ρU)=0, ∂ρU∂t +∇⋅(ρUU+p)=−Srp, ∂E∂t +∇⋅((E+p)U)=−Sre, 1c∂Iν∂t +Ω⋅∇Iν=Qν,

where , , and are the density, velocity, total energy and pressure of the matter, respectively. is the radiation momentum source, is the radiation energy source, is the speed of light, is radiation intensity, is the photon direction of flight and is the angle and frequency dependent radiation source representing the radiation-matter interaction. The radiation-transport equation is essentially the conservation of the photon number, which reveals the relationship between the photon free transport and radiation-matter interaction, i.e. photon emission, photon absorption and photon scatter. The source terms and can be written as the zeroth and first order frequency-integrated angular moments of , respectively,

 Sre≡∫4π∫∞0QνdΩdν=∂ε∂t+∇⋅F, Srp≡1c∫4π∫∞0ΩQνdΩdν=1c2∂F∂t+∇⋅P,

 ε =aRT4, F =−κ∇T4+43UaRT4, P =13aRT4,

where is temperature, is diffusion constant and is radiation constant representing the ratio of the radiation energy to the material energy. The time-derivative of radiation flux in the total-momentum can be dropped in accordance with the diffusion approximation. Therefore, the equation of radiation hydrodynamics in the equilibrium-diffusion limit can be written as

 ∂ρ∂t+∇⋅(ρU)=0,∂ρU∂t+∇⋅(ρUU+p∗)=0,∂E∗∂t+∇⋅((E∗+p∗)U)=∇⋅(κ∇T4), (1)

where the total energy and the total pressure are given as

 E∗= 12ρU2+pγ−1+aRT4, p∗= p+13aRT4.

The polytropic ideal gas is considered, and the equation of state is given by

 p=(γ−1)ρe=(γ−1)ρcvT,

where is the specific internal energy, is the heat capacity at constant volume and is the specific heat ratio.

### 2.2 BGK model

In this paper, a high-order gas-kinetic scheme will be presented for two-dimensional flows. To recover the macroscopic equation of radiation hydrodynamics Eq.(1), the modified BGK equation BGK-1 ; BGK-2 can be written as

 ft+ufx+vfy=g−fτ, (2)

where is the gas distribution function, is the equilibrium distribution, is the particle velocity and is the collision time. For the equation of radiation hydrodynamics, a modified equilibrium state function is introduced radiative-GKS-2

 g(x,u,t)=ρ(λ1λ2(λ1+λ2)π)d/2(λ1π)K1/2(λ2π)K2/2e−(λ1λ2λ1+λ2(u−U)2+λ1ξ21+λ2ξ22),

where for two dimensional system, is the macroscopic velocity, the internal variables are defined as and

and the internal degrees of freedom of

and satisfy

 K1+d=2/(γ−1),  K2+d=6.

The parameters can be given by

 ρ2λ1=p,  ρ2λ2=13aRT4.

The collision term also satisfies the compatibility condition

 ∫ψg−fτdΞ=0,

where and .

According to the Chapman-Enskog expansion, the Euler and Navier-Stokes equations can be derived form the BGK equation GKS-Xu1 ; GKS-Xu2 . Similarly, the macroscopic equations Eq.(1) can be also derived from Eq.(2). Taking zeroth-order Chapman-Enskog expansion, i.e. with , and taking moments of BGK equation Eq.(2), we have

 ∫ψ(gt+ugx+vgy)dΞ=0,

where

 Q=∫ψgdΞ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ρρUρV12ρ(U2+V2+K1+d2λ1+K2+d2λ2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

and

 F(Q)=∫ψugdΞ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ρUρU2+ρ2λ1+ρ2λ2ρUV12ρU(U2+V2+K1+d+22λ1+K2+d+22λ2)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

According to the definition of and the relation of , the hyperbolic part of macroscopic equations Eq.(1

) can be recovered and the vector form is used in the following sections

 ∂Q∂t+∂F(Q)∂x+∂G(Q)∂y=S(Q),

where is the conservative variable, and are the fluxes in and directions and is the source term for radiation diffusion.

## 3 High-order gas-kinetic scheme

### 3.1 Temporal discretization

Recently, based on the time-dependent flux function of the generalized Riemann problem solver (GRP) GRP-high-1 ; GRP-high-2 and gas-kinetic scheme (GKS) GKS-high-1 ; GKS-high-2 , a two-stage fourth-order time-accurate discretization was developed for Lax-Wendroff type flow solvers, particularly applied for the hyperbolic conservation laws. Considering the following time-dependent equation with the initial condition

 dQijdt=L(Qij),

where is an operator for spatial derivative of flux. Introducing an intermediate state at , the two-stage temporal discretization can be written as

 Q∗ij=Qnij+12ΔtL(Qnij)+18Δt2∂∂tL(Qnij),Qn+1ij=Qnij+ΔtL(Qnij)+16Δt2(∂∂tL(Qnij)+2∂∂tL(Q∗ij)). (3)

It can be proved that for the hyperbolic equations the two-stage time stepping method Eq.(3) provides a fourth-order time accurate solution for at . Based on the high-order spatial reconstruction WENO-JS ; WENO-Z , successes have also been achieved for the construction of high-order gas-kinetic scheme for Euler and Navier-Stokes equations GKS-high-1 ; GKS-high-2 . The two-stage method provides a reliable framework to develop high-order scheme with the implementation of second-order flux function. Most importantly, due to the use of both flux function and its temporal derivative, this scheme is robust and works perfectly from the subsonic to the hypersonic flows.

In this paper, the high-order gas-kinetic scheme will be developed for the radiation hydrodynamics equation Eq.(1) as well. For simplicity, the two-dimensional uniform mesh is used. Taking moments of the BGK equation Eq.(2) and integrating with respect to space for the cell , the semi-discrete finite volume scheme can be written as

 dQijdt =L(Qij)+S(Qij), (4)

where is the cell averaged conservative variables over the cell . The operator for hydrodynamic part is given by

 L(Qij)=−1Δx(Fi+1/2,j(t)−Fi−1/2,j(t))−1Δy(Gi,j+1/2(t)−Gi,j−1/2(t)), (5)

and the operator is given by

 S(Qij)=1ΔxΔyS(Qij),

where and are the cell size, and are the time dependent numerical fluxes at cell interfaces in and directions and is the source for radiation diffusion. Without considering the source term, the explicit scheme is used for the Euler and Navier-Stokes equations GKS-high-1 . However, the time scales of radiation diffusion and fluid advection are different. It will make the time step of an explicit scheme very small, and the explicit scheme can be only used efficiently for the mildly non-relativistic regime for RHE. To improve the efficiency, the RHE should be discretized in an implicit-explicit (IMEX) procedure, i.e. the fluid advection term is treated explicitly and the radiation component is treated implicitly. Similar with Eq.(3), introducing an intermediate state at , the two-stage temporal discretization for Eq.(4) is given as follows

 Q∗ij=Qnij+Δt2L(Qnij)+Δt28Lt(Qnij)+Δt4(S(Qnij)+S(Q∗ij)),Qn+1ij=Qnij+ΔtL(Qnij)+Δt26(Lt(Qnij)+2Lt(Q∗ij))+Δt6(S(Qnij)+4S(Q∗ij)+S(Qn+1ij)). (6)

Compared with the original two-stage method Eq.(3), the trapezoid integration is used for the source terms.

It can be proved that Eq.(6) provides a third-order accurate approximation for the system with sources Eq.(4) at . Integrating Eq.(4) on the time interval , we have

 Qn+1ij−Qnij =∫tn+Δttn(L+S)(Qij(t))dt.

To prove this proposition above, the following Taylor expansion need to be satisfied

 ∫tn+Δttn(L+S)(Qij(t))dt= Δt(L+S)(Qnij)+Δt22(L+S)t(Qnij) (7) + Δt36(L+S)tt(Qnij)+O(Δt4).

According to Eq.(4) and Cauchy-Kovalevskaya method, the temporal derivatives can be given by

 Lt=LQ(L+S), St=SQ(L+S),

and

For the operator , we have the following expansion up to the corresponding order

 S(Q∗ij)=S(Qnij)+SQ(Q∗ij−Qnij)+SQQ2(Q∗ij−Qnij)2+O(Q∗ij−Qnij)3,S(Qn+1ij)=S(Qnij)+SQ(Qn+1ij−Qnij)+SQQ2(Qn+1ij−Qnij)2+O(Qn+1ij−Qnij)3. (8)

Substituting Eq.(8) into Eq.(6), we have

 Q∗ij−Qnij= = Δt2(L(Qnij)+S(Qnij))+Δt28(Lt(Qnij)+St(Qnij))+O(Δt3), (9)

and

 Qn+1ij−Qnij= Δt(L(Qnij)+S(Qnij))+Δt26(Lt(Qnij)+2Lt(Q∗ij))) + Δt6(SQ(Qn+1ij−Qnij)+4SQ(Q∗ij−Qnij)) + Δt12(SQQ(Qn+1ij−Qnij)2+4SQQ(Q∗ij−Qnij)2). (10)

For the operator , we have the following expansion up to the corresponding order as well

 L(Q∗ij)=L(Qnij)+LQ(Q∗ij−Qnij)+LQQ2(Q∗ij−Qnij)2+O(Q∗ij−Qnij)2,LQ(Q∗ij)=LQ(Qnij)+LQQ(Q∗ij−Qnij)+LQQQ2(Q∗ij−Qnij)2+O(Q∗ij−Qnij)2. (11)

Substituting Eq.(3.1) and Eq.(11) into Eq.(3.1), it is easy to verify

 Qn+1ij−Qnij =Δt(L+S)(Qnij)+Δt22(LQ+SQ)(L+S)(Qnij) +Δt36((LQQ+SQQ)(L+S)2+(LQ+SQ)2(L+S))(Qnij)+O(Δt4).

Therefore, the two-stage method Eq.(6) provides a third-order temporal discretization for radiation hydrodynamic equations.

### 3.2 Discretization for hydrodynamic part

In the following subsections, the implementation of the hydrodynamic part and radiative part will be given for Eq.(6). For the hydrodynamic part, the numerical flux in -direction can be given by Gaussian quadrature

 Fi+1/2,j(t)=1Δy∫yj+1/2yj−1/2Fi+1/2(y,t)dy=2∑ℓ=1ωℓ∫ψuf(xi+1/2,yjℓ,t,u,v,ξ)dΞ, (12)

where is the Gaussian quadrature point and are quadrature weights. To construct the gas distribution function at the cell interface, the integral solution of BGK equation Eq.(2) is used

 f(xi+1/2,yjℓ,t,u,v,ξ)=1τ∫t0g(x′,y′,t′,u,v,ξ)e−(t−t′)/τdt′+e−t/τf0(−ut,−vt,u,v,ξ), (13)

where is the initial gas distribution function, is the corresponding equilibrium state, and and are the trajectory of particles. Similar with the gas-kinetic scheme for Euler and Navier-Stokes equations, the second-order gas-kinetic solver GKS-Xu2 can be written as follows

 f(xi+1/2,yjℓ,t,u,v,ξ)= (1−e−t/τ)g0+((t+τ)e−t/τ−τ)(¯¯¯au+¯¯bv)g0 + (t−τ+τe−t/τ)¯Ag0 + e−t/τgr[1−(τ+t)(aru+brv)−τAr)](1−H(u)) + e−t/τgl[1−(τ+t)(alu+blv)−τAl)]H(u). (14)

The coefficients in Eq.(3.2) can be determined by the reconstructed derivatives and compatibility condition

 ⟨al,r⟩=∂Ql,r∂x,⟨bl,r⟩=∂Ql,r∂y,⟨al,ru+bl,rv+Al,r⟩=0,

and

 ⟨¯¯¯a⟩=∂Q0∂x,⟨¯¯b⟩=∂Q0∂y,⟨¯¯¯au+¯¯bv+¯¯¯¯A⟩=0,

where the moments of the equilibrium are defined by

 ⟨...⟩=∫g(...)ψdΞ.

Compared with the Euler and Navier-Stokes equations, the procedure for RHE is a little more complicated. As an example, the spatial derivative of can be written as

 ∂g∂x=ag, (15)

where

 a=a⋅~ψ=a1+a2u+a3v+12a4(u2+v2)+12a5ξ21+12a6ξ22.

Taking moments of Eq.(15), we have

 ∂Q∂x=∫∂g∂xψ% dΞ=∫agψdΞ,

To obtain the connections , the equation can be given as a linear system

 Ma=∂Q∂x, (16)

where

 M=∫ψ⊗~ψgdΞ.

The detailed formulation of can be found in Appendix. The system Eq.(16) seems to be under determined, but we can also find the unique solution.

According to the chain rule, the derivative of

can be also written as

 ∂g∂x=∂g∂ρ∂ρ∂x+∂g∂U∂U∂x+∂g∂V∂V∂x+∂g∂λ1∂λ1∂x+∂g∂λ2∂λ2∂x.

It’s easily to verify that

 ∂g∂λ1=C1g−ξ21g, ∂g∂λ2=C2g−ξ22g,

where and are the functions without and , respectively. Comparing the coefficients of and with Eq.(15), we have

 12a5ξ21g=−ξ21g∂λ1∂x, 12a6ξ22g=−ξ22g∂λ2∂x.

According to the definition of and , we have

 a5 =−2∂λ1∂x=TxRT2, a6 =−2∂λ2∂x=−ρxT−4ρTxaRT5/3,

where

 Tx =1ρcv+4aRT3((E∗x−(ρU)(ρU)x+(ρV)(ρV)xρ)+((ρU)2+(ρV)2)ρx2ρ2−ρxcvT).

With the solution for and , the unique solution of the linear system is given by

 a4 =(1ρ(A−UB−VC)−K18λ21a5−K28λ22a6)/L2, a3 =1ρLC−Va4, a2 =1ρLB−Ua4, a1 =1ρ∂ρ∂x−Ua2−Va3−B′1a4−K14λ1a5−K24λ2a6,

where

 A =∂E∗∂x−12(U2+V2+K1+22λ1+K2+22λ2)∂ρ∂x, B =∂(ρU)∂x−U∂ρ∂x, C =∂(ρV)∂x−V∂ρ∂x, L =12λ1+12λ2.

Similarly, the coefficients for temporal derivative can be also determined. Thus, the gas distribution function Eq.(3.2) and the numerical fluxes Eq.(12) can be fully constructed.

In order to utilize the two-stage temporal discretization, the temporal derivatives of the flux function need to be determined. The flux in the time interval is expanded as the following linear form

 Fi+1/2,j(t)=Fni+1/2,j+∂tFni+1/2,j(t−tn). (17)

where the coefficients and can be fully determined by solving the linear system

 Fni+1/2,jΔt +12∂tFni+1/2,jΔt2=∫tn+ΔttnFi+1/2,j(t)dt, 12Fni+1/2,jΔt +18∂tFni+1/2,jΔt2=∫tn+Δt/2tnFi+1/2,j(t)dt.

Similarly, the coefficients and corresponding to the flux in -direction can be constructed as well. According to Eq.(5), and its temporal derivative at can be given by

 L(Qnij) =−1Δx(Fni+1/2,j(t)−Fni−1/2,j(t))−1Δy(Gni,j+1/2(t)−Gni,j−1/2(t)), ∂tL(Qnij)= −1Δx(∂tFni+1/2,j(t)−∂tFni−1/2,j(t))−1Δy(∂tGni,j+1/2(t)−∂tGni,j−1/2(t)).

With the procedure at the intermediate state, can be constructed as well.

### 3.3 Discretization for radiative part

According to the definition of source terms, we only need to solve energy equations at two-stages for Eq.(6), which are nonlinear systems with respect to the temperature . At each stage, the energy equations can be simplified as the nonlinear system

 F(T)=0, (18)

where is a nonlinear function from to and is the number of cells. To obtain an approximate solution of the nonlinear system, the above system is rewritten as

 Jδ=−F,

where is the Jacobian of . In order to avoid calculating the Jabobian directly, can be approximated by the -derivative as follows

 J(T)⋅δ=F(T+σδ)−F(T)σ, (19)

where is a small scalar. In fact, it is not trivial to choose , which has great influence on the stability of algorithm and improper values will cause blowing up. According to GMRES-2 , a special is given as follows

 σ=√η∥δ∥22max{|T⋅δ|, typT⋅|δ|}⋅sign(T⋅δ),

where is the machine epsilon, , , is the typical size of and the typical size of real number equals to its order of magnitude plus one. With the initial guess for Eq.(18) and for Eq.(19), the initial residual is

 r(0)=−F−Jδ(0),

and its generated Krylov subspace is

 Km≡span{r(0),Jr(0),…,Jm−1r(0)}.

According to the nonlinear Newton-GMRES method GMRES-2 , an approximate solution can be given. Besides, in order to improve computation efficiency, a restart algorithm is applied in orthogonalization process GMRES-1 . In actual calculation, the restart times can be no more than 5 times if the restart step and convergence condition are set appropriately.

## 4 Numerical tests

In this section, the one and two dimensional tests cases are provided to the correctness and robustness of current scheme. In the computation, the collision time takes

 τ=ϵΔt+C|pl−prpl+pr|Δt,

where and denote the pressure on the left and right sides of the cell interface. In smooth flow regions, the collision term takes and the gas distribution function Eq.(3.2) will reduce to

 f=g0(1+¯At).

To achieve the spatial accuracy, the classical multidimensional fifth-order WENO reconstruction is adopted and more details can be found in GKS-high-1 ; WENO-Z . In order to eliminate the spurious oscillation and improve the stability, the WENO reconstruction is performed for the characteristic variables, and the detailed analysis for characteristics can be found in Radiation-Diffusion-4 . Without special statement, the gas with is used in the following numerical examples.

### 4.1 Accuracy test

The advection of density perturbation problem for Euler equations is extended to test the order of accuracy for RHE. In this test case, the initial condition is set as follows

 ρ(x)=1+0.2sin(πx), u(x)=1, p(x)=1,

where the computational domain is , the periodic boundary condition is applied and .

As reference, the case with and is tested. With these parameters, the RHE degenerates into the Euler equation and the exact solution can be given by

 ρ(x,t)=1+0.2sin(π(x−t)), u(x,t)=1, p