A L2-norm regularized incremental-stencil WENO scheme for compressible flows

TL;DR

This paper introduces a novel L2-norm regularized incremental-stencil WENO scheme that enhances robustness and accuracy in simulating compressible flows with complex structures by adaptively modulating stencil weights.

Contribution

It proposes a new L2-norm regularization approach integrated into WENO weighting, improving accuracy and efficiency in flow simulations over existing methods.

Findings

Achieves high robustness in complex flow simulations.

Resolves fine flow structures effectively.

Demonstrates improved computational efficiency.

Abstract

For the simulation of compressible flow with a broadband of length scales and discontinuities, the WENO schemes using incremental stencil sizes other than uniform ones are promising for more robustness and less numerical dissipation. However, in smooth region, large weights may be assigned to smaller stencils due to the lack of high-order derivatives in the smoothness indicator compared with that of larger stencils, and may degrade the order accuracy and lead too much numerical dissipation to resolve fine flow structures. In order to cope with this drawback, based on the stencil selection of WENO-IS [Wang et al., IJMF 104 (2018): 20-31], we propose a L2-norm regularized incremental-stencil WENO scheme. In this method, a new L2-norm regularization is introduced into the WENO weighting strategy to modulate the weights of incremental-width stencils by taking account the L2-norm error term. A high-order non-dimensional discontinuity detector is then utilized as the regularization parameter for adaptive control. In addition, a hybrid method is adopted to further improve the performance and the computational efficiency. A number of benchmark cases suggest that the present scheme achieves very good robustness and fine-structure resolving capabilities.

Tables

Table 1. Table 1: Computational time for double Mach reflection of a strong shock on the 512 × 128 512 128 512\times 128 grid.

	WENO-JS	WENO-HY	WENO-Z	Present
CPU time (s)	$1721$	$1146$	$1986$	$1446$

Table 2. Table 2: Computational time for viscous shock tube problem on the 250 × 125 250 125 250\times 125 grid.

	WENO-JS	WENO-HY	WENO-Z	Present
CPU time (s)	$1690$	$1176$	$1740$	$1343$

Equations

$$\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}=0,$$

$$\frac{d{u}_j}{dt}=-\frac{1}{\Delta x}\left(h_{j+1/2}-h_{j-1/2}\right),$$

$$f\left[u\left(x\right)\right]=\frac{1}{\Delta x}{\int }_{x-\Delta x/2}^{x+\Delta x/2}h(\xi )d\xi .$$

$$\frac{d{u}_j}{dt}\approx -\frac{1}{\Delta x}\left({\hat{f}}_{j+1/2}-{\hat{f}}_{j-1/2}\right).$$

$${\hat{f}}_{j+1/2}=\sum _{k=0}^3{\omega }_k{\hat{f}}_{k,j+1/2},$$

$$\begin{array}{rl}& {\hat{f}}_{0,j+1/2}=\frac{1}{2}{f}_j+\frac{1}{2}{f}_{j+1},\\ & {\hat{f}}_{1,j+1/2}=-\frac{1}{2}{f}_{j-1}+\frac{3}{2}{f}_j,\\ & {\hat{f}}_{2,j+1/2}=\frac{1}{3}{f}_j+\frac{5}{6}{f}_{j+1}-\frac{1}{6}{f}_{j+2},\\ & {\hat{f}}_{3,j+1/2}=\frac{1}{3}{f}_{j-2}-\frac{7}{6}{f}_{j-1}+\frac{11}{6}{f}_j.\end{array}$$

$${\hat{f}}_{k,j+1/2}=h_{j+1/2}+A_k\Delta x^{r_k}+O\left(\Delta x^{r_k+1}\right),$$

$${\omega }_k=\frac{{\alpha }_k}{{\sum }_{s=0}^3{\alpha }_s},{\alpha }_k=d_k\left(1+\frac{{\tau }_5}{{\beta }_k+\epsilon }\right),$$

$${\beta }_k=\sum _{l=1}^{r_k-1}\Delta x^{2l-1}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left(\frac{d^l}{d{x}^l}{\hat{f}}_k\left(x\right)\right)}^{2}dx.$$

$${\tau }_5=\sum _{l=3}^4\Delta x^{2l-1}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left(\frac{d^l}{d{x}^l}{\hat{f}}_k\left(x\right)\right)}^{2}dx,$$

$$\begin{array}{rl}& {\beta }_0={\left(f_{j+1}-f_j\right)}^2,\\ & {\beta }_1={\left(f_j-f_{j-1}\right)}^2,\\ & {\beta }_2=\frac{13}{12}{\left(f_j-2f_{j+1}+f_{j+2}\right)}^2+\frac{1}{4}{\left(3f_j-4f_{j+1}+f_{j+2}\right)}^2,\\ & {\beta }_3=\frac{13}{12}{\left(f_{j-2}-2f_{j-1}+f_j\right)}^2+\frac{1}{4}{\left(f_{j-2}-4f_{j-1}+3f_j\right)}^2,\\ & {\tau }_5=\frac{13}{12}{\left(f_{j+2}-4f_{j+1}+6f_j-4f_{j-1}+f_{j-2}\right)}^2+\frac{1}{4}{\left(f_{j+2}-2f_{j+1}+2f_{j-1}-f_{j-2}\right)}^2,\end{array}$$

$$\begin{array}{rl}& {\beta }_0={f^{\prime }}_j^2\Delta x^2+{f^{\prime }}_j{f^{\prime }{}^{\prime }}_j\Delta x^3+\left(\frac{1}{4}{f^{\prime }{}^{\prime }}_j^2+\frac{1}{3}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^4+\left(\frac{1}{12}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }{}^{\prime }}_j+\frac{1}{6}{f^{\prime }{}^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^5+O\left(\Delta x^6\right),\\ & {\beta }_1={f^{\prime }}_j^2\Delta x^2-{f^{\prime }}_j{f^{\prime }{}^{\prime }}_j\Delta x^3+\left(\frac{1}{4}{f^{\prime }{}^{\prime }}_j^2+\frac{1}{3}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^4-\left(\frac{1}{12}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }{}^{\prime }}_j+\frac{1}{6}{f^{\prime }{}^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^5+O\left(\Delta x^6\right),\\ & {\beta }_2={f^{\prime }}_j^2\Delta x^2+\left(\frac{13}{12}{f^{\prime }{}^{\prime }}_j^2-\frac{2}{3}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^4+\left(-\frac{1}{2}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }{}^{\prime }}_j+\frac{13}{6}{f^{\prime }{}^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^5+O\left(\Delta x^6\right),\\ & {\beta }_3={f^{\prime }}_j^2\Delta x^2+\left(\frac{13}{12}{f^{\prime }{}^{\prime }}_j^2-\frac{2}{3}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^4-\left(-\frac{1}{2}{f^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }{}^{\prime }}_j+\frac{13}{6}{f^{\prime }{}^{\prime }}_j{f^{\prime }{}^{\prime }{}^{\prime }}_j\right)\Delta x^5+O\left(\Delta x^6\right),\\ & {\tau }_5={f^{\prime }{}^{\prime }{}^{\prime }}_j\Delta x^6+\frac{13}{12}{f^{\prime }{}^{\prime }{}^{\prime }{}^{\prime }}_j\Delta x^8+O\left(\Delta x^{10}\right).\end{array}$$

$$\frac{{\tau }_5}{{\beta }_k+\epsilon }=O\left(\Delta x^4\right),k=0,1,2,3.$$

$${\alpha }_0=d_0\left(1+\frac{{\tau }_5}{{\beta }_0+\epsilon }\cdot \frac{{\tau }_5}{{\beta }_{01}+\epsilon }\right),{\alpha }_1=d_1\left(1+\frac{{\tau }_5}{{\beta }_1+\epsilon }\cdot \frac{{\tau }_5}{{\beta }_{01}+\epsilon }\right),$$

$${\beta }_{01}=\frac{13}{12}{\left(f_{j-1}-2f_j+f_{j+1}\right)}^2+\frac{1}{4}{\left(f_{j-1}-f_{j+1}\right)}^2$$

$$\frac{{\tau }_5}{{\beta }_0+\epsilon }\cdot \frac{{\tau }_5}{{\beta }_{01}+\epsilon }\gg \frac{{\tau }_5}{{\beta }_2+\epsilon },\frac{{\tau }_5}{{\beta }_1+\epsilon }\cdot \frac{{\tau }_5}{{\beta }_{01}+\epsilon }\gg \frac{{\tau }_5}{{\beta }_2+\epsilon },$$

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}\left(\mathbf{U}\right)}{\partial x}=0,$$

$$\bm{U}=\left[\begin{array}[]{c}\rho\\ \rho u\\ {E}\end{array}\right],\bm{F}\left(\bm{U}\right)=\left[\begin{array}[]{c}\rho u\\ \rho u^{2}+p\\ {\left(E+p\right)u}\end{array}\right].\\$$

$$\frac{d{\mathbf{U}}_j}{dt}=-\frac{1}{\Delta x}\left({\hat{\mathbf{F}}}_{j+1/2}-{\hat{\mathbf{F}}}_{j-1/2}\right),$$

$$v_m^s={\mathbf{L}}_{j+1/2}^s\cdot {\mathbf{U}}_m,g_m^s={\mathbf{L}}_{j+1/2}^s\cdot {\mathbf{F}}_m,$$

$$f_m^{s,\pm }=\frac{1}{2}\left(g_m^s\pm {\alpha }^s{v}_m^s\right),$$

$$f_{j+1/2}^{s,+}=\sum _{k=0}^3{\omega }_k^{+}{f}_{k,j+1/2}^{s,+},f_{j+1/2}^{s,-}=\sum _{k=0}^3{\omega }_k^{-}{f}_{k,j+1/2}^{s,-},$$

$${\hat{f}}_{j+1/2}^s=f_{j+1/2}^{s,+}+f_{j+1/2}^{s,-}.$$

$${\hat{\mathbf{F}}}_{j+1/2}=\sum _{s=0}^2{\mathbf{R}}_{j+1/2}^s{\hat{f}}_{j+1/2}^s.$$

$$E_k=\frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left({\hat{f}}_k\left(x\right)-f_k\left(x\right)\right)}^{2}dx,$$

$$\begin{array}{rl}& E_0=\frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left({\hat{f}}_0\left(x\right)-f_{01}\left(x\right)\right)}^{2}dx,\\ & E_1=\frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left({\hat{f}}_1\left(x\right)-f_{01}\left(x\right)\right)}^{2}dx,\\ & E_2=\frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left({\hat{f}}_2\left(x\right)-f_2\left(x\right)\right)}^{2}dx,\\ & E_3=\frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left({\hat{f}}_3\left(x\right)-f_3\left(x\right)\right)}^{2}dx,\end{array}$$

$$\begin{array}{rl}& E_0=\frac{1}{45}\left(f_{j+1}^2-4f_j{f}_{j+1}+2f_{j-1}{f}_{j+1}+4f_j^2-4f_{j-1}{f}_j+f_{j-1}^2\right),\\ & E_1=\frac{1}{45}\left(f_{j+1}^2-4f_j{f}_{j+1}+2f_{j-1}{f}_{j+1}+4f_j^2-4f_{j-1}{f}_j+f_{j-1}^2\right),\\ & E_2=0,\\ & E_3=0.\end{array}$$

$${\omega }_k^r=\frac{{\alpha }_k^r}{{\sum }_{s=0}^3{\alpha }_s^r},{\alpha }_k^r=d_k\left(1+\frac{\lambda {\tau }_5}{\lambda {\beta }_k+E_k+\epsilon }\right),$$

$$\begin{array}{rl}& E_0=E_1=\frac{1}{45}\left({f^{\prime \prime }}_j\Delta x^4+\frac{1}{6}{f^{\prime \prime }}_j{f^{\prime \prime \prime \prime }}_j\Delta x^6+O\left(\Delta x^6\right)\right),\\ & E_2=E_3=0.\end{array}$$

$${\sigma }^s={\left(\frac{\Delta v_{j+\frac{1}{2}}^s}{\tilde{\rho }}\right)}^2,$$

Full text

††journal: Journal of Computational Physics

1 Introduction

The simulation of compressible flow is challenging when there present a broadband of length scales and discontinuities. The numerical schemes must be of low dissipation in smooth region without damping the fine flow structures, and be able to introduce sufficient dissipation in vicinity of discontinuities for suppressing non-physical oscillations. In order to address this issue, many high-order shock capturing schemes, like total variation diminishing (TVD)[1], essentially non-oscillatory (ENO)[2] and weighted essentially non-oscillatory (WENO) [3] schemes, have been developed. Among them, the WENO-JS scheme proposed by Jiang and Shu [4] exhibits high-order accuracy and robust shock-capturing properties and have been extensively employed in direct numerical simulations (DNS) and large eddy simulations (LES). In spite of these advantages, the classic WENO-JS scheme still faces several shortcomings. It suffers from order-degeneration, e.g. not able to recover the formal order of accuracy near critical points where the first or higher derivatives vanish [5], and is too dissipative for sustaining fine structures in turbulence or aeroacoustic flows [6]. Furthermore, the computation may fail when very strong discontinuities exist or multiple discontinuities are too close to each other [7].

In order to overcome the drawbacks of order-degeneration and excessive dissipation of the classical WENO-JS scheme, one promising approach is to modify the weighting strategy, e.g. the adaptation mechanism, by increasing the weight of the optimal upwind linear scheme in relatively smooth regions. Following this idea, various schemes have been developed, such as WENO-M scheme [5] with a corrective mapping applied to the classical WENO weights, WENO-Z scheme [8, 9] with a new high-order smoothness indicator and WENO-RL or WENO-RLTV [10] scheme switching on the optimal weights when the magnitude of smoothness indicators are comparable. Another approach is to introduce the contribution of the downwind stencil and use optimal weights for central scheme instead of upwind scheme as pointed in Refs. [6, 11, 12]. Take WENO-CU6 scheme [12] as an example, Hu et al. employed the 6th-order central scheme as the optimal linear scheme and introduced a new smooth indicator for the adaptation between central and upwind stencils. Later, this scheme was extended to WENO-CU6-M with a physically-motivated scale separation approach for implicit large eddy simulation (ILES) [13]. While WENO-CU6-M works well for incompressible and compressible turbulence simulation with low dissipation and maintains the shock-capturing capability, it produces spurious waves because of the non-dissipative optimal scheme. In addition, as demonstrated in Ref. [14], one version of WENO-CU6-M scheme fails in the benchmark test cases ”Interacting blast waves” and ”Double Mach reflection” and suggests inferior robustness than the classic WENO-JS scheme. Fu et al. [14] proposed a family of TENO scheme which dynamically assembles a set of low-order candidate stencils with incrementally increasing stencil width. The combination of low-order approximation polynomials with incremental stencils makes TENO scheme considerably more robust. However, positivity-preserving methods, such as those in Ref. [15, 16], are still required for TENO schemes when high density ratio or strong shock waves present.

Recently, Wang et al. [17] developed a 5th-order incremental-stencil WENO (denote as WENO-IS) scheme for multi-phase flows. Since both 2- and 3-point stencils are used, in the region with strong discontinuities, the reconstruction would prefer to choose the 2-point stencils, which makes the scheme, compared to the classic WENO-JS schemes, very robust for material interface with high density ratio and strong shock waves. In their weighting strategy, a slight modification is introduced to decrease the weights for 2-point stencils in relatively smooth region and to avoid excessive dissipation and order degeneration near critical points. However, this scheme is still too dissipative to capture fine flow structures, as also shown in Section 4, due to the too much weights for the 2-point stencils.

In this paper, we propose a L2-norm regularized incremental-stencil WENO scheme following the idea of neural-network optimization [18]. Similar to WENO-IS scheme, both incremental 2- and 3-point stencils are used for the 5th-order reconstruction. The difference is that a L2-norm regularization is introduced to address the problem of order degeneration and excessive dissipation. First, the L2-norm errors between the approximations from 2-point stencils and a classical 3-point stencil are utilized to modulate the weights of the incremental stencils. Then, a high-order non-dimensional discontinuity detector is adopted as the regularization parameter to control the influence of the error term. With a further hybridization with a linear scheme, less numerical dissipation in very smooth region and considerable higher computational efficiency are achieved. The reminder of the paper is organized as follows. Section 2 gives a brief review of the WENO-IS scheme following Wang et al. [17] and the application of WENO-IS to Euler equations is given. In Section 3, the L2-norm regularization and control parameter are proposed. Benchmark test cases with strong discontinuities and fine flow structures are presented in Section 4, and brief concluding remarks are given in the last Section 5.

2 WENO-IS scheme

In this section, we briefly review the 5th-order WENO-IS scheme for hyperbolic conservation laws following Ref. [17]. First, we illustrate the details of WENO-IS scheme based on the discretization of one-dimensional scalar hyperbolic conservation law. Then, its application in Euler equations is given.

2.1 Semi-discretization of a hyperbolic conservation law

For simplicity, we consider the following one-dimensional scalar hyperbolic conservation law

[TABLE]

where $u$ denotes the conservative variable and $f\left(u\right)$ is the flux function. Here, the characteristic velocity is assumed to be positive with $\frac{\partial{f\left({u}\right)}}{{\partial u}}>0$.

On a uniform grid, we denote $x_{j}=j\Delta x$, where $\Delta x$ is the grid spacing. The quantities at the point $x_{j}$ are identified by the subscript $j$. The semi-discretized form of Eq. (1) using the method of lines can be approximated by a conservative finite difference scheme as

[TABLE]

where $h_{j\pm 1/2}=h\left(x_{j\pm 1/2}\right)$ and is implicitly defined by

[TABLE]

Numerically, $h_{j\pm 1/2}$ are approximated by high-order polynomial interpolations $\widehat{f}_{j\pm 1/2}$ or numerical fluxes and Eq. (2) can be expressed as

[TABLE]

2.2 Candidate stencils with incremental width

For 5th-order WENO-IS scheme, the 5-point stencil $S_{5}$ as shown in Fig. 1 is served as the full stencil for the evaluation of numerical fluxes $\widehat{f}_{j+1/2}$. Different from the classic WENO-JS scheme, which subdivides the full stencil into three stencils with the same width, the WENO-IS scheme uses four stencils with incremental width including two 2-point and two 3-point stencils, e.g. $S_{0}$, $S_{1}$, $S_{2}$ and $S_{3}$ as shown in Fig. 1, for reconstruction.

The numerical flux $\widehat{f}_{j+1/2}$ is computed by a convex combination of the four candidate stencils fluxes as

[TABLE]

where $\widehat{f}_{k,j+1/2}$ and $\omega_{k}$ are the reconstructed fluxes from candidate stencils and their non-linear weights, respectively. These reconstructed fluxes are

[TABLE]

From Taylor series expansion, we can obtain

[TABLE]

where $A_{k}$, $k=0,1,2,3$ are independent of $\Delta x$ and $r_{k}$ denote the points of the candidate stencil.

The non-linear weights $\omega_{k}$ are given by

[TABLE]

where coefficients $d_{0}=\frac{4}{10},d_{1}=\frac{2}{10},d_{2}=\frac{3}{10},d_{3}=\frac{1}{10}$ are the optimal weights and $\varepsilon$ is a small positive parameter to prevent division by zero and we choose $\varepsilon=10^{-20}$ in this paper. $\beta_{k}$ are the smoothness indicators of each candidate stencil, while $\tau_{5}$ is the global reference smoothness indicator on the full stencil. Following Jiang and Shu  [4], the smoothness indicators at $S_{k}$ are evaluated by

[TABLE]

The global reference smoothness indicator is obtained by

[TABLE]

which only contains the high-order component of the smoothness indicator for the full stencil [14, 17]. The formulations of $\beta_{k}$ and $\tau_{5}$ in terms of cell-centered fluxes $f_{j}$ can be expressed as

[TABLE]

and their Taylor series expansions at $x_{j}$ are

[TABLE]

Substituting Eq. (12) to Eq. (8), we can observe that

[TABLE]

Therefore, the condition of $\omega_{k}-d_{k}=O\left(\Delta x^{4}\right)$ is satisfied at the region without critical points, implying that the WENO-IS scheme has the expected 5th-order convergence from the analysis in Ref. [17].

2.3 Weighting strategy

From Eq. (9) and Eq. (11), we can observe that the smoothness indicators of 2-point stencils are evaluated from the 1st-order derivatives of linear functions, while these of 3-point stencils are from the 1st and 2nd-order derivatives of parabolic functions. Due to the lack of 2nd-order derivative term, the smoothness indicator of 2-point stencils can be much smaller than that of stencils with 3-point especially near critical point and larger weights will be posed on 2-point stencils, which may cause order degeneration eventually. In order to address this issue, Wang et al. [17] made a slight modification of the weighting strategy in Eq. (8) as

[TABLE]

where

[TABLE]

is the smoothness indicator of $S_{01}$, an original 3-point stencil of WENO-JS as shown in Fig. 1, and is obtained from 1st and 2nd-order derivatives. This modification decreases the weights of 2-point stencils considerably and avoids order degeneration at critical point.

In relatively smooth region, the effectiveness of this modification has been confirmed by a number of test cases in Ref. [17]. The WENO-IS scheme achieves the formal order of accuracy and produces less errors than classical WENO-JS scheme. However, in the region with shocks or high wavenumber waves, this modification may not decrease the weights of 2-point stencils sufficiently, which can still lead to excessive numerical dissipation. Suppose there is a shock locating in stencil $S_{3}$ and the solutions in stencils $S_{0}$, $S_{1}$ and $S_{2}$ are smooth (see Fig. 1). Under this condition, $\tau_{5}$ is comparable to $\beta_{3}$ and $\tau_{5}=O\left(1\right)$, while $\beta_{k}=O\left(\Delta x^{2}\right)$(k=0,1,2,01) are all small numbers. It is straightforward to show that

[TABLE]

which means larger weights will be posed on 2-point stencils. Same conclusion can be drawn when the shock locates in other stencils. Thus, this modification enlarges the numerical dissipation near the discontinuities. Although larger numerical dissipation improves the robustness of the scheme, this modification imposes too much excessive numerical dissipation to resolve the fine flow structures with broadband length scales well (also see the cases in section 4).

Note that TENO scheme  [14] avoid the problem of large weight posed on low-order stencils by a ENO-like strategy. However, the numerical stability of TENO scheme is inferior to WENO-JS hence also WENO-IS scheme.

2.4 Application in Euler equations

For simplicity, the detailed procedure on the application of WENO-IS for Euler equations is illustrated based on the one-dimensional Euler equations of gas dynamics

[TABLE]

where

[TABLE]

Here, $\rho$ is density, $p$ is pressure and $u$ denotes the velocity in $x-$direction. $E=p/\left(\gamma-1\right)+1/2\rho u^{2}$ is total energy per unit volume with ideal-gas equation of states used.

Within a uniform grid, at the point $x_{j}$, the semi-discretized form of Eq. (17) gives

[TABLE]

where $\hat{\bm{F}}_{j\pm 1/2}$ are the numerical fluxes at $x_{j\pm 1/2}$. With a typical characteristic-wise method, the numerical fluxes are reconstructed within the local characteristic flied. The procedure to obtain $\hat{\bm{F}}_{j+1/2}$ includes the following steps:

(1) At $x_{j+1/2}$, evaluate the Jacobian matrix $\bm{A}=\partial\bm{F}/\bm{U}$ at a Roe-average state.

(2) Compute the left and right eigenvector matrix $\bm{L}_{j+1/2}^{s}$, $\bm{R}_{j+1/2}^{s}\left(s=1,2,3\right)$ as well as the eigenvalues $\lambda_{j+1/2}^{s}\left(s=1,2,3\right)$ of the Roe-average Jacobian matrix $\bm{A}$ .

(3) Within each candidate stencils, transform the physical fluxes and conservative variables into characteristic filed as

[TABLE]

where $j-2<m<j+3$ and $s=1,2,3$.

(4) Carry out flux splitting in characteristic space. We can obtain

[TABLE]

where $\alpha^{s}=\left|\lambda_{j+1/2}^{s}\right|$ for Roe flux (RF) splitting. Otherwise, $\alpha^{s}=max\left|\lambda_{l}^{s}\right|$ for Lax-Friedrichs flux (LF) splitting with $l$ representing the entire computational domain or for local Lax-Friedrichs flux (LLF) splitting when $l$ represents the local stencil.

(5) Then, reconstruct the numerical flux with WENO-IS scheme, which gives

[TABLE]

where $f_{k,j+1/2}^{s,+}$ and $\omega_{k}^{+}$ are the positive fluxes by candidate stencils and their non-linear weights as defined in Eq. (6) and Eq. (8), respectively. The negative term $f_{k,j+1/2}^{s,-}$ and $\omega_{k}^{-}$ can be computed in the similar way according to the symmetry at the point $x_{j+1/2}$. The numerical flux in each characteristic space can be finally evaluated by

[TABLE]

(6) At last, transform the numerical flux back into physical space as

[TABLE]

3 A L2-norm regularized WENO-IS scheme

In Ref. [18] (their section 3.2.1), in order to find a way to reduce overfitting of the neural-network, an optimization technique, known as L2-norm regularization, is introduced. The essential idea of L2-norm regularization is to add an extra term, called the regularization term, to the original cost function which only contains the L2-norm error term, with a hyper-parameter to control the relative amount of regularization. This technique can be viewed as a way of compromising between the original L2-norm error term and the L2-norm regularization term, which makes the neural-network better at generalizing beyond the training data.

In the classical WENO schemes, as an analog of optimization technique, we can find that only the regularization terms based on the smoothness indicators are utilized. The L2-norm error term is neglected due to the same order of accuracy for each candidate stencil. However, for WENO-IS scheme, the reconstruction accuracy of 2-point stencil is different from that of the 3-point stencil, therefore, neglecting the L2-norm error terms is not suitable anymore. To cope with the discrepancy of numerical accuracy between 2-point and 3-point stencils, an optimization strategy, L2-norm regularization, is introduced to the WENO weighting strategy.

3.1 L2-norm regularization

For the conservation law in Eq. (1), a form of L2-norm error term is

[TABLE]

which is obtained by the difference between reconstruct flux and reference flux. For the 5-th order WENO-IS scheme, the three reconstructed fluxes of classic WENO-JS scheme are chosen as the reference fluxes. Therefore, the L2-norm error term of each stencil for WENO-IS scheme can be evaluated by

[TABLE]

where $f_{01}\left(x\right)$, $f_{2}\left(x\right)$ and $f_{3}\left(x\right)$ are the reconstructed fluxes of the corresponding stencils $S_{01}$, $S_{2}$ and $S_{3}$ as shown in Fig. 1. Since the reconstructed fluxes and the reference fluxes of stencil $S_{2}$ and $S_{3}$ are the same, their L2-norm error terms equal to zero accordingly. Therefore, the above error terms yield

[TABLE]

With L2-norm regularization, the weighting strategy of WENO-IS scheme is modified as

[TABLE]

where $\lambda$ is known as the regularization parameter and we will discuss it in later subsection. The Taylor series expansions at $x_{j}$ of the L2-norm error term $E_{k}$ as shown in Eq. (27) are

[TABLE]

It is straightforward to check from Eq. (12) that $\omega_{k}^{r}-d_{k}=O\left(\Delta x^{4}\right)$ is still satisfied for the present weighting strategy. Hence, with L2-norm regularization the modified WENO-IS scheme also has the expect 5th-order convergence at non-critical points. Note that, this modification doesn’t solve the order degeneration problem at critical points. Although some weighting strategies can be applied to cope with this case as noted by Borges et al. [8], it will increase numerical dissipation inevitably. Alternatively, in the present paper, we introduce a hybrid scheme to resolve the order degeneration problem, which will be discussed in section 3.3.

3.2 An adaptive weights control

It is obvious that the effect of the regularization is to make up the discrepancy of numerical accuracy between 2-point and 3-points stencils so as to reduce the weights of 2-point stencils. The relative importance of the regularization term and the added L2-norm error term depends on the value of $\lambda$. Here, we introduce a non-dimensional discontinuity detector following Hu et al. [19] to devise an adaptive regularization parameter. Due to the design of non-dimensional discontinuity detector is based on Euler equations, in this section, the detailed procedure is illustrated using one dimensional Euler equations as presented in section 2.3.

The non-dimensional discontinuity detector is defined by

[TABLE]

where

[TABLE]

and ${\tilde{\rho}}$ is the Roe-average density corresponding to the Jacobian matrix $\bm{A}$ at $x_{j+\frac{1}{2}}$.

We define a regularization parameter by

[TABLE]

where $c$ is a positive constant and we set $c=1$ in this work. When the full 6-point stencil contains discontinuities, the regularization parameter becomes a large number, which makes the original smoothness indicator the leading term and 2-point stencils will be assigned larger weights. Note that, the relatively large weights of 2-point stencils evaluated by the present method are less than these of the original WENO-IS scheme under this condition. As will be shown in section 4, with this treatment the scheme is still robust even for strong discontinuities and less numerical dissipation is imposed. Otherwise, the reconstruction prefer to choose 3-point stencils in the flow field with fine flow structures or the smooth region where $\lambda$ is small.

3.3 The hybrid scheme

Although the added L2-norm error term and the designed adaptive parameter have the advantages as discussed above, the modified WENO-IS scheme is not suitable for smooth region due to the preference of 3-point stencil of the reconstruction when $\lambda$ is very small. Furthermore, compared with its optimal linear scheme, the calculation of the modified WENO-IS scheme is time-consuming. To overcome these limitations, with the discontinuity detector presented in above section, it is easy to introduce a hybrid scheme flowing Hu et al. [19] to achieve better accuracy in the smooth region and improve computational efficiency. The numerical flux of this hybrid scheme is switching from that of the nonlinear WENO-IS scheme and its optimal linear upwind scheme, which is given by

[TABLE]

where $\sigma_{j+1/2}$ equals to one in smooth region otherwise zero in vicinity of discontinuities. Here, the numerical flux $\hat{\bm{F}}_{j+1/2}^{UPS}$ of the optimal linear upwind scheme is evaluated by

[TABLE]

It is obvious that this linear scheme omits the computation of the nonlinear weights as presented in Eq. (8) and  (22), and decreases by $5/6$ the characteristic-projection operations of Eq. (20) as noted in Ref. [19]. Therefore, the hybridization will save much computational time.

The threshold is defined as

[TABLE]

where $L$ is the characteristic length scale of the problem, $\alpha$ is a positive integer and $C$ is a positive constant. In this paper, we choose $C=1$ and $\alpha=3$. When $\sigma^{s}<\epsilon$, we think the flow field is smooth enough and numerical flux is evaluated by the optimal linear upwind scheme. Otherwise, it is obtained by the regularized WENO-IS scheme.

4 Numerical cases

In order to assess the performance of the proposed scheme, various test cases are considered, including shock tube problems, shock/entropy wave interaction problems and a problem involving very strong discontinuities. In this section, the local Lax-Friendrichs flux is adopted for flux splitting in one dimensional cases, while Lax-Friendrichs flux is used for two dimensional cases if not mentioned otherwise. The third order TVD Runge-Kutta scheme is utilized for time integration with a CFL number of 0.5. In the following, ”WENO-JS” denotes the classical 5th-order WENO scheme [4], ”WENO-Z” represents the scheme devised in Ref. [8], ”WENO-IS” denotes the original scheme in Ref. [17], ”WENO-HY” is the hybrid method devised by Hu et al. [19], ”Present” is the proposed scheme of this paper, and ”Exact” denotes the theoretical or convergent solution.

4.1 Propagation of broadband sound waves

This case is taken form Ref. [11]. A propagation of sound wave packet which contains various length scales of acoustic turbulent structures is computed. The initial condition is given as

[TABLE]

where

[TABLE]

is the energy spectrum which reaches its maximum at $k=k_{0}$. Here, $\phi_{k}$ is a random number between [math] and $1$, $\varepsilon=0.001$, $\gamma=1.4$ and the sound speed $c=\sqrt{\gamma p/\rho}$. The computational domain is $x\in\left[0,1\right]$ and period boundary conditions are imposed at the boundaries. The computing is executed on a 128-point grid for one period of time with CFL number of $0.2$. In this case, we choose $k_{0}=12$, which means most part of the energy is concentrated on high wavenumbers [11]. The numerical result of pressure distribution is shown in Fig. 2. We can find that the original WENO-IS scheme produces more excessive numerical dissipation for the high wavenumber waves although it has better numerical stability property as stated in Ref. [17]. However, with the modification of this paper, the present method achieves better resolution properties than WENO-Z and WENO-JS schemes, especially in the region near critical points.

4.2 Shock-tube problems

In this test case, we show the proposed regularized WENO-IS scheme passing the shock tube problems, which include the Sod problem [20], the Lax problem [21] and the 1-2-3 problem [22]. These computations are executed on a 200-point grid.

The initial condition for Sod problem is

[TABLE]

and the final time is $t=2$.

The initial condition for Lax problem is

[TABLE]

and the final time is $t=0.14$.

The initial condition for 123 problem is

[TABLE]

and the final time is $t=1$.

The comparison of the numerical results of density and velocity distribution is shown in Fig. 3. The ”Exact” denotes the theoretical solution. It is obvious that the present method performs as well as other WENO schemes.

4.3 Interaction blast waves

A two-blast-wave interaction problem  [23] is considered and the initial condition is

[TABLE]

The Roe flux is adopted for flux splitting and a reflection boundary condition is applied at $x=0$ and $x=1$. The computing is carried on up to $t=0.038$ on a 400-point grid.

Fig. 4 gives the comparison of density profile and a enlarged part. The ”Exact” denotes the convergent result which is computed by WENO-JS at the resolution of $N=2500$. It can be observed that all methods capture the strong shock waves very well. However, the original WENO-IS scheme show considerate larger numerical dissipation and the present method exhibits comparatively less dissipation than others, especially near the right peak of the density profile as shown in Fig. 4.

4.4 The Shu-osher problem

This case gives a Mach 3 shock wave interaction with a sine entropy wave and the initial condition is

[TABLE]

The computing is carried on a 200-point grid with the domain of $x\in\left[0,10\right]$ and the finial time is $t=1.8$. The density profiles computed by different schemes are shown in Fig. 5. The reference ”Exact” denotes a high-resolution result solved by WENO-JS on a 3200-point grid. It can be clearly noticed that all schemes can capture the shock waves well. However, in the region downstream the shock wave, the present method shows superior resolution than WENO-Z, WENO-JS and WENO-IS schemes in reproducing the fine flow structures. Also, compared with TENO5 and TENO5-opt scheme as shown in  [14](their Fig. 15), the density profile of the present method is more close to the convergent density profile especially in the region downstream the shock wave.

4.5 Double Mach reflection of a strong shock

This is a two-dimensional test case taken from Woodward and Colella [23] on the double Mach reflection of a strong shock. The initial condition is

[TABLE]

and the computation is solved on the domain of $\left[0,4\right]\times\left[0,1\right]$. This case represents a right moving shock wave with Mach 10 initially locates at $x=0.1667$, $y=0$ and makes an angle of $60^{\circ}$ with $x-$axis. For the lower boundary, the reflection wall condition is imposed for $0.1667<x<4$, while the post shock condition is applied for $0<x<0.1667$. The computation is carried on two grids with solutions of $512\times 128$ and $1024\times 256$ and the final time is $t=0.2$.

The density profiles on $512\times 128$ and $1024\times 256$ grids are shown in Fig. 6 and Fig. 7, respectively. The close-up view of the ”blow up” region is presented to compare the fine flow structure resolving abilities of different schemes. It can be noticed that all methods can capture the main flow features, such as the Mach stem and the near wall jet. However, the fine structures around the slip line and the near wall jet are considerably different in their solutions. Among these schemes, the present method resolves the finest flow structures and strongest near wall jet. Although WENO-Z scheme resolves nearly the same amount of roll-up structures with the present method, it is inferior in terms of the direction and the shape of the near wall jet compared with the results with very high resolutions, e.g. Fig.8 in Ref. [24], Fig.2.1 in Ref. [25] and Fig.20 in Ref. [7]. Table 1 gives the CPU time for the entire computing, which counts from the initialization to output, of double Mach reflection problem with different methods on the $512\times 128$ grid. We can find that the present method is computational efficient and only costs less than three fourth of WENO-Z.

4.6 Two-dimensional viscous shock-tube problem

This viscous shock tube problem is taken from Ref. [26] and has been investigated in Refs. [27, 28]. In this problem, the shock wave reflected from the end wall interacts with the boundary layer on the side wall induced by the incident shock. These interactions result in a complex system of vortices, shock wave bifurcation and other various flow structures. The flow filed containing various fine structures makes it a perfect case for testing high-order schemes. The initial condition is

[TABLE]

where $\gamma=1.4$ with ideal gas used. The non-slip adiabatic condition is applied to all boundary conditions of the tube. Since the configuration is symmetric about the line $y=0.5$, only half of the domain is computed, which is $\left[0,1\right]\times\left[0,0.5\right]$. Here, we set the Prandtl number $Pr=0.73$ and the final time is $1.0$. The viscosity is assumed to be constant and the Reynolds number of 200 and 1000 are considered.

The density profiles of the problem with the Reynolds number of 200 are shown in Fig. 8. The grid converged density profiles in Sjögreen and Yee [27](their Fig.2) and Zhou et al. [28](their Fig.3) obtained with much higher grid resolutions are taken as the reference for comparison. Considering the position of the triple point, the shape and the height of the primary vortex and the orientation of the long axis of the primary vortex, it is obvious that the present method gives considerable better result than other schemes. The comparison of the entire computational time of this case by various methods is shown in Table  2. We can see that the present method is more efficient than WENO-JS and WENO-Z scheme and only costs nearly three fourth of WENO-Z in this case.

Fig. 9 gives the density profile of the problem with the Reynolds number of 1000. The grid convergent result can be found in Ref. [28](their Fig.6) and Ref. [14](their Fig.21). With the Reynolds number increased to 1000, more fine flow structures appear in the flow field. As shown in Fig. 9, the present method performs well in terms of shock capturing and fine-structure resolution. We can see that the lambda-shape shocks captured by the present method and the WENO-Z scheme are more sharply than that of WENO-JS and WENO-HY. The big rotating structures at the lower right corner obtained by WENO-JS, WENO-HY, WENO-Z and TENO5 (see Fig.22 in Ref. [14]) are obviously different from the reference solution. However, the one produced by the present method fits the reference well. Furthermore, the primary vortex, the adjacent small vortices as well as the jet beneath the oblique shock are predicted quite well with the present method.

4.7 Mach 2000 jet flows

In order to verify the outstanding numerical stability of the present method, we consider the Mach 2000 jet problem, which has been computed in Zhang and Shu [16] with positive preserving limiters. The computing is carried out on the domain $\left[0,1\right]\times\left[0,0.25\right]$ with the resolution of $640\times 160$. The initial condition is $\left(\rho,u,v,p\right)=\left(0.5,0,0,0.4127\right)$ for the entire domain. A reflective condition is applied at the bottom boundary and an outflow condition is used at the right and top boundaries. For the left boundary, an inflow condition is adopted with states $\left(\rho,u,v,p\right)=\left(5,800,0,0.4127\right)$ when $y<0.05$ and $\left(\rho,u,v,p\right)=\left(0.5,0,0,0.4127\right)$ otherwise. We set $\gamma=5/3$ and the speed of the jet is 800, which gives about Mach 2100 with respect to the sound speed in the jet gas. In the computing, the CFL number is set to $0.25$ and the final time is $0.001$ as in the reference [15]. Note that, without additional positive preserving method applied, the WENO-JS, WENO-HY and WENO-Z schemes all blow out during the computing. However, WENO-IS and the present method perform well. The computed density and pressure profiles of WENO-IS and the present method are shown in Fig. 10. We can observe that these results are in good agreement with those in Zhang and Shu [16] (their Fig.4.6) and Hu et al. [15](their Fig.4). Furthermore, the present method shows better resolution than WENO-IS scheme and the classical WENO scheme combined with positive preserving flux limiters [16, 15] with more fine flow structures.

5 Conclusions

In this paper, a L2-norm regularized incremental-stencil WENO scheme is proposed for compressible flows. Following the stencil selection of WENO-IS scheme, we make the following modifications: (a) a L2-norm error term is introduced in the weighting strategy of WENO-IS scheme to cope with the discrepancy of the smoothness indicators between 2- and 3-point stencils; (b) a suitable non-dimensional regularization parameter is chosen to adaptively control the contribution of the error and regularized terms; (c) a hybridization with the optimal linear scheme is introduced to improve the performance of the scheme at smooth region and further improve computational efficiency. These modifications make the method can predict more fine flow structures by reducing numerical dissipation, and the remaining of 2-point stencils helps to compute the flow field with strong discontinuities. A set of test cases reveal that the proposed method has good fine-structure resolving capability and keeps the superior robustness.

Acknowledgements

The first author is partially supported by Xidian University (China) and the project of Natural Science Foundation of Shaanxi Province (Grant No:2019JM-186).