03964.com

文档资料库 文档搜索专家

文档资料库 文档搜索专家

JOURNAL OF COMPUTATIONAL PHYSICS ARTICLE NO.

126, 202–228 (1996)

0130

Ef?cient Implementation of Weighted ENO Schemes*

GUANG-SHAN JIANG?

AND

CHI-WANG SHU

Division of Applied Mathematics, Brown University, Providence, Rhode Island 02912 Received August 14, 1995; revised January 3, 1996

In this paper, we further analyze, test, modify, and improve the high order WENO (weighted essentially non-oscillatory) ?nite difference schemes of Liu, Osher, and Chan. It was shown by Liu et al. that WENO schemes constructed from the rth order (in L1 norm) ENO schemes are (r 1)th order accurate. We propose a new way of measuring the smoothness of a numerical solution, emulating the idea of minimizing the total variation of the approximation, which results in a ?fth-order WENO scheme for the case r 3, instead of the fourth-order with the original smoothness measurement by Liu et al. This ?fth-order WENO scheme is as fast as the fourth-order WENO scheme of Liu et al. and both schemes are about twice as fast as the fourth-order ENO schemes on vector supercomputers and as fast on serial and parallel computers. For Euler systems of gas dynamics, we suggest computing the weights from pressure and entropy instead of the characteristic values to simplify the costly characteristic procedure. The resulting WENO schemes are about twice as fast as the WENO schemes using the characteristic decompositions to compute weights and work well for problems which do not contain strong shocks or strong re?ected waves. We also prove that, for conservation laws with smooth solutions, all WENO schemes are convergent. Many numerical tests, including the 1D steady state nozzle ?ow problem and 2D shock entropy wave interaction problem, are presented to demonstrate the remarkable capability of the WENO schemes, especially the WENO scheme using the new smoothness measurement in resolving complicated shock and ?ow structures. We have also applied Yang’s arti?cial compression method to the WENO schemes to sharpen contact discontinuities. 1996 Academic Press, Inc.

1. INTRODUCTION

In this paper, we further analyze, test, modify, and improve the WENO (weighted essentially non-oscillatory) ?nite difference schemes of Liu, Osher, and Chan [9] for the approximation of hyperbolic conservation laws of the type ut div f(u) 0, (1.1)

* Research supported by ARO Grant DAAH04-94-G-0205, NSF Grants ECS-9214488 and DMS-9500814, NASA Langley Grant NAG1-1145 and Contract NAS1-19480 while the second author was in residence at ICASE, NASA Langley Research Center, Hampton, VA 236810001, and AFOSR Grant 95-1-0074. ? Current address: Department of Mathematics, UCLA, Los Angeles, CA 90024. 202

0021-9991/96 $18.00 Copyright ? 1996 by Academic Press, Inc. All rights of reproduction in any form reserved.

or perhaps with a forcing term g(u, x, t) on the right-hand side. Here u (u1 , ..., um), f (f1 , ..., fd), x (x1 , ..., xd) and t 0. WENO schemes are based on ENO (essentially nonoscillatory) schemes, which were ?rst introduced by Harten, Osher, Engquist, and Chakravarthy [5] in the form of cell averages. The key idea of ENO schemes is to use the ‘‘smoothest’’ stencil among several candidates to approximate the ?uxes at cell boundaries to a high order accuracy and at the same time to avoid spurious oscillations near shocks. The cell-averaged version of ENO schemes involves a procedure of reconstructing point values from cell averages and could become complicated and costly for multi-dimensional problems. Later, Shu and Osher [14, 15] developed the ?ux version of ENO schemes which do not require such a reconstruction procedure. We will formulate the WENO schemes based on this ?ux version of ENO schemes. The WENO schemes of Liu et al. [9] are based on the cell-averaged version of ENO schemes. For applications involving shocks, second-order schemes are usually adequate if only relatively simple structures are present in the smooth part of the solution (e.g., the shock tube problem). However, if a problem contains rich structures as well as shocks (e.g., the shock entropy wave interaction problem in Example 4, Section 8.3), high order shock capturing schemes (order of at least three) are more ef?cient than low order schemes in terms of CPU time and memory requirements. ENO schemes are uniformly high order accurate right up to the shock and are very robust to use. However, they also have certain drawbacks. One problem is with the freely adaptive stencil, which could change even by a round-off perturbation near zeroes of the solution and its derivatives. Also, this free adaptation of stencils is not necessary in regions where the solution is smooth. Another problem is that ENO schemes are not cost effective on vector supercomputers such as CRAY C-90 because the stencil-choosing step involves heavy usage of logical statements, which perform poorly on such machines. The ?rst problem could reduce the accuracy of ENO schemes for certain functions [12]; however, this can be remedied by embedding certain parameters (e.g., threshold and biasing factor) into the

WEIGHTED ENO SCHEMES

203

stencil choosing step so that the preferred linearly stable stencil is used in regions away from discontinuities. See [1, 3, 13]. The WENO scheme of Liu, Osher, and Chan [9] is another way to overcome these drawbacks while keeping the robustness and high order accuracy of ENO schemes. The idea is the following: instead of approximating the numerical ?ux using only one of the candidate stencils, one uses a convex combination of all the candidate stencils. Each of the candidate stencils is assigned a weight which determines the contribution of this stencil to the ?nal approximation of the numerical ?ux. The weights can be de?ned in such a way that in smooth regions it approaches certain optimal weights to achieve a higher order of accuracy (an rth-order ENO scheme leads to a (2r 1)th-order WENO scheme in the optical case), while in regions near discontinuities, the stencils which contain the discontinuities are assigned a nearly zero weight. Thus essentially non-oscillatory property is achieved by emulating ENO schemes around discontinuities and a higher order of accuracy is obtained by emulating upstream central schemes with the optimal weights away from the discontinuities. WENO schemes completely remove the logical statements that appear in the ENO stencil choosing step. As a result, the WENO schemes run at least twice as fast as ENO schemes (see Section 7) on vector machines (e.g., CRAY C-90) and are not sensitive to round-off errors that arise in actual computation. Atkins [1] also has a version of ENO schemes using a different weighted average of stencils. Another advantage of WENO schemes is that its ?ux is smoother than that of the ENO schemes. This smoothness enables us to prove convergence of WENO schemes for smooth solutions using Strang’s technique [18]; see Section 6. According to our numerical tests, this smoothness also helps the steady state calculations, see Example 4 in Section 8.2. In [9], the order of accuracy shown in the error tables (Table 1–5 in [9]) seemed to suggest that the WENO schemes of Liu et al. are more accurate than what the truncation error analysis indicated. In Section 2, we carry out a more detailed error analysis for the WENO schemes and ?nd that this ‘‘super-convergence’’ is indeed super?cial: the ‘‘higher’’ order is caused by larger error on the coarser grids instead of smaller error on the ?ner grids. Our error analysis also suggests that the WENO schemes can be made more accurate than those in [9]. Since the weight on a candidate stencil has to vary according to the relative smoothness of this stencil to the other candidate stencils, the way of evaluating the smoothness of a stencil is crucial in the de?nition of the weight. In Section 3, we introduce a new way of measuring the smoothness of the numerical solution which is based upon minimizing the L2 norm of the derivatives of the reconstruction polynomials, emulating the idea of minimizing

the total variation of the approximations. This new measurement gives the optimal ?fth-order accurate WENO scheme when r 3 (the smoothness measurement in [9] gives a fourth-order accurate WENO scheme for r 3). Although the WENO schemes are faster than ENO schemes on vector supercomputers, they are only as fast as ENO schemes on serial computers. In Section 4, we present a simpler way of computing the weights for the approximation of Euler systems of gas dynamics. The simpli?cation is aimed at reducing the ?oating point operations in the costly but necessary characteristic procedure and is motivated by the following observation: the only nonlinearity of a WENO scheme is in the computation of the weights. We suggest using pressure and entropy to compute the weights, instead of the local characteristic quantities. In this way one can exploit the linearity of the rest of the scheme. The resulting WENO schemes (r 3) is about twice as fast as the original WENO scheme which uses local characteristic quantities to compute the weights (see Section 7). The same idea can also be applied to the original ENO schemes. Namely, we can use the undivided differences of pressure and entropy to replace the local characteristic quantities to choose the ENO stencil. This has been tested numerically but the results are not included in this paper since the main topic here is the WENO schemes. WENO schemes have the same smearing at contact discontinuities as ENO schemes. There are mainly two techniques for sharpening the contact discontinuities for ENO schemes. One is Harten’s subcell resolution [4] and the other is Yang’s arti?cial compression (slope modi?cation) [20]. Both were introduced in the cell average context. Later, Shu and Osher [15] translated them into the point value framework. In one-dimensional problems, the subcell resolution technique works slightly better than the arti?cial compression method. However, for two or higher dimensional problems, the latter is found to be more effective and handy to use [15]. We will highlight the key procedures of applying the arti?cial compression method to the WENO schemes in Section 5. In Section 8, we test the WENO schemes (both the WENO schemes of Liu et al. and the modi?ed WENO schemes) on several 1D and 2D model problems and compare them with ENO schemes to examine their capability in resolving shock and complicated ?ow structures. We conclude this paper by a brief summary in Section 9. The time discretization of WENO schemes will be implemented by a class of high order TVD Runge–Kutta-type methods developed by Shu and Osher [14]. To solve the ordinary differential equation

du dt

L(u),

(1.2)

204

JIANG AND SHU

where L(u) is a discretization of the spatial operator, the third-order TVD Runge–Kutta is simply u(1) u(2) u

n 1

TABLE I Coef?cients ar k,l

r k 0 1 0 1 2 l 0 1/2 1/2 1/3 1/6 1/3 l 1 3/2 1/2 7/6 5/6 5/6 11/6 1/3 1/6 l 2

un un u

n

tL(un) u(1) u

(2)

tL(u(1)) tL(u ).

(2)

(1.3)

2

3

Another useful, although not TVD, fourth-order Runge– Kutta scheme is u(1) u(2) u(3) un

1

un un un ( un

tL(un) tL(u(1)) tL(u(2)) u(1) 2u(2) u(3)) tL(u(3)). (1.4) where df (u) /du one can de?ne f (u) f (u) f (u), (2.4) 0. For example, (2.5) 0 and df (u) /du f (u) ( f (u) u),

This fourth-order Runge–Kutta scheme can be made TVD by an increase of operation counts [14]. We will mainly use these two Runge–Kutta schemes in our numerical tests in Section 8. The third-order TVD Runge–Kutta scheme will be referred to as ‘‘RK-3’’ while the fourth-order (nonTVD) Runge–Kutta scheme will be referred to as ‘‘RK-4.’’

2. THE WENO SCHEMES OF LIU, OSHER, AND CHAN

where max f (u) and the maximum is taken over the whole relevant range of u. This is the global Lax–Friedrichs (LF) ?ux splitting. For other ?ux splitting, especially the Roe ?ux splitting with entropy ?x (RF); see [15] for details. Let ? j 1/2 and f j 1/2 be, resp. the numerical ?uxes obtained f from the positive and negative parts of f (u), we then have

? fj

1/2

?j f

1/2

?j f

1/2 .

(2.6)

In this section, we use the ?ux version of ENO schemes as our basis to formulate WENO schemes of Liu et al. and analyze their accuracy in a different way from that used in [9]. We use one-dimensional scalar conservation laws (i.e., d m 1 in (1.1)) as an example: ut f (u)x 0. (2.1)

Let us discretize the space into uniform intervals of size j x. Various quantities at xj will be x and denote xj identi?ed by the subscript j. The spatial operator of the WENO schemes, which approximates f (u)x at xj , will take the conservative form L 1 ? ( fj x

? Here we will only describe how f j 1/2 is computed in [9] on the basis of ?ux version of ENO schemes. For simplicity, we will drop the ‘‘ ’’ sign in the superscript. The formulas for the negative part of the split ?ux are symmetric (with respect to xj 1/2) and will not be shown. As we know, the rth-order (in L1 sense) ENO scheme chooses one ‘‘smoothest’’ stencil from r candidate stencils and uses only the chosen stencil to approximate the ?ux hj 1/2 . Let’s denote the r candidate stencils by Sk , k 0, 1, ..., r 1, where

Sk (xj

k r 1,

xj

k r 2,

..., xj k).

1/2

?j f

1/2),

(2.2)

If the stencil Sk happens to be chosen as the ENO interpolation stencil, then the rth-order ENO approximation of hj 1/2 is

where the numerical ?ux ?j 1/2 approximates hj 1/2 f h(xj 1/2) to a high order with h(x) implicitly de?ned by [15] f (u(x)) 1 x

x x x/2 x/2

?j f

where

1/2

qr ( fj k

k r 1,

..., fj k),

(2.7)

r 1

h( ) d .

(2.3)

qr (g0 , ..., gr 1) k

l 0

a r gl . k,l

(2.8)

We can actually assume f (u) 0 for all u in the range of our interest. For a general ?ux, i.e., f (u) 0, one can split it into two parts either globally or locally,

k, l r 1, are constant coef?cients. For Here a r , 0 k,l later use, we provide these coef?cients for r 2, 3, in Table I.

WEIGHTED ENO SCHEMES

205

TABLE II Optimal Weights C r k

To just use the one smoothest stencil among the r candidates for the approximation of hj 1/2 , is very desirable near discontinuities because it prohibits the usage of information on discontinuous stencils. However, it is not so desirable in smooth regions because all the candidate stencils carry equally smooth information and thus can be used together to give a higher order (higher than r, the order of the base ENO scheme) approximation to the ?ux hj 1/2 . In fact, one could use all the r candidate stencils, which all together contain (2r 1) grid values of f to give a (2r 1)th-order approximation of hj 1/2 :

Cr k r r 2 3

k

0

k

1

k

2

1/3 1/10

2/3 6/10

— 3/10

?j f

1/2

q2r 11( fj r

r 1

r 1,

..., fj

r 1)

(2.12) (

k

?j f

1/2

q2r 11( fj r

r 1,

..., fj

r 1)

(2.9)

C r )qr ( fj k k

k r 1,

..., fj k).

k 0

which is just the numerical ?ux of a (2r 1)th-order upstream central scheme. As we know, high order upstream central schemes (in space), combined with high order Runge–Kutta methods (in time), are stable and dissipative under appropriate CFL numbers and thus are convergent, according to Strang’s convergence theory [18] when the solution of (1.1) is smooth (see Section 6). The above facts suggest that one could use the (2r 1)th-order upstream central scheme in smooth regions and only use the rthorder ENO scheme near discontinuities. As in (2.7), each of the stencils can render an approximation of hj 1/2 . If the stencil is smooth, this approximation is rth-order accurate; otherwise, it is less accurate or even not accurate at all if the stencil contains a discontinuity. One could assign a weight k to each candidate stencil Sk , k 0, 1, ..., r 1, and use these weights to combine the r different approximations to obtain the ?nal approximation of hj 1/2 as

Recalling (2.9), we see that, the ?rst term on the righthand side of the above equation is a (2r 1)th-order r 1 1, if we require approximation of hj 1/2 . Since k 0 C r k r 1 1, the last summation term can be written as k k 0

r 1

(

k 0

k

C r )(qr ( fj k k

k r 1,

..., fj k)

hj

1/2).

(2.13)

Each term in the last summation can be made O(h2r 1) if

k

Cr k

O(hr 1)

(2.14)

for k 0, 1, ..., r 1. Here, h x. Thus C r will bear k the name of optimal weight. The question now is how to de?ne the weight such that (2.14) is satis?ed in smooth regions while essentially nonoscillatory property is achieved. In [9], the weight k for stencil Sk is de?ned by

k k 0 r 1

? fj

r 1 1/2 k 0 r kqk( fj k r 1 ,

..., fj k),

(2.10) where

,

(2.15)

where qr ( fj k r 1 , ..., fj k) is de?ned in (2.8). To achieve k essentially non-oscillatory property, one then requires the weight to adapt to the relative smoothness of f on each candidate stencil such that any discontinuous stencil is effectively assigned a zero weight. In smooth regions, one can adjust the weight distribution such that the resulting ? approximation of the ?ux fj 1/2 is as close as possible to that given in (2.9). Simple algebra gives the coef?cients C r such that k

r 1

k

(

Cr k , k ISk)p

0, 1, ..., r

1.

(2.16)

q2r 11( fj r

r 1

r 1,

..., fj

r 1) k 0

C r qr ( fj k k

k r 1,

..., fj k)

(2.11)

and k 0 C r 1 for all r 2. For r 2, 3, these coef?cients k are given in Table II. Comparing (2.11) with (2.10), we get

Here is a positive real number number which is introduced to avoid the denominator becoming zero (in our later tests, we will take 10 6. Our numerical tests indicate that the result is not sensitive to the choice of , as long as it is in the range of 10 5 to 10 7); the power p will be discussed in a moment; ISk in (2.16) is a smoothness measurement of the ?ux function on the kth candidate stencil. r 1 It is easy to see that k 0 k 1. To satisfy (2.14), it suf?ces to have (through a Taylor expansion analysis) ISk D(1 O(hr 1)) (2.17)

206

JIANG AND SHU

for k 0, 1, ..., r 1, where D is some non-zero quantity independent of k. As we know, an ENO scheme chooses the ‘‘smoothest’’ ENO stencil by comparing a hierarchy of undivided differences. This is because these undivided differences can be used to measure the smoothness of the numerical ?ux on a stencil. In [9], ISk is de?ned as

r 1 r l

IS1

(( f h ( f h2)2

f h2)2 O(h5) f h2)2 O(h5).

(f h

f h2)2) (2.24)

IS2

(( f h ( f h2)2

(f h

f h2)2) (2.25)

ISk

l 1 i 1

(f[j

k r

i l

r, l])2

,

(2.18)

where f [ , ] is the lth undivided difference: f [ j, 0] f [ j, l] fj f[j 1, l 1] f [ j, l 1], k 1, ..., r 1.

For example, when r ISk When r ISk (f[j k

2, we have 1, 1])2, k 0, 1. (2.19)

3, (2.18) gives (( f [ j (f[j k k 2, 1])2 (f[j k 0, 1, 2. 1, 1])2) (2.20)

2, 2])2, k

In smooth regions, Taylor expansion analysis of (2.18) gives ISk ( f h)2(1 O(h)), k 0, 1, ..., r 1,

r 1

We can see that the second-order terms are different from stencil to stencil. Thus (2.22) is no longer valid at critical points of f (u(x)) which implies that the WENO scheme of Liu et al. for r 3 is only third-order accurate at these points. In fact, the weights computed from the smoothness measurement (2.18) diverge far away from the optimal weights near critical points (see Fig. 1 in the next section) on coarse grids (10 to 80 grid points per wave). But on ?ne grids, since the smoothness measurements ISk for all k are relatively smaller than the non-zero constant in (2.16), the weights become close to the optimal weights. Therefore the ‘‘super-convergence’’ phenomena appeared in Tables 1–5 in [9] are caused by large error commitment on coarse grids and less error commitment on ?ner grids when using the errors of the ?fth-order central scheme as reference (see Tables III and IV). At discontinuities, it is typical that one or more of the r candidate stencils reside in smooth regions of the numerical solution while other stencils contain the discontinuities. The size of the discontinuities is always O(1) and does not change when the grid is re?ned. So we have for a smooth stencil Sk , ISk O(h2p) (2.26)

(2.21)

where f f (uj). Note the O(h) term is not O(h ) that we would want to have (see (2.17)). Thus in smooth monotone regions, i.e., f 0, we have

k

and for a non-smooth stencil Sl , ISl O(1). (2.27)

Cr k

O(h), k

0, 1, ..., r

1.

(2.22)

Recalling (2.12), we see that the WENO schemes with the smoothness measurement given by (2.18) is (r 1)th-order accurate in smooth monotone regions of f (u(x)). This result was proven in [9] using a different approach. For r 2, this is optimal in the sense that the third-order upstream central scheme is approximated in most smooth regions. However, this is not optimal for r 3, for which this measurement can only give fourth-order accuracy while the optimal upstream central scheme is ?fth-order accurate. We will introduce a new measurement in the next section which will result in an optimal order accurate WENO scheme for the r 3 case. When r 3, Taylor expansion of (2.20) gives IS0 (( f h ( f h2)2 f h2)2 O(h5) (f h f h2)2) (2.23)

From the de?nition of the weights (2.15), we can see that, for this non-smooth stencil Sl , the corresponding weight l satis?es

l

O(h2p).

(2.28)

Therefore for small h and any positive integer power p, the weight assigned to the non-smooth stencil vanishes as h 0. Note, if there is more than one smooth stencils in the r candidates, from the de?nition of the weights in (2.15), we expect each of the smooth stencils will get a weight which is O(1). In this case, the weights do not exactly resemble the ‘‘ENO digital weights.’’ However, if a stencil is smooth, the information that it contains is useful and should be utilized. In fact, in our extensive numerical experiments, we ?nd the WENO schemes in [9] work very well at shocks. We also ?nd that p 2 is adequate to obtain essentially non-oscillatory approximations at least

WEIGHTED ENO SCHEMES

207

1

for r 2, 3, although it is suggested in [9] that p should be taken as r, the order of the base ENO schemes. We will use p 2 for all our numerical tests. Notice that, in discussing accuracy near discontinuities, we are simply concerned with spatial approximation error. The error due to time evolution is not considered. In summary, WENO schemes of Liu et al. de?ned by (2.10), (2.15), and (2.18) have the following properties: 1. They involve no logical statements which appear in the base ENO schemes. 2. The WENO scheme based on the rth-order ENO scheme is (r 1)th-order accurate in smooth monotone regions, although this is still not as good as the optimal order (2r 1). 3. They achieve essentially non-oscillatory property by emulating ENO schemes at discontinuities. 4. They are smooth in the sense that the numerical ?ux ?j 1/2 is a smooth function of all its arguments (for a general f ?ux, this is also true if a smooth ?ux splitting method is used, e.g., global Lax–Friedrichs ?ux splitting).

3. A NEW SMOOTHNESS MEASUREMENT

IS1 IS2

( fj ( fj

2 fj 2 fj

1

fj 1)2 fj 2)

2

( fj

1

fj 1)2 4 fj

1

(3.3) fj 2) .

2

(3 fj

(3.4)

In smooth regions. Taylor expansion of (3.2)–(3.4) gives, respectively, IS0 IS1 IS2 where f ( f h2)2 ( f h2)2 (f h )

2 2

(2 f h (2 f h (2 f h 0, then

f h3)2 f h3)2 f h)

3 2

O(h6) O(h6) O(h ),

6

(3.5) (3.6) (3.7)

f (uj). If f ISk ( f h)2(1

O(h2)), k

0, 1, 2,

(3.8)

which means the weights (see (2.15)) resulting from this measurement satisfy (2.17) for r 3; thus we obtain a ?fth-order (the optimal order for r 3) accurate WENO scheme. Moreover, this measurement is also more accurate at critical points of f (u(x)). When f 0, we have ISk ( f h2)2(1 O(h2)), k 0, 1, 2, (3.9)

In this section, we present a new way of measuring the smoothness of the numerical solution on a stencil which can be used to replace (2.18) to form a new weight. As we 1)thknow, on each stencil Sk , we can construct a (r order interpolation polynomial, which if evaluated at x xj 1/2 , renders the approximation of hj 1/2 given in (2.7). Since total variation is a good measurement for smoothness, it would be desirable to minimize the total variation for the approximation. Consideration for a smooth ?ux and for the role of higher order variations leads us to the following measurement for smoothness: let the interpolation polynomial on stencil Sk be qk(x); we de?ne

r 1

which implies that the weights resulting from the measurement (3.1) is also ?fth-order accurate at critical points. To illustrate the different behavior of the two measurements (i.e., (2.18) and (3.1)) for r 3 in smooth monotone regions, near critical points or near discontinuities, we compute the weights 0 , 1 , and 2 for the function sin 2 xj 1 if 0 xj xj 0.5, 1.

fj

sin 2 xj if 0.5

(3.10)

ISk

l 1

xj 1/2 xj 1/2

h

2l 1

(l) (qk )2

dx,

(3.1)

(l) where qk is the lth-derivative of qk(x). The right-hand side of (3.1) is just a sum of the L2 norms of all the derivatives of the interpolation polynomial qk(x) over the interval (xj 1/2 , xj 1/2). The term h2l 1 is to remove h-dependent factors in the derivatives of the polynomials. This is similar to, but smoother than, the total variation measurement based on the L1 norm. It also renders a more accurate WENO scheme for the case r 3, when used with (2.15) and (2.16). When r 2, (3.1) gives the same measurement as (2.18). However, they become different for r 3. For r 3, (3.1) gives

IS0

( fj

2

2 fj

1

fj)2

( fj

2

4 fj

1

3 fj)2

(3.2)

at all half grid points xj 1/2 , where xj j x, xj 1/2 xj x/2, and x . We display the weights 0 and 1 in Fig. 1 ( 2 1 0 1 is omitted in the picture). Note the optimal weight for 0 is C 3 0.1 and for 1 is C 3 0 1 0.6. We can see that the weights computed with (2.18) (referred to as the original measurement in Fig. 1) are far less optimal than those with the new measurement, especially around the critical points x , . However, near the discontinuity x , the two measurements behave similarly; the discontinuous stencil always gets an almostzero weight. Moreover, for the grid point immediately left to the discontinuity, 0 and 1 , which means, when only one of the three stencils is non-smooth, the other two stencils get O(1) weights. Unfortunately, these weights do not approximate a fourth-order scheme at this point. A similar situation happens to the point just to the right of the discontinuity. For simplicity of notation, we use WENO-X-3 to stand

208

JIANG AND SHU

FIG. 1. A comparison of the two smoothness measurements.

for the third-order WENO scheme (i.e., r 2, for which the original and new smoothness measurement coincide) (where X LF, Roe, RF refer, respectively, to the global Lax–Friedrichs ?ux splitting, Roe’s ?ux splitting, and Roe’s ?ux splitting with entropy ?x). The accuracy of this scheme has been tested in [9]. We will use WENO-X-4 to represent the fourth-order WENO scheme of Liu et al. (i.e., r 3 with the original smoothness measurement of Liu et al.) and WENO-X-5 to stand for the ?fth-order WENO scheme resulting from the new smoothness measurement. In later sections, we will also use ENO-X-Y to denote conventional ENO schemes of ‘‘Y’’th order with ‘‘X’’ ?ux splitting. We caution the reader that the orders here are in L1 sense. So ENO-RF-4 in our notation refers to the same scheme as ENO-RF-3 in [15]. In the following we test the accuracy of WENO schemes on the linear equation: ut ux 0, 1 x 1, (3.11) (3.12)

u0(x) sin4( x). Again we see that WENO-RF-4 is more accurate than WENO-RF-5 on the coarsest grid (N 20) but becomes less accurate than WENO-RF-5 on ?ner grids. The order of accuracy for WENO settles down later than in the previous example. Notice that this is the example for which ENO schemes lose their accuracy [12].

4. A SIMPLE WAY FOR COMPUTING WEIGHTS FOR EULER SYSTEMS

For system (1.1) with d 1, the derivatives dfi /dxi , i 1, ..., d, are approximated dimension by dimension; for

TABLE III Accuracy on ut

Method WENO-RF-4 N 10 20 40 80 160 320 10 20 40 80 160 320 10 20 40 80 160 320

ux

0 with u0(x)

L order — 2.13 2.81 3.05 3.37 3.86 — 4.36 4.99 4.95 5.07 5.03 — 4.96 4.99 5.00 5.00 5.00

sin( x)

L1 error 7.93e-3 1.32e-3 1.56e-4 1.13e-5 6.88e-7 2.74e-8 1.60e-2 7.41e-4 2.22e-5 6.91e-7 2.17e-8 6.79e-10 3.07e-3 9.92e-5 3.14e-6 9.90e-8 3.11e-9 9.73e-11 L1 order — 2.59 3.08 3.79 4.04 4.65 — 4.43 5.06 5.01 4.99 5.00 — 4.95 4.98 4.99 4.99 5.00

L error 1.31e-2 3.00e-3 4.27e-4 5.17e-5 4.99e-6 3.44e-7 2.98e-2 1.45e-3 4.58e-5 1.48e-6 4.41e-8 1.35e-9 4.98e-3 1.60e-4 5.03e-6 1.57e-7 4.91e-9 1.53e-10

u(x, 0)

u0(x), periodic.

In Table III, we show the errors of the two schemes at t 1 for the initial condition u0(x) sin( x) and compare them with the errors of the ?fth-order upstream central scheme (referred to as CENTRAL-5 in the following tables). We can see that WENO-RF-4 is more accurate than WENO-RF-5 on the coarsest grid (N 10) but becomes less accurate than WENO-RF-5 on the ?ner grids. Moreover, WENO-RF-5 gives the expected order of accuracy starting at about 40 grid points. In this example and the one for Table IV, we have adjusted the time step to t ( x)5/4 so that the fourth-order Runge–Kutta in time is effectively ?fth-order. In Table IV, we show errors for the initial condition

WENO-RF-5

CENTRAL-5

WEIGHTED ENO SCHEMES

209 1, are the weights in the sth

TABLE IV Accuracy on ut

Method WENO-RF-4 N 20 40 80 160 320 640 20 40 80 160 320 640 20 40 80 160 320 640

ux

0 with u0(x)

L order — 1.56 2.43 3.68 4.08 3.81 — 3.60 2.31 3.88 4.80 5.48 — 4.40 4.89 4.97 5.00 5.00

sin ( x)

L1 error 3.29e-2 9.99e-3 1.44e-3 8.31e-5 3.06e-6 9.57e-8 4.91e-2 3.64e-3 5.00e-4 2.17e-5 6.17e-7 1.57e-8 3.35e-2 1.52e-3 5.09e-5 1.60e-6 4.99e-8 1.56e-9 L1 order — 1.72 2.79 4.12 4.76 5.00 — 3.75 2.86 4.53 5.14 5.30 — 4.46 4.90 4.99 5.00 5.00

4

Here k,s , k 0, 1, ..., r characteristic ?eld,

k,s k(ls

L error 7.31e-2 2.48e-2 4.60e-3 3.59e-4 2.12e-5 1.51e-6 1.08e-1 8.90e-3 1.80e-3 1.22e-4 4.37e-6 9.79e-8 5.23e-2 2.47e-3 8.32e-5 2.65e-6 8.31e-8 2.60e-9

fj

r 1,

..., ls fj

r 1)

(4.3)

which is a nonlinear function ( k is de?ned by (2.15)). The numerical ?uxes obtained in each characteristic ?eld can then be projected back to the component space by ? fj

m 1/2 s 1

WENO-RF-5

?j f

1/2,srs .

(4.4)

CENTRAL-5

example, when approximating df1 /dx1 , one ?xes xl , l 1, and uses a one-dimensional approximation in the direction of x1 . In the following, we only discuss how to approximate df1 /dx1 and will drop the index ‘‘1’’ for simplicity. We will also assume that all the eigenvalues of the Jacobian df/du are nonnegative (a condition identical to f 0 in the scalar equation). For general ?ux, one can split it locally into positive and negative parts just as in the scalar case. The formulas for the negative part of the ?ux will be omitted due to symmetry. ? For systems of equations, the ?ux fj 1/2 are usually approximated in the (local) characteristic ?elds. Let us take Aj 1/2 to be some average Jacobian at xj 1/2 , e.g., the arithmetic mean Aj f u (4.1)

u (uj uj 1)/2

1/2

or for Euler systems, the Roe’s mean matrix [11]. We denote by rs (column vector) and ls (row vector) the sth right and left eigenvectors of Aj 1/2 , respectively. Then the scalar WENO scheme can be applied to each of the characteristic ?elds. For example, (2.10) becomes

Because of the nonlinearity of the weights (see (4.3)), the above procedure involves many local projections (or vector vector products). In fact, these projections are responsible for most of the ?oating point operations of WENO schemes (true also for ENO schemes). Moreover, these projections cannot be avoided if the weights are to be computed from the projected quantities. However, if the weights can be computed from other quantities, we then can exploit the linearity of the rest of the scheme (e.g., the linearity of qr ) to reduce the number of ?oating point operations bek cause the only nonlinear part of WENO schemes is in the calculation of the weights. The question then is what quantities can serve as replacements of the projected values. Obviously for each characteristic ?eld, the replacing quantity must indicate the jump discontinuities in that ?eld. Although such quantities are yet to be discovered for general systems of equations, we ?nd, after an extensive searching and trial, that pressure and entropy are good replacements for the projected values when Euler systems are concerned, at least for problems without strong shocks and re?ective waves. Namely, we will use pressure to compute the weights in the genuinely nonlinear characteristic ?elds (s 1, m) and use entropy for the linearly degenerate ?eld(s) (1 s m). The motivation: (1) The pressure does not jump at contact discontinuities but always jumps at shocks; (2) The entropy jumps at contact discontinuities but jumps only slightly at a weak shock. Since the pressure and entropy can be obtained independent of the characteristic projection procedure, we can reformulate the WENO schemes to take advantage of the linearity of the rest of the scheme. Let us de?ne

r 1 r k,sqk(fj k r 1 , k 0

Fj

1/2,s

..., fj k), s

1, ..., m.

(4.5)

?j f

r 1 1/2,s k 0 r k,sqk(ls

fj

k r 1,

..., ls fj k)

(4.2)

which gives the numerical ?ux in the sth characteristic ?eld.

For Euler systems, the sth (1 s m) characteristic ?eld is linearly degenerate. These ?elds have the same characteristic speed (eigenvalue) and the weights are all computed from the entropy. So we have for all 1 s m,

210

k,s k,2

JIANG AND SHU

k

0, ..., r

1

m(a1 , ..., an) s min ai , if s

1 i n

and, therefore, Fj

1/2,s

sign(a1)

sign(an) (5.3)

Fj

1/2,2

0, and

j

otherwise;

is given by fj

j 1

for all l s m. Combine (4.2) and (4.4) and use the linearity of qr to take out ls ; we get k fj ? fj

m 1/2 s 1

2fj fj fj

fj

1

2

1

fj

,

1

(5.4)

(ls Fj (l1 ( Fj

1/2,s)rs

1/2,1

Fj

1/2,2))r1

(4.6) Fj

1/2,2 .

(lm ( Fj

1/2,m

Fj

1/2,2))rm

As we can see, we only need two projections from component space to characteristic space and two inverse projections, plus the few operations for computing Fj 1/2,s , s 1, 2, m. We will denote, by WENO-LF-5-PS, the WENO scheme for the case r 3, which uses pressure and entropy for weight computation in conjunction with the new smoothness measurement (3.1), the weights (2.15), and global Lax–Friedrichs ?ux splitting (according to our numerical tests, the original smoothness measurement of Liu et al. does not perform well at shocks when combined with the above way of computing weights). The accuracy of WENO-LF-4, WENO-LF-5, and WENO-LF-5-PS on the 1D Euler system is tested using an initial condition which produces a smooth solution, the same example used in Section 6. The result is similar to the scalar case in Table III and thus will not be shown.

5. SHARPENING OF CONTACT DISCONTINUITIES

where is a positive parameter. We will use 33 as suggested by Yang [20] in all our tests in Section 8, although this parameter can be tuned to optimize the results for individual problems. The case of a 0 can be treated symmetrically and the generalization to variable coef?cient or nonlinear problems is rather straightforward. See [15] for details. We would like to point out that, in smooth regions, Yang’s ACM adds to the regularly computed ?uxes a term which is of the same order of the truncation error size; thus no accuracy loss is expected. We will apply the above sharpening technique only to contact discontinuities or contact characteristic ?led(s) in case of Euler systems. A scheme which uses the above arti?cial technique will be denoted by adding to its name the suf?x ‘‘-A’’, e.g., WENO-LF-5-A.

6. CONVERGENCE FOR SMOOTH SOLUTIONS

As we can see from the previous sections, the WENO schemes are smooth in the sense that the spatial operator L, L L( fj r , fj

r 1,

..., fj

r 1),

(6.1)

For a linear, constant coef?cient problem ( f (u) au in (2.1)), Yang’s arti?cial compression method, when applied to the WENO schemes is simply (assuming a 0)

? jA 1/2 f

?j f

1/2

cj

1/2 ,

(5.1)

? where fj 1/2 is the ?ux obtained by one of the methods introduced in the previous three sections and

cj m

j

1/2

2

?j m( f R 1/2 ? R 1/2 fj

? fj

fj

1/2 ,

? R 1/2 fj

?j f

1/2), fj 1

is in?nitely differentiable to any of its arguments (see (2.2), (2.10), (2.15), (2.16), and (2.18) or (3.1)). Here r 2 is the L1 order of the base ENO scheme. In case of a general ?ux, if a smooth ?ux splitting is used (e.g., the global Lax– Friedrichs ?ux splitting), the smoothness of the WENO schemes is unchanged. Strang’s theorem (Theorem I in [18]) implies that, for a conservation law whose ?ux function and solution have enough continuous derivatives, a smooth, consistent scheme is convergent if its ?rst variation (see [18] for the de?nition) is l2-stable. It is easy to see that, for the scalar one-dimensional conservation law (2.1) with f 0, the spatial operator of ? WENO schemes has the following simple ?rst variation L,

(5.2)

1

? L

j r 1 l j r

?j f

1/2 ,

.

L (uj , ..., uj )ul ul

r 1 1,

? Here ? R 1/2 is obtained by the same method for fj 1/2 , prefj tending a 0; m is the usual minmod function de?ned by

f (uj ) 2r (q r 1(uj x

..., uj

r 1)

(6.2)

q 2r 11(uj r , ..., uj r

r 2))

WEIGHTED ENO SCHEMES

211

TABLE V CFL Numbers

n r r 2 3 3 n 4

because ( k / ul )(uj , ..., uj ) 0 and k(uj , ..., uj ) Cr k for all k 0, ..., r 1 and l j r, ..., j r 1; (6.2) can be rewritten into a summation of a (2r 2)th-order central difference D2r 2 and a (2r 1)th-order upwind biased difference,

? L

f (uj ) 2r 2 (D (uj x ( 1)r

1 r 2r 1

1.625 1.434

1.745 1.731

r 1,

..., uj

r 1)

(6.3)

Note. n

order of the Runge–Kutta scheme.

uj r ), 1D, 2D, and 3D Euler systems are solved. The 3D Euler system is (1.1) with d 3, m 5, and u f(u) g(u) h(u) ( , u, v, w, E)T; ( u, P u 2, uv, uw, u(E v , vw, v(E w 2, w(E

2

where r (r 1)!(r 1)!/(2r 1)! 0 and 2r 1 is the (2r 1)th-order forward difference operator. Applying the classical Fourier analysis to the ?rst variation, we see that the (2r 2)th-order central difference renders purely imaginary spectrum while the second term in (6.3), which is just a (2r 1)th-order upwind biased difference, has a spectrum of the form

2r 1

(7.1) P))T; P)) ; P))T,

T

(7.2) (7.3) (7.4)

( v, vu, P

( w, wu, wv, P

2

2r 1

r

sin

2

sin

2

i cos

2

(6.4) where P ( 1)(E (u 2 v2 w 2)).

2 ; (6.3) and (6.4) together imply that the where 0 ? spectrum of the operator L lies fully on the left half of the complex plane. Therefore, with appropriately chosen CFL number, the ?rst variation of the WENO schemes are l2-stable when the third or higher order Runge–Kutta time discretization is used. Let us de?ne by u(x0 t0 , x) the numerical solution at (x0 , t0) Rd R for grid size x and ?xed CFL number. For general scalar conservation laws, the same analysis gives the following. THEOREM 6.1. For the initial value problem of (1.1) with m 1 (i.e., scalar conservation laws), (x0 , t0) Rd R , if the exact solution v and df/dv, g have r [(d 1)/2] q0 2 continuous derivatives in the domain of dependence of (x0 , t0) as de?ned in [18], the WENO schemes using a smooth ?ux splitting and a nth-order Runge–Kutta scheme (n max(r, 3)) satisfy u(x0 , t0 , x) v(x0 , t0) O( x r ) (6.5)

The initial condition is 1 0.2 sin( (x y z)), u v w 1, P 1.

for an appropriately chosen CFL number. Here q0 is a small constant integer (see [18]). For a few special cases, we list the CFL numbers in Table V.

7. EFFICIENCY COMPARISON

In this section, we compare the ef?ciency of WENOLF-4, WENO-LF-5, WENO-LF-5-PS, and ENO-LF-4 on a vector supercomputer (CRAY C-90) and two serial workstations (SUN Sparc10 and SGI Indigo2).

Here we use f, g, h, x, y, z, instead of f1 , f2 , f3 , x1 , x2 , x3 . The 1D and 2D Euler systems and their initial conditions can be deduced from the above 3D problem by removing the extra degree of freedom(s). We display the CPU time of ENO-LF-4, WENO-LF4, WENO-LF-5, and WENO-LF-5-PS (all with RK-4) on CRAY C-90, Sparc10, and SGI Indigo2 in Table VI. We observe that WENO-LF-4 and WENO-LF-5 are at least twice as fast as ENO-LF-4 on CRAY C-90 and WENOLF-5-PS is 2.5 times as fast as ENO-LF-4 for 1D Euler problems, 3.2 times as fast for 2D Euler problems, and 3.9 times as fast for 3D Euler systems. On the workstations, WENO-LF-4 and WENO-LF-5 are a bit faster than ENOLF-4 on SUN Sparc10 but a bit slower on SGI Indigo2. WENO-LF-5-PS is 1.5 to 2.2 times as fast as ENO-LF-4 on SUN Sparc10 and on SGI Indigo2. As a reference, we also include the CPU times of a typical second-order scheme [8] (the positive scheme, Van Leer’s limiter with second-order Runge–Kutta scheme in time, our own implementation) in the following tables. We can see the second-order scheme is about 10 times as fast as ENO-LF-4 on CRAY C-90, 4.5 times as fast on SUN Sparc10 and 3.5 times as fast on SGI Indigo2. In Table VII, the number of ?oating point operations and the MFlops (million ?oating-point operations per sec-

212

JIANG AND SHU

TABLE VI CPU Time in Seconds

d 1 2 3 N 1600 200 60 2nd-order 1.75 13.13 15.48 ENO-LF-4 16.67 122.52 171.42 WENO-LF-4 7.44 63.93 76.79 WENO-LF-5 7.45 60.84 78.89 WENO-LF-5-PS 6.29 37.67 43.47

1 2 3

1600 200 40

69.43 512.33 178.95

SUN Sparc 10 (66 MHz, HyperSparc), compiled with ‘‘-r8 -O4’’ 311.22 317.55 319.02 2582.25 2132.50 2116.53 807.75 716.05 754.88 SGI Indigo2 (75 MHz R8000), compiled with ‘‘-r8 -O3’’ 66.21 73.88 555.51 564.54 626.92 699.58

d

215.95 1163.72 389.77

1 2 3

1600 200 60

21.03 151.26 167.44

77.01 578.22 715.91

58.14 347.48 366.29

Note. N points in each spatial dimension; 104

iterations for the d-dimensional system.

ond) are given for the second-order scheme, ENO-Roe-4, ENO-LF-4, WENO-LF-5, WENO-LF-5-PS, ENO-Roe-4A, and WENO-LF-5-A. The operation count and MFlops for WENO-LF-4 is about the same as those for WENOLF-5, thus omitted in the table. We can see all the WENO schemes achieve the speed of about 500 MFlops, which is 50% of the peak speed of CRAY C-90. The decrease of MFlops for high dimensions is because of the shorter array length N used in our tests. Notice also that the operation count per grid point per Runge–Kutta stage, of the full

characteristic based fourth- or ?fth-order ENO schemes using Lax–Friedrichs building blocks, is about 3 to 4 times that of the second-order schemes. This ratio actually decreases to only about 1.5 if the Roe building block is used instead; i.e., f (u) and f (u) are not approximated separately. This is somewhat surprising, as it was commonly believed that high order methods are much more expensive than lower order ones. When Roe building block is used, Yang’s arti?cial compression causes 40% increase in operation count for 1D Euler systems and 65% increase for 2D

TABLE VII Number of Operations per Runge–Kutta Stage per Grid Point and MFlops on CRAY C-90

Scheme 2nd-order d 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 1 2 x y x y x/y 9 22 39 3 6 9 3 6 9 27 70 129 13 26 39 4 10 37 110 x 8 20 36 19 50 93 39 102 189 3 6 9 3 6 9 33 106 19 70 sign(x) 3 8 15 0 0 0 0 0 0 0 0 0 0 0 0 6 24 12 48 x y, 3 6 9 3 6 9 3 6 9 3 6 9 4 8 12 3 6 3 6 x MFlops 478 400 350 179 191 — 223 219 190 557 503 442 474 453 357 164 178 526 482

ENO-Roe-4

ENO-LF-4

WENO-LF-5

WENO-LF-5-PS

ENO-Roe-4-A WENO-LF-5-A

82 239 476 102 309 663 244 791 1751 235 703 1466 145 341 576 144 511 375 1379

83 248 506 98 304 656 233 766 1718 284 838 1718 129 315 579 135 484 447 1654

Note. d

the spatial dimension.

WEIGHTED ENO SCHEMES

213

Euler systems as we can see from the operation counts for ENO-Roe-4 and ENO-Roe-4-A. When a Lax–Friedrichs building block is used, the increase of operation count is 65% for 1D Euler systems and 100% for 2D Euler systems as shown by the operation counts for WENO-LF-5 and WENO-LF-5-A. The increase in CPU time is well re?ected by the above percentages. We note that ef?ciency comparison is a highly circumstantial practice. It depends on the hardware, the compiler options, and the coding of the speci?c schemes, as well as on the choice of speci?c scheme features, such as the time discretization. We have tried to detail all the circumstances of our tests as much as possible.

8. NUMERICAL RESULTS

best at all waves. Note, we have adjusted the CFL number for WENO-Roe-5-A from 0.4 to 0.2. For CFL 0.4, using 40 in (5.4) gives similar results. EXAMPLE 2 (Non-convex problems). We test the thirdand ?fth-order WENO schemes on the Buckley–Leverett problem whose ?ux is f (u) 4u 2 (1 (8.1)

4u

2

u)2

8.1. Scalar Conservation Laws in One Dimension EXAMPLE 1 (Linear equation). We solve the linear equation: ut ux 0, 1 x 1,

u(x, 0) where

uu0(x), periodic,

with initial data u 1 in [ , 0] and u 0 elsewhere. (For the numerical results of the fourth-order WENO scheme of Liu et al., see [9]). The exact solution is a shock– rarefaction–contact discontinuity mixture. The results obtained by WENO-RF-3 (with RK-3) and WENO-RF-5 (with RK-4) are shown in Fig. 3. We can see both schemes converge to the correct entropy solution and give sharp shock pro?le. Note that, around discontinuities, WENO schemes are simulating the base ENO schemes. Therefore the sharpness of the shock pro?le obtained by the WENO schemes are only expected to be as good as that obtained by the base ENO schemes. However, in terms of this sharpness, our tests show that the third-order WENO scheme is comparable to the third-order ENO

(G(x, , z 1, u0(x) 1 10(x

)

G(x, , z

)

4G(x, , z)),

0.8 0.4

x x 0.2; x 0.6;

0.6; 0.2;

0.1) , ) F(x, , a ) 4F(x, , a)),

0 0.4

x

(F(x, , a 0,

otherwise.

G(x, , z) F(x, , a)

e

(x z)2

max(1

2

(x

a)2, 0).

schemes, instead of the base second-order ENO scheme, and the ?fth-order WENO scheme is comparable to the fourth-order ENOscheme. 8.2. Euler System in One Dimension EXAMPLE 1 (1D Riemann problems). We consider here two well-known problems which have the following Riemann type initial conditions: uL if x uR if x 0 0.

The constants are taken as a 0.5, z 0.7, 0.005, 10, and log 2/36 2. The solution contains a smooth but narrow combination of Gaussians, a square wave, a sharp triangle wave, and a half ellipse. We compute the solution up to t 8 with 200 points. The results are shown in Fig. 2. We observe that both WENO-Roe-4 and WENO-Roe-5 perform better than ENO-Roe-4 for all the four types of waves in the initial condition. WENO-Roe-4 does better than WENO-Roe-5 at acute turns in the solution curve (or spike-like peaks), but WENO-Roe-5 does better for the square wave and at obtuse turns in the solution curve. With Yang’s arti?cial compression technique, WENO-Roe-5-A performs the

u(x, 0)

The ?rst one is Sod’s problem [17]. The initial data are (

L,

qL , PL)

(1, 0, 1), (

R,

qR , PR)

(0.125, 0, 0.1).

214

JIANG AND SHU

FIG. 2. Linear equation; third-order Runge–Kutta in time; 200 points, CFL 4; (c) WENO-Roe-5; (d) WENO-Roe-5-A.

0.4 (0.2 for (d) only), T

8: (a) ENO-Roe-4; (b) WENO-Roe-

The second one is the Riemann problem proposed by Lax [7], ( (

L, R,

qL , PL) qR , PR)

(0.445, 0.698, 3.528), (0.5, 0, 0.571).

The numerical results are presented in Fig. 4. We can

see that all schemes give a correct solution with good resolution. WENO-RF-5 is better than WENO-LF-5-PS which is, in turn, better than WENO-RF-3. We note that Figs. 4b and 4c (Figs. 4e and 4f) are comparable, respectively, to Fig. 10 (Fig. 11) in [15]. Also see Fig. 9a (Fig. 10a) in [9]. WENO-LF-5-A does much better than all other schemes at the contact discontinuities. We would like to point out

WEIGHTED ENO SCHEMES

215

FIG. 3. The Buckley–Leverett problem: (a) WENO-RF-3; (b) WENO-RF-5.

that, according to our experience with extensive numerical testing, these two problems, especially the Lax’s problem, are tough test cases for non-characteristic-based schemes of order at least three. Oscillations can easily appear for such schemes. Here WENO-LF-5-PS performs well in these two cases. To save space, the ?gures for velocity and pressure are not shown. EXAMPLE 2 (1D shock entropy wave interaction). In this example, we test the WENO schemes on a model that involves a moving shock interacting with an entropy wave of small amplitude. On a domain [0, 5], the initial condition is 3.85714; e

sin (kx)

u ; u

2.629369; P 0; P

10.33333; when x 1; when x

0.5 0.5,

where and k are the amplitude and wave number of the entropy wave, respectively. The mean ?ow is a pure rightmoving Mach 3 shock. If is small compared to the shock strength, the shock will march to the right at approximately the non-perturbed shock speed and generate a sound wave which travels along with the ?ow behind the shock. At the same time, the perturbing entropy wave, after ‘‘going through’’ the shock, is compressed and ampli?ed and travels approximately at the speed of u c, where u and c are the velocity and speed of the sound of the mean ?ow left to the shock. The ampli?cation factor for the entropy wave can be obtained by linear analysis. See [10, 21] for details. In order to get rid of the transient waves due to the non-numerical initial shock pro?le, we let the shock move up to x 4.5 and then shuf?e it back to x 0.5. The solution is examined when the shock reaches x 4.5 the second time. The goal of this test is to examine the stability and accuracy of the WENO schemes at the presence of the shock. Since the entropy wave here is set to be very weak relative to the shock, any excessive oscillation could pollute

the generated waves (e.g., the sound waves) and the ampli?ed entropy waves. In our tests, we take 0.01 and k 13. The amplitude of the ampli?ed entropy waves predicted by the linear analysis is 0.08690716 (shown in the following ?gures as horizontal solid lines). First we use 800 points which is effectively 20 points in each wave length of generated entropy wave. Since the generated sound waves (or pressure wave) are of lower frequency than the ampli?ed entropy waves, they are much better resolved by this grid size. Therefore we only display the entropy component of the numerical solutions. WENO-LF-4, WENO-LF-5, and WENO-LF-5-PS are used in our tests and the results are shown in Fig. 5 (the mean ?ow has been subtracted from the numerical solution). We see that all three schemes catch the ampli?ed entropy waves quite well. WENO-LF-4 performs the best on this grid and WENO-LF-5 ranks second. In order to examine the relative performance of WENO-LF-4 and WENO-LF-5, we run the same test on a grid of 1200 points. The results for these two schemes are displayed in Fig. 6. We can see that on this grid (approximately 30 points per wave length), WENO-LF-5 is as accurate as WENO-LF-4. In fact, on ?ner grids, WENO-LF-5 becomes more accurate than WENO-LF-4. This is in good agreement with our accuracy test in Section 4. For the purpose of comparison with low order schemes, we also include the entropy computed by a typical second-order scheme [8] (half Van Leer’s limiter, half Superbee limiter with second-order Runge–Kutta scheme in time, 2000 points). The advantage of using higher order schemes for this example is apparent. EXAMPLE 3. (Two interacting blast waves). We consider here the interaction of two blast waves. The initial data are uL u(x,0) if 0 x x x 0.1 0.9 1,

uM if 0.1 uR if 0.9

216

JIANG AND SHU

FIG. 4. Density: (a)–(d), Sod’s problem; (e)–(h), Lax’s problem.

where

L M R

1, uL 10 2, PR

uM

uR

0,

PL

103, PM

102.

The re?ective boundary condition is applied at both x 0 and x 1. See [19] for a detailed discussion of this problem. Three grids are used: 199, 399, 799 points. We examine our numerical solutions at t 0.038. The ‘‘exact’’ solution (solid lines in all the pictures) are computed by ENO-RF-

WEIGHTED ENO SCHEMES

217

FIG. 5. 1D shock entropy wave interaction: entropy,

t

0.6

x.

218

JIANG AND SHU

FIG. 6. 1D shock entropy wave interaction (continued): Entropy,

t

0.6

x.

5 with 1600 points. In Fig. 7, we show the density computed by WENO-RF-3 (with RK-3), WENO-RF-4, WENO-RF5, and WENO-RF-5-A (with RK-4). We observe that the fourth-order and ?fth-order WENO schemes are much better than the third-order WENO scheme and the results are comparable with those obtained by the unmodi?ed ENO-RF-4 (see Fig. 12 in [15]. Note, the fourth-order ENO scheme in L1 norm was denoted as ENO-RF-3 there). WENO-RF-4 is slightly better than WENO-RF-5 on the medium grid while on the ?ne grid WENO-RF-5 seems to be better. The results of WENORF-5-A on the coarse and medium grids are nearly as good as WENO-RF-5 on, respectively, medium and ?ne grids. EXAMPLE 4 (Quasi-one-dimensional nozzle ?ow). In this example, we use the WENO schemes to solve the steady state quasi-1D nozzle ?ow. The governing equation of the quasi-1D nozzle ?ow is the 1D Euler system with the forcing term g(u, x) Ax ( u, u2, u(E A P))T,

and 1.8 at x 1 (the exit). The exit ?ow condition is then decided by the prescribed shock position, which is x 0.5 in our test. In Fig. 8, we display the density computed by WENORoe-4 and WENO-Roe-5 with 34 points. We can see both schemes converge nicely to the exact solution (solid line in the pictures). The residue computed with both schemes settles down to 10 7 for this grid, and to a smaller number for a ?ner grid. This example shows that WENO has its advantage for steady state calculations. 8.3. Euler System in Two Dimensions EXAMPLE 1. (Oblique Sod’s problem). The purpose of this test is to analyze the capability of WENO schemes in resolving waves that are oblique to the computational mesh. The 2D Sod’s problem is solved where the initial /2). If jump makes an angle against x-axis (0 /2, we have the one-dimensional Sod’s problem. If 0 /2, all the waves produced will be oblique to the rectangular computational mesh. We take our computational domain to be [0, 6] [0, 1] and position the initial jump at (x, y) (2.25, 0). The physical domain varies with and is taken as [0, 6/sin ] [0, 1/sin ]. The scaling factor 1/sin is to ensure the same grid resolution normal to the wave propagation on a given mesh at some ?xed

where A A(x) is the cross area function of the nozzle and Ax dA/dx. The nozzle here is unit length long, whose shape is determined by assuming a linear, isentropic Mach number distribution, which is 0.8 at x 0 (the entrance)

WEIGHTED ENO SCHEMES

219

FIG. 7. Two interacting blast waves.

time for all choices of . See [3] for details. We take to be arc tan 1, arc tan 2, and arc tan 4. The solution is computed up to t 1.2 on a 96 16 mesh. WENO-LF-4, WENO-LF-5, WENO-LF-5-A, and WENO-LF-5-PS are used in our tests and the results are

compared with that obtained by ENO-LF-4 (all with RK4 and t 0.6 x). For the case arc tan 1, we display the density contours obtained by WENO-LF-5-PS in Fig. 9a; in Fig. 9b, we show the densities at y 0 obtained by all four schemes. We can see that all WENO schemes

220

JIANG AND SHU

FIG. 8. Density: Steady quasi-1D nozzle ?ow; 34 points, RK-3 in time.

are doing well in resolving the oblique waves and their differences from the ENO-LF-4 (except WENO-LF-5-A) are barely noticeable. WENO-LF-5-A gives sharp pro?le of the contact discontinuity as expected. In Figs. 9c–f a more quantitative study is carried out. Namely, for each scheme, we measure the differences between oblique cases and the one-dimensional case. We can see that WENOLF-4 and WENO-LF-5 perform similarly as ENO-LF-4 does while WENO-LF-5-PS gives slightly larger deviation near the contact discontinuity. However, this small difference can be regarded as negligible. We want to note that WENO-LF-5-PS performs well at the shock. EXAMPLE 2 (A Mach 3 wind tunnel with a step). This model problem has been carefully examined in [19]. The setup of the problem is the following: The wind tunnel is 1 length unit wide and 3 length units long. The step is 0.2 length units high and is located 0.6 length units from the left-hand end of the tunnel. The problem is initialized by a right-going Mach 3 ?ow. Re?ective boundary conditions are applied along the walls of the tunnel and in-?ow and out-?ow boundary conditions are applied at the entrance (left-hand end) and the exit (right-hand end). For the treatment of the singularity at the corner of the step, we adopt the same technique used in [19], which is based on the assumption of a nearly steady ?ow in the region near the corner. WENO-LF-4 and WENO-LF-5 are used in our tests and the results are compared with those obtained by ENO-LF4 (all with RK-4 and t 0.6 x). Two grids are used: 122 39 and 242 79. They correspond respectively to the medium and ?ne grids in [19]. In Figs. 10 to 11 we show the density component obtained by all three schemes on the two grids. We can see that all the schemes perform well with good resolution. Relatively speaking, WENO-LF-4 and WENO-LF-5 have slightly better resolution at the contact line (originated from the Mach step) and contain less visible ‘‘bumps,’’ which are indeed small numerical oscillations, than ENOLF-4.

WENO-LF-5-PS does not work for this problem because of the strong re?ecting waves. EXAMPLE 3 (Double Mach re?ection of a strong shock). The computational domain for this problem is chosen to be [0, 4] [0, 1]. The re?ecting wall lies at the bottom of the computational domain starting from x . Initially a right-moving Mach 10 shock is positioned at x ,y 0, and makes a 60 angle with the x-axis. For the bottom boundary, the exact postshock condition is imposed for the part from x 0 to x and a re?ective boundary condition is used for the rest. At the top boundary of our computational domain, the ?ow values are set to describe the exact motion of the Mach 10 shock. See [19] for a detailed description of this problem. Two grids have been used in our tests: 240 59 and 480 119. They correspond to the medium and ?ne grids in [19], respectively. We will only show the solutions on part of the domain: [0, 3] [0, 1], where most of the ?ow features are located. We use WENO-LF-4, WENO-LF-5, and ENO-LF-4 (all with RK-3 and t 0.6 x) in our tests. We show the density contours obtained by these three schemes; see Figs. 12 and 13. We see that all three schemes resolve the two Mach stems well. Again WENO-LF-5-PS does not work because of the strong re?ecting wave pattern in this problem. EXAMPLE 4 (2D shock entropy wave interaction). In this example, we test the WENO schemes on a 2D model that involves the interaction between a normal shock and a weak entropy wave which makes an angle r (0, /2) against the x-axis. If r 0, we have essentially the 1D problem (see Example 2 in Section 8.2). Since the weak entropy waves are oblique to the shock, the waves generated by the interaction are much more dif?cult to resolve than in the 1D case. Our goal here is to further examine the capability of the WENO scheme in capturing such small scale waves at the presence of shock. See [21, 16] for detailed discussions on this subject. The setup of the

WEIGHTED ENO SCHEMES

221

FIG. 9. Oblique Sod’s problem; (a) density contours, WENO-LF-5-PS, 0: (c) ENO-LF-4; (d) WENO-LF-4; (e) WENO-LF-5; (f) WENO-LF-5-PS.

arc tan 1; (b) density, y

0,

arc tan 1. (c)–(f)

1D ,

y

problem is the following: for a right moving normal shock of Mach number M, we add a small entropy wave to the ?ow on the right of the shock which is equivalent to changing only the density of the ?ow on the right of the shock to:

re

r(sin r)/Pr

,

where r and Pr are respectively the density and pressure kr (x cos r y sin of the right state of the shock, r r) and kr is the entropy wave number. In order to enforce periodic boundary conditions in the y-direction, we take the computation domain to be [0, 5] [0, 2 /kr sin r]. We initially position the normal shock at x 0.5 and allow

222

JIANG AND SHU

FIG. 10. Flow past a forward facing step. Density on medium grid, 122

39: (a) WENO-LF-4; (b) WENO-LF-5; (c) ENO-LF-4.

it to move up to x 4.5 and then shuf?e the data back to x 0.5. We extract the data at the time when the shock moves up to x 4.5 again. See [16] for a similar implementation. In our tests, we take ur , the velocity on the right of 0.01, kr the shock, to be 0 and we set M 3, r 15, r 30 . We measure the performance by comparing the amplitude of the ampli?ed entropy waves, which is computed by a Fourier analysis in the y direction for all ?xed grid values x [3.4, 4.4]. Both ENO-LF-4 and ENO-LF-5 suffer a loss of accuracy if not modi?ed (see [16]). WENO-LF-4 and WENO-LF5 work nicely without any modi?cation. The loss of accuracy of the ENO schemes can be ?xed easily by the techniques introduced in [1, 13], which effec-

tively force an upstream-centered stencil to be chosen away from the shock and a free-adapted stencil to be used near the shock. The techniques can also be adapted to enhance the performance of WENO schemes by modifying the weights as C r , if ISl k wk ,

1

?k

for any l

0, ..., r

1,

(8.2)

otherwise.

0, ..., r 1, are the where is taken as 2 and wk , k regularly computed weights. Equation (8.2) leads to the optimal weights being used for stencils away from the shock and regularly computed weights being used near the shock. We denote the modi?ed ENO-LF-4, ENO-LF-5, WENOLF-4, and WENO-LF-5 to be, respectively, ENO2-LF-4,

WEIGHTED ENO SCHEMES

223

FIG. 11. Flow past a forward facing step (continued). Density on ?ne grid, 242

79: (a) WENO-LF-4; (b) WENO-LF-5; (c) ENO-LF-4.

ENO2-LF-5, WENO2-LF-4, and WENO2-LF-5. Note that, with the modi?ed weights, WENO-LF-4 becomes ?fthorder accurate in smooth regions. In all tests, RK-4 is used with t 0.6 x. First we use 800 points in the x–direction and 20 points in the y–direction, which give approximately 20 points per entropy wave length in both directions. In Fig. 14, the amplitude of the ampli?ed entropy waves obtained by all aforementioned schemes are displayed. The solid horizontal line is the amplitude predicted by linear analysis which is 0.08744786 (see [10, 21]). We see that the modi?ed schemes generally perform better in terms of accuracy and decaying rate than the unmodi?ed schemes. As a reference, we have also included the amplitudes obtained by a typical second-order scheme [8] (the positive scheme, half Van Leer’s limiter, half Superbee limiter with a second-order

Runge–Kutta scheme in time, 800 points). We can see that high order schemes perform much better than the secondorder schemes, in terms of accuracy and decaying rate for this problem. WENO-LF-5-PS does not perform well for this problem even with the remedy above. This indicates that the pressure–entropy combination is not good enough to indicate precisely the smoothness of the numerical solution. This causes oscillations generated at shocks and thus destroys the accuracy of the scheme in resolving the waves which have ‘‘undergone’’ the interaction with the shock. Remark. We have seen that WENO-LF-5-PS does not work for the step problem and the double Mach re?ection problem because it cannot handle the re?ective boundary properly. This can be explained by the following: the usual

224

JIANG AND SHU

FIG. 12. Double Mach re?ection. Density on medium grid, 240

59: (a) WENO-LF-4; (b) WENO-LF-5; (c) ENO-LF-4.

way of imposing re?ective boundary condition1 is to reverse the normal velocity at the grid points which are symmetric with respect to that boundary while setting other ?ow quantities (density, pressure, and tangential velocity) to be the same; in particular, the pressure and entropy at each pair of symmetric grid points are identical. Therefore neither the pressure nor the entropy can indicate possible jumps in the normal velocity. This failure will result in an unstable weight distribution in the normal direction near the re?ective boundary and cause fatal errors such as density becoming negative. An immediate ‘‘?x’’ seems to be

We assume here the physical boundary is ?at, as is the case in the aforementioned problems.

1

using the normal velocity to compute the weights for one of the linearly degenerate ?elds. Unfortunately, the jump in the normal velocity is not like a contact discontinuity, which belongs solely to one of the characteristic ?elds. While this might cure the ill distribution of the weights in the ?eld in which the velocity is used, it cannot cure this ill distribution in other ?elds. However, WENO-LF-5-PS can be applied to problems where the re?ective boundary is not playing a vital role. As an example, we look at the following model problem. EXAMPLE 5 (Shock vortex interaction). This model problem describes the interaction between a stationary shock and a vortex. The computational domain is taken to be [0, 2] [0, 1]. A stationary Mach 1.1 shock is posi-

WEIGHTED ENO SCHEMES

225

FIG. 13. Double Mach re?ection (continued). Density on ?ne grid, 480

119: (a) WENO-LF-4; (b) WENO-LF-5; (c) ENO-LF-4.

tioned at x 0.5 and normal to the x-axis. Its left state is ( , u, v, P) (1, , 0, 1). A small vortex is superposed to the ?ow left to the shock and centers at (xc , yc) (0.25, 0.5). We describe the vortex as a perturbation to the velocity (u, v), temperature (T P/ ), and entropy (S ln (P/ ) of the mean ?ow and we denote it by the tilde values: u ? v ? e

(1

2

)

sin

2

(8.3) (8.4)

(1

2

e ( 0,

(1

)

cos

)

? T ? S

1) 2e2 4

(8.5) (8.6)

where r/rc and r (x xc)2 ( y yc)2. Here indicates the strength of the vortex, controls the decay rate of the vortex, and rc is the critical radius for which the vortex has the maximum strength. In our tests, we choose 0.3, rc 0.05, and 0.204. The above de?ned vortex is a steady state solution to the 2D Euler equation. We use a grid of 251 100 which is uniform in y but re?ned in x around the shock using a Roberts transform (see [2] and the references there). The upper and lower boundaries are intentionally set to be re?ective. The pressure contours obtained by WENO-LF-5-PS at t 0.05, t 0.20, and t 0.35 are shown in Figs. 15(a–c). We can see that for this problem, where the re?ective boundary is nonessential, WENO-LF-5-PS works nicely. To appreciate this further, we look at the

226

JIANG AND SHU

FIG. 14. 2D shock entropy wave interaction. Amplitude of ampli?ed entropy waves. 800 points (about 20 points per entropy wave length).

FIG. 15. 2D shock vortex interaction. Pressure, (a)–(d) WENO-LF-5-PS. Thirty contours: (a) t 0.05; (b) t (e)–(g) t 0.60. Ninety contours from 1.19 to 1.37: (e) WENO-LF-5-PS; (f) WENO-LF-5; (g) ENO-LF-4.

0.20; (c) t

0.35; (d) t

0.80;

WEIGHTED ENO SCHEMES

227

solution at t 0.8. By this time one branch of the shock bifurcations has reached the top boundary and been re?ected. The pressure contours obtained by WENOLF-5-PS at this moment are shown in Fig. 15d. We see that the re?ection is well captured. In Figs. 15(e–g), we compare the results obtained by WENO-LF-5-PS, WENO-LF-5 and ENO-LF-4. Ninety contours are drawn for the pressure component in the range of (1.19, 1.37). We see that the three methods give approximately the same resolution. A careful examination reveals that WENO schemes are slightly better in the sense that less numerical noise is generated. Between the two WENO schemes, WENO-LF-5 seems a little better for the same reason. For a qualitative comparison, see also [2]. Note that a different vortex is used there. EXAMPLE 6 (Flow past a cylinder). In this test, we use the WENO schemes to simulate the supersonic ?ow past a cylinder. In the physical space, a cylinder of unit radius is positioned at the origin on a x–y plane. The computational domain is chosen to be [0, 1] [0, 1] on – plane. The mapping between the computational domain and the physical domain is x y (Rx (Ry (Rx (Ry 1) ) cos( (2 1) ) sin( (2 1)) 1)), (8.7) (8.8)

FIG. 16. Flow past a cylinder: (a) physical grid; (b) pressure. WENOLF-5 with RK-4, 20 contours.

where we take Rx 3, Ry 6, and 5 /12. See [16] for the eigenvalues and eigenvectors of 2D Euler systems on general structured grids. A uniform mesh of 60 80 in the computational domain is used. For an illustration of the mesh in the physical space (drawing every other grid line), see Fig. 16a. The problem is initialized by a Mach 3 shock moving toward the cylinder from the left. The re?ective boundary condition is imposed at the surface of the cylinder, i.e., 1, the in?ow boundary condition is applied at 0, and the out?ow boundary condition is applied at 0, 1. The pressure contour obtained by WENO-LF-5 with RK-4 and t 0.6 x is shown in Fig. 16b. Similar results can be obtained by WENO-LF-4 and ENO-LF-4.

9. CONCLUSION

With the new smooth measurement, which is based on minimizing the L2 norm of the derivatives of the interpolation polynomials, the WENO schemes formulated from the rth-order ENO schemes can be made (2r 1)th-order accurate in smooth regions of the ?ux function (in spatial variable), at least for r 2, 3. However, at discontinuities, all WENO schemes are just rth-order accurate (r is the order of the base ENO scheme).

The fourth-order WENO scheme of Liu et al. and the ?fth-order WENO scheme resulting from the new smoothness measurement are found to be at least twice as fast as the fourth-order ENO schemes on vector supercomputers (e.g., CRAY C-90) and as fast on serial machines (therefore on parallel machines as well). Many 1D and 2D numerical tests suggest that both WENO schemes are very robust for shock calculations. The fourth-order WENO scheme of Liu et al. is slightly more accurate than ?fth-order WENO scheme on coarse grids (20 points or less per wave length) but it becomes less accurate on ?ner grids. For Euler systems, we also suggest computing the weights from pressure and entropy, instead of the projected values. The resulting WENO schemes are about twice as fast as the WENO schemes which use the projected values to compute weights and work well for problems which do not contain strong shocks or strong re?ected waves. More detailed numerical results for WENO schemes can be found in [6]. We have also adopted the arti?cial compression method of Yang [20] to enhance the performance of WENO schemes at contact discontinuities. However, the CPU cost is increased by as much as 100% when a Lax– Friedrichs building block is used. We believe the idea of the arti?cial compression method can be adapted directly into the weight de?nition to achieve the sharpening effect at a much lower expense. This will be investigated in the future.

ACKNOWLEDGMENTS

We thank Jay Casper, Andrew Godfrey, Xu-Dong Liu, Eitan Tadmor, and Huanan Yang for helpful discussions. We also want to thank Stanley Osher for valuable suggestions.

228

REFERENCES

JIANG AND SHU 10. J. McKenzie and K. Westphal, Phys. Fluids 11, 2350 (1968). 11. P. Roe, J. Comput. Phys. 27, 1 (1978). 12. A. Rogerson and E. Meiberg, J. Sci. Comput. 5, 151 (1990). 13. C-W. Shu, J. Sci. Comput. 5, 127 (1990). 14. C-W. Shu and S. Osher, J. Comput. Phys. 77, 439 (1988). 15. C-W. Shu and S. Osher, J. Comput. Phys. 83, 32 (1989). 16. C-W. Shu, T. A. Zang, G. Erlebacher, D. Whitaker, and S. Osher, Appl. Numer. Math. 9, 45 (1992). 17. G. Sod, J. Comput. Phys. 27, 1 (1978). 18. G. Strang, Numer. Math. 6, 37 (1964). 19. P. Woodward and P. Colella, J. Comput. Phys. 54, 115 (1984). 20. H. Yang, J. Comput. Phys. 89, 125 (1990). 21. T. A. Zang, M. Y. Hussaini, and D. M. Bushnell, AIAA J. 22, 13 (1984).

1. H. Atkins, AIAA Paper 91-1557, in AIAA 10th Computational Fluid Dynamics Conference, Honolulu, Hawaii, June 1991. 2. J. Casper, AIAA J. 30, 2829 (1992). 3. J. Casper, C-W Shu, and H. L. Atkins, AIAA J. 32, 1970 (1994). 4. A. Harten, J. Comput. Phys. 83, 148 (1989). 5. A. Harten, B. Engquist, S. Osher, and S. Chakravarthy, J. Comput. Phys. 71, 231 (1987). 6. G-S. Jiang, Ph.D. thesis, Division of Applied Mathematics, Brown University, 1995. 7. P. D. Lax, Commun. Pure Appl. Math. 7, 159 (1954). 8. X-D. Liu and P. D. Lax, CFD J., to appear. 9. X-D. Liu, S. Osher, and T. Chan, J. Comput. Phys. 115, 200 (1994).

赞助商链接

相关文章:

更多相关标签: