Alireza Mazaheri
email: alireza.mazaheri@nasa.gov

Chi-Wang Shu
email: shu@dam.brown.edu

Vincent Perrier
email: vincent.perrier@inria.fr

Bounded and compact weighted essentially nonoscillatory limiters for discontinuous Galerkin schemes: Triangular elements

Keywords:

 to high-order unstructured DG schemes. A compact positivity preserving limiter is applied to the solutions to ensure pressure and density remain within physical ranges at all time. It is then verified that the bounded WENO-DG maintains the formal order of accuracy of the underlying DG schemes in the smooth regions. The performance of the proposed WENO-DG is also demonstrated with inviscid test cases including the classical Riemann problems, shock-turbulence interaction, scramjet, blunt body flows, and the double Mach Reflection problems.

Introduction

Consider the unsteady compressible Euler equations for a perfect gas flow w t + ∇ • f(w) = 0, w(x, t = 0) = w 0 (x), [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF] where

w =   ρ ρu ρ E   , f(w) =   ρu ρuu T + pI ρuH   , H = E + p ρ , (2)
compact WENO-DG limiters presented here, therefore, require a stencil only as large as (nf aces) (0) + 1 number of elements, where (nf aces) (0) corresponds to the number of faces of the elements in which the limiter is being applied. For example, for triangles, the compact WENO stencil has only up to 4 elements (see also Fig. 1). In the first proposed approach, the compact stencil is considered, and a series of linear polynomials is constructed by solving a (dof × dof ) linear system, where dof denotes degrees of freedom. In the second approach, a linear polynomial is constructed by solving a series of constrained least-squares (LSQ) minimization problems within the compact stencil. In both approaches, a WENO polynomial limiter is obtained with a convex combination of the original polynomial and the constructed linear polynomials that are added using nonlinear weights. The numerical fluxes are evaluated with the Local Lax-Friedrichs (LLxF) flux, and a compact positivity preserving limiter [START_REF] Wang | Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations[END_REF] is also applied to the solutions to ensure pressure and density remain bounded and physical at all time.

The paper is organized as follows. Section 2 describes in detail the proposed compact WENO polynomial limiters for unstructured DG schemes with simplex elements. The positivity preserving limiter for the Euler equations is briefly discussed in Section 3. A few sample problems, ranging from the classical Sod and Lax problems to blunt body flows and double Mach reflection test case are presented in Section 4. The concluding remarks are given in Section 5.

Proposed WENO-DG polynomial limiters

Two new classes of compact WENO polynomial limiters are proposed here for unstructured DG schemes. These schemes are in part inspired by the recent successful WENO finite difference and finite volume procedures of Zhu and Shu [START_REF] Zhu | On the convergence to steady state solution of a new class of high order WENO schemes[END_REF][START_REF] Zhu | Numerical study on the convergence to steady state solutions of a new class of finite volume WENO schemes: triangular meshes[END_REF], which are extensions and steady state computations of the earlier works reported by Zhu and Qiu [START_REF] Zhu | A new fifth order finite difference WENO scheme for solving hyperbolic conservation laws[END_REF][START_REF] Zhu | A new type of finite volume WENO schemes for hyperbolic conservation laws[END_REF][START_REF] Zhu | New finite volume weighted essentially non-oscillatory schemes on triangular meshes[END_REF][START_REF] Zhu | On the convergence to steady state solution of a new class of high order WENO schemes[END_REF] and Levy et al. [START_REF] Levy | Central WENO schemes for hyperbolic systems of conservation laws[END_REF]. The proposed WENO procedures are designed for simplex elements, and are compact, requiring information only from the immediate neighbors. For each presented approach, a detailed step-by-step procedure is outlined that may be used for implementation of the WENO as a post-processing step to the DG solution after each time iteration.

Consider a compact triangular stencil {0, 1, 2, 3} as shown in Fig. 1, where the element 0 is the target element for which the limiter is applied. The goal is to reconstruct a new candidate polynomial for the target element, p (0)new , that is nonoscillatory in the presence of discontinuities, and dependent only on the available information within the compact stencil. We remark that the polynomials constructed throughout this document are defined with a set of basis functions expressed in a reference element.

The original polynomial on the target element, p (0) , may be expressed in the following identity relation as

p(0) := p (0) = γ 0 1 γ 0 p (0) - i γ i γ 0 p(i) + i γ i p(i) , (4)
which holds for arbitrary positive coefficients γ j satisfying 3 j=0 γ j = 1, where, γ j , are the linear weights defined as γ 0 = 0.997, γ j = 0.001 ( j = 0), and the p(i) are the to-be-reconstructed polynomials for the elements i = 1, 2, 3 (i.e., the immediate neighbors of the target element). The p(i) reconstruction procedures are discussed in detail later. The linear weights are then replaced with the nonlinear weights, ω, and the following reconstructed candidate polynomial is sought for the target element4

p (0)new = ω 0 1 γ 0 p (0) - i γ i γ 0 p(i) + i ω i p(i) , (5)
where the nonlinear weights are computed by adopting the high-order WENO-FD [START_REF] Zhu | On the convergence to steady state solution of a new class of high order WENO schemes[END_REF] and WENO-FV [START_REF] Zhu | Numerical study on the convergence to steady state solutions of a new class of finite volume WENO schemes: triangular meshes[END_REF] techniques of Zhu and Shu for the present WENO-DG schemes. These techniques are based on the WENO-Z strategy given in Refs. [START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF][START_REF] Castro | High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws[END_REF][START_REF] Don | Accuracy of the weighted essentially non-oscillatory conservative finite difference schemes[END_REF].

The nonlinear weights are then defined as

ω i = ω i j ω j , ω j = γ j 1 + τ (ǫ + β j ) 2 , τ = i |β 0 -β i | n 2 , ( 6
)
where n is the total number of reconstructed polynomials used in constructing the WENO polynomial p (0)new as given in Eq. (5) (i.e., maximum integer value of i) plus the polynomial of the target element 0 . Here ǫ = 10 -12 is a small number to avoid zero denominator, and β is the classical smoothness indicator defined as (see e.g., Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] and Kolb [START_REF] Kolb | On the full and global accuracy of a compact third order weno scheme[END_REF])

β i = k |l|=1 | 0 | (|l|-1) 0 ∂ |l| ∂x l 1 ∂ y l 2 p(i) 2 , l = (l 1 , l 2 ), (7)
where k denotes the polynomial order, and | 0 | is the volume of the target element. The above procedure is complete with reconstruction of the polynomials p(i) . This is discussed next where two compact reconstruction approaches are presented in detail.

Polynomial reconstruction approach 1

In this section, a procedure to reconstruct a set of polynomials, p(i) , is described in detail. The reconstructed polynomials are used for the WENO polynomial limiter given by Eq. ( 5). First, a procedure for a scalar case is given followed by a system case.

Scalar case

1. Compute the means of the polynomials on the compact WENO stencil {0, 1, 2, 3}:

u (i) = 1 | i | i p (i) , i = 0, 1, 2, 3, (8)
where i denotes neighboring elements with corresponding volume of | i |.

| i | i p( j) = u (i) , (9)
where i corresponds to the element number in the compact WENO stencil (see Fig. 1), and j refers to the candidate linear polynomial number. For a triangle, we have three such candidates, and for each candidate, we use the elements in the ith set to construct a linear polynomial p( j) :

( j, i) := (1, {0, 1, 2}); (2, {0, 1, 3}); (3, {0, 2, 3}). ( 10
)
Remark. All the reconstructed linear polynomials are defined within the target element, with the basis functions, φ, also defined on the target element. Such basis functions are denoted as φ(0) i . For example, consider the j = 1 polynomial, p(1) . This is defined as

p(1) = 2 i=0 ũi φ(0) i . ( 11
)
This polynomial is then used for evaluation of the integral given in Eq. ( 9) with appropriate quadrature rules. The elements in the ith set {0, 1, 2} are needed in the construction of the polynomial p(1) . Therefore, a 3 × 3 linear system is solved to obtain the three coefficients (degrees of freedom) of the reconstructed linear polynomial. Thus, for the j and ith set used in this example, the following linear system is solved:

     1 | 0 | q w q ϕ (0) 0 | q 1 | 0 | q w q ϕ (0) 1 | q 1 | 0 | q w q ϕ (0) 2 | q 1 | 1 | q w q ϕ (0→1) 0 | q 1 | 1 | q w q ϕ (0→1) 1 | q 1 | 1 | q w q ϕ (0→1) 2 | q 1 | 2 | q w q ϕ (0→2) 0 | q 1 | 2 | q w q ϕ (0→2) 1 | q 1 | 2 | q w q ϕ (0→2) 2 | q        ũ0 ũ1 ũ2   =   u (0) u (1) u (2)   , (12)
where subscript q denotes summation over the element quadrature points, w q corresponds to the quadrature weights, and the arrow indicates the extrapolation of the basis onto the neighbor elements 1 and 2. The above system is valid for interior elements as shown in Fig. 1a. For elements that have faces on physical boundaries, the above system (12) takes a different form. For instance, consider ith set {1, ∂ 1 , ∂ 2 } where ∂ denotes the corresponding boundary face of the target element that is on a boundary. In this case, the system (12) reduces to

     1 | 0 | q w q ϕ (0) 0 | q 1 | 0 | q w q ϕ (0) 1 | q 1 | 0 | q w q ϕ (0) 2 | q 1 |∂ 1 | ∂q w ∂q ϕ (0) 0 | ∂q 1 |∂ 1 | ∂q w ∂q ϕ (0) 1 | ∂q 1 |∂ 1 | ∂q w ∂q ϕ (0) 2 | ∂q 1 |∂ 2 | ∂q w ∂q ϕ (0) 0 | ∂q 1 |∂ 2 | ∂q w ∂q ϕ (0) 1 | ∂q 1 |∂ 2 | ∂q w ∂q ϕ (0) 2 | ∂q        ũ0 ũ1 ũ2   =   u (0) u ∂ 1 u ∂ 2   , (13)
where ∂ denotes the length of the corresponding boundary face, ∂q indicates summation over the corresponding boundary face quadratures, and u ∂ l is the mean value evaluated on the corresponding boundary face of the target element after the boundary condition is applied.

The same process is repeated to construct the other two linear polynomials p(2) and p(3) .

3. Project the three basis functions of the constructed linear polynomials onto the basis space of the original polynomial of the target element, p (0) ; i.e., the (k

+ d)!/(k d!) basis functions,
where k is the polynomial order, and d is the element dimension.

4. Reconstruct the new WENO limiter as given in Eq. ( 5) with i = 1, 2, 3.

System case

1. Follow steps 1 and 2 given for the scalar case.

2. Project p( j) , j = 0, 1, 2, 3 onto the characteristic fields

p( j) i = L i • p( j) , (14)
where L i are the left eigenvectors based on the mean values of the target cell and the normal direction i. For a triangle, there are three normal directions corresponding to each face of the element. Note: p(0) := p (0) .

U N C O R R E C T E D P R O O F
•• (••••) •••-•••
••• (••••) •••-•••

Use

p j i , and follow steps 3 and 4 given for the scalar case, and evaluate p(0) ne w i for each direction i; i.e., p(0

) ne w i = ω 0 γ 0 p(0) i + j ω j -ω 0 γ j γ 0 p( j) i , j = 1, 2, 3, i = 1, 2, 3 . (15)
p (0) ne w i = R i • p(0) ne w i , (16)
where R i are the right eigenvectors based on the mean values of the target cell and the normal direction i. 5. Obtain the final reconstructed WENO polynomial for the target cell with weighted averaging:

p (0) ne w = i p (0) ne w i | i | i | i | , ( 17
)
where | i | corresponds to the volume of the neighboring element for face i of the target element. If the face is a boundary face, use the volume of the target element, | 0 |.

U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••

Polynomial reconstruction approach 2

In this section, a second approach in constructing the polynomials p(i) is presented. In this approach, a constraint minimization problem is solved using the least squares (lsq) approach. The proposed lsq linear polynomial reconstruction procedure for unstructured DG schemes is also compact, and its implementation is nonintrusive similar to the first presented approach.

The aim here is to construct a set of linear polynomials, plsq , on the target cell such that arg min

ũ l∈S    1 | l | l plsq ( ũ, x) dx -u (l)    2 , S = {1, 2, 3} , (18)
where l corresponds to the immediate neighbors of the target element, is satisfied subject to a constraint that the means of the constructed polynomials remain the same as the means of the original polynomials on the target element; i.e., 0

plsq dx = u (0) . (19)
Thus, the following minimization problem is sought by employing a penalization parameter (called Lagrange multiplier) λ: Test case where ũi are the coefficients of the to-be-constructed linear polynomial for ith degrees of freedom, dof , and | l | is the volume of the neighbor l. The constraint minimization problem may be expanded and expressed with employing proper quadrature rules as

E( ũi , λ) = arg min ũi    l∈S    1 | l | l plsq ( ũi , x) dx -u (l)    2 + λ    1 | 0 | 0 plsq ( ũi , x) dx -u (0)       , (20)
U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
ρ ρu ρ v ρ E ρ ρu ρ v ρ E 20×20 2.66E-1 3.679E-4 2.752E-4 2.795E-4 1.153E-3 - - - - 40×40 1.40E-1 4.803E-5 3.425E-5 3.455E-
ρ ρu ρ v ρ E ρ ρu ρ v ρ E 20×20 2.66E-1 3.679E-4 2.752E-4 2.795E-4 1.153E-3 - - - - 40×40 1.40E-1 4.803E-5 3.425E-5 3.455E-
ρ ρu ρ v ρ E ρ ρu ρ v ρ E 20×20 2.66E-1 7.134E-6 7.199E-6 7.190E-6 2.358E-5 - - - - 40×40 1.40E-
ρ ρu ρ v ρ E ρ ρu ρ v ρ E 20×20 2.66E-1 2.285E-7 1.849E-7 1.879E-7 7.260E-7 - - - - 40×40 1.40E-
ρ ρu ρ v ρ E ρ ρu ρ v ρ E 20×20 2.66E-1 2.285E-7 1.849E-7 1.879E-7 7.260E-7 - - - - 40×40 1.40E-
••• (••••) •••-•••
ρ L u L p L ρ R u R p R x d t f inal
E( ũi , λ) = arg min ũi   l∈S 1 | l | i q ũi φ(0→l) i,q ω (l) q | det J | (l) -u (l) 2 + λ 1 | 0 | i q ũi φ(0) i,q ω (0) q | det J | (0) -u (0)   , (21)
where ω (l) q is the weight of the quadrature for the qth quadrature point of the lth neighbor,

| det J | (l) = | 0 | (l)
is the determinant of the Jacobian of the transformation from the reference to the physical element of neighbor l, and ϕ (0→l) i,q is the ith basis function in the reference element of the target element extrapolated to the lth neighbor, and evaluated at the corresponding qth quadrature point.

The minimization problem is solved by setting the derivatives of E w.r.t. the jth dof coefficient ũ j and λ to zero; i.e., + λ

∂ E ∂ ũ j = l 1 | l | q φ(0→l) j,q ω (l) q | det J | (l) 1 | l | i q ũi φ(0→l) i,q ω (l) q | det J | (l) -u (l) U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
1 | 0 | q φ(0) j,q ω (0) q | det J | (0) = 0, (22)
∂ E ∂λ = 1 | 0 | i q ũi φ(0) i,q ω (0) q | det J | (0) -u (0) = 0. ( 23
)
This may be expressed in a matrix form as

l 1 | l | 2 l φ(0→l) j l φ(0→l) i 1 | 0 | 0 φ(0) j sym. 0 A ji ũi λ U i = l u (l) | l | | l | φ(0→l) j u (0) B j , ( 24
)
where A ji is a symmetric matrix, and λ = λ/2.

It is important to note that whenever the integration is performed on the neighbor element l, the linear polynomial basis defined on the target element, φ(0) , must be extrapolated to the neighboring element. This is emphasized in Eq. ( 24) by employing the notation ϕ (0→l) j . For convenience and simplicity in the discussion, however, this notation is removed in the subsequent text, and such extrapolation is therefore assumed to be understood whenever the integration is performed on the neighboring elements.

For interior triangular elements similar to one shown in Fig. 1a, the matrix A ji is a symmetric 4×4 matrix, and U i and B j are vectors of length four (three for dof + one for the constraint). They take the following forms

U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
A ji =         l 1 | l | 2 l φ(0) 0 2 l 1 | l | 2 l φ(0) 0 | l | φ(0) 1 l 1 | l | 2 l φ(0) 0 | l | φ(0) 2 1 | 0 | 0 φ(0) 0 S ym. l 1 | l | 2 l φ(0) 1 2 l 1 | l | 2 l φ(0) 1 l φ(0) 2 1 | 0 | 0 φ(0) 1 S ym. S ym. l 1 | l | 2 l φ(0) 2 2 1 | 0 | 0 φ(0) 2 S ym. S ym. S ym. 0         , ( 25
)
U i =     ũ0 ũ1 ũ2 λ     , B j =         l u (l) | l | l φ(0) 0 l u (l) | l | l φ(0) 1 l u (l) | l | l φ(0) 2 u (0)         . ( 26
)
For triangular elements whose faces may be on boundaries, such as the one depicted in Fig. 1b, where two of the target element faces are on the boundary, the set of neighbors in the linear reconstruction is therefore l = {1, ∂ 1 , ∂ 2 }, where ∂ denotes the corresponding boundary face of the target element that is on a boundary. For instance, consider the triangle shown in Fig. 1b, for which the following expressions for the A 11 and its corresponding right-hand-side B 1 are obtained

A 11 = 1 | 1 | 2    1 φ(0) 0    2 + 1 |∂ 1 | 2    ∂ 1 φ(0) 0    2 + 1 |∂ 2 | 2    ∂ 2 φ(0) 0    2 , (27)
U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
| 1 | | 1 | φ(0) 0 + u ∂ 1 |∂ 1 | ∂ 1 φ(0) 0 + u ∂ 2 |∂ 2 | ∂ 2 φ(0) 0 , (28)
where u ∂ l is the mean value evaluated on the corresponding boundary face of the target element after the boundary condition is applied. Note that ∂ l denotes that the integration is performed on the corresponding boundary face of the target element and thus, no extrapolation is needed for the last two terms of the expression [START_REF] Levy | Central WENO schemes for hyperbolic systems of conservation laws[END_REF]; extrapolation must be applied to the first term of the expression. The complete step-by-step procedure with the second presented approach is outlined below.

Scalar case

1. Compute the means of the polynomials on the compact WENO stencil {0, 1, 2, 3}. This is the same as the step 1 of the approach 1.

2. Construct the plsq polynomial as outlined above.

3. Use plsq , and follow steps 3 and 4 of approach 1, and evaluate p(0) ne w ; i.e., 2. Project p(0) and plsq onto the characteristic field by multiplying them with the left eigenvectors

p (0) ne w = ω 0 γ 0 p (0) + ω lsq -ω 0 γ lsq γ 0 plsq . ( 29
••• (••••) •••-•••
p(0) i = L i • p(0) , (30)
plsq i = L i • plsq , (31)
where L i are the left eigenvectors based on the mean values of the target cell and the normal direction i. For triangle, there are three normal directions corresponding to each face of the element. Note: p(0) := p (0) .

Evaluate

p(0) ne w i for each direction i using p(0) i and plsq i ; i.e., p(0

) ne w i = ω 0 γ 0 p(0) i + ω lsq -ω 0 γ lsq γ 0 plsq i . (32)
p (0) ne w i = R i • p(0) ne w i , (33)
where R i are the right eigenvectors based on the mean values of the target cell and the normal direction i. 5. Obtain the final reconstructed WENO polynomial for the target cell with weighted averaging:

p (0) ne w = i p (0) ne w i | i | i | i | , (34)

Positivity preserving limiter (PPL) for Euler

The proposed WENO-DG does not necessary bound the density and pressure within their physical values and thus, a proper bound preserving strategy must be adopted to avoid occurrences of unphysical quantities in time and/or space. The bound preserving limiter must also preserve the formal order of accuracy of the underlying DG scheme. Here, a positivity preserving algorithm of Wang et al. [START_REF] Wang | Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations[END_REF] is adopted and applied to the polynomials in conjunction with the WENO polynomial limiters. The presented PPL limiter for Euler equations, Eq. ( 1), is completely local requiring only information within each element, and therefore, enabling an extremely efficient parallel implementation. Similar to the presented WENO-DG, the given PPL procedure is also nonintrusive and may be applied as a postprocessing step to the polynomials after each time iteration. The step-by-step PPL procedure for DG methods is:

1. Compute the means of the polynomials w = (ρ, ρu, ρ E) T on each element

w i = 1 w i , (35)
where, i = 0 . . . (d + 1) denotes the indices of the vector of conservative variables w, and d refers to the dimension.

Construct a new polynomial for density as follows

ρ = θ (1) (ρ -ρ) + ρ, θ (1) = min x∈q 1, ρ -ǫ ρ -ρ q (x) , ( 36
)
where ǫ is a small number, which is taken as 10 -13 in this study, and subscript q denotes values defined on Gauss-Lobatto quadrature points. This step enforces positivity for density.

U N C O R R E C T E D P R O O F

••• (••••) •••-•••
4. Compute the following scaling factor for every quadrature point, q

θ q (x) = 1 : p( ŵ) ≥ 0 p(w)/ p(w) -p( ŵq (x)) : otherwise . ( 38
)
5. Limit the ŵ polynomials to obtain w w = θ (2) ŵw + w, θ (2) = min x∈q θ q (x) .

(39)

6. Replace w with the limited polynomials w after each time iteration. For explicit time schemes such as SSP RK, this procedure is applied after each RK stage.

Numerical results

The limiters' abilities to achieve the desired order of accuracy is verified first using an inhomogeneous Euler system. A few sample examples are then presented to assess the performance of the proposed limiters. These include the classical Riemann problems (Sod and Lax), the Shu-Osher shock-density interactions, the Mach 3 wind tunnel case of Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF], a Mach 3 two-strut scramjet, Mach 3 and Mach 10 blunt body problems, and the Mach 10 double Mach reflection problem. In all the test cases presented here, the corresponding domains are discretized with irregular triangular elements and simulations are carried with a CFL condition that is formally proven to be stable [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] under the following condition

U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
||λ|| dt h ≤ 1 2(2k + 1) , ( 40
)
where ||λ|| is the magnitude of the largest characteristic quantity of the hyperbolic system, dt is the time step, h is the element size, and k is the polynomial order. For each test case, the proposed WENO-DG schemes are applied and the predicted results are shown for second-, third-, fourth-, and fifth-order solutions. The computed nonlinear weights profiles are also presented.

Verification -inhomogeneous Euler

Consider the two-dimensional Euler equations [START_REF] Eberle | Generalized fluxvectors for hypersonic shock capturing[END_REF] which is imposed on the domain boundaries. A series of randomly generated irregular triangular grids is considered (see Fig. 2), and steady state solutions of the DG ( P k), k = 1, 2, 3, 4, scheme with and without the presented WENO and PPL limiters are obtained. The convergence history plots for third-, fourth-, and fifth-order cases shown in Figs. 3 and4 indicate that the residuals have settled down to 10 -12 -10 -15 values (machine zero) for both presented WENO approaches. Fig. 5 shows the difference in residuals history for both explicit SSP (3,3) RK and implicit Euler backward time discretization schemes.

w t + ∇ • f(w) = S, w(x, t = 0) = w 0 (x),
ρ = 1 + 0.2 sin(x + y), u = 1, v = 1, p = 1 + 0.2 sin(x + y),
The predicted high-order DG ( P k) solutions are also compared with the exact solution on the given grids by computing the L 2 error. The resulting L 2 errors are tabulated in Tables 1-4 and shown in Fig. 6. The L 2 errors for both WENO limiters are identical to the decimal point shown and therefore only one set of values is given. The WENO and positivity preserving limiters are both applied to all the elements within the domain. The (k + 1)th order of accuracy for both DG and the WENO-DG schemes are verified, confirming that the proposed WENO maintains the expected order of accuracy of the underlying DG ( P k) scheme. It is also interesting to note that the presented error magnitudes are order of magnitude smaller than the corresponding results reported with the third-, fourth-, and fifth-order WENO-FV schemes in Ref. [START_REF] Zhu | A new third order finite volume weighted essentially non-oscillatory scheme on tetrahedral meshes[END_REF].

Riemann problems

Here, a two-dimensional irregular triangular domain is considered to test the WENO scheme against the classical one-dimensional Sod and Lax problems. The computational domain (x, y) ∈ [-0.5, -0.05] × [0.5, 0.05] is discretized with triangular elements with a characteristic mesh size of h/100; see Fig. 7. The two states of the gas (γ = 1.4) are separated 5. The proposed WENO-DG is then applied to these problems, and the simulations are continued until the t f inal is reached.

U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
Figs. 8-13 show the predicted solutions, (ρ, u, p), against the exact Sod and Lax values. The predicted fifth-order WENO-DG results are comparable with the fifth-order ADER-CWENO scheme of Dumbser et al. [START_REF] Dumbser | Central WENO schemes for hyperbolic conservation laws on fixed and moving unstructured meshes[END_REF].

Shu-Osher problem

The one-dimensional idealization of the shock-turbulence interaction suggested by Shu-Osher [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes II[END_REF] is considered in a two-dimensional framework. The goal of this test is to assess the proposed limiters in capturing the shock wave and its interactions with the unsteady low frequency density fluctuations and its waves propagations downstream of the shock. The

This corresponds to a Mach 3 shock (γ = 1.4) interacting with the sine waves density field. The results at t = 1.8 are presented in Figs. [START_REF] Qiu | Hermite WENO schemes and their application as limiters for Runger-Kutta discontinuous Galerkin method: one dimensional case[END_REF][START_REF] Qiu | Runge-Kutta discontinuous Galerkin method using WENO limiters[END_REF][START_REF] Qiu | Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method II: two dimensional case[END_REF]. The fine structured of the shock-density wave interactions are clearly captured by the proposed limiter. to remove the known artifacts that is present in the second-order result. These artifacts are largely removed in the higher-order results.

Mach 3 forward facing step

Consider the forward facing step problem that was originally proposed by Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. This test is often referred to as the Mach 3 wind tunnel test. The computational domain is shown in Fig. 17. The initial conditions correspond to a uniform flow moving to the right with Mach 3, (ρ, u, v, p) = (1.4, 3.0, 0.0, 1.0). The inflow condition is imposed to the left boundary while a do-nothing boundary condition is set to the right boundary. A reflecting boundary condition is applied to other surfaces, and the flow is initialized with the left boundary values. The domain (see Fig. 17) is discretized using irregular triangular elements with characteristic element size of h/160. The high-order solutions at t = 4.0 are presented in Fig. 18. Note that no modification, neither to the scheme nor to the grid resolution, is applied to the corner singularity. The corner step singularity is known to produce an erroneous entropy layer in addition to spurious Mach stem at the bottom wall. These artifacts are clearly present in the presented second-order result; the predicted higher order solutions are almost free from these artifacts.

Scramjet

Consider a two-strut scramjet (see Fig. 19 and Table 6) with a Mach 3 inflow imposed on the left surface boundary. For this test, the second proposed approach produced nearly machine zero residuals as illustrated in Fig. 20. The profiles of the computed nonlinear weights along the centerline of the two-strut scramjet are shown in Fig. 21. The corresponding density contours in the range of 1.5 and 8.0 are shown in Fig. 22. These results are in agreement with the previously published results [START_REF] Kumar | Numerical Analysis of Scramjet-Inlet Flow Field by Using Two-Dimensional Navier-Stokes Equations[END_REF][START_REF] Eberle | Generalized fluxvectors for hypersonic shock capturing[END_REF][START_REF] Alauzet | Estimateur dérreur géométrique et métriques anisotropes pour l'adaptation de maillage: Partie II: exemples d'applications[END_REF] on a similar configuration; the exact geometrical information used in generating the previously published articles was not known to the authors. 6.

Table 6

The coordinates of the lower half portion of the two-strut scramjet geometry shown in Fig. 19. The strong shocks are correctly captured but there exists some waviness in the density contour lines, particularly when the polynomial order and the intensity of the shock are increased. However, these striation-like patterns, which are due to irregularity of the grid elements in the vicinity of the shock, appear to have minimal effects on the surface quantities as illustrated in Fig. 26. Further improvement to the proposed WENO is needed to reduce these striation-like patterns in the post shock regions. This will be reported in future studies.

Double Mach reflection test case

The double Mach reflection problem is originally proposed by Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] as a benchmark test for Euler solvers. The problem consists of a Mach 10 shock front that meets a 30-deg inclined ramp. The shock front initially sits at x 0 = 1/6, and makes a 60-deg angle with the reflecting wall (x-axis). The properties of the undisturbed air (γ = 1.4) ahead of the shock are ρ = 1.4 and p = 1. Hence, the shock speed is |u s | = 10. As it is a customary, this problem is solved in a computational domain (x, y) ∈ [0, 4] × [0, 1] with a coordinate system that is aligned with the ramp; see Fig. 27. The domain is discretized with irregular triangular elements and a mesh size of h = 1/200.

U N C O R R E C T E D P R O O F Please cite
••• (••••) •••-•••
Considering the ramp-aligned coordinate system, the problem is therefore setup with the following initial condition:

(ρ, u, v, p)| t=0 = (8.0, +8.25 cos(π /6), -8.25 sin(π /6), 116.5), x < x s , (1.4, 0.0, 0.0, 1.0), x ≥ x s ,

where the shock position is computed as,

x s (y) = x 0 + y tan(π /6).

For the boundary conditions, the post-shock values are imposed for the short region from x = 0 to x s along the lower boundary at y = 0, while for the rest of the lower boundary, a reflecting boundary condition is imposed. This ensures that the reflecting shock is attached to the wall. The initial post-shock condition is also assigned at the left boundary along the x = 0 axis. The boundary condition on the upper boundary along the y = |u s |/ cos(π /6) t. The readers are referred to Ref. [START_REF] Kemm | On the proper setup of the double Mach reflection as a test case for the resolution of gas dynamics codes[END_REF] for discussions on the effects of the upper boundary condition and the computational domain size on numerical artifacts in the predicted solutions; such effects are not investigated in this study.

••• (••••) •••-••• 25
The computed nonlinear weights along the y = 0.3 is plotted in Fig. 28. The density contours for a truncated domain (x, y) = [0, 3.2] × [0, 1] (the computational domain extends to x = 4) with 41 equally spaced isolines from 1.5 to 22.5 at t = 0.2 are shown in Fig. 29 for the second-, third-, fourth-, and fifth-order WENO-DG. The closeup views of the triple Mach points region, (x, y) ∈ [2, 2.9] × [0, 0.6], are shown in Fig. 30. Clearly, the complexity of the triple Mach points is captured by the proposed limiters on the irregular triangular mesh. The quality of the solution improves with increasing the polynomial order on the same identical irregular triangular mesh. The resolution of the curled flow structures along the primary slip line, which is caused by its interactions with the secondary reflected shock emanating from the secondary triple point, is often used for judging the quality of the numerical scheme (see Ref. [START_REF] Kemm | On the proper setup of the double Mach reflection as a test case for the resolution of gas dynamics codes[END_REF] for definitions of the primary and the secondary triple points and slip lines). The presented results provide a remarkable curled flow structures compared to the solutions reported by Hu and Shu [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF], Zhu et al. [START_REF] Zhu | Runge-Kutta discontinuous Galerkin method with a simple and compact Hermite WENO limiter[END_REF], and Dumbser et al. [START_REF] Dumbser | Central WENO schemes for hyperbolic conservation laws on fixed and moving unstructured meshes[END_REF]. Further improvement in the WENO could reduce the noise in the predicted contour lines. Grid adaptation could also enhance the results further.

Concluding remarks

Two compact WENO limiters were proposed for DG schemes for irregular triangular elements, and detailed step-by-step construction procedures were outlined. The WENO-DG limiters were bounded with a compact bound-preserving limiter for the Euler equations. It was verified that the proposed bounded WENO-DG polynomial limiters preserve the expected order of accuracy of the underlying DG schemes when they are applied to smooth regions. A few test cases, including the classical Riemann problems, Shu-Osher shock-turbulence interaction, scramjet, blunt body, and double Mach reflection problems, involving strong shocks were performed, and the results for second-, third-, fourth-, and fifth-order compact and bounded WENO-DG are presented. It was shown that the second proposed WENO-DG is more effective for some problems in

Highlights

• Compact high-order WENO for DG methods on triangular elements.

• Bounded polynomial limiter with a positivity preserving.

• (k + 1)th order of accuracy on irregular triangular elements.

• Nonintrusive and simple implementation procedures.

View publication stats View publication stats

Fig. 1 .

 1 Fig. 1. The compact WENO-DG stencil for triangular elements.

  For boundary elements, first obtain the right values u R by applying the boundary condition along the boundary face, and then integrate u R to get the mean values. If1

Fig. 2 .

 2 Fig. 2. A sample of randomly generated irregular triangular elements.

Fig. 4 .

 4 Fig. 4. Verification -Residual history for the finest mesh.

Fig. 5 .

 5 Fig. 5. Example 1 -Residual history for a) third-order strong stability preserving (SSP) Runge-Kutta time scheme, and b) Backward Euler time scheme.

Fig. 6 .

 6 Fig. 6. Verification -Order of accuracy plots for the conserved density, x and y momentums, and energy for DG and the proposed WENO-DG on irregular triangular elements.

Fig. 7 . 5

 75 Fig. 7. Schematic of the Riemann problem, and the irregular triangular grid with h = 1/100. The two states of the gas are separated at x d = 0.

Fig. 8 .

 8 Fig. 8. Predicted density profile for the Sod problem with h/100 irregular triangular elements at t = 0.2, and comparison against the exact solution.

Fig. 9 .

 9 Fig. 9. Predicted velocity profile for the Sod problem with h/100 irregular triangular elements at t = 0.2, and comparison against the exact solution.

Fig. 10 .

 10 Fig. 10. Predicted pressure profile for the Sod problem with h/100 irregular triangular elements at t = 0.2, and comparison against the exact solution.

Fig. 11 .

 11 Fig. 11. Predicted density profile for the Lax problem with h/100 irregular triangular elements at t = 0.14, and comparison against the exact solution.

Fig. 12 .

 12 Fig. 12. Predicted velocity profile for the Lax problem with h/100 irregular triangular elements at t = 0.14, and comparison against the exact solution.

Fig. 13 .

 13 Fig. 13. Predicted pressure profile for the Lax problem with h/100 irregular triangular elements at t = 0.14, and comparison against the exact solution.3. Create a new vector of conservative variables ŵ with the new density polynomial ρ obtained in the previous step

Fig. 14 .

 14 Fig. 14. Predicted density profile for the Shu-Osher problem with h/40 irregular triangular elements at t = 1.8.

( 41 )

 41 with the vector of conservative variables w and the source S, cos(x + y) 0.6 cos(x + y) 0.6 cos(x + y) 1.8 cos(x + y)

Fig. 15 .

 15 Fig. 15. Predicted velocity profile for the Shu-Osher problem with h/40 irregular triangular elements at t = 1.8.

Fig. 16 .

 16 Fig. 16. Predicted entropy profile, ln(p/ρ γ ), for the Shu-Osher problem with h/40 irregular triangular elements at t = 1.8.

  computational domain (x, y) ∈ [-5.0, -0.1] × [5.0, 0.1] is discretized with triangular elements and a characteristic mesh size of h/40. The domain is initialized as (ρ, u, v, p)| t=0 = (3.857143, 2.629369, 0., 10.33333), x < 4.0, (1. + 0.2 sin(5x), 0., 0., 1.), x ≥ 4.0.

Fig. 17 .

 17 Fig. 17. The geometry of the Mach 3 wind tunnel forward facing step with a sample irregular grid.

Fig. 18 .

 18 Fig. 18. High-order DG with proposed WENO and positivity preserving with irregular triangular elements (h = 1/160) for the forward facing step problem. Shown are 30 equidistance density contour lines at t = 4.0 from 0.32 to 6.15. The corner step singularity neither in the scheme nor by the grid is treated

UFig. 19 .

 19 Fig. 19. The geometry of the two-strut scramjet test case with a sample irregular grid with 10,000 vertices. Coordinates of the lower half portion of the geometry is given in Table6.

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Absolute value of the elemental residuals obtained with the second proposed WENO approach for the Mach 3 two-strut scramjet test case.

Fig. 22 .Fig. 23 .

 2223 Fig.22. Hundred equally spaced density contours in the range of 1.5 and 8.0 using the second presented WENO approach for the Mach 3 two-strut scramjet test case with 10,000 unadapted irregular triangular elements.

Fig. 24 .

 24 Fig. 24. Nonlinear weights along the stagnation line of the Mach 10 blunt body problem.

4. 6 .Fig. 25 .

 625 Fig. 25. Thirty equally spaced density contours in the range of 1 and 8 for the Mach 3 (top row) and Mach 10 (bottom row) blunt body test cases.

Fig. 26 .

 26 Fig. 26. Top row: Mach 3. Bottom row: Mach 10; a,d) density profile along a line normal to the stagnation point, b,e) surface density, and c,f) surface temperature. Data are extracted from a low-order visualization.

Fig. 27 .

 27 Fig. 27. Schematic of the double Mach reflection, and the computational domain (x, y) ∈ [0, 4] × [0, 1] with the ramp-aligned coordinate system.

Fig. 28 .

 28 Fig. 28. Nonlinear WENO weights, ω j , along the y = 0.3 for the double Mach reflection case.

UFig. 29 .

 29 Fig. 29. High-order DG with proposed WENO and positivity preserving with irregular triangular elements (h = 1/200) in (x, y) ∈ [0, 4] × [0, 1]. Shown are 41 equidistance density contour lines at t = 0.2 from 1.5 to 22.5.

Fig. 30 . 28 A

 3028 Fig. 30. High-order DG with proposed WENO and positivity preserving with irregular triangular elements (h = 1/200) in (x, y) ∈ [0, 4] × [0, 1]. Shown are zoom-in views around the Mach stems at t = 0.2 with 41 equidistance density contour lines from 1.5 to 22.5.

  this article in press as: A. Mazaheri et al., Bounded and compact weighted essentially nonoscillatory limiters for discontinuous Galerkin schemes: Triangular elements, J. Comput. Phys. (2019), https://doi.org/10.1016/j.jcp.2019.06.023

	JID:YJCPH AID:8747 /FLA	[m3G; v1.260; Prn:17/06/2019; 14:18] P.7 (1-28)
	A. Mazaheri et al. / Journal of Computational Physics

Table 1

 1 Verification -Order of accuracy verification for the second-order DG with and without WENO on irregular triangular elements.

	Grid	h size	L 2 Error: DG (P 1)				Order: DG (P 1)
			ρ	ρu	ρ v	ρ E	ρ	ρu	ρ v	ρ E
	20×20	2.66E-1	6.015E-3	6.323E-3	6.328E-3	2.056E-2	-	-	-	-
	40×40	1.40E-1	1.452E-3	1.520E-3	1.524E-3	4.960E-3	2.21	2.22	2.21	2.22
	60×60	8.93E-2	5.950E-4	6.279E-4	6.270E-4	2.056E-3	1.99	1.97	1.98	1.96
	80×80	6.93E-2	3.422E-4	3.602E-4	3.595E-4	1.176E-3	2.18	2.19	2.19	2.20
	100×100	5.61E-2	2.170E-4	2.255E-4	2.260E-4	7.449E-4	2.17	2.23	2.21	2.17
	Grid	h size	L 2 Error: WENO-DG (P 1)			Order: WENO-DG (P 1)
			ρ	ρu	ρ v	ρ E	ρ	ρu	ρ v	ρ E
	20×20	2.66E-1	6.015E-3	6.323E-3	6.328E-3	2.056E-2	-	-	-	-
	40×40	1.40E-1	1.452E-3	1.520E-3	1.524E-3	4.960E-3	2.21	2.22	2.21	2.22
	60×60	8.93E-2	5.950E-4	6.279E-4	6.270E-4	2.056E-3	1.99	1.97	1.98	1.96
	80×80	6.93E-2	3.422E-4	3.602E-4	3.595E-4	1.176E-3	2.18	2.19	2.19	2.20
	100×100	5.61E-2	2.170E-4	2.255E-4	2.260E-4	7.449E-4	2.17	2.23	2.21	2.17

Table 2

 2 Verification -Order of accuracy verification for the third-order DG with and without WENO on irregular triangular elements.

	Grid	h size	L 2 Error: DG (P 2)	Order: DG (P 2)

Table 3

 3 Verification -Order of accuracy verification for the fourth-order DG with and without WENO on irregular triangular elements.

	Grid	h size	L

[START_REF] Gottlieb | High order strong stability preserving time discretizations[END_REF]

Error: DG (P 3) Order: DG (P 3)

Table 4

 4 Verification -Order of accuracy verification for the fifth-order DG with and without WENO on irregular triangular elements.

	Grid	h size	L 2 Error: DG (P 4)	Order: DG (P 4)

  this article in press as: A. Mazaheri et al., Bounded and compact weighted essentially nonoscillatory limiters for discontinuous Galerkin schemes: Triangular elements, J. Comput. Phys. (2019), https://doi.org/10.1016/j.jcp.2019.06.023 JID:YJCPH AID:8747 /FLA [m3G; v1.260; Prn:17/06/2019; 14:18] P.24 (1-28) 24 A. Mazaheri et al. / Journal of Computational Physics

Classical trouble cell indicators are not employed in the presented examples and therefore, the proposed WENO is applied to all the computational elements. Hence, all the elements in the computational domain are considered as target elements. This approach is beneficial in ensuring that the proposed WENO does not affect the desired order of accuracy of the DG scheme even if it is artificially applied to smooth region.

(V. Perrier 1 Research supported by the NASA Langley CIF/IRAD project, and the NASA STMD ESM/GCD program. 2 Research supported by ARO grant W911NF-15-1-0226 and NSF grant DMS-1719410. 3 Research supported by NASA Langley CIF/IRAD project while visiting NASA Langley from January -March 2017.

U N C O R R E C T E D P R O O F

Please cite this article in press as: A. Mazaheri