f: flux at the interface computed using Rusanov scheme
function roe(nx,gamma,uL,uR,f,fL,fR)
dd = Array{Float64}(undef,3)
dF = Array{Float64}(undef,3)
V = Array{Float64}(undef,3)
gm = gamma-1.0
for i = 1:nx+1
#Left and right states:
rhLL = uL[i,1]
uuLL = uL[i,2]/rhLL
eeLL = uL[i,3]/rhLL
ppLL = gm*(eeLL*rhLL - 0.5*rhLL*(uuLL*uuLL))
hhLL = eeLL + ppLL/rhLL
rhRR = uR[i,1]
```
*Fluids* **2019**, *4*, 159

```
uuRR = uR[i,2]/rhRR
eeRR = uR[i,3]/rhRR
ppRR = gm*(eeRR*rhRR - 0.5*rhRR*(uuRR*uuRR))
hhRR = eeRR + ppRR/rhRR
alpha = 1.0/(sqrt(abs(rhLL)) + sqrt(abs(rhRR)))
uu = (sqrt(abs(rhLL))*uuLL + sqrt(abs(rhRR))*uuRR)*alpha
hh = (sqrt(abs(rhLL))*hhLL + sqrt(abs(rhRR))*hhRR)*alpha
aa = sqrt(abs(gm*(hh-0.5*uu*uu)))
D11 = abs(uu)
D22 = abs(uu + aa)
D33 = abs(uu - aa)
beta = 0.5/(aa*aa)
phi2 = 0.5*gm*uu*uu
R11, R21, R31 = 1.0, uu, phi2/gm
R12, R22, R32 = beta, beta*(uu + aa), beta*(hh + uu*aa)
R13, R23, R33 = beta, beta*(uu - aa), beta*(hh - uu*aa)
L11, L21, L31 = 1.0-phi2/(aa*aa), phi2 - uu*aa, phi2 + uu*aa
L12, L22, L32 = gm*uu/(aa*aa), aa - gm*uu, -aa - gm*uu
L13, L23, L33 = -gm/(aa*aa), gm, gm
for m = 1:3
V[m] = 0.5*(uR[i,m]-uL[i,m])
end
dd[1] = D11*(L11*V[1] + L12*V[2] + L13*V[3])
dd[2] = D22*(L21*V[1] + L22*V[2] + L23*V[3])
dd[3] = D33*(L31*V[1] + L32*V[2] + L33*V[3])
dF[1] = R11*dd[1] + R12*dd[2] + R13*dd[3]
dF[2] = R21*dd[1] + R22*dd[2] + R23*dd[3]
dF[3] = R31*dd[1] + R32*dd[2] + R33*dd[3]
for m = 1:3
f[i,m] = 0.5*(fR[i,m]+fL[i,m]) - dF[m]
end
end
end
```
We test the Riemann solver using Sod shock tube problem [34]. The time evolution of Sod shock tube problem is governed buy Euler equations. This problem is used for testing numerical schemes to study how well they can capture and resolve shocks and discontinuities and their ability to produce correct density profile. The initial conditions for the Sod shock tube problem are

✝ ✆

$$
\begin{pmatrix} \rho\_L \\ p\_L \\ u\_L \end{pmatrix} = \begin{pmatrix} 1.0 \\ 1.0 \\ 0.0 \end{pmatrix} \quad \text{when } \mathbf{x} < 0.5; \qquad \begin{pmatrix} \rho\_R \\ p\_R \\ u\_R \end{pmatrix} = \begin{pmatrix} 0.125 \\ 0.1 \\ 0.0 \end{pmatrix} \quad \text{when } \mathbf{x} > 0.5. \tag{76}
$$

We use third-order Runge-Kutta numerical scheme for time integration from *t* = 0 to *t* = 0.2. We use two grid resolutions to see the capability of Riemann solver to capture the shock. The Sod shock tube problem has Dirichlet boundary condition at the left and right boundary. We use linear interpolation to find the values of *q* at ghost points. Figure 12 shows the density, velocity, pressure, and energy at final time *t* = 0.2 using Roe solver. We consider the solution obtained with high-grid resolution as the true solution and compare it with the low-resolution results. From Figure 12, we observe that there are small oscillations near discontinuities at low-resolution. We use WENO reconstruction to compute the flux at interface and WENO scheme is not strictly non-oscillatory. However, WENO scheme does not allow oscillations to amplify to large values. To remove these oscillations further, for example, we can use characteristic-wise reconstructions [35].

**Figure 12.** Evolution of density, velocity, energy, and pressure at *t* = 0.2 for the shock tube problem computed using Riemann solver based on Roe averaging. The true solution is calculated with *N* = 8192 grid resolution with Δ*t* = 0.00005 and the low-resolution results are for *N* = 256 grid resolution with Δ*t* = 0.0001.

#### *5.2. HLLC Riemann Solver*

We present one of the most widely used approximate Riemann solver based on HLLC scheme [30,36]. They used lower and upper bounds on the characteristics speeds in the solution of Riemann problem involving left and right states. These bounds are approximated as

$$S\_L = \min(\mu\_{L\prime}, \mu\_R) - \max(a\_{L\prime}, a\_R), \tag{77}$$

$$S\_R = \min(\mu\_{L\prime}, \mu\_R) + \max(a\_{L\prime}, a\_R)\_\prime \tag{78}$$

where *SL* and *SR* are the lower and upper bound for the left and right state characteristics speed. For HLLC scheme we also include middle wave of speed *S*<sup>∗</sup> given by

$$S\_{\*} = \frac{p\_R - p\_L + \rho\_L u\_L (S\_L - u\_L) \rho\_R u\_R (S\_R - uR)}{\rho\_L (S\_L - u\_L) - \rho\_R (S\_R - u\_R)}. \tag{79}$$

The mean pressure is given by

$$P\_{LR} = \frac{1}{2} \left[ p\_L + p\_R + \rho\_L (S\_L - \mu\_L)(S\_\ast - \mu\_l) + \rho\_R (S\_R - \mu\_R)(S\_\ast - \mu\_R) \right]. \tag{80}$$

The fluxes are computed as

$$F\_{i+1/2} = \begin{cases} F^L, & \text{if } S\_L \ge 0\\ F^R, & \text{if } S\_R \le 0\\ \frac{S\_\*(S\_L u\_L F^L) + S\_L P\_{LR} D\_\*}{S\_L - S\_\*}, & \text{if } S\_L \le 0 \text{ and } S\_\* \ge 0\\ \frac{S\_\*(S\_R u\_R F^R) + S\_R P\_{LR} D\_\*}{S\_R - S\_\*}, & \text{if } S\_R \ge 0 \text{ and } S\_\* \le 0 \end{cases} \tag{81}$$

where *<sup>D</sup>*<sup>∗</sup> = [0, 1, *<sup>S</sup>*∗]. Once the flux is available at the interface, we can integrate the solution in time using Equation (66). The implementation of HLLC scheme is given in Listing 16. We perform the time integration using third-order Runge-Kutta numerical scheme. Figure 13 shows the density, velocity, pressure, and energy at final time *t* = 0.2 computed using Riemann solver based on HLLC scheme. The oscillations in the numerical solution calculated by the HLLC scheme is slightly less than those calculated by the Roe's Riemann solver. The true solution is computed using *N* = 8192 grid points and Δ*t* = 0.00005 and is compared with the low-resolution results for *N* = 256 grid points and Δ*t* = 0.0001.

**Listing 16.** Implementation of Julia function to calculate the interface flux using Riemann solver based on HLLC scheme. ✞ ☎

```