s,e: stat and end index of the matrix
function ctdms(a,b,c,alpha,beta,r,x,s,e)
bb = Array{Float64}(undef, e)
u = Array{Float64}(undef, e)
z = Array{Float64}(undef, e)
gamma = -b[s]
bb[s] = b[s] -gamma
bb[e] = b[e] - alpha*beta/gamma
for i = s+1:e-1
bb[i] = b[i]
end
tdms(a,bb,c,r,x,s,e) # call regular Thomas~algorithm
u[s] = gamma
u[e] = alpha
for i = s+1:e-1
u[i] = 0.0
end
tdms(a,bb,c,u,z,s,e) # call regular Thomas~algorithm
```

```
fact = (x[s] + beta*x[e]/gamma)/(1.0 + z[s] + beta*z[e]/gamma)
for i = s:e
x[i] = x[i] - fact*z[i]
end
end
✝ ✆
```
We test the CRWENO-5 method for the same problem as the WENO-5 scheme. The evolution of shock generated from sine wave with time is shown in Figure 9 for both Dirichlet and periodic boundary conditions computed using CRWENO-5 method. We use the same parameter as the WENO-5 scheme. We store snapshots of the solution field at different time steps to see the evolution of the shock. It can be seen from Figure 9 that the CRWENO-5 scheme captures the shock formed at final time *t* = 0.25 accurately.

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

#### **4. Inviscid Burgers Equation: Conservative Form**

In this section, we explain two approaches for solving the inviscid Burgers equation in its conservative form. The inviscid Burgers equation can be represented in the conservative form as below

$$\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0,\qquad\text{where}\quad f = \left(\frac{u^2}{2}\right). \tag{51}$$

The *f* is termed as the flux. We need to modify the grid arrangement to solve the inviscid Burgers equation in its conservative form. The modified grid for one-dimensional problem is shown in Figure 10. Here, <sup>Δ</sup>*<sup>x</sup>* <sup>=</sup> *xi*<sup>+</sup>1/2 <sup>−</sup> *xi*−1/2. The solution field is stored at the center of each cell and this type of grid is primarily used for finite volume discretization. We explain two approaches to compute the nonlinear term in the inviscid Burgers equation.

**Figure 10.** Finite difference grid for one-dimensional problems where the the solution is stored at the center of the cells. The flux is also computed at left and right boundary. The stencil required for reconstruction of left and right side flux is highlighted with blue rectangle. Ghost points are shown by red color.

#### *4.1. Flux Splitting*

The conservative form of the inviscid Burgers equation allows us to use the flux splitting method and WENO reconstruction to compute the flux at the interface. In this method, we first compute the flux *f* at the nodal points. The nonlinear advection term of the inviscid Burgers equation is computed as

$$\left.\frac{\partial f}{\partial x}\right|\_{x=x\_i} = \frac{f\_{i+1/2} - f\_{i-1/2}}{\Delta x},\tag{52}$$

$$\left. \frac{\partial f}{\partial x} \right|\_{x=x\_i} = \frac{f\_{i+1/2}^L - f\_{i-1/2}^L}{\Delta x} + \frac{f\_{i+1/2}^R - f\_{i-1/2}^R}{\Delta x} \tag{53}$$

at a particular node *i*. The superscripts *L* and *R* refer to the positive and negative flux component at the interface and the subscripts *i* − 1/2 and *i* + 1/2 stands for the left and right side interface of the nodal point *i*. We use WENO-5 reconstruction discussed in Section 3.1 to compute the left and right side flux at the interface. First, we compute the flux at nodal points and then split it into positive and negative components. This process is called as Lax-Friedrichs flux splitting and the positive and negative component of the flux at a nodal location is calculated as

$$f\_i^+ = \frac{1}{2}(f\_i + au\_i), \quad f\_i^- = \frac{1}{2}(f\_i - au\_i), \tag{54}$$

where *α* is the absolute value of the flux Jacobian (*α* = | *∂ f <sup>∂</sup><sup>u</sup>* |). We chose the values of *α* as the maximum value of *ui* over the stencil used for WENO-5 reconstruction, i.e.,

$$\mathfrak{a} = \mathfrak{m} \text{ax}(|\mathfrak{u}\_{i-2}|, |\mathfrak{u}\_{i-1}|, |\mathfrak{u}\_i|, |\mathfrak{u}\_{i+1}|, |\mathfrak{u}\_{i+2}|). \tag{55}$$

Once we split the nodal flux into its positive and negative component we reconstruct the left and right side fluxes at the interface using WENO-5 scheme using below formulas [24]

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

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

where the nonlinear weights are computed using Equations (36) and (37) . Once we have left and right side fluxes at the interface we use Equation (53) to compute the flux derivative. We demonstrate the flux splitting method for a shock generated by sine wave similar to Section 3. We use the third-order Runge-Kutta numerical scheme for time integration and periodic boundary condition for the domain. For the conservative form, we need three ghost points on left and right sides, since we compute the flux at the left and right boundary of the domain also. The periodic boundary condition for fluxes at ghost points is given by

$$f\_{-3} = f\_{\text{N}-2}, \quad f\_{-2} = f\_{\text{N}-1}, \quad f\_{-1} = f\_{\text{N}}, \quad f\_{\text{N}+1} = f\_1, \quad f\_{\text{N}+2} = f\_2, \quad f\_{\text{N}+3} = f\_3. \tag{58}$$

The Julia function to compute the right hand side of the inviscid Burgers equation is detailed in Listing 12. The implementation of flux reconstruction and weight computation in Julia is similar to Listings 6 and 8.

**Listing 12.** Implementation of Julia function to calculate the right hand side of inviscid Burgers equation using flux splitting method. ✞ ☎

```