tceic.com

学霸学习网 这下你爽了

学霸学习网 这下你爽了

- 学校后勤管理效能提升浅析
- Convergence of adaptive control schemes using least-squares parameter estimates
- 1 On Sliding Window Schemes For Discrete Least-Squares Approximation By Trigonometric Polyn
- Data distribution and communication schemes for least squares problems on massively distrib
- A Second Order Learning Scheme based on Iteratively Reweighted Least Squares
- Generalization based on least squares adjustment
- Least-squares meshes
- COMPACT RECONSTRUCTION SCHEMES WITH WEIGHTED ENO LIMITING FOR HYPERBOLIC CONSERVATION LAWS
- Least-squares for second order elliptic problems
- Cholesky-based methods for sparse least squares the benefits of regularization
- A Unified Analysis of A Weighted Least Squares Method for First-Order Systems
- Weighted ENO schemes for Hamilton-Jacobi equations
- Adaptive and iterative least squares support vector regression based on quadratic Renyi entropy
- [2015][Song][TNNLS]Kernel-Based Least Squares Temporal Difference with Gradient Correction
- A Novel Feature Fusion Method Based on Partial Least Squares Regression

High-Order ENO Schemes for Unstructured Meshes Based on Least-Squares Reconstruction

Carl F. Ollivier-Gooch Mathematics and Computer Science Division Argonne National Laboratory

1 Introduction

High-order accurate schemes for conservation laws for unstructured meshes are not nearly so well advanced as such schemes for structured meshes. Consequently, little or nothing is known about the possible practical advantages of high-order discretization on unstructured meshes. This article is part of an ongoing e ort to develop high-order schemes for unstructured meshes to the point where meaningful information can be obtained about the trade-o s involved in using spatial discretizations of higher than second-order accuracy on unstructured meshes. This article describes a high-order accurate ENO reconstruction scheme, called DD-L2-ENO, for use with vertex-centered upwind ow solution algorithms on unstructured meshes. The solution of conservation equations in this context can be broken naturally into three phases: 1. Solution reconstruction, in which a polynomial approximation of the solution is obtained in each control volume. 2. Flux integration around each control volume, using an appropriate ux function and a quadrature rule with accuracy commensurate with that of the reconstruction. 3. Time evolution, which may be implicit, explicit, multigrid, or some hybrid. This article focuses primarily on solution reconstruction. A new high-order ENO reconstruction

Postdoctoral Researcher. Currently: Assistant Professor, Department of Mechanical Engineering, University of British Columbia, 2324 Main Mall, Vancouver, BC V6T 1Z4, Canada. Voice: (604) 822-1854. Fax: (604) 822-2403. E-mail: cfog@mech.ubc.ca. This work was supported in part by the Mathematical, Information, and Computational Sciences Division subprogram of the O ce of Computational and Technology Research, U.S. Department of Energy, under Contract W-31-109-Eng-38. This article is declared a work of the U.S. Government and is not subject to copyright protection in the United States.

technique for unstructured meshes is presented. The scheme is uniformly accurate for smooth functions, even near extrema. Near discontinuities, the scheme gracefully reduces the order of accuracy to control overshoots. Because the scheme is based on leastsquares reconstruction, implementation on unstructured meshes is straightforward. Finally, the present scheme has better convergence behavior than stencilsearching ENO schemes because the ENO property is obtained by the use of weights which vary smoothly with the data rather than by switching. Section 2 summarizes of existing reconstruction techniques. Section 3 discusses the new reconstruction scheme, with some attention given to implementation details. Section 4 gives examples of the capabilities of the reconstruction scheme. Section 5 discusses a number of technical points concerning construction of a ow solver compatible with the new scheme. Section 6 gives several example solutions to the Euler equations. Finally, Section 7 gives some conclusions from the present work and discusses some issues still remaining for high-order accurate solution of the Euler equations on unstructured meshes.

2 Overview of Reconstruction Techniques

The use of high-order spatial discretization on structured meshes is commonplace. The development of MUSCL schemes 1] focused on attaining high-order accuracy for smooth solutions, with a drop in accuracy near discontinuities and near extrema in the solution. More recently, essentially non-oscillatory (ENO) schemes have been developed to ensure uniformly high-order accuracy for all points with a smooth neighborhood. Early work in one dimension 2, 3, 4] demonstrated the feasibility of this reconstruction scheme, which searches for the smoothest stencil for reconstruction in each control volume. Extensions to multiple space dimensions soon fol1

2 lowed 5, 6, 7, 8]. Stencil-searching ENO schemes share the problem that small changes in the solution from time step to time step causes stencil \switching" and prevents convergence to a numerical steady-state. This problem has recently been addressed by a new family of weighted ENO (WENO) schemes 9, 10]. These schemes use a weighted sum of all possible stencils rather than searching for the smoothest possible stencil. Stencils containing non-smooth data are not excluded by theses schemes, but instead are given weights on the order of truncation error. Because the weights vary smoothly with the data, these schemes should converge well (although we are unaware of any studies of this issue). These schemes are one dimensional and are applied direction-by-direction for reconstruction on multidimensional structured meshes. Several extensions of ENO schemes to unstructured meshes have been made based on stencilsearching approaches 11, 12]. These schemes, like their structured-mesh counterparts, are guaranteed to reconstruct based on smooth data when this is possible. They have the same convergence di culties as structured stencil-searching schemes. A more common approach to reconstruction on unstructured meshes is to use least-squares reconstruction followed by some limiting procedure to eliminate overshoots 13, 14]. The reconstruction scheme described here is an extension of previous work 16] on the use of data-dependent least-squares reconstruction to produce ENO schemes. Previously, high-order reconstruction was demonstrated in one dimension, and second-order reconstruction and ow solutions were shown in two dimensions. The present work discusses high-order reconstruction and ow solution in two dimensions. The reconstruction scheme uses a xed stencil similar to that used by a typical k-exact least-squares scheme. An initial data-independent least-squares reconstruction is computed, and the smoothness of the data on the stencil is inferred from the results. Data-dependent weights are applied to the least-squares problem to virtually eliminate the in uence of non-smooth data, and a data-dependent reconstruction is computed. As in WENO schemes, the data-dependent weights are chosen to satisfy the ENO property of Liu, Osher, and Chan 9]. An important di erence between the present scheme and WENO schemes is their behavior when a smooth stencil does not exist. Where enough smooth data exists | near a single discontinuity or boundary,

Similar in intent are SLIP schemes 15], which are closely related to FCT schemes. SLIP schemes are local extremum diminishing and can in principle be extended to higher than second-order accuracy.

for example | each scheme reconstructs to the nominal order of accuracy. Where there are not enough smoothly connected neighbors, the scheme automatically reduces the order of accuracy of the reconstruction locally rather than contaminating the reconstruction. The result is a reconstruction scheme that is well behaved near multiple discontinuities while retaining high-order accuracy elsewhere.

3 Data-dependent Least-Squares Reconstruction

Consider a domain that has been tessellated; the tessellation has a characteristic length scale x, at least locally. The median dual of the tessellation de nes for each vertex vi a surrounding control volume Vi . For any function u(~ x) de ned on and its control-volume averaged values ui, the DD-L2-ENO will compute an expansion Ri(~ x?~ xi) about vi that conserves the mean; has compact support; reconstructs exactly polynomials of degree k ? (equivalently, Ri(~ x?~ xi ) ? u(~ x) = O xk+1 ); and satis es the ENO property of Liu, Osher, and Chan 9]. The remainder of this section describes this process in two dimensions. Reduction to one dimension, extension to three dimensions, and application to structured meshes are all straightforward variations on the theme.

3.1 Conservation of the Mean

Z

Vi

Conservation of the mean within a control volume requires that Ri (~ x?~ xi )dA =

Z

This can be accomplished by using zero-mean polynomials in expanding about vi , that is, by writing @u (x ? x ? x ) Ri(~ x?~ xi) = ui + @x i i i + @u @y (y ? yi ? y i ) u (x ? xi)2 ? x2i +@ @x2 2

2

Vi

u(~ x)dA:

(1)

i

i

3 @ 2 u ((x ? x )(y ? y ) ? xy ) + @x@y i i i

i

Reconstruction CV Added for second order Added for third order Added for fourth order

1 Z (x ? x )n (y ? y )m dA: (3) i i Ai Vi By inspection, the expansion of Equation 2 satis es Equation 1. As a practical matter, the integral of Equation 3 is most easily computed by using Green's theorem to convert it to a boundary integral around Vi . Z xnym = (n +11)A (x ? xi )n+1 (y ? yi )m dy (4) i @Vi Figure 1: Reconstruction Stencils for DD-L2-ENO This integral may be evaluated exactly by using a Gaussian quadrature of appropriate order along the Equivalently, when reconstructing a general function boundary of the control volume. u(~ x) from control volume averages, ? xk+1 : R ( ~ x ? ~ x ) = u ( ~ x ) + O (6) i i 3.2 Compact Support practice, this accuracy requirement means that the Compact support implies that the reconstruction Ri In modi Taylor series expansion of Ri given in Equawill only use data from a stencil fVj gi whose members tion 2 ed must be carried out through the kth derivaare both physically near ~ xi and topologically near tives. To compute derivatives, we seek to mincontrol volume i. The size of the compact stencil imize the error in these predicting mean value of the is determined by the number of required derivatives. function for control volumes inthe the fVj gi. The In practice, including additional neighbors allows lee- error for a single control volume is stencil given by way for ignoring some non-smooth data while retain1 Z R (~ ing high-order accurate reconstruction. The use of 3, (7) xi )dA ? uj : Ej;i = A i x?~ 8, and 14 neighbors for second-, third-, and fourthj Vj order accuracy (to compute 2, 5, and 9 derivatives, Because this integral can be expressed in terms of gerespectively) seems to be su cient. quantities and approximations to derivatives Stencils are determined iteratively. The initial ometric of u at x ~ i , the errors Ej;i can be used to formulate stencil consists of rst vertex neighbors. For con- a least-squares for the derivatives. The retrol volumes that need a larger stencil, second vertex mainder of this problem subsection develops this formulation neighbors are added. Additional layers of neighbors in detail. are added until a large enough stencil has been found The mean value, for a single control volume V , of j for each control volume; this process is illustrated in the reconstructed function R is i Figure 1. Nearly all interior points use rst neighZ bors for second order and add second neighbors for 1 x?~ xi )dA = ui (8) third and fourth order. Boundary points often add Aj Vj Ri(~ ! another layer of neighbors. Stencils are computed in Z 1 @u a pre-processing step and stored for later use. (x ? xi )dA ? xi + @x A j Vj i ! Z 3.3 Accuracy for Smooth Functions 1 @u + @y A (y ? yi )dA ? y i j Vj Accuracy of a reconstruction for smooth functions i ! Z can be stated in two equivalent ways. The recon2 1 1 @ u 2 2 struction can be said to be k-exact for some k if, + @x2 2A (x ? xi) dA ? 2 x i j Vj i when reconstructing P (~ x) fxm yn : m + n kg from ! Z control volume averages, 2 1 @ u + @x@y A (x ? xi)(y ? yi )dA ? xyi j Vj Ri(~ x?~ xi) P (~ x): (5) i xn y m i

, where

2 u (y ? yi )2 ? y2 i + +@ @y2 i 2

(2)

4

2 u 1 Z (y ? y )2 dA ? 1 y2 + +@ i @y2 i 2Aj Vj 2 i

!

squares problem for the derivatives.

To avoid computing moments of each control volume in fVj gi about vi , replace x ? xi and y ? yi with (x ? xj )+(xj ? xi ) and (y ? yj )+(yj ? yi ), respectively. Expanding and integrating, we obtain

2 6 6 6 6 6 4

i 1 Z R (~ @u (x + (x ? x ) ? x ) (12) u + x ? x ~ ) = i j j i i i i Aj Vj @x i where ?y + (y ? y ) ? y c 2 + @u c y ij wij yb2 ij Lij = wij x bij wij y bij wij x i ij wij x @y i j j i (13) 2 2 + 2x (x ? x ) ? (x ? x )2 ? x2 x @ u j j j i j i i and + @x2 1 : 2 i wij = (14) 2 ? j ~ xj ? ~ x i j2 @ u + @x@y xy j + xj (yj ? yi ) + (xj ? xi)y j i 3.4 Solution of the Least-Squares +(xj ? xi )(yj ? yi ) ? xyi ) Problem 2 2 2 2 u y j + 2y j (yj ? yi ) ? (yj ? yi ) ? y i +@ In the present work, the least-squares problem of @y2 i 2 + : (9) Equation 12 is solved by using Householder transformations to reduce the left-hand side of Equation 12 to upper-triangular form. After the upper trianguThe geometric terms in this equation are of the gen- larization is complete, back-substitution yields the eral form required derivatives. There are several good reasons to use this approach instead of the simpler normal equation solution to the least-squares problem. 1 Z ((x ? x ) + (x ? x ))n nym ij xd j j i Aj Vj Using Householder transformations gives a more accurate solution to the least-squares problem ((y ? yj ) + (yj ? yi ))m dA ? xnym i m X n m n than using normal equations, especially for illX k = ( x ? x ) conditioned matrices. The error in the solution is j i k l=0 k=0 l O( K ) using Householder transformations and ? O K 2 for normal equations, where K is the (yj ? yi )l xn?k ym?l j ? xn ym i :(10) condition number of the non-square matrix and is machine precision 17]. This also implies In these terms, we can write greater robustness. As a further improvement in robustness, the 1 Z R (~ @u @u Householder transform approach can detect sinxi ) = ui + @x x bij + @y y bij Aj Vj i x ? ~ gular and nearly singular matrices on the y. If i i the least-squares problem is (nearly) singular, a b 2 2 2 2 @ u y ij x2ij + @ u x u c column with (nearly) zero elements on and below c y + +@ @x2 i 2 @x@y i ij @y2 i 2 the diagonal will be encountered during House+ : (11) holder triangularization. This failure occurs because the stencil is inadequate to support the requested number of derivatives. To resolve this, Equation 11 evaluates the mean value of the reeither more points must be added to the reconconstruction Ri (~ x?~ xi ) for a control volume j , given struction stencil or the reconstruction must be the low-order derivatives of the solution at vi and modi ed to include fewer derivatives. The latter low-order moments of the control volumes. With this course is adopted in this work. Derivatives are in hand, we can easily write down a weighted leastcomputed only to the highest order for which all

Li1 Li2 Li3 .. . LiN

0 @u @x 3B @u B @y B 7 @u B 7 B @x 7 B 7 @u B @x@y 7 5B B @u B @ @y .. .

1 2

2

2

2

1 2

2

2

1 C C C C C C C C C C A

0 w (u ? u ) i i B wi (u ? ui) B wi (u ? ui) =B B B .. @ .

1 1 2 3 2 3

wiN (uN ? ui )

1 C C C ; C C A

5 derivatives can be computed; the additional incomplete set of derivatives is discarded, because no increase in order of accuracy is possible by retaining them. After the upper triangularization of the leastsquares problem is complete, the residual for the solution is available at virtually no cost. Before back substitution, the problem looks like the following. 0 1 2x x x 3 x 0 @u 1 B r1 C @x C B r2 C 6 x x x7 B @u C B . C 6 7 B 6 @y C B . x x7 B . C 6 7 1 @ 2u C B B C 6 7 . 2 B C B C 2 @x . r 6 7 m . x B C B C 2 = 6 7 @ u B C B C r 6 7 m +1 C @x@y x B C B 6 7 2 B C B .. C 1@ u 6 7 B B 2 @y 2 C 6 7 . C B C @ A 4 5 .. 0 @ A r n ? 1 . i rn (15) If we seek m derivatives using an n point stencil, the rst m equations will be satis ed exactly. The remaining n ? m equations will not be; the ^ , which is the same as the residual for residual R the original problem, is with geometric weights and preventing overshoots by heuristically limiting, or reducing, the derivatives (e.g., 13, 14]). While this approach is not unsuccessful, it provides only a mechanical solution to an underlying theoretical problem. Speci cally, the stencil for a control volume i near a discontinuity will include control volumes j that lie on the opposite side of the discontinuity. Because the function is not smooth, approximating data in Vj by a modi ed Taylor series around vi is inappropriate. Ignoring this mathematical fact causes the unphysically large derivatives that limiting seeks to reduce. A better alternative is to reconstruct using only data that is smoothly connected to data in i. This approach is taken directly by ENO schemes, which by design search for a smooth stencil and completely exclude non-smooth data from the reconstruction. WENO schemes work less directly, using all possible stencils and weighting those containing non-smooth data with a weight that is of the order of the truncation error. In the present weighted least-squares context, weights are assigned to control volumes rather than to stencils. Nevertheless, ? the goal is to weight nonsmooth data with O xk+1 to satisfy the ENO condition of Liu, Osher, and Chan 9]. We seek to construct a data-dependent weight that will multiply the previously calculated geometric weight. This construction is based on two observations. 1. If the function is non-smooth and the neighborhood of Vi crosses a discontinuity, then a modi ed Taylor expansion does not adequately describe the function locally and there will be one or more control volumes j for which 1 Z R (~ xi )dA ? uj = O (1) : (18) Aj Vj i x ? ~ This means that the residual R of the leastsquares problem will be O (1). On the other hand, for stencils that cover only smooth regions of the function, ? Ri(~ x?~ xi) ? u(~ x) = O xk+1 (19) ? and R is also O xk+1 . Therefore, R can be used as a gauge of the smoothness of the data for the entire stencil fVj gi . 2. Whether the data in Vj is smoothly connected to data in Vi can be determined asymptotically by evaluating 8 O (1) connected ? x?1 smoothly uj ? u i = < O not smoothly j~ xj ? ~ xi j : connected (20)

v u n uX ^ rl : R=t

2

l=m+1

(16)

Scaling this by the RMS value of the geometric weight removes local mesh scale e ects: ^ (17) R = qPR 2 : wij =n

R has several uses. Within the context of this work, R will be used to compute data-dependent weightings. R also is a good measure of how

well the solution is approximated locally, making it yet another candidate for use as a re nement measure.

3.5 Reconstruction of Non-smooth Functions

The reconstruction scheme described above is designed for smooth functions. For non-smooth functions | those with O (1) discontinuities | such a scheme allows overshoots of O (1). This is not desirable for either function approximation or scienti c computation, where such overshoots can easily produce aphysical values. This problem has typically been addressed by performing a reconstruction

6 We seek a data-dependent weighting that uses R j ?ui to detect stencils with non-smooth data and ju xj ?~ ~ xi j to determine which data within that stencil should be excluded and which included. One appropriate weighting is 1 DD = Wij ; (21) (k +1) j ?ui 1 + kR ju xj ?~ ~ xi j

( +1)

and

0 = wij

j ?ui is a smoothness indicator analwhere R ju xj ?~ ~ xi j ogous to the divided di erence indicators of 9, 10] but appropriate for unstructured meshes and leastsquares reconstruction. R is computed from a data-independent least-squares reconstruction, as described above. We are concerned only with the DD , so the value of k is asymptotic behavior of Wij not critical; 0.1 seems to be a good choice. Asymptotically, the behavior of this weighting for the three important cases is 8 1 + O ? x(k+1) smooth stencil > < smooth data in stencil DD Wij = > O (1) w= nonsmooth data : O ? x?(k+1) nonsmooth data (22) That is, for stencils containing only smoothly connected data, the data-dependent weights are all approximately 1, ensuring that the good qualities of the data-independent reconstruction will be preserved for smooth functions. For stencils that are not entirely smooth, the data-dependent weight for non-smooth data is smaller than that for smooth data by a factor of the order of truncation error. These weightings satisfy the ENO condition of Liu, Osher, and Chan 9]. The data-dependent least-squares problem is closely related to the data-independent problem. The j th row in the least-squares problem of Equation 12 is modi ed by scaling with the data-dependent weight. 3.6 Summary With no further computation, the least-squares prob- The data-dependent least-squares approach can be lem becomes used to produce function reconstructions that sat0 @u 1 isfy the ENO condition. The least-squares heritage @x C 0 1 2 L0 3 B 0 of these DD-L2 -ENO schemes allows them to be ap@u wi1 (u1 ? ui ) C B i 1 @y C B plied easily to function reconstruction on unstruc0 0 2 B 6 wi2 (u2 ? ui ) C Li2 7 1@ u C B C B 7 6 2 C B 2 @x tured meshes in multiple dimensions. The algorithm 0 0 C B 6 Li3 7 B 2 u C = B wi3 (u3 ? ui ) C (23) 7 6 @ C B can be summarized as follows: C 6 .. @x@y C B A @ 5B 4 ... 7 . B 1 @2u C 2 C 0 (uN ? ui ) wiN L0iN B Input. A computational mesh, structured or @ 2 @y .. A unstructured. . i The average value of a function to be reconwhere structed for each control volume. y 0x 0c 0y 0x 0y Because data-dependent weights should all be O (()1) for c b2ij 2 L0ij = wij bij wij bij wij y ij wij ; smooth ij wij x data, a cuto value of 0.1 is used to determine whether (24) a weight is \high".

k

As data-dependent weights are assigned, the number that are numerically large is counted. If too few control volumes are assigned high data-dependent weights, insu cient smoothly connected data is available to compute the desired number of derivatives.y Similar to the data-independent case, this results in a lowering of the nominal order of accuracy of the reconstruction locally. This loss of accuracy is most likely to occur for control volumes straddling discontinuities; near the intersection of two discontinuities; or near the impingement of a discontinuity on a boundary. Because no non-smooth data has a signi cant impact on the reconstruction, however, large overshoots in the reconstruction are not expected. This would not be true if non-smooth data were included. In contrast, stencil-searching ENO schemes such as Abgrall's 12] can sometimes avoid local degradation of accuracy near a discontinuity by extending the stencil farther away from the discontinuity. However, for cases with medium resolution in which some control volumes do not have enough smoothly connected neighbors, stencil-searching schemes will use non-smooth (and therefore irrelevant) data. This approach can lead to large overshoots for such cases. Also, consider a control volume Vi that is divided by a discontinuity and therefore has an averaged function value that lies between the values on either side of the discontinuity. Here, reconstruction makes little physical sense because no smoothly connected data exists. The present scheme can detect such a situation and choose to reconstruct the solution in Vi as piecewise constant, whereas stencil searching schemes will still seek a high-order polynomial reconstruction.

j~ xj ? ~ xi j2

1

DD : Wij

(25)

7 tion, in the form of a modi ed Taylor series expansion about each vertex, valid within the control volume surrounding the vertex. Preprocessing. For each control volume Vi , nd a su ciently large set of nearby control volumes fVj gi for reconstruction of the desired order. Preprocessing. Compute all control volume moments that will be needed in the least-squares problem of Equation 12. Reconstruction. For each control volume each time a function is reconstructed: 1. Compute geometric weights, and construct the arrays needed for the data-independent least-squares problem of Equation 12. 2. Solve the least-squares problem using Householder transformations. The solution algorithm should not destroy the original arrays and should return both the solution and the residual of the least-squares problem, along with information about how many derivatives were actually calculable. 3. Using the residual of the least-squares problem, compute a data-dependent weight for each control volume in the stencil and multiply the appropriate row of the original least squares problem by this weight. Keep track of how many high weights there are, since this limits the number of derivatives that can be plausibly calculated. 4. Solve the least-squares problem again. The derivatives computed give an ENO reconstruction when substituted into Equation 2.

Output. An ENO reconstruction of the func-

10

-1

L1 Norm of Reconstruction Error

10

-2

10

-3

10

-4

10

-5

Second Order Third Order Fourth Order 100 1000 10000 Number of Vertices 100000

10

-6

10

-1

L2 Norm of Reconstruction Error

10

-2

10

-3

10

-4

10

-5

Second Order Third Order Fourth Order 100 1000 10000 Number of Vertices 100000

10

-6

Convergence Order of Error Norms Nominal Order 2 3 4 L1 2.23 2.90 4.12 Norm L2 2.16 2.89 3.99

2: Convergence of k-Exact Reconstruction of The reconstruction should be uniformly high-order Figure a Smooth Function accurate for smooth functions. To show this, we have reconstructed the function (26) function within each control volume. Sixth-order quadrature is used to compute the error norms in on the square 0; 1] 0; 1] using a series of unstruc- order eliminate quadrature error in the computatured triangular meshes ranging in size from 491 to tion oftothe norms. The accompanying table veri es 26651 vertices. The L1 and L2 norms of the error that the expected asymptotic convergence rates are in reconstruction of this function are shown in Fig- achieved. ure 2 for second- through fourth-order accuracy. Error norms are computed by using the di erence be- The new reconstruction scheme was also tested for tween the analytic function and the reconstructed non-smooth function reconstruction in two dimenu(x; y) = cos( x2 + 4 y)

4 Function Reconstruction

8

1.0 First Order Second Order

0.5

0.0

-0.5

-1.0 -1.0

-0.5

0.0

0.5

1.0

Figure 5: Control Volumes That Fail to Attain ThirdFigure 3: Contours of Function De ned by Equa- Order Accuracy (k = 2) tion 27.

1.0 1.0 First Order Second Order Third Order 0.5 0.5

0.0

0.0

-0.5

-0.5

-1.0 -1.0

-0.5

0.0

0.5

1.0

-1.0 -1.0

-0.5

0.0

0.5

1.0

Figure 4: Control Volumes That Fail to Attain Figure 6: Control Volumes That Fail to Attain Fourth-Order Accuracy (k = 3) Second-Order Accuracy (k = 1) sions. The function chosen is that of Abgrall 18].z p x cos2 y (x ? cot p 2 y) u(x; y) = f f (x + cot 2 y) + cos (2 y) x > cos2 y (27) with 8 ?r sin ? 3 r2 r ?1 < 2 3 (28) jrj < 1 f (r) = : jsin (2 r)j 3 2r ? 1 + 1 sin (3 r) r 1 6 3 Contour plots of this function are shown in Figure 3 Because the reconstruction scheme is k-exact (with or without data-dependent weighting), we know that the accuracy of the reconstruction in smooth regions of the function will be of order k +1. Norms of the error in the reconstruction are meaningless in this context, because the di erence between the actual function and its reconstruction is guaranteed to be O (1) in control volumes that are crossed by a discontinuity. Instead, Figures 4{6 show the control volumes for which the nominal order of accuracy was not obtained on an isotropic triangular mesh with 26651 vertices for second through fourth order. The discontinuities in the function are clearly visible in these gures as control volumes with reduced reconstruction accuracy. In all control volumes away from the z Several typographical errors in the de nition of this func- discontinuities, the nominal order of accuracy is attion in 18] cause a mismatch with the plotted function there; tained. Table 1 shows the distribution of actual rethe function shown here matches the plots in 18]. construction accuracy for each case.

9 Table 1: Actual Accuracy of Reconstruction for Abgrall's Function Nominal Achieved Order Order 1 2 3 4 2 245 26406 | | 3 233 82 26336 | 4 290 108 300 25953

Canonical Second Order Theoretical Second Order

Canonical Third and Fourth Order

Theoretical Third and Fourth Order

5 Flow Solver Implementation

In addition to reconstruction, there are several other details in the construction of a high-order accurate Figure 7: Choices for Gauss Integration Points for ow solver which require careful attention. the Median Dual

5.1 Flux Quadrature

After the solution has been reconstructed from control-volume averages to a polynomial in each control volume, the second phase of the ow solution procedure is computation of a ux integral, or residual, for each control volume. This integration must be done to at least the same order of accuracy as the solution reconstruction to obtain high-order accuracy. In the present work, the integration is performed by Gaussian quadrature 19] around the boundary of each control volume. Gaussian quadrature has the property that an N -point quadrature along a line segment is 2N -order accurate. Accordingly, rst- and second-order accurate schemes (k = 0; 1) use N = 1, and third- and fourth-order schemes (k = 2; 3) use N = 2. Unfortunately, for the median dual this is not quite the full story, since the control volume boundary separating two adjacent volumes is not a single segment but two segments, each extending from the middle of an edge of the mesh to the centroid of a triangle. For rst and second order, this is not a problem; a single point can still be used. Although in practice using the mid-edge seems adequate, in principle this point should be the centroid of the pair of line segments. For third and fourth order, two points on each part of the segment are, in principle, required. These options are illustrated in Figure 7. In both cases, the present work sticks to the letter of the law, using the quadrature points shown on the right side of the gure. As we shall see, computational time is dominated by reconstruction for high-order, so the additional ux computations are not a severe time penalty. At each quadrature point, the ux is evaluated using Roe's approximate Riemann solver 20].

A piecewise linear boundary description, as given by a polygonal input set, gives a second-order accurate representation of the shape of a smoothly curved boundary. For higher-order schemes, a more accurate boundary representation is required. The approach taken here is to specify, for each boundary segment, the inward normal at each end, as shown in Figure 8. This is equivalent to specifying slopes

5.2 Boundary Description and Implementation

Figure 8: Schematic of High-Order Boundary Description at the end of each segment and allows the exibility to include sharp corners, because normals are de ned edgewise rather than pointwise. This slope information is used to construct a piecewise cubic representation of the surface, which is fourth-order accurate for smoothly curved boundaries. This cubic representa-

10 tion is used to determine the Gauss integration data and control volume moments for control volumes on the boundary. Ultimately the question of practicality of high-order schemes will depend on trade-o s between accuracy and computational requirements. As preliminary data towards settling this question, requirements for CPU time and memory for the present implementation are tabulated here. Table 2: Resource Requirements for High-Order Euler Solution on Two-Dimensional Unstructured Meshes

Resource Order of Accuracy Usage 2 3 4 a CPU Time ( sec / vertex / evaluation) Reconstruction 255 1580 2180 Interior Flux Quad 130 573 649 Boundary Flux Quad 5.5 12.4 12.4 Memory (words / vertex)b Solution 4f 4f 4f Derivatives 8f 20f 36f Neighbor Informationc 7i 19i 19i Gauss Point Locationsd 6f 24f 24f Gauss Point Normalsz 6f 12f 12f Gauss Weightsz 3f 6f 6f Total 27f 84f 100f +7i +19i +19i

Error in Computed Density

-1

10

5.3 Computational Resources

10

-2

10

-3

10

-4

L1, Second L2, Second L1, Third L2, Third L1, Fourth L2, Fourth

10

-5

10

100 1000 Number of Vertices

10000

Convergence Order of Error Norms Nominal Order 2 3 4 L1 1.79 2.88 3.64 Norm L2 1.78 2.87 3.61

a On a 110 MHz SPARC 5 b f = oating point, i = integer c Average for interior vertices d Assuming three times as many edges

Figure 9: Error in Density for Supersonic Vortex Flow used for this case, containing 73, 253, 941, and 4398 vertices. Error norms are computed from the solution in the same way as for reconstruction. Figure 9 shows the L1 and L2 norms of the error in density for this problem, and the accompanying table indicates the asymptotic order of accuracy achieved. To verify that the scheme behaves well for ows with weak shocks, AGARD test case 1 21] was computed on a mesh with 4156 vertices, shown in Figure 10. Because the current scheme is an ENO scheme, we are interested not only in the solutions but also in the convergence behavior, to verify that it is possible to converge to machine zero. Figure 11 shows the surface pressure coe cients for this case for second- and third-order accuracy. There is little visible di erence in the solutions at this scale. Each solution shows a mild oscillation near the shock. Of note also is the behavior near the leading edge stagnation point. In

as vertices

6 Flow Solutions

6.1 Supersonic Vortex Flow

A smooth ow was calculated to validate the accuracy of the high-order boundary shape de nition and ux quadrature. The ow chosen was a supersonic vortex ow in a quarter-circular annulus. This is an isentropic ow in which the velocity is inversely proportional to radius and the density is given by = i 1 M 2 1 ? ri2 1+ ? 2 i r2

1 ? 1

6.2 Transonic Airfoil Flow

;

(29)

where the subscript i denotes quantities at the inner boundary. For this case, ri was chosen to be 2, ro = 3, and Mi = 2. Four isotropic triangular meshes were

11

Figure 10: Mesh for AGARD Test Case 1 (4156 ver- Figure 12: Stagnation Pressure Coe cient Near tices) Leading Edge for AGARD Test Case 1 (Second Order)

-1.6

-1.2

Surface Pressure Coefficient

-0.8

-0.4

0.0

Second Order Third Order

0.4

0.8

1.2 0.0

0.2

0.4 x/c

0.6

0.8

1.0

Figure 11: Surface Pressure Coe cient for AGARD Figure 13: Stagnation Pressure Coe cient Near Test Case 1 Leading Edge for AGARD Test Case 1 (Second Order) this region, stagnation pressure loss is often relatively high because of poor resolution of the rapidly turning ow. Figures 12 and 13 demonstrate a marked value for second order is about 0.072, while for third improvement in the losses at the leading edge when order it is about 0.036. going from second to third order. In each gure, con- Figure 14 shows the convergence history for these tour levels are separated by 0.04 and the contours cases. Multigrid W-cycles were used in conjunction away from the body are at 0 or 0.04. The peak with local time step and a three-stage time-stepping

12

0.8650

0

10 10 10 10 Maximum Density Residual 10 10 10 10 10 10 10 10 10

-1

-2

0.3986 0.3460 0.2811 0.1211

-3

-4

-5

-6

0

1.0880

2.1071 2.1998

2.9410

3.3940

3.9899

-7

-8

Figure 15: Scramjet Geometry

Second Order Third Order Fourth Order

-9

-10

-11

-12

0

100

200 MG Cycles

300

400

Figure 14: Convergence Histories, AGARD Test Case 1 scheme. No serious e ort has been made to optimize the multistage scheme or the CFL number. The gure shows that convergence for the third-order scheme is not as good as for second order, but both converge well with none of the hang-ups often seen with stencilsearching ENO schemes. The fourth-order scheme does hang up, with oscillations in shock position and shape being the culprit. Further study is needed to eliminate this problem. The scramjet con guration introduced by Kumar 22] was computed at M=5, = 0, solely as a robustness demonstration. The geometry for this case is shown in Figure 15. A nearly uniform, isotropic triangular mesh with 8841 vertices was generated. Figure 16 shows density contours for a second-order accurate solution of this problem. Much of the detail of the ow is missing because of poor resolution (for example, the throat is only nine cells across), but the broad outlines of the shock re ections and interactions are present. Clearly, local re nement would be of tremendous bene t in resolving this ow.

Figure 16: Scramjet Density Contours ENO scheme. Reconstruction accuracy for smooth functions has been veri ed, and the reconstruction scheme behaves gracefully near function discontinuities. Accuracy for smooth ows has been established by comparing computation with an analytic solution. The scheme has been shown to be capable of converging to machine zero. Finally, the scheme's ability to handle complex, high Mach number ows robustly has been demonstrated. A number of open questions remain. Extension to a three-dimensional ow solver will likely require switching to a cell-centered scheme because of the complex shape of vertex-centered control volumes. Might it also be advantageous to use a cell-centered scheme in two dimensions? In smooth regions of the ow, data-dependent reconstruction is unnecessary; an excellent representation of the solution is provided by the dataindependent reconstruction. Can the scheme recognize and exploit this feature without causing convergence to hang up by switching? In very at regions of the ow, the use of fourthorder accuracy is, in a practical sense, overkill. Second- or even rst-order accuracy may be sufcient to represent the solution. How must a p-re nement scheme be designed to choose optimal accuracy locally? Some convergence and robustness questions re-

6.3 Scramjet Con guration

7 Conclusions

This article has demonstrated the feasibility of computing high-order accurate solutions to the Euler equations on unstructured triangular meshes using an

13 main. Certain cases violate positivity, and con- 11] Durlofsky, L. J., Enquist, B., and Osher, S., vergence rates are not particularly good. \Triangle Based Adaptive Stencils for the Solution of Hyperbolic Conservation Laws," Journal of Computational Physics, vol. 98, pp. 64{73, Jan. 1992. 1] van Leer, B., \Towards the Ultimate Conserva- 12] Abgrall, R., \On Essentially Non-oscillatory tive Di erence Scheme. V. A Second-Order SeSchemes on Unstructured Meshes: Analysis quel to Godunov's Method," Journal of Compuand Implementation," Journal of Computational tational Physics, vol. 32, pp. 101{136, 1979. Physics, vol. 114, no. 1, pp. 45{58, 1994. 2] Harten, A. and Osher, S., \Uniformly High- 13] Barth, T. J., \Recent Developments in High Order Accurate Nonoscillatory Schemes," SIAM Order K-Exact Reconstruction on Unstructured Journal on Numerical Analysis, vol. 24, pp. 279{ Meshes." AIAA paper 93-0668, Jan. 1993. 309, Apr. 1987. 14] Venkatakrishnan, V., \On the Accuracy of Lim3] Harten, A., Osher, S., Engquist, B., and iters and Convergence to Steady-State SoluChakravarthy, S. R., \Some Results on Unitions." AIAA paper 93-0880, Jan. 1993. formly High-Order Accurate Essentially NonOscillatory Schemes," Applied Numerical Math- 15] Jameson, A., \Arti cial Di usion, Upwind Biasing, Limiters and their E ect on Accuracy ematics, vol. 2, pp. 347{377, 1986. and Multigrid Convergence in Transonic and Hy4] Harten, A., Enquist, B., Osher, S., and personic Flows." AIAA paper 93-3359-CP, July Chakravarthy, S. R., \Uniformly High Order 1993. Accurate Essentially Non-oscillatory Schemes, III," Journal of Computational Physics, vol. 71, 16] Ollivier-Gooch, C. F., \Quasi-ENO Schemes for Unstructured Meshes Based on Unlimpp. 231{303, Aug. 1987. ited Data-Dependent Least-Squares Reconstruc5] Shu, C.-W. and Osher, S., \E cient Impletion," Journal of Computational Physics. To apmentation of Essentially Non-oscillatory Shockpear. Capturing Schemes, II," Journal of Computa17] Golub, G. H. and van Loan, C. F., Matrix Comtional Physics, vol. 83, pp. 32{78, 1989. putations. Baltimore, Maryland: The Johns 6] Shu, C.-W., Zang, T. A., Erlebacher, G., Hopkins University Press, 1983. Whitaker, D., and Osher, S., \High-order ENO Schemes Applied to Two- and Three- 18] Abgrall, R., \Design of an Essentially Nonoscillatory Reconstruction Procedure on FiniteDimensional Compressible Flow," Applied NuElement Type Meshes." ICASE Report No. merical Mathematics, vol. 9, no. 1, pp. 45{71, 91-84, NASA Langley Research Center, 1991. 1992. NASA CR 189574. 7] Casper, J., An Extension of Essentially NonOscillatory Shock Capturing Schemes to Multi- 19] Stroud, A. H. and Secrest, D., Gaussian QuadraDimensional Systems of Conservation Laws. ture Formulas. Englewood Cli s, N.J.: PrenticePhD thesis, Old DominionUniversity, Dec. 1990. Hall, 1966. 8] Godfrey, A. G., Mitchell, C. R., and Wal- 20] Roe, P. L., \Approximate Riemann Solvers, Paters, R. W., \Practical Aspects of Spatially rameter Vectors, and Di erence Schemes," JourHigh-Order Accurate Methods," AIAA Journal, nal of Computational Physics, vol. 43, pp. 357{ vol. 31, pp. 1634{1642, Sept. 1993. 372, 1981. 9] Liu, X.-D., Osher, S., and Chan, T., \Weighted 21] AGARD Fluid Dynamics Panel, Test Cases for Essentially Non-oscillatory Schemes," Journal of Inviscid Flow Field Methods. AGARD, May Computational Physics, vol. 115, pp. 200{212, 1985. AGARD Advisory Report AR-211. 1994. 22] Kumar, A., \Numerical Analysis of the 10] Jiang, G.-S. and Shu, C.-W., \E cient ImpleScramjet-Inlet Flow Field by Using Twomentation of Weighted ENO Schemes," Journal Dimensional Navier-Stokes Equations," Techniof Computational Physics, vol. 126, pp. 202{228, cal Paper 1940, NASA, December 1981. June 1996.

References