Table Of Content47thAIAAAeroscpaceSciencesMeeting,Orlando,FL,USA
High-Order Energy Stable WENO Schemes
Nail K. Yamaleev(cid:3)
North Carolina A&T State University, Greensboro, 27411, USA
Mark H. Carpentery
NASA Langley Research Center, Hampton, VA 23681, USA
A third-order Energy Stable Weighted Essentially Non{Oscillatory (ESWENO) (cid:12)nite
di(cid:11)erence scheme developed by Yamaleev and Carpenter was proven to be stable in the
energy norm for both continuous and discontinuous solutions of systems of linear hyper-
bolic equations. Herein, a systematic approach is presented that enables \energy stable"
modi(cid:12)cations for existing WENO schemes of any order. The technique is demonstrated
by developing a one-parameter family of (cid:12)fth-order upwind-biased ESWENO schemes;
ESWENO schemes up to eighth order are presented in the appendix. New weight func-
tions are also developed that provide (1) formal consistency, (2) much faster convergence
for smooth solutions with an arbitrary number of vanishing derivatives, and (3) improved
resolution near strong discontinuities.
Nomenclature
a = constant wave speed
d(r) = target value for w arising from stencil S
(r)
D = matrix de(cid:12)ning the discrete derivative operator
D = arti(cid:12)cial dissipation matrix operator
ad
D = skew-symmetric portion of derivative matrix
skew
D = symmetric skew-symmetric portion of derivative matrix
sym
D(cid:22) = energy-stable derivative operator
D = discrete undivided di(cid:11)erence operator
1
f = continuous (cid:13)ux function (linear or nonlinear)
f^ = discrete (cid:13)ux function built at position j+1=2
j+1=2
f(cid:22) = projection of continuous (cid:13)ux onto grid (linear or nonlinear)
f0 = derivative of the (cid:13)ux function
h(x) = numerical (cid:13)ux function implicitly de(cid:12)ning f(x)
j = grid index
P = symmetric positive de(cid:12)nite matrix de(cid:12)ning discrete norm
Q = skew-symmetric matrix de(cid:12)ning dispersive portion of derivative
R = symmetric matrix de(cid:12)ning dissipative portion of derivative
S , S = Stencil shifted one cell to the left, right
L R
u = solution
(cid:3)AssociateProfessor,DepartmentofMathematics,E-mail: nkyamale@ncat.edu, MemberAIAA.
ySeniorResearchScientist,ComputationalAerosciencesBranch,E-mail: mark.h.carpenter@nasa.gov,MemberAIAA.
1of30
AmericanInstitute ofAeronautics andAstronautics
w(r) = Weight function for stencil S
(r)
w(r) = Weight function for stencil S , everywhere in the computational domain
: (r)
x = spatial coordinate
(cid:11) = solution dependent component of weight function arising from stencil S
r (r)
(cid:12)(r) = classical smoothness indicator for stencil r
(cid:14) = nonzero parameter used to prevent division by zero (ESWENO)
(cid:1)x = grid spacing
(cid:15) = nonzero parameter used to prevent division by zero (WENO)
(cid:21)5 = diagonal matrix resulting from elemental decomposition of 5th-order R
0
(cid:21) = local component resulting from elemental decomposition of R
i0
(cid:21)(cid:22) = strictly positive local dissipation term
i0
(cid:27) = (cid:12)fth-order solution dependent component of weight function
5
(cid:27)~ = (cid:12)fth-order solution dependent function proposed by Borgeset al.8
5
’ = parameter in one parameter family of 5th-order schemes
I. Introduction
Anew,third-orderweightedessentiallynonoscillatoryscheme(calledEnergy-StableWENO)hasrecently
been proposed and developed by the authors of the present paper.1 In reference,1 we prove that the third-
orderESWENOschemeisenergystable,thatis,stableinanL energynorm,forsystemsoflinearhyperbolic
2
equationswithbothcontinuousanddiscontinuoussolutions. Stabilityisexplicitlyachieved(byconstruction)
by requiring that the ESWENO scheme satis(cid:12)es a nonlinear summation-by-parts (SBP) condition at each
instant in time. Thus, L strict stability is attained without the need for a total variation bounded (TVB)
2
(cid:13)uxreconstructionoralarge-time-stepconstraint,23 and.4 Herein, wegeneralizeandextendthe third-order
ESWENO methodology1 to an arbitrary order of accuracy. Similar to the third-order ESWENO scheme,
the new families of higher order (up to eighth order) ESWENO schemes are provably stable in the energy
norm and retain the underlying WENO characteristics of the background schemes. Numerical experiments
demonstrate that the new family of ESWENO schemes provides the design order of accuracy for smooth
problems and delivers stable essentially nonoscillatorysolutions for problems with strong discontinuities.
Another issue that is addressed in this paper is the consistency of the new class of ESWENO schemes.
The consistency of any WENO-type scheme fully depends on a proper choice of the weight functions. On
one hand, for smooth solutions the weights should provide a rapid convergence of the WENO scheme to
the corresponding underlying linear scheme. On the other hand, the weights should e(cid:11)ectively bias the
stencil away from strong discontinuities. The high-order upwind-biased WENO schemes with conventional
smoothness indicators that are presented in reference2 are too dissipative for solving problems with a large
amount of structure in the smooth part of the solution, such as direct numerical simulations of turbulence,
or aeroacoustics,5.6 Furthermore, as has been shown in references7 and,8 the classical weight functions of
the (cid:12)fth-order WENO scheme fail to provide the design order of convergence near smooth extrema, where
the (cid:12)rst derivative of the solution becomes equal to zero. New approaches are proposed in references7 and8
to improve the error convergence near the critical points. Although these new weight functions recover the
(cid:12)fth order of convergenceof the WENO scheme near smooth extrema, the problem persists if the (cid:12)rst- and
second-orderderivatives vanish simultaneously.8 An attempt to resolve this loss of accuracy is presented in
reference.8 This proposed resolutionprovidesonly a partialremedy for the problem; the same degeneration
inthe orderofconvergenceoccursifatleastthe (cid:12)rstthreederivativesbecomeequaltozero. Tofully resolve
thisproblem,weproposenewweightstoprovidefastererrorconvergencethanthosepresentedinreference,8
andimposesomeconstraintsontheweightparameterstoguaranteethattheWENOandESWENOschemes
aredesign-orderaccurateforsu(cid:14)cientlysmoothsolutionswithanarbitrarynumberofvanishingderivatives.
This paper is organized as follows. In section 2, energy estimates for the continuous and corresponding
discretewaveequationsarepresented. Insection3,wepresentaone-parameterfamilyof (cid:12)fth-orderWENO
2of30
AmericanInstitute ofAeronautics andAstronautics
schemes; one value of the parameter yields a central scheme that converges with sixth-order accuracy.
In section 4, we present a systematic methodology for constructing ESWENO schemes of any order and
demonstrate the methodology by transforming the family of WENO schemes presented in section 3, into a
familyof(cid:12)fth-orderESWENOschemes. Insection5,weanalyzetheconsistencyofthenewclassofESWENO
schemesandwederivesu(cid:14)cientconditionsfortheweightsfunctionsthatensurethattheESWENOschemes
are design-order accurate regardless of the number of vanishing derivatives in the solution. The tuning
parameters in the weight functions are also optimized. In section 6, we present numerical experiments that
corroborateour theoretical results. We summarize and draw conclusions in section 7.
II. Energy Estimates
Consider a linear, scalar wave equation
@u + @f =0; f =au; t(cid:21)0; 0(cid:20)x(cid:20)1;
@t @x (1)
u(0;x)=u (x)
0
where a is a constant and u (x) is a bounded piecewise continuous function. Without loss of generality,
0
assumethata(cid:21)0,andfurther assumethat the problem isperiodic onthe interval0(cid:20)x(cid:20)1. Applying the
energy method to equation (1) leads to
d
kuk2 =0 (2)
dt L2
where k(cid:1)k is the continuous L norm. Thus, the continuous problem de(cid:12)ned in equation (1) is neutrally
L2 2
stable.
Wenowdevelopusingmimetic techniques(see1 or9) aclassof discretespatialoperatorsthatisneutrally
stable or dissipative. The continuous target operator used for this development is the following singular
perturbed wave equation:
N
@u + @f = ((cid:0)1)n(cid:0)1 @n (cid:22)@nu ; f =au; t(cid:21)0; 0(cid:20)x(cid:20)1;
@t @x @xn @xn (3)
n=1
u(0;x)=u P(x) ; (cid:0) (cid:1)
0
where (cid:22) = (cid:22)(u) is a non-negative C1 function of u. As before, we assume that equation (3) is subject to
periodic boundary conditions. Our goal is to match each spatial term in equation (3) with an equivalent
discrete term that maintains neutral stability (or dissipates) of the discrete energy norm.
Webeginbyshowingthatthetermsontheright-handsideofequation(3)aredissipativetherebyensuring
stability. Multiplying equation (3) by u and integrating it overthe entire domain yields
1
1 d 1 1 N @n @nu
kuk2 + au2 = ((cid:0)1)n(cid:0)1u (cid:22) dx : (4)
2dt L2 2 @xn @xn
(cid:12)0 n=1Z (cid:18) (cid:19)
(cid:12) X0
(cid:12)
(cid:12)
Integratingeachtermontheright-handsidebypartsandaccountingforperiodicboundaryconditionsyields
the following energy estimate:
1
d N @nu 2
kuk2 =(cid:0)2 (cid:22)(u) dx(cid:20)0 : (5)
dt L2 @xn
n=1Z (cid:18) (cid:19)
X0
All the perturbation terms included in equation (3) provide dissipation of energy.
Turning now to the discrete case, we de(cid:12)ne a uniform grid x = j(cid:1)x, j = 0;J, with (cid:1)x = 1=J.
j
On this grid, we de(cid:12)ne a (cid:13)ux (cid:22)f = au(cid:22) and its derivative (cid:22)f = au(cid:22) , where u(cid:22) = [u(x ;t);:::;u(x ;t)]T
x x 0 J
and u(cid:22) = [u (x ;t);:::; u (x ;t)]T are projections of the continuous solution and its derivative onto the
x x 0 x J
3of30
AmericanInstitute ofAeronautics andAstronautics
computationalgrid. Next,wede(cid:12)neapth-orderapproximationforthe(cid:12)rst-orderderivativeterminequation
(1) as
@(cid:22)f
=D(cid:22)f +O((cid:1)xp) : (6)
@x
Placing a mild restriction on the generality of the derivative operator (see1 or9), the matrix D can be
expressed in the following form:
D=P(cid:0)1[Q+R] ; Q+QT =0
R=RT ; vTRv(cid:21)0 (7)
P =PT ; vTPv>0
for any real vector v 6= 0. By choosing the matrix R as a discrete analog of the dissipation operator in
equation (3), we have
N
R= D n(cid:3)[D n]T ; (8)
1 1
n=1
X
where (cid:3) is a diagonal positive semide(cid:12)nite matrix and D is the di(cid:11)erence matrix
1
.. 0
.
D1 =2 (cid:0)1 1 3 :
0 ..
6 . 7
6 7
4 5
By using the SBP operators[eqs. (6){(8)], the semi-discrete counterpart of equation (1) becomes
N
@u
+P(cid:0)1Qf =(cid:0) P(cid:0)1D n(cid:3)[D n]Tf; (9)
1 1
@t
n=1
X
where f =au, u=[u (t);u (t);:::;u (t);]T is the discrete approximation of the solution u of equation (1),
0 1 J
and Q and (cid:3) arenonlinearmatrices (i.e., Q=Q(u) and (cid:3)=(cid:3)(u). To showthat the above (cid:12)nite-di(cid:11)erence
scheme is stable, the energy method is used. Multiplying equation (9) with uTP yields
N
1 d kuk2 +auTQu=(cid:0)a [D n]Tu T (cid:3)[D n]Tu; (10)
2dt P 1 1
n=1
X(cid:0) (cid:1)
where k(cid:1)k is the P norm (i.e., kuk2 =uTPu). Adding equation (10) to its transpose yields
P P
N
d kuk2 +auT Q+QT u=(cid:0)2a [D n]Tu T (cid:3)[D n]Tu : (11)
dt P 1 1
n=1
(cid:0) (cid:1) X(cid:0) (cid:1)
If we account for periodic boundary conditions and the skew-symmetry of Q, then the second term on the
left-hand side vanishes, and the energy estimate becomes
N
d kuk2 =(cid:0)2a [D n]Tu T (cid:3)[D n]Tu(cid:20)0 : (12)
dt P 1 1
n=1
X(cid:0) (cid:1)
The right-hand side of equation (12) is nonpositive because the diagonal matrix (cid:3) is positive semide(cid:12)nite
[vT(cid:3)v(cid:21)0 for all real v of length (J +1)] and a(cid:21)0; thus the stability of the (cid:12)nite-di(cid:11)erence scheme given
by equation (9) is assured. This result can be summarized in the following theorem:
Theorem 1. The approximation [eq. (9)] of the problem [eq. (1)] is stable if equations [(6){(8)] hold.
4of30
AmericanInstitute ofAeronautics andAstronautics
Remark 1. Despite the fact that the initial boundary value problem [eq. (1)] is linear, the (cid:12)nite-di(cid:11)erence
scheme [eq. (9)] constructed for approximation of equation (1) is nonlinear, because the matrices Q and (cid:3)
(and in principle P) are assumed to depend on the discrete solution u.
Remark 2. The onlyconstraintsthat areimposed onthe matrixQand the diagonalmatrix(cid:3)arethe skew-
symmetry of the former and positive semide(cid:12)niteness of the later. No other assumptions have been made
about a speci(cid:12)c form of the matrices Q and (cid:3) to guarantee the stability of the (cid:12)nite-di(cid:11)erence scheme [eq.
(9)].
Remark 3. The discrete operators that are de(cid:12)ned by equations [(6){(8)] are similar in form to those that
are used for conventional SBP operators(See refs.9,10). What is new, however,is the fact that the matrices
Q and (cid:3) depend on u.
f
j+1/2
S6
x x x x x x
j-2 j-1 j j+1 j+2 j+3
S x
LL j+1/2
S
L
S
R
S
RR
Figure1. Extendedsix-pointstencilS6,andcorrespondingcandidatestencilsSLL;SL;SR;andSRR forone-parameter
family of(cid:12)fth-order WENO schemes.
Next, a new one-parameter family of (cid:12)fth-order WENO schemes is developed, and then is used as the
starting point in the development of a family of \energy-stable"WENO schemes. Great careis exercised in
thedevelopingtheWENOschemestoensurethatdesign-orderaccuracyisachievedinthevicinityofsmooth
extrema.
III. Fifth- and Sixth-order WENO Schemes
Anyconventionalhigh-orderWENO(cid:12)nite-di(cid:11)erenceschemeforthescalarone-dimensionalwaveequation
(1) can be written in the following semidiscrete form:
duj + f^j+12 (cid:0)f^j(cid:0)21 =0: (13)
dt (cid:1)x
For the (cid:12)fth-order WENO scheme that is presented in reference,2 the numerical (cid:13)ux f^ is computed
j+1
2
as a convex combination of three third-order (cid:13)uxes de(cid:12)ned on the following three-point stencils: S =
LL
fx ;x ;x g, S = fx ;x ;x g, and S = fx ;x ;x g. (See (cid:12)gure 1.) Note that this set of
j(cid:0)2 j(cid:0)1 j L j(cid:0)1 j j+1 R j j+1 j+2
stencils is not symmetric with respect to the (j + 1) point; thus, the (cid:12)fth-order WENO scheme is biased
2
in the upwind direction. A central WENO scheme can be constructed from the conventional (cid:12)fth-order
WENO scheme by including an additional downwind candidate stencil S = fx ;x ;x g, so that
RR j+1 j+2 j+3
5of30
AmericanInstitute ofAeronautics andAstronautics
the collection of allfour stencils issymmetric with respect to point (j+ 1). a The WENO (cid:13)ux, constructed
2
in this manner, is given by
f^ =wLL fLL +wL fL +wR fR +wRR fRR ; (14)
j+1 j+1=2 j+1=2 j+1=2 j+1=2 j+1=2 j+1=2 j+1=2 j+1=2
2
where f(r) , r=fLL;L;R;RRgare third-order (cid:13)uxes de(cid:12)ned on these four stencils:
j+1=2
f(u )
j(cid:0)2
fLL(u ) 2 (cid:0)7 11 0 f(u )
j+1=2 0 j(cid:0)1 1
fL(u ) (cid:0)1 5 2 f(u )
0 j+1=2 1 = 10 1B j C (15)
fR(u ) 6 2 5 (cid:0)1 B f(u ) C
BBB fRR(ujj++11==22) CCC BBB 0 11 (cid:0)7 2 CCCBBBB f(ujj++21) CCCC
@ A @ AB f(uj+3) C
B C
@ A
andwLL,wL,wR,andwRR areweightfunctionsthatareassignedtothefourstencilsS ,S ,S ,andS ,
LL L R RR
respectively. The terms wLL, wL, wR, and wRR in equation (14) are nonlinear weight functions. These
have both preferred values that are derived from an underlying linear scheme as well as solution-dependent
components. The preferred values are given by
dLL = 1 (cid:0)’ ; dL = 6 (cid:0)3’ ; dR = 3 +3’ ; dRR = ’ ; (16)
10 10 10
where ’ is a parameter. The convergence rate of the scheme [eqs. (13){(16)] with the preferred values
w(r) =d(r) andr =fLL;L;R;RRgisgreaterthanorequalto5forallvaluesoftheparameter’inequation
(16). For the speci(cid:12)c value ’ = 1 , the (cid:12)fth-order term vanishes and the convergence rate increases to 6.
c 20
The classical (cid:12)fth-order upwind-biased WENO scheme of Jiang and Shu is obtained for ’=0.
The weight functions wLL, wL, wR, and wRR needed in equation (14) are given by
w(r) = (cid:11)r ;
j+21 Pl (cid:11)l (17)
where
(cid:28)
(cid:11) =d(r) 1+ p ; r=fLL;L;R;RRg; (18)
r (cid:15)+(cid:12)(r)
(cid:18) (cid:19)
and (cid:15) (cid:20) O((cid:1)x2). The functions (cid:12)(r) are the classical smoothness indicators. [See equation (59) for the
general expression.] For (cid:12)fth-order schemes, they are given by
(cid:12)LL = 13(u (cid:0)2u +u )2 + 1(u (cid:0)4u +3u )2
12 j(cid:0)2 j(cid:0)1 j 4 j(cid:0)2 j(cid:0)1 j
(cid:12)L = 13(u (cid:0)2u +u )2 + 1(u (cid:0)u )2
12 j(cid:0)1 j j+1 4 j(cid:0)1 j+1 (19)
(cid:12)R = 13(u (cid:0)2u +u )2 + 1(3u (cid:0)4u +u )2
12 j j+1 j+2 4 j j+1 j+2
(cid:12)RR = 13(u (cid:0)2u +u )2 + 1((cid:0)5u +8u (cid:0)3u )2:
12 j+1 j+2 j+3 4 j+1 j+2 j+3
The (cid:28) also varies with discretization order; the expression for (cid:28) is given by
p 5
((cid:0)f +5f (cid:0)10f +10f (cid:0)5f +f )2 ; for ’6=0
(cid:28) = j(cid:0)2 j(cid:0)1 j j+1 j+2 j+3 (20)
5 ( (fj(cid:0)2(cid:0)4fj(cid:0)1+6fj (cid:0)4fj+1+fj+2)2 ; for ’=0
Expressionsfor the fourth-, seventh-, and eighth-order WENO schemes appear in the appendix.
aTheaboveapproachtoconstructingcentralWENOschemesisproposedinreference.6 Ingeneral,centralhigh-orderWENO
schemesthatarebuiltinthismanner,areunstablewhenunresolvedfeaturesorstrongdiscontinuitiesexistinthecomputational
domain. The generalized energy-stable methodology presented in the next section guarantees the stability of the even order
approximationswhilemaintainingtheirnonoscillatoryproperties.
6of30
AmericanInstitute ofAeronautics andAstronautics
The WENO mechanics expressed in equations (17){(20) represents a signi(cid:12)cant departure from the
mechanics used in the original algorithm by Jiang and Shu.2 For example, equation (18) replaces the
expression
d(r)
(cid:11) = (21)
r ((cid:15)+(cid:12)(r))2
that is usedin the weight functions givenin reference.2 Furthermore,(cid:15)(cid:20)O((cid:1)x2) replaces(cid:15)=O(10(cid:0)6) in.2
Equations (17){(20) more closely resemble the weight and smoothness indicators proposed by Borges et al.
in8 There are di(cid:11)erences, however, in how the parameters (cid:28) and (cid:15) are chosen; di(cid:11)erences motivated by the
p
following observations.
Borgeset al.8 proposes the following function (cid:28)~ :
5
(cid:28)~ =j(cid:12)LL(cid:0)(cid:12)Rj (22)
5
anda(cid:12)xedvaluefortheparameter(cid:15). Althoughtheaboveweightsandsmoothnessindicators[eqs. (17){(19)]
with the value of (cid:28)~ given by equation (22) signi(cid:12)cantly outperform the conventional weights of Jiang and
5
Shu2 for both continuous and discontinuous solutions, three remaining shortcomingsinclude:
1. Thevalueof(cid:28)~ givenbyequation(22)isnotsmoothbecauseanabsolutevalue functionisused, which
5
may reduce accuracy at points where (cid:12)LL(cid:0)(cid:12)R changes sign.
2. This approach currently does not generalize to WENO schemes with design orders other than 5. For
example, neither (cid:28)~ =j(cid:12)LL(cid:0)(cid:12)Rj nor (cid:28)~ =j(cid:12)LL(cid:0)(cid:12)RRj provides the design order of accuracy for the
6 6
central sixth-order WENO scheme.
3. Near critical points where f0, f00, and f000 vanish simultaneously, the modi(cid:12)ed weights [eqs. (17){
(19), and (22)] fail to provide the design order of convergence for the (cid:12)fth-order WENO scheme if no
constraints are imposed on the parameter (cid:15) other than (cid:15) > 0. An example that demonstrates this
property is presented in Section 6.1.
Selecting(cid:28) inequation(20)tobeaquadraticfunctionofthe(cid:12)fth-degreeundivideddi(cid:11)erencede(cid:12)nedon
5
the entiresix-pointstencil, circumventsthe (cid:12)rst,and seconddrawbacksencounteredwhenusing the scheme
suggested in reference8 . Indeed, (cid:28) is a C1 function in its arguments, and can readily be generalized to
5
WENO schemes of order p by choosing (cid:28) to be the pth-degree undivided di(cid:11)erence that is de(cid:12)ned on the
p
entire (p+1){point stencil. In contrast to the (cid:28)~ found in reference,8 which is of order O((cid:1)x5) for smooth
5
solutions, the proposed function (cid:28) as given by equation (20) is of order O((cid:1)x8) and O((cid:1)x10) for ’ = 0
5
and ’ 6= 0, respectively; thus much faster convergence of the (cid:12)fth- and sixth-order WENO schemes to the
corresponding underlying linear schemes is achieved. b
A remedy for the third shortcoming requires an additional modi(cid:12)cation to the approach proposed in
reference.8 Boththe(cid:12)fth-andsixth-orderWENOschemeswiththenewweightsthataregivenbyequations
(17{20)aredesign-orderaccurateforsmoothsolutions,includingpointsatwhichthe (cid:12)rst-andsecond-order
derivatives of the solution vanish simultaneously. However, if all derivatives up to fourth are equal to zero,
then the (cid:12)fth- and sixth-order WENO schemes locally become only third-order accurate. To fully resolve
this issue, the constraint (cid:15)(cid:20)O((cid:1)x2) is imposed herein.
Equations(13){(20)describeanewfamilyof(cid:12)fth-orderWENOschemes. Theprimarydi(cid:11)erencebetween
existingWENOschemesandthisnewfamily,isinthestencilbiasingmechanicsdescribedbyequations(17){
(20). Detailed theoretical justi(cid:12)cations for the choice of (cid:28) and (cid:15) used in equations (17){(20) are presented
p
in sections 5.2, 5.3. There, and again in the results section, it is shown that these parameters provide fast
convergence of the new WENO and ESWENO schemes to the corresponding underlying linear schemes for
smooth solutions, and deliver improved shock-capturingcapabilities near unresolved features.
bThe conventional (cid:12)fth-order WENO scheme that corresponds to ’ = 0, does not include the downwind stencil SRR;
therefore, the fourth-degree undivided di(cid:11)erence built on the available (cid:12)ve-point stencil, is included in equation (20). This
approachtohandlethenarrowstencilencountered with’=0,alsogeneralizestootherordersofaccuracy.
7of30
AmericanInstitute ofAeronautics andAstronautics
IV. A General Approach to Constructing High-Order Energy-Stable Schemes
Algorithm (1) transforms an existing WENO scheme into an ESWENO scheme that is characterized
by the following four properties: 1) a bounded energy estimate for arbitrary nonsmooth initial data, 2)
conservation, 3) design order accuracy for su(cid:14)ciently smooth data, and 4) discontinuity (shock) capturing
capabilities that are similar to those of the base WENO scheme.
Algorithm 1 ( Transformation from WENO to ESWENO).
1: Express the base WENO derivative operator as matrix D.
2: Decompose D into symmetric and skew-symmetric components as D=D +D :
sk sym
3: Add an arti(cid:12)cial dissipation operator D such that the modi(cid:12)ed symmetric matrix
ad
D +D is positive semide(cid:12)nite.
sym ad
3a: Form the decomposition D = s [D i](cid:3) [D i]T,
sym i=0 1 i 1
2s+1 is the bandwidth of the matrix D.
P
The existence of this decomposition is established in the appendix.
3b: Modify the diagonal terms of (cid:3) such that they are smoothly positive.
i
One such approach is
((cid:21)(cid:22) ) = 1 ((cid:21) )2 +(cid:14)2 + ((cid:21) ) (23)
i j;j 2 i j;j i i j;j
where (cid:14) , 0(cid:20)hiq(cid:20)s are small positive cionstants that may depend on (cid:1)x.
i
3c: Form the symmetric, positive semide(cid:12)nite matrix
D = P(cid:0)1 s [D i](cid:3)(cid:22) [D i]T (24)
ad i=0 1 i 1
4: Form the energy-stable operator
P
D(cid:22) = D + D (25)
ad
Byconstruction,algorithm(1)isapplicableto1-Dperiodicschemesofanyorderandproducesamodi(cid:12)ed
schemethatautomaticallysatis(cid:12)esthe(cid:12)rstandsecondproperties: boundedenergyestimateandconservation
(See reference1 for details.) Likewise, as shown in the next section, the additional terms added in the
ESWENO formulation do not degrade the formal accuracy of the original WENO discretization (property
three). However, the degree to which the two formulations di(cid:11)er for unresolved data is not clear. Thus, the
properties of the new ESWENO scheme must be tested to ensure it retains the desirable properties of the
original formulation.
ESWENO Schemes of Fifth- and Sixth-order
We now apply algorithm (1) to the one-parameter family of (cid:12)fth-order WENO schemes [eqs. (17){(20)] for
the scalar 1-D wave equation (1) with a(cid:21)0. Combining equation (15) with equation (14) and substituting
the resulting WENO (cid:13)ux f^ into equation (13) produces a stencil for the jth grid point of the form
j+1
2
0 0 2wjL+L1=2 (cid:0)7wjL+L1=2 11wjL+L1=2 0 0 0 1
Dj5;k = 6(cid:1)1xBBBBBBBBBBBBBBBBBBBB@ (cid:0)2wjL000000(cid:0)L1=2 17wwjjLL(cid:0)(cid:0)00000L11==22 (cid:0)(cid:0)(cid:0)(cid:0)12511wwwwjjLRjL000jL(cid:0)(cid:0)+(cid:0)L1111====2222 (cid:0)(cid:0)(cid:0)25125ww1wwwjjRLjjLR++00jR(cid:0)(cid:0)(cid:0)11R11==1==22=222 117251wwwwwjjjjRRLRjR(cid:0)(cid:0)++00+RR11111=====22222 (cid:0)(cid:0)(cid:0)217wwwjRjjRR0000(cid:0)++RR111===222 2wjR+000000R1=2 CCCCCCCCCCCCCCCCCCCCA0BBBBBBBBBBBBB@ fffffff((((((uuuuuu(ujjjjjj(cid:0)(cid:0)(cid:0)+++j123321))))))) 1CCCCCCCCCCCCCA (26)
with k on the interval j(cid:0)3 (cid:20) k (cid:20) j+3. An explicit expression for the di(cid:11)erentiation matrix D5 follows
immediately from equation (26).
8of30
AmericanInstitute ofAeronautics andAstronautics
The derivative matrix D5 is now decomposed into symmetric and skew-symmetric parts as
D5 = D5 + D5 :
skew sym
As with the third-order case (ref.1), the skew-symmetric component of D5 and the norm, take the forms
D5 = P(cid:0)1Q ; Q +QT =0 ; P =(cid:1)xI :
skew 5 5 5
while the matrix D5 is expressed as
sym
D5 =P(cid:0)1 D3(cid:3)5[D3]T + D2(cid:3)5[D2]T + D1(cid:3)5[D1]T + D0(cid:3)5[D0]T : (27)
sym 1 3 1 1 2 1 1 1 1 1 0 1
(cid:16) (cid:17)
The matrices (cid:3)5 are diagonal with expressions for the jth element de(cid:12)ned by
j
((cid:21)5) = 1 wLL (cid:0)wRR (28)
3 j;j 6 j+5=2 j+1=2
h i
+wLL (cid:0)4wLL +wL
((cid:21)5) = 1 j+3=2 j+5=2 j+3=2 (29)
2 j;j 12" (cid:0)wR +4wRR (cid:0)wRR #
j+1=2 j(cid:0)1=2 j+1=2
3wLL (cid:0)5wLL +2wLL
j+1=2 j+3=2 j+5=2
((cid:21)5) = 1 2 +wjL+1=2 (cid:0)wjL+3=2 3 (30)
1 j;j 12 +wR (cid:0)wR
6 j(cid:0)1=2 j+1=2 7
6 7
6 (cid:0)2wRR +5wRR (cid:0)3wRR 7
6 j(cid:0)3=2 j(cid:0)1=2 j+1=2 7
4 5
(cid:0)wLL (cid:0)wL (cid:0)wR (cid:0)wRR
((cid:21)5) = 1 j(cid:0)1=2 j(cid:0)1=2 j(cid:0)1=2 j(cid:0)1=2 =0 (31)
0 j;j 2" +wLL +wL +wR +wRR #
j+1=2 j+1=2 j+1=2 j+1=2
Because w(r) =1,(cid:3)5 is alwaysequal tozero,andis notincluded inD5 . The signof the diagonalterms
r : 0 ad
((cid:21) 5) ; i=1;2;3;could be either positive or negative;thus, the conventional (cid:12)fth- and sixth-orderWENO
i j;j P
schemes may become locally unstable. If we de(cid:12)ne ((cid:21)(cid:22)5) , i=1;2;3; to be smoothly positive
i j;j
1
((cid:21)(cid:22)5) = [ ((cid:21)5)2 +(cid:14)2 + ((cid:21)5) ] ;
i j;j 2 i j;j i i j;j
q
then the additional arti(cid:12)cial dissipation operator becomes D(cid:22)5 = P(cid:0)1 s [D i]((cid:3)(cid:22)5)[D i]T, and the re-
ad i=1 1 i 1
sulting energy-stable scheme is obtained by adding the additional dissipation term to the original WENO
P
scheme. That is,
D(cid:22)5 = D5 + D(cid:22)5 : (32)
ad
By construction, D5 satis(cid:12)es allof the conditions of Theorem1, therebyprovidingstabilityof the (cid:12)fth- and
sixth-order ESWENO schemes that are de(cid:12)ned by equation (32).
V. Consistency Analysis
V.A. Necessary and Su(cid:14)cient Conditions for Consistency of WENO Schemes
We now derive the necessary and su(cid:14)cient conditions for the weight functions w(r) for a family of pth-
:
order WENO schemes to attain the design order of accuracy. Similar conditions have been obtained for the
conventional (cid:12)fth-order WENO scheme in reference.7 With the same approach discussed in section 3 for
p=5, aone-parameterfamilyof pth-orderWENO(cid:13)uxescanbeconstructedbyusingaconvexcombination
(r)
of (s+1) (cid:13)uxes f of sth order as
:
f^ = w(r) f(r) ; (33)
j(cid:6)1=2 j(cid:6)1=2 j(cid:6)1=2
r
X
9of30
AmericanInstitute ofAeronautics andAstronautics
with
p
f(r) =h(x )+ c(r)(cid:1)xl+O((cid:1)xp+1); (34)
j(cid:6)1=2 j(cid:6)1=2 l
l=s
X
where h(x) is the numerical (cid:13)ux function that is implicitly de(cid:12)ned as
x+(cid:1)x
2
1
f(x)= h((cid:17))d(cid:17) (35)
(cid:1)x
Z
x(cid:0)(cid:1)x
2
and c(r) are constants that do not depend on (cid:1)x.
l
Thecorrespondingpth-orderWENOoperatorthatapproximatesthe(cid:12)rst-orderspatialderivativeisgiven
by
w(r) f(r) (cid:0)w(r) f(r)
f^ (cid:0)f^ j+1=2 j+1=2 j(cid:0)1=2 j(cid:0)1=2
[Dpf]j = j+1=2(cid:1)x j(cid:0)1=2 = Pr (cid:16) (cid:1)x (cid:17); (36)
where[(cid:1)] isajthcomponentofavector,theindexrinequation(36)sweepsoverall(s+1)stencils,andw(r)
j
isanonlinearweightfunctionthatisassignedtothecorrespondings-pointstencilS . Forsu(cid:14)cientlysmooth
r
solutions, the weights w(r) approach their preferred values d(r), so that the WENO operator converges to
:
the target linear operator DTarget as
fTarget(cid:0)fTarget d(r)fj(+r)1=2(cid:0)d(r)fj((cid:0)r)1=2
DTargetf = j+1=2 j(cid:0)1=2 = r (cid:16) (cid:17); (37)
j (cid:1)x P (cid:1)x
(cid:2) (cid:3)
where d(r) is a one-parameter family of constants chosen to ensure pth-order convergence of the target
operator to the exact value of the (cid:12)rst-order derivative at x . That is,
j
@f
DTargetf = +O((cid:1)xp): (38)
j @x
(cid:12)x=xj
(cid:2) (cid:3) (cid:12)
(cid:12)
The coe(cid:14)cients d(r) for one-parameter families of linea(cid:12)r target schemes up to eighth order are given in the
appendix. All target linear schemes in the family are (2s(cid:0)1)th-order accurate, except one central scheme,
which is (2s)th-order accurate.
Subtracting equation (37) from equation (36) and using equation (34), we have
[Dpf] (cid:0) DTargetf = Prh(wj(r+)1=2(cid:0)d(r))fj(+r)1=2(cid:0)(wj(r(cid:0))1=2(cid:0)d(r))fj((cid:0)r)1=2i
j j (cid:1)x
= (cid:1)1x (cid:2)(wj(r+)1=2(cid:0)(cid:3) d(r))h(r)(xj+1=2)(cid:0)(wj(r(cid:0))1=2 (cid:0)d(r))h(r)(xj(cid:0)1=2) (39)
r
p h i
+ P(cid:1)xl(cid:0)1c(r) (w(r) (cid:0)d(r))(cid:0)(w(r) (cid:0)d(r)) +O((cid:1)xp)
l j+1=2 j(cid:0)1=2
l=s r
h i
PP
From equation (39) it immediately follows that to retain pth-order accuracy, the weights of the WENO
operator Dp should satisfy the following necessary and su(cid:14)cient conditions:
w(r) (cid:0)d(r) =O((cid:1)xp+1)
j(cid:6)1=2
r
Phc(r) (w(r) (cid:0)id(r))(cid:0)(w(r) (cid:0)d(r)) =O((cid:1)xp(cid:0)s+1)
r s j+1=2 j(cid:0)1=2 (40)
h i
P:::
c(r) (w(r) (cid:0)d(r))(cid:0)(w(r) (cid:0)d(r)) =O((cid:1)x)
p j+1=2 j(cid:0)1=2
r
h i
P
Here, we use the following properties of h(x) and d(r) : f0(x)= h(xj+1=2)(cid:0)h(xj(cid:0)1=2), and d(r) =1.
(cid:1)x
r
P
10of30
AmericanInstitute ofAeronautics andAstronautics
