NASA/TM–2009–216166

Boundary Closures For Fourth-order

Energy Stable Weighted Essentially

Non-Oscillatory Finite Difference

Schemes

Travis C. Fisher

Purdue University, West Lafayette, Indiana

Mark H. Carpenter

Langley Research Center, Hampton, Virginia

Nail K. Yamaleev

North Carolina A&T State University, Greensboro, North Carolina

Steven H. Frankel

Purdue University, West Lafayette, Indiana

December 2009
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated

to the advancement of aeronautics and space

science. The NASA scientific and technical

information (STI) program plays a key part

in helping NASA maintain this important

role.

The NASA STI Program operates under the

auspices of the Agency Chief Information

Officer. It collects, organizes, provides for

archiving, and disseminates NASA’s STI.

The NASA STI Program provides access to

the NASA Aeronautics and Space Database

and its public interface, the NASA Technical

Report Server, thus providing one of the

largest collection of aeronautical and space

science STI in the world. Results are

published in both non-NASA channels and

by NASA in the NASA STI Report Series,

which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results

of NASA programs and include extensive

data or theoretical analysis. Includes

compilations of significant scientific and

technical data and information deemed to

be of continuing reference value. NASA

counterpart of peer-reviewed formal

professional papers, but having less

stringent limitations on manuscript length

and extent of graphic presentations.

• TECHNICAL MEMORANDUM.

Scientific and technical findings that are

preliminary or of specialized interest, e.g.,

quick release reports, working papers, and

bibliographies that contain minimal

annotation. Does not contain extensive

analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION.

Collected papers from scientific and

technical conferences, symposia, seminars,

or other meetings sponsored or

co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to

NASA’s mission.

Specialized services also include creating

custom thesauri, building customized

databases, and organizing and publishing

research results.

For more information about the NASA STI

Program, see the following:

• Access the NASA STI program home page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA STI Help

Desk at 443-757-5803

• Phone the NASA STI Help Desk at

443-757-5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076–1320
NASA/TM–2009–216166

Boundary Closures For Fourth-order

Energy Stable Weighted Essentially

Non-Oscillatory Finite Difference

Schemes

Travis C. Fisher

Purdue University, West Lafayette, Indiana

Mark H. Carpenter

Langley Research Center, Hampton, Virginia

Nail K. Yamaleev

North Carolina A&T State University, Greensboro, North Carolina

Steven H. Frankel

Purdue University, West Lafayette, Indiana

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23681-2199

December 2009
Acknowledgments

Work performed in partial fulfillment of LaRC LARSS summer fellowship program, Technical

monitor Robert Baurle. Special thanks are extended to the faculty and staff of Technical

University of Delft, Delft, The Netherlands, for many stimulating discussions on this and other

topics.

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not

constitute an offical endorsement, either expressed or implied, of such products or manufacturers by the

National Aeronautics and Space Administration.

Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076-1320

443-757-5802
Abstract

A general strategy exists for constructing Energy Stable Weighted Essentially Non–

Oscillatory (ESWENO) finite difference schemes up to eighth-order on periodic domains. These ESWENO schemes satisfy an energy norm stability proof for both

continuous and discontinuous solutions of systems of linear hyperbolic equations.

Herein, boundary closures are developed for the fourth-order ESWENO scheme that

maintain wherever possible the WENO stencil biasing properties, while satisfying

the summation-by-parts (SBP) operator convention, thereby ensuring stability in

an L2 norm. Second-order, and third-order boundary closures are developed that

achieve stability in diagonal and block norms, respectively. The global accuracy

for the second-order closures is three, and for the third-order closures is four. A

novel set of non-uniform flux interpolation points is necessary near the boundaries

to simultaneously achieve 1) accuracy, 2) the SBP convention, and 3) WENO stencil

biasing mechanics.

Contents

1 Introduction 2

2 Definitions 4

2.1 Complementary Grids . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.2 Summation-By-Parts Operators . . . . . . . . . . . . . . . . . . . . . 5

2.2.1 Central Fourth-order SBP Operators . . . . . . . . . . . . . . 6

3 Expressing Summation-By-Parts Operators in Flux Form 8

3.1 General Structure of ESWENO Schemes . . . . . . . . . . . . . . . . 8

3.2 Nonuniform Flux-Points . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.3 Interpolation Operators and Preferred Weights . . . . . . . . . . . . 11

4 Energy Stabilization Terms 11

5 Numerical Tests 12

5.1 Eigenvalues of the ESWENO operator . . . . . . . . . . . . . . . . . 13

5.2 Linear Wave Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 14

5.2.1 Sine Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

5.2.2 Square Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

5.3 Euler Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

5.3.1 Subsonic One-Dimensional Nozzle . . . . . . . . . . . . . . . 17

5.3.2 Transonic One-Dimensional Nozzle . . . . . . . . . . . . . . . 17

5.3.3 Two Dimensional Inviscid Vortex . . . . . . . . . . . . . . . . 18

5.4 Supersonic Shear Layers . . . . . . . . . . . . . . . . . . . . . . . . . 19

6 Conclusions 21

A Smoothness Indicators 23

1
B Recipe 25

C Fourth-order Diagonal-Norm Operator: D2-4-2 27

C.1 Differentiation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . 27

C.2 Interpolation Operators . . . . . . . . . . . . . . . . . . . . . . . . . 28

C.3 Target Weights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

C.4 Energy Stable Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

D Fourth-order Block-Norm Operator: D3-4-3 31

D.1 Differentiation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . 31

D.2 Interpolation Operators . . . . . . . . . . . . . . . . . . . . . . . . . 32

D.3 Target Weights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

D.4 Energy Stable Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

E References 34

1 Introduction

High order Weighted Essentially Non–Oscillatory (WENO) methods are ideally

suited for high fidelity simulations of complex physics containing discontinuities.

An example of a typical application is a canonical simulation of sound generated by

a shock-vortex interaction [1,2]. The high-order nature of WENO efficiently resolves

the subtle details of the sound generation and propagation, while the stencil biasing

mechanics ensures robust, high resolution properties in the vicinity of the shock.

Although ideal for infinite or periodic domains, conventional WENO formulations encounter serious challenges at nodes that are near domain boundaries. To

illustrate this, first realize that all high-order finite-difference (HOFD) formulations

use “inward biased” stencils near boundaries that maintain the accuracy and stability of the interior scheme. Next, note that WENO schemes invoke stencil biasing

mechanics throughout the entire domain, including nodes adjoining the boundaries

(e.g. fourth-order WENO schemes test three candidate stencils at each node). Thus,

not only must the boundary closures used for WENO schemes be 1) inward-biased,

2) stable and 3) accurate, they must also be stable for any possible combination

of candidate stencils occurring anywhere in the domain. Simultaneously satisfying

these constraints is a remote possibility if left strictly to chance.

Many ad hoc boundary closures have been proposed for conventional WENO

schemes, but most implementations still rely on low order boundary closures to

achieve stability and robustness. For flow near a wall, Shen et al. [3] constructed the

domain such that the no slip condition is enforced at the WENO flux point. Pressure

was extrapolated to the wall using a third order scheme, and the next WENO flux

cell was treated with a third order reconstruction. This procedure reduces the global

order of accuracy of their 5th order WENO scheme. Alternatively, ghost points can

be used in combination with physical boundary conditions to specify data in stencils

that extend outside of the domain [4]. This reduces the generality of the boundary

2
treatment. Another strategy is to simply discard stencils that extend outside of the

domain [4].

The principle objective of this work is to develop boundary closures to complement the periodic domain Energy Stable WENO finite difference methodology

developed and reported in Refs. [5] and [6]. The new interior/boundary ESWENO

schemes (henceforth referred to as “finite-domain ESWENO”) retain all the salient

features of the original periodic schemes, including: 1) conservation and L2-energy

stability for constant coefficient (linear) hyperbolic systems, including those with

discontinuous initial or boundary data, 2) design order accuracy throughout the

entire domain, especially regions near the boundaries or near smooth extrema, and

3) full stencil biasing mechanics (Left, Central, Right) at all possible points.

The finite-domain ESWENO schemes are constructed by adding a “special” nonlinear artificial dissipation term to a conventional WENO scheme (e.g. those found

in Refs. [7, 8]). The additional term is design-order accurate for smooth solutions,

including smooth extrema, and is constructed such that the resulting ESWENO

scheme is stable in the L2-energy norm for both continuous and discontinuous solutions.1

In this initial work, we focus exclusively on the central fourth-order class

of ESWENO schemes, leaving the development of 6th- and 8th-order schemes to a

later date.

The strategy for developing finite-domain ESWENO schemes is for the most

part equivalent to that used to build the periodic domain ESWENO schemes. The

essential elements are summarized as follows:

• Develop a 4th-order finite-domain target scheme that is stable, conservative

and accurate for smooth flows. This task is accomplished using the summationby-parts (SBP) [6, 9] framework (as was done in the periodic case). The SBP

framework requires the construction of accurate discrete schemes that satisfy

basic matrix conditions, that if met guarantee discrete stability and conservation.

• The target scheme is recast into the “dual-grid” framework of the conventional

WENO approach, whereby the solution is stored/advanced at the gridpoints,

while interface fluxes are constructed at the “half-points”.

• The stability of the resulting stencil is tested, and if unstable, is modified with

a special artificial dissipation term. This ensures that the discrete eigenvalues

of the spatial scheme are bounded within the left-half of the complex plane.

The organization of the paper is as follows. Section 2 introduces a dual-grid

nomenclature and summarizes the necessary SBP formulas including a discrete stability proof. Section 3 summarizes the implementation of the ESWENO scheme,

and then expresses the baseline ESWENO operator in SBP form. Section 4 formulates the artificial dissipation term necessary for the ESWENO approach. Section

1The conventional WENO scheme, does not have such an energy estimate. Thus, while the

spectrum of the symmetric part of the ESWENO scheme is always located in the left half of the

complex plane, the symmetric part of the WENO scheme could have positive eigenvalues.

3
5 presents a summary of numerical test cases. Conclusions and directions are presented in Section 6. An appendix is included with all the implementation details.

2 Definitions

2.1 Complementary Grids

The implementation of the ESWENO (WENO) stencil biasing mechanics, as a matter of convenience, uses two complementary sets of gridpoints, differing in dimension

by one. As such, two discrete sets of points

xN = [x1, x2, · · · , xN−1, xN ]

T

; x ̄M = [ ̄x0, x ̄1, · · · , x ̄N−1, x ̄N ]

T

are defined on the interval 0 ≤ x ≤ 1. Furthermore, a projection is defined of any

function φ onto xN and x ̄M respectively as

φ = [φ(x1), φ(x2), · · · , φ(xN−1), φ(xN )]T; φ ̄ = [φ( ̄x0), φ( ̄x1), · · · , φ( ̄xN−1), φ( ̄xN )]T

.

The xN is a uniform set of “solution points” and stores the discrete solution u.

Surrounding each solution point, xN j, a control volume is constructed, whereby the

governing conservation laws are enforced. A second set of “flux points”, x ̄M, is situated at the interfaces between adjoining control volumes, where data interpolated

from u is used to construct interface fluxes  ̄f and establish a simple flux balance.

Figure (1) illustrates the relationship between the solution-mesh xN and the

flux-mesh x ̄M.

f+ f−

x1 x2 x3 x4 xN−3 xN−2 xN−1 xN

xi

ˆfi+1/2

fi−1 fi fi+1 fi+2

SL SC SR

xf1 xf2 xf3 xf4 xfi xfN−4 xfN−3 xfN−2 xfN−1

xf0 xfN

xi+1

ˆfi−1/2

Figure 1. The computational domain is depicted. The ESWENO solution and flux

points, and the smoothness indicator stencils S

L, SC and SR are shown.

The biasing mechanics move data via interpolation or differentiation between

the two sets of points using non-square matrix operations. For example, a flux

constructed on the high-dimensional set from data on the low-dimensional set, is

represented as

 ̄f = MIN f (1)

where M = N + 1 and MIN is an [(N + 1) × N] matrix.

4
2.2 Summation-By-Parts Operators

Consider a linear, scalar wave equation

∂u

∂t +

∂f

∂x = 0, f = au, t ≥ 0, 0 ≤ x ≤ 1,

u(x, 0) = u0(x) ; u(0, t) = g(t)

(2)

where u is the continuous solution, a is a constant and u0(x) is a bounded piecewise

continuous function. Without loss of generality, assume that a ≥ 0 on the interval

0 ≤ x ≤ 1. Applying the energy method to equation (2) leads to

d

dtkuk

2

L2 + a[u

2

(1, t) − g

2

(t)] = 0 (3)

where k·kL2is the continuous L2 norm. Thus, equation (2) admits an energy solution

that is strictly dissipative (if g(t) = 0, a ≥ 0) for all time. We now develop using

mimetic techniques (see Ref. [6] or Ref. [10]) a class of discrete spatial operators for

which the discrete energy is bounded from above by the continuous target estimate

provided by equations (2) and (3).

Define the low dimensional space xN to be a uniform grid xj+1 = j∆x, j =

0, · · · , N − 1, with ∆x = 1/(N − 1). On this grid, we define a flux f = au and

its derivative fx = aux, where u = [u(x1, t), . . . , u(xN , t)]T and ux = [ux(x1, t), . . . ,

ux(xN , t)]T are projections of the continuous solution u and its derivative ux onto the

computational grid xN . Next, define a fourth-order approximation for the first-order

derivative term in equation (2) as

∂f

∂x = Df + O(∆x

4

) (4)

Placing a mild restriction on the generality of the derivative operator (see Ref. [10]

or Ref. [6]), the matrix D is expressed in the following form:

D = P

−1

[Q + R] ; Q + QT = Diag[−1, 0, · · · , 0, 1]

R = RT; v

T Rv ≥ 0

P = P

T

; v

TPv > 0

(5)

for any real vector v 6= 0. The matrix R is defined by the expression

R = Λ0 + D1Λ1[D1]

T + D1[D1]TΛ2D1[D1]T + D1[D1]T D1Λ3[D1]T D1[D1]T

(6)

which allows a constructive means of forming a symmetric semi-definite discrete

operator. Herein we define the matrices Λi to be diagonal positive semidefinite

matrices of the appropriate size: i.e.

Λe = Diag[λ1, · · · , λN ] ; λj ≥ 0 , j = 1, N ; e = 2p, p = 0, 1

Λo = Diag[λ0, · · · , λN ] ; λj ≥ 0 , j = 0, N ; o = 2p + 1, p = 0, 1

5
Using the definitions of the SBP operators [eqs. (4)–(6)], and defining on xN

the discrete approximate solution u and flux f = au, the semi-discrete counterpart

of equation (2) becomes

∂u

∂t + P

−1Q f =

−P−1

h

Λ0 + D1Λ1[D1]

T + D1[D1]TΛ2D1[D1]T + D1[D1]T D1Λ3[D1]T D1[D1]T

i

f

(7)

To show that the above finite-difference scheme is stable, the energy method is

used. Multiplying equation (7) with u

TP, adding the result to its transpose, and

rearranging the expression using equation (5), yields the semi-discrete energy equation2

d

dtkuk

2

P + a[u

2

(N, t) − g

2

(t)] = −2a

X

q

i=0



[D1

i

]

Tu

T

Λi[D1

i

]

Tu ≤ 0 (8)

The right-hand side of equation (8) is non-positive because the diagonal matrix

Λn is positive semidefinite [v

TΛnv ≥ 0 for all real v of length N] and a ≥ 0; thus

the stability of the finite-difference scheme given by equation (7) is assured. This

result can be summarized in the following theorem:

Theorem 1 The approximation [eq. (7)] of the problem [eq. (2)] is stable if equations [(4)–(6)] hold.

Remark. The finite-difference scheme [eq. (7)] constructed as an approximation

of equation (2) is inherently nonlinear, despite the linearity of (2). Indeed, the

matrices Q and Λn are nonlinear matrices [i.e., Q = Q(u) and Λn = Λn(u)].

Remark. Although the matrices Q and Λn are constrained in form by Q + QT =

Diag[−1, 0, · · · , 0, 1] and v

TΛnv ≥ 0, no specific assumptions are made about the

values of the coefficients.

Remark. All matrices present in equation (7) need not be square. Indeed, it is

advantageous to define the [N × (N + 1)] dimensional matrix D1 to be

D1 =





−1 1 0 0 0 0

0 −1 1 0 0 0

0 0 .

.

.

.

.

. 0 0

0 0 0 −1 1 0

0 0 0 0 −1 1





(9)

to properly account for the difference between the dimension of the xN and x ̄M

grids.

2.2.1 Central Fourth-order SBP Operators

The target operator used herein for the ESWENO scheme is a fourth-order nondissipative (i.e. R = 0) central SBP operator, for which equation (5) reduces

2The nomenclature D1i

, i = 1, q used in equation (8) is not valid for nonsquare matrices. Herein,

define D1

2p

to mean {D1[D1]

T

}

p

; p = 1, q.

6
to D = P

−1Q subject to the constraints Q + QT = Diag[−1, 0, · · · , 0, 1] and

P = P

T

; v

TPv > 0.

The matrices P and Q are composed of the conventional banded operators in

the interior of the domain:

PI = ∆x I ; QI = Pentadiagonal

1

12

,

−8

12

, 0,

8

12

,

−1

12 

.

Boundary closures satisfying the SBP convention require auxiliary stencils at Nb

points near each boundary. Elsewhere, it is shown [11] that Nb ≥ 4 is required to

satisfy the necessary accuracy, symmetry and skew-symmetry constraints. Herein,

the Nb = 4 is assumed, whereby the matrices P and Q take the form

P = ∆x





P0 0 0

0 I 0

0 0 (P0)

P T



 ; Q =





Q0 Qd 0

−QT

d QI Qd

0 −QT

d − (Q0)

P T



 (10)

with

P0 =





p11 p12 p13 p14

p12 p22 p23 p24

p13 p23 p33 p34

p14 p24 p34 p44





; Q0 =





−1

2

q12 q13 q14

−q12 0 q23 q24

−q13 −q23 0 q34

−q14 −q24 −q34 0





;

Qd =





0 · · ·

0 · · ·

−1

12 0 · · ·

8

12

−1

12 0 · · ·





(11)

where P T denotes the per-symmetric transpose.

Two types of P-norms are typically used to prove stability with SBP discretizations: diagonal, and block. The diagonal-norm is advantageous when considering

equations of more complexity than the constant coefficient, linear equation considered herein. Indeed, SBP proofs using the diagonal-norm extend to variablecoefficient problems [12], as well as certain nonlinear formulations [13,14], a property

not shared by the block-norm boundary closures. Diagonal-norm boundary closures

can not preserve the full design-order potential of the interior operator. Closures

satisfying the matrix constraints in equation (8) with a fourth-order interior scheme,

can have accuracy no higher than second-order. Thus, the global accuracy of this

combination is only third-order [15, 16], or in general p + 1-order despite 2p-order

accuracy in the interior. (The nomenclature p − 2p − p; p = 1, · · · is sometimes used

to describe the diagonal-norm formulations. See Ref. [9] and Ref. [11] for further

discussion.)

Although block-norm closures are less robust than their diagonal-norm counterparts, for variable-coefficient linear or nonlinear equations, they enable the use

of 2p − 1-order stencils, and preserve the formal accuracy of the 2p-order interior

operator. For example, third-order boundary closures preserve the formal accuracy

of the fourth-order interior scheme. Furthermore, the additional parameters in the

7
P-norm (six, in the (4 × 4) case) can be used to optimize other properties of the

third-order closure. Indeed, two free parameters (herein denoted θl, l = 1, 2) remain,

after satisfying all accuracy, and SBP matrix constraints, when using a (4×4) block

norm to close a fourth-order interior. Both parameters, as will be shown in the

next section, are needed to develop stable, and accurate fourth-order ESWENO

boundary closures.

3 Expressing Summation-By-Parts Operators in Flux

Form

The next objective is to recast the target fourth-order central SBP operator (D =

P

−1Q defined on the xN grid) into the “dual-grid” framework of the conventional

ESWENO/WENO approach. Of particular interest is extending the ESWENO’s

“interior” stencil biasing mechanics, to those stencils at or near domain boundaries.

We begin with a cursory description of the fourth-order ESWENO scheme, followed

by its representation in matrix form.

3.1 General Structure of ESWENO Schemes

The conventional high-order WENO/ESWENO finite difference scheme for the scalar

1-D wave equation (Eq. (2)) can be written component-wise in the following semidiscrete form:

duj

dt +

 ̄fj+ 1

2

−  ̄fj− 1

2

∆ ̄xMj

= 0 (12)

The numerical flux  ̄fj+ 1

2

for the central fourth-order WENO scheme is given by

the expression

 ̄fj+ 1

2

=  ̄w

L

j+1/2

 ̄f

L

j+1/2 +  ̄w

C

j+1/2

 ̄f

C

j+1/2 +  ̄w

R

j+1/2

 ̄f

R

j+1/2

, (13)

where  ̄f

(r) = MI

(r)

N

f, r = {L, C, R} are second order fluxes constructed from data

interpolated using three distinct stencils S

L, SC, and SR.

3

The nonlinear weight functions  ̄w

(r)

j+1/2

, r = {L, C, R} used in equation (13),

are constructed such that the interface flux  ̄fj+ 1

2

is a convex combination of the

three candidate second-order fluxes. The weight functions are constructed using the

expressions

w ̄

(r)

j+1/2 =

α ̄

(r)

X

j+1/2

r

α ̄

(r)

j+1/2

;  ̄α

(r)

j+1/2 =  ̄d

(r)

j+1/2



1 +

τ ̄j+1/2

 + β ̄

(r)

j+1/2



 , r = L, C, R (14)

3Figure (1) shows a schematic of the interpolation process. Each second-order flux is constructed from data obtained using only two points: S

L = {xj−1, xj}, SC = {xj , xj+1}, and

S

R = {xj+1, xj+2}. The r = {L, C, R} nomenclature identifies the origin of the data relative to

the interface position. The interpolates need not be limited to three; indeed S

LL = {xj−2, xj−1},

and S

RR = {xj+2, xj+3} can be included near boundaries to increase the fidelity of the boundary

closures.

8
where  ̄d

(r)

j+1/2

is the target weight of the candidate stencil. If no stencil biasing occurs

(a well resolved solution), the target central difference operators are recovered, and

the functions  ̄w

(r) used as weights in equation (13) reduce to

ω ̄

(r)

j+1/2 =

 ̄d

(r)

X

j+1/2

r

 ̄d

(r)

j+1/2

=  ̄d

(r)

j+1/2

.

Definitions of the parameters  ̄τj+1/2, , β ̄

(r)

j+1/2

, r = L, C, R, as well as implementation

details are included in the appendix.

The flux derivative expressed in equation (12) in component form can be written

as

 ̄fj+ 1

2

−  ̄fj− 1

2

∆ ̄xMj

→



[δx ̄M]

−1D1

 ̄f



j

(15)

where [δx ̄M] is the diagonal matrix [δx ̄M] = Diag[D1x ̄M] that accounts for variable

cell sizes. Further expanding the flux function  ̄f in equation (15) in terms of matrix

interpolation operators, and equating the result with the target SBP form, yields

the expression

[δx ̄M]

−1D1

 ̄f = [δx ̄M]

−1D1



ω ̄

L

(

 ̄d) MI

L

N +  ̄ω

C(

 ̄d) MI

C

N +  ̄ω

R(

 ̄d) MI

R

N



f = P

−1Qf

(16)

The solution to equation (16) is the central objective of this section; specifically, to

express a target fourth-order SBP discretization P

−1Q, in terms of the conventional

1) interpolation operators MIN , 2) weight functions  ̄ω, and 3) difference operator

D1 needed to implement the ESWENO formulation.

Remark. There is no guarantee that a solution exists to Equation (16). The

conventional periodic domain operators MI

r

N , ω ̄

r

; r = L, C, R on uniformly spaced

grids xN , x ̄M, match the fourth-order centered P

−1Q identically in the interior of

the domain. The task at hand is to determine boundary operators for MI

r

N , ω ̄

r

; r =

L, C, R and the grids xN , x ̄M that match P

−1Q at boundaries.

Remark. Equation (16) is a nonlinear function of the flux grid x ̄M. Indeed,

assuming a uniform solution grid xN , the functional dependency of the matrix operators are

h

d ̄(r) = d ̄(r)(x ̄M) ; MI

(r)

N = MI

(r)

N

(x ̄M)

i

, r = L, C, R ;

P = P(θl) , Q = Q(θl) , l = 1, 2

3.2 Nonuniform Flux-Points

An analytic solution to equation (16) is not immediately forthcoming using Mathematica 7.0 [17], if the solution xN and flux x ̄M points are uniformly spaced on the

domain 0 ≤ x ≤ 1:

xj+1 = j∆x , j = 0, N − 1 , ∆x = 1/(N − 1) ;

x ̄k =

(xk+xk+1)

2

, k = 1, N − 1 , x ̄0 = 0, x ̄N = 1,

9
for either the diagonal-norm or the block-norm SBP operators P

−1Q. Increasing the

boundary block sizes to (6×6), while providing numerous additional free parameters

θl, does not remove an apparent inconsistency in the constraint equations.

Multiple solutions exist for equation (16), for either the diagonal- or block-norms,

if the second, third, and fourth flux points x ̄M are shifted to a nonuniform set of

points near the boundaries. (Figure (1) illustrates the relationship between the

solution-mesh xN and flux-mesh x ̄M, including the shifted flux points near each

boundary.) The location of the flux points is precisely related to the P-norm used to

prove stability for the SBP operator. The following lemma is a central contribution

of this work:

Lemma 1 Assume N uniformly distributed solution points xN . There exists an N+

1 dimensional set of flux points x ̄M for which the central, fourth-order, ESWENO

interpolation operators, can be expressed in SBP form. The discrete set of points

x ̄M satisfies the following expression:

P

−1D1x ̄M = 1 (17)

with 1 the unit vector.

Proof: The proof is constructive. Given x ̄M, the general solutions are determined

using Mathematica. Owing to the complexity of the general expressions, only two

optimized solutions are included herein: one derived using a diagonal-norm, the

other a block-norm. Note that D1x ̄M = P1 is an equivalent statement because P

is nonsingular. 

Although a mathematical motivation for Lemma (1) has not been found, a dimensional argument can be used to explain the relationship between P and x ̄M.

First, note that the D1 matrix is an undivided difference, and does not contain

any reference to the precise location of the flux-points  ̄f, nor to their separation;

D1 is dimensionless and simply forms the difference between the left and right cell

interface states. The dimensionality (or inverse weight) in the derivative operation is provided by the term [δx ̄M]

−1

. Left multiplying the flux difference D1f

by 1

T

[δx ̄M] (the matrix equivalent of integration), yields the domain flux balance

1

T

[δx ̄M]D1

 ̄f = f(1) − f(0). Clearly, the individual terms in the matrix [δx ̄M] are

the quadrature weights needed for discrete integration. Thus, the matrix [δx ̄M] is

filling precisely the same role in equation (16), as the P matrix in the SBP context

D = P

−1Q. It is not surprising that P and [D1x ̄M] are closely related.

The next two examples demonstrate the validity of Lemma (1) as well as provide the flux points x ̄M needed for the diagonal- and block-norm fourth-order SBP

operators. First, consider the [N × N] diagonal P-norm Pd defined by

Pd = ∆x Diag 

17

48

59

48

43

48

49

48 1 · · · 1

49

48

43

48

59

48

17

48 

Distribute the flux points x ̄M as

x ̄M =



0, x ̄1, x ̄2, x ̄3,

7∆x

2

, · · · ,



1 −

7∆x

2



,(1 − x ̄3),(1 − x ̄2),(1 − x ̄1), 1

T

(18)

10
where

x ̄1 =  ̄x0 +

17∆x

48

;  ̄x2 =  ̄x1 +

59∆x

48

;  ̄x3 =  ̄x2 +

43∆x

48

;  ̄x4 =  ̄x3 +

49∆x

48

=

7∆x

2

.

Then, the vector D1x ̄M is precisely the reciprocal of the diagonal element of Pd.

The product P

−1

d D1x ̄M is trivially the identity vector 1, [the condition of Lemma

(1)] since Pd is a diagonal matrix.

Lemma (1) generalizes immediately to block-norm P cases. Recall that P is not

uniquely defined in the fourth-order, centered, block-norm case. The matrix P is a

function of two independent parameters (i.e. P = P(θl), l = 1, 2), after all accuracy,

and matrix constraints are satisfied. In this specific case (but not in general) the

matrix-vector product P1, used in Lemma (1) to determine x ̄M, is independent of

the parameters θl. Thus, the following set of points x ̄M

x ̄M =



0,

43∆x

144 ,

61∆x

36 ,

349∆x

144 ,

7∆x

2

, · · · ,



1 −

7∆x

2



,



1 −

349∆x

144 

,



1 −

61∆x

36 

,



1 −

43∆x

144 

, 1

T

(19)

satisfied the conditions of Lemma (1), independent of the values of the available

parameters θl, l = 1, 2.

Remark. Existence of a set of points x ̄M does not imply that x ̄M is unique,

although no other set has been identified.

Remark. The solution- (xN ) and flux- (x ̄M) grids are interdigitated, the exceptions being the collocated endpoints. Given this set of flux points x ̄M, the interpolation stencils S

r

, r = {L, C, R} are consistently defined for all points  ̄xMj

, j = 2, N−2.

The flux-points next to each boundary  ̄xMj

, j = 1 and j = N −1 only have a centered

and an inward-biased stencil. The first and last flux-points, in addition to inward

biased interpolation stencils, have an “exact” interpolant derived from the collocated solution-point. This property is shared with conventional WENO boundary

operators.

3.3 Interpolation Operators and Preferred Weights

An analytic solution to equation (16) exists for either the diagonal-norm or the

block-norm case, when x ̄M is given by equation (18) or (19), respectively. The

diagonal-norm case required additional interpolants at the boundaries, to become

realizable. Specifically, the interpolation operators I

LL, IRR, and weights  ̄wLL, w ̄RR,

were added in equations (13) and (16). The additional terms were only required at

the boundaries.

Interpolation operators for the WENO/ESWENO3-4-3 and WENO/ESWENO2-4-2

operators are given in the appendix.

4 Energy Stabilization Terms

When the nonlinear weights used in equation (13) satisfy the target conditions,

w ̄

(r) =  ̄ω(r)

, r = {L, C, R}, then all stencils are stable without the dissipation term

P

−1R. With nonlinear weights, however, there is no guarantee that the baseline

11
WENO operator D is stable. To ensure stability, an artificial dissipation term Dad is

added to the WENO operator D, such that the symmetric portion of the resulting

(ESWENO) operator D ̄ is positive-definite. Specifically, D ̄ satisfies the following

matrix conditions: D ̄ = D + Dad subject to Dsym + Dad ≥ 0.

The matrix Dad is constructed by expanding the symmetric matrix Dsym into

its elemental components using the decomposition4

Dsym = P

−1

h

Λ0 + D1Λ1[D1]

T + D1[D1]TΛ2D1[D1]T + D1[D1]T D1Λ3[D1]T D1[D1]T

i

.

The existence of this decomposition is established in Ref. [6].

Next, the diagonal terms of Λi are modified to be smoothly positive using the

expression:

λ ̄

(i)

j =

1

2

 r

λ

(i)

j

2

+ δ

2

i − λ

(i)

j

!

(20)

where δ is dependent on the grid spacing.

δ1 = (∆ξ)

3

δ2 = (∆ξ)

2

δ3 = (∆ξ)

2

(21)

The artificial dissipation operator takes the final form

Dadf = P

−1Rf =

P

−1

h

Λ ̄

0 + D1Λ ̄1[D1]

T + D1[D1]TΛ ̄

2D1[D1]

T + D1[D1]T D1Λ ̄

3[D1]

T D1[D1]T

i

f

(22)

The dissipation matrix Λ0 (at least herein) is identically zero for both the diagonal and block norm schemes; thus it is convenient to define the vector

ψ ̄ =

h

Λ ̄

1[D1]

T + [D1]TΛ ̄

2D1[D1]

T + [D1]T D1Λ ̄

3[D1]

T D1[D1]T

i

f (23)

This definition of ψ ̄ allows the artificial dissipation term to be combined with the

original flux  ̄f given by equation (13) to yield the following:

df

dx = P

−1D1



 ̄f + ψ ̄



+ O(∆x

4

) (24)

5 Numerical Tests

The accuracy and stability characteristics of the new ESWENO framework are tested

using the one dimensional advection equation and the quasi-one dimensional and

two dimensional Euler equations. As a more practical test, the framework was

implemented in a two dimensional reacting Navier-Stokes solver and was used to

simulate supersonic non-reacting and reacting shear layers. All simulations are

integrated in time or pseudo-time using the same five-step, fourth-order Runge Kutta

solver [18].

4Note that the dimensions of the even indexed Λi are (N × N), while the odd indexed Λi have

dimensions (N + 1 × N + 1). This is the result of the difference matrix D1 having dimensions

[N × (N + 1)].

12
5.1 Eigenvalues of the ESWENO operator

To demonstrate that the finite-domain ESWENO operator is positive semidefinite,

an eigenvalue analysis is conducted using the one dimensional scalar advection equation with discontinuous initial data. We begin by rearranging equation (2) into the

form

P

d ̄u

dt = −a (Q + R)  ̄u (25)

Setting  ̄u = ˆue

kt (see Ref. [19]) in equation (25) yields

− λPˆu = (Q + R) ˆu, λ =

k

a

(26)

Stability is guaranteed if all real eigenvalues satisfy the condition λr ≤ 0 since a > 0.

ESWENO/WENO operators depend on the calculated solution, so eigenvalues were

calculated over ten characteristic times using N = 100 grid cells to ensure consistency of the analysis. The boundary condition used in the simulation was g(t) = 0.

Eigenvalues were computed for the ESWENO/WENO3-4-3 scheme applied to a

convected square wave. The resulting maximum real eigenvalues are plotted against

time in Figure (2).

Figure 2. Maximum real eigenvalues for a convected square wave are plotted as a

function of time for the WENO3-4-3 operator with and without the energy stable

addition.

It is clear that the bounded WENO3-4-3 operator has eigenvalues in the right

half of the complex plane. This is fixed by the addition of the energy stable term in

the ESWENO3-4-3 scheme. This result is similar to that observed in the periodic

case [5, 6].

13
5.2 Linear Wave Equation

Numerical solutions of the linear advection equation were examined to test the accuracy and dissipation characteristics of the finite-domain ESWENO schemes. The

convergence rate of the ESWENO2-4-2 and ESWENO3-4-3 operators are calculated

for a smooth solution. To evaluate the robustness of the schemes, a wave solution

with a time-dependent, discontinuous inlet condition is examined. It is important

to note that since the ESWENO3-4-3 operator is constructed using a block P-norm,

simply assigning the boundary solution to be the boundary data (injection boundary conditions) does not maintain energy stability. Instead, both ESWENO schemes

utilize the Simultaneous Approximation Term (SAT) penalty method to enforce the

boundary conditions [20].

5.2.1 Sine Wave

To test the global order of convergence, the ESWENO scheme is used to approximate

the solution to the linear advection equation with the exact solution

u(x, t) = sin(x − at), a = 1, x ∈ (−π, π) (27)

Ten characteristic times were simulated. The L2 and L∞ errors are tabulated along

with the convergence rate in Tables (1) and (2), respectively, for the full block-norm

ESWENO operator and are compared to the errors from the block-norm linear centered difference operator. For the diagonal norm operators, the errors are tabulated

in Tables (3) and (4). It is clear that the ESWENO operator yields the design order of convergence for the linear advection equation, and provides almost the same

accuracy as the linear counterpart.

Table 1. L2 error and convergence rate of a smooth solution to the 1D advection

equation for the linear and ESWENO operators with the block norm boundary

closure.

Linear Block Norm ESWENO3-4-3

Number of Cells L2 Error L2 Rate L2 Error L2 Rate

100 1.60e-06 - 1.65e-06 -

200 1.01e-07 3.99 1.04e-07 3.99

400 6.31e-09 4.00 6.55e-09 3.99

800 3.95e-10 4.00 4.10e-10 4.00

1600 2.52e-11 3.97 2.61e-11 3.97

5.2.2 Square Wave

A square wave advection problem with the exact solution

u(x, t) = tanh (b (x − at + 2π)) − tanh 

b



x − at +

3π

2



2

, a = 1, x ∈ (−π, π)

(28)

14
Table 2. L∞ error and convergence rate of a smooth solution to the 1D advection

equation for the linear and ESWENO operators with the block norm boundary

closure.

Linear Block Norm ESWENO3-4-3

Number of Cells L∞ Error L∞ Rate L∞ Error L∞ Rate

100 3.24e-06 - 3.43e-06 -

200 2.07e-07 3.97 2.18e-07 3.97

400 1.30e-08 3.99 1.38e-08 3.99

800 8.17e-10 3.99 8.65e-10 3.99

1600 5.20e-11 3.97 5.51e-11 3.97

Table 3. L2 error and convergence rate of a smooth solution to the 1D advection

equation for the linear and ESWENO operators with the block norm boundary

closure.

Linear Diagonal Norm ESWENO2-4-2

Number of Cells L2 Error L2 Rate L2 Error L2 Rate

100 4.20e-05 - 2.05e-05 -

200 5.21e-06 3.01 2.30e-06 3.16

400 6.48e-07 3.01 2.72e-07 3.08

800 8.09e-08 3.00 3.30e-08 3.04

1600 1.01e-08 3.00 4.07e-09 3.02

was tested, where b = 200 caused a very sharp gradient that could not be resolved.

The numerical solution was calculated for up to t = 10.0, and N = 200 grid cells were

used. The initial condition was u0(x, 0) = 0 in the domain. The wave was passed

into the domain, advected through the interior, and then passed out of the domain.

This problem tests the dissipation and oscillatory characteristics of the boundary

treatment as well as the interior treatment. A comparison of the numerical solution

to the exact solution is shown for the two different boundary norms in Figure (3)

at two different times. The solutions shown in Figure (3) prevent oscillations, and

smearing of the discontinuity caused by the boundary treatment, particularly at

the inlet, is reasonable. The solution using the ESWENO2-4-2 operator in the left

figure (t = 2.0) exhibits slightly more dissipative error than the solution using the

ESWENO3-4-3 operator.

5.3 Euler Equations

The quasi-one dimensional Euler equations given by

∂U

∂t +

∂F

∂x =

Ax

A

G

U =





ρ

ρu

ρE



 F =





ρu

ρuu + p

(ρE + p)u



 G =





ρu

ρuu

(ρE + p)u





(29)

15
Table 4. L∞ error and convergence rate of a smooth solution to the 1D advection

equation for the linear and ESWENO operators with the block norm boundary

closure.

Linear Diagonal Norm ESWENO2-4-2

Number of Cells L∞ Error L∞ Rate L∞ Error L∞ Rate

100 7.79e-05 - 6.60e-05 -

200 9.96e-06 2.97 8.16e-06 3.02

400 1.25e-06 3.00 1.07e-06 2.93

800 1.57e-07 3.00 1.37e-07 2.97

1600 1.96e-08 3.00 1.73e-08 2.99

Figure 3. Numerical solutions using the ESWENO2-4-2 and ESWENO3-4-3 operator compared to the exact solution of an advected square wave at two different

times.

were used to test the extension of the ESWENO operators to non-linear hyperbolic

systems of equations. In order to properly calculate stencil biasing parameters,

the conservative fluxes in equation (29) are transformed into characteristic form

and split into forward and backward propagating waves using global Lax-Friedrichs

splitting.

f

± =

1

2

(Fc ± ΛmaxUc) (30)

where Fc and Uc represent the characteristic fluxes and variables, respectively. The

transformation from conservative to characteristic variables is performed at each flux

point using transformation matrices constructed from the Roe averaged variables at

the flux point. The ESWENO reconstruction is used to calculate interpolated fluxes,

and then the interpolated fluxes are transformed back to physical space, where the

gradient is calculated.

The solutions to the quasi-1D Euler equations utilize the SAT penalty boundary

16
conditions in order to ensure stability of the solution and consistency of the boundary

closure. The details of the SAT penalty method for the Euler and Navier-Stokes

equations can be found in Sv ̈ard, Carpenter, and Nordstr ̈om [21].

5.3.1 Subsonic One-Dimensional Nozzle

The first nonlinear test problem was a steady state, subsonic flow in a one-dimensional

nozzle. The area was defined as

A(x) = 1 − 0.8x(1 − x), x ∈ (0, 1) (31)

The inlet Mach number was set to M = 0.5 and the outlet pressure was set to

be equivalent to the inlet. The exact solution is plotted in Figure (4). The L2

error norm of the density and the convergence rate for both ESWENO schemes is

tabulated in Table (5). The calculated convergence rates match the design order.

Remark. The design order of convergence is achieved despite the fact that the

transformation matrices are calculated from second order accurate Roe averages.

This is because the transformation is only used to transform to wave space and

back to physical space in order to facilitate flux splitting. The accuracy of this

transformation does not affect the order of accuracy of the interpolation. Instead,

the accuracy of the transformation relates to how accurately the physical fluxes can

be decomposed into waves.

Table 5. L2 error and convergence rate of a smooth solution to the quasi-one dimensional Euler equations for ESWENO2-4-2 and ESWENO3-4-3.

ESWENO2-4-2 ESWENO3-4-3

Number of Cells L2 Error L2 Rate L2 Error L2 Rate

25 6.55E-05 - 2.13E-06 -

50 6.14E-06 3.42 1.14E-07 4.22

100 8.05E-07 2.93 6.43E-09 4.15

200 1.03E-07 2.97 3.82E-10 4.07

5.3.2 Transonic One-Dimensional Nozzle

The next test of the boundary-closed ESWENO scheme was a steady state transonic

one-dimensional nozzle. The cross sectional area of the transonic nozzle is

A(x) = 1.398 + 0.347 tanh(0.8x − 4), x ∈ (0, 10) (32)

The inlet Mach number was set to M = 1.5, and the shock through the overexpanded nozzle occurred at x = 5.0. The resulting density profile from each

ESWENO scheme is compared to the exact solution in Figure (5). The solution

is non-oscillatory, and spreads the shock over a few cells. The results from the different boundary treatments do not differ appreciably, which can be attributed to

the zero gradient at the boundary.

17
Figure 4. Solution for the quasi-1D subsonic nozzle problem is shown.

Remark. The transonic one-dimensional nozzle is a problem where boundary

effects do not significantly influence the solution, but the problem demonstrates the

stability of the ESWENO boundary closure for a steady state problem. The behavior

of the solution to this problem in the current work exhibits the same characteristics

that were observed in the periodic case [5]. Similarly, Sod’s shock tube problem and

the sine-shock wave interaction problems exhibited the same characteristics. The

results for those problems are not repeated here.

5.3.3 Two Dimensional Inviscid Vortex

An inviscid two-dimensional vortex convection problem was used to test the accuracy of the finite-domain ESWENO framework in multiple dimensions. The exact

solution to the convected vortex is

f(x, y, t) = 1 −



(x − x0 − U∞t)

2 + (y − y0)2



,

ρ(x, y, t) = 1 − 

2

γ − 1

8π

2

exp (f(x, y, t)) 1

γ−1

, p =

ρ

γ

γ

u(x, y, t) = U∞ − 

y − y0

2π

exp 

f(x, y, t)

2



,

v(x, y, t) = 

x − x0 − U∞t

2π

exp 

f(x, y, t)

2



,

U∞ = M∞c∞, x ∈ (0, 10), y ∈ (−5, 5), (x0, y0) = (5, 0), t ≥ 0,

(33)

where  = 1.0, M∞ = 0.5, and γ = 1.4 were used. Density contours from the initial

condition and the approximate solution at t = 10 for the ESWENO3-4-3 operator

with N = 200 × 200 grid cells are shown in Figure (6).

18
Figure 5. Density is plotted for Sod’s shock tube problem at t = 0.2 using N = 100

cells.

Both boundary closures were tested for the two-dimensional vortex and compared to corresponding linear operators. The accuracy is affected not only by the

boundary closure of the finite difference scheme, but also by the treatment of the

boundary conditions. In this case, the SAT penalty method is used to specify

the exact boundary data in a well posed manner. This is meant to limit reflections from the boundaries to be influenced only by the boundary closure and the

SAT penalty term. A convergence study would not yield meaningful results using

more common boundary condition treatments, such as Navier-Stokes Characteristic

Boundary Conditions.

The L2 and L∞ convergence rates for the ESWENO3-4-3 operator are tabulated

in Tables (6) and (7), respectively. The convergence rates for the ESWENO2-4-2

operator are given in Tables (8) and (9). The results show that the linear operators

recover the design order accuracy. The ESWENO operators yield better than design

order convergence rates for lower resolutions, but achieve design order convergence

as resolution increases. The ESWENO operators exhibit larger errors for lower

resolution than the corresponding linear operator. However, as the critical regions

are better resolved, the linear operators and ESWENO operators yield similar L∞

errors and nearly identical L2 errors. This is the expected behavior, since the nonlinear weights in the ESWENO operator approach the target weights as resolution

is improved, recovering the linear operator.

5.4 Supersonic Shear Layers

A more practical test of ESWENO is the two dimensional non-reacting and reacting supersonic hydrogen-air shear layers. These cases test the behavior of the new

19
(a) t = 0 (b) t = 10

Figure 6. Density contours are shown for a two dimensional isentropic vortex solved

using ESWENO3-4-3

Table 6. L2 error and convergence rates for the ESWENO3-4-3 and linear block

norm centered difference operators for the two dimensional inviscid vortex.

Linear Block Norm ESWENO3-4-3

Number of Cells L2 Error L2 Rate L2 Error L2 Rate

50 × 50 2.49E-05 - 5.32E-05 -

100 × 100 1.64E-06 3.92 2.25E-06 4.57

200 × 200 1.04E-07 3.98 1.08E-07 4.37

400 × 400 6.57E-09 3.99 6.62E-09 4.04

scheme where shocks, high gradients, and physical instabilities are present and interact. A successful scheme must provide enough dissipation to damp oscillations

near shocks and unresolved species and temperature gradients, but not be so dissipative that it damps the physical Kelvin-Helmholtz instability, which mixes the fuel

and oxidizer streams. In addition to shocks and high gradients, the cases presented

here include a discontinuous inlet condition and cannot be solved using standard

centered finite differences without the addition of filtering and limiters.

The hydrogen air mixing layer used for this study does not include any forcing

of the inlet conditions. The top stream has molar species fractions XH2 = 0.5 and

XN2 = 0.5, with a temperature of T = 300 K. The bottom stream is composed of

pure air, with inlet temperature of T = 300 K for non-reacting flow and T = 1500 K

for reacting flow. The Mach number of the top stream is M = 2.0, while the bottom

stream is M = 1.2. The full two dimensional compressible Navier-Stokes equations

were solved. ESWENO was used to calculate convective gradients, while a 4th

order, wide stencil SBP centered difference was used to calculate gradients of the

viscous terms. The species were assumed to be thermally perfect ideal gases, and

the open-source thermodynamics library Cantera [22] was used to calculate thermodynamic and transport properties. A hydrogen air chemical kinetic mechanism

20
Table 7. L∞ error and convergence rates for the ESWENO3-4-3 and linear block

norm centered difference operators for the two dimensional inviscid vortex.

Linear Block Norm ESWENO3-4-3

Number of Cells L∞ Error L∞ Rate L∞ Error L∞ Rate

50 × 50 2.46E-04 - 4.91E-04 -

100 × 100 1.68E-05 3.87 3.31E-05 3.89

200 × 200 1.07E-06 3.97 1.53E-06 4.44

400 × 400 6.71E-08 4.00 9.04E-08 4.08

Table 8. L2 error and convergence rates for the ESWENO2-4-2 and linear diagonal

norm centered difference operators for the two dimensional inviscid vortex.

Linear Diagonal Norm ESWENO2-4-2

Number of Cells L2 Error L2 Rate L2 Error L2 Rate

50 × 50 3.49E-05 - 6.48E-05 -

100 × 100 3.05E-06 3.52 4.43E-06 3.87

200 × 200 3.40E-07 3.17 3.64E-07 3.60

400 × 400 4.13E-08 3.04 4.16E-08 3.13

with 13 chemical species and 25 reactions was used.

The density contours for the supersonic shear layers are shown in Figure (7). In

both cases, the Kelvin-Helmholtz instability is clearly predicted, and the solution

is stable despite the presence of shocks and unresolved species and temperature

gradients. This suggests that the new ESWENO scheme exhibits the robustness

required for chemically reacting calculations without adding excessive numerical

dissipation.

6 Conclusions

A general strategy is presented in Refs. [5] and [6], for constructing Energy Stable

Weighted Essentially Non–Oscillatory (ESWENO) finite difference schemes on periodic

domains. ESWENO schemes up to eighth-order are developed and proved to be stable in the energy norm for both continuous and discontinuous solutions of systems

of linear hyperbolic equations. Herein, boundary closures are developed for the

fourth-order ESWENO scheme that maintain wherever possible the WENO stencil

biasing properties, while satisfying the summation-by-parts (SBP) operator convention, thereby ensuring stability in an L2 norm. A novel set of non-uniform flux

interpolation points is necessary near the boundaries to simultaneously achieve 1)

accuracy, 2) the SBP convention, and 3) WENO stencil biasing mechanics. The

novelty lies in the recognition that the discrete set of flux points x ̄M must be consistent with the stability norm P used in the SBP formulation, and in fact are derived

from the norm itself. Using the new flux-points x ̄M, second-order, and third-order

boundary closures are developed that achieve stability in a diagonal- and block-norm

sense, respectively.

21
Table 9. L∞ error and convergence rates for the ESWENO2-4-2 and linear diagonal

norm centered difference operators for the two dimensional inviscid vortex.

Linear Diagonal Norm ESWENO2-4-2

Number of Cells L∞ Error L∞ Rate L∞ Error L∞ Rate

50 × 50 3.17E-04 - 5.48E-04 -

100 × 100 2.42E-05 3.71 4.46E-05 3.62

200 × 200 2.28E-06 3.41 3.11E-06 3.84

400 × 400 2.42E-07 3.23 2.74E-07 3.51

Extensive numerical validation is presented to assess the efficacy of the new

boundary closures; the test problems included 1) unsteady 1-D linear wave equation (hyperbolic), 2) steady quasi one-dimensional nozzle flow solving the nonlinear Euler equations, 3) unsteady propagation of a 2-D Euler vortex, and 4) unsteady convection-diffusion-reaction of a supersonic hydrogen-air mixing layer. The

global accuracy of the second-order diagonal-norm operators and third-order closures block-norm operators, is shown to be “design-order”; i.e. three and four,

respectively, for all test cases. Furthermore, it is shown that the accuracy of the

new ESWENO operators is very close to that of the linear target operator. Finally,

an eigenvalue solver is used to show that the ESWENO discrete operator is stable

for smooth and discontinuous initial data.

Four appendices are included to provide detailed instructions on the implementation of the new ESWENO scheme. The appendices include sections describing the

1) smoothness indicators used by ESWENO, including near-wall stencil definitions,

2) a numerical recipe section that includes pseudo-code for the implementation of

ESWENO on nonlinear hyperbolic wave equations, (i.e. Euler), 3) the diagonalnorm interpolation and weights formulae, and the dissipation coefficients needed to

stabilize the baseline WENO approach, and 4) the block-norm interpolation and

weights formulae, and the dissipation coefficients needed to stabilize the baseline

WENO approach.

22
(a) Non-Reacting Flow

(b) Reacting Flow

Figure 7. Density contours for a non-reacting and reacting supersonic shear layer

utilizing fourth order ESWENO.

Appendix A

Smoothness Indicators

The smoothness indicators, β ̄

(r)

j+1/2

, are calculated at all flux points using the

expressions

β ̄L

j+1/2 = (fj − fj−1)

2

β ̄C

j+1/2 = (fj+1 − fj )

2

β ̄R

j+1/2 = (fj+2 − fj+1)

2

β ̄LL

j+1/2 = (fj−1 − fj−2)

2

β ̄RR

j+1/2 = (fj+3 − fj+2)

2

(A1)

(Note that the indicators β ̄LL and β ̄RR are only used near the boundaries, and then

only when a diagonal matrix norm is being used to prove stability.)

To guarantee that the downwind stencil weight does not exceed that of the

central or upwind weights, the downwind smoothness indicator is modified using

23
the expression

β ̄d

j+1/2 =

 

1

3

X

r

[β ̄

(r)

j+1/2

]

k

!1/k

, (A2)

where k is an even integer.

An additional stencil biasing parameter,  ̄τj+1/2, is needed for the ESWENO

scheme. Here,  ̄τj+1/2is a quadratic function of the highest undivided difference

available on the stencil for  ̄fj+1/2. From flux points xf2to xfN−2[see Fig. (1)],

τ ̄j+1/2 = (−fj−1 + 3fj − 3fj+1 + fj+2)

2

(A3)

At the point xf1,  ̄τj+1/2is biased toward the interior of the domain,

τ ̄j+1/2 = (−fj + 3fj+1 − 3fj+2 + fj+3)

2

, (A4)

and at xfN−1,  ̄τj+1/2is biased in a mirrored fashion.

The parameter  is a function of number of points in the discretization.

 = max(kf0k, kf

0

0k)x6=xd

(∆ξ)

3

, ∆ξ =

1

N

, (A5)

where kf0k and kf0k

0

represent a norm of the flux and the gradient of the flux from

the initial condition, where points near a discontinuity are excluded.A1

A1For applications where the solution changes drastically during simulation,  can be rescaled as

necessary.

24
Appendix B

Recipe

The recipe for calculating the gradient using the fourth-order ESWENO schemes

is summarized below.

1. Construct the convective flux f at solution points (xN ) according to the governing differential equation(s).

2. For systems of equations (e.g. Euler), construct the flux Jacobian matrices:

A =

∂F

∂U at the flux points (x ̄M) using Roe average variables formed from

nearest neighbor solution point data. Form the eigenvector decomposition

A = SΛS

−1

.

3. At the flux point  ̄xj+ 1

2

, identify all solution points xN that contribute to the

interpolants MI

(r)

N

, r = L, C, R. Use the eigenvector matrix S

−1

j+ 1

2

to rotate the

solution and flux at these points, into characteristic form C = S

−1U ; Fc(C) =

S

−1F(U).

4. Form the Lax-Friedrichs characteristic fluxes fc

± =

1

2

(fc ± λmaxc).

5. Perform interpolations on each candidate stencil,  ̄f

(r) = MI

(r)

N

fc.

6. Calculate the stencil biasing parameters:

(a) Calculate β ̄rfor each flux point according to Eq. (A1), except at the

endpoints.

(b) Calculate  ̄τj+1/2 according to Eq. (A3) and Eq. (A4).

7. Calculate and normalize the weights using the stencil biasing parameters and

the target weights in Eq. (14).

8. Calculate Λi from the weights using Eqs. (C3), (C4), and (C5) for the diagonalnorm scheme or Eqs. (D5) and (D6) for the block-norm.

9. Modify the diagonal of Λi to be smoothly positive according to Eq. (20).

10. Calculate the WENO flux from the weights and candidate interpolations using

Eq. (13).

11. Calculate the Energy Stable flux, ψ ̄, from Eq. (23).

12. Reconstruct the fluxes in characteristic space,  ̄fj+1/2 =  ̄f

+

j+1/2 +  ̄f

−

j+1/2

and

ψ ̄

j+1/2 = ψ ̄+

j+1/2 + ψ ̄−

j+1/2

, and transform back to physical space.

13. Calculate the gradient using the inverse of the P-norm, as in Eq. (24).

25
Remark. It is important to note that all equations presented are for the forward

propagating waves,  ̄f

+

j+1/2

. The equations for interpolation of backward propagating

waves (  ̄f

−

j+1/2

) are found simply by flipping r = (LL, L, C, R, RR) about C. For the

Λi terms, after flipping r about C, the diagonal matrix is multiplied by −1. Eq.

(23) is also multiplied by −1 to calculate ψ ̄−

j+1/2

.

26
Appendix C

Fourth-order Diagonal-Norm Operator: D2-4-2

C.1 Differentiation Matrix

The target SBP differentiation operator D2-4-2 can be written in the form

D =

1

∆x





−

24

17

59

34 −

4

17 −

3

34 0 0 0 0 0 0 0

−

1

2

0

1

2

0 0 0 0 0 0 0 0

4

43 −

59

86 0

59

86 −

4

43 0 0 0 0 0 0

3

98 0 −

59

98 0

32

49 −

4

49 0 0 0 0 0

0 0 1

12 −

2

3

0

2

3 −

1

12 0 0 0 0

0 0 0 .

.

.

.

.

.

.

.

.

.

.

.

.

.

. 0 0 0

0 0 0 0 1

12 −

2

3

0

2

3 −

1

12 0 0

0 0 0 0 0 4

49 −

32

49 0

59

98 0 −

3

98

0 0 0 0 0 0 4

43 −

59

86 0

59

86 −

4

43

0 0 0 0 0 0 0 0 −

1

2

0

1

2

0 0 0 0 0 0 0 3

34

4

17 −

59

34

24

17





where the diagonal mass matrix P is

P = ∆xDiag 

17

48

59

48

43

48

49

48 1 · · · 1

49

48

43

48

59

48

17

48 

The skew-symmetric matrix “Q” follows immediately from P D.

27
C.2 Interpolation Operators

The WENO/ESWENO2-4-2 requires the inclusion of two extra stencils near the

boundary. The interpolation coefficients are presented below.

ILL =





0 0

0 0

0 0

−

71

48

119

48

0 0

.

.

.

.

.

.

0 0

−

73

48

121

48

−

17

12

29

12

−

79

48

127

48

0 0





IL =





0 0

0 0

−

7

12

19

12

−

23

48

71

48

−

1

2

3

2

.

.

.

.

.

.

−

1

2

3

2

−

25

48

73

48

−

5

12

17

12

−

31

48

79

48

0 0





IC =





1 0

31

48

17

48

5

12

7

12

25

48

23

48

1

2

1

2

.

.

.

.

.

.

1

2

1

2

23

48

25

48

7

12

5

12

17

48

31

48

0 1





IR =





0 0

79

48 −

31

48

17

12 −

5

12

73

48 −

25

48

3

2 −

1

2

.

.

.

.

.

.

3

2 −

3

2

71

48 −

23

48

19

12 −

7

12

0 0

0 0





IRR =





0 0

127

48 −

79

48

29

12 −

17

12

121

48 −

73

48

0 0

.

.

.

.

.

.

0 0

119

48 −

71

48

0 0

0 0

0 0





(C1)

C.3 Target Weights

The target weights for the five interpolation stencils are:

d242 =





0 0 1 0 0

0 0 24

31

1013

4898

3

158

0

11

56

51

70

3

40 0

3

142

357

3266

408

575

4

25 0

0

1

6

2

3

1

6

0

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

0

1

6

2

3

1

6

0

0

4

25

408

575

357

3266

3

142

0

3

40

51

70

11

56 0

3

158

1013

4898

24

31 0 0

0 0 1 0 0





(C2)

28
C.4 Energy Stable Terms

The nonzero terms of the diagonal matrices, Λi, for the ESWENO2-4-2 are

Diag (Λ1) =





0

1

96 

14 − 28w

(L)

2 − 71w

(LL)

3 − 31w

(R)

1 − 31w

(RR)

1



1

96 

−8 + 28w

(L)

2 − 23w

(L)

3 + 48w

(LL)

3 + 31w

(R)

1 − 20w

(R)

2 − 48w

(RR)

1



1

96 

2 + 23w

(L)

3 − 24w

(L)

4 + 23w

(LL)

3 + 20w

(R)

2 − 25w

(R)

3 + 79w

(RR)

1



1

96 

24w

(L)

4 − 24w

(L)

5 + 25w

(R)

3 − 24w

(R)

4



.

.

.

1

4



w

(L)

i − w

(L)

i+1 + w

(R)

i−1 − w

(R)

i



.

.

.

1

96 

24w

(L)

N−4 − 25w

(L)

N−3 + 24w

(R)

N−5 − 24w

(R)

N−4



1

96 

−2 + 25w

(L)

N−3 − 79w

(LL)

N−1 − 20w

(L)

N−2 + 24w

(R)

N−4 − 23w

(R)

N−3 − 23w

(RR)

N−3



1

96 

8 − 31w

(L)

N−1 + 48w

(LL)

N−1 + 20w

(L)

N−2 + 23w

(R)

N−3 − 48w

(RR)

N−3 − 28w

(R)

N−2



1

96 

−14 + 31w

(L)

N−1 + 31w

(LL)

N−1 + 71w

(RR)

N−3 + 28w

(R)

N−2



0





(C3)

Diag (Λ2) =





0

1

96 

28w

(L)

2 + 142w

(LL)

3 − 31w

(R)

1 − 110w

(RR)

1



1

96 

23w

(L)

3 + 94w

(LL)

3 − 20w

(R)

2 − 158w

(RR)

1



1

96 

24w

(L)

4 − 25w

(R)

3



.

.

.

1

4



w

(L)

i − w

(R)

i−1



.

.

.

1

96 

25w

(L)

N−3 − 24w

(R)

N−4



1

96 

158w

(LL)

N−1 + 20w

(L)

N−2 − 23w

(R)

N−3 − 94w

(RR)

N−3



1

96 

31w

(L)

N−1 + 110w

(LL)

N−1 − 142w

(RR)

N−3 − 28w

(R)

N−2



0





(C4)

29
Diag (Λ3) =





0

0

1

96 

−71w

(LL)

3 + 79w

(RR)

1



.

.

.

0

.

.

.

1

96 

−79w

(LL)

N−1 + 71w

(RR)

N−3



0

0





(C5)

The matrix Diag (Λ0) is identically zero.

30
Appendix D

Fourth-order Block-Norm Operator: D3-4-3

D.1 Differentiation Matrix

The target SBP differentiation operator D3-4-3 satisfies the expression D3-4-3 =

P

−1

3-4-3Q3-4-3 [see equation (10)]. The inverse P-norm for the ESWENO3-4-3 scheme

is given by

P

−1

3-4-3 =

1

∆x





P

−1

0

0 0

0 I 0

0 0 P

−1

0

P T



 ; P

−1

0 =





71043003

15979414 −

7964853

15979414

4820919

15979414

2199519

15979414

−

7964853

15979414

14093655

15979414

331299

15979414 −

1445721

15979414

4820919

15979414

331299

15979414

18932367

15979414

253203

15979414

2199519

15979414 −

1445721

15979414

253203

15979414

15938451

15979414





(D1)

where P T denotes the per-symmetric transpose. The Q3-4-3 matrix is given by

Q3-4-3 =





−

1

2

101

144 −

29

144 0 0 0 0 0 0 0 0

−

101

144 0

17

24 −

1

144 0 0 0 0 0 0 0

29

144 −

17

24 0

85

144 −

1

12 0 0 0 0 0 0

0

1

144 −

85

144 0

2

3 −

1

12 0 0 0 0 0

0 0 1

12 −

2

3

0

2

3 −

1

12 0 0 0 0

0 0 0 .

.

.

.

.

.

.

.

.

.

.

.

.

.

. 0 0 0

0 0 0 0 1

12 −

2

3

0

2

3 −

1

12 0 0

0 0 0 0 0 1

12 −

2

3

0

85

144 −

1

144 0

0 0 0 0 0 0 1

12 −

85

144 0

17

24 −

29

144

0 0 0 0 0 0 0 1

144 −

17

24 0

101

144

0 0 0 0 0 0 0 0 29

144 −

101

144

1

2





(D2)

31
D.2 Interpolation Operators

The interpolation coefficients for the ESWENO3-4-3 are

IL =





0 0

0 0

−

25

36

61

36

−

61

144

205

144

−

1

2

3

2

.

.

.

.

.

.

−

1

2

3

2

−

83

144

227

144

−

11

36

47

36

−

101

144

245

144

0 0





IC =





1 0

101

144

43

144

11

36

25

36

83

144

61

144

1

2

1

2

.

.

.

.

.

.

1

2

1

2

61

144

83

144

25

36

11

36

43

144

101

144

0 1





IR =





0 0

245

144 −

101

144

47

36 −

11

36

227

144 −

83

144

3

2 −

1

2

.

.

.

.

.

.

3

2 −

1

2

205

144 −

61

144

61

36 −

25

36

0 0

0 0





(D3)

D.3 Target Weights

The target weights for each stencil in ESWENO3-4-3 are given below for each flux

point.

d343 =





0 1 0

0

72

101

29

101

29

100

189

275

1

44

1

61

4248

5063

12

83

1

6

2

3

1

6

.

.

.

.

.

.

.

.

.

1

6

2

3

1

6

12

83

4248

5063

1

61

1

44

189

275

29

100

29

101

72

101 0

0 1 0





(D4)

32
D.4 Energy Stable Terms

The nonzero terms of the diagonal matrices, Λi, for the ESWENO3-4-3 are

Diag (Λ1) =





0

1

288 

58 − 100w

(L)

2 − 101w

(R)

1



1

288 

−56 + 100w

(L)

2 − 61w

(L)

3 + 101w

(R)

1 − 44w

(R)

2



1

288 

22 + 61w

(L)

3 − 72w

(L)

4 + 44w

(R)

2 − 83w

(R)

3



1

288 

72w

(L)

4 − 72w

(L)

5 + 83w

(R)

3 − 72w

(R)

4



.

.

.

1

4



w

(L)

i − w

(L)

i+1 + w

(R)

i−1 − w

(R)

i



.

.

.

1

288 

72w

(L)

N−4 − 83w

(L)

N−3 + 72w

(R)

N−5 − 72w

(R)

N−4



1

288 

−22 + 83w

(L)

N−3 − 44w

(L)

N−2 + 72w

(R)

N−4 − 61w

(R)

N−3



1

288 

56 − 101w

(L)

N−1 + 44w

(L)

N−2 + 61w

(R)

N−3 − 100w

(R)

N−2



1

288 

−58 + 101w

(L)

N−1 + 100w

(R)

N−2



0





(D5)

Diag (Λ2) =





0

1

288 

100w

(L)

2 − 101w

(R)

1



1

288 

61w

(L)

3 − 44w

(R)

2



1

288 

72w

(L)

4 − 83w

(R)

3



.

.

.

1

4



w

(L)

i − w

(R)

i−1



.

.

.

1

288 

83w

(L)

N−3 − 72w

(R)

N−4



1

288 

44w

(L)

N−2 − 61w

(R)

N−3



1

288 

101w

(L)

N−1 − 100w

(R)

N−2



0





(D6)

33
Appendix E

References

1. Erlebacher, G.; Hussaini, M.; and Shu, C.-W.: Interaction of a shock with a

longitudinal vortex. Journal of Fluid Mechanics, vol. 337, 1997, pp. 129–153.

2. Zhang, S.; Jiang, S.; Zhang, Y.; and Shu, C.: The mechanism of sound generation in the interaction between a shock wave and two counter-rotating vortices.

Physics of Fluids, vol. 21(7), 2009.

3. Shen, Y.; Wang, B.; and Zha, G.: Implicit WENO Scheme and High Order

Viscous Formulas for Compressible Flows. 25th AIAA Applied Aerodynamics

Conference, no. 2007-4431, 2007.

4. Cockburn, B.; Johnson, C.; Shu, C.-W.; and Tadmor, E.: Advanced Numerical

Approximation of Nonlinear Hyperbolic Equations. Springer-Verlag, 1998.

5. Yamaleev, N. K.; and Carpenter, M. H.: Third-order Energy stable WENO

Scheme. Journal of Computational Physics, vol. 228, 2009, pp. 3025–3047.

6. Yamaleev, N. K.; and Carpenter, M. H.: A systematic methodology for constructing high-order energy stable WENO schemes. Journal of Computational

Physics, vol. 228, 2009, pp. 4248–4272.

7. Jiang, G.; and Shu, C.-W.: Efficient implementation of weighted ENO schemes.

Journal of Computational Physics, vol. 126, 1996, pp. 202–228.

8. Martin, M.; Taylor, E.; Wu, M.; and Weris, V.: A bandwidth-optimized WENO

scheme for the effective direct numerical simulation of compressible turbulence.

Journal of Computational Physics, vol. 220, 2006, pp. 270–289.

9. Kreiss, H.; and Scherer, G.: Mathematical Aspects of Finite Elements in Partial

Differential Equations, Academic Press, New York, Finite element and finite

difference methods for hyperbolic partial differential equations. 1974.

10. Carpenter, M. H.; Nordstr ̈om, J.; and Gottlieb, D.: A stable and conservative interface treatment of arbitrary spatial accuracy. Journal of Computational

Physics, vol. 148, 1999, pp. 341–365.

11. Strand, B.: Summation by parts for finite difference approximations for d/dx.

Journal of Computational Physics, vol. 110(1), 1994, pp. 47–67.

12. Nordstr ̈om, J.: Conservative finite difference formulations, variable coefficients,

energy estimates and artificial dissipation. Journal of Scientific Computing,

vol. 29, 2006, pp. 375–404.

34
13. Kitson, A.; McLachlan, R. I.; and Robidoux, N.: Skew-adjoint finite difference

methods on nonuniform grids. New Zealand Journal of Mathematics, vol. 32,

no. 2, 2003, pp. 139–159.

14. McLachlan, R. I.; and Robidoux, N.: Antisymmetry, pseudospectral methods,

and conservative PDEs, World Scientific, EQUADIFF 99, B. Fiedler, K. Groger,

and J. Sprekels, eds. 2000, pp. 994–999.

15. Gustafsson, B.: The convergence rate for difference approximations to mixed

initial boundary value problems. Mathematics of Computation, vol. 29, 1975,

pp. 396–406.

16. Sv ̈ard, M.; and Nordstr ̈om, J.: On the order of accuracy for difference approximations of initial-boundary value problems. Journal of Computational Physics,

vol. 218, 2006, pp. 333–352.

17. http://www.wolfram.com/products/mathematica/index.html, 2009.

18. Carpenter, M. H.; and Kennedy, C. A.: Fourth-Order 2N-Storage Runge-Kutta

Schemes. TM 109112, NASA, 1994.

19. Kim, J. W.: Optimized boundary compact finite difference schemes for computational aeroacoustics. Journal of Computational Physics, vol. 225, 2007,

pp. 995–1019.

20. Carpenter, M. H.; Gottlieb, D.; and Abarbanel, S.: Time-Stable Boundary Conditions for Finite-Difference Schemes Solving Hyperbolic Systems: Methodology

and Application to High-Order Compact Schemes. Journal of Computational

Physics, vol. 111, 1994, pp. 220–236.

21. Sv ̈ard, M.; Carpenter, M. H.; and Nordstr ̈om, J.: A stable high-order finite difference scheme for the compressible Navier-Stokes equations, far-field boundary

conditions. Journal of Computational Physics, vol. 225, 2007, pp. 1020–1038.

22. http://www.cantera.org, 2009.

35
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704–0188

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection

of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports

(0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be

subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18

1. REPORT DATE (DD-MM-YYYY)

01-12-2009

2. REPORT TYPE

Technical Memorandum

3. DATES COVERED (From - To)

10/2009–12/2009

4. TITLE AND SUBTITLE

Boundary Closures For Fourth-order Energy Stable Weighted Essentially

Non-Oscillatory Finite Difference Schemes

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

599489.02.07.07.03.13.01

6. AUTHOR(S)

Travis C. Fisher, Mark H. Carpenter, Nail K. Yamaleev, Steven H. Frankel

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, Virginia 23681-2199

8. PERFORMING ORGANIZATION

REPORT NUMBER

L–19792

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

10. SPONSOR/MONITOR’S ACRONYM(S)

NASA

11. SPONSOR/MONITOR’S REPORT

NUMBER(S)

NASA/TM–2009–216166

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Category 64

Availability: NASA CASI (443) 757-5802

13. SUPPLEMENTARY NOTES

An electronic version can be found at http://ntrs.nasa.gov.

14. ABSTRACT

A general strategy exists for constructing Energy Stable Weighted Essentially Non–Oscillatory (ESWENO) finite difference

schemes up to eighth-order on periodic domains. These ESWENO schemes satisfy an energy norm stability proof for both

continuous and discontinuous solutions of systems of linear hyperbolic equations. Herein, boundary closures are developed for

the fourth-order ESWENO scheme that maintain wherever possible the WENO stencil biasing properties, while satisfying the

summation-by-parts (SBP) operator convention, thereby ensuring stability in an L2 norm. Second-order, and third-order

boundary closures are developed that achieve stability in diagonal and block norms, respectively. The global accuracy for the

second-order closures is three, and for the third-order closures is four. A novel set of non-uniform flux interpolation points is

necessary near the boundaries to simultaneously achieve 1) accuracy, 2) the SBP convention, and 3) WENO stencil biasing

mechanics.

15. SUBJECT TERMS

high-order finite difference methods, weighted essentially nonoscillatory,schemes, energy estimate, numerical stability,

artificial dissipation

16. SECURITY CLASSIFICATION OF:

a. REPORT

U

b. ABSTRACT

U

c. THIS PAGE

U

17. LIMITATION OF

ABSTRACT

UU

18. NUMBER

OF

PAGES

40

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802