Fully discrete WENO with double entropy condition for hyperbolic conservation laws

ABSTRACT

This paper put forward a new fully discrete scheme construction method – double entropy condition solution formula method. With that, we turn the state-of-the-art semi-discrete WENO + RK scheme into a fully discrete scheme, which is named as Full-WENO. A major difficulty of this work is that we lack exact solution expressions for nonlinear equations in general cases. A feasible way we can go is to linearize equations and get quasi-exact solution formulas. The critical challenge is keeping both accuracy and efficiency in a scheme. Then, we get a class of new high-order schemes far better than traditional WENO schemes in the following aspects: (1) One-step to consistent high accuracy order in both space and time; (2) Resolution improves with the increasing CFL number; (3) Less CPU time and memory space, 1/s times of WENO with s-stage RK method in theory; (4) Excellent entropy condition satisfying property. Compared with  our original work , the new method applies the more sophisticated WENO reconstruction and solves the resolution loss problems in multi-dimensional cases. The numerical tests show that the new scheme is equipped with the merits of high efficiency, high resolution and low dissipation, especially for long-time nonlinear problems.

KEYWORDS:

• solution formula method
• fully discrete
• full-WENO
• Euler equations
• hyperbolic conservation laws
• unsteady compressible flow

1. Introduction

Aerodynamic designs of vehicles strongly depend on solutions of fluid equations, however, exact solutions can hardly be achieved due to nonlinearity of equations. At present stage, solutions are usually obtained in two ways: numerically through computational fluid dynamics (CFD) and experimentally through wind tunnels. Wind tunnel experiments are based on real fluids and thus the measurement results are credible, but also costly and difficult to obtain complete information of flow. While CFD can obtain more detailed flow information, there are no guarantees that the results are reliable for all flow conditions. The two main factors that contribute to the credibility of CFD are algorithms and turbulence, and this paper is concerned with the former.

From its earliest days, CFD has struggled with various flaws in its algorithms. One of the most serious flaws is the numerical oscillation of classical schemes caused by shocks. There have been ample research efforts aimed at dealing with this problem, with partial success, using techniques such as conservation forms, upwind difference, total variation diminishing (TVD) condition (Harten, 1983), limiters (Davis, 1987; Harten, 1983; Sweby, 1984), entropy condition for discontinuities (Lax, 1971), Riemann solvers (such as Roe, HLL, HLLC, etc.) (Toro, 2009), etc. With the pursuit of high-fidelity numerical solutions to meet engineering needs, a large number of distinctive high-order schemes were constructed with different ideas. By applying a first-order TVD scheme to equations with modified fluxes, Harten constructed a second-order TVD (Harten, 1983) scheme. By using limiters designed with the monotonicity preserving condition, VanLeer constructed the second-order- MUSCL scheme along the lines of the Godunov scheme, and now has been widely used in different fields (Sohn, 2005; Zhao et al., 2019). Since the TVD condition is too strict, it is not suitable for construction of higher-order schemes. Some researchers also constructed higher-order schemes which originated by Harten and developed by Osher and Harten (1987) and Shu (1997) with a more relaxed essentially non-oscillatory (ENO) condition. Several years later, a compact semi-discrete weighted essentially non-oscillatory (WENO) scheme came out associated with Runge-Kutta (RK) method which originated by LIU and Osher and improved by Jiang and Shu (1996). Soon, it became the mainstream in high-order schemes due to its excellent properties in shock-capturing. Under the impetus of Shu, etc., WENO schemes have been rapidly and widely used in many industries. Now, it has emerged a variety of WENO + RK schemes, such as MPWENO (Balsara & Shu, 2000), multistep WENO (Shen et al., 2014; Shen & Zha, 2014), multi-resolution WENO (Wang et al., 2021) and in addition with Lax-Wendroff type WENO schemes (LW-WENO) (Li & Du, 2016; Qiu, 2007; Zorío et al., 2017), ADER-WENO (Dumbser, 2014) schemes united with Arbitrary high-order DErivative Riemann (ADER) approaches (Titarev & Toro, 2002; Toro & Hidalgo, 2009), and Semi-Lagrangian/Eulerian-Lagrangian (SL/EL) WENO schemes (Huang et al., 2016; Huang & Arbogast, 2017).

Now we discuss the mentioned WENO schemes in detail according to their time discretization method: (1) WENO + RK schemes are semi-discrete, which use multi-stage RK method to enhance time accuracy order and avoid spurious oscillations. High-order RK methods not only occupy a heavy burden in both computing time and memory space, but also difficult to guarantee TVD and robustness properties of the scheme. Therefore, researchers usually apply the third-order TVD-RK (Jiang & Shu, 1996) method with a small CFL number. However, for some applications, such as numerical simulation of compressible turbulence and wave propagation problems involving long-time evolution it would be beneficial to use schemes which converge with higher order both in time and space. (2) LW-WENO schemes are fully discrete, which realize time discretization by replacing all time derivatives with space derivatives through LW procedure (Lax & Wendroff, 1960). Some can even improve their efficiency to twice of WENO + RK. Nevertheless, LW type schemes are extremely complex and not easy to implement in high-order situations, even need a multistep strategy (Li & Du, 2016). (3) ADER-WENO schemes are one-step and fully discrete, which realize time discretization by LW procedure and simplified generalized Riemann problem (GRP) solver. They decompose the difficult problem into a sequence of m conventional Riemann problems and finally achieve m-th accuracy order (it can be arbitrary m-th accuracy order in theory). However, this may be a little costly when applied to Euler equations. Titarev and Toro (2002) told that, under the same CFL number, ADER schemes can be 2–3 times faster than WENO + RK for linear equations with constant coefficients, but only 50% faster for 1D Euler equations; (4) SL/EL-WENO schemes are fully discrete, which realize time discretization by tracing characteristic lines. SL/EL-WENO can maintain its robustness in a relax CFL number, even nearly free from the limitation of CFL in cases of scalar nonlinear conservation laws (Huang et al., 2016). For scalar nonlinear cases, they apply RK method to trace the characteristic line, which leads to a very complex procedure compared with normal semi-discrete WENO + RK. And for 1D Euler equations, their method does not work very well due to there are three characteristic lines, even computational cost is tripled (Huang & Arbogast, 2017).

Different from aforementioned methods, this paper constructs a one-step fully discrete scheme through solution formula method (Dong et al., 2002). The main point of solution formula method is to construct a quasi-exact solution formula for the partial differential equation and then discretize it into a difference scheme. Since systems of conservation laws are nonlinear hyperbolic type and their solutions may contain discontinuities, it is not easy to construct exact solution formulas for general initial values. So, first we integrate them once in space and obtain the Hamilton-Jacobi (HJ) equations. Although HJ equations are still nonlinear (in flux), their solutions are continuous and thus easier for discretization. Via linearizing the flux of nonlinear HJ equations, we can obtain the quasi-exact solution formula. Applying Newton interpolation with limiters to the (quasi-) exact solution, Dong et al. (2002) constructed a second-order fully discrete entropy condition (EC) scheme, and then also develop it into a high-order version (Zhou & Dong, 2021, 2022). Since the solution formula of nonlinear equations containing discontinues is not easy to obtain, EC schemes inlay the entropy condition of discontinuities into the flux and thus get a quasi-exact solution after a flux linearization technique. This method can be generalized to systems, using discontinuity entropy conditions to get local quasi-exact solutions for systems (This method is easy to extend to systems, simply after applying the entropy condition of systems you can get the local quasi-exact solution formula). In this paper, we leverage the solution formula method and obtain a one-step fully discrete scheme by operating Newton interpolation with WENO construction, which is named as Full-WENO. Moreover, we find that it may result in some resolution loss or serve oscillation in multi-dimensional problems when distinguishing the wave evolution of Euler equations by the single entropy condition adopted in EC schemes. The reason is that it is only guided by velocity, however, the projected velocity in multi-dimensional problems may be identified as different properties in different directions when applying Strang split technique (Strang, 1968). So, we design a more accurate and reliable double entropy condition, which applies selective flux reconstruction according to both velocity and pressure.

Main framework of this paper is as follows: the second section introduces general construction steps of fully discrete schemes based on solution formula method; the third section proposes the concrete construction process of Full-WENO; the fourth section extends new scheme to three-dimensional Euler equations in curvilinear coordinates; the fifth section is numerical experiments, includes accuracy order tests, tests of scalar equation, one- and two-dimensional Euler equations, and sonic point test; the sixth section concludes.

2. Construction of fully discrete schemes via solution formula method

This section considers the construction of schemes for following one-dimensional Euler equations of perfect gas. (1) $\begin{array}{rl}& {\mathbf{\text{u}}}_t+\mathbf{\text{f}}(\mathbf{\text{u}}{)}_x=0\iff {\left(\begin{array}{l}\rho \\ \rho u\\ E\end{array}\right)}_t+{\left(\begin{array}{l}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right)}_x=0,\\ & p=(\gamma -1)(E-\frac{\rho u^2}{2}),c=\sqrt{\gamma \frac{p}{\rho }}.\end{array}$(1) Where ρ is density, u is velocity, p is pressure, E is total energy, γ is specific heat ratio, c is sound speed. The Jacobian matrix of flux A = f_u(u) can be written into the form containing only velocity u and sound speed c, namely A(u,c). This matrix has complete left/right eigenvectors and eigenvalues, which can also be rewritten into the function of u and c: L(u,c), R(u,c), λ(u,c).

The Euler equations are actually conservation laws (CL). First, we integrate the conservation laws once into HJ equations. And then, we construct solution formulas of HJ equations which will be used to obtain the numerical flux of conservative schemes. During this, there are two key techniques: (1) entropy condition linearization; (2) non-oscillation Newton interpolation methods. Introducing the space integration of conserved variables, we obtain HJ form of Euler equations. (2) $\begin{array}{rl}& \left(\text{CL}\right)\{\begin{array}{l}{\mathbf{\text{u}}}_t+\mathbf{\text{ f}}(\mathbf{\text{u}}{)}_x=0,\\ \mathbf{\text{u}}(x,t^n)={\mathbf{\text{u}}}^n\left(x\right),\end{array}\\ & \int \mathbf{u}dx\equiv \mathbf{v}\downarrow \uparrow \mathbf{u}\equiv {\mathbf{v}}_x\\ & \left(\text{HJ}\right)\{\begin{array}{l}{\mathbf{\text{v}}}_t+\mathbf{\text{f}}\left({\mathbf{\text{v}}}_x\right)=0,\\ \mathbf{\text{v}}(x,t^n)={\mathbf{\text{v}}}^n\left(x\right).\end{array}\end{array}$(2) Notice that the weak solution of CL systems is corresponding to the continuous solution of HJ systems, and it’s easier to get formulas of continuous solution. By using Newton interpolation, we can reconstruct solutions in continuous fields from discrete data (numerical solution). This method greatly simplifies the derivation of schemes, which is easier for readers to understand the essence of the algorithm. By the way, the solution formula in this paper all indicates to the solution of HJ equation in Equation (2). Note: ‘≡’ means ‘define’ all over this paper.

2.1. Construction of solution formula

Differential equation is an algebraic relation containing derivatives, while its solution is an algebraic relation containing no derivatives. Every part of scalar equation or matrix equations has same structure and is self-consistent in algebraic relationships, which means one can directly construct their solution formulas. However, Euler equations are vector systems, so we should use eigen-matrix to transform vectors into scalars and get solution formulas. There are five steps for constructing the solution formula of HJ systems: (1) scalarize (diagonalize) the system with eigenvector matrix; (2) linearize the flux via first-order Taylor expansion or linear interpolation; (3) construct scalarized solution formula; (4) compose solution formula for systems; (5) write general numerical flux by solution formula. Details of every step are as follows:

1. Premultiply the HJ systems by local constant left eigenvectors (3) $\begin{array}{rl}& \{\begin{array}{l}{\mathbf{\text{v}}}_t+\mathbf{\text{f}}\left({\mathbf{\text{v}}}_x\right)=0,\\ \mathbf{\text{v}}(x,t^n)={\mathbf{\text{ v}}}^n\left(x\right),\end{array}\\ & \downarrow \mathbf{f}\left(\mathbf{u}\right)⟹\text{diagonlization}{\phi }^k\left(\mathbf{u}\right)\equiv {\overline{\mathbf{L}}}^k\mathbf{f}\left(\mathbf{u}\right)\downarrow \\ & \{\begin{array}{l}{\mathbf{\text{L}}}^k{\mathbf{\text{v}}}_t+{\mathbf{\text{L}}}^k\mathbf{\text{f}}\left({\mathbf{\text{v}}}_x\right)=0,\\ {\mathbf{\text{L}}}^k\mathbf{\text{v}}(x,t^n)={\mathbf{\text{L}}}^k{\mathbf{\text{v}}}^n\left(x\right).\end{array}\end{array}$(3)

Where L^k is k-th row of matrix L, v and f are vectors, L^kv and L^kf are scalars, so Equation (3) shows the way from system (vector equations) to a group of scalar equations. Note: the initial value vⁿ(x) starts from n-th timestep tⁿ for the convenience of constructing numerical schemes. For one-dimensional Euler equations, k = 1, 2, 3, which means the k-th characteristic field.

1. Rewrite the flux into linear form to linearize the scalarized equations. (4) $\begin{array}{rl}& \{\begin{array}{l}{\mathbf{\text{L}}}^k{\mathbf{\text{v}}}_t+{\mathbf{\text{L}}}^k\mathbf{\text{f}}\left({\mathbf{\text{v}}}_x\right)=0,\\ {\mathbf{\text{L}}}^k\mathbf{\text{v}}(x,t^n)={\mathbf{\text{L}}}^k{\mathbf{\text{v}}}^n\left(x\right),\end{array}\stackrel{{\phi }^k\left(\mathbf{\text{u}}\right)\stackrel{\text{linearlization}}{⟶}{\lambda }^k{\mathbf{\text{L}}}^k\mathbf{\text{u}}-{\phi }^{\ast k}}{⟶}\\ & \{\begin{array}{l}{\mathbf{\text{L}}}^k{\mathbf{\text{v}}}_t+{\lambda }^k{\mathbf{\text{L}}}^k{\mathbf{\text{v}}}_x-{\phi }^{\ast k}=0,\\ {\mathbf{\text{L}}}^k\mathbf{\text{v}}(x,t^n)={\mathbf{\text{L}}}^k{\mathbf{\text{v}}}^n\left(x\right).\end{array}\end{array}$(4)

Where φ^k(u) is a nonlinear function, λ^kL^ku-φ*^k is a linear function, so Equation (3) shows the way from nonlinear equation to linear equation. Note: L^ku is a scalar variable, λ^k and φ*^k are constants.

1. The solution formula for the linearized scalarized equations can be given (5) $\begin{array}{rl}& \{\begin{array}{l}{\mathbf{\text{L}}}^k{\mathbf{\text{v}}}_t+{\lambda }^k{\mathbf{\text{L}}}^k{\mathbf{\text{v}}}_x-{\phi }^{\ast k}=0,\\ {\mathbf{\text{L}}}^k\mathbf{\text{v}}(x,t^n)={\mathbf{\text{L}}}^k{\mathbf{\text{v}}}^n\left(x\right).\end{array}\stackrel{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}}{⟶}\\ & \times {\mathbf{\text{L}}}^k{\mathbf{\text{v}}}^{n+1}\left(x\right)={\mathbf{\text{L}}}^k{\mathbf{\text{v}}}^n(x-{\lambda }^k\tau )+\tau {\phi }^{\ast k}.\end{array}$(5)

If we denote υ≡L^kv, a≡λ^k, f*≡φ*^k, the equation of Equation (5) can be written as υ_t+ aυ_x – f* = 0, this linear equation has an exact solution expression as υ(x,t^n + 1) = υⁿ(x-aτ)+τf*, just the same as the solution expression of Equation (5), where τ≡t^n + 1-tⁿ is time step size.

1. The solution formula for systems can be composed by right eigenvectors (6) $\begin{array}{rl}{\mathbf{\text{v}}}^{n+1}\left(x\right)& ={\mathbf{\text{R}}}^1\left({\mathbf{\text{L}}}^1{\mathbf{\text{v}}}^n\right(x-{\lambda }^1\tau )+\tau {\phi }^{\ast 1})\\ & +{\mathbf{\text{R}}}^2\left({\mathbf{\text{L}}}^2{\mathbf{\text{v}}}^n\right(x-{\lambda }^2\tau )+\tau {\phi }^{\ast 2})\\ & +{\mathbf{\text{R}}}^3\left({\mathbf{\text{L}}}^3{\mathbf{\text{v}}}^n\right(x-{\lambda }^3\tau )+\tau {\phi }^{\ast 3}).\end{array}$(6)

These four steps above can reach the solution formula expression of HJ systems. Where L^k, R^k, λ^k, φ^*k are all undetermined local constants, which will be given by discontinuity entropy condition in next section. Note: R^k is the k-th column of matrix R.

1. From solution formula Equation (6), we can obtain the general expression of numerical flux for conservative scheme (7) $\begin{array}{rl}{\mathbf{\text{u}}}_j^{n+1}& ={\mathbf{\text{u}}}_j^n-\frac{\tau }{h}({\widetilde{\mathbf{\text{f}}}}_{j+\frac{1}{2}}-{\widetilde{\mathbf{\text{f}}}}_{j-\frac{1}{2}})\stackrel{{\mathbf{\text{u}}}_j^n\equiv \frac{{\mathbf{\text{v}}}_{j+\frac{1}{2}}^n-{\mathbf{\text{v}}}_{j-\frac{1}{2}}^n}{h}}{⟵}\\ {\mathbf{\text{v}}}_{j+\frac{1}{2}}^{n+1}& ={\mathbf{\text{v}}}_{j+\frac{1}{2}}^n-\tau {\widetilde{\mathbf{\text{f}}}}_{j+\frac{1}{2}}\\ {\widetilde{\mathbf{\text{f}}}}_{j+\frac{1}{2}}& =\frac{{\mathbf{\text{v}}}_{j+\frac{1}{2}}^n-{\mathbf{\text{v}}}_{j+\frac{1}{2}}^{n\text{ + 1}}}{\tau },\\ {\mathbf{\text{v}}}_{j+\frac{1}{2}}^{n\text{ + 1}}& =\sum _{k=1}^3{\mathbf{\text{R}}}^k\left({\mathbf{\text{L}}}^k{\mathbf{\text{ v}}}^n\right(x_{j+\frac{1}{2}}-{\lambda }^k\tau )+\tau {\phi }^{\ast k}),\\ & \Updownarrow \\ {\widetilde{\mathbf{\text{f}}}}_{j+\frac{1}{2}}& =\sum _{k=1}^3{\mathbf{\text{R}}}^k({\lambda }^k{\mathbf{\text{L}}}^k{\overline{\mathbf{\text{u}}}}_{j+\frac{1}{2}}-{\phi }^{\ast k}),\\ {\mathbf{\text{ L}}}^k{\overline{\mathbf{\text{u}}}}_{j+\frac{1}{2}}& \equiv {\mathbf{\text{L}}}^k\frac{{\mathbf{\text{v}}}_{j+\frac{1}{2}}^n-{\mathbf{\text{v}}}^n(x_{j+\frac{1}{2}}-{\lambda }^k\tau )}{{\lambda }^k\tau }.\end{array}$(7)

Note: Equation (7) is just the quasi-exact solution in a numerical scheme form, it is not the final discrete scheme yet. ${\mathbf{\text{L}}}^k{\overline{\mathbf{\text{u}}}}_{j+\frac{1}{2}}$ in above equations can be simplified into scalar form due to L^k are local constants. And for each k the equations are all same, the k can be omitted. (8) $\begin{array}{rl}{\mathbf{L}}^k{\overline{\mathbf{u}}}_{j+{\displaystyle \frac{1}{2}}}={\overline{u}}_{j+\frac{1}{2}}& \equiv {\displaystyle \frac{{\upsilon }_{j+\frac{1}{2}}^n-{\upsilon }^n(x_{j+\frac{1}{2}}-\lambda \tau )}{\lambda \tau }}\\ \downarrow \nu & \equiv \frac{\lambda \tau }{h}\downarrow \\ {\overline{u}}_{j+\frac{1}{2}}& \equiv {\displaystyle \frac{{\upsilon }_{j+\frac{1}{2}}^n-{\upsilon }_{j+\frac{1}{2}-\nu }^n}{\nu h}}.\end{array}$(8) Where the values with subscript j+1/2+ν are generally not at the grid nodes, so we need an interpolation method to construct the initial value υⁿ(x). These are the whole idea of solution formula and general expression of numerical flux. Then we just need to solve two key techniques: (1) Discontinuity entropy condition linearization technique to decide local constants in solution formulas; (2) Non-oscillation Newton interpolation method to discrete numerical solution and to reconstruct the initial value υⁿ(x), and this paper we apply the thought of WENO reconstruction.

2.2. Double entropy condition linearization method

There exist global linearization methods for scalar nonlinear function, and thus exact solution formula can be constructed. However, the global linearization for nonlinear vector function cannot be achieved easily, so we try to construct quasi-exact formula by local linearization. Local linearization is adopted at the interface of grid cells [uⁿ_j, uⁿ_j+1], which is a discontinuity formed by bilateral values. By the way, the linearization method for nonlinear vector is not unique, now we apply following double entropy condition to determine the local constants (coefficients and constants of linear functions) in Equation (7):

u-entropy condition: We apply velocity entropy condition to obtain the local constants L, R, λ due to they are matrix as a whole that do not split with characteristic fields. Specifically, the local constants L^k, R^k, λ^k are constructed with the guide of bilateral velocity relationship u_L≷u_R at arbitrary discontinuities [u_L, u_R], so that we can scalarize the system.

(9) $\begin{array}{rl}u& \equiv \{\begin{array}{ll}\frac{\sqrt{{\rho }_L}{u}_L+\sqrt{{\rho }_R}{u}_R}{\sqrt{{\rho }_L}+\sqrt{{\rho }_R}}& u_L>u_R\\ \frac{{\rho }_L{u}_L+{\rho }_R{u}_R}{{\rho }_L+{\rho }_R}& u_L\le u_R\end{array},\\ H& \equiv \{\begin{array}{ll}\frac{\sqrt{{\rho }_L}{H}_L+\sqrt{{\rho }_R}{H}_R}{\sqrt{{\rho }_L}+\sqrt{{\rho }_R}}& u_L>u_R\\ H\left(\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}\right)& u_L\le u_R\end{array},\\ c& =\sqrt{(\gamma -1)(H-\frac{u^2}{2})}\to \{\begin{array}{l}\mathbf{\text{A}}\equiv \mathbf{\text{ A}}(u,c),\\ \mathbf{\text{L}}\equiv \mathbf{\text{L}}(u,c),\\ \mathbf{\text{R}}\equiv \mathbf{\text{R}}(u,c),\\ \left\{{\lambda }^k\right\}\equiv \mathit{\lambda }(u,c).\end{array}\end{array}$(9)

Concretely, we apply Roe means for u_L> u_R, while the characteristic lines tend to converge into a shock. And we apply arithmetic means for u_L≤ u_R, while the characteristic line tends to diverge into rarefaction waves. With this procedure, we can diagonalize the system. And here λ^k will be used in following steps.

p-entropy condition: We apply pressure entropy condition to obtain the local constants φ*^k due to they can be customized according to whether the characteristic fields are rarefaction waves or compression waves. Specifically, the local constants φ*^k are constructed with the guide of pressure relationship $p\gtrless p_{L/R}$ at arbitrary discontinuities [u_L, u_R], so that we can linearize the scalarized flux. Where p is an approximate middle pressure that can be obtained by discontinuity decomposition from bilateral status. (10) $\begin{array}{rl}p& ={(max\left(0,\frac{\frac{u_L-u_R}{2}+\frac{c_L}{\gamma -1}+\frac{c_R}{\gamma -1}}{\frac{c_L}{\gamma -1}{\left(\frac{1}{p_L}\right)}^{\frac{\gamma -1}{2\gamma }}+\frac{c_R}{\gamma -1}{\left(\frac{1}{p_R}\right)}^{\frac{\gamma -1}{2\gamma }}}\right))}^{\frac{2\gamma }{\gamma -1}},\\ {\phi }^{\ast 1}& \equiv {\mathbf{\text{ L}}}^1\left(\begin{array}{l}{\lambda }^1\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}-\frac{\mathbf{\text{f}}\left({\mathbf{\text{u}}}_L\right)+\mathbf{\text{f}}\left({\mathbf{\text{u}}}_R\right)}{2},p>p_L,\\ {\lambda }^1\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}-\mathbf{\text{ f}}\left(\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}\right),p\le p_L,\end{array}\right),\\ {\phi }^{\ast 2}& \equiv {\mathbf{\text{ L}}}^2\left({\lambda }^2\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}-\mathbf{\text{ f}}\left(\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}\right)\right),\\ {\phi }^{\ast 3}& \equiv {\mathbf{\text{L}}}^3\left(\begin{array}{l}{\lambda }^3\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}-\frac{\mathbf{\text{f}}\left({\mathbf{\text{u}}}_L\right)+\mathbf{\text{f}}\left({\mathbf{\text{u}}}_R\right)}{2},p>p_R,\\ {\lambda }^3\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}-\mathbf{\text{ f}}\left(\frac{{\mathbf{\text{u}}}_L+{\mathbf{\text{u}}}_R}{2}\right),p\le p_R,\end{array}\right).\end{array}$(10) Concretely, we apply Roe means for p > p_L/R, since these compression waves may evolve into shocks. Otherwise, we apply arithmetic means. This step is to linearize the scalarized equations. Note: Equations (7) and (8) and Equations (9) and (10) with [u_L, u_R] ≡ [uⁿ_j, uⁿ_j+1] make up the so-called Solution Formula Method of this paper.

Above techniques which determine the local constants by predicting the discontinuities evolve into shocks or rarefaction waves are named as entropy condition. Similar to the statement in (Zhou & Dong, 2021, 2022), Equations (9) and (10) that provides basic second order in temporal discretization for schemes in solution formula method can be named as baseline double entropy condition linearization. If we replace the arithmetic means (u_L+ u_R)/2 into a reconstructed high-order version ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$, it can provide higher-order nonlinear temporal accuracy order for scheme, details are shown in next section.

By the way, compared to the single entropy condition guided only by velocity applied in (Zhou & Dong, 2021, 2022), the double entropy condition guided by velocity and pressure in union will be more accurate, especially for multi-dimensional problems. The reason is the projected velocity may be identified as different properties (compression wave or rarefaction wave) in different directions during applying Strang split technique (Strang, 1968), especially at the slip lines (contact discontinuities in 1D problems). Concretely, the scheme may apply Roe mean for compression waves in one direction and apply reconstructed ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$ in the other direction for rarefaction waves, which may lead to a resolution loss or severe oscillation, such as that shown in section 5.4, Double Mach reflection problem. With double entropy condition, we can more precisely distinguish compression waves and rarefaction waves. Moreover, we customize the reconstruction proposal for three fields of Euler equations. Due to the second field of 1D Euler equations (second, fourth and fifth fields in 3D Euler equations) is equipped with linear features, we safely apply the high-order reconstructed ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$.

3. Construction of full-WENO via solution formula method

3.1. Initial value reconstruction

Using a Newton interpolation with WENO reconstruction (Jiang & Shu, 1996), we construct the initial value function L^kvⁿ(x) from initial data L^kv_iⁿ, i = j + $\frac{1}{2}$. And for writing convenience, we denote υⁿ(x)≡L^kvⁿ(x). (2r−1)-th Newton interpolation will occupy 2r−1 grids, consider the upwind condition and two numerical fluxes, the scheme will occupy 2r + 1 grids in total. According to WENO reconstruction, we weight several r-th stencils by performing with smoothness indicators, which not only helps to reach (2r−1)-th order at smooth regions but also can avoid the oscillation caused by high-order interpolation at risk regions. For example, a fifth-order stencil can be weighted by three third-order stencils, a seventh-order stencil can be weighted by four fourth-order stencils. Totally, this process will provide basic (2r−1)-th order in space for Full-WENO scheme. We give the specific expressions now. (Note: we give fifth and seventh schemes as examples here because it’s not easy to give the simple general expressions)

3.1.1. Initial value reconstruction for fifth order full-WENO (Full-WENO5)

The fully discrete interpolation and weighted coefficients of Full-WENO5 can be given as.

(11) $\begin{array}{rl}{\overline{u}}_{j+\frac{1}{2}}& =\{\begin{array}{ll}{u}_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+({\mathrm{\Delta }}_j^2+({\mathrm{\Delta }}_{j-\frac{1}{2}}^3+{\mathrm{\Delta }}_j^4\frac{3-\nu }{5})\frac{-1-\nu }{4})\frac{2-\nu }{3})\frac{1-\nu }{2},& a\ge 0,\\ u_{j+1}+({\mathrm{\Delta }}_{j+\frac{3}{2}}+({\mathrm{\Delta }}_{j+1}^2+({\mathrm{\Delta }}_{j+\frac{3}{2}}^3+{\mathrm{\Delta }}_{j+1}^4\frac{-3-\nu }{5})\frac{1-\nu }{4})\frac{-2-\nu }{3})\frac{-1-\nu }{2},& a\le 0,\end{array}\begin{array}{c}\stackrel{\text{weighted by 3 third order stencils}}{⟶}\\ \underset{\text{finial fifth order stencils}}{⟵}\end{array}\\ {\overline{u}}_{j+\frac{1}{2}}& =\{\begin{array}{cc}{\gamma }_1^{-}{u}_{j+\frac{1}{2}}^{\left(1\right)-}+{\gamma }_2^{-}{u}_{j+\frac{1}{2}}^{\left(2\right)-}+{\gamma }_2^{-}{u}_{j+\frac{1}{2}}^{\left(3\right)-},& a\ge 0,\\ {\gamma }_1^{+}{u}_{j+\frac{1}{2}}^{\left(1\right)+}+{\gamma }_2^{+}{u}_{j+\frac{1}{2}}^{\left(2\right)+}+{\gamma }_2^{+}{u}_{j+\frac{1}{2}}^{\left(3\right)+},& a\le 0.\end{array}\\ & \stackrel{a\ge 0}{\overbrace{\{\begin{array}{l}{u}_{j+\frac{1}{2}}^{\left(1\right)-}=u_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+{\mathrm{\Delta }}_{j-1}^2\frac{2-\nu }{3})\frac{1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(2\right)-}=u_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+{\mathrm{\Delta }}_j^2\frac{2-\nu }{3})\frac{1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(3\right)-}=u_j+({\mathrm{\Delta }}_{j+\frac{1}{2}}+{\mathrm{\Delta }}_{j+1}^2\frac{-1-\nu }{3})\frac{1-\nu }{2},\end{array}\{\begin{array}{l}{\gamma }_1^{-}=\frac{1}{20}(1+v)(2+v),\\ {\gamma }_2^{-}=\frac{1}{10}(3-v)(2+v),\\ {\gamma }_2^{-}=\frac{1}{20}(3-v)(2-v),\end{array}}}\\ & \stackrel{a\le 0}{\overbrace{\{\begin{array}{l}{u}_{j+\frac{1}{2}}^{\left(1\right)+}=u_{j+1}+({\mathrm{\Delta }}_{j+\frac{1}{2}}+{\mathrm{\Delta }}_j^2\frac{1-\nu }{3})\frac{-1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(2\right)+}=u_{j+1}+({\mathrm{\Delta }}_{j+\frac{3}{2}}+{\mathrm{\Delta }}_{j+1}^2\frac{-2-\nu }{3})\frac{-1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(3\right)+}=u_{j+1}+({\mathrm{\Delta }}_{j+\frac{3}{2}}+{\mathrm{\Delta }}_{j+2}^2\frac{-2-\nu }{3})\frac{-1-\nu }{2},\end{array}\{\begin{array}{c}{\gamma }_1^{+}=\frac{1}{20}(2+v)(3+v),\\ {\gamma }_2^{+}=\frac{1}{10}(2-v)(3+v),\\ {\gamma }_2^{+}=\frac{1}{20}(2-v)(1-v).\end{array}}}\end{array}$(11)

Where ν = aτ/h, a = λ^k, Δ^m is the m-th order difference, ${\mathrm{\Delta }}_{j+\frac{1}{2}}=u_{j+1}^n-u_j^n,{\mathrm{\Delta }}_j^2={\mathrm{\Delta }}_{j+\frac{1}{2}}-{\mathrm{\Delta }}_{j-\frac{1}{2}},{\mathrm{\Delta }}_{j+\frac{1}{2}}^3={\mathrm{\Delta }}_{j\text{ + 1}}^2-{\mathrm{\Delta }}_j^2$, etc. For nonlinear cases, we need to consider the accuracy of λ, and we achieve it through the flux reconstruction technique in section 3.2. The expression of ${\overline{u}}_{j+\frac{1}{2}}$ in Equation (11) is deduced via Newton interpolation of υⁿ(x), see Appendix A for detailed process.

Smoothness indicator for Full-WENO5, from (Jiang & Shu, 1996) (12) $\begin{array}{rl}& \stackrel{a\ge 0}{\overbrace{\{\begin{array}{l}{\text{β}}_1^{-}=\frac{13}{12}(u_{j-2}-2u_{j-1}+u_j{)}^2\\ +\frac{1}{4}(u_{j-2}-4u_{j-1}+3u_j{)}^2,\\ {\text{β}}_2^{-}=\frac{13}{12}(u_{j-1}-2u_j+u_{j+1}{)}^2\\ +\frac{1}{4}(u_{j-1}-u_{j+1}{)}^2,\\ {\text{β}}_3^{-}=\frac{13}{12}(u_j-2u_{j+1}+u_{j+2}{)}^2\\ +\frac{1}{4}(3u_j-4u_{j+1}+u_{j+2}{)}^2,\end{array}}}\\ & \stackrel{a\le 0}{\overbrace{\{\begin{array}{l}{\text{β}}_1^{+}=\frac{13}{12}(u_{j-1}-2u_j+u_{j+1}{)}^2\\ +\frac{1}{4}(u_{j-1}-4u_j+3u_{j+1}{)}^2,\\ {\text{β}}_2^{+}=\frac{13}{12}(u_j-2u_{j+1}+u_{j+2}{)}^2+\frac{1}{4}(u_j-u_{j+2}{)}^2,\\ {\text{β}}_3^{+}=\frac{13}{12}(u_{j+1}-2u_{j+2}+u_{j+3}{)}^2\\ +\frac{1}{4}(3u_{j+1}-4u_{j+2}+u_{j+3}{)}^2.\end{array}}}\end{array}$(12) Final interpolation of Full-WENO5 can be written as (13) $\begin{array}{r}{\overline{u}}_{j+\frac{1}{2}}=\{\begin{array}{c}\frac{\frac{{\gamma }_1^{-}}{{(ϵ+{\text{β}}_1^{-})}^2}}{\frac{{\gamma }_1^{-}}{{(ϵ+{\text{β}}_1^{-})}^2}+\frac{{\gamma }_2^{-}}{{(ϵ+{\text{β}}_2^{-})}^2}+\frac{{\gamma }_3^{-}}{{(ϵ+{\text{β}}_3^{-})}^2}}{u}_{j+\frac{1}{2}}^{\left(1\right)-}\\ +\frac{\frac{{\gamma }_2^{-}}{{(ϵ+{\text{β}}_2^{-})}^2}}{\frac{{\gamma }_1^{-}}{{(ϵ+{\text{β}}_1^{-})}^2}+\frac{{\gamma }_2^{-}}{{(ϵ+{\text{β}}_2^{-})}^2}+\frac{{\gamma }_3^{-}}{{(ϵ+{\text{β}}_3^{-})}^2}}{u}_{j+\frac{1}{2}}^{\left(2\right)-}\\ +\frac{\frac{{\gamma }_3^{-}}{{(ϵ+{\text{β}}_3^{-})}^2}}{\frac{{\gamma }_1^{-}}{{(ϵ+{\text{β}}_1^{-})}^2}+\frac{{\gamma }_2^{-}}{{(ϵ+{\text{β}}_2^{-})}^2}+\frac{{\gamma }_3^{-}}{{(ϵ+{\text{β}}_3^{-})}^2}}{u}_{j+\frac{1}{2}}^{\left(3\right)-},& a\ge 0,\\ \frac{\frac{{\gamma }_1^{+}}{{(ϵ+{\text{β}}_1^{+})}^2}}{\frac{{\gamma }_1^{+}}{{(ϵ+{\text{β}}_1^{+})}^2}+\frac{{\gamma }_2^{+}}{{(ϵ+{\text{β}}_2^{+})}^2}+\frac{{\gamma }_3^{+}}{{(ϵ+{\text{β}}_3^{+})}^2}}{u}_{j+\frac{1}{2}}^{\left(1\right)+}\\ +\frac{\frac{{\gamma }_2^{+}}{{(ϵ+{\text{β}}_2^{+})}^2}}{\frac{{\gamma }_1^{+}}{{(ϵ+{\text{β}}_1^{+})}^2}+\frac{{\gamma }_2^{+}}{{(ϵ+{\text{β}}_2^{+})}^2}+\frac{{\gamma }_3^{+}}{{(ϵ+{\text{β}}_3^{+})}^2}}{u}_{j+\frac{1}{2}}^{\left(2\right)+}\\ +\frac{\frac{{\gamma }_3^{+}}{{(ϵ+{\text{β}}_3^{+})}^2}}{\frac{{\gamma }_1^{+}}{{(ϵ+{\text{β}}_1^{+})}^2}+\frac{{\gamma }_2^{+}}{{(ϵ+{\text{β}}_2^{+})}^2}+\frac{{\gamma }_3^{+}}{{(ϵ+{\text{β}}_3^{+})}^2}}{u}_{j+\frac{1}{2}}^{\left(3\right)+},& a\le 0.\end{array}\end{array}$(13) Where ϵ is a positive real number introduced to avoid the denominator becoming zero, and we will take ϵ = 10⁻⁶ in later tests. Equations (12) and (13) come from reference (Jiang & Shu, 1996).

3.1.2. Initial value reconstruction for seventh order full-WENO (Full-WENO7)

The fully discrete interpolation and weighted coefficients of Full-WENO7 can be given as.

(14) $\begin{array}{rl}{\overline{u}}_{j+\frac{1}{2}}& =\{\begin{array}{ll}{u}_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+({\mathrm{\Delta }}_j^2+({\mathrm{\Delta }}_{j-\frac{1}{2}}^3+({\mathrm{\Delta }}_j^4+({\mathrm{\Delta }}_{j-\frac{1}{2}}^5+{\mathrm{\Delta }}_j^6\frac{4-\nu }{7})\frac{-2-\nu }{6}\left)\frac{3-\nu }{5}\right)\frac{-1-\nu }{4}\left)\frac{2-\nu }{3}\right)\frac{1-\nu }{2},& a\ge 0,\\ u_{j+1}+({\mathrm{\Delta }}_{j+\frac{3}{2}}+({\mathrm{\Delta }}_{j+1}^2+({\mathrm{\Delta }}_{j+\frac{3}{2}}^3+({\mathrm{\Delta }}_{j+1}^4+({\mathrm{\Delta }}_{j+\frac{3}{2}}^5+{\mathrm{\Delta }}_{j+1}^6\frac{-4-\nu }{7})\frac{2-\nu }{6}\left)\frac{-3-\nu }{5}\right)\frac{1-\nu }{4}\left)\frac{-2-\nu }{3}\right)\frac{-1-\nu }{2},& a\le 0,\end{array}\\ & \begin{array}{c}\stackrel{\text{weighted by 4 fourth order stencils}}{⟶}\\ \underset{\text{finial seventh order stencils}}{⟵}\end{array}\\ {\overline{u}}_{j+\frac{1}{2}}& =\{\begin{array}{cc}\sum _{k=1}^4{\gamma }_k^{-}{u}_{j+\frac{1}{2}}^{\left(k\right)-},& a\ge 0,\\ \sum _{k=1}^4{\gamma }_k^{+}{u}_{j+\frac{1}{2}}^{\left(k\right)+},& a\le 0.\end{array}\\ & \stackrel{a\ge 0}{\overbrace{\begin{array}{l}\{\begin{array}{l}{u}_{j+\frac{1}{2}}^{\left(1\right)-}=u_j^n+({\mathrm{\Delta }}_{j-\frac{1}{2}}+({\mathrm{\Delta }}_{j-1}^2+{\mathrm{\Delta }}_{j-\frac{3}{2}}^3\frac{3-\nu }{4}\left)\frac{2-\nu }{3}\right)\frac{1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(2\right)-}=u_j^n+({\mathrm{\Delta }}_{j-\frac{1}{2}}+({\mathrm{\Delta }}_j^2+{\mathrm{\Delta }}_{j-\frac{1}{2}}^3\frac{-1-\nu }{4}\left)\frac{2-\nu }{3}\right)\frac{1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(3\right)-}=u_j^n+({\mathrm{\Delta }}_{j-\frac{1}{2}}+({\mathrm{\Delta }}_j^2+{\mathrm{\Delta }}_{j+\frac{1}{2}}^3\frac{-1-\nu }{4}\left)\frac{2-\nu }{3}\right)\frac{1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(4\right)-}=u_j^n+({\mathrm{\Delta }}_{j+\frac{1}{2}}+({\mathrm{\Delta }}_{j+1}^2+{\mathrm{\Delta }}_{j+\frac{3}{2}}^3\frac{-2-\nu }{4}\left)\frac{-1-\nu }{3}\right)\frac{1-\nu }{2},\end{array}\\ \{\begin{array}{l}{\gamma }_1^{-}=\frac{1}{210}(1+v)(2+v)(3+v),\\ {\gamma }_2^{-}=\frac{1}{70}(4-v)(2+v)(3+v),\\ {\gamma }_3^{-}=\frac{1}{70}(4-v)(3-v)(3+v),\\ {\gamma }_4^{-}=\frac{1}{210}(4-v)(3-v)(2-v),\end{array}\end{array}}}\\ & \stackrel{a\le 0}{\overbrace{\begin{array}{l}\{\begin{array}{l}{u}_{j+\frac{1}{2}}^{\left(1\right)+}=u_{j+1}^n+({\mathrm{\Delta }}_{j+\frac{1}{2}}+({\mathrm{\Delta }}_j^2+{\mathrm{\Delta }}_{j-\frac{1}{2}}^3\frac{2-\nu }{4}\left)\frac{1-\nu }{3}\right)\frac{-1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(2\right)+}=u_{j+1}^n+({\mathrm{\Delta }}_{j+\frac{3}{2}}+({\mathrm{\Delta }}_{j+1}^2+{\mathrm{\Delta }}_{j+\frac{1}{2}}^3\frac{1-\nu }{4}\left)\frac{-2-\nu }{3}\right)\frac{-1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(3\right)+}=u_{j+1}^n+({\mathrm{\Delta }}_{j+\frac{3}{2}}+({\mathrm{\Delta }}_{j+1}^2+{\mathrm{\Delta }}_{j+\frac{3}{2}}^3\frac{1-\nu }{4}\left)\frac{-2-\nu }{3}\right)\frac{-1-\nu }{2},\\ u_{j+\frac{1}{2}}^{\left(4\right)+}=u_{j+1}^n+({\mathrm{\Delta }}_{j+\frac{3}{2}}+({\mathrm{\Delta }}_{j+2}^2+{\mathrm{\Delta }}_{j+\frac{5}{2}}^3\frac{-3-\nu }{4}\left)\frac{-2-\nu }{3}\right)\frac{-1-\nu }{2},\end{array}\\ \{\begin{array}{l}{\gamma }_1^{+}=\frac{1}{210}(2+v)(3+v)(4+v),\\ {\gamma }_2^{+}=\frac{1}{70}(3-v)(3+v)(4+v),\\ {\gamma }_3^{+}=\frac{1}{70}(3-v)(2-v)(4+v),\\ {\gamma }_4^{+}=\frac{1}{210}(3-v)(2-v)(1-v).\end{array}\end{array}}}\end{array}$(14)

The expression of ${\overline{u}}_{j+\frac{1}{2}}$ in Equation (14) is deduced via Newton interpolation of υⁿ(x), the detailed process is the same as fifth order case and is omitted here.

Smoothness indicator for Full-WENO7, from (Balsara & Shu, 2000) (15) $\begin{array}{rl}& \stackrel{a\ge 0}{\overbrace{\{\begin{array}{l}{\text{β}}_1=u_{j-3}(547u_{j-3}-3882u_{j-2}-4642u_{j-1}\\ -1854u_j)+u_{j-2}(7043u_{j-2}-1724u_{j-1}+7042u_j)\\ +u_{j-1}(11003u_{j-1}-9402u_j)+2107u_j^2,\\ {\text{β}}_2=u_{j-2}(267u_{j-2}-1642u_{j-1}+1602u_j-494u_{j+1})\\ +u_{j-1}(2843u_{j-1}-5966u_j+1922u_{j+1})\\ +u_j(3443u_j-2522u_{j+1})+547u_{j+1}^2,\\ {\text{β}}_3=u_{j-1}(547u_{j-1}-2522u_j+1922u_{j+1}-494u_{j+2})\\ +u_j(3443u_j-5966u_{j+1}+1602u_{j+2})\\ +u_{j+1}(2843u_{j+1}-1642u_{j+2})+267u_{j+2}^2,\\ {\text{β}}_4=u_j(2107u_j-9402u_{j+1}+7042u_{j+2}-1854u_{j+3})\\ +u_{j+1}(11003u_{j+1}-17246u_{j+2}+4642u_{j+3})\\ +u_{j+2}(7043u_{j+2}-3882u_{j+3})+547u_{j+3}^2,\end{array}}}\\ & \stackrel{a\le 0}{\overbrace{\{\begin{array}{l}{\text{β}}_1=u_{j+1}(2107u_{j+1}-9402u_j+7042u_{j-1}\\ -1854u_{j-2})+u_j(11003u_j-17246u_{j-1}\\ +4642u_{j-2})+u_{j-1}(7043u_{j-1}-3882u_{j-2})\\ +547u_{j-2}^2,\\ {\text{β}}_2=u_{j+2}(547u_{j+2}-2522u_{j+1}+1922u_j\\ -494u_{j-1}\left)+u_{j+1}\right(3443u_{j+1}-5966u_j+1602u_{j-1})\\ +u_j(2843u_j-1642u_{j-1})+267u_{j-1}^2,\\ {\text{β}}_3=u_{j+3}(267u_{j+3}-1642u_{j+2}+1602u_{j+1}\\ -494u_j)+u_{j+2}(2843u_{j+2}-5966u_{j+1}+1922u_j)\\ +u_{j+1}(3443u_{j+1}-2522u_j)+547u_j^2,\\ {\text{β}}_4=u_{j+4}(547u_{j+4}-3882u_{j+3}-4642u_{j+2}\\ -1854u_{j+1})+u_{j+3}(7043u_{j+3}\\ -1724u_{j+2}+7042u_{j+1})\\ +u_{j+2}(11003u_{j+2}-9402u_{j+1})+2107u_{j+1}^2.\end{array}}}\end{array}$(15) Final interpolation of Full-WENO7 can be written as (16) ${\overline{u}}_{j+\frac{1}{2}}=\{\begin{array}{cc}\sum _{k=1}^4\frac{\frac{{\gamma }_k^{-}}{{(ϵ+{\text{β}}_k^{-})}^2}}{\sum _{k=1}^4\frac{{\gamma }_k^{-}}{{(ϵ+{\text{β}}_k^{-})}^2}}{u}_{j+\frac{1}{2}}^{\left(k\right)-},& a\ge 0,\\ \sum _{k=1}^4\frac{\frac{{\gamma }_k^{+}}{{(ϵ+{\text{β}}_k^{+})}^2}}{\sum _{k=1}^4\frac{{\gamma }_k^{+}}{{(ϵ+{\text{β}}_k^{+})}^2}}{u}_{j+\frac{1}{2}}^{\left(k\right)+},& a\le 0.\end{array}$(16) Equations (15) and (16)come from reference (Balsara & Shu, 2000).

3.2. Flux reconstruction for full-WENO

This process is to reach designed accuracy for nonlinear equations, to be exact, we need to obtain the local constants L^k, R^k, λ^k, φ*^k for the composition of final numerical flux. We have given the baseline entropy condition linearization in section 2.2 and guarantee that the designed scheme is fully high order for linear equations. Now we show the fully high-order version in nonlinear cases. If we want high-order nonlinear accuracy, we need to linearize the flux with exact solution at time t^n + 1 for rarefaction waves and compression waves before shocks formed. However, we cannot get the nonlinear exact solution, so we need to reconstruct the quasi-exact solution ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$ by the data at time tⁿ to replace the arithmetic means in Equations (9) and (10). More specifically, the quasi-exact solution ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$ origins from solution formula method for conserved variables, which actually is the first-order derivative of ${\mathbf{\text{v}}}_{j+\frac{1}{2}}^{n\text{ + 1}}$ in Equation (7): (17) $\begin{array}{rl}{\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}& ={\mathbf{\text{R}}}^1\left({\mathbf{\text{L}}}^1{\mathbf{\text{v}}}_x^n\right(x_{j+\frac{1}{2}}-{\lambda }^1\tau \left)\right)+{\mathbf{\text{R}}}^2\left({\mathbf{\text{L}}}^2{\mathbf{\text{v}}}_x^n\right(x_{j+\frac{1}{2}}-{\lambda }^2\tau \left)\right)\\ & +{\mathbf{\text{R}}}^3\left({\mathbf{\text{L}}}^3{\mathbf{\text{v}}}_x^n\right(x_{j+\frac{1}{2}}-{\lambda }^3\tau \left)\right).\end{array}$(17) Here the initial L^k, R^k, λ^k in ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$ can be acquired by baseline u-entropy condition in Equation (9)

Then the reconstructed L^k, R^k, λ^k can be given by u-entropy condition and φ*^k can be given p entropy condition. (18) $\begin{array}{rl}u& \equiv \{\begin{array}{ll}\frac{\sqrt{{\rho }_j^n}{u}_j^n+\sqrt{{\rho }_{j+1}^n}{u}_{j+1}^n}{\sqrt{{\rho }_j^n}+\sqrt{{\rho }_{j+1}^n}}& u_j^n>u_{j+1}^n\\ u_{j+\frac{1}{2}}^{n+1}& u_j^n\le u_{j+1}^n\end{array},\\ H& \equiv \{\begin{array}{ll}\frac{\sqrt{{\rho }_j^n}{H}_j^n+\sqrt{{\rho }_{j+1}^n}{H}_{j+1}^n}{\sqrt{{\rho }_j^n}+\sqrt{{\rho }_{j+1}^n}}& u_j^n>u_{j+1}^n\\ H\left({\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}\right)& u_j^n\le u_{j+1}^n\end{array},\\ c& =\sqrt{(\gamma -1)(H-\frac{u^2}{2})}\to \{\begin{array}{l}\mathbf{\text{A}}\equiv \mathbf{\text{A}}(u,c),\\ \mathbf{\text{L}}\equiv \mathbf{\text{L}}(u,c),\\ \mathbf{\text{R}}\equiv \mathbf{\text{R}}(u,c),\\ \mathit{\lambda }\equiv \mathit{\lambda }(u,c),\end{array}\\ p& ={(max\left(0,\frac{\frac{u_j^n-u_{j+1}^n}{2}+\frac{c_j^n}{\gamma -1}+\frac{c_{j+1}^n}{\gamma -1}}{\frac{c_j^n}{\gamma -1}{\left(\frac{1}{p_j^n}\right)}^{\frac{\gamma -1}{2\gamma }}+\frac{c_{j+1}^n}{\gamma -1}{\left(\frac{1}{p_{j+1}^n}\right)}^{\frac{\gamma -1}{2\gamma }}}\right))}^{\frac{2\gamma }{\gamma -1}},\\ {\phi }^{\ast 1}& \equiv {\mathbf{\text{L}}}^1\left(\begin{array}{ll}{\lambda }^1\frac{{\mathbf{\text{u}}}_j^n+{\mathbf{\text{u}}}_{j+1}^n}{2}-\frac{\mathbf{\text{f}}\left({\mathbf{\text{u}}}_j^n\right)+\mathbf{\text{f}}\left({\mathbf{\text{u}}}_{j+1}^n\right)}{2}& p>p_j^n\\ {\lambda }^1{\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}-\mathbf{\text{f}}\left({\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}\right)& p\le p_j^n\end{array}\right),\\ {\phi }^{\ast 2}& \equiv {\mathbf{\text{L}}}^2\left({\lambda }^2{\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}-\mathbf{\text{f}}\left({\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}\right)\right),\\ {\phi }^{\ast 3}& \equiv {\mathbf{\text{L}}}^3\left(\begin{array}{ll}{\lambda }^3\frac{{\mathbf{\text{u}}}_j^n+{\mathbf{\text{u}}}_{j+1}^n}{2}-\frac{\mathbf{\text{f}}\left({\mathbf{\text{u}}}_j^n\right)+\mathbf{\text{f}}\left({\mathbf{\text{u}}}_{j+1}^n\right)}{2}& p>p_{j+1}^n\\ {\lambda }^3{\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}-\mathbf{\text{f}}\left({\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}\right)& p\le p_{j+1}^n\end{array}\right).\end{array}$(18) Due to the accuracy order of L^k, R^k, λ^k getting from baseline u-entropy condition is not equal to the reconstructed ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$, it will influence the accuracy to some content. So, theoretically, iteration is needed for getting a more accurate ${\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}$. However, we do not recommend this process owing to the complexity of Euler equations and consideration of computing burden. Actually, the numerical experiments tell that the error gently corrected by iteration contributes few influences in the vast majority of tests, except for testing nonlinear accuracy order. Moreover, the greatest influence factors on resolutions are the reconstructed λ^k φ*^k, so we even do not need to reconstruct L^k, R^k for almost all tests. It means we can directly use the L^k, R^k in Equation (17) which is calculated by baseline entropy condition. By the way, we also keep in line with the above statement in following numerical experiments. (Note: we can get the entropy condition linearization formula for the scalar equation if we remove the eigenvectors LR and use the eigenvalue for determining the selection between Roe means and reconstruction, see details in Zhou & Dong, 2021, 2022)

Denote ${\mathbf{\text{L}}}^k{\mathbf{\text{u}}}_{j+\frac{1}{2}}^{n+1}\equiv u_{j+\frac{1}{2}}^{\ast }$, according to the accuracy order analysis in (Zhou & Dong, 2021, 2022), final accuracy order of scheme via solution formula method = min (accuracy order of initial value reconstruction, 2 ×accuracy order of flux reconstruction). So, we can obtain designed order by r-th interpolation of u*_j+${}_{\frac{1}{2}}$. While doing r-th Newton interpolation for initial value υⁿ(x) in section 3.1, we can obtain their corresponding first-order derivative function u ⁿ(x). Inserting x = x${}_{j+\frac{1}{2}}$-νh, we get the final formula. And we use the smoothness indicators in section 3.1 to maintain the non-oscillation property. Details are shown as follows.

3.2.1. Flux reconstruction for r = 3, full-WENO5

(19) $\begin{array}{rl}& if({\text{β}}_k^2=min({\text{β}}_1^2,{\text{β}}_2^2,{\text{β}}_3^2\left)\right)\Rightarrow u_{j+\frac{1}{2}}^{\ast }=u_k^{\ast },\\ & \left[\begin{array}{l}{u}_1^{\ast }\\ u_2^{\ast }\\ u_3^{\ast }\end{array}\right]\equiv \{\begin{array}{ll}\left[\begin{array}{l}{u}_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}\\ +{\mathrm{\Delta }}_{j-1}^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}\\ +{\mathrm{\Delta }}_j^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}\\ u_j+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}\\ +{\mathrm{\Delta }}_{j+1}^2\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)}{3!}\end{array}\right],& a\ge 0,\\ \left[\begin{array}{l}{u}_{j+1}+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}\\ +{\mathrm{\Delta }}_j^2\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)}{3!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}\\ +{\mathrm{\Delta }}_{j+1}^2\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)}{3!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}\\ +{\mathrm{\Delta }}_{j+2}^2\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)}{3!}\end{array}\right],& a\le 0.\end{array}\end{array}$(19) The expression of $u_{j+\frac{1}{2}}^{\ast }$ in Equation (19) is deduced via Newton interpolation of υⁿ(x), see Appendix A for detailed process. The expansion form of Equation (19) is in Appendix A

3.2.2. Flux reconstruction for r = 4, full-WENO7

(20) $\begin{array}{rl}& if({\text{β}}_k^2=min({\text{β}}_1^2,{\text{β}}_2^2,{\text{β}}_3^2,{\text{β}}_4^2\left)\right)\Rightarrow u_{j+\frac{1}{2}}^{\ast }=u_k^{\ast },\\ & \left[\begin{array}{l}{u}_1^{\ast }\\ u_2^{\ast }\\ u_3^{\ast }\\ u_4^{\ast }\end{array}\right]\equiv \{\begin{array}{ll}\left[\begin{array}{l}{u}_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_{j-1}^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}+{\mathrm{\Delta }}_{j-\frac{3}{2}}^3\frac{d}{dv}\frac{-(-v)(1-v)(2-v)(3-v)}{4!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_j^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}+{\mathrm{\Delta }}_{j-\frac{1}{2}}^3\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)(2-v)}{4!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_j^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}+{\mathrm{\Delta }}_{j+\frac{1}{2}}^3\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)(2-v)}{4!}\\ u_j+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)}{3!}+{\mathrm{\Delta }}_{j+\frac{3}{2}}^3\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)(1-v)}{4!}\end{array}\right],& a\ge 0,\\ \left[\begin{array}{l}{u}_{j+1}+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}+{\mathrm{\Delta }}_j^2\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)}{3!}+{\mathrm{\Delta }}_{j-\frac{1}{2}}^3\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)(2-v)}{4!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)}{3!}+{\mathrm{\Delta }}_{j+\frac{1}{2}}^3\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)(1-v)}{4!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)}{3!}+{\mathrm{\Delta }}_{j+\frac{3}{2}}^3\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)(1-v)}{4!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{d}{dv}\frac{-(-1-v)(-v)}{2!}+{\mathrm{\Delta }}_{j+2}^2\frac{d}{dv}\frac{-(-2-v)(-1-v)(-v)}{3!}+{\mathrm{\Delta }}_{j+\frac{5}{2}}^3\frac{d}{dv}\frac{-(-3-v)(-2-v)(-1-v)(-v)}{4!}\end{array}\right],& a\le 0.\end{array}\end{array}$(20)

Where the smoothness indicators β_k in flux reconstruction can also be calculated by Equation (12) (for Full-WENO5) and Equation (15) (for Full-WENO7). With the minimum squared value of β_k, we can get the smoothest stencil for flux reconstruction and avoid the suspicious oscillations at risk regions as much as possible. Moreover, the β_k we get in this step can be directly applied in initial value reconstruction for saving computing time, which means the β_k applied in initial value reconstruction and flux reconstruction are same. The expansion form of Equation (20) is in Appendix A.

Additionally, for some cases with high demand in robustness, such as Test 12 Double Mach reflection in section 5.4, flux reconstruction requires an order reduction. We’d like to do as follows.

(21) $u_{j+\frac{1}{2}}^{\ast }=\{\begin{array}{lll}{u}_j,& if({\mathrm{\Delta }}_{j-\frac{1}{2}}{\mathrm{\Delta }}_{j+\frac{1}{2}}<0\cap {\mathrm{\Delta }}_{j-1}^2{\mathrm{\Delta }}_j^2<0\cap {\mathrm{\Delta }}_j^2{\mathrm{\Delta }}_{j+1}^2<0),& a\ge 0,\\ u_{j+1},& if({\mathrm{\Delta }}_{j+\frac{1}{2}}{\mathrm{\Delta }}_{j+\frac{3}{2}}<0\cap {\mathrm{\Delta }}_j^2{\mathrm{\Delta }}_{j+1}^2<0\cap {\mathrm{\Delta }}_{j+1}^2{\mathrm{\Delta }}_{j+2}^2<0),& a\le 0.\end{array}$(21)

Equation (21) means we do the order reduction while middle convexity is different with two ends of it (nonconvex-convex-nonconvex, convex-nonconvex-convex). This method balances robustness and resolution to some extent.

Note 1: The entropy condition of solution formula method and smoothness indicator of WENO guarantee high resolution and non-oscillation properties together. Then Full-WENO can be composed as follows.

Full-WENO5 is composed of Equations (7–13) and (18 and 19).

Full-WENO7 is composed of Equations (7–10) and (14–16) and (18 and 20).

Note 2: If we remove the flux reconstruction, take ν = 0 and associate it with RK method, the scheme in this paper will degrade to a normal semi-discrete WENO + RK scheme.

Note 3: If we remove the eigenvectors, solution formula, and entropy condition, then use RK method to trace the characteristic line, we can get the SL/EL-WENO (Huang et al., 2016) for scalar nonlinear equation. What we mean is Full-WENO may quite different from SL/EL-WENO (Huang et al., 2016; Huang & Arbogast, 2017), especially for systems.

Note 4: In theory, with same computational condition, the computing speed of Full-WENO will be s times faster than same order WENO with s-stage RK method. For example, for six-stage RK5 and nine-stage RK7 (Butcher, 2016) (see Appendix B), the computing speed of Full-WENO5 will be 6 times faster than WENO5 + RK5, and Full-WENO7 will be 9 faster than WENO7 + RK7.

4. Multi-dimensional cases

For three-dimensional Euler equations in curvilinear coordinates, we use dimensional splitting method (Strang, 1968) and turn it into three one-dimensional equations. (22) $\begin{array}{rl}& {\mathbf{\text{u}}}_t+\sum _{e=\xi ,\eta ,\zeta }{\mathbf{\text{f}}}^e{\left(\mathbf{\text{u}}\right)}_e=0\iff {\left(\begin{array}{l}\rho \\ \rho u\\ \rho v\\ \rho w\\ E\end{array}\right)}_t\\ & +\sum _{e=\xi ,\eta ,\zeta }{\left(\begin{array}{l}\rho (e_{x}u+e_{y}v+e_{z}w)\\ \rho (e_{x}u+e_{y}v+e_{z}w)u+e_{x}p\left(\mathbf{\text{u}}\right)\\ \rho (e_{x}u+e_{y}v+e_{z}w)v+e_{y}p\left(\mathbf{\text{u}}\right)\\ \rho (e_{x}u+e_{y}v+e_{z}w)w+e_{z}p\left(\mathbf{\text{u}}\right)\\ (e_{x}u+e_{y}v+e_{z}w)(E+p(\mathbf{\text{u}}\left)\right)\end{array}\right)}_e=0,\\ & \Downarrow \text{dimensional splitting}\\ & {\mathbf{\text{u}}}_t+{\mathbf{\text{f}}}^e(\mathbf{\text{u}}{)}_e=0\iff {\left(\begin{array}{l}\rho \\ \rho u\\ \rho v\\ \rho w\\ E\end{array}\right)}_t\\ & +{\left(\begin{array}{l}\rho (e_{x}u+e_{y}v+e_{z}w)\\ \rho (e_{x}u+e_{y}v+e_{z}w)u+e_{x}p\left(\mathbf{\text{u}}\right)\\ \rho (e_{x}u+e_{y}v+e_{z}w)v+e_{y}p\left(\mathbf{\text{u}}\right)\\ \rho (e_{x}u+e_{y}v+e_{z}w)w+e_{z}p\left(\mathbf{\text{u}}\right)\\ (e_{x}u+e_{y}v+e_{z}w)(E+p(\mathbf{\text{u}}\left)\right)\end{array}\right)}_e=0,\\ & (e=\xi ,\eta ,\zeta ).\end{array}$(22) Equation (1) and the last equation of Equation (22) is absolutely the same form except that the number of components is different. So it is easy to generalize all formulas of schemes (Full-WENO5 and Full-WENO7) for Equation (22) from Equation (1). For example, the general expression of numerical flux for conservative scheme in 3D can be written as follows (23) $\begin{array}{rl}{\mathbf{\text{u}}}_j^{n+1}& ={\mathbf{\text{u}}}_j^n-\frac{\tau }{h_e}({\widetilde{\mathbf{\text{f}}}}_{j+\frac{1}{2}}^e-{\widetilde{\mathbf{\text{f}}}}_{j\text{ - }\frac{1}{2}}^e),e=\xi ,\eta ,\zeta ,\\ k& =1,2,3,4,5,\\ {\widetilde{\mathbf{\text{f}}}}_{j+\frac{1}{2}}^{\text{e}}& =\sum _{k=1}^5{\mathbf{\text{R}}}^{e,k}({\lambda }^{e,k}{\mathbf{\text{L}}}^{e,k}{\overline{\mathbf{\text{u}}}}_{j+\frac{1}{2}}-{\phi }^{\ast e,k}),\\ {\mathbf{\text{L}}}^{e,k}{\overline{\mathbf{\text{u}}}}_{j+\frac{1}{2}}& \equiv {\mathbf{\text{L}}}^{e,k}\frac{{\mathbf{\text{v}}}_{j+\frac{1}{2}}^n-{\mathbf{\text{v}}}^n(e-{\lambda }^{e,k}\tau )}{{\lambda }^{e,k}\tau }.\end{array}$(23) More detailed information on three-dimensional cases can be found in our former work (Zhou & Dong, 2021, 2022).

5. Numerical experiments

This section we do numerical experiments for Full-WENO, which contains accuracy order tests, scalar equation tests, Euler equation tests, and sonic point test. For all tests, we compare our new schemes with frequently used WENO-RK schemes. If not specifically stated WENO-RK uses TVD-RK3 (Jiang & Shu, 1996) (see Appendix B) and Local Lax-Friedrichs (LLF) flux (Svenn Tveit, 2011). All numerical simulations are calculated by our own codes developed in Fortran95 with Microsoft Visual Studio 2013 and no specific data is used.

5.1. Tests of accuracy order

This subsection we test accuracy order of new scheme, which contains linear equation, Burgers equation, and point-by-point accuracy order tests.

5.1.1. Test 1 Accuracy order test, linear equation

This is an accuracy order test of linear equation, initial value is given as follows. (24) $\{\begin{array}{ll}{u}_t+u_x=0,& t=2,\\ u_0\left(x\right)=\sin \left(\pi x\right),& x\in [-1,1].\end{array}$(24)

5.1.2. Test 2 Accuracy order test, Burgers equation

This is an accuracy order test of nonlinear equation, initial value is given as follows. (25) $\{\begin{array}{ll}{u}_t+{\left(\frac{u^2}{2}\right)}_x=0,& t=0.15,\\ u_0\left(x\right)=0.5+\sin \left(\pi x\right),& x\in [0,2].\end{array}$(25) As shown in Figures 1 and 2, both Full-WENO and WENO + RK almost achieve designed accuracy order. Full-WENO reaches similar accuracy order compared with the same order WENO + RK, however, with smaller errors. While WENO + RK3 only shows about third order under CFL = 0.5 due to the space discrete does not dominate the accuracy order at this CFL number. In test 1 Full-WENO7 and WENO7 + RK7 only reach about sixth order and they will reach about eighth order if we continue to refine the grids. After further research, we found this phenomenon may cause by the positive real number ϵ in smoothness indicator, which has been discussed in reference (Henrick et al., 2005).

Figure 1. Accuracy order test, linear equation. (A) 5-th order schemes (B) 7-th order schemes.

Figure 2. Accuracy order test, Burgers equation. (A) 5-th order schemes (B) 7-th order schemes.

5.1.3. Test 3 Point-by-point accuracy order (PbP order)

Here we construct a new method for testing accuracy order. As the name suggests, point-by-point accuracy order is to reckon the accuracy order at every point by applying the theoretical error and accuracy. This method can visually locate positions of order loss and is helpful for researchers checking and analyzing bugs of algorithms. It’s meaningful for research of smoothness indicators. Consider following linear equation. (26) $\{\begin{array}{ll}{u}_t+a{u}_x=0,& t\le T,\\ a=1,& T=2.\end{array}$(26) point-by-point accuracy order can be given as (27) $q\left(x\right)=\ln \frac{R\left(x\right)}{C\left(x\right)}/\ln h.$(27) Where $C\left(x\right)=C{u}^{(q+1)}=C\frac{|u^{(q+1)}\left(x\right)|}{\overline{|u^{(q+1)}|}}$, and theoretical error R = Ch^q has been given in Table 1, C is theoretical error coefficient, q is theoretical accuracy order, |u^(q + 1)(x)| is q-th derivative of initial value, $\overline{|u^{(q+1)}|}$ is the average value of q-th derivative, h is space step. Here we test an infinitely smooth case, initial value can be found in Equation (24), and for this case we have $u^{(q+1)}=\{\begin{array}{l}\frac{{\pi }^{q+1}\cos \left(\pi x\right)}{2{\pi }^q},q=2k-1,\\ \frac{{\pi }^{q+1}\sin \left(\pi x\right)}{2{\pi }^q},q=2k,\end{array}k=1,2,3\dots$

Table 1. Accuracy order test, theoretical error of schemes.

	Theoretical error R (a = 1, T = 2, v = CFL)	CFL = 0.1	CFL = 0.5	CFL = 0.9	CFL = 1.0
Full-WENO5	$R_{\text{Full - WENO5}}=\left(aT\right|\frac{(1-{\nu }^2)(2^2-{\nu }^2)(3-\nu )}{6!}\left|2{\pi }^5\right)h^5$	19.475 h^5	11.954 h^5	2.164 h^5	0
Full-WENO7	$R_{\text{Full - WENO7}}=\left(aT\right|\frac{(1-{\nu }^2)(2^2-{\nu }^2)(3^2-{\nu }^2)(4-\nu )}{8!}\left|2{\pi }^7\right)h^7$	41.497 h^7	25.808 h^7	4.611 h^7	0
WENO5	$R_{\text{WENO5}}=\left(aT\right|\frac{2^2\times 3}{6!}\left|2{\pi }^5\right)h^5$			20.401 h^5
WENO7	$R_{\text{WENO7}}=\left(aT\right|\frac{2^2\times 3^2\times 4}{8!}\left|2{\pi }^7\right)h^7$			43.147 h^7

Table 1 shows that, for linear cases, the theoretical error of WENO + RK is not relevant to CFL number without the consideration of the error origin from RK method. When CFL→1 the theoretical error of Full-WENO tends to be zero, and when CFL→0 it tends to be same as semi-discrete WENO (for nonlinear cases Full-WENO needs to consider the error of ν = λτ/h). Figure 3 shows the point-by-point accuracy order, where optimal value is calculated directly by optimal weight without smoothness indicator. Figure 3 also tells that Full-WENO achieves similar accuracy order compared with same order WENO + RK, however, WENO + RK3 only achieves third to fourth order. And we can also find, it may introduce new errors after applying smoothness indicator and thus trigger an overall accuracy reduction, which means there is still room for improvement in this smoothness indicator. The accuracy order test method in Test 1 may mask this problem due to the division of two errors at different grids, and finally get designed accuracy.

Figure 3. Accuracy order test, point-by-point accuracy. (A) 5-th order schemes (B) 7-th order schemes.

5.2. Scalar equation

5.2.1. Test 4 Linear scalar equation with multiple extremes

This test is composited by a series of smooth and unsmooth functions, which contains a Gaussian, a square, a triangle, and a semi-ellipse. It is widely used to test the performance of discontinuity capturing and dissipation. Initial value is given as follows.

(28) $\{\begin{array}{l}{u}_0\left(x\right)=\{\begin{array}{ll}\frac{1}{6}\left(G\right(x,\text{β},z-\delta )+G(x,\text{β},z+\delta )+4G(x,\text{β},z\left)\right),& -0.8\le x\le -0.6,\\ 1,& -0.4\le x\le -0.2,\\ 1-|10(x-0.1)|,& 0\le x\le 0.2,\\ \frac{1}{6}\left(F\right(x,\alpha ,a-\delta )+F(x,\alpha ,a+\delta )+4F(x,\alpha ,a\left)\right),& 0.4\le x\le 0.6,\\ 0,& \text{else,}\end{array}\\ \begin{array}{l}G=e^{-\text{β}{(x-z)}^2},F(x,\alpha ,a)=\sqrt{max(1-{\alpha }^2{(x-a)}^2,0)},\\ z=-0.7,\delta =0.005,\text{β}=(\log 2)/\left(36{\delta }^2\right),a=0.5,\alpha =10.\end{array}\end{array}$(28)

The result of Figure 4 is solved with 200 points and t = 8, which shows the quite different solution properties between these two schemes. Combining Table 1, we can know the resolution of Full-WENO for this linear case tend to be similar with semi-discrete WENO + RK when CFL→0 and to be exact solution when CFL→1. On the contrary, due to RK3 method needs a small CFL number to maintain its robustness and accuracy, the resolution of WENO + RK3 performs worse with the increasing CFL number

Figure 4. Linear scalar equation, multiple extremes case. (A) Full-WENO5, overall view (B) Full-WENO5, local view (C) WENO5 + RK3, overall view (D) WENO5 + RK3, local view (E) Full-WENO7, overall view (F) Full-WENO7, local view (G) WENO7 + RK3, overall view (H) WENO7 + RK3, local view.

5.2.2. Test 5 Burgers equation

This is a nonlinear case, the initial value has been given in Equation (25). Forty points and period boundary are used, and we output the solution at t = 0.5 and t = 20 respectively. And we set CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3. Figure 5 shows that Full-WENO performs better at the discontinuity and still keeps a high resolution after a long-time evolution, while WENO + RK3 obviously loses its resolution. This test tells that Full-WENO is also equipped with an excellent low dissipative property in nonlinear cases, which is mainly attributed to the consistent high-order spatial and temporal accuracy.

Figure 5. Nonlinear scalar equation, Burgers equation. (A) 5-th order schemes, t = 0.5 (B) 7-th order schemes, t = 0.5 (A) 5-th order schemes, t = 20 (B) 7-th order schemes, t = 20.

5.3. One-dimensional Euler equations

5.3.1. Test 6 Sod shock tube problem

This is a typical Riemann problem for 1D Euler equations. Initial value is given as follows. (29) $(\rho ,u,p)=\{\begin{array}{ll}(1,0,1),& 0\le x\le 0.5,\\ (0.125,0,0.1),& 0.5\le x\le 1.\end{array}$(29) The solutions at t = 0.2 with 200 points, CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3 and reconstructed eigenvectors are plotted in Figure 6. It can be seen that, owing to the high order flux reconstruction for rarefaction waves and linear field, Full-WENO resolve the rarefaction waves and contract discontinuity much better than WENO + RK3. Also, owing to Roe means are applied for the compression waves that may evolve into shocks, the performance of Full-WENO at the shock is also better than WENO + RK3 with LLF flux.

Figure 6. 1D Euler equations, Sod shock tube problem. (A) 5-th order schemes, overall view (B) 5-th order schemes, local view (C) 7-th order schemes, overall view (D) 7-th order schemes, local view.

5.3.2. Test 7 Blast-wave problem

We now consider the interaction of two blast waves. Initial value is given as follows. (30) $(\rho ,u,p)=\{\begin{array}{ll}(1,0,{10}^3),& 0\le x<0.1,\\ (1,0,{10}^{-2}),& 0.1\le x<0.9,\\ (1,0,{10}^2),& 0.9\le x<1.\end{array}$(30) A reflective boundary condition is imposed at both ends. The simulation is performed with 500 points until final time t = 0.038. The reference solution is obtained by second-order TVD scheme (Harten, 1983) with 10000 points. It is shown in Figure 7 that Full-WENO resolves the density profile much better than same order WENO + RK3, especially near the valley x = 0.75 and the right peak x = 0.78, which means that Full-WENO exhibits less dissipation than WENO + RK3.

Figure 7. 1D Euler equations, Blast-Wave problem. (A) 5-th order schemes (B) 7-th order schemes.

5.3.3. Test 8 Shu-Osher problem

This is a typical example for testing the performance of high-order schemes when the solution contains both shocks and complex smooth region structures. In this case, a one-dimensional Mach 3 shock wave interacts with a perturbed density field generating both small-scale structures and discontinuities, hence it is selected to validate shock-capturing and wave-resolution capability. Initial value is given as follows. (31) $\begin{array}{r}(\rho ,u,p)=\{\begin{array}{ll}(3.857,2.629,10.333),& -5\le x<-4,\\ (1+0.2\sin (5x),0,1),& -4\le x<5.\end{array}\end{array}$(31) Here we solve the case with 200 and 400 grids respectively, and we set CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3. The reference solution is also obtained by second-order TVD scheme with 10,000 points. Figure 8 is the solution output at t = 1.8, which tells that Full-WENO performs much better than same order WENO + RK3 both at coarse and fine grids. This test well confirms that Full-WENO is less dissipation than WENO + RK3, which mainly attributes to its high-order accuracy in time evolution.

Figure 8. 1D Euler equations, Shu-Osher problem. (A) 5-th order schemes, 200 points (B) 7-th order schemes, 200 points (C) 5-th order schemes, 400 points (D) 7-th order schemes, 400 points.

5.4 Two-dimensional Euler equations

5.4.1. Test 9 Two-dimensional Riemann problem

This is a Riemann problem with interaction of planar shocks. Initial value is given as follows. (32) $\begin{array}{rl}& \left(p,u,v,p\right)\\ & =\{\begin{array}{ll}\left(0.138,1.206,1.206,0.029\right),& x<0.8,y<0.8,\\ \left(0.5323,0,1.206,0.3\right),& x>0.8,y<0.8,\\ \left(1.5,0,0,1.5\right),& x>0.8,y>0.8,\\ \left(0.5323,1.206,0,0.3\right),& x<0.8,y>\mathrm{0.8.}\end{array}\end{array}$(32) We solve the 2D Euler equations in a computational domain of (x, y)$\in$[0,1] × [0,1] with the mesh resolution 400 × 400 until final time t = 0.8. We set CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3. Figure 9 shows that Full-WENO resolves much better than same order WENO + RK3 in predicting the small-scale structures

Figure 9. 2D Euler equations, density distribution of 2D Riemann problem, 30 contours from 0.12 to 1.76. (A) Full-WENO5 (B) WENO5 (C) Full-WENO7 (D) WENO7.

5.4.2. Test 10 Shock vortex interaction problem

This model problem describes the interaction between a stationary shock and a vortex. At the initial moment, a stationary Mach 1.1 shock is positioned at x = 0.5 and normal to the x-axis. Its left state is $\left(p,u,v,p\right)=(1,\sqrt{\gamma },0,1)$, and its right state can be obtained through Hugoniot relation. A small vortex is superposed to the flow left to the shock and centers at (x_c, y_c) = (0.25, 0.5). The state of vortex can be described as a perturbation to velocity (u, v), temperature (T = p/ρ) and entropy (S = ln(p/ρ^γ)) of the mean flow, then we can denote it by the tilde values. (33) $\begin{array}{l}\tilde{u}=ϵ\tau e^{\alpha \left(1-{\tau }^2\right)}\sin \theta ,\\ \tilde{v}=-ϵ\tau e^{\alpha \left(1-{\tau }^2\right)}\cos \theta ,\\ \tilde{T}=-\frac{(\gamma -1){ϵ}^2{e}^{2\alpha \left(1-{\tau }^2\right)}}{4\alpha \gamma },\\ \tilde{S}=0,\\ \tau =r/r_c,r=\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}.\end{array}$(33) Where ϵ indicates vortex strength, α controls the decay rate of vortex, r_c is the critical radius to control the maximum vortex strength. Here we set ϵ = 0.3, r_c= 0.05, α = 0.204. The computational domain is taken as (x, y)$\in$[0, 2] × [0, 1], and the mesh of 200 × 80 is used. We set CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3. Figures 10 and 11 are the solution output at t = 0.6.

Figure 10. 2D Euler equations, shock vortex interaction, density at y = 0.5. (A) 5-th order schemes (B) 7-th order schemes.

Figure 11. 2D Euler equations, density distribution of shock vortex interaction, 30 contours from 1 to 1.22. (A) Full-WENO5 (B) WENO5 (C) Full-WENO7 (D) WENO7.

Figure 10 is the comparison of the density along the center line of y = 0.5, which shows Full-WENO has less numerical dissipation according to the performance of moving vortex. And Full-WENO also shows its excellent capturing capacity in Figures 10 and 11, for which the shock of Full-WENO is much sharper than WENO + RK3. It is mainly because the Roe means inlaid in flux linearization is friendly to shock. However, LLF flux applied in WENO + RK3 may introduce some dissipation. Moreover, if use Roe flux for WENO + RK3, some entropy fix may be needed to maintain its robustness and correctness which may also introduce some dissipation, see analysis in section 5.6.

5.4.3. Test 11 Rayleigh-Taylor instability problem

This problem illustrates the flow caused by the instability in domain (x, y)$\in$ [0, 0.25] × [0, 1]. Initial value is given as follows. (34) $\begin{array}{rl}& \left(\rho ,u,v,p\right)\\ & =\{\begin{array}{l}\left(2,0,-0.025a\cdot \cos \left(8\pi x\right),2y+1\right),0\le y<0.5,\\ \left(1,0,-0.025a\cdot \cos \left(8\pi x\right),y+1.5\right),0.5\le y\le 1.\end{array}\end{array}$(34) Where, $a=\sqrt{\gamma p/\rho }$ $\left(\gamma \text{ = }\frac{5}{3}\right)$ is the sound speed. The boundary conditions are set as follows: the initial condition at the bottom boundary is (ρ, u, v, p) = (2, 0, 0, 1); while at the top boundary (ρ, u, v, p) = (1, 0, 0, 2.5) is assigned; at the left and right boundaries, the reflecting boundary conditions are used. We solve the two-dimensional Euler equations by adding the source term (0, 0, ρ, ρv). Mesh 240 × 960 is used in this test, and we set CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3. We output the result at t = 1.95. As described in Figure 12, Full-WENO has resolved much richer vortical structures than same order WENO. This is because the uniform high order accuracy in time and space of Full-WENO will lead to less numerical dissipation. Moreover, a more relaxed CFL lead to fewer time steps for Full-WENO, which also helps to reduce the numerical dissipation.

Figure 12. 2D Euler equations, density distribution of Rayleigh-Taylor instability problem, 50 contours from 0.9 to 2.2. (A) Full-WENO5 (B) WENO5 (C) Full-WENO7 (D) WENO7.

5.4.4. Test 12 Double Mach reflection problem

This is a classic test for investigating high-resolution schemes. We solve the 2D Euler equations in a computational domain of (x, y)$\in$[0, 4] × [0, 1], and the mesh of 1600 × 400 is adopted. Here we also set CFL = 1 for Full-WENO, CFL = 0.4 for WENO + RK3. At the bottom of computational domain, it is a reflecting wall. Initially, a strong shock of Mach 10 makes a 60° angle with the bottom wall which begins at x = 1/6 and extends to the top of the domain at y = 0. Ahead of the shock, the initial condition is ρ = 1.4, p = 1.0, u = 0, v = 0, and the exact post-shock condition is imposed by Hugoniot relation. As shown in Figure 13, Full-WENO shows its merits of high resolution and low dissipation by resolving much more small structures compared with WENO + RK3.

Figure 13. 2D Euler equations, density distribution of Double Mach reflection problem, 30 contours from 2 to 22. (A) Full-WENO5 (double entropy condition) (B) Full-WENO5 (single entropy condition) (C) WENO5 (D) Full-WENO7 (double entropy condition) (E) Full-WENO7 (single entropy condition) (F) WENO7.

We also provide the original single entropy condition which is guided only by velocity. In Figure 13, we can find the resolution of double entropy condition improves much compared to that of single entropy condition (the reason has been explained in section 2.2).

5.5. CPU time test

This section we test the CPU time (in sec.) of 1D and 2D tests respectively. All tests are obtained by same workstation (AMD Ryzen Threadripper 3790X 32-Core Processor, 3.69 GHz) and same working environment.

5.5.1. Test 13 CPU time test, one-dimensional Euler equations

Here we compare the efficiency of fully discrete Full-WENO and semi-discrete WENO + RK. Here we choose the Sod shock tube, initial value sees Equation (29). In theory, the computing speed of Full-WENO should be s/CFL_semi/fully (CFL_semi/fully = (CFL of WENO + RK)/(CFL of Full-WENO)) times faster than same order WENO with s-stage RK method. As described in Figure 14, for 3-stage RK3, 6 stage-RK5, and 9-stage RK7, the CPU time tests all basically meet this law. Taking RK3 as an example: (1) With same CFL number, the computing speed of Full-WENO achieve about 2.8 times compared with same order WENO + RK3; (2) Considering WENO + RK3 requires a small CFL number for its resolution, with different CFL numbers, the computing speed of Full-WENO (CFL = 1) achieves about 6.8 times compared with WENO + RK3 (CFL = 0.4).

Figure 14. CPU time test, one-dimensional Euler equations, Sod shock tube problem. (A) 5-th order schemes (B) 7-th order schemes.

More specifically, here we extract a data from Figure 14 (800 points) to conclude the computing speed comparison in same order. Figure 15 tells that compared to the same time-space accuracy WENO + RK, Full-WENO is extremely equipped with competition in efficiency. The computing speed of Full-WENO5 can reach about 5.5 times as much as WENO5 + RK5 and Full-WENO7 can reach about 8 times as much as WENO7 + RK7 under same computing conditions, which also basically conforms to the theoretical speed (namely 6 and 9 times, respectively). Moreover, the storage memory of RK method is also more than s times compared to Full-WENO.

Figure 15. CPU time test, computing speed comparison in same order.

5.5.2. Test 14 CPU time test, two-dimensional Euler equations

The CPU time test of all 2D cases has been given in the corresponding test under section 5.4 (in sec). From Tables 2–5, we can obviously find the computing speed of Full-WENO achieve about 2.8 times as much as WENO + RK3 under CFL = 1, and achieve about 6.8 times if WENO + RK3 uses CFL = 0.4 with consideration of robustness and resolution. Because the entropy condition applies different strategies for compression waves and rarefaction waves (less computing burden for compression waves), some differences may happen in CPU time rate. (Note: We do not consider the CPU time of source terms in Rayleigh-Taylor instability problem)

Table 2. CPU time test, Two-dimensional Riemann problem.

Two-dimensional Riemann problem, grid number = 400 × 400, t = 0.8, 28 threads
scheme	CFL	CPU time (Rate)	CFL	CPU time (Rate)
r = 3
Full-WENO5	0.4	74.7 (2.44)	1.0	30.6 (1.00)
WENO5 + RK3	0.4	206.1 (6.74)	1.0	84.0 (2.75)
r = 4
Full-WENO7	0.4	88.4 (2.44)	1.0	36.2 (1.00)
WENO7 + RK3	0.4	244.1 (6.75)	1.0	100.8 (2.79)

Table 3. CPU time test, Shock vortex interaction.

Shock vortex interaction, grid number = 200 × 80, t = 0.6, 1 thread
scheme	CFL	CPU time (Rate)	CFL	CPU time (Rate)
r = 3
Full-WENO5	0.4	8.7 (2.23)	1.0	3.9 (1.00)
WENO5 + RK3	0.4	24.9 (6.38)	1.0	11.2 (2.87)
r = 4
Full-WENO7	0.4	11.2 (2.15)	1.0	5.2 (1.00)
WENO7 + RK3	0.4	32.8 (6.31)	1.0	15.0 (2.92)

Table 4. CPU time test, Rayleigh-Taylor instability problem.

Rayleigh-Taylor instability problem, grid number = 240 × 960, t = 1.95, 28 threads
scheme	CFL	CPU time (Rate)	CFL	CPU time (Rate)
r = 3
Full-WENO5	0.4	279.1 (2.42)	1.0	115.1 (1.00)
WENO5 + RK3	0.4	792.1 (6.88)	1.0	330.2 (2.87)
r = 4
Full-WENO7	0.4	335.1 (2.41)	1.0	139.1 (1.00)
WENO7 + RK3	0.4	935.0 (6.72)	1.0	388.6 (2.79)

Table 5. CPU time test, Double Mach reflection problem.

Double Mach reflection problem, grid number = 1600 × 400, t = 0.2, 28 threads
scheme	CFL	CPU time (Rate)	CFL	CPU time (Rate)
r = 3
Full-WENO5	0.4	297.4 (2.42)	1.0	123.9 (1.00)
WENO5 + RK3	0.4	846.2 (7.07)	1.0	356.8 (2.88)
r = 4
Full-WENO7	0.4	383.0 (2.40)	1.0	159.6 (1.00)
WENO7 + RK3	0.4	1090.1 (6.83)	1.0	448.5 (2.81)

5.6. Sonic point test

For schemes do not satisfy the entropy condition may render a non-physical solution at the rarefaction wave with a sonic point, such as WENO with Roe flux (just like that Roe scheme does not satisfy entropy condition). However, Full-WENO can efficiently avoid the non-physical solution by formula solution method constructed with entropy condition linearization.

5.6.1. Test 15 Sonic point test, reversed shock problem

We design a reversed shock problem that contains a sonic point at rarefaction. Five hundred points are used, initial value is given as follows. (35) $(\rho ,u,p)=\{\begin{array}{ll}(5.4,\text{7/9,}\text{31/3}),& -0.5\le x\le 0,\\ (1.4,3.0,1.0),& 0\le x\le \mathrm{1.5.}\end{array}$(35) In this test, we output the result at t = 0.3 and use the fifth-order scheme as an example. WENO5 + RK3 with Roe flux is set as the comparison, and is imposed with an entropy fix Q(λ) (Harten, 1983) ( = an extra artificial viscosity at the sonic point) (36) $Q\left(\lambda \right)\equiv \{\begin{array}{cc}|\lambda |,& |\lambda |>ϵ,\\ \frac{{\lambda }^2+{ϵ}^2}{2ϵ},& |\lambda |\le ϵ,\end{array}0\le ϵ\le 1.$(36) As shown in Figure 16, WENO5-Roe has a large unphysical discontinuity and tends to be normal only if an extra artificial viscosity is used. While Full-WENO5 has well resolved the solution owing to its entropy condition solution formula without any artificial viscosity.

Figure 16. Sonic point test, reversed shock problem.

5.7. Summary for numerical experiments

Overall, the numerical tests show that Full-WENO is equipped with the merits of high efficiency including high accuracy order, high resolution to discontinuities, and high speed of computation, when compared to traditional WENO-RK. The accuracy order tests verify that new schemes reach the designed accuracy order. Scalar tests including multiple extremes test and nonlinear test show that new scheme’s resolution to discontinuities and extremes is obviously higher than WENO-RK. One-dimensional Euler equation tests including Sod shock tube, Blast-wave, and Shu-Osher, prove that the solution formula method is feasible in fluid flow simulations, and keeps the priority to WENO-RK. Two-dimensional Euler equation tests including 2D-Reimann, Shock-vortex, RT-instabilities, and Double Mach reflection show more cases of high-quality results of new schemes using double entropy conditions, such as capturing much more fine structures than WENO + RK. CPU time test shows that new scheme has a speed of computation many times faster than WENO-RK. Sonic point test shows that new scheme can effectively avoid the unphysical discontinuity without any extra artificial viscosity, thanks to entropy condition inlaid in flux linearization.

6. Conclusions

The solution formula method is a construction method for difference schemes which we are in favor of. It is equipped with following merits: clear philosophy, easy derivation, natural upwind, and consistent space-time accuracy achievable. In this paper, we construct a Full-WENO scheme by solution formula method combined with some techniques of WENO reconstruction. The new scheme has following five highlights: (1) One-step and fully discrete: Full-WENO gets rid of RK method and achieves consistent high accuracy order both in space and time. And Full-WENO can perform robustly under CFL = 1; (2) The excellent performance along with CFL number: When CFL→1, Full-WENO not only tends to be exact in linear cases, but also to be more accurate in nonlinear cases. This is named the traveling wave solution property. However, all semi-discrete schemes propelled with RK method do not have this merit, which indicates that Full-WENO can reserve more accurate information of solutions. (3) High resolution: Full-WENO resolves better than same order WENO + RK for all CFL numbers. More importantly, the larger the CFL number, the more obvious the performance gap will be; (4) High computing speed: The computing speed of Full-WENO can reach about s times as much as semi-discrete WENO + s-stage RK under same computational condition. Specifically, Full-WENO5 is nearly 5∼6 times faster than WENO5 + RK5, Full-WENO7 is nearly 8∼9 times faster than WENO7 + RK7; (5) Double entropy condition linearization: Full-WENO uses the newly designed double entropy condition linearization guided by both velocity and pressure. Compared with single entropy condition only guided by velocity, the new method is more accurate and reliable. It remedies the accuracy loss caused by rough prediction of previous method (especially for multi-dimensional problems using Strang split). Moreover, compared with WENO + RK which needs entropy fix for Roe flux (or using LLF flux which has more numerical viscosity) to avoid non-physical solutions, Full-WENO embedded with entropy condition automatically avoids non-physical solutions without extra artificial numerical viscosity.

The solution formula method is valuable for engineering applications, particularly for long-time evolution problems. The entropy condition linearization is a critical technique in solution formula method which largely determines the resolution and robustness in nonlinear situations. However, there are still some flaws in current work: (1) How to remove the influence of numerical disturbances on the entropy condition selection which may result in a resolution loss due to incorrectly choose the Roe means; (2) How to exactly predict when the compression waves evolve into shocks (we can also apply the high-order flux reconstruction for compression waves before they evolve into shocks), etc. By solving these, the quality and robustness of numerical solutions may be lifted to a new level. Therefore, in the future, we plan to optimize the solution formula method from these two aspects. Moreover, we also hope to develop this into a more efficient large time step scheme, which has been simply tried in (Zhou & Dong, 2022). If realized, it may greatly improve engineering efficiency.

Note

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendices

Appendix A: Discretization process of Newton Interpolation in Full-WENO

We give the detailed Newton interpolation process for both initial value reconstruction and flux reconstruction of solution formula method. Consider Full-WENO5, for initial value reconstruction in case a ≥ 0, there are three third-order Newton interpolations of υ at [x${}_{j-\frac{1}{2}}$, x${}_{j+\frac{1}{2}}$] with six data around it, the middle one is as follows (A1) $\begin{array}{rl}\upsilon \left(x\right)={\upsilon }_{j+\frac{1}{2}}^n& +\frac{{\upsilon }_{j+\frac{1}{2}}-{\upsilon }_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}}(x-x_{j+\frac{1}{2}})+\frac{\frac{{\upsilon }_{j+\frac{3}{2}}-{\upsilon }_{j+\frac{1}{2}}}{x_{j+\frac{3}{2}}-x_{j+\frac{1}{2}}}-\frac{{\upsilon }_{j+\frac{1}{2}}-{\upsilon }_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}}}{x_{j+\frac{3}{2}}-x_{j-\frac{1}{2}}}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\\ & +\frac{\frac{\frac{{\upsilon }_{j+\frac{3}{2}}-{\upsilon }_{j+\frac{1}{2}}}{x_{j+\frac{3}{2}}-x_{j+\frac{1}{2}}}-\frac{{\upsilon }_{j+\frac{1}{2}}-{\upsilon }_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}}}{x_{j+\frac{3}{2}}-x_{j-\frac{1}{2}}}-\frac{\frac{{\upsilon }_{j+\frac{1}{2}}-{\upsilon }_{j-\frac{1}{2}}}{x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}}-\frac{{\upsilon }_{j-\frac{1}{2}}-{\upsilon }_{j-\frac{3}{2}}}{x_{j-\frac{1}{2}}-x_{j-\frac{3}{2}}}}{x_{j+\frac{1}{2}}-x_{j-\frac{3}{2}}}}{x_{j+\frac{3}{2}}-x_{j-\frac{3}{2}}}\times (x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\\ ={\upsilon }_{j+\frac{1}{2}}^n& +u_j(x-x_{j+\frac{1}{2}})+\frac{u_{j+1}-u_j}{2h}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})+\frac{\frac{u_{j+1}-u_j}{2h}-\frac{u_j-u_{j-1}}{2h}}{3h}(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\\ ={\upsilon }_{j+\frac{1}{2}}^n& +u_j(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\end{array}$(A1) We write the three Newton interpolations in square brackets and deduce the expression of ${\overline{u}}_{j+\frac{1}{2}}$ as follows (A2) $\begin{array}{rl}\upsilon \left(x\right)& =\left[\begin{array}{c}{\upsilon }_{j+\frac{1}{2}}^n+u_j(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-1}^2}{3!h^2}(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\\ {\upsilon }_{j+\frac{1}{2}}^n+u_j(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\\ {\upsilon }_{j+\frac{1}{2}}^n+u_j(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}}{2h}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j+1}^2}{3!h^2}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})(x-x_{j+\frac{3}{2}})\end{array}\right],\\ x_{j+\frac{1}{2}}& =(j+\frac{1}{2})h,x=(j+\frac{1}{2}-\nu )h\Rightarrow x-x_{j-\frac{3}{2}}=(2-\nu )h,x-x_{j-\frac{1}{2}}=(1-\nu )h,x-x_{j+\frac{1}{2}}=(-\nu )h,\\ & x-x_{j+\frac{3}{2}}=(-1-\nu )h,\\ {\upsilon }_{j+\frac{1}{2}}^{n+1}& \equiv \upsilon \left(\right(j+\frac{1}{2}-\nu \left)h\right)=\left[\begin{array}{l}{\upsilon }_{j+\frac{1}{2}}^n+u_j(-\nu )h+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(1-\nu )(-\nu )h^2+\frac{{\mathrm{\Delta }}_{j-1}^2}{3!h^2}(2-\nu )(1-\nu )(-\nu )h^3\\ {\upsilon }_{j+\frac{1}{2}}^n+u_j(-\nu )h+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(1-\nu )(-\nu )h^2+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}(2-\nu )(1-\nu )(-\nu )h^3\\ {\upsilon }_{j+\frac{1}{2}}^n+u_j(-\nu )h+\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}}{2h}(1-\nu )(-\nu )h^2+\frac{{\mathrm{\Delta }}_{j+1}^2}{3!h^2}(1-\nu )(-\nu )(-1-\nu )h^3\end{array}\right],\\ {\overline{u}}_{j+\frac{1}{2}}& \equiv \frac{{\upsilon }_{j+\frac{1}{2}}^n-{\upsilon }_{j+\frac{1}{2}}^{n+1}}{\nu h}=\left[\begin{array}{l}{u}_j+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2}(1-\nu )+\frac{{\mathrm{\Delta }}_{j-1}^2}{3!}(2-\nu )(1-\nu )\\ u_j+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2}(1-\nu )+\frac{{\mathrm{\Delta }}_j^2}{3!}(2-\nu )(1-\nu )\\ u_j+\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}}{2}(1-\nu )+\frac{{\mathrm{\Delta }}_{j+1}^2}{3!}(1-\nu )(-1-\nu )\end{array}\right]=\left[\begin{array}{l}{u}_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+{\mathrm{\Delta }}_{j-1}^2\frac{2-\nu }{3})\frac{1-\nu }{2}\\ u_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+{\mathrm{\Delta }}_j^2\frac{2-\nu }{3})\frac{1-\nu }{2}\\ u_j+({\mathrm{\Delta }}_{j+\frac{1}{2}}+{\mathrm{\Delta }}_{j+1}^2\frac{-1-\nu }{3})\frac{1-\nu }{2}\end{array}\right].\end{array}$(A2) For a ≥ 0, there is one fifth-order Newton interpolation at [x${}_{j-\frac{1}{2}}$, x_j${}_{j+\frac{1}{2}}$] with six data around it (A3) $\begin{array}{l}\upsilon \left(x\right)={\upsilon }_{j+\frac{1}{2}}^n+u_j(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\\ +\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}^3}{4!h^3}(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})(x-x_{j+\frac{3}{2}})+\frac{{\mathrm{\Delta }}_j^4}{5!h^4}(x-x_{j-\frac{5}{2}})(x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})(x-x_{j+\frac{3}{2}}),\\ x_{j+\frac{1}{2}}=(j+\frac{1}{2})h,x=(j+\frac{1}{2}-\nu )h\Rightarrow x-x_{j-\frac{5}{2}}=(3-\nu )h,x-x_{j-\frac{3}{2}}=(2-\nu )h,x-x_{j-\frac{1}{2}}=(1-\nu )h,\\ x-x_{j+\frac{1}{2}}=(-\nu )h,x-x_{j+\frac{3}{2}}=(-1-\nu )h,x-x_{j+\frac{5}{2}}=(-2-\nu )h,\\ {\upsilon }_{j+\frac{1}{2}}^{n+1}\equiv \upsilon \left(\right(j+\frac{1}{2}-\nu \left)h\right)={\upsilon }_{j+\frac{1}{2}}^n+u_j(-\nu )h+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}(1-\nu )(-\nu )h^2+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}(2-\nu )(1-\nu )(-\nu )h^3\\ +\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}^3}{4!h^3}(2-\nu )(1-\nu )(-\nu )(-1-\nu )h^4+\frac{{\mathrm{\Delta }}_j^4}{5!h^4}(3-\nu )(2-\nu )(1-\nu )(-\nu )(-1-\nu )h^5,\\ {\overline{u}}_{j+\frac{1}{2}}\equiv \frac{{\upsilon }_{j+\frac{1}{2}}^n-{\upsilon }_{j+\frac{1}{2}}^{n+1}}{\nu h}=u_j+({\mathrm{\Delta }}_{j-\frac{1}{2}}+({\mathrm{\Delta }}_j^2+({\mathrm{\Delta }}_{j-\frac{1}{2}}^3+{\mathrm{\Delta }}_j^4\frac{3-\nu }{5})\frac{-1-\nu }{4}\left)\frac{2-\nu }{3}\right)\frac{1-\nu }{2}.\end{array}$(A3) For flux reconstruction in case a ≥ 0, we need the derivatives of three third-order Newton interpolations as follows (A4) $\begin{array}{rl}u\left(x\right)& ={\upsilon }_x\left(x\right)=\left[\begin{array}{l}{u}_j\frac{d}{dx}(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}\frac{d}{dx}\left((x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\right)+\frac{{\mathrm{\Delta }}_{j-1}^2}{3!h^2}\frac{d}{dx}\left((x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\right)\\ u_j\frac{d}{dx}(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}\frac{d}{dx}\left((x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\right)+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}\frac{d}{dx}\left((x-x_{j-\frac{3}{2}})(x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\right)\\ u_j\frac{d}{dx}(x-x_{j+\frac{1}{2}})+\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}}{2h}\frac{d}{dx}\left((x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})\right)+\frac{{\mathrm{\Delta }}_{j+1}^2}{3!h^2}\frac{d}{dx}\left((x-x_{j-\frac{1}{2}})(x-x_{j+\frac{1}{2}})(x-x_{j+\frac{3}{2}})\right)\end{array}\right],\\ u_{j+\frac{1}{2}}^{n+1}& \equiv u(j+\frac{1}{2}-\nu )=\left[\begin{array}{l}{u}_j\frac{d}{dv}(-(-\nu \left)\right)+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}\frac{d}{dv}(-(1-\nu \left)\right(-\nu \left)\right)h+\frac{{\mathrm{\Delta }}_{j-1}^2}{3!h^2}\frac{d}{dv}(-(2-\nu \left)\right(1-\nu \left)\right(-\nu \left)\right)h^2\\ u_j\frac{d}{dv}(-(-\nu \left)\right)+\frac{{\mathrm{\Delta }}_{j-\frac{1}{2}}}{2h}\frac{d}{dv}(-(1-\nu \left)\right(-\nu \left)\right)h+\frac{{\mathrm{\Delta }}_j^2}{3!h^2}\frac{d}{dv}(-(2-\nu \left)\right(1-\nu \left)\right(-\nu \left)\right)h^2\\ u_j\frac{d}{dv}(-(-\nu \left)\right)+\frac{{\mathrm{\Delta }}_{j+\frac{1}{2}}}{2h}\frac{d}{dv}(-(1-\nu \left)\right(-\nu \left)\right)h+\frac{{\mathrm{\Delta }}_{j+1}^2}{3!h^2}\frac{d}{dv}(-(1-\nu \left)\right(-\nu \left)\right(-1-\nu \left)\right)h^2\end{array}\right],\\ & =\left[\begin{array}{l}{u}_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_{j-1}^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_j^2\frac{d}{dv}\frac{-(-v)(1-v)(2-v)}{3!}\\ u_j+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{d}{dv}\frac{-(-v)(1-v)}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{d}{dv}\frac{-(-1-v)(-v)(1-v)}{3!}\end{array}\right].\end{array}$(A4) And the expansion form of flux reconstruction in Equation (19) (for Full-WENO5) can be given as follows (A5) $\left[\begin{array}{l}{u}_1^{\ast }\\ u_2^{\ast }\\ u_3^{\ast }\end{array}\right]\equiv \{\begin{array}{ll}\left[\begin{array}{l}{u}_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_{j-1}^2\frac{3v^2-6v+2}{3!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_j^2\frac{3v^2-6v+2}{3!}\\ u_j+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{3v^2-1}{3!}\end{array}\right],& a\ge 0,\\ \left[\begin{array}{l}{u}_{j+1}+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_j^2\frac{3v^2-1}{3!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{3v^2+6v+2}{3!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_{j+2}^2\frac{3v^2+6v+2}{3!}\end{array}\right],& a\le 0.\end{array}$(A5) The expansion form of flux reconstruction in Equation (20) (for Full-WENO7) can be given as follows (A6) $\left[\begin{array}{l}{u}_1^{\ast }\\ u_2^{\ast }\\ u_3^{\ast }\\ u_4^{\ast }\end{array}\right]\equiv \{\begin{array}{ll}\left[\begin{array}{l}{u}_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_{j-1}^2\frac{3v^2-6v+2}{3!}+{\mathrm{\Delta }}_{j-\frac{3}{2}}^3\frac{-4v^3+18v^2-22v+6}{4!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_j^2\frac{3v^2-6v+2}{3!}+{\mathrm{\Delta }}_{j-\frac{1}{2}}^3\frac{-4v^3+6v^2-2v-2}{4!}\\ u_j+{\mathrm{\Delta }}_{j-\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_j^2\frac{3v^2-6v+2}{3!}+{\mathrm{\Delta }}_{j+\frac{1}{2}}^3\frac{-4v^3+6v^2-2v-2}{4!}\\ u_j+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{-2v+1}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{3v^2-1}{3!}+{\mathrm{\Delta }}_{j+\frac{3}{2}}^3\frac{-4v^3-6v^2+2v+2}{4!}\end{array}\right],& a\ge 0,\\ \left[\begin{array}{l}{u}_{j+1}+{\mathrm{\Delta }}_{j+\frac{1}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_j^2\frac{3v^2-1}{3!}+{\mathrm{\Delta }}_{j-\frac{1}{2}}^3\frac{-4v^3+6v^2-2v-2}{4!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{3v^2+6v+2}{3!}+{\mathrm{\Delta }}_{j+\frac{1}{2}}^3\frac{-4v^3-6v^2+2v+2}{4!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_{j+1}^2\frac{3v^2+6v+2}{3!}+{\mathrm{\Delta }}_{j+\frac{3}{2}}^3\frac{-4v^3-6v^2+2v+2}{4!}\\ u_{j+1}+{\mathrm{\Delta }}_{j+\frac{3}{2}}\frac{-2v-1}{2!}+{\mathrm{\Delta }}_{j+2}^2\frac{3v^2+6v+2}{3!}+{\mathrm{\Delta }}_{j+\frac{5}{2}}^3\frac{-4v^3-18v^2-22v-6}{4!}\end{array}\right],& a\le 0.\end{array}$(A6)

Appendix B: Runge-Kutta method used in WENO-RK

TVD third-order three-stage RK method

(B1) $\{\begin{array}{l}{u}_t=L\left(u\right),u^{n+1}=u^n+\tau (\frac{\text{1}}{\text{6}}{L}_1+\frac{\text{1}}{\text{6}}{L}_{\text{2}}+\frac{\text{4}}{\text{6}}{L}_3)\\ L_1\equiv L\left(u^n\right),n\\ L_2\equiv L(u^n+\tau L_1),n+1\\ L_3\equiv L(u^n+\tau (\frac{1}{\text{4}}{L}_1+\frac{1}{\text{4}}{L}_2\left)\right),n+\frac{1}{2}\end{array}$(B1)

Nystrom fifth-order six-stage RK method

(B2) $\{\begin{array}{l}{u}_t=L\left(u\right),u^{n+1}=u^n+\tau (\frac{23}{192}{L}_1+\frac{125}{192}{L}_3+\frac{-81}{192}{L}_5+\frac{125}{192}{L}_6)\\ L_1\equiv L\left(u^n\right),n\\ L_2\equiv L(u^n+\frac{\tau }{3}{L}_1),n+\frac{1}{3}\\ L_3\equiv L(u^n+\tau (\frac{4}{25}{L}_1+\frac{6}{25}{L}_2\left)\right),n+\frac{2}{5}\\ L_4\equiv L(u^n+\tau (\frac{1}{4}{L}_1+\frac{-12}{4}{L}_2+\frac{15}{4}{L}_3\left)\right),n+1\\ L_5\equiv L(u^n+\tau (\frac{6}{81}{L}_1+\frac{90}{81}{L}_2+\frac{-50}{81}{L}_3+\frac{8}{81}{L}_4\left)\right),n+\frac{2}{3}\\ L_6\equiv L(u^n+\tau (\frac{6}{75}{L}_1+\frac{36}{75}{L}_2+\frac{10}{75}{L}_3+\frac{8}{75}{L}_4\left)\right),n+\frac{4}{5}\end{array}$(B2)

Butcher seventh-order nine-stage RK method

(B3) $\{\begin{array}{l}{u}_t=L\left(u\right),u^{n+1}=u^n+\tau (\frac{32}{105}{L}_4+\frac{1771561}{6289920}{L}_5+\frac{243}{2560}{L}_6+\frac{16807}{74880}{L}_7+\frac{77}{1440}{L}_8+\frac{11}{270}{L}_9)\\ L_1\equiv L\left(u^n\right),n\\ L_2\equiv L(u^n+\frac{\tau }{6}{L}_1),n+\frac{1}{6}\\ L_3\equiv L(u^n+\frac{\tau }{3}{L}_2),n+\frac{1}{3}\\ L_4\equiv L(u^n+\tau (\frac{1}{8}{L}_1+\frac{3}{8}{L}_3\left)\right),n+\frac{1}{2}\\ L_5\equiv L(u^n+\tau (\frac{148}{1331}{L}_1+\frac{150}{1331}{L}_3+\frac{-56}{1331}{L}_4\left)\right),n+\frac{2}{11}\\ L_6\equiv L(u^n+\tau (\frac{-404}{243}{L}_1+\frac{-170}{27}{L}_3+\frac{4024}{1701}{L}_4+\frac{10648}{1701}{L}_5\left)\right),n+\frac{2}{3}\\ L_7\equiv L(u^n+\tau (\frac{2466}{2401}{L}_1+\frac{1242}{343}{L}_3+\frac{-19176}{16807}{L}_4+\frac{-51909}{16807}{L}_5+\frac{1053}{2401}{L}_6\left)\right),n+\frac{6}{7}\\ L_8\equiv L(u^n+\tau (\frac{5}{154}{L}_1+\frac{96}{539}{L}_4+\frac{-1815}{20384}{L}_5+\frac{-405}{2464}{L}_6+\frac{49}{1144}{L}_7\left)\right),n\\ L_9\equiv L(u^n+\tau (\frac{-113}{32}{L}_1+\frac{-195}{22}{L}_3+\frac{32}{7}{L}_4+\frac{29403}{3584}{L}_5+\frac{-729}{512}{L}_6+\frac{1029}{1408}{L}_7+\frac{21}{16}{L}_8\left)\right),n+1\end{array}$(B3)

According to (Butcher, 2016): above order 4, it is no longer possible to obtain order s with just s stages. For order 5, six stages are required, and for order 7, nine stages are required.