un: matrix for storing solution at every time step
s, e = 1, nx
for k = 2:nt+1
i=1 # left boundary
a[i], b[i], c[i], d[i] = 0.0, 1.0, 0.0, 0.0
for i = 2:nx # interior points
a[i] = 12.0/(dx*dx)-2.0/(alpha*dt)
b[i] = -24.0/(dx*dx) - 20.0/(alpha*dt)
c[i] = 12.0/(dx*dx)-2.0/(alpha*dt)
d[i] = -2.0/(alpha*dt)*(un[k-1,i+1] + 10.0*un[k-1,i] + un[k-1,i-1])
end
```
*Fluids* **2019**, *4*, 159

```
i = nx+1 # right boundary
a[i], b[i], c[i], d[i] = 0.0, 1.0, 0.0, 0.0
un[k,:] = tdma(a,b,c,d,s,e) # thomas algorithm function
end
```
✝ ✆

**Figure 5.** Comparison of exact solution and numerical solution computed using implicit compact Pade numerical scheme. The numerical solution is computed using Δ*t* = 0.0025 and Δ*x* = 0.025.

#### **3. Inviscid Burgers Equation: Non-Conservative Form**

In this section, we discuss the solution of the inviscid Burgers equation which is a nonlinear hyperbolic partial differential equation. The hyperbolic equations admit discontinuities, and the numerical schemes used for solving hyperbolic PDEs need to be higher-order accurate for smooth solutions, and non-oscillatory for discontinuous solutions [19]. The inviscid Burgers equation is given below

$$
\mu \frac{\partial u}{\partial t} + \mu \frac{\partial u}{\partial x} = 0.\tag{29}
$$

First we use FTCS numerical scheme for discretization of the inviscid Burgers equation as follow

$$\frac{u\_i^{(n+1)} - u\_i^{(n)}}{\Delta t} + u\_i^{(n)} \frac{u\_{i+1}^{(n)} - u\_{i-1}^{(n)}}{2\Delta x} = 0. \tag{30}$$

We use the shock generated by a sine wave to demonstrate the capability of FTCS numerical scheme to compute the numerical solution. The initial condition is *u*<sup>0</sup> = sin(2*πx*). We integrate the solution from time *t* = 0 to *t* = 0.25 with Δ*t* = 0.0001 and divide the computational domain into 200 grids. The numerical solution computed by FTCS at 10 equally spaced time steps between initial and final time is shown in Figure 6. It can be observed that the FTCS scheme does not capture the shock formed at *x* = 0.5 and produce oscillations in the solution. It can be demonstrated that FTCS scheme is numerically unstable for advection equations.

There are several ways in which the discontinuity in the solution (e.g., shock) can be handled. There are classical shock-capturing methods like MacCormack method, Lax-Wanderoff method, and Beam-Warming method as well as flux limiting methods [20,21]. We can also use higher-order central difference schemes with artificial dissipation [22]. Low pass filtering might also help to diminish the growth of the instability modes [23]. There is also a class of modern shock capturing methods called as high-resolution schemes. These high-resolution schemes are generally upwind biased and take the direction of the flow into account. Here, in particular, we discuss higher-order weighted essentially non-oscillatory (WENO) schemes [24] that are used widely as shock-capturing numerical methods. We also use compact reconstruction weighted essentially non-oscillatory (CRWENO) schemes [25] which have lower truncation error compared to WENO schemes.

**Figure 6.** Evolution of shock with time for the initial condition *u*<sup>0</sup> = sin(2*πx*) using FTCS scheme.

#### *3.1. WENO-5 Scheme*

We refer readers to a text by Shu [26] for a comprehensive overview of the development of WENO schemes, their mathematical formulation, and the application of WENO schemes for convection dominated problems. We will use WENO-5 scheme for one-dimensional inviscid Burgers equation with the shock developed by a sine wave. We will discuss how one can go about applying WENO reconstruction for finite difference schemes. The finite difference approximation to the inviscid Burgers equation in non-conserved form can be written as

$$\frac{\partial u\_i}{\partial t} + u\_i \frac{u\_{i+\frac{1}{2}}^L - u\_{i-\frac{1}{2}}^L}{\Delta x}, \quad \text{if} \quad u\_i > 0 \tag{31}$$

$$\frac{\partial u\_i}{\partial t} + u\_i \frac{u\_{i+\frac{1}{2}}^R - u\_{i-\frac{1}{2}}^R}{\Delta x}, \quad \text{otherwise}, \tag{32}$$

which can be combined into a more compact upwind/downwind notation as follows

$$\frac{\partial u\_i}{\partial t} + u\_i^+ \frac{u\_{i+\frac{1}{2}}^L - u\_{i-\frac{1}{2}}^L}{\Delta x} + u\_i^- \frac{u\_{i+\frac{1}{2}}^R - u\_{i-\frac{1}{2}}^R}{\Delta x} = 0,\tag{33}$$

where we define *u*<sup>+</sup> *<sup>i</sup>* = max(*ui*, 0) and *u*<sup>−</sup> *<sup>i</sup>* = min(*ui*, 0). the reconstruction of *<sup>u</sup><sup>L</sup> i*+ <sup>1</sup> 2 uses a biased stencil with one more point to the left, and that for *u<sup>R</sup> i*+ <sup>1</sup> 2 uses a biased stencil with one more point to the right. The stencil used for reconstruction of left and right side state (i.e., *u<sup>L</sup> <sup>i</sup>*+1/2 and *<sup>u</sup><sup>R</sup> <sup>i</sup>*+1/2) is shown in Figure 7. The WENO-5 reconstruction for left and right side reconstructed values is given as [24]

$$u\_{i+\frac{1}{2}}^L = w\_0^L \left(\frac{1}{3}u\_{i-2} - \frac{7}{6}u\_{i-1} + \frac{11}{6}u\_i\right) + w\_1^L \left(-\frac{1}{6}u\_{i-1} + \frac{5}{6}u\_i + \frac{1}{3}u\_{i+1}\right) + w\_2^L \left(\frac{1}{3}u\_i + \frac{5}{6}u\_{i+1} - \frac{1}{6}u\_{i+2}\right), \tag{34}$$

$$u\_{l-\frac{1}{2}}^R = w\_0^R \left( -\frac{1}{6}u\_{l-2} + \frac{5}{6}u\_{l-1} + \frac{1}{3}u\_l \right) + w\_1^R \left( \frac{1}{3}u\_{l-1} + \frac{5}{6}u\_l - \frac{1}{6}u\_{l+1} \right) + w\_2^R \left( \frac{11}{6}u\_l - \frac{7}{6}u\_{l+1} + \frac{1}{3}u\_{l+2} \right), \tag{35}$$

where nonlinear weights are defined as

$$w\_k^L = \frac{n\_k}{n\_0 + n\_1 + n\_2}, \quad u\_k = \frac{d\_k^L}{(\beta\_k + \epsilon)^2}, \quad k = 0, 1, 2\tag{36}$$

*Fluids* **2019**, *4*, 159

$$w\_k^R = \frac{a\_k}{a\_0 + a\_1 + a\_2}, \quad a\_k = \frac{d\_k^R}{(\beta\_k + \varepsilon)^2}, \quad k = 0, 1, 2\tag{37}$$

in which the smoothness indicators are defined as

$$\beta\_0 = \frac{13}{12}(u\_{i-2} - 2u\_{i-1} + u\_i)^2 + \frac{1}{4}(u\_{i-2} - 4u\_{i-1} + 3u\_i)^2,\tag{38}$$

$$\beta\_1 = \frac{13}{12}(u\_{i-1} - 2u\_i + u\_{i+1})^2 + \frac{1}{4}(u\_{i-1} - u\_{i+1})^2,\tag{39}$$

$$\beta\_2 = \frac{13}{12}(u\_i - 2u\_{i+1} + u\_{i+2})^2 + \frac{1}{4}(3u\_i - 4u\_{i+1} + 3u\_{i+2})^2. \tag{40}$$

We use *d<sup>L</sup>* <sup>0</sup> = 1/10, *<sup>d</sup><sup>L</sup>* <sup>1</sup> = 3/5, and *<sup>d</sup><sup>L</sup>* <sup>2</sup> = 3/10 in Equation (36) to compute nonlinear weights for calculation of left side state using Equation (34) and *d<sup>R</sup>* <sup>0</sup> = 3/10, *<sup>d</sup><sup>R</sup>* <sup>1</sup> = 3/5, and *<sup>d</sup><sup>R</sup>* <sup>2</sup> = 1/10 in Equation (37) to compute nonlinear weights for calculation of right side state using Equation (35). We set <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>6</sup> to avoid division by zero. To avoid repetition, we provide only the implementation of Julia function to compute the left side state in Listing 6. The right side state is computed in a similar way with coefficients for *d<sup>R</sup> <sup>k</sup>* in Equation (37).

We perform the time integration using third-order Runge-Kutta numerical scheme as discussed in Section 2.2 and its implementation is similar to Listing 2. The right hand side term of the inviscid Burgers equation is computed using the point value *ui* and the derivative term is evaluated using the left and right side states at the interface. Based on the direction of flow i.e., the sign of *ui*, we either use left side state or right side state to compute the derivative. The implementation of right hand side state calculation in Julia for the inviscid Burgers equation is given in Listing 7.

**Figure 7.** Finite difference grid for one-dimensional problems where the the solution is stored at the interface of the cells. The first and last points lie on the left and right boundary of the domain. The stencil required for reconstruction of left and right side state is highlighted with blue rectangle. Ghost points are shown by red color.

**Listing 6.** Implementation of left side state calculation function in Julia. ✞ ☎

```