Two Dimensional Analysis and Optimization of Hybrid MDCD-TENO Schemes

In this article, a quasi-linear semi-discrete analysis of shock capturing schemes in two dimensional wavenumber space is proposed. Using the dispersion relation of the two dimensional advection and linearized Euler equations, the spectral properties of a spatial scheme can be quantiﬁed in two dimensional wavenumber space. A hybrid scheme (HYB-MDCD-TENO6) which combines the merits of the Minimum Dispersion and Controllable Dissipation (MDCD) scheme with the Targeted Essentially Non-Oscillatory (TENO) scheme was developed and tested. Using the two dimensional analysis framework, the scheme was spectrally optimized in such a way that the linear part of the scheme can be separately optimized for its dispersion and dissipation properties. In order to compare its performance against existing schemes, the proposed scheme as well as the baseline schemes were tested against a series of benchmark test cases. It was found that the HYB-MDCD-TENO6 scheme provides similar or better resolution as compared to the baseline TENO6 schemes for the same grid size.

to their superior spectral properties, spectral methods and compact schemes are readily used. However, these methods are mainly limited to flows without discontinuities in the solution (such as shockwaves) and are found to cause non-physical oscillations when applied to flows with discontinuities. On the other hand, the use of Essentially Non Oscillatory (ENO) and Weighted Essentially Non Oscillatory (WENO) type schemes [14] are found to be robust for flows with discontinuities. However, these schemes in their original form are found to be too dissipative for turbulent flows. In addition, their dispersion properties may not be optimal. Research efforts documented in the literature have concentrated on the development of numerical schemes with high resolution and good shock capturing capabilities.
In order to develop and improve such schemes, it is important to have a framework to analyse and characterise a non-linear scheme. For this purpose, the Approximate Dispersion Relation (ADR) method was introduced in [17]. Subsequently, Li [12] introduced a simplified form of non-linear spectral analysis method based on the ADR method. The key idea of this method is to consider the semi-discretized one dimensional advection equation and to use only the fundamental modes obtained from a Fourier transform of the semi-discretized equation to deduce the modified wavenumber. A further observation concerning such an analysis is that most of the analysis technique are usually based on the semi-discretized dispersion relation of the one dimensional advection equation. However, when one considers multi dimensional problems, schemes that are analysed and optimized in one dimensional wavespace may not have isotropic behaviour. Therefore, the main advantage of a multi dimensional wavenumber analysis and optimization technique is the ability to quantify and reduce the amount of isotropic error. In this work, the method outlined in [12] is used and then extended to two dimensional wavenumber analysis.
Recent development of low dissipation shock capturing schemes can be found in [6], where Fu has developed a modified version of the WENO scheme (known as the TENO scheme) based on the use of a cut-off parameter which completely removes a sub schemes in its flux reconstruction when the normalized smoothness indicator exceeds a certain value. As compared to the WENO scheme, the TENO scheme is able to achieve better resolution. A natural choice to achieve better resolution properties is to combine the TENO/WENO/ENO scheme with another scheme which possesses spectral-like resolution to form a hybrid scheme. The development of such hybrid schemes can be found in [1,3,16,20]. Most of these work concerns the hybridization of the WENO scheme with either a compact, central or spectral scheme. One of the key issues with hybrid scheme concerns the smooth transition from one scheme to another as the flow field changes from a smooth to a region where the solution is discontinuous. Ren [18] improved the hybrid compact-WENO scheme of Pirozzoli [16] by designing a continuous weighting function to avoid the abrupt transition from the compact to WENO scheme. A key issue concerning hybrid scheme relates to the proper design of a hybridization strategy which ensures the robustness of the numerical solver in one dimensional and multi dimensional problems.
Other than the hybridization of schemes, another approach which can improve the spectral properties of a shock capturing schemes is to use optimization technique. Optimization technique for explicit and implicit finite difference schemes can be found in [7,11,19,21,22]. Most of these technique are usually based on the pioneering work of Tam and Webb who devised the dispersion-relationpreserving (DRP) schemes for computational acoustics. Within the context of a non-linear shock capturing schemes, the methods used in [21,26] are based on optimizing the weighting of each candidate stencils. In [7], the underlying linear sub scheme of the Targeted Essentially Non Oscillatory (TENO) scheme is optimized based on re-expressing the linear scheme into a separate set of finite difference coefficients that govern the dispersion and dissipation property of the numerical scheme separately. The dispersion properties are optimized by solving a set of linear algebraic equation based on the condition that the integrated error across the wavenumber range is minimized. No free parameter is involved. In this article, the method used in [21] is further explored in connection with the Targeted Essentially Non Oscillatory (TENO) scheme [6]. In [21], a class of finite difference scheme with minimized dissipation and controllable dispersion (MDCD) was developed. The key concept of the MDCD scheme is that the coefficients of the finite difference scheme can be expressed in terms of two free independent parameters which allow for separate optimization of its dispersion and dissipation properties. The concept of the MDCD scheme is then coupled with the TENO scheme to allow for the weighting function of each candidate stencil to be expressed in term of two free independent parameters. This leads to the development of the MDCD-TENO scheme. This approach differs from the approach used in [7], where the coefficients of the underlying linear sub scheme is being optimized instead. Since the nonlinear mechanism of the shock capturing scheme is needed only in the vicinity of a discontinuity and that it can drastically worsen the spectral properties of the scheme in smooth flow, a hybrid scheme which blends the linear MDCD scheme with the non-linear MDCD-TENO scheme is proposed. A robust blending method based on the inherent normalized smoothness indicator of the TENO scheme is proposed. The hybrid MDCD-TENO scheme is then spectrally optimized in two dimensional wavenumber space for its dispersion property.
Concerning the desired spectral properties of a finite difference schemes, it is generally accepted that the dispersion error should be minimized based on some chosen criteria such as the dispersion relation of the one dimensional advection equation. As for the dissipation of the numerical scheme, there exists no established guidelines for how the dissipation should be optimized. In practice, a scheme with too much dissipation makes it unsuitable for DNS. On the other hand, the use of non dissipative central difference scheme directly on the Navier Stokes equation can lead to numerical instabilities associated with the aliasing error arising from the non-linear advection terms [17]. As pointed out by Pirozzoli [17], a small amount of dissipation is useful since there exists significant dispersion error in the high wavenumber range and that it would be desirable to damp out the waves propagating at an incorrect speed. The main dilemma concerning the optimization of the dissipation is that the optimal dissipation is problem dependent. In this work, the developed approach allow the user to adjust the numerical dissipation through a free parameter. In the work of Fu [7], the choice of the leading finite difference coefficient that govern the dissipation property has a direct effect on numerical dissipation. However, the choice of the stencil length would impose additional constraint on its value, such that it must fall within specified range that satisfy the monotonicity condition to ensure that anti dissipation behaviour does not appear at intermediate wavenumber. A final note concerning the optimal dissipation can be found in [8], where the authors have proposed an approximate dispersion-dissipation relation for finite-difference schemes for estimating the dissipation required to damp spurious high-wavenumber waves based on a relation between the group velocity and the dissipation rate.
The paper is organized in the following manner. Firstly, a review of the existing WENO and TENO shock capturing schemes are presented. This is followed by a review of one dimensional non-linear spectral analysis method [12] and the Minimum Dispersion and Controllable Dissipation (MDCD) scheme [21]. Thereafter, the one dimensional non-linear spectral analysis method is extended to two dimensional wavenumber space. A hybrid scheme was then developed based on the combination of the Minimum Dispersion and Controllable Dissipation (MDCD) scheme and the Targeted Essentially Non Oscillatory (TENO) scheme. The two dimensional wavenumber analysis framework was then used to optimize the dispersion property of the hybrid scheme in two dimensional wavenumber space. Finally, the spectrally optimized hybrid schemes are compared against existing baseline schemes through a series of benchmark test cases.

Key concepts of WENO and TENO schemes
In order to describe the WENO scheme, a general hyperbolic conservation law in 1D is considered. The one dimensional conservation law can be described as follows: u(x, t) is the conservative variable and f (u) is the flux function. By considering a uniformly discretized space, equation (1) can be written in a semi-discrete form using a conservative finite difference scheme: In order to introduce the correct upwinding, the numerical flux at the cell interface is generally split into two parts as follow;f i+1/2 =f + i+1/2 +f − i+1/2 . In this paper, the flux vector splitting approach of Lax-Friedreich is used [9]. For the sake of brevity, only the positive part of the split numerical flux is described below. The stencils of the negative split numerical flux can be obtained through symmetry at the cell interface, x i+1/2 . The flux at the cell interface is interpolated with discrete f i of the surrounding nodal values using a convex combination of r candidatestencil fluxes. Mathematically, the numerical flux,f i+1/2 can be expressed as: The non-linear weights in equation (3) are defined as: represents a small number to avoid division by zero. For illustration purposes, the fifth order WENO scheme of Jiang and Shu [9] will be presented. r is equal to 3 for the WENO5-JS scheme. is set to 1.0e-6 and p is set to 2.0. The linear weights are given by c 0 = 0.1, c 1 = 0.6, c 2 = 0.3. The sub schemes to be used in equation (3) are given as: The smoothness indicator determines the weighting contribution of each sub schemes towards the flux reconstruction at the cell interface. The smoothness indicators, β i of the WENO5-JS scheme are given by: Details of the 7th order WENO-JS scheme can be found in [9,14].
Following this, the TENO5 and TENO6 scheme described in [4,5] will be introduced. The key differences between the TENO and WENO type schemes lies in the use of a cut off parameter to remove a linear sub scheme if it crosses a discontinuity. The other major difference between the WENO and TENO scheme lies in the sub stencils used in the flux reconstruction. While the sub schemes used within the WENO scheme has the same stencil width, the sub stencils within the TENO schemes are constructed with incremental width. For example, the 7th order WENO is compose of four 4 point stencils while the 7th order TENO is compose of three 3 point stencils and two 4 point stencils. In the case of having two closely located discontinuity at both end of the main scheme's stencil width, the TENO scheme is able to guarantee the ENO property since there exist one centrally located sub scheme which can avoid crossing with discontinuities located at either end of the main scheme stencil width. As for the WENO scheme, the sub stencils are located chronologically from the most upwind biased location to the most downwind biased location. This ensures that all WENO sub schemes crosses with the discontinuities located at either end of the main scheme stencil width. In this paper, only the TENO5 scheme is presented for completeness. Further details of the TENO6 scheme can be found in [5]. The numerical flux of the TENO scheme is expressed similarly as equation (3), with the exception of the weighting definition. The weighting of the TENO scheme is defined as: where δ k represents a sharp cut off function and is defined as: C T is taken to be 10 −6 and 10 −7 for the TENO5 and TENO6 schemes respectively in this paper. The normalised smoothness indicator, χ k is defined as: where β k is given as: For the TENO5 scheme, the sixth order τ 5 is given as: The linear sub schemes of the TENO5 scheme are given:

1D non-linear spectral analysis and MDCD scheme
In this section, the WENO and TENO schemes introduced earlier will be analysed using non-linear spectral analysis method. In addition, the concept of the Minimum Dispersion and Controllable Dispersion (MDCD) scheme will be introduced. Before presenting an analysis of these schemes, a short review of the non-linear spectral analysis method is presented.

1D non-linear spectral analysis
Unlike a linear scheme, the modified wavenumber of a non-linear scheme cannot be easily derived analytically and therefore a numerical approach was introduced to deduce the spectral properties of the numerical scheme. One of the first non-linear spectral analysis approach was the ADR method introduced in [17]. The essential idea of ADR is based on the ratio of the spectra solution at time, τ of the primary Fourier mode with respect to its initial complex amplitude of Fourier mode. For a single mode initial condition, the solution spectra at time τ obtained through the use of a non-linear shock capturing scheme contains multiple Fourier modes instead of a single mode as in the case of a linear scheme. Furthermore, the different Fourier modes may interact at later time stages which affect the evolution of each Fourier mode. The framework of the ADR method may be further simplified by neglecting the temporal discretization effect in the case of a very small time step and performing a Fourier analysis of the semi-discretized system instead [12]. A numerical flux, F j at node j can be written in terms of its spatial derivative as: In the context of a non-linear scheme, the flux function must contain the harmonic modes in addition to the fundamental mode, e ik 0 x j , that is whereF k 0,1,2... is the complex amplitude of the Fourier modes. This is in contrast to linear schemes, for which the flux function contains only the fundamental mode. Similar to the ADR approach, a quasi linear analysis accounts for only the leading order Fourier modes of the non-linear effects embodied within a shock capturing scheme. Therefore, only the primary mode is considered for further analysis. The complex amplitude of the leading order Fourier mode,F k 0 can thus be obtained through a Discrete Fourier Transform (DFT), F j is the flux obtained numerically through the non-linear differentiation given in equation (14). The real part of the modified wavenumber relates to the dispersion property and the imaginary part relates to the dissipation property. Using an initial condition, f j = e ik 0 x j , which is a sinusoidal function with wavenumber k 0   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 and only the leading order Fourier mode in equation (15), the numerical flux can be expressed as: Using equation (17), the modified wavenumber of the non-linear scheme can thus be expressed in terms of the leading order complex Fourier amplitude. The real and imaginary part of the modified wavenumber correponds to the dispersion and dissipation property of the numerical scheme.
The leading order complex Fourier mode is obtained through equation (16). In order to obtain the dispersion or dissipation curve along the entire wavenumber spectrum, the wavenumber of the initial condition is varied for the entire spectrum so as to produce a list of fundamental complex Fourier modes.

MDCD Scheme
The concept of the Minimum Dispersion and Controllable Dissipation (MDCD) scheme can be traced back to [21] where the authors have introduced a novel approach for one dimensional optimization of linear finite difference schemes based on two free parameter. It was shown that the dispersion and dissipation properties of the linear scheme can be separately controlled through two independent parameters; γ disp and γ diss . As compared to existing optimization technique, the main advantage of this technique lies in the modification of the dispersion and dissipation property separately without affecting one another. The general form of the MDCD schemes is given as [21]: where The binomial coefficient n k is defined by n k = n! k!(n−k)! . γ disp = α+β and γ disp = α − β are two independent parameters which can be used to separately optimize the dispersion and dissipation property of the scheme separately. It was shown mathematically that the dispersion and dissipation property of such a scheme can be expressed as [21] :   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 65 which proves that γ disp and γ diss are two independent parameter. In order to facilitate the construction of the MDCD scheme and the MDCD-TENO scheme which will be introduced in later sections, the spectral properties of the MDCD scheme will be studied in detail for the case of r = 3. When r = 3, the expanded form of the fourth order MDCD scheme is given as [21]: The real and imaginary part of the modified wavenumber of equation (21) can be written as:

Spectral properties of WENO, TENO and MDCD schemes
A comparison of the dispersion and dissipation property of the baseline WENO, TENO and MDCD schemes are presented in figure 1. For the purpose of comparison sake, various values of γ disp and γ diss are considered for the MDCD scheme. It can be seen that the TENO schemes have slightly better dispersion and lower dissipation as compared to WENO schemes having the same stencil width. For the fourth order MDCD scheme, it can be seen that a value of γ disp = 0.0545 gives the optimal dispersion property while a value of γ disp = 0.035 gives a less than optimal dispersion property. The dissipation of the MDCD scheme is negative as long as γ diss > 0. For a value of γ diss = 0, the MDCD scheme will have no numerical dissipation. With the use of increasing positive value of γ diss , the amount of numerical dissipation in the high wavenumber region increases.

Two dimensional optimization of non-linear schemes
In this part, we present a framework for which the non-linear scheme can be optimized in two dimensional wavenumber space based on the dispersion relation of the two dimensional advection equations. The basic idea is to optimize the weighting contribution of the linear part of the nonlinear scheme based on the dispersion relation of the two dimensional advection equation equations. In particular, the TENO scheme is split into a linear and nonlinear part from which the weighting 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 associated with the linear part of the TENO scheme can be expressed in terms of the two free parameters of the MDCD schemes. Expressing the linear scheme in terms of the weighting parameters allows for the specific tailoring of the dispersion and dissipation property of the numerical scheme. The only drawback of such an approach is that the schemes are optimized for their dissipation and dispersion properties at the expense of their formal order of accuracy. Different from [21], the cost function used in this optimization is the minimization of the integrated sum of dispersive error, E disp in two-dimensional wavenumber space. Since the dispersion and dissipation property of the numerical scheme can be separately optimized, the cost function only considers the dispersion error. The dispersion error for every reduced wavenumber combination is multiplied by an exponential term, e β(2π−k x ∆x−k y ∆y) so as to ensure a larger weighting for error closer to the lower wavenumber range. β is a parameter which controls the magnitude of the weighting.
(25) It should be noted that ∆ = ∆x = ∆y. r 1 and r 2 are set to 0.1 and 0.8 respectively. They represent the considered wavenumber range for which the error is integrated for the optimization process. The free parameters to be optimized in this context is the values of γ disp and γ diss used in the weighting of the sub scheme of the TENO6 scheme. The remaining parameters are kept the same as the original TENO6 scheme. By weighting the integrated error more in the low wavenumber range, the error in the low wavenumber range can be reduced more effectively through the optimization procedure.

The MDCD-TENO scheme
The semi-discrete finite difference scheme in conservative form can be written as: wheref i+1/2 is the numerical flux. The numerical flux of the TENO scheme can be rewritten in terms of a combination of a linear and nonlinear part: where   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 The linear part mainly affects the formal order of accuracy while the nonlinear part is responsible for capturing the discontinuities. For the nonlinear part, ω k is the weighting related to the relative smoothness of f on each set of candidate stencils. It is defined similarly as the original TENO scheme given in equations (7) to (10). For the linear part, C r k is the optimal weight which can be written in terms of the free parameters. The optimal weights are derived by comparing the expanded form of the linear part of equation (27) to the form of the minimized dispersion and controllable dissipation (MDCD) schemes given in equation (19). The resulting optimal weight is shown in table 1. Table 1: Optimal weights C k for MDCD-TENO6 scheme, r = 3 The key advantage of this framework is that this allow for specific tailoring of the dispersion or dissipation of the nonlinear scheme based on the absolute value of |γ disp | and |γ diss | respectively.

The hybrid scheme
Since the nonlinear mechanism of the shock capturing scheme may deteriorate the spectral properties of the scheme even in the smooth region, a hybrid scheme which combines the MDCD-TENO6 scheme with its linear counterpart is proposed. The reason for the use of the fourth order MDCD scheme lies in the fact that the stencil width of the linear MDCD scheme is the same as that of the TENO6 scheme. The hybridization strategy of the numerical flux is inherently based on the smoothness indicator of the TENO6 scheme. The numerical flux of the hybrid scheme is expressed aŝ where σ i+1/2 is the blending factor.f M DCD  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 where r c is a threshold parameter used to control the contribution of the linear scheme to the overall spectral property. r c is set to 0.03 based on numerical experiments and the optimal spectral properties desired. A lower value of r c would imply that a higher percentage of the numerical flux would be computed by the linear MDCD scheme. However, setting too low value of r c would lead to a loss of shock capturing properties. r i+1/2 is the blending function computed from the smoothness indicator, β k of the TENO scheme; The algorithm works as follows: n is set to 4 and 8 for one dimensional and multi-dimensional problems respectively. In equation (32), σ i+1/2 is taken to be 1 if |β k,max − β k,min | is less than ∆x n .

Spectral properties of the hybrid scheme
In this subsection, the spectral properties of the hybrid scheme with respect to the original TENO6 scheme will be analysed in both one dimensional and two dimensional wavenumber space. The proposed hybrid scheme will be denoted as the HYB-MDCD-TENO6 here-after. Since the dispersion and dissipation properties of the numerical scheme can be tailored separately, the value of γ disp will first be optimized for its dispersion properties. In the following part, the integrated error, E disp is plotted as a function of γ disp for different β values considered. This is shown in figure 8. After that, the spectral properties of the HYB-MDCD-TENO6 scheme for different values of γ disp is compared in figure 9.  Using the framework of the MDCD scheme ensures that both the dispersion and dissipation property can be separately controlled through γ disp and γ diss respectively. For stability reason, the dissipation should be non negative for all wavenumber. It can be seen in figure 9 that an optimal γ disp value of 0.0568 gives the optimal dispersion property for the hybrid scheme. In principle, the dissipation must be as small as possible but sufficiently large to ensure the stability of the simulation and the damping of the spurious numerical oscillation. Consequently, the optimal value of dissipation is problem dependent. An advantage of using the MDCD scheme is 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 as a function of γ diss for different β values  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 that this allows for the flexibility in adjusting the dissipation without affecting the dispersion property. In the instance that there is no numerical stability issues or spurious numerical oscillation, the dissipation of the hybrid scheme can be almost eliminated by choosing a very small |γ diss |. The flexibility in adjusting the dissipation makes it possible to develop a robust code for the simulation of problems with discontinuities. In principle, γ diss can be chosen in an adaptive way. For the purpose of comparison, the value of γ diss will be set to 0.0001 and 0.03. In two dimensional wavenumber space, the phase and amplitude error of the hybrid scheme is illustrated in figures 11 and 12.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 where the vectors U , F , G, H are expressed as:  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 65 The viscous stress tensor and heat flux are given by: (37) Other derived quantities can be calculated from the vector U as follows: The viscosity is computed based on Sutherlaw law: T ∞ is defined as the dimensional free stream temperature. For these cases, the third order Total Variation Diminishing (TVD) RK3 scheme is used for time advancement :   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 The spatial discretization methods employed for these cases are the WENO-JS, TENO or the hybrid MDCD-TENO6 schemes for the convective terms. For the use of the hybrid MDCD-TENO6 schemes, γ diss is set to 0.0001 for all numerical tests. The viscous stress and heat fluxes are computed through a generic 4th order central difference scheme. In order to avoid numerical oscillation, the local characteristic decomposition method is employed. For the construction of the left and right eigenvector, the Roe-averaged primitive variables are used. The constructed eigenvectors of the Euler equations are then used to transform the numerical flux from the physical to characteristic space and vice versa. Flux splitting and flux reconstruction are performed in the characteristic space. Finally, the reconstructed characteristic fluxes are then transformed back into physical space for the calculation of the derivatives. Further details can be found in [15,21]. Numerical results of the different schemes are compared with the original WENO5-JS scheme under a coarse grid. A fine grid case based on the WENO5-JS scheme with uniform grid in x and y is used for reference comparison. For each test cases, the boundary conditions are imposed by specifying appropriate values at the ghost points and the boundary nodes (if required) relevant to the specific boundary condition. With the use of ghost points, numerical schemes of the same stencil width can be used throughout the computational domain even at the boundary nodes. All two and three dimensional computations are performed using a parallel hybrid MPI-OMP Euler/Navier Stokes code while the one dimensional computations are performed using a serial Euler code.

One dimensional cases
In this section, several numerical tests for the one dimensional Euler equations are considered for the purpose of evaluating the performance of the non-linear scheme. For the one dimensional Euler equation, local characteristic decomposition method is employed, from which the eigenvectors at the cell interface are computed from the Roe-averaged primitive variables. Numerical results are compared between the WENO5-JS, WENO7-JS TENO5 and TENO6 and the hybrid MDCD-TENO6 schemes.

Sod problem
The initial condition of the Sod problem is given by:  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 The results are plotted in figure 13. It can be seen that the results of the hybrid scheme are closer to the reference results as compared to the baseline TENO6 scheme. Also, it can be seen that the both the TENO5 and TENO6 give closer results to the reference case than the baseline WENO5-JS and WENO7-JS scheme. As previously mentioned, the blending factor is a measure of the usage of the linear MDCD scheme. It can be clearly seen that σ remain 1.0 when the flow region is smooth and become less than 1.0 when the flow display a sharp change in gradient. By utilizing only the linear MDCD scheme in the smooth region of the flow, better spectral properties of the spatial discretization scheme can be obtained. In order to have a fair comparison of the numerical schemes in terms of their efficiency, the 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 performance metric used in [13] is adopted in this study. The efficiency metric is defined as follows: where is a measure of the error metric and T is the CPU time needed for the simulation to reach the solution at t = 0.14. The L 1 error norm definition is used for the error metric. The sum of the absolute difference between the reference solution and the solution for each scheme across all grid points is considered. Since E ∆ is not a constant for different sets of and T combination, the error metric is kept to a similar range of value for the different schemes by refining the grid when necessary. For each scheme, the grid is successively refined until the L 1 error is about 3.5 × 10 −1 and the results are summarized in   In table 3, it can be seen that the HYB-MDCD-TENO6 scheme is the most efficient scheme among the different spatial discretization schemes. The HYB-MDCD-TENO6 scheme can achieve similar magnitude of L 1 error with lesser grid points.

Lax problem
The initial condition of the Lax problem is given by: The results are plotted in figure 14. It can be seen that the results of the hybrid scheme are slightly closer to the reference results as compared to the baseline TENO6 scheme. Also, it can be seen that the both the TENO5 and TENO6 schemes give better results than the baseline WENO5-JS and WENO7-JS scheme. For the plot of the blending factor, it can be clearly seen that σ remain 1.0 when the flow region is smooth and become less than 1.0 when the flow display a significant change in gradient .   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65

Shu Osher problem
This test problem represents an interaction between a Mach 3 shockwave and a sine entropy wave which generates a flow field with both small scale structures and discontinuities. The initial condition of the Lax problem is given by:  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 TENO6 scheme. Also, it can be seen that the both the TENO5 and TENO6 schemes give better results than the baseline WENO5-JS and WENO7-JS scheme. The plot of the blending factor display similar trends as compared to previous cases. σ remain 1.0 when the flow region is smooth and become less than 1.0 when the flow display a significant change in gradient.

2D Double Mach Reflection of a strong shock
The double Mach reflection test is a close representation of planar shock reflection in the air from a wedge. Instead of orientating the geometry to create an angled surface, the shock is orientated at an angle and the wall is kept parallel to the bottom domain. For this case, only the two dimensional Euler equation is considered. It is a widely used benchmark to test the ability of the numerical scheme to capture shock and resolve the small scale structures. For the present simulation, the computational domain is taken to be [0,4] x [0,1]. From x >= 1 6 , the lower boundary is set to be a reflecting free slip wall. At the initial state, a right moving Mach 10 shock is placed at x = 1 6 with an incidence angle of 60 • to the x-axis. For the top boundary, the fluid variables are defined to exactly follow the evolution of the Mach 10 shockwave. Mathematically the top boundary is defined as: (45) As the shock propagates from left to right, a self similar shock structure with two triple points evolve. At the primary triple point, the incident shock, reflected shock and Mach stem intersects. A secondary triple points is formed from the intersection of the reflected shock, bowed Mach stem and a secondary reflected shock [10]. A primary slip line emanates from both triple points. In figure 16, it can be seen that the cases using the TENO5 and TENO6 scheme are able to produce finer scale structure along the primary slip line (curled flow structure) as compared to the baseline WENO5-JS scheme. As compared to the fine grid WENO5-JS case, both the TENO5 and TENO6 cases are able to produce either similar or finer scale structures despite using a coarser grid. The hybrid schemes are able to produce even finer scale structures at the primary slip line as compared to the baseline TENO6 scheme. A plot of the blending factor in the x and y direction for the hybrid scheme is shown in figure 17. It can be seen that at region where there is sharp change in the gradient of the flow, the blending factor approaches zero. At region where the gradient of the flow is smooth, the blending factor approaches one. At region of intermediate gradient, the blending factor lies between zero and one.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65

2D Euler -Shock Vortex Interaction
This problem describes the interaction of a stationary shock and a vortex. Details of the initial and boundary conditions were obtained from [15]. The computational domain is set to [0,2] x [0,1]. A Mach 1.1 shock is positioned at x = 0.5 and normal to the x− axis. Its left and right states are given as (1.1691, 1.1133, 0, 1.245), otherwise.
A small vortex is superposed on the left state region and its center is located at (x c , y c ) = (0.25, 0.5). The initial condition to the left of the normal shock is given by: where τ = r/r c and r = (x − x c ) 2 + (y − y c ) 2 . = 0.3 is a measure of the strength of the vortex. α = 0.204 controls the decay rate of the vortex and r c = 0.05 is the critical radius at which the vortex has the maximum strength. A grid of (N x , N y ) = (250, 100) is used. The reference results is computed on a refined grid of (N x , N y ) = (1000, 400). Both inflow and outflow values at the ghost points are linearly extrapolated from the interior points. Top and bottom sides are imposed symmetry boundary condition.
The results are compared at t = 0.6 and they are illustrated in figure 18. Density distribution of different results at y=0.5 are illustrated in figure 19. It can be seen that the different spatial discretization methods give very similar results although the TENO5, TENO6 and HYB-MDCD-TENO6 schemes has slightly more numerical oscillation at region close to the normal shock. It can be seen that all methods capture the shock and vortex well. It can be seen that the HYB-MDCD-TENO6 scheme is able to achieve sharper normal shock profile and lower vortex core density as compared to the baseline TENO6 schemes. This is because of the lower numerical dissipation. Also, it can be seen that both the TENO5 and TENO6 cases give better results than the WENO5-JS and WENO7-JS cases at these regions. A plot of the blending factor in the x and y direction for the hybrid scheme is shown in figure 20.
An equally spaced grid of (N x , N y ) = (500, 100) is used. The reference result is calculated on a refined grid of (N x , N y ) = (2000, 400). Density contours are shown in figure 21 at t = 120. It can be observed that the WENO5-JS, TENO5, TENO6 and HYB-MDCD-TENO6 cases obtain similar results which is very comparable to the reference result. However, it should be noted that the self similar structure of the downstream vortices of the reference case showcases smaller scale structure close to its vortex core. As compared to the baseline TENO6 case, the HYB-MDCD-TENO6 case showcases presence of smaller scale structure within its downstream vortex core from X 160 onward. Also, the size of the downstream vortices are comparatively larger. A plot of the blending factors, σ x and σ y of the first flux component in each spatial directions is shown in figure 22. The slanted 'white line' region where the blending factor is close to zero corresponds to the location of the impinging oblique shock interacting with the mixing layer.
6.4 Three dimensional case