Interface fluxes (Rusanov)
for i = 1:nx+1
f[i] = 0.5*(fR[i]+fL[i]) - 0.5*ps[i]*(uR[i]-uL[i])
end
end
function fluxes(nx,u,f)
for i = 1:nx+1
f[i] = 0.5*u[i]*u[i]
end
end
✝ ✆
```
**Figure 11.** Evolution of shock from sine wave (*u*<sup>0</sup> = sin(2*πx*)) computed using conservative form of inviscid Burgers equation. The reconstruction of the field is done using WENO-5 scheme, and flux splitting and Riemann solver approach is used to compute the nonlinear term.

#### **5. One-Dimensional Euler Solver**

In this section, we present a solution to one-dimensional Euler equations. Euler equations contain a set of nonlinear hyperbolic conservation laws that govern the dynamics of compressible fluids neglecting the effect of body forces, and viscous stress [30]. The one-dimensional Euler equations in its conservative form can be written as

$$\frac{\partial q}{\partial t} + \frac{\partial F}{\partial \mathbf{x}} = 0,\tag{61}$$

where

$$q = \begin{pmatrix} \rho \\ \rho u \\ \rho \varepsilon \end{pmatrix}, \quad F = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho u h \end{pmatrix}.$$

in which

$$h = \varepsilon + \frac{p}{\rho}, \quad p = \rho(\gamma - 1)\left(\varepsilon - \frac{1}{2}\mu^2\right). \tag{62}$$

Here, *ρ* and *p* are the density, and pressure respectively; *u* is the horizontal component of velocity; *e* and *h* stand for the internal energy and static enthalpy, and *γ* is the ratio of specific heats. Equation (61) can also be written as

$$A\frac{\partial q}{\partial t} + A\frac{\partial q}{\partial \mathbf{x}} = 0,\tag{63}$$

where *A* = *<sup>∂</sup><sup>F</sup> <sup>∂</sup><sup>q</sup>* is the convective flux Jacobian matrix. The matrix *A* is given as

$$A = \frac{\partial F}{\partial q} = \begin{pmatrix} 0 & 1 & 0 \\ \phi^2 - u^2 & (3 - \gamma)u & \gamma - 1 \\ (\phi^2 - h)u & h - (\gamma - 1)u^2 & \gamma u \end{pmatrix},\tag{64}$$

where *φ*<sup>2</sup> = <sup>1</sup> <sup>2</sup> (*<sup>γ</sup>* <sup>−</sup> <sup>1</sup>)*u*2. We can write the matrix *<sup>A</sup>* as *<sup>A</sup>* <sup>=</sup> *<sup>R</sup>*Λ*<sup>L</sup>* in which <sup>Λ</sup> is the diagonal matrix of real eigenvalues of *A*, i.e., Λ = diag[*u*, *u* + *a*, *u* − *a*]. Here, *L* and *R* are the matrices of which the columns are the left and right eigenvectors of *A*; *L* = *R*−1. The matrices *L* and *R* are given below [31]

$$L = \begin{pmatrix} 1 - \frac{\phi^2}{a^2} & (\gamma - 1)\frac{u}{a^2} & -\frac{(\gamma - 1)}{a^2} \\ \phi^2 - ua & a - (\gamma - 1)u & \gamma - 1 \\ \phi^2 + ua & -a - (\gamma - 1)u & \gamma - 1 \end{pmatrix}, \quad R = \begin{pmatrix} 1 & \beta & \beta \\ u & \beta(u + a) & \beta(u - a) \\ \frac{\phi^2}{(\gamma - 1)} & \beta(h + ua) & \beta(h - ua) \end{pmatrix}, \tag{65}$$

where *a* is the speed of sound and is given by *a*<sup>2</sup> = *γp*/*ρ*, and *β* = 1/(2*a*2).

We use the grid similar to the one used for the conservative form of the inviscid Burgers equation and is shown in Figure 10. The semi-discrete form of the Euler equations can be written as

$$\frac{\partial q\_i}{\partial t} + \frac{F\_{i+1/2} - F\_{i-1/2}}{\Delta x} = 0,\tag{66}$$

where *qi* is the cell-centered values stored at nodal points and *Fi*±1/2 are the fluxes at left and right cell interface. We use third-order Runge-Kutta numerical method for time integration. We use WENO-5 reconstruction to compute the left and sight states of the fluxes at the interface. We include three different Riemann solvers to compute the flux *Fi*±1/2 at interfaces to be used in Equation (66).

#### *5.1. Roe's Riemann Solver*

First, we present the Roe solver for one-dimensional Euler equations. According to the Gudanov theorem [32], for a hyperbolic system of equations, if the Jacobian matrix of the flux vector is constant (i.e., *A* = *<sup>∂</sup><sup>F</sup> <sup>∂</sup><sup>q</sup>* = constant), the exact values of fluxes at the interfaces can be calculated as

$$F\_{i+1/2} = \frac{1}{2} \left( F\_{i+1/2}^R + F\_{i+1/2}^L \right) - \frac{1}{2} R |\Lambda| L \left( q\_{i+1/2}^R - q\_{i+1/2}^L \right),\tag{67}$$

where |Λ| is the diagonal matrix consisting of absolute values eigenvalues of the Jacobian matrix and the matrices *L* and *R* are given in Equation (65). However, Jacobian matrix is not constant (i.e., *F* = *F*(*q*)). The Roe solver [33] is an approximate Riemann solver and it uses below formulation to find the numerical fluxes at the interface

$$F\_{i+1/2} = \frac{1}{2} \left( F\_{i+1/2}^R + F\_{i+1/2}^L \right) - \frac{1}{2} \mathcal{R} |\Lambda| L \left( q\_{i+1/2}^R - q\_{i+1/2}^L \right),\tag{68}$$

where the bar represents the Roe average between the left and right states. In order to derive *A*¯ we need to know *u*¯, ¯ *h*, and *a*¯ and they are computed using Roe averaging procedure. In order to derive *u*¯, ¯ *h*, and *a*¯ Roe used Taylor series expansion of *F* around *q<sup>L</sup>* and *q<sup>R</sup>* points.

$$F(q) = F(q^L) + A(q^L)(q - q^L),\tag{69}$$

$$F(q) = F(q^R) + A(q^R)(q - q^R). \tag{70}$$

Neglecting higher-order terms and equation above two equations

$$F(q^L) - F(q^R) = A(q^L - q^R) + \left(A(q^R) - A(q^L)\right)q. \tag{71}$$

Roe [20] derived approximate matrix *<sup>A</sup>*¯ such that it satisfies *<sup>A</sup>*¯(*qL*, *<sup>q</sup>R*) <sup>→</sup> *<sup>A</sup>*(*q*) as *<sup>q</sup>*¯ <sup>→</sup> *<sup>q</sup>*. Therefore, we get the following set of equations (written in vector form)

$$F(q^L) - F(q^R) = \bar{A}(q^L - q^R). \tag{72}$$

Using Equation (72) we can compute *u*¯, ¯ *h*, and *a*¯. First we reconstruct the left and right states of *q* at the interface (i.e., *q<sup>R</sup> <sup>i</sup>*+1/2 and *<sup>q</sup><sup>L</sup> <sup>i</sup>*+1/2) using WENO-5 scheme. Then, we can compute the left and right state of the fluxes (i.e., *F<sup>L</sup>* and *FR*) using Equation (61). The Roe averaging formulas to compute approximate values for constructing *R*¯ and *L*¯ are given below

$$\vec{u} = \frac{\mu\_R \sqrt{\rho\_R} + \mu\_L \sqrt{\rho\_L}}{\sqrt{\rho\_R} + \sqrt{\rho\_L}},\tag{73}$$

$$\mathcal{H} = \frac{h\_{\mathcal{R}}\sqrt{\rho\_{\mathcal{R}}} + h\_{L}\sqrt{\rho\_{L}}}{\sqrt{\rho\_{\mathcal{R}}} + \sqrt{\rho\_{L}}},\tag{74}$$

where the left and right states are computed using WENO-5 reconstruction. The speed of the sound is computed using the below equation

$$d = \sqrt{(\gamma - 1)\left[\hbar - \frac{1}{2}\vec{n}^2\right]}.\tag{75}$$

The eigenvalues of the Jacobian matrix are *λ*<sup>1</sup> = *u*¯, *λ*<sup>2</sup> = *u*¯ + *a*, and *λ*<sup>3</sup> = *u*¯ − *a*. The implementation of Roe's Riemann solver is given in Listing 14.

**Listing 14.** Implementation of Julia function to calculate the right hand side of Euler equations using Roe's Riemann solver. ✞ ☎

```