Finite volume HWENO schemes for nonconvex conservation laws1footnote 11footnote 1Research was supported by NSFC grants 11571290, 91530107, Air Force Office of Scientific Research FA9550-16-1-0179 and NSF DMS-1522777.

# Finite volume HWENO schemes for nonconvex conservation laws111Research was supported by NSFC grants 11571290, 91530107, Air Force Office of Scientific Research FA9550-16-1-0179 and NSF DMS-1522777.

Xiaofeng Cai222Department of Mathematics, University of Houston, Houston, 77204. E-mail: xfcai@math.uh.edu., Jianxian Qiu333School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modeling & High-Performance Scientific Computing, Xiamen University, Xiamen, Fujian, 361005, P.R. China. E-mail: jxqiu@xmu.edu.cn., and Jingmei Qiu444Department of Mathematics, University of Houston, Houston, 77204. E-mail: jingqiu@math.uh.edu.

Abstract: Following the previous work of Qiu and Shu [SIAM J. Sci. Comput., 31 (2008), 584-607], we investigate the performance of Hermite weighted essentially non-oscillatory (HWENO) scheme for nonconvex conservation laws. Similar to many other high order methods, we show that the finite volume HWENO scheme performs poorly for some nonconvex conservation laws. We modify the scheme around the nonconvex regions, based on a first order monotone scheme and a second entropic projection, to ensure entropic convergence. Extensive numerical tests are performed. Compare with the earlier work of Qiu and Shu which focuses on 1D scalar problems, we apply the modified schemes (both WENO and HWENO) to one-dimensional Euler system with nonconvex equation of state and two-dimensional problems.

Keywords: Nonconvex conservation laws, Finite volume HWENO scheme, Entropy solution, Entropic projection.

## 1 Introduction

In this paper, we consider the Cauchy problem for nonconvex hyperbolic conservation laws:

 {qt+∇⋅F(q)=0,q(x,0)=q0(x), (1.1)

whose entropy solutions may admit composite waves which involve a sequence of shocks and rarefaction waves and are difficult to be resolved numerically. Such examples include scalar conservation laws with nonconvex flux functions and hyperbolic systems such as the Euler system and magnetohydrodynamics system with a nonconvex equation of state (EOS) [6, 17, 7, 12].

It is well known that first order monotone schemes converge to entropy solutions of both convex and nonconvex conservation laws [3], but with a relatively slow convergence rate. It has also been known [5, 11] that there are some nonconvex conservation laws, for which high order schemes such as the ones with weighted essentially non-oscillatory (WENO) reconstructions [14] and discontinuous Galerkin methods [2] would fail to converge to the entropy solution. There have been great research effort in ensuring entropic convergence for general nonlinear conservation laws, for example by adding entropy viscosity [4] and by modifying reconstruction operators. Examples for the latter approach include the computationally inexpensive strategy proposed in [5] on an adaptive choice between a low order dissipative reconstruction and a high order central WENO scheme, as well as low order modifications around nonconvex regions to ensure entropic convergence proposed in [11]. Compare with the work in [5], with more computational effort, the second order entropic convergence of the schemes can be rigorously proved [1].

This paper is a natural extension of our earlier work in [11]. We investigate the performance of the finite volume Hermite WENO (HWENO) scheme for nonconvex conservation laws and apply the corresponding modifications as being done in [11]. In addition to the scalar examples discussed in [11], we investigate the performance of modified WENO and HWENO scheme for 2D problems, as suggested in [5]. The FV HWENO scheme was originally proposed in [8, 9]. The key idea of the scheme is to evolve more pieces of information, i.e. functions and their spatial gradients, per computational cell. With such mechanism, the HWENO scheme has relatively compact stencils, hence it is easier to handle boundary conditions compared with the traditional WENO scheme [14]. Moreover, with the same formal accuracy, compact stencils are known to exhibit better resolution of small scale structures by improving dispersive and dissipative properties.

An outline of this paper is as follows. Section 2 describes the high order FV HWENO scheme. In Section 3, FV HWENO schemes with a first order monotone modification and a second order modification using an entropic projection around nonconvex regions are proposed for nonconvex conservation laws. In Section 4, numerical examples are shown to demonstrate the effectively of proposed schemes. Concluding remarks are given in Section 5.

## 2 Description of FV HWENO schemes

We briefly review the FV HWENO scheme for solving conservation laws [8, 9, 20]. The idea of HWENO method is to numerically evolve both the function and its spatial gradients, and use these information in the reconstruction process. Thus it leads to a more compact reconstruction stencil compared with the traditional WENO scheme [13, 14].

General scheme formulation of HWENO. Taking the gradient with respect to spatial variables in (1.1), we obtain the evolution equation for function’s gradients,

 (∇q)Tt+∇⋅(∇⊗F(q))=0, (2.1)

where is a tensor product. The FV HWENO scheme is defined for the equations:

 Ut+∇⋅F(U)=0, (2.2)

where and . We integrate the system (2.2) on a control volume , which is an interval for 1D cases or a rectangle for 2D cases. The integral form of (2.2) reads,

 ddt¯¯¯¯¯UΩk=−1|Ωk|∫∂ΩkF(U)⋅nds (2.3)

where is the volume of and represents the outward unit normal vector to . The line integral in (2.3) can be approximated by a -point Gaussian quadrature on each side of :

 ∫∂ΩkF(U)⋅nds≈S∑s=1|∂Ωks|L∑l=1ωlF(U(Gsl,t))⋅n, (2.4)

where and are Gaussian quadrature points on and weights respectively. is evaluated by a numerical flux (approximate or exact Riemann solvers). We adopt the Lax-Friedrichs flux in this paper, which is given by

 F(U(Gsl,t))⋅n≈12[F(U−(Gsl,t))+F(U+(Gsl,t))]⋅n−α(U+(Gsl,t)−U−(Gsl,t)),

where is taken as an upper bound for eigenvalues of the Jacobian along the direction , and and are the reconstructed values of at Gaussian point inside and outside . Finally, the semi-discretization HWENO scheme (2.3) can be written in the following ODE form:

 ddt¯¯¯¯¯U=L(¯¯¯¯¯U). (2.5)

The ODE system (2.5) is then discretized in time by a strong stability preserving Runge-Kutta (RK) method in [15]. The following third-order version is used in this paper,

 ¯¯¯¯¯U(1)=¯¯¯¯¯Un+ΔtL(¯¯¯¯¯Un),¯¯¯¯¯U(2)=34¯¯¯¯¯Un+14(¯¯¯¯¯U(1)+ΔtL(¯¯¯¯¯U(1))),¯¯¯¯¯Un+1=13¯¯¯¯¯Un+23(¯¯¯¯¯U(2)+ΔtL(¯¯¯¯¯U(2))). (2.6)

A scalar 1D example. As an example, we consider a scalar 1D equation,

 qt+f(q)x=0. (2.7)

Taking the derivative of (2.7), we obtain the equation for the derivative,

 ξt+H(q,ξ)x=0, (2.8)

where and Let and denote approximation to cell averages of and over cell respectively, the semi-discrete FV HWENO scheme is designed by approximating spatial derivatives in equation (2.7) and (2.8) with the following flux difference form,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩d¯qjdt=−ˆf⎛⎝q−j+12,q+j+12⎞⎠−ˆf⎛⎝q−j−12,q+j−12⎞⎠Δx,d¯ξjdt=−ˆH⎛⎝q−j+12,q+j+12;ξ−j+12,ξ+j+12⎞⎠−ˆH⎛⎝q−j−12,q+j−12;ξ−j−12,ξ+j−12⎞⎠Δx, (2.9)

where and are reconstructed with high order from neighboring cell averages and . The details of such reconstruction procedures can be found in [8]. is a monotone numerical flux (non-decreasing in the first argument and non-increasing in the second argument), and is non-decreasing in the third argument and non-increasing in the fourth argument. In this paper, we use the Lax-Friedrichs flux in [8],

 ˆf(a,b)=12[f(a)+f(b)−α(b−a)],ˆH(a,b;c,d)=12[H(a,c)+H(b,d)−α(d−c)], (2.10)

where . For the first order “building block” of the HWENO scheme with the Lax-Friedrichs flux, the total variation stability is proved in [8].

A scalar 2D example. We consider a 2D problem on a rectangular domain :

 qt+f(q)x+g(q)y=0. (2.11)

We consider a set of uniform mesh with ,

 a=x12
 c=y12

Let , and , , be cell averages. Taking spatial derivatives of (2.11), we obtain

 ξt+H(q,ξ)x+R(q,ξ)y=0, (2.12) ηt+K(q,η)x+S(q,η)y=0, (2.13)

where . A semi-discrete FV HWENO discretization is given by

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ddt¯¯¯qij=−1Δx(ˆfi+12,j−ˆfi−12,j)−1Δy(ˆgi,j+12−ˆgi,j−12),ddt¯¯¯ξij=−1Δx(ˆHi+12,j−ˆHi−12,j)−1Δy(ˆRi,j+12−ˆRi,j−12),ddt¯¯¯ηij=−1Δx(ˆKi+12,j−ˆKi−12,j)−1Δy(ˆSi,j+12−ˆSi,j−12). (2.14)

We define

 ˆfi+12,j=1Δy∫yj+12yj−12f(q(xi+12,y))dy≈L∑igωigˆf(q−i+12,ig,q+i+12,ig), (2.15) ˆHi+12,j=1Δy∫yj+12yj−12f′(q(xi+12,y))ξdy≈L∑igωigˆH(q−i+12,ig,q+i+12,ig,ξ−i+12,ig,ξ+i+12,ig), (2.16) ˆKi+12,j=1Δy∫yj+12yj−12f′(q(xi+12,y))ηdy≈L∑igωigˆK(q−i+12,ig,q+i+12,ig,η−i+12,ig,η+i+12,ig), (2.17)

as the average of fluxes over the right boundary of cell , where the integrations are evaluated by applying the -point Gauss quadrature. The flux functions , and are taken as the Lax-Friedrich flux as in the 1D, see eq. (2.10), and are reconstructed with high order by neighboring cell averages of , and . In the next paragraph, we briefly describe such reconstruction procedure. , and in eq. (2.14) are evaluated in a similar fashion as the average of fluxes over the top boundary of a cell.

We only review the fourth order reconstruction in 2D and refer to [20] for more details. We relabel the cell and its neighboring cells as as shown in Figure 2.1, where is relabeled as . We construct the quadratic polynomials in the following stencils, , to approximate . For instance, a quadratic polynomial can be reconstructed based on the information in the stencil . Such reconstruction will reconstruct a quadratic polynomial on . Similar reconstructions can be done for stencil , and . For stencil to , only cell averages are used in the reconstruction process. We remark that other combination of information are possible for reconstructing 2D quadratic polynomial. The one we just mentioned seems to be very robust and is implemented in our numerical experiments. If we choose the linear weights denoted by such that

 q(Gl)=8∑n=1γ(l)npn(Gl) (2.18)

is valid for any polynomial of degree at most 3, leading to a fourth-order approximation of at the point for all sufficiently smooth functions . Notice that (2.18) holds for any polynomial of degree at most 2 if . There are four additional constraints on the linear weights so that (2.18) holds for and . The rest of free parameters are determined by a least square procedure to minimize .

As for the derivatives (e.g. ), a third-order approximation in each stencil is enough to obtain the fourth-order approximation to . For instance, a cubic polynomial on can be reconstructed based on the information in the stencil . Similar reconstructions can be done for stencil , and . The information in the stencil is adopted to approximate . Finally, can be chosen. The nonlinear weights of 2D HWENO reconstruction can be designed by following the way of the WENO method.

## 3 The modified FV HWENO schemes for nonconvex conservation laws

Although FV HWENO schemes can be successfully applied in many applications [8, 9, 20, 10, 19], they perform poorly for some nonconvex conservation laws as shown below. To remedy this, we propose a first order monotone modification and a second order modification with an entropic projection around nonconvex regions.

### 3.1 An example of nonconvex conservation laws with poor performance for the FV HWENO scheme

We first show a nonconvex conservation law, for which the FV HWENO scheme performs poorly in converging to the entropy solution. We consider the scalar equation (2.7) with the nonconvex flux and the initial condition

 q0(x)={π/64,if x<0,255π/64,if x≥0. (3.1)

It is shown in Figure 3.2, that the numerical solution of the high order FV HWENO scheme does not converge to the entropy solution (solid black lines given by the first order Godunov scheme with a very refined mesh). One of the rarefaction waves in the compound wave is missing.

### 3.2 First order monotone modification

In this subsection, we propose a first order modification to the FV HWENO scheme for 1D nonconvex conservation laws following a similar idea in [11]. The scheme can be summarized as follows, after a suitable initialization to obtain and .

1. Perform the HWENO reconstruction [8].

At each cell interface, say , reconstruct the point values and derivative values using neighboring cell average and respectively by the fifth order HWENO reconstruction procedure in Section 2.

2. Identify the troubled cell boundary .

Criterion I: A cell boundary is good, if , and all fall into the same linear, convex or concave region of the flux function . Otherwise, it is defined to be a troubled cell boundary.

3. At troubled cell boundaries, modify the numerical flux and with a discontinuity indicator in [18]. Specifically, the discontinuity indicator is defined as

 ϕj=βjβj+γj (3.2)

where

 αj=|¯¯¯qj−1−¯¯¯qj|2+ε, τj=|¯¯¯qj+1−¯¯¯qj−1|2+ε, βj=τjαj−1+τjαj+2, γj=(qmax−qmin)2αj.

Here is a small positive number taken as in the code, and and are the maximum and minimum values of over all cells. The discontinuity indicator has the property that

• is on the order of in smooth regions.

• is close to near a strong discontinuity.

Let and where

 qm,−j+12=(1−ϕ2j)q−j+12+ϕ2j¯¯¯qj,qm,+j+12=(1−ϕ2j)q+j+12+ϕ2j¯¯¯qj+1, (3.3)
 ξm,−j+12=(1−ϕ2j)ξ−j+12+ϕ2j¯¯¯ξj,ξm,+j+12=(1−ϕ2j)ξ+j+12+ϕ2j¯¯¯ξj+1, (3.4)

with defined by (3.2), if is a troubled cell boundary. Otherwise, at good cell boundaries, and .

4. Evolve the cell averages and by (2.9).

###### Remark 1.

When a troubled cell boundary is at a strong discontinuity, , hence , and indicating a first order monotone scheme is taking effect around a nonconvex discontinuous region. When a troubled cell boundary is in a smooth region, the modification is obtained with the magnitude at most of the size

 ϕ2jmax(∣∣qj−q−j+12∣∣,∣∣qj+1−q+j+12∣∣)∼O(Δx5),
 ϕ2jmax(∣∣ξj−ξ−j+12∣∣,∣∣ξj+1−ξ+j+12∣∣)∼O(Δx5),

hence it does not affect the fifth order accuracy of the scheme.

It is natural to extend the above first order modification to 2D problems. With the system, the HWENO reconstructions are performed in local characteristics directions [8], then a first order monotone modification in the form of (3.3) and (3.4) is applied. For 2D problems, we identify trouble cell boundaries, that is to check if the convexity fails, via Gaussian points along cell boundaries, e.g. in equation (2.15). Similarly the first order modification is performed with respect to these Gaussian points along cell boundaries.

### 3.3 Second order modification with an entropic projection

A MUSCL type method with an entropic projection is proposed in [1]. It is proved in the same paper that schemes with such entropic projection enjoy cell entropy inequalities for all convex entropy functions. In the following, we apply such entropic projection as a modification to the HWENO scheme around nonconvex regions to ensure entropic convergence.

#### 3.3.1 Review of the MUSCL method satisfying all the numerical entropy inequalities.

The procedure of the MUSCL scheme satisfying all entropy inequalities can be summarized below. Let the numerical solution at time level be written as with over the cell . Initially, , with where the minmod function is defined as follows,

 minmod(a,b)=⎧⎨⎩0,if ab<0,min(a,b), if a,b≥0,max(a,b), if a,b≤0.

It consists of two steps to evolve from to .

1. Exact evolution : Evolve (2.7) exactly for a time step , to obtain a solution , which in general is not a piecewise linear function anymore.

2. An entropic projection : Find a second order approximation to by a piecewise linear function , satisfying

 ∫IjU(qn+1(x))dx≤∫IjU(˜qn+1(x))dx, ∀j (3.5)

for all convex entropy function . Second order reconstruction satisfying (3.5) can be obtain by setting the cell average as

 ¯¯¯qn+1j=1Δxj∫Ij˜qn+1(x)dx (3.6)

and the slope as

 sn+1j=D˜qn+1|Ij=minmodIjζ(y) (3.7)

where

 ζ(y)=2Δxj⎛⎝1xj+12−y∫xj+12y˜qn+1(x)dx−1y−xj−12∫yxj−12˜qn+1(x)dx⎞⎠ (3.8)

The minmod function of on the interval is defined as

 minmod(a,b)g(x)=⎧⎪ ⎪⎨⎪ ⎪⎩0,if ∃y1,y2∈(a,b), s.t. g(y1)g(y2)≤0,min(a,b)g(y),if g(y)>0, ∀y∈(a,b),max(a,b)g(y),if g(y)<0, ∀y∈(a,b). (3.9)

In summary, the scheme can be written out in the following abstract form

 qn+1=P1∘TΔt(qn)≐Q1(Δt)(qn). (3.10)

It enjoys the following convergence theorem as proved in [1].

###### Theorem 1.

[1]. Let , be the exact entropy solution to (2.7) with the initial data then there exists a constant , such that

 ∥Q1(Δt)nq0−q(⋅,T)∥L1≤C(f∞√TΔt+Δx√T/Δt). (3.11)

The second order MUSCL scheme with the entropic projection (3.10) converges to the unique entropy solution, when .

#### 3.3.2 Second order modification to the fifth order FV HWENO schemes

At each time step evolution, , , over the cell is updated in each RK stage. For instance, at the initial stage, is obtain by the HWENO reconstruction from and . and refer to approximations with entropic projection to the left and right boundaries of when the cell is detected as a trouble cell. Initially, and come from a MUSCL scheme with a minmod reconstruction. To show the idea of second order modification to the fifth order FV HWENO schemes, we present the procedure to update from .

Step 1.

Compute by performing the HWENO reconstruction from .

Step 2.

Update and :

1. Identify the troubled cell boundaries, for which we refer to Criterion I in section 3.2 for the details.

2. Only at trouble cell boundaries, modify numerical fluxes and as follows. We let and , where

 qm,−j+12=(1−ϕ2j)qn,−j+12+ϕ2jqn,rj,qm,+j+12=(1−ϕ2j)qn,+j+12+ϕ2jqn,lj+1, (3.12)
 ξm,−j+12=(1−ϕ2j)ξn,−j+12+ϕ2j¯¯¯ξnj,ξm,+j+12=(1−ϕ2j)ξn,+j+12+ϕ2j¯¯¯ξnj+1, (3.13)

with defined by (3.2) at the troubled cell boundary.

3. Update the solution at the first RK stage as follows,

 ¯¯¯q(1)j=¯¯¯qnj−ΔtΔx(ˆfj+12−ˆfj−12),
 ¯¯¯ξ(1)j=¯¯¯ξnj−ΔtΔx(ˆHj+12−ˆHj−12).
Step 3.

Update and for a nonconvex troubled cell : we first identify nonconvex troubled cells by the following criterion.

Criterion II: A cell is called a good cell, if with , fall into the same linear, convex or concave region of the flux function . Otherwise, it is defined to be a nonconvex troubled cell.

Trouble cells.

At a nonconvex troubled cell , we apply a first order scheme on a refined mesh by evolving a time step . Specifically, we evolve equation (2.7) with the initial condition

 ¯¯¯qnl+snlσl, for x∈Il, l=j−1,j,j+1 (3.14)

where and