qll v '?'% _ --- NASA Contractor Report 187502 ICASE Report No. 91-8 ICASE RECENT DEVELOPMENTS IN SHOCK-CAPTURING SCHEMES Ami Harten Contract No. NAS 1-18605 January 1991 Institute for Computer Applications in Science and Engineering NASA Langley Research Center Hampton, Virginia 23665-5225 Operated by the Universities Space Research Association IW A National AeronaUlic.q and Spare Adminisfration I__ngley Rese_rcli Center I"larnplon, Virginia 23665-5225 "(_ASA-C;'-187502) R_c,'r:NT DEVELOPMENTS IN N91-17_£9 SNOCK-CAPTURING SCHEMES Fin,_l Report (ICASF) ]z, p _. CSCL ]2A UncI as G3/04 0332310 Recent Developments in Shock-Capturing Schemes Ami Harten 1 School of Mathematical Sciences Tel- Aviv University Tel-Aviv ISRAEL and Department of Mathematics University of California Los Angeles, CA 90024-1555 ABSTRACT In this paper we review the development of the shock-capturing methodology, paying special attention to the increasing nonlinearity in its design and its relation to interpolation. It is well-known that high-order approximations to a discontinuous function generate spurious oscillations near the discontinuity (Gibbs phenomenon). Unlike standard finite-difference methods which use a fixed stencil, modern shock- capturing schemes use an adaptive stencil which is selected according to the local smoothness of the solution. Near discontinuities this technique automatically switches to one-sided approximations, thus avoiding the use of discontinuous data which brings about spurious oscillations. 1Research was supported by the National Aeronautics and Space Administration under NASA Contract No. NAS1-18605 while the author was in residence at the Institute for Computer Ap- plications in Science and Engineering (ICASE), NASA Langley Research Center_ Hampton, VA 23665. 1 Introduction In this paper, we describe and analyze numerical techniques that are designed to approximate weak solutions of hyperbolic systems of conservation laws in several space dimensions. For sake of exposition, we shall describe these methods as they apply to the pure initial value problem (IVP) for a one-dimensional scalar conservation law u,+/(u)_=0, u(x,0)=uo(z). (i.i) To further simplify our presentation, we assume that the flux/(u) is a convex function, i.e., f"(u) > 0 and that the initial data u0(x) are piecewise smooth functions which are either periodic or of compact support. Under these assumptions, no matter how smooth uo is, the solution u(x,t) of the IVP (1.1) becomes discontinuous at some finite time t = t=. In order to extend the solution for t > t=, we introduce the notion of weak solutions, which satisfy d/b (i.2.) -gi ,_d:_+f(,.,(b,t))- f(Ka,t))=0 for all b > a and t > 0. ielation (i.2a) implies that u(z, t) satisfies the PDE in (i.1) wherever it is smooth, and the iankine-Hugoniot jump relation f(,_(y+o,t))- f(_,(y- o,t))= [,_(y+ o,t)- ,_(y- o,t)]-_ (1.2b) across curves z = y(t) of discontinuity. It is well-known that weak solutions are not uniquely determined by their initial data. To overcome this difficulty, we consider the IVP (1.1) to be the vanishing viscosity limit e I 0 of the parabolic problem (,_'),+f(,_'),,=_,L _,o(x,o=) _,o(_.), (l.3a) and identify the unique "physically relevant" weak solution of (1.1) by limut The limit solution (1.3) can be characterized by an inequality that the values uL = u(y- O,t),uR = u(y + O,t) and s = dy/dt have to satisfy; this inequality is called an entropy condition; admissible discontinuities are called shocks. When f(u) is convex, this inequality is equivalent to Lax's shock condition (1.4) where a(u) = f'(u) is the characteristic speed (see [81for more details). Weturn nowtodescribefinitedifferenceapproximationsforthe numericalsolution of the IVP (1.1). Let v_'denote the numerical approximation to u(x_,t,,) where x_= jh, t. = nT";let Vh(X,t) be a globally defined numerical approximation associated with the discrete values {v_'},oo < j < co,n >_0. The classical approach to the design of numerical methods for partial differential equations is to obtain a solvable set of equations for {v_'} by replacing derivatives in the PDE by appropriate discrete approximations. Therefore, there is a conceptual difficulty in applying classical methods to compute solutions which may become dis- continuous. Lax and Wen&off [9] overcame this difculty by considering numerical approximations to the weak formulation (1.2a) rather than to the PDE (1.1). For this purpose, they have introduced the notion of schemes in conservation form: = - (Tj+-½7j_ -) :)j; (l.5a) V?+I here A = -r/h and f_+½ denotes 7,+½= v,"+k); (1.hb) 7(wl,..., w2k) is a numerical flux function which is consistent with the flux f(u), in the sense that 7(u,,,,...u) =/(,,); (1.5c) Eh denotes the numerical solution operator. Lax and Wendroff proved that if the numerical approximation converges boundedly almost everywhere to some function u, then u is a weak solution of (1.1), i.e., it satisfies the weak formulation (1.2a). Consequently discontinuities in the limit solution automatically satisfy the Rankine- Hugoniot relation (1.2b). We refer to this methodology as shock-capturing (a phrase coined by H. Lomax). In the following, we list the numerical flux function of various 3-point schemes (k = l in (1.5b)): (i) The Lax-Friedrichs scheme [7] (1.6) -f(w,,w2) = 2[f(wl) + f(w2) - _(w2 - wl)] (ii) Godunov's scheme [1] 7(w_,w2) f(V(O; w,,w2)); (1.7a) here V(m/t;w,, w2) denotes the self-similar solution of the IVP (1.1) with the initial data ?AO(X) L w_ m>O (iii) The Cole-Murman scheme [12]: (1.8a) f(wl,w_) = _[/(Wl) + f(w_) - [_(wl,w2)](_2 - Wl)] where if wl # w2 (1.8b) { $(,,,=)-1(_,,) if wl = w2 2 (iv) The Lax-Wendroffscheme[9]: f(wl,w2) = _{f(wl) -Ff(w2) --ha wl --Fw_ [f(w_)--f(wl)]}. (1.9) 2 Let E(_) denote the evolution operator of the exact solution of (1.1) and let Eh denote the numerical solution operator defined by the RHS of (1.5a). We say that the numerical scheme is r-th order accurate (in a pointwise sense) if its local truncation error satisfies E(r).u- Zh.u = 0(h "+1) • (1.10) for all sufficiently smooth u; here r = 0(h). If r > 0, we say that the scheme is consistent. The schemes of Lax-Friedrichs (1.6), Godunov (1.7), and Cole-Murman (1.8) are first order accurate; the scheme of Lax-Wendroff (1.9) is second order accurate. We remark that the Lax-Wendroff theorem states that if the scheme is convergent, then the limit solution satisfies the weak formulation (1.2b); however, it need not be the entropy solution of the problem (see [4]). It is easy to see that the schemes of Cole-Murman (1.8) and Lax-Wendroff (1.9) admit a stationary "expansion shock" (i.e., f(uL) = f(uR) with a(uL) < a(uR)) as a steady solution. This problem can be easily rectified by adding sufficient numerical dissipation to the scheme (see [11] and [3]). 2 Interpolatory Schemes and Linear Discontinu- ities Let us consider the constant coefficient case f(u) = au, a-" const, in (1.1), i.e., + = 0, 0)= (2.1a) the solution to which is In this case the schemes (1.6)- (1.9) take the form K + v'_+I = __, Cl(v)v']l+- (Eh" v")j (2.2) l=-K- where v = )_a is the CFL number. The coefficients Cl(v) are independent of the numerical solution v"; this makes Eh a linear operator. We say that the numerical scheme Eh is (linearly) stable if [l(Eh)'_l] < C for 0_ n_" < T, r = O(h). (2.3a) In the constant coefficient case the scheme is stable if and only if it satisfies von Neumann's condition K + [ _ Ct(v)l_[ <1 for all 0<_<_r. (2.3b) l=-K- 3 It is easy to see that all the 3 point schemes (1.6) - (1.9) are stable under the CFL condition Iv[--l_al __ I. (2.3c) | The notion of stability (2.3a) is related to convergence through Lax's equivalence theorem, which states that a consistent linear scheme is convergent if and only if it is stable (see [13] for more details). Let us denote by SI" the stencil of (r + 1) successive points starting with zi (2.4.) _'i = {xi, xi+l,... ,xi+,}, let P(x; S_, u) denote the unique polynomial of degree r interpolating the (r+ 1) vvaalulueses of u on this stencil and let Q(x; u) denote the piecewise polynomial interpolation_ion ooff (2.4b) = S_(j), xj-1 _ _ We refer to the numerical scheme v$+1= Q(_j-,,;v") (2.4c) as interpolatory scheme. Clearly, the interpolatory scheme (2.4) is r-th order accumratea.te. When Q(x; v") is the piecewise linear interpolation of v" (i.e., r = 1, i(j) = j - 1 in (2.45)) then (2.4c) is the first-order accurate upwind scheme; in the conasstan;atnt coefficient case this scheme is identical to those of Godunov (1.7) and Cole-Murmirlnaan, n (1.8). Next let us assume a > 0 and consider the second order case r = 2 in which eiii:c_iiwo::asicohf i Q(x; v") is a piecewise-parabolic interpolation of v'_. There are two different choices of stencil in (2.4): Taking Q in Ixj_ 1, x_] to be the parabola through S]_1 = {X#_I, al#, ali+l } (i.e., i(j) = j - 1) results in the Lax-Wendroff scheme (1.9); taking Q in [zi_l,xj] to be the parabola through S]_2 = {xj__,x__l,zj} (i.e., (i(j) = j - 2) results in the second-order upwind scheme. We turn now to consider the application of these schemes to the step function H(z) {0 j<0 1 _>o ,H_= 1 j>l_ For the first order upwind scheme we get that ; (2.55) Q(x;H)= /h O < :r. < h i xh<<0_ for the Lax-Wendroff scheme x<-h lz x -h<z<O O(x; H) = 2h(i + X) ; (2.5c) 1- ½(1- _)(2- _) O<z<h 0 1 h<z 4 for the secondorderupwind schemewegetthat _(1+ _) 0 < x < h (2.5d) _- 1)(2- • h < 0 z_<._ Q(z; H) = 11-4-](_ X) 2h 2h We observe that Q in (2.5b) is a monotone function of x; consequently the numerical solution by Godunov's scheme to these data is also monotone. On the other hand Q for the second order schemes (2.5c) - (2.5d) is not a monotone function. For the Lax-Wendroff scheme Q is negative in -h < x < 0 and has a minimum of -0.125; similarly for the second order upwind scheme Q is larger than 1 in h < z < 2h with a maximum of 1.125. This observation explains the Gibbs-like phenomenon of generating spurious oscillations in calculating discontinuous data with these second order schemes. We say that the scheme Eh is monotonicity preserving if v monotone =_ Eh. v monotone. (2.5) Clearly the numerical solution of a monotonicity preserving scheme to initial data of a step-function is always monotone and therefore the discontinuity propagates without generating spurious oscillations. Godunov has shown that the linear scheme (2.2) is monotonicity preserving if and only if Cl(v)_>0, -g_ <l<K+; (2.7) this implies that a monotonicity-preserving scheme which is linear is necessarily only first-order accurate. It took some time to realize the Godunov's monotonicity the- orem does not mean that there are no high-order accurate monotonicity preserving schemes; it only means that there are no such linear ones. Hence high-order accurate monotonicity-preserving schemes are nonlinear in an essential way. The second-order accurate schemes mentioned above are linear because the choice of the stencil (2.4) is fixed. Let us consider now a piecewise-quadratic interpolation which is made nonlinear by an adaptive selection of the stencil in (2.4b). For the interval [__1, _j] let us consider the two stencils S]_ 2 = {xj_2, zj_,, x_} and S]_ 1 = (Zi_l, z j, zj+l }7and select the one in which the interpolant is smoother. If we measure the smoothness of u by the second derivative of the corresponding parabola we select d_ j - 2 if [_-_P(z; __2,u)] < [-_P(x;S]_I,U)[ (2.s.) i(j) i j-1 otherwise When we apply this selection of stencil to the step-function H(z) (2.5a) we get that for[_-1,_0]wechoosethestencilS___= {___,___,_0)forwhichP(_;S___H,) =0;for the interval [zl, x2] we choose the stencil S_ = {x_, z2, x3} for which P(z; S_, H) _ 1. As is evident from comparing (2.5c) and (2.5d) it does not matter which stencil we assignto [x0,xl]sinceboth parabolae aremonotone there;with (2.8a)we select S2 -1 for[x0,xl].Thus we get in (2.4) Q(x;H)= _1K=(1+_) 0<x-<h - (2.86) 0 x_<0 1 h_<x which is a monotone function of x although it is actually a plecewise-quadratic poly- i nomial. The use of an adaptive stencil is the main idea behind the Essentially Non- Oscillatory (ENO) schemes to be described later in this paper. It extends to high order of accuracy in a straightforward manner: For r-th order accu:ruacyacy wwee c,:oonrsidsiedrer for [zi_l,zj] the r stencils S" S r S;__. We choose i(j) iinn ((22..44bb)) tto(, bbe_ tthhee j--r, i-r+l,'" ", one which minimizes Id- P(x; ,u)l fori=j-r,...,j- 1. (2.9) 3 Total Variation Stability and TVD Schemes An immense body of work has been done to find out whether stability of constant coefficient scheme with respect to all "frozen coefficients" associated with the problem, implies convergence in the variable coefficient case and in the nonlinear case. In the variable coefficient case, where the numerical solution operator is linear and Lax's equivalence theorem holds, it comes out that the stability of the variable coefficient scheme depends strongly on the dissipativity of the constant coefficient one, i.e., on the particular way it damps the high-frequency components in the Fourier representation of the numerical solution. In the nonlinear case, under assumptions of sufficient smoothness of the PDE, its solution and the functional definition of the numericai scheme, Strang proved that linear stability of the first variation of the scheme implies its convergence; we refer the reader to [13] for more details. In the case of discontinuous solutions of nonlinear problems, linearly stable schemes are not necessarily convergent; when such a scheme fails to converge, we refer to this case as "nonlinear instability." The occurrence of a nonlinear instability is usually associated with insufficient numerical dissipation which triggers exponential growth of the h{gh-frequency components of the numerical solution. The following theorem states that a stronger sense of stability, namely uniform boundedness of the total variation of the numerical solution, does imply convergence to a weak solution. Theorem 3.1. Let Vh be a numerical solution of a conservative scheme (1.5). (i) If TV(vh(.,t)) < C. TV(uo) (3.1) where TV( ) denotes the total variation in x and C is a constant independent of h for 0 < t < T, then any refinement sequence h _ 0 with T = 0(h) has a convergent subsequence hi --+ 0 that converges in --lr)°¢to a weak solution of (1.1). 6