Discontinuity Preserving Scheme

Non-linear schemes are widely used in high-speed flows to capture the discontinuities present in the solution. Limiters and weighted essentially non-oscillatory schemes (WENO) are the most common non-linear numerical schemes. Most of the high-resolution schemes use the piecewise parabolic reconstruction (PPR) technique for their robustness. However, it may be impossible to achieve non-oscillatory and dissipation-free solutions with the PPR technique without non-linear switches. Most of the shock-capturing schemes use excessive dissipation to suppress the oscillations present in the discontinuities. To eliminate these issues, an algorithm is proposed that uses the shock-capturing scheme (SCS) in the first step, and then the result is refined using a novel scheme called the Discontinuity Preserving Scheme (DPS). The present scheme is a hybrid shock capture-fitting scheme. The present scheme has outperformed other schemes considered in this work, in terms of shock resolution in linear and non-linear test cases. The most significant advantage of the present scheme is that it can resolve shocks with three grid points. KeywordsShock capturing scheme, WENO, High-resolution schemes, Conservative schemes, Finite volume method.


Introduction
In the finite volume framework, shock-capturing scheme (SCS) is often implemented using a flux difference splitting approach. The two essential components are suitable reconstruction operator to determine approximate left and right states, and an approximate Riemann solver to obtain the numerical flux. The numerical dissipation is present in both these steps. Unlike in the case of a smooth solution, obtaining a higher-order result using a linear reconstruction is not possible in the presence of discontinuities. The order barrier theorem of Godunov (1959) stated that all schemes higher than the first-order are non-monotonous. The use of first-order schemes would result in excessive dissipation and smearing of shocks. Van Leer (1979) extended the SCS to second-order for shock problems using the concept of Total Variation Diminishing (TVD). The TVD schemes make use of nonlinear switching functions. However, they reduce to lower-order at the shocks and suffer from loss of accuracy in non-smooth extrema.
Because TVD schemes introduce excessive dissipation and reduce to lower-order at shocks, they give diffused results. These limitations were partially overcome by Harten et al. (1987). By the introduction of essentially non-oscillatory (ENO) scheme. Under this scheme, the solution's smoothness is evaluated on several stencils, and then the flux is determined based on the smoothest stencil that eliminates interpolating through discontinuity. These schemes are nonlinear as the interpolation coefficients depend on the solution. This idea was extended by Liu et al. (1994) through the WENO scheme. Here, instead of using the smoothest stencil, the weighted average of the fluxes on each stencil is taken. The weights are based on the smoothness of the solution in each stencil. WENO scheme is less diffusive than the limiters. However, the FVM codes using WENO reconstruction schemes cannot accurately produce step-like solutions because of dissipation introduced by parabolic basis, directional and artificial dissipation introduced by Riemann solver.
Most of the SCS are based on PPR (Collella and Woodward, 1984). It is known that the PPR without a non-linear switch suffers from Runge's phenomenon, so Marquina (1994) investigated non-parabolic reconstruction. Deng et al. (2018) and Xiao et al. (2005) use a reconstruction technique that significantly reduces the numerical dissipation present in the reconstruction basis. Their scheme named THINC (tangent of hyperbola for interface) uses a hyperbolic tangent function, which fits a step-like discontinuity well. Also, this function is differentiable and monotone. However, this procedure cannot reduce the artificial dissipation caused by the Riemann solver. It is impossible to solve shock problems without Riemann solver or artificial dissipation using existing algorithms. Eliminating the Riemann solver or reduction of artificial dissipation may reduce the robustness of the solver. It is common practice to refine the grid or try to achieve a higher-order reconstruction to get an accurate shock structure. The former increases the cost of computation, and the latter leads to oscillations in the solution. To overcome this issue, we have proposed a scheme to remove the numerical dissipation present in the schemes.
The conservative equations are solved using SCS in this method. WENO (Shu, 2009) is being used in the present study for the Euler equation, and backward in space discretization is used for the convective equation. To preserve the discontinuities present in the solution, DPS is applied to the solution of SCS. This approach is similar to the JST-scheme introduced by Jameson (2017), where artificial dissipation or filtering is carried out on a solution obtained from the regular SCS to eliminate the oscillations present in the solution. We have validated our scheme for linear convection equation and non-linear Euler equations and found some significant improvement over others in terms of resolution. The present scheme can resolve the discontinuity with three grid points, irrespective of the SCS and Riemann solver used. This method is similar to the THINC scheme of Deng et al. (2018), but THINC is a shock-capturing scheme, and the present scheme is the hybrid shock capturing-fitting scheme.

Methodology
The details of the conservative equation and numerical schemes used are described in this section. We have used the WENO scheme (Jiang and Shu, 1996) for interpolation. As the exact Riemann solvers are computationally expensive, Rusanov's scheme (also referred to as Local Lax-Fredrich / LLF) is used as an approximate Riemann solver. The values at the cell-face are calculated by interpolating conservative variables.

Discretization of the Euler Equation
The Euler equation locally produces step-function like discontinuity and has the exact solution for some cases, we have tested our schemes on the Euler equation. One-dimensional Euler equation is: Using conservative discretization and method of lines it becomes Once the right-hand side of (2) is calculated, time-marching is performed using the HRK31 method (Neelan and Nair, 2018) for the Euler equation and the forward Euler method for the linear convection equation.

WENO Algorithm
The flux ( F ) is a function of the conservative variable U so the face values (i.e.) 0.5 i U  are calculated using the WENO scheme and its basis are (Biswas and Dubey, 2018)   0 0.5 2 1 1 2 7 11 6 The ideal weights of this formulation are: 0 1 10   , 1 6 10   and 2 3 10 The final reconstruction function is: where,  is a small number added to the denominator to avoid getting an unphysical number. j  is the smoothness indicator used in Appadu and Peer (2013).

Discontinuity Preserving Scheme
Once the conservative equation is solved using SCS, the solution is refined using DPS. The DPS first detects the location of the shock by calculating the gradient. The gradient is calculated using the following formulations: The position of the critical points (shocks) is obtained by calculating the first gradient of  using local maxima calculating procedure. Once the shock locations are identified, a hyperbolic-tangent function is used to fit the shock. The function is similar to the function proposed by Deng et al. (2018) that is Let 1 k be the location of the discontinuity. It takes a value 0 i  (center of the discontinuity) locally and t s is the span of smoothing function. In most of the cases, 5 t s  to 15 and give a good result. If two shocks are very close, Eq. (4) is applied to the solution of SCS over the interval 1 t ks  to 1 t ks  . This means DPS is only applied at the location of discontinuities.

Numerical Result and Discussion
To check the robustness of the present algorithm, we have tested the scheme on linear convection equation and the Euler equation.

Linear Convection Equation
One-dimensional linear convection equation is: The initial condition used is: The problem is solved in the domain 22 to   . Forward Euler method is used for time integration, and the backward Euler method is used for spatial discretization (FTBS). The number of grid points () n used is 101 and solved using 0.8 CFL  and simulated up to 1.0053s flow time. Figure 1 and Figure 2 demonstrate the convection equation solution. The current scheme is less dissipative than FTBS. It is worth noting that the current scheme (FTBS-DPS) used three points to resolve the discontinuity, but FTBS used twelve points to resolve it.

Sod Shock-Tube Problem
The present scheme is tested on Sod shock tube problem (Sod, 1978). The problem is solved using 400 grid points with CFL number 0.5, up to flow time = 0.1 s. The initial condition of the problem is: The Euler equation is solved using the finite volume method (FVM), interpolation is performed using WENO and Rusanov's scheme is used to calculate interface flux. The velocity plots of the problem are shown in Figure 3 and Figure 4. It is evident from Figure 4 that the current model has a better high-resolution property than the standard SCS, and it is comparable with the analytical solution. To resolve the discontinuity, the WENO scheme used eight grid points for this problem, THINC used five grid points, but the proposed scheme used three points.  (1, 3.55, 1) 0.5. The discretization procedure used here is the same as that for the Sod shock tube problem. Figure  5 and Figure 6 show the velocity plots. The present scheme shows a better shock resolving property than the WENO scheme. Incorrect shock strength is due to an incorrect weighting of the WENO scheme. It originated from the reconstruction of density in the left-state (i-2 to i+2), and it needs some improvement. In this case, with some dissipation, the limiters can capture the shock strength correctly.

Right Expansion and Left Strong Shock
The problem is solved using 400 grid points with CFL = 0.8 up to flow time 0.01 s. The Riemann solver used here is Rusanov's scheme. The initial condition used is:   (7, 0, 1) 0.5 ,, (10, 0, 1) 0.5 The solution to this problem using WENO, WENO-DPS and the exact solution is shown in Figure 7 and Figure 8. Here, WENO scheme used 13 grid points to resolve the discontinuity but WENO-DPS used 3 grid points.

Lax Test Case
Lax test case (Lax and Liu, 1998) is solved using a total number of The Lax test case solution can be found in Figure 9, Figure 10, Figure 11 and Figure 12. It is evident from Figure 10 that WENO-DPS cannot completely remove the oscillations in the solution but can be cured by carefully adjusting the span of WENO-DPS shown in Figure 12, where 20 t s  it is used. WENO used twelve points to resolve the discontinuity in the convection equation, but the current algorithm (WENO-DPS) used three points to achieve the same.
The density solution to this problem is shown in Figure 13 and Figure 14. For this problem, WENO-DPS did not show significant improvement in the solution over WENO and THINC.
In the Sod shock tube problem, WENO scheme used eight points to resolve the discontinuity and THINC used five grid points, but the proposed algorithm used three points. Similarly, for the Mach 3 test case, the current scheme was capable of resolving shocks with three points, but the WENO scheme used twelve points, and THINC used five grid points. So, the current scheme could be bundled with SCS to avoid numerical dissipation of shock. This algorithm results in a fifth-order convergence in smooth regions and a third-order convergence in a non-smooth region, similar to the WENO scheme. This can be checked by expanding the Taylor series.

Two-Dimensional DPS
Flow past a 10-degree wedge with an incoming Mach number of 6 is considered. The spatial terms are solved using the Roe-Riemann solver along with an edge-based limiting technique using minmod limiter. Time advancement is achieved using the implicit Euler method. Unstructured quadrilateral grid is used with a node-centred, second-order finite-volume method. The domain [ 2,3] [0, 5]  is discretized using 100 100  quadrilateral cells. A second-order Gaussian kernel, moving least square procedure is employed to obtain gradient values. CFL ramping is used to improve stability at the start of time advancement.
We could split the 2-dimensional problem using two 1-dimensional problems and apply the DPS algorithm, but it may lead to the stair-case effect. Figure 15 and Figure 16 show the solution of the above test case post-processed using 1-D DPS. In this case, DPS algorithm performed worse than the standard shock capturing algorithm. We need some modifications in the present 1-D DPS algorithm so an alternative 2-D DPS algorithm is presented. Here, the location of discontinuity is calculated using (3) in horizontal and vertical directions. Once the discontinuities are located, we fit the following function using moving least square algorithm. The function used is: Where, a, b, c and d are calculated using the Levenberg-Marquardt algorithm (Kanzow et al., 2005). The solution of DPS algorithm is shown in Figure 17. In post-processing, the geometry is removed from the control volume to make the interpolation simple. Figure 18 shows the pressure plot of SCS and SCS-DPS at y = 0.5 line. In this case, DPS able to resolve the shock with three grid points but SCS took 12 points to resolve the discontinuity.

Advantages and Disadvantages of the Present Scheme
It has an excellent shock resolving capacity, regardless of grid size/ CFL/ Riemann solver etc. It does not allow any numerical dissipation or dispersion at the shocks. It is computationally more economical because it can retain the shock structure even on a coarse grid. Some of the Riemann solvers may not have enough artificial viscosity to suppress the oscillations that can be cured using the present approach. In addition to that, a well-tuned DPS can produce solution matches with the analytical solution. The limitations of the method are cut-off should be tuned based on the problem, discretization and grid size. It is not recommended for very diffusive schemes, especially for the first-order schemes. If not tuned well, it may remove physical oscillations present in the vicinity of the shocks. It may not work well on multi-dimensional unstructured grids. If the grid is not aligned with the shock, it may lead to a stair-case effect. The Levenberg-Marquardt algorithm can also be used in the one-dimensional test cases, which could improve the results when the oscillations present in the vicinity of shocks. However, it is computationally expensive than the presented algorithm for the one-dimensional case.

Conclusions
A discontinuity preserving scheme (DPS) is introduced in the present study. The objective is to preserve the discontinuities present in the solution that could otherwise be dissipated by the numerical dissipation of shock capture schemes (SCS). This algorithm is similar to the artificial dissipation method and the spatial filters used in large eddy simulation to remove certain or unwanted oscillations present in the solution. Instead of applying DPS at each iteration, it can be applied at a few time iterations or at the end of time iteration to reduce the computational cost. For linear and non-linear test cases, the current scheme shows better shock-resolving properties than SCS. It can resolve the shock with a lesser number of grids than standard SCS. This algorithm offers a cost-effective way to solve the shocks with fewer grid points than the regular SCS. This scheme is a hybrid shock fitting shock capturing algorithm. The standard shock fitting algorithm relies on analytical expression, but the present scheme relies on numerical solution obtained using SCS. It is much simpler than the standard shock-fitting technique.