A hybrid Hermite WENO scheme for hyperbolic conservation laws

TL;DR

This paper introduces a hybrid HWENO scheme combining shock detection and adaptive reconstruction for hyperbolic conservation laws, improving accuracy and efficiency near discontinuities while maintaining robustness and compactness.

Contribution

The paper presents a novel hybrid HWENO scheme that adaptively switches between high-order linear and HWENO reconstructions based on shock detection, enhancing performance over existing methods.

Findings

Demonstrates high accuracy and resolution in benchmark tests.

Shows robustness and efficiency in handling discontinuities.

Maintains scheme compactness and reduces oscillations.

Abstract

In this paper, we propose a hybrid finite volume Hermite weighted essentially non-oscillatory (HWENO) scheme for solving one and two dimensional hyperbolic conservation laws. The zeroth-order and the first-order moments are used in the spatial reconstruction, with total variation diminishing Runge-Kutta time discretization. The main idea of the hybrid HWENO scheme is that we first use a shock-detection technique to identify the troubled cell, then, if the cell is identified as a troubled cell, we would modify the first order moment in the troubled cell and employ HWENO reconstruction in spatial discretization; otherwise, we directly use high order linear reconstruction. Unlike other HWENO schemes, we borrow the thought of limiter for discontinuous Galerkin (DG) method to control the spurious oscillations, after this procedure, the scheme would avoid the oscillations by using HWENO reconstruction nearby discontinuities and have higher efficiency for using linear approximation straightforwardly in the smooth regions. In addition, the hybrid HWENO scheme still keeps the compactness. A collection of benchmark numerical tests for one and two dimensional cases are performed to demonstrate the numerical accuracy, high resolution and robustness of the proposed scheme.

Tables

Table 1. Table 3.1: 1D-Burgers’ equation: initial data u ​ ( x , 0 ) = 0.5 + s ​ i ​ n ​ ( π ​ x ) 𝑢 𝑥 0 0.5 𝑠 𝑖 𝑛 𝜋 𝑥 u(x,0)=0.5+sin(\pi x) . HWENO and WENO schemes. T = 0.5 / π 𝑇 0.5 𝜋 T=0.5/\pi . L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} errors and orders. Uniform meshes with N 𝑁 N cells.

	HWENO scheme				WENO scheme
$N$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
10	1.21E-02		1.00E-01		1.90E-02		7.46E-02
20	1.06E-03	3.52	1.00E-02	3.32	2.06E-03	3.20	1.23E-02	2.60
40	4.23E-05	4.65	5.25E-04	4.26	1.22E-04	4.08	1.05E-03	3.55
80	1.24E-06	5.09	1.70E-05	4.95	4.36E-06	4.80	4.78E-05	4.46
160	4.26E-08	4.87	4.84E-07	5.13	1.64E-07	4.74	1.41E-06	5.09
320	1.13E-09	5.24	1.43E-08	5.08	4.78E-09	5.10	7.35E-08	4.26
	Hybrid HWENO scheme				Hybrid WENO scheme
$N$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
10	1.18E-03		6.00E-03		1.44E-02		7.32E-02
20	4.18E-05	4.82	3.69E-04	4.02	1.58E-03	3.19	1.47E-02	2.31
40	8.51E-07	5.62	1.14E-05	5.02	9.45E-05	4.06	1.29E-03	3.51
80	1.46E-08	5.87	2.26E-07	5.65	2.39E-06	5.30	3.11E-05	5.38
160	2.66E-10	5.78	3.59E-09	5.98	7.15E-08	5.06	9.40E-07	5.05
320	5.65E-12	5.56	5.93E-11	5.92	2.12E-09	5.08	2.82E-08	5.06

Table 2. Table 3.2: 1D-Euler equations: initial data ρ ​ ( x , 0 ) = 1 + 0.2 ​ sin ⁡ ( π ​ x ) 𝜌 𝑥 0 1 0.2 𝜋 𝑥 \rho(x,0)=1+0.2\sin(\pi x) , μ ​ ( x , 0 ) = 1 𝜇 𝑥 0 1 \mu(x,0)=1 and p ​ ( x , 0 ) = 1 𝑝 𝑥 0 1 p(x,0)=1 . HWENO and WENO schemes. T = 2 𝑇 2 T=2 . L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} errors and orders. Uniform meshes with N 𝑁 N cells.

	HWENO scheme				WENO scheme
$N$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
10	3.98E-03		6.25E-03		1.13E-02		1.66E-02
20	1.39E-04	4.84	2.50E-04	4.64	6.26E-04	4.17	9.94E-04	4.06
40	4.00E-06	5.12	8.18E-06	4.93	2.04E-05	4.94	3.72E-05	4.74
80	1.22E-07	5.04	2.43E-07	5.08	6.45E-07	4.98	1.21E-06	4.94
160	3.73E-09	5.03	6.71E-09	5.18	2.01E-08	5.01	3.67E-08	5.05
320	1.11E-10	5.07	1.91E-10	5.13	6.09E-10	5.04	1.01E-09	5.19
	Hybrid HWENO scheme				Hybrid WENO scheme
$N$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
10	1.82E-06		2.82E-06		2.55E-03		4.25E-03
20	3.71E-08	5.62	5.73E-08	5.62	8.94E-05	4.83	1.47E-04	4.86
40	1.02E-09	5.18	1.60E-09	5.16	2.91E-06	4.94	4.67E-06	4.97
80	3.10E-11	5.05	4.86E-11	5.04	9.22E-08	4.98	1.47E-07	4.99
160	9.61E-13	5.01	1.51E-12	5.01	2.90E-09	4.99	4.59E-09	5.00
320	3.00E-14	5.00	4.71E-14	5.00	9.10E-11	5.00	1.43E-10	5.00

Table 3. Table 3.3: 2D-Burgers’ equation: initial data u ​ ( x , y , 0 ) = 0.5 + s ​ i ​ n ​ ( π ​ ( x + y ) / 2 ) 𝑢 𝑥 𝑦 0 0.5 𝑠 𝑖 𝑛 𝜋 𝑥 𝑦 2 u(x,y,0)=0.5+sin(\pi(x+y)/2) . Hybrid HWENO and WENO schemes. T = 0.5 / π 𝑇 0.5 𝜋 T=0.5/\pi . L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} errors and orders. Uniform meshes with N x × N y subscript 𝑁 𝑥 subscript 𝑁 𝑦 N_{x}\times N_{y} cells.

	Hybrid HWENO scheme				WENO scheme
$N_x\times N_y$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
$40\times 40$	7.07E-05		6.32E-04		8.20E-05		6.74E-04
$80\times 80$	3.95E-06	4.16	4.28E-05	3.88	4.06E-06	4.34	3.91E-05	4.11
$120\times 120$	7.31E-07	4.16	7.62E-06	4.26	6.30E-07	4.60	5.67E-06	4.76
$160\times 160$	2.19E-07	4.19	2.30E-06	4.16	1.66E-07	4.64	1.42E-06	4.81
$200\times 200$	8.67E-08	4.15	8.96E-07	4.23	5.65E-08	4.82	4.96E-07	4.72
$240\times 240$	4.07E-08	4.14	4.19E-07	4.17	2.28E-08	4.99	2.22E-07	4.40

Table 4. Table 3.4: 2D-Euler equations: initial data ρ ​ ( x , y , 0 ) = 1 + 0.2 ​ sin ⁡ ( π ​ ( x + y ) ) 𝜌 𝑥 𝑦 0 1 0.2 𝜋 𝑥 𝑦 \rho(x,y,0)=1+0.2\sin(\pi(x+y)) , μ ​ ( x , y , 0 ) = 1 𝜇 𝑥 𝑦 0 1 \mu(x,y,0)=1 , ν ​ ( x , y , 0 ) = 1 𝜈 𝑥 𝑦 0 1 \nu(x,y,0)=1 and p ​ ( x , y , 0 ) = 1 𝑝 𝑥 𝑦 0 1 p(x,y,0)=1 . Hybrid HWENO and WENO schemes. T = 2 𝑇 2 T=2 . L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} errors and orders. Uniform meshes with N x × N y subscript 𝑁 𝑥 subscript 𝑁 𝑦 N_{x}\times N_{y} cells.

	Hybrid HWENO scheme				WENO scheme
$N_x\times N_y$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
$40\times 40$	5.66E-06		8.90E-06		5.11E-06		1.14E-05
$80\times 80$	1.86E-07	4.92	2.93E-07	4.93	1.19E-07	5.43	3.12E-07	5.18
$120\times 120$	2.61E-08	4.85	4.10E-08	4.85	1.38E-08	5.31	3.79E-08	5.20
$160\times 160$	6.66E-09	4.75	1.05E-08	4.75	3.09E-09	5.21	7.58E-09	5.59
$200\times 200$	2.36E-09	4.65	3.71E-09	4.65	9.73E-10	5.18	2.03E-09	5.91
$240\times 240$	1.03E-09	4.56	1.62E-09	4.55	3.79E-10	5.17	7.21E-10	5.67

Table 5. Table 3.5: Isentropic vortex test. Hybrid HWENO and WENO schemes. T = 10 𝑇 10 T=10 . L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} errors and orders. Uniform meshes with N x × N y subscript 𝑁 𝑥 subscript 𝑁 𝑦 N_{x}\times N_{y} cells.

	Hybrid HWENO scheme				WENO scheme
$N_x\times N_y$ cells	$L^1$ error	order	$L^{\infty }$error	order	$L^1$ error	order	$L^{\infty }$ error	order
$40\times 40$	1.82E-04		4.53E-03		1.39E-03		2.73E-02
$80\times 80$	8.18E-06	4.47	1.34E-04	5.08	7.00E-05	4.32	1.73E-03	3.98
$120\times 120$	1.16E-06	4.82	2.01E-05	4.67	1.09E-05	4.58	1.89E-04	5.46
$160\times 160$	2.82E-07	4.92	5.01E-06	4.83	2.70E-06	4.85	3.86E-05	5.52
$200\times 200$	9.37E-08	4.93	1.68E-06	4.89	9.12E-07	4.87	1.66E-05	3.78
$240\times 240$	3.81E-08	4.93	6.88E-07	4.90	3.77E-07	4.85	7.51E-06	4.36

Equations

$$\left\{\begin{array}[]{ll}u_{t}+f(u)_{x}=0,\\ u(x,0)=u_{0}(x).\\ \end{array}\right.$$

$$\{\begin{array}{rl}\frac{d{\overline{u}}_i(t)}{dt}& =-\frac{1}{\Delta x}\left({\hat{f}}_{i+1/2}-{\hat{f}}_{i-1/2}\right),\\ \frac{d{\overline{v}}_i(t)}{dt}& =-\frac{1}{2\Delta x}\left({\hat{f}}_{i-1/2}+{\hat{f}}_{i+1/2}\right)+\frac{1}{\Delta x}{F}_i(u),\end{array}$$

$${\hat{f}}_{i+1/2}=\frac{1}{2}\left(f(u_{i+1/2}^{-})+f(u_{i+1/2}^{+})\right)-\frac{\alpha }{2}\left(u_{i+1/2}^{+}-u_{i+1/2}^{-}\right),$$

$$F_i(u)=\frac{1}{\Delta x}{\int }_{I_i}f(u)dx\approx \sum _{l=1}^4{\omega }_{l}f(u(x_l^G,t)).$$

$$x_1^G=x_{i-1/2},x_2^G=x_{i-\sqrt{5}/10},x_3^G=x_{i+\sqrt{5}/10},x_4^G=x_{i+1/2},$$

$$\frac{\left|{\int }_{\partial I_i^{-}}(u^h{|}_{I_i}-u^h{|}_{I_{n_i}})ds\right|}{h_i^{\frac{k+1}{2}}|\partial I_i^{-}|||u^h{|}_{I_i}||}>1,$$

$$p_i^3(x)=u_i^{(0)}+u_i^{(1)}\left(\frac{x-x_i}{\Delta x}\right)+u_i^{(2)}\left[{\left(\frac{x-x_i}{\Delta x}\right)}^2-\frac{1}{12}\right]+u_i^{(3)}\left[{\left(\frac{x-x_i}{\Delta x}\right)}^3-\frac{3}{20}\left(\frac{x-x_i}{\Delta x}\right)\right],$$

$$\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_i^3(x)dx={\overline{u}}_{i+j},j=-1,0,1,\frac{1}{\Delta x}{\int }_{I_i}{p}_i^3(x)\frac{x-x_i}{\Delta x}dx={\overline{v}}_i.$$

$$u_i^{(0)}={\overline{u}}_i,u_i^{(1)}=12{\overline{v}}_i,u_i^{(2)}=\frac{1}{2}({\overline{u}}_{i-1}-2{\overline{u}}_i+{\overline{u}}_{i+1}),u_i^{(3)}=-\frac{5}{11}({\overline{u}}_{i-1}+24{\overline{v}}_i-{\overline{u}}_{i+1}).$$

$$\begin{array}{rl}\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_1(x)dx& ={\overline{u}}_{i+j},j=-1,0,\frac{1}{\Delta x}{\int }_{I_{i-1}}{p}_1(x)\frac{x-x_{i-1}}{\Delta x}dx={\overline{v}}_{i-1},\\ \frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_2(x)dx& ={\overline{u}}_{i+j},j=-1,0,1,\\ \frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_3(x)dx& ={\overline{u}}_{i+j},j=0,1,\frac{1}{\Delta x}{\int }_{I_{i+1}}{p}_3(x)\frac{x-x_{i+1}}{\Delta x}dx={\overline{v}}_{i+1},\end{array}$$

$$\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_0(x)dx={\overline{u}}_{i+j},j=-1,0,1,\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_0(x)\frac{x-x_{i+j}}{\Delta x}dx={\overline{v}}_{i+j},j=-1,1.$$

$$\begin{array}{rl}\frac{1}{\Delta x}{\int }_{I_i}{p}_1(x)\frac{x-x_i}{\Delta x}dx& =\frac{1}{6}{\overline{u}}_i-\frac{1}{6}{\overline{u}}_{i-1}-{\overline{v}}_{i-1},\\ \frac{1}{\Delta x}{\int }_{I_i}{p}_2(x)\frac{x-x_i}{\Delta x}dx& =\frac{1}{24}{\overline{u}}_{i+1}-\frac{1}{24}{\overline{u}}_{i-1},\\ \frac{1}{\Delta x}{\int }_{I_i}{p}_3(x)\frac{x-x_i}{\Delta x}dx& =\frac{1}{6}{\overline{u}}_{i+1}-\frac{1}{6}{\overline{u}}_i-{\overline{v}}_{i+1},\\ \frac{1}{\Delta x}{\int }_{I_i}{p}_0(x)\frac{x-x_i}{\Delta x}dx& =\frac{5}{76}{\overline{u}}_{i+1}-\frac{5}{76}{\overline{u}}_{i-1}-\frac{11}{38}{\overline{v}}_{i-1}-\frac{11}{38}{\overline{v}}_{i+1}.\end{array}$$

$$\frac{1}{\Delta x}{\int }_{I_i}{p}_0(x)\frac{x-x_i}{\Delta x}dx=\frac{1}{\Delta x}\sum _{n=1}^3{\gamma }_n{\int }_{I_i}{p}_n(x)\frac{x-x_i}{\Delta x}dx,$$

$${\beta }_n=\sum _{\alpha =1}^r{\int }_{I_i}{\Delta x}^{2\alpha -1}(\frac{d^{\alpha }{p}_n(x)}{d{x}^{\alpha }}{)}^{2}dx,n=1,2,3.$$

$$\{\begin{array}{rl}{\beta }_1& =4({\overline{u}}_{i-1}-{\overline{u}}_i+6{\overline{v}}_{i-1}{)}^2+\frac{13}{3}({\overline{u}}_{i-1}-{\overline{u}}_i+12{\overline{v}}_{i-1}{)}^2,\\ {\beta }_2& =\frac{1}{4}({\overline{u}}_{i-1}-{\overline{u}}_{i+1}{)}^2+\frac{13}{12}({\overline{u}}_{i-1}-2{\overline{u}}_i+{\overline{u}}_{i+1}{)}^2,\\ {\beta }_3& =4({\overline{u}}_i-{\overline{u}}_{i+1}+6{\overline{v}}_{i+1}{)}^2+\frac{13}{3}({\overline{u}}_i-{\overline{u}}_{i+1}+12{\overline{v}}_{i+1}{)}^2.\end{array}$$

$${\omega }_n=\frac{{\bar{\omega }}_n}{{\sum }_{l=1}^3{\bar{\omega }}_l},\text{mbox}with{\bar{\omega }}_n=\frac{{\gamma }_n}{({\beta }_n+\epsilon {)}^2},n=1,2,3,$$

$${\overline{v}}_i=\frac{1}{\Delta x}\sum _{n=1}^3{\omega }_n{\int }_{I_i}{p}_n(x)\frac{x-x_i}{\Delta x}dx.$$

$$\begin{array}{rl}\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_1(x)dx& ={\overline{u}}_{i+j},\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_1(x)\frac{x-x_{i+j}}{\Delta x}dx={\overline{v}}_{i+j},j=-1,0,\\ \frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_2(x)dx& ={\overline{u}}_{i+j},j=-1,0,1,\frac{1}{\Delta x}{\int }_{I_i}{p}_2(x)\frac{x-x_i}{\Delta x}dx={\overline{v}}_i,\\ \frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_3(x)dx& ={\overline{u}}_{i+j},\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_3(x)\frac{x-x_{i+j}}{\Delta x}dx={\overline{v}}_{i+j},j=0,1,\end{array}$$

$$\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_0(x)dx={\overline{u}}_{i+j},\frac{1}{\Delta x}{\int }_{I_{i+j}}{p}_0(x)\frac{x-x_{i+j}}{\Delta x}dx={\overline{v}}_{i+j},j=-1,0,1.$$

$$\begin{array}{rl}{p}_1(x_{i+1/2})& =\frac{3}{4}{\overline{u}}_{i-1}+\frac{1}{4}{\overline{u}}_i+\frac{7}{2}{\overline{v}}_{i-1}+\frac{23}{2}{\overline{v}}_i,\\ p_2(x_{i+1/2})& =\frac{2}{33}{\overline{u}}_{i-1}+\frac{5}{6}{\overline{u}}_i+\frac{7}{66}{\overline{u}}_{i+1}+\frac{60}{11}{\overline{v}}_i,\\ p_3(x_{i+1/2})& =\frac{1}{2}{\overline{u}}_i+\frac{1}{2}{\overline{u}}_{i+1}+2{\overline{v}}_i-2{\overline{v}}_{i+1},\\ p_0(x_{i+1/2})& =\frac{13}{108}{\overline{u}}_{i-1}+\frac{7}{12}{\overline{u}}_i+\frac{8}{27}{\overline{u}}_{i+1}+\frac{25}{54}{\overline{v}}_{i-1}+\frac{241}{54}{\overline{v}}_i-\frac{28}{27}{\overline{v}}_{i+1}.\end{array}$$

$$p_0(x_{i+1/2})=\sum _{n=1}^3{\gamma }_n{p}_n(x_{i+1/2}),$$

$$\{\begin{array}{rl}{\beta }_1& =\frac{1}{16}({\overline{u}}_{i-1}-{\overline{u}}_i+6{\overline{v}}_{i-1}+54{\overline{v}}_i{)}^2+\frac{13}{48}(15{\overline{u}}_{i-1}-15{\overline{u}}_i+66{\overline{v}}_{i-1}+114{\overline{v}}_i{)}^2+\\ & \frac{3905}{16}({\overline{u}}_{i-1}-{\overline{u}}_i+6{\overline{v}}_{i-1}+6{\overline{v}}_i{)}^2,\\ {\beta }_2& =\frac{1}{484}({\overline{u}}_{i-1}-{\overline{u}}_{i+1}-240{\overline{v}}_i{)}^2+\frac{13}{12}({\overline{u}}_{i-1}-2{\overline{u}}_i+{\overline{u}}_{i+1}{)}^2+\\ & \frac{355}{44}({\overline{u}}_{i-1}-{\overline{u}}_{i+1}+24{\overline{v}}_i{)}^2,\\ {\beta }_3& =\frac{1}{16}({\overline{u}}_i-{\overline{u}}_{i+1}+54{\overline{v}}_i+6{\overline{v}}_{i+1}{)}^2+\frac{13}{48}(15{\overline{u}}_i-15{\overline{u}}_{i+1}+114{\overline{v}}_i+66{\overline{v}}_{i+1}{)}^2+\\ & \frac{3905}{16}({\overline{u}}_i-{\overline{u}}_{i+1}+6{\overline{v}}_i+6{\overline{v}}_{i+1}{)}^2,\end{array}$$

$${\omega }_n=\frac{{\bar{\omega }}_n}{{\sum }_{l=1}^3{\bar{\omega }}_l},\text{mbox}with{\bar{\omega }}_n=\frac{{\gamma }_n}{({\beta }_n+\epsilon {)}^2},n=1,2,3.$$

$$u_{i+1/2}^{-}=\sum _{n=1}^3{\omega }_n{p}_n(x_{i+1/2}).$$

$$u_{i-1/2}^{+}=p_0(x_{i-1/2})=\frac{8}{27}{\overline{u}}_{i-1}+\frac{7}{12}{\overline{u}}_i+\frac{13}{108}{\overline{u}}_{i+1}+\frac{28}{27}{\overline{v}}_{i-1}-\frac{241}{54}{\overline{v}}_i-\frac{25}{54}{\overline{v}}_{i+1},$$

$$u_{i+1/2}^{-}=p_0(x_{i+1/2})=\frac{13}{108}{\overline{u}}_{i-1}+\frac{7}{12}{\overline{u}}_i+\frac{8}{27}{\overline{u}}_{i+1}+\frac{25}{54}{\overline{v}}_{i-1}+\frac{241}{54}{\overline{v}}_i-\frac{28}{27}{\overline{v}}_{i+1}.$$

$$\begin{array}{rl}{u}_{i-\sqrt{5}/10}=p_0(x_{i-\sqrt{5}/10})& =-(\frac{101}{5400}\sqrt{5}+\frac{1}{24}){\overline{u}}_{i-1}+\frac{13}{12}{\overline{u}}_i+(\frac{101}{5400}\sqrt{5}-\frac{1}{24}){\overline{u}}_{i+1}-\\ & (\frac{3}{20}+\frac{841}{13500}\sqrt{5}){\overline{v}}_{i-1}-\frac{10289}{6750}\sqrt{5}{\overline{v}}_i+(\frac{3}{20}-\frac{841}{13500}\sqrt{5}){\overline{v}}_{i+1},\end{array}$$

$$\begin{array}{rl}{u}_{i+\sqrt{5}/10}=p_0(x_{i+\sqrt{5}/10})& =(\frac{101}{5400}\sqrt{5}-\frac{1}{24}){\overline{u}}_{i-1}+\frac{13}{12}{\overline{u}}_i-(\frac{101}{5400}\sqrt{5}+\frac{1}{24}){\overline{u}}_{i+1}+\\ & (\frac{841}{13500}\sqrt{5}-\frac{3}{20}){\overline{v}}_{i-1}+\frac{10289}{6750}\sqrt{5}{\overline{v}}_i+(\frac{3}{20}+\frac{841}{13500}\sqrt{5}){\overline{v}}_{i+1}.\end{array}$$

$$\displaystyle\left\{\begin{array}[]{lll}u^{(1)}&=&u^{n}+\Delta tL(u^{n}),\\ u^{(2)}&=&\frac{3}{4}u^{n}+\frac{1}{4}u^{(1)}+\frac{1}{4}\Delta tL(u^{(1)}),\\ u^{(n+1)}&=&\frac{1}{3}u^{n}+\frac{2}{3}u^{(2)}+\frac{2}{3}\Delta tL(u^{(2)}).\end{array}\right.$$

$$\left\{\begin{array}[]{ll}u_{t}+f(u)_{x}+g(u)_{y}=0,\\ u(x,y,0)=u_{0}(x,y).\\ \end{array}\right.$$

Full text

**A hybrid Hermite WENO scheme for hyperbolic conservation laws 111The research is partly supported by Science Challenge Project, No. TZ2016002, NSAF grant U1630247 and NSFC grant 11571290. **

Zhuang Zhao222School of Mathematical Sciences, Xiamen University, Xiamen, Fujian 361005, P.R. China. E-mail: zzhao@stu.xmu.edu.cn., Yibing Chen333Institute of Applied Physics and Computational Mathematics, Beijing 100094, China. E-mail: chen_yibing@iapcm.ac.cn. and Jianxian Qiu444School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modeling and High-Performance Scientific Computing, Xiamen University, Xiamen, Fujian 361005, P.R. China. E-mail: jxqiu@xmu.edu.cn.

Abstract

In this paper, we propose a hybrid finite volume Hermite weighted essentially non-oscillatory (HWENO) scheme for solving one and two dimensional hyperbolic conservation laws, which would be the fifth order accuracy in the one dimensional case, while is the fourth order accuracy for two dimensional problems. The zeroth-order and the first-order moments are used in the spatial reconstruction, with total variation diminishing Runge-Kutta time discretization. Unlike the original HWENO schemes [28, 29] using different stencils for spatial discretization, we borrow the thought of limiter for discontinuous Galerkin (DG) method to control the spurious oscillations, after this procedure, the scheme would avoid the oscillations by using HWENO reconstruction nearby discontinuities, and using linear approximation straightforwardly in the smooth regions is to increase the efficiency of the scheme. Moreover, the scheme still keeps the compactness as only immediate neighbor information is needed in the reconstruction. A collection of benchmark numerical tests for one and two dimensional cases are performed to demonstrate the numerical accuracy, high resolution and robustness of the proposed scheme.

Key Words: Hermite WENO scheme, hyperbolic conservation laws, discontinuous Galerkin method, limiter

AMS(MOS) subject classification: 65M60, 35L65

1 Introduction

In this paper, we design a hybrid Hermite weighted essentially non-oscillatory (HWENO) scheme in the finite volume framework, which is the fifth order accuracy in the one dimensional case and the fourth order accuracy for two dimensional problems. The HWENO scheme was derived from essentially non-oscillatory (ENO) and weighted essentially non-oscillatory (WENO) schemes, which have been widely applied for nonlinear hyperbolic conservation laws in recent decades. In 1985, Harten and Osher [11] constructed a weaker version of the total variation diminishing (TVD) criterion [10], which gave a framework for the reconstruction to design higher order ENO schemes. Then, Harten et al. [13] developed the finite volume ENO schemes for solving one dimensional problems, in which they selected stencil adaptively in terms of the local smoothness, and in [12], Harten extended the finite volume ENO schemes to two dimensional hyperbolic conservation laws. In 1994, the first WENO scheme was constructed by Liu, Osher and Chan [20] mainly based on ENO scheme, where they used a nonlinear convex combination of all the candidate stencils to obtain higher order accuracy in smooth regions, and it was a third-order finite volume method in the one dimensional case. In 1996, Jiang and Shu [18] proposed the third and fifth-order finite difference WENO schemes in multi-space dimensions, which gave a general framework to design the smoothness indicators and nonlinear weights. Ever since then, the WENO schemes have been further developed in the finite difference and finite volume frameworks presented in the literature [14, 21, 32, 42, 44], and more detailed review for WENO schemes can refer to [34].

However, if we want to achieve higher order accuracy for WENO schemes, we need to enlarge the stencil for the spatial reconstruction. To make the stencil more compact, Qiu and Shu [28, 29] developed the WENO methodology, which were first taken as limiters for Runge-Kutta discontinuous Galerkin methods, termed as Hermite WENO (HWENO) schemes. After this, many HWENO schemes were developed for solving hyperbolic conservations laws [39, 7, 16, 35, 24, 43, 36, 8, 25]. The HWENO schemes can achieve higher order accuracy than standard WENO schemes on the same stencils. On the other hand, to avoid spurious oscillations, all reconstructions of WENO and HWENO schemes are based on local characteristic decompositions for the systems. Actually, the cost to compute the nonlinear weights and local characteristic decompositions is very high for WENO and HWENO schemes. To increase the efficiency, Pirozzoli [27] designed an efficient hybrid compact-WENO scheme, on which they chose compact up-wind schemes in the smooth regions, while used WENO schemes in the discontinuous regions. Hill and Pullin [15] developed a hybrid scheme, which combined the tuned center-difference schemes with WENO schemes. The purpose was to expect the nonlinear weights would be achieved automatically in the smooth regions away from shocks, but a switching principle was still necessary. Next, Li and Qiu [23] studied the hybrid WENO scheme using different switching principles [30], which illustrated that the troubled-cells indicator introduced by Krivodonova et al.[19] (KXRCF) has the ability to identify the discontinuities well. Other different schemes introduced by [40, 38, 45, 17] for hyperbolic conservation laws also showed the good performances of the KXRCF troubled-cell indicator. In this paper, we would choose it as the indicator to identify the troubled-cell where the solutions may be discontinuous. The main idea for the hybrid WENO schemes [23, 46] was that they both used non-linear WENO reconstruction near discontinuities, while employed upwind linear approximation directly in the smooth regions. We also notice that using upwind linear approximation in the smooth regions is one of the choice, and the other methods can also be used, such as B. Costa and W.S. Don used spectral method [5] and central finite difference scheme [6] in the smooth regions, which also can increase the efficiency obviously.

The hybrid HWENO scheme in this paper is different from the original HWENO schemes [28, 29] as we take the thought of limiter for discontinuous Galerkin (DG) method to control the spurious oscillations. Since the solutions of nonlinear hyperbolic conservation laws often contain discontinuities, the derivative equations for the HWENO schemes need to deal with the derivatives or the first order moments, which would be relatively large nearby discontinuities. Therefore, the HWENO schemes listed by [28, 29, 39, 35, 24, 43, 36, 8, 25] all used different stencils to discretize the space for the original equations and the derivative equations, respectively. The variables of the derivative equations for the hybrid HWENO scheme are the first order moments, which also can be seen in the discontinuous Galerkin (DG) methods [1, 2, 3, 4] and other HWENO schemes [41, 26, 36]. In one sense, these HWENO schemes can be seen as an extension by DG methods, and Dumbser et al. [9] gave a general and unified framework to define the numerical scheme extended by DG method, termed as $P_{N}P_{M}$ method, in which $P_{N}$ represents a piecewise polynomial of degree $N$ used as test functions in DG method and $P_{M}$ is a polynomial of degree $M$ reconstructed by the test functions of degree $N$ for computing the numerical fluxes, and $M\geq N$. It is well known that DG methods use limiters to modify the first order or higher order moments in the discontinuous regions, therefore, we adopt this thought by adding a HWENO limiters [28, 29, 22] to modify the first order moments nearby discontinuities, and use HWENO procedure to reconstruct the point values on the interface of troubled cell.

The main procedures of the hybrid HWENO scheme are given as follows. At first, we use the KXRCF troubled-cell indicator [19] to identify troubled cells, then, we modify the first order moments in the troubled cells by the HWENO limiters [28, 29, 22]. After we modify the first order moments for all troubled cells in the computing domain, we would employ HWENO reconstruction at the points on the interface of the troubled cell, but use linear approximation at the internal points for spatial discretization, otherwise we directly use high order linear approximation. For the systems, all HWENO reconstructions are based on local characteristic decompositions to avoid spurious oscillations just like the classical WENO scheme [18]. Compared with other HWENO schemes [28, 29, 39, 7, 16, 35, 24, 43, 36, 8, 25], we borrow the idea of limiter for discontinuous Galerkin (DG) method to control the spurious oscillations, which would have two advantages. The one is to control the oscillations, and another one is to increase the efficiency for we directly use linear approximation in the smooth regions, where the limiter doesn’t work under this circumstance. In short, the hybrid HWENO scheme avoids the spurious oscillations well and has higher efficiency, while it still keeps the compactness as only immediate neighbor information is needed in the reconstruction.

The organization of the paper is as follows: in Section 2, we present the construction and implementation of the finite volume hybrid HWENO scheme in the one and two dimensional cases in detail. In Section 3, extensive numerical tests are performed to demonstrate the numerical accuracy, efficiency and robustness of the proposed scheme. Concluding remarks are given in Section 4.

2 Hybrid Hermite WENO scheme

In this section, we will introduce the detailed implementation procedures of the hybrid HWENO scheme for one and two dimensional hyperbolic conservational laws, which would be the fifth order accuracy in the one dimensional case, while is the fourth order accuracy for two dimensional problems.

2.1 One dimensional case

We first consider one dimensional scalar hyperbolic conservation laws

[TABLE]

For simplicity, we divide the spatial domain with uniform meshes by $I_{i}=[x_{i-1/2},x_{i+1/2}]$, where the cell center is $x_{i}=\frac{x_{i-1/2}+x_{i+1/2}}{2}$, and the mesh size is set as $\Delta x=x_{i+1/2}-x_{i-1/2}$.

To design a HWENO scheme, we multiply the governing equation (2.1) by $\frac{1}{\Delta x}$ and $\frac{x-x_{i}}{{(\Delta x)}^{2}}$, respectively, and integrate them over $I_{i}$, then, apply the integration by parts and use the numerical flux to approximate the values of the flux at the interface of $I_{i}$, finally, the semi-discrete finite volume HWENO scheme is

[TABLE]

with initial conditions $\overline{u}_{i}(0)=\frac{1}{\Delta x}\int_{I_{i}}u_{0}(x)dx$ and $\overline{v}_{i}(0)=\frac{1}{\Delta x}\int_{I_{i}}u_{0}(x)\frac{x-x_{i}}{\Delta x}dx$, where $\overline{u}_{i}(t)$ is the cell average as $\frac{1}{\Delta x}\int_{I_{i}}u(x,t)dx$ and $\overline{v}_{i}(t)$ is the first order moment as $\frac{1}{\Delta x}\int_{I_{i}}u(x,t)\frac{x-x_{i}}{\Delta x}dx$. Here, $\hat{f}_{i+1/2}$ is the numerical flux which is the approximation to the values of the flux $f(u)$ at the interface point $x_{i+1/2}$ and $F_{i}(u)$ is the numerical integration for the flux $f(u)$. We adopt the Lax-Friedrichs numerical flux method to define the $\hat{f}_{i+1/2}$:

[TABLE]

where $\alpha=\max_{u}|f^{\prime}(u)|$. The numerical integration $F_{i}(u)$ is approximated by four-point Gauss-Lobatto quadrature formula, and the specific expression is given as follows,

[TABLE]

Here, $\omega_{1}=\omega_{4}=\frac{1}{12}$ and $\omega_{2}=\omega_{3}=\frac{5}{12}$. The quadrature points on the cell $I_{i}$ are

[TABLE]

where $x_{i+a}$ is defined as $x_{i}+a\Delta x$.

The general frameworks for the hybrid HWENO scheme are given as follows. In Steps 1 and 2, we’ll introduce the procedures of the spatial reconstruction for the semi-discrete scheme (2.2). In Step 3, the equations (2.2) is discretized in time by the third order TVD Runge-Kutta methodology [31].

Step 1. Identify the troubled-cell and modify the first order moment in the troubled-cell.

Step 1.1. Identify the troubled-cell.

Troubled-cell means that the solution in the cell may be discontinuous, and in [30], Qiu and Shu investigated different troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods. As suggested in [30], we use KXRCF troubled-cell indicator by Krivodonova et al. [19] to identify the discontinuities. We first divide the interface of the cell $I_{i}$ into two parts $\partial I_{i}^{-}$ and $\partial I_{i}^{+}$, in which the flow is into ($\overrightarrow{v}\cdot\overrightarrow{n}<0$, $\overrightarrow{n}$ is the normal vector to $\partial I_{i}$) and out ($\overrightarrow{v}\cdot\overrightarrow{n}>0$) of $I_{i}$, respectively. The cell $I_{i}$ is finally identified as a troubled cell, if:

[TABLE]

where $h_{i}$ is the radius of the circumscribed circle in the cell $I_{i}$, $I_{n_{i}}$ is the neighbor of $I_{i}$ on the side of $\partial I_{i}^{-}$, the norm is $L_{\infty}$ norm in the one dimensional case and $k$ is the degree of the polynomial $u_{h}$ approximating to $u(x)$, and we take $k=2$ in this paper. We should reconstruct the polynomial $u_{h}$ which is used only in the troubled-cell indicator (2.3) to identify troubled cell, not in the reconstruction procedure for solution.

We use the information $\overline{u}_{i-1}$, $\overline{u}_{i}$, $\overline{u}_{i+1}$ and $\overline{v}_{i}$ to reconstruct a cubic polynomial $p_{i}^{3}(x)$ on the orthogonal basis function space $\left\{1,\frac{x-x_{i}}{\Delta x},\left(\frac{x-x_{i}}{\Delta x}\right)^{2}-\frac{1}{12},\left(\frac{x-x_{i}}{\Delta x}\right)^{3}-\frac{3}{20}\left(\frac{x-x_{i}}{\Delta x}\right)\right\}$, and the expressions is

[TABLE]

satisfying

[TABLE]

We have:

[TABLE]

$u_{h}$ is taken as $u_{i}^{(0)}+u_{i}^{(1)}\left(\frac{x-x_{i}}{\Delta x}\right)+u_{i}^{(2)}\left[\left(\frac{x-x_{i}}{\Delta x}\right)^{2}-\frac{1}{12}\right]$, which is the projection of $p_{i}^{3}(x)$ in the quadratic orthogonal function space $\left\{1,\frac{x-x_{i}}{\Delta x},\left(\frac{x-x_{i}}{\Delta x}\right)^{2}-\frac{1}{12}\right\}$. Dropping the cubic term doesn’t affect the accuracy for $u_{h}$ is used only in the troubled-cell indicator (2.3) to identify troubled cell, we use next high order HWENO methodology to reconstruct the first order moment in the troubled cell.

**Step 1.2. Modify the first order moment in the troubled-cell. **

If the cell $I_{i}$ is identified as a troubled cell, we would modify the first order moment $\overline{v}_{i}$. The procedure to modify the first order moment is the same as that HWENO limiter [28]. At first, we give three small stencils $S_{1}=\{I_{i-1},I_{i}\}$, $S_{2}=\{I_{i-1},I_{i},I_{i+1}\}$, $S_{3}=\{I_{i},I_{i+1}\}$, and a large stencil $S_{0}=\{S_{1},S_{2},S_{3}\}=S_{2}$, then, we obtain three quadratic Hermite polynomials $p_{1}(x),p_{2}(x),p_{3}(x)$ on $S_{1},S_{2},S_{3}$, respectively, as

[TABLE]

and get a quartic polynomial $p_{0}(x)$ on $S_{0}$, satisfying

[TABLE]

Then, we use these polynomials to reconstruct $\overline{v}_{i}$, and their explicit results based on the moments $\{\overline{u}_{i},\overline{v}_{i}\}_{i}$ are

[TABLE]

The linear weights $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ can be obtained easily, just following as

[TABLE]

which leads to $\gamma_{1}=\frac{11}{38}$, $\gamma_{2}=\frac{8}{19}$ and $\gamma_{3}=\frac{11}{38}$, then, we compute the smoothness indicators $\beta_{n}$, which measure how smooth the functions $p_{n}(x)$ in the target cell $I_{i}$, and we use the same definition for the smoothness indicators as in [18, 34],

[TABLE]

Here, $r=2$ is the degree of the polynomials $p_{n}(x)$, and their explicit expressions are shown as

[TABLE]

Then the nonlinear weights are computed by:

[TABLE]

where $\varepsilon$ is a small positive number to avoid the denominator by zero, and we take $\varepsilon=10^{-6}$ in this paper. Hence, the first order moment $\overline{v}_{i}$ is finally modified by

[TABLE]

Step 2. The reconstruction procedure for Gauss-Lobatto points values.

In this subsection, we would give the details of reconstruction procedure for Gauss-Lobatto points values $u^{\pm}_{i\mp 1/2}$ and $u_{i\pm\sqrt{5}/10}$ from $\{\overline{u}_{i},\overline{v}_{i}\}_{i}$. Similarly to Step 1, we first give the stencils $S_{1},S_{2},S_{3}$ and $S_{0}$. If one of the cells in stencil $S_{0}$ is identified as a troubled cell, we would use the HWENO method described in Step 2.1 to reconstruct the $u^{\pm}_{i\mp 1/2}$; otherwise we use the upwind linear approximation method described in Step 2.2 to reconstruct the $u^{\pm}_{i\mp 1/2}$. And the reconstruction of $u_{i\pm\sqrt{5}/10}$ is described in Step 2.3.

Step 2.1. The HWENO reconstruction for $u^{-}_{i+1/2}$.

If one of the cells in stencil $S_{0}$ is identified as a troubled cell, $u^{-}_{i+1/2}$ is reconstructed by HWENO procedure. We reconstruct three cubic polynomials $p_{1}(x),p_{2}(x),p_{3}(x)$ on $S_{1},S_{2},S_{3}$, respectively, such that:

[TABLE]

and reconstruct a quintic polynomial $p_{0}(x)$ on $S_{0}$, as

[TABLE]

Based on (2.8) and (2.9), we compute the approximations of $u^{-}_{i+1/2}$ by these polynomials at the point $x_{i+1/2}$, and their explicit expressions are

[TABLE]

Then we get the linear weights $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$, according to

[TABLE]

and we have $\gamma_{1}=\frac{25}{189}$, $\gamma_{2}=\frac{22}{63}$ and $\gamma_{3}=\frac{14}{27}$. Then, we compute the smoothness indicators $\beta_{n}$, which measure how smooth the functions $p_{n}(x)$ in the cell $I_{i}$. Again, we use the formula (2.6) to compute the smoothness indicators, where $r=3$, and we have the explicit expressions for the smoothness indicators:

[TABLE]

and the nonlinear weights are computed by:

[TABLE]

Here $\varepsilon$ is a small positive number taken as $10^{-6}$. Hence, the final value of $u^{-}_{i+1/2}$ is reconstructed by

[TABLE]

The reconstruction to $u^{+}_{i-1/2}$ is mirror symmetric with respect to $x_{i}$ of the above procedure.

Step 2.2. The linear approximation for $u^{\mp}_{i\pm 1/2}$.

If nether cell in $S_{0}$ is identified as troubled cell, then we will use upwind linear reconstruction for $u^{\mp}_{i\pm 1/2}$, that is we use $p_{0}(x)$ to approximate $u$ directly, and we have:

[TABLE]

and

[TABLE]

Step 2.3. The linear reconstruction for $u_{i\pm\sqrt{5}/10}$.

We would like to use linear reconstruction for $u_{i\pm\sqrt{5}/10}$ in all cells, then, $u_{i\pm\sqrt{5}/10}$ are finally approximated by the following expressions, respectively,

[TABLE]

and

[TABLE]

Step 3. When we have finished the spatial discretization following Steps 1 and 2, the semi-discrete schemes (2.2) are discretized in time by the third order TVD Runge-Kutta method [31]:

[TABLE]

Remark 1: For one dimensional scalar equation, the solution $u$ is taken as the indicator variable and $\overrightarrow{v}$ is defined as $f^{\prime}(u)$; while for one dimensional Euler equations, the density $\rho$ and the energy $E$ are taken as the indicator variables, respectively, and $\overrightarrow{v}$ is set as the velocity $\mu$ of the fluid.

Remark 2: For the systems, such as the one dimensional compressible Euler equations, in order to achieve better qualities at the price of more complicated computations, the HWENO approximation is always used with a local characteristic field decomposition seen in e.g. [33, 34] for details, while the linear approximation is used in component by component.

2.2 Two dimensional case

Similarly to one dimensional case, we first consider two dimensional scalar hyperbolic conservation laws

[TABLE]

For simplicity of presentation, the computing domain is divided by uniform meshes. The mesh sizes are $\Delta x=x_{i+1/2}-x_{i-1/2}$ in the $x$ direction and $\Delta y=y_{j+1/2}-y_{j-1/2}$ in the $y$ direction, and each cell of the mesh $I_{i,j}$ is taken as $[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$ with its cell center $(x_{i},y_{j})=(\frac{x_{i-1/2}+x_{i+1/2}}{2},\frac{y_{j-1/2}+y_{j+1/2}}{2})$. In the next procedures, $x_{i}+a\Delta x$ is defined as $x_{i+a}$, while $y_{j}+b\Delta y$ is set as $y_{j+b}$.

To design a HWENO scheme, we multiply the equation (2.15) by $\frac{1}{\Delta x\Delta y}$, $\frac{x-x_{i}}{(\Delta x)^{2}\Delta y}$ and $\frac{y-y_{j}}{(\Delta y)^{2}\Delta x}$ on both sides, respectively, and we integrate them over $I_{i,j}$, then, apply the integration by parts and employ the numerical flux to approximate the values of the flux at the points on the interface of $I_{i,j}$, lastly, we get the semi-discrete finite volume HWENO scheme, and the explicit formulas are given as follows,

[TABLE]

with initial conditions \overline{u}_{i,j}(0)=$$\frac{1}{\Delta x\Delta y}\int_{I_{i,j}}u_{0}(x,y)dxdy, \overline{v}_{i,j}(0)=$$\frac{1}{\Delta x\Delta y}\int_{I_{i,j}}u_{0}(x,y)\frac{x-x_{i}}{\Delta x}dxdy and \overline{w}_{i,j}(0)=$$\frac{1}{\Delta x\Delta y}\int_{I_{i,j}}u_{0}(x,y)\frac{y-y_{j}}{\Delta y}dxdy. $\overline{u}_{i,j}(t)$ is the cell average set as \frac{1}{\Delta x\Delta y}$$\int_{I_{i,j}}u(x,y,t)dxdy; $\overline{v}_{i,j}(t)$ is the first order moment in the $x$ direction defined as \frac{1}{\Delta x\Delta y}$$\int_{I_{i,j}}u(x,y,t)\frac{x-x_{i}}{\Delta x}dxdy and $\overline{w}_{i,j}(t)$ is the first moment in the $y$ direction taken as \frac{1}{\Delta x\Delta y}$$\int_{I_{i,j}}u(x,y,t)\frac{y-y_{j}}{\Delta y}dxdy. $\hat{f}(u(x_{i+1/2},y))$ is a numerical flux which is the approximation to the values of the numerical flux $f(u)$ at the interface point $(x_{i+1/2},y)$ and $\hat{g}(u(x,y_{j+1/2}))$ is a numerical flux to approximate the values of $g(u)$ at the interface point $(x,y_{j+1/2})$.

Just as in the one dimensional case, we will approximate the integral terms of equations (2.16) by employing numerical integration. Since we construct a fourth-order accuracy scheme, 2-point Gaussian will be used in each numerical quadrature, then, these approximated formulas for the integral terms are given as follows,

[TABLE]

[TABLE]

[TABLE]

and the approximated expressions for other integral terms are similar. Here, $\omega_{1}=\frac{1}{2}$ and $\omega_{2}=\frac{1}{2}$ are the quadrature weights, and the coordinates of the Gaussian points over the cell $I_{i,j}$ are

[TABLE]

In the two dimensional case, two fluxes in the $x$ direction and in the $y$ direction are approximated by Lax-Friedrichs numerical flux:

[TABLE]

and

[TABLE]

where $\alpha=\max_{u}|f^{\prime}(u)|$, $\beta=\max_{u}|g^{\prime}(u)|$, and $G_{b}$ is the Gaussian point on the interface of the cell $I_{i,j}$.

The general frameworks for the hybrid HWENO scheme are: in Steps 4 and 5, we present the spatial reconstruction for the semi-discrete scheme (2.16). In Step 6, the equations (2.16) is discretized by the third order TVD Runge-Kutta method [31] in time.

Step 4. Identify the troubled-cell and modify the first order moments in the troubled-cell.

We also use the KXRCF troubled-cell indicator by Krivodonova et al. [19] (KXRCF) to identify the discontinuities, which has been introduced in the one dimensional problems, and its explicit expression can be seen in (2.3). In particular, the troubled-cell indicator works separately in the $x$ and $y$ directions for two dimensional case, and the norm is still $L_{\infty}$ norm in the two dimensional case. The cell $I_{i,j}$ is finally identified as a troubled-cell, if it is identified either in $x$ direction or $y$ direction. In addition, if the cell $I_{i,j}$ is identified as a troubled cell, we mark the cells $I_{i-1,j}$, $I_{i+1,j}$, $I_{i,j-1}$ and $I_{i,j+1}$ as troubled cells in practice, as the spatial reconstruction for the neighbor cells also need to use the information of $I_{i,j}$.

If the cell $I_{i,j}$ is identified as a troubled cell, we would modify the first order moments $\overline{v}_{i,j}$ and $\overline{w}_{i,j}$ following as next procedures. We modify the first order moments in the troubled cells using dimensional by dimensional manner, and the modification procedures are the same as the one dimensional case. More explicitly, if the cell $I_{i,j}$ is identified as a troubled cell, we use these information $\overline{u}_{i-1,j}$, $\overline{u}_{i,j}$, $\overline{u}_{i+1,j}$, $\overline{v}_{i-1,j}$, $\overline{v}_{i+1,j}$ to reconstruct the value of $\overline{v}_{i,j}$, and the procedures are the same as the expressions introduced in Step 1 for one dimensional case, and the procedures for the modification of $\overline{w}_{i,j}$ are similar.

Step 5. Reconstruct the point values of the solutions $u$ at the specific points.

This step is to reconstruct the point values of $u^{+}(x_{i-1/2},y_{j\pm\sqrt{3}/6})$, $u^{-}(x_{i+1/2},y_{j\pm\sqrt{3}/6})$, $u^{+}(x_{i\pm\sqrt{3}/6},y_{j-1/2})$, $u^{-}(x_{i\pm\sqrt{3}/6},y_{j+1/2})$ and $u(x_{i\pm\sqrt{3}/6},y_{j\pm\sqrt{3}/6})$ in the cell $I_{i,j}$. If the cell is identified as a troubled cell, in Step 5.1, the interface points values of the cell $I_{i,j}$ are reconstructed by HWENO methodology but the internal points values of $I_{i,j}$ are approximated by linear approximation, respectively; otherwise, we directly use the linear approximation presented in Step 5.2.

Step 5.1. Reconstruct the point values of the solutions $u$ at the interface points by HWENO approximation and approximate the internal points values using linear approximation.

We first give the big stencil $S_{0}$ in Figure 2.1. For simplicity, we rebel the cell $I_{i,j}$ and its neighboring cells as $I_{1},...,I_{9}$. Particularly, the new label of the cell $I_{i,j}$ is $I_{5}$ and the symbols ”*” in $I_{5}$ represent the locations of the solutions $u$ where we need to reconstruct. We also give eight small stencils: $S_{1},S_{2},...,S_{8}$ shown in Figure 2.2, which were first introduced by Qiu and Shu in [29]. Noticed that we reconstruct the point values of solutions $u$ in the cell $I_{i,j}$, then, we would like to use more information in the cell $I_{i,j}$, such as the first order moments $\overline{v}_{i,j}$ and $\overline{w}_{i,j}$, so we construct eight incomplete cubic reconstruction polynomials, and these polynomials have the expressions as follows

[TABLE]

satisfying as

[TABLE]

for

[TABLE]

Then, we combine these eight incomplete cubic polynomials to obtain a fourth-order approximation for the reconstruction of these points values of the solutions $u$. For simplicity, we use $G_{k}$ to represent the specific points, where we want to reconstruct. At first, we use the linear weights $\gamma_{1}^{k},...,\gamma_{8}^{k}$ to joint these eight small polynomials, satisfying as

[TABLE]

If $\sum_{n=1}^{8}\gamma_{n}^{k}=1$, the requirement (2.18) are always satisfied for any incomplete cubic polynomial $u$, and the form of the polynomial is presented in (2.17), but we still have two other constraints on the linear weights to hold the requirement (2.18) for $u=x^{2}y$ and $xy^{2}$. Subject to these three constraints listed above, it leaves five free parameters to calculate the linear weights, and we can obtain $\gamma_{1}^{k},...,\gamma_{8}^{k}$ easily and uniquely by minimizing $\sum_{n=1}^{8}(\gamma_{n}^{k})^{2}$. In fact, the linear weights $\gamma_{1}^{k},...,\gamma_{8}^{k}$ determined by this least square methodology are all positive in the implementation. For simplicity, we only present the eight linear weights to reconstruct $u^{-}(x_{i+1/2},y_{j+\sqrt{3}/6})$ at the interface point, and the values are $\frac{3533+351\sqrt{3}}{37040}$, $\frac{5727+351\sqrt{3}}{37040}$, $\frac{3533-351\sqrt{3}}{37040}$, $\frac{5727-351\sqrt{3}}{37040}$, $\frac{10599-1867\sqrt{3}}{111120}$, $\frac{17181-415\sqrt{3}}{111120}$, $\frac{10599+1867\sqrt{3}}{111120}$, $\frac{17181+415\sqrt{3}}{111120}$, respectively, and the linear weights for other points on the interface can be obtained by symmetry. In addition, it is interesting that the linear weights to reconstruct $u(x_{i\pm\sqrt{3}/6},y_{j\pm\sqrt{3}/6})$ are all $\frac{1}{8}$.

Similarly as in the one dimensional problems, if $G_{k}$ is inside of $I_{i,j}$, we directly use linear approximation to reconstruct $u(G_{k})$ as $\sum_{n=1}^{8}\gamma_{n}^{k}p_{n}(G_{k})$, and we’d better simplify the formals in advance, instead of calculating in the codes; while $G_{k}$ is located on the interface of the cell $I_{i,j}$, we need to employ the next HWENO reconstruction procedures, then, we first compute the smoothness indicators $\beta_{n}$, which measure how smooth the function $p_{n}(x,y)$ in the cell $I_{i,j}$. The formula was listed by [14], given as follows,

[TABLE]

where $l=(l_{1},l_{2})$, $|l|=l_{1}+l_{2}$, then, we can get the non-linear weights using the linear weights and the smoothness indicators, having

[TABLE]

where $\varepsilon$ is set as $10^{-6}$ just as in one dimensional case. Finally, the approximation for the point values of the solutions $u$ at the interface point $G_{k}$ is given by

[TABLE]

where ”” is ”+” when $G_{k}$ is located on the left or bottom interface of the cell $I_{i,j}$, while ”” is ”-” on the right or top interface of $I_{i,j}$.

Step 5.2. Reconstruct the point values of the solutions $u$ at the specific points by linear approximation straightforwardly.

In this step, we’ll use the same polynomials and linear weights introduced in Step 5.1, then, the linear approximation of the solutions $u$ at reconstructed point $G_{k}$ can be directly taken as

[TABLE]

If $G_{k}$ is located on the interface of the cell $I_{i,j}$, ”” has the same meaning just as in Step 5.1; otherwise, ”” will be blank. Similarly, we also can obtain the simplified formulas easily for the linear approximation of $u^{*}(G_{k})$ in advance, instead of calculating over and over again in the codes.

Noticed that the reconstruction for the points values of the solutions $u$ only has the fourth order accuracy for the information we used here is not enough to reconstruct interpolation polynomial with degree 4, therefore, the scheme only is the fourth order in the two dimensional case.

Step 6. Discretize the semi-discrete scheme (2.16) in time by the third order TVD Runge-Kutta method [31].

When we have finished Steps 4 and 5, the semi-discrete scheme (2.16) is discretized in time by the third order TVD Runge-Kutta method, and the explicit expression has been presented in (2.14) for the one dimensional case.

Remark 3: The KXRCF indicator is satisfying for two dimensional hyperbolic conservation laws. For two dimensional scalar equation, the solution $u$ is set as the indicator variable. $\overrightarrow{v}$ is taken as $f^{\prime}(u)$ in the $x$ direction, while it is defined as $g^{\prime}(u)$ in the $y$ direction; for two dimensional Euler equations, the density $\rho$ and the energy $E$ are taken as the indicator variables, respectively. $\overrightarrow{v}$ is defined as the velocity $\mu$ in the $x$ direction of the fluid, while it is set as the velocity $\nu$ in the $y$ direction of the fluid.

Remark 4: For the systems, such as the two dimensional compressible Euler equations, similarly as in the one dimension case, we first use the KXRCF indicator to identify the troubled cell in Step 4. If the cell $I_{i,j}$ is identified as a troubled cell, we modify the first order moments for each component. For the modification, it is different from the one dimensional case for it has two first order moments and two flux functions, so we modify the moments $\overline{v}_{i,j}$ in the $x$ direction in terms of the local characteristic direction provided by $f(u)$, while reconstruct $\overline{w}_{i,j}$ in the $y$ direction based on the local characteristic direction of $g(u)$. For each local characteristic direction, we follow the procedures of Step 4 to reconstruct the first order moments in the troubled cells. For the reconstruction for the point values of the solutions $u$, all HWENO procedures are performed on the local characteristic decompositions, and the linear approximations are based on component by component.

3 Numerical tests

In this section, we perform the numerical results of the hybrid HWENO scheme in the one and two dimensional cases, which is outlined in the previous section. If no otherwise specified, the CFL number is set as 0.6 for one dimensional tests and 0.45 for two dimensional examples.

3.1 Accuracy tests

For simplicity, Hybrid HWENO scheme is denoted as the hybrid HWENO scheme introduced in the previous section, while HWENO scheme is represented as that we modify the first order moments for every cell and employ HWENO reconstruction at the interface points for the spatial discretization on the basis of the hybrid HWENO scheme. WENO scheme was listed by Jiang and Shu [18], while Hybrid WENO scheme using KXRCF troubled-cell indicator was introduced by Li and Qiu [23]. Since the four schemes all have fifth order accuracy in one dimensional problems, we will compare their performance in the one dimensional accuracy tests. For two dimensional smooth tests, as the Hybrid HWENO scheme is based on the finite volume framework, we also make the comparisons with the classical fifth order finite volume WENO scheme narrated in [33].

Example 3.1. We solve the one dimensional Burgers’ equation:

[TABLE]

with the initial condition $u(x,0)=0.5+sin(\pi x)$ and periodic boundary condition. We present the numerical results at $t=0.5/\pi$ when the solution is still smooth, then, the numerical errors and orders are shown in Table 3.1 with $N$ uniform meshes for HWENO and WENO schemes. From the table, we can see that all four schemes have fifth order accuracy. Firstly, we know the hybrid schemes have less errors than the original schemes, meanwhile, we also can find that two HWENO schemes have less errors than corresponding WENO schemes with the same number of cells. In Figure 3.1, we show numerical errors against CPU times by using four different schemes, which illustrates Hybrid HWENO scheme has much higher efficiency than other three schemes, meanwhile, the two HWENO schemes only need three cells while the two WENO schemes need five cells for the spatial reconstruction.

Example 3.2. We consider one dimensional Euler equations:

[TABLE]

in which $\rho$ is the density, $\mu$ is the velocity, $E$ is the total energy and $p$ is the pressure. The initial conditions are set to $\rho(x,0)=1+0.2\sin(\pi x)$, $\mu(x,0)=1$, $p(x,0)=1$ and $\gamma=1.4$. The computing domain is $x\in[0,2\pi]$. Periodic boundary condition is applied here. The exact solution is $\rho(x,t)=1+0.2\sin(\pi(x-t))$, $\mu(x,0)=1$, $p(x,0)=1$ and the final computing time is $T=2$. The numerical errors and orders of the density for the HWENO and WENO schemes are given in Table 3.2, which shows four schemes achieve the designed fifth order accuracy. Similarly, the hybrid schemes have less errors than the original schemes and the two HWENO schemes have less errors than corresponding WENO schemes with the same number of cells. In addition, Figure 3.2 represents numerical errors against CPU times using four different schemes, which shows Hybrid HWENO scheme has higher efficiency than other three schemes, and the HWENO schemes are more compact than the WENO schemes.

Example 3.3. Consider the two dimensional Burgers’ equation:

[TABLE]

with the initial condition $u(x,y,0)=0.5+sin(\pi(x+y)/2)$ and periodic boundary condition in each direction. We perform the numerical results at $t=0.5/\pi$. In this time, the solution is still smooth, then, the numerical errors and order are shown in Table 3.3 for Hybrid HWENO and WENO schemes. We can see that both schemes achieve the designed order accuracy. In Figure 3.3, we present their numerical errors against CPU times, which illustrates Hybrid HWENO scheme has higher efficiency than WENO scheme, meanwhile, the hybrid HWENO scheme is more compact for only immediate neighbor information is needed in the spatial reconstruction.

Example 3.4. We consider two dimensional Euler equations:

[TABLE]

Here $\rho$ is the density, $(\mu,\nu)$ is the velocity, $E$ is the total energy, and $p$ is the pressure. The initial conditions are taken as $\rho(x,y,0)=1+0.2\sin(\pi(x+y))$, $\mu(x,y,0)=1$, $\nu(x,y,0)=1$, $p(x,y,0)=1$ and $\gamma=1.4$. The computing domain is $(x,y)\in[0,2]\times[0,2]$. Periodic boundary conditions are applied in $x$ and $y$ directions. The exact solution of $\rho$ is $\rho(x,y,t)=1+0.2\sin(\pi(x+y-2t))$, $\mu(x,y,0)=1$, $\nu(x,y,0)=1$, $p(x,y,0)=1$ and the computing time is up to $T=2$. Table 3.4 gives the numerical errors and orders of the density for the hybrid HWENO and WENO schemes, and we can know both two schemes achieve the designed fourth and fifth order accuracy, respectively. In addition, we also present their numerical errors against CPU times in Figure 3.4, which shows Hybrid HWENO scheme has higher efficiency than WENO scheme

Example 3.5. Isentropic vortex test [33] for two dimensional Euler equations. A isentropic vortex is added to the mean flow $(\rho,\mu,\nu,p,\gamma)$ $=$ $(1,1,1,1,1.4)$ with perturbations in ($\mu,\nu$), the temperature $T=p/\rho$ and no perturbation in the entropy $S=p/\rho^{\gamma}$, so that the initial conditions are

[TABLE]

in which $\epsilon$ represents the vortex strength taken as 5 here and $r^{2}=x^{2}+y^{2}$. The computational domain is $[-5,5]\times[-5,5]$, and periodic boundary conditions are applied in $x$ and $y$ directions. The vortex movement is aligned with the diagonal of the computational domain, and it returns to the initial position with time period 10. we present the numerical errors and orders of the density for the hybrid HWENO and WENO schemes in Table 3.5, then we can see both two schemes achieve the designed accuracy, respectively. In addition, we also give their numerical errors against CPU times in Figure 3.5, which illustrates Hybrid HWENO scheme has higher efficiency than WENO scheme with smaller numerical errors and less CPU times.

3.2 Non-smooth tests

We present the results of the hybrid HWENO scheme for the non-smooth tests. To make comparisons with WENO scheme, we also give the numerical results for some of one dimensional non-smooth tests. In addition, we have computed these tests for HWENO scheme likewise, but the results will not be presented here for saving space. Actually, they have similar performance for the problems with discontinuities, but HWENO scheme uses much more computing time. Moreover, we also test the non-smooth problems by the method that we don’t modify the first order moments of any cells and use HWENO reconstruction at the interface points of each cell for the spatial discretization. Unfortunately, all non-smooth tests have obvious oscillations near discontinuities, even for the one-dimensional Burgers’ equation with smaller CFL number, and some tests don’t work in this case even though using smaller CFL number, such as Shu-Osher and two blast waves problems, and so on, which illustrate that the modification for these first order moments in the troubled cells for the hybrid HWENO scheme is significant to avoid oscillations and keep the hybrid HWENO scheme be robust. In the implementation, there are only a small part of cells in which we need to modify their first order moments for the hybrid HWENO scheme.

Example 3.6. We solve the one-dimensional Burgers’ equation (3.1) as in Example 3.1. The same initial and boundary conditions are applied here. The computing time is up to $t=1.5/\pi$, and the solution is discontinuous by this time, then, we plot the numerical solution of the hybrid HWENO scheme and the exact solution in Figure 3.6, and we also test this example by the method that we don’t modify any first order moments and directly use HWENO reconstruction at the interface points of each cell for the spatial discretization. We set the CFL number as 0.45 in this case for the code doesn’t work with the original CFL number, and its numerical solution is also presented in Figure 3.6. From this figure, we can see that if we don’t modify the first order moments in the troubled cell and directly use HWENO procedures at the interface points of each cell for spatial discretization, which would have obvious oscillations even with smaller CFL number. However, we also find that the hybrid HWENO scheme can avoid the non-physical oscillations and has a good resolution, which shows that the modification for the first order moments near the discontinuities is a significant and essential procedure. In addition, there are only 4.52% cells where we need to modify their first order moments by calculating.

Example 3.7. We now consider a one dimensional nonlinear non-convex scalar Buckley-Leverett problem

[TABLE]

with the initial condition: $u=1$ for $-\frac{1}{2}\leq x\leq 0$ and $u=0$ elsewhere. Inflow and outflow conditions are applied at left and right boundary, respectively, and the computing time is up to $t=0.4$. The exact solution contains both shock and rarefaction, moreover, some high-order schemes may non-converge to the right entropy solution for this test. We present the numerical results in Figure 3.7, and we can find that the hybrid HWENO scheme performs well for capturing the correct entropy solution and has a good resolution.

Example 3.8. We solve the 1D Euler equations Riemann initial condition for the Lax problem

[TABLE]

The final computing time is up to $T=0.16$, and we first present the performances of the exact solution and the computed density $\rho$ obtained with the hybrid HWENO and WENO schemes by using 200 uniform cells in Figure 3.8. The zoomed in picture and the time history of the cells where we modify the first order moments in the hybrid HWENO scheme are also given in Figure 3.8. In this test case, there are about 10.71% cells in which we modify their first order moments, which means most regions use linear approximation and have non-decomposition in the characteristic direction, therefore, the hybrid HWENO scheme saves about 62.5% computational time than the HWENO scheme, meanwhile, the modification in the troubled cells is very important to control the oscillations. The hybrid HWENO scheme also keeps good resolution, and the hybrid HWENO and WENO schemes have similar numerical performances, but the hybrid HWENO scheme only needs the immediate neighbor information.

Example 3.9. The Shu-Osher problem, which describes shock interaction with entropy waves [34], and the initial condition is

[TABLE]

As we know, when the solutions both contains shocks and complex smooth region structures, a good scheme would simulate it well. Actually, this test case is a typical example with a moving Mach=3 shock interacting with sine waves in density. We solve this problem up to $T=1.8$. In Figure 3.9, we present the computed density $\rho$ by the hybrid HWENO scheme and WENO schemes against the referenced ”exact” solution, the zoomed in picture and the time history of the cells where we modify the first order moments in the hybrid HWENO scheme. The referenced ”exact” solution is computed by the fifth order finite difference WENO scheme [18] with 2000 grid points. We find that there are only 2.42% cells where we need to modify their first order moments by calculating, which saves near 66.7% CPU time than HWENO scheme, but the modification in the troubled cells is very significant to make the scheme be robust. In addition, we also see that the hybrid HWENO scheme has higher resolution than WENO scheme shown in Figure 3.9. Particularly, the hybrid HWENO scheme only needs three cells while the WENO scheme needs five cells for the spatial reconstruction.

Example 3.10. We now consider the interaction of two blast waves, and the initial conditions are:

[TABLE]

The final computing time $T=0.038$ and the reflective condition is applied in two boundaries. In Figure 3.10, we present the computed density by the hybrid HWENO and WENO schemes against the reference ”exact” solution, the zoomed in picture and the time history of the cells where we modify the first order moments in the hybrid HWENO scheme. Particularly, the reference ”exact” solution is a converged solution computed by the fifth order finite difference WENO scheme [18] with 2000 grid points. In the implementation, we find that the hybrid HWENO scheme saves about 58.5% computational time as there are almost 11.31% cells where we need to modify their first order moments. The modification for the first order moment in the troubled cells makes the hybrid HWENO scheme be robust, meanwhile, the hybrid HWENO scheme has higher resolution than WENO scheme.

Example 3.11. We consider the two-dimensional Burgers’ equation (3.3) seen in Example 3.3 with the same initial and boundary conditions here. The finial computing time is $t=1.5/\pi$, and the solution is discontinuous, then, we plot the numerical solution and the exact solution in Figure 3.11 with $80\times 80$ uniform meshes. Similarly as in one dimensional Burgers’ equation, the hybrid HWENO scheme has good resolution nearby discontinuities.

Example 3.12. Double Mach reflection problem [37], which is modeled by the two-dimensional Euler equations (3.4). The computational domain is $[0,4]\times[0,1]$. A reflection wall lies at the bottom, starting from $x=\frac{1}{6}$, $y=0$, making a $60^{o}$ angle with the $x$-axis. For the bottom boundary, the exact post-shock condition is imposed for the part from $x=0$ to $x=\frac{1}{6}$ and the rest is reflection boundary condition, while it is the exact motion of the Mach 10 shock for the top boundary. $\gamma=1.4$ and the final computing time is set as $T=0.2$. In Figure 3.12, we present the numerical results in region $[0,3]\times[0,1]$, the cells where we modify the first order moments in the hybrid HWENO scheme at the final time and the blow-up region around the double Mach stems. Again, the hybrid HWENO scheme works well for this test case, meanwhile, there are only 3.55% cells in which we need to modify their first order moments, which saves about 68.2% computational time as most regions directly use linear approximation, and it shows the hybrid HWENO scheme has higher efficiency than HWENO scheme.

Example 3.13. The final example is about a Mach 3 wind tunnel with a step modeled by the two-dimensional Euler equations (3.4), which is also originally from [37]. The setup of the problem is as follows. The wind tunnel is 1 length unit wide and 3 length units long. The step is 0.2 length units high and is located 0.6 length units from the left-hand end of the tunnel. The problem is initialized by a right-going Mach 3 flow. Reflective boundary conditions are applied along the wall of the tunnel and inflow/outflow boundary conditions are applied at the entrance/exit. The corner of the step is a singular point and we treat it as the same way in [37], which is based on the assumption of a nearly steady flow in the region near the corner. The final time is $T=4$. In Figure 3.13, we present the results for the hybrid HWENO scheme with $960\times 320$ uniform cells and the cells where we modify the first order moments in the hybrid HWENO scheme at the last time step. We can notice that the hybrid HWENO scheme gets good resolutions in the non-smooth regions, moreover, we find that there are only 11.68% cells where we need to modify their first order moments in our implementation, which shows most regions directly use high order linear approximation and we don’t need to modify their first order moments, and we also know the hybrid HWENO scheme has higher efficiency than the HWENO scheme for saving near 64.2% CPU time by calculating.

4 Concluding remarks

In this paper, a hybrid finite volume Hermite weighted essentially non-oscillatory (HWENO) scheme is designed for solving hyperbolic conservation laws. Compared with other HWENO methods [28, 29, 39, 7, 16, 35, 24, 43, 36, 8, 25] for hyperbolic conservation laws, we bring the idea of limiter for discontinuous Galerkin (DG) method to modify the first order moments in the troubled cells, meanwhile, we find that there are only a small part of cells where we need to modify their first order moments, and we also reconstruct the interface point values for the solutions using HWENO approximation in these troubled cells. In other words, there are many cells in which we reconstruct the point values for the solutions directly employing linear approximation, which makes the hybrid HWENO scheme be higher efficiency. In general, the modification for the first order moments in the troubled cells is significant to avoid oscillations and keep the scheme be robust for non-smooth numerical tests, and these smooth and non-smooth numerical results all illustrate the good performances of the hybrid HWENO scheme.