# Full text of "Finite Volume Methods: Foundation and Analysis"

## See other formats

Finite volume methods: foundation and analysis Timothy Barth 1 and Mario Ohlberger 2 1 NASA Ames Research Center, Information Sciences Directorate, Moffett Field, California, 94035, USA 2 Institute of Applied Mathematics, Freiburg University, Hermann-Herder-Str. 10, 79104 Freiburg, Germany and CSCAMM, University of Maryland, 4127 CSIC Building, College Park, Maryland, 20742-3289, USA ABSTRACT Finite volume methods axe a class of discretization schemes that have proven highly successful in approximating the solution of a wide variety of conservation law systems. They are extensively used in fluid mechanics, porous media flow, meteorology, electromagnetics, models of biological processes, semi-conductor device simulation and many other engineering areas governed by conservative systems that can be written in integral control volume form. This article reviews elements of the foundation and analysis of modern finite volume methods. The primary advantages of these methods are numerical robustness through the obtention of discrete maximum (minimum) principles, applicability on very general unstructured meshes, and the intrinsic local conservation properties of the resulting schemes. Throughout this article, specific attention is given to scalar nonlinear hyperbolic conservation laws and the development of high order accurate schemes for discretizing them. A key tool in the design and analysis of finite volume schemes suitable for non-oscillatory discontinuity capturing is discrete maximum principle analysis. A number of building blocks used in the development of numerical schemes possessing local discrete maximum principles are reviewed in one and several space dimensions, e.g. monotone fluxes, E-fluxes, TVD discretization, non- oscillatory reconstruction, slope limiters, positive coefficient schemes, etc. When available, theoretical results concerning a priori and a posteriori error estimates are given. Further advanced topics are then considered such as high order time integration, discretization of diffusion terms and the extension to systems of nonlinear conservation laws. KEY words: finite volume methods, conservation laws, non-oscillatory approximation, discrete maximum principles, higher order schemes Contents 1 Introduction: Scalar nonlinear conservation laws 2 1.1 Characteristics of scalar conservation laws 3 1.2 Weak solutions 4 1.3 Entropy weak solutions and vanishing viscosity 5 1.4 Measure- valued or entropy process solutions 6 2 Finite volume (FV) methods for nonlinear conservation laws 7 2.1 Godunov finite volume discretizations 8 2.1.1 Monotone schemes 10 Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Ren6 de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 2 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 2.1.2 E-flux schemes 11 2.2 Stability, convergence and error estimates 12 2.2.1 Discrete maximum principles and stability 12 2.2.2 Convergence results 13 2.2.3 Error estimates and convergence rates 16 2.2.4 A posteriori error estimate 17 2.2.5 A priori error estimate 17 2.2.6 Convergence proofs via the streamline diffusion discontinuous Galerkin finite element method 18 3 Higher order accurate FV generalizations 19 3.1 Higher order accurate FV schemes in 1-D 19 3.1.1 TVD schemes 20 3.1.2 MUSCL schemes 22 3.1.3 ENO/WENO schemes 24 3.1.4 Reconstruction via primitive function 25 3.1.5 ENO reconstruction .25 3.1.6 WENO reconstruction 26 3.2 Higher order accurate FV schemes in multi-dimensions 27 3.2.1 Positive coefficient schemes on structured meshes 28 3.2.2 FV schemes on unstructured meshes utilizing linear reconstruction. ... 30 3.2.3 Linear reconstruction operators on simplicial control volumes 35 3.2.4 Linear reconstruction operators on general control volumes shapes. ... 36 3.2.5 General p-exact reconstruction operators on unstructured meshes. ... 38 3.2.6 Positive coefficient schemes on unstructured meshes 39 4 Further Advanced Topics 41 4.1 Higher order time integration schemes 41 4.1.1 Explicit SSP Runge-Kutta methods 42 4.1.2 Optimal second and third order nonlinear SSP Runge-Kutta methods. . 43 4.2 Discretization of elliptic problems 43 4.3 Conservation laws including diffusion terms 45 4.3.1 Choices of the numerical diffusion flux djk 45 4.3.2 Note on stability, convergence and error estimates 46 4.4 Extension to systems of nonlinear conservation laws 47 4.4.1 Numerical flux functions for systems of conservation laws 48 5 Concluding Remarks 51 1. Introduction: Scalar nonlinear conservation laws Many problems from physics, chemistry, biology, mechanics, and gas dynamics lead to the study of nonlinear hyperbolic conservation laws. As a prototype conservation law, consider the Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 3 Cauchy initial value problem d t u + V • f(u) = in R d x K + , (la) u(x,0) = u {x) inK d . (lb) Here u(x,t) : R d x R + ->■ R denotes the dependent solution variable, f(u) € C l (R) denotes the flux function, and uo(x) : R d -> R the initial data. The function u is a classical solution of the scalar initial value problem if u € C 1 (R d x K + ) satisfies (la-lb) pointwise. An essential feature of nonlinear conservation laws is that, in general, gradients of u blow up in finite time, even when the initial data uo is arbitrarily smooth. Beyond some critical time to classical solutions of (la-lb) do not exist. This behavior will be demonstrated using the method of characteristics. By introducing the notion of weak solutions of (la-lb) together with an entropy condition, it then becomes possible to define a class of solutions where existence and uniqueness is guaranteed for times greater than to. These are precisely the solutions that are numerically sought in the finite volume method. 1.1. Characteristics of scalar conservation laws Let u be a classical solution of (la) subject to initial data (lb). Further, define the vector a(u)=f'(u) = (f[(u);---,f' d (u)) T . A characteristic T y is a curve (x(t),t) such that x'(t) = a(u(x(t),t)) for*>0 , x(0) = y . Since u is assumed to be a classical solution, it is readily verified that —u(x(t),t) = d t u + x'(t)Vu = d t u + a(u)Wu = d t u + V ■ f(u) = . at Therefore, u is constant along a characteristic curve and T y is a straight line since x'(t) = a(u(x(t),t)) = a{u(x{0),0)) = a(u(y,0)) = a(u (y)) = const . In particular x(t) is given by x(t) = y + ta(u (y)) ■ (2) This important property may be used to construct classical solutions. If a; is fixed and y determined as a solution of (2) , then u(x,t) -u {y) . This procedure is the basis of the so-called method of characteristics. On the other hand, this construction shows that the intersection of any two straight characteristic lines leads to a contradiction in the definition of u(x, t). Thus, classical solutions can only exist up to the first time to at which any two characteristics intersect. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, RenS de Borst and Thomas J.R. Hughes. © 2004 John Wiley h Sons, Ltd. ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 1.2. Weak solutions Since, in general, classical solutions only exist for a finite time t , it is necessary to introduce the notion of weak solutions that are well-defined for times t > to- Definition 1.1 (Weak solution) Let uo € L°°(E d ). Then, u is a weak solution of (la-lb) ifuQ I/°°(E d x M + ) and (la-lb) hold in the distributional sense, i.e. J J (ud t (f> + /(u) • V(/>) dt dx + J u 4>(x, 0)dx=0 for all <j> € Ctf (R d x E+ ) . (3) R d R+ R d Note that classical solutions are weak solutions and weak solutions that lie in C 1 (E d x 1K + ) satisfy (la-lb) in the classical sense. It can be shown (see Kruzkov, 1970; Oleinik, 1963) that there always exists at least one weak solution to (la-lb) if the flux function / is at least Lipschitz continuous. Nevertheless, the class of weak solutions is too large to ensure uniqueness of solutions. An important class of solutions are piecewise classical solutions with discontinuities separating the smooth regions. The following lemma gives a necessary and sufficient condition imposed on these discontinuities such that the solution is a weak solution (see for example Godlewski and Raviart, 1991; Kroner, 1997). Later a simple example is given where infinitely many weak solutions exist. Lemma 1.2 (Rankine-Hugoniot jump condition) Assume that R d x R + is separated by a smooth hyper surface S into two parts fij and Cl r . Furthermore, assume u is a C 1 -function on fit and Q r , respectively. Then, u is a weak solution of (la-lb) if and only if the following two conditions hold: i) u is a classical solution in Cli and Cl T , respectively. ii) u satisfies the Rankine-Hugoniot jump condition, i.e. [u]s = [f(u)] -v onS . (4) Here, (v, —s) T denotes a unit normal vector for the hypersurface S and [w] denotes the jump in w across the hypersurface S. In one space dimension, it may be assumed that S is parameterized by (cr(t),t) such that s = <r'(t) and u = 1. The Rankine-Hugoniot jump condition then reduces to . = m 0«S. (5) [u] Example 1.3 (Non-uniqueness of weak solutions) Consider the one- dimensional Burg- ers' equation, f(u) = u 2 /2, with Riemann data: Uo(x) = ui for x < and uo{x) = u T for x > 0. Then, for any a > max(ui, —u T ) a function u given by u(x,t) — < Ui, x < Sit —a, sit < x < , fi . a, < x < s 2 t ^ ' u r , s^t < x is a weak solution if s± = (ui — a)/2 and s 2 = (a + u T )/2. This is easily checked since u is piecewise constant and satisfies the Rankine-Hugoniot jump condition. This elucidates a one- parameter family of weak solutions. In fact, there is also a classical solution whenever ui < u T . Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 5 In this case, the characteristics do not intersect and the method of characteristics yields the classical solution {U[, x < Ult x/t, mt < x < u T t . (7) u r , u r t < x This solution is the unique classical solution but not the unique weak solution. Consequently, additional conditions must be introduced in order to single out one solution within the class of weak solutions. These additional conditions give rise to the notion of a unique entropy weak solution. 1.3. Entropy weak solutions and vanishing viscosity In order to introduce the notion of entropy weak solutions, it is useful to first demonstrate that there is a class of additional conservation laws for any classical solution of (la). Let u be a classical solution and n : M -> E a smooth function. Multiplying (la) by tj'(u), one obtains = r)'{u)d t u + r)'(u)V ■ f(u) = d t r](u) + V • F{u) (8) where F is any primitive of n'f. This reveals that for a classical solution u, the quantity rj(u), henceforth called an entropy function, is a conserved quantity. Definition 1.4 (Entropy - entropy flux pair) Lett] : R — ► K be a smooth convex function and F:I-4la smooth function such that F' = r,'f (9) in (8). Then (n, F) is called an entropy - entropy flux pair or more simply an entropy pair for the equation (la). Note 1.5 (Kruzkov entropies) The family of smooth convex entropies n may be equivalently replaced by the non-smooth family of so-called Kruzkov entropies, i.e. n K {u) = \u - k\ for all /tSl (see Kroner, 1997). Unfortunately, the relation (8) can not be fulfilled for weak solutions in general, as it would lead to additional jump conditions which would contradict the Rankine-Hugoniot jump condition lemma. Rather, a weak solution may satisfy the relation (8) in the distributional sense with inequality. To see that this concept of entropy effectively selects a unique, physically relevant solution among all weak solutions, consider the viscosity perturbed equation a t u, + V-/K) = eAu £ (10) with e > 0. For this parabolic problem, it may be assumed that a unique smooth solution u € exists. Multiplying by rj and rearranging terms yields the additional equation d t n{u t ) + V • F(u ( ) = eAry(u £ ) - e7/"K)|Vw| 2 . Furthermore, since n is assumed convex (n" > 0), the following inequality is obtained d t n(u e ) + V • F(u e ) < eAn(u e ) . Taking the limit e — ► establishes (see Malek, Necas, Rokyta and Ruzicka, 1996) that u e converges towards some u a.e. in M d x R + where u is a weak solution of (la-lb) and satisfies the entropy condition d t r)(u) + V-F(u)<0 (11) Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Ren6 de Borst and Thomas J.R. Hughes. © 2004 John Wiley &: Sons, Ltd. 6 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS in the sense of distributions on R d x R+ . By this procedure, a unique weak solution has been identified as the limit of the approximating sequence u e . The obtained solution u is called the vanishing viscosity weak solution of (la-lb). Motivated by the entropy inequality (11) of the vanishing viscosity solution, it is now possible to introduce the notion of entropy weak solutions. This notion is weak enough for the existence and strong enough for the uniqueness of solutions to (la-lb). Definition 1.6 (Entropy weak solution) Let u be a weak solution of (la-lb). Then, u is called an entropy weak solution if u satisfies for all entropy pairs (n, F) f f (n(u)d t <t> + F(u) ■ V<p) dtdx + f r)(u )4>(x,0) dx>0 for all <f> € C^(R d x E+,E + ) . (12) R<fR + R d Prom the vanishing viscosity method, it is known that entropy weak solutions exist. The following L 1 contraction principle guarantees that entropy solutions are uniquely defined (see Kruzkov, 1970). Theorem 1.7 (L 1 contraction principle) Let u and v be two entropy weak solutions of (la-lb) with respect to initial data uo and v<$. Then, the following L 1 contraction principle holds IK->*) -*>(•> *)IIli(r") < IK-follii(R<') (13) for almost every t > 0. This principle demonstrates a continuous dependence of the solution on the initial data and consequently the uniqueness of entropy weak solutions. Finally, note that an analog of the Rankine-Hugoniot condition exists (with inequality) in terms of the entropy pair for all entropy weak solutions [r)(u)] s > [F(u)\ -v on S . (14) 1.4- Measure-valued or entropy process solutions The numerical analysis of conservation laws requires an even weaker formulation of solutions to (la-lb). For instance, the convergence analysis of finite volume schemes makes it necessary to introduce so called measure-valued or entropy process solutions (see DiPerna, 1985; Eymard, Galluoet and Herbin, 2000). Definition 1.8 (Entropy process solution) A function /j,(x,t,a) £ I/°°(R d x E + x (0,1)) is called an entropy process solution of (la-lb) if u satisfies for all entropy pairs (t],F) i J J J T](fi)d t {4> + F{/i) ■ V</>) da dtdx+ J n(u )<p{x, 0) dx > for all <p € Cl (R d xl + ,l+) . R<IR+ U d The most important property of such entropy process solutions is the following uniqueness and regularity result (see Eymard, Galluoet and Herbin, 2000 [Theorem 6.3]). Theorem 1.9 (Uniqueness of entropy process solutions) Let uq £ I/°°(E d ) and f € C 1 (K). The entropy process solution fi of problem (la-lb) is unique. Moreover, there exists a function u € L°°(E' i x E + ) such that u{x,t) = fj,(x,t,a) a.e. for {x,t,a) e R d x E+ x (0,1) and u is the unique entropy weak solution of (la- lb). Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 7 2. Finite volume (FV) methods for nonlinear conservation laws In the finite volume method, the computational domain, Q C M. d , is first tessellated into a collection of non overlapping control volumes that completely cover the domain. Notationally, let T denote a tessellation of the domain ft with control volumes T € T such that UxeTT = Ci. Let hr denote a length scale associated with each control volume T, e.g. hr = diam(T). For two distinct control volumes T* and Tj in T, the intersection is either an oriented edge (2-D) or face (3-D) e^ with oriented normal i/;j or else a set of measure at most d — 2. In each control volume, an integral conservation law statement is then imposed. Definition 2.1 (Integral conservation law) An integral conservation law asserts that the rate of change of the total amount of a substance with density u in a fixed control volume T is equal to the flux f of the substance through the boundary dT d_ dt I udx+ f f(u)-dv = . (15) JT JdT This integral conservation law statement is readily obtained upon spatial integration of the divergence equation (la) in the region T and application of the divergence theorem. The choice of control volume tessellation is flexible in the finite volume method. For example, Fig. • storage location 4fo control volume a. Cell-centered b. Vertex-centered Figure 1. Control volume variants used in the finite volume method: (a) cell-centered and (b) vertex-centered control volume tessellation. 1 depicts a 2-D triangle complex and two typical control volume tessellations (among many others) used in the finite volume method. In the cell-centered finite volume method shown in Fig. la, the triangles themselves serve as control volumes with solution unknowns stored on a per triangle basis. In the vertex-centered finite volume method shown in Fig. lb, control volumes are formed as a geometric dual to the triangle complex and solution unknowns stored on a per triangulation vertex basis. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 8 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 2.1. Godunov finite volume discretizations Fundamental to finite volume methods is the introduction of the control volume cell average for each Tj £ T i r u dx . (16) -f 3 ~ \T } For stationary meshes, the finite volume method can be interpreted as producing an evolution equation for cell averages / udx^lTjl^uj . (17) dt 'dt u > Godunov, 1959 pursued this interpretation in the discretization of the gas dynamic equations by assuming piecewise constant solution representations in each control volume with value equal to the cell average. However, the use of piecewise constant representations renders the numerical solution multivalued at control volume interfaces thereby making the calculation of a single solution flux at these interfaces ambiguous. The second aspect of Godunov's scheme and subsequent variants was the idea of supplanting the true flux at interfaces by a numerical flux function, g(u,v) :lxR4l, a Lipschitz continuous function of the two interface states u and v. A single unique numerical flux was then calculated from an exact or approximate local solution of the Riemann problem in gas dynamics posed at these interfaces. Figure 2 depicts a representative 1-D solution profile in Godunov's method. For a given control volume Tj = [xj-i/2,Xj+i/2], Riemann problems are solved at each interface Xj±i/ 2 - For example, at the interface Xj +1 / 2 the Riemann problem counterpart of (la-lb) d T W j+1/2 {Z,T) + c /(u; i+1/2 (£,T)) - in K. x R+ for u>j+i/2 (£>?■) € IK with initial data Wj+i/ 2 (£,0) = = / U; \ u 3+l if£<0 if£>0 is solved either exactly or approximately. From this local solution, a single unique numerical flux at Xj + i/ 2 is computed from g(uj,Uj + i) = /(w., +1 /2(0,R + )). This construction utilizes the fact that the solution of the Riemann problem at £ = is a constant for all time r > 0. U(x,t) Xj-l/2 Xj + l/2 Xj +3 Q Figure 2. 1-D control volume, Tj = [xj-i/2,Xj+i/2], depicting Godunov's interface Riemann problems, Wj±i/2 (£; T )> from piecewise constant interface states. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 9 In higher space dimensions, the flux integral appearing in (15) is similarly approximated by /. f{u)-dvx, ^ 9jk{uj,u k ) (18) dTj Ve jk €dTj where the numerical flux is assumed to satisfy the properties: • (Conservation) This property insures that fluxes from adjacent control volumes sharing a mutual interface exactly cancel when summed. This is achieved if the numerical flux satisfies the identity 9jk(u,v) = -gkj(v,u) . (19a) • (Consistency) Consistency is obtained if the numerical flux with identical state arguments reduces to the true flux of that same state, i.e. g jk (u,u)= f f(u)-dv . (19b) Combining (17) and (18) yields perhaps the simplest finite volume scheme in semi-discrete form. Let V£ denote the space of piecewise constants, i.e. V° = {v\v\ T e x (T), VTeT} (20) with x(T) a characteristic function in the control volume T. Definition 2.2 (Semi-discrete finite volume method) The semi-discrete finite volume approximation of (la-lb) utilizing continuous in time solution representation, t € [0, t], and piecewise constant solution representation in space, Uh(t) € V^, such that u j(*) = W1 u h(x,t) dx with initial data %(0) = WH u o(x) dx and numerical flux function gjk(uj,Uk) is given by the following system of ordinary differential equations u i + vf~\ Yl 9jk(uj,u k ) = , VTj &T . (21) dt ' IT,-, This system of ordinary differential equations can be marched forward using a variety of explicit and implicit time integration methods. In Sect. 4.1, time integration schemes that preserve properties of the spatial discretization are considered in more detail. Let u" denote a numerical approximation of the cell average solution in the control volume Tj at time t n = nAi. A particularly simple time integration method is the forward Euler scheme d < +1 - «7 dt 3 At thus producing the fully-discrete finite volume form. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 10 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Definition 2.3 (Fully-discrete finite volume method) The fully-discrete finite volume approximation of (la- lb) for the time slab interval [t n ,t n + At] utilizing the piecewise constant solution representation in space, u% € Vfi, such that u» = -L f <(x) dx with initial data u j °- ■ - ' -'^ dx and numerical flux function ffjfc(uj ,u£) is given by the following fully-discrete system U " +1=U "~#I £ #*(«>*)> VTjZT. (22) In subsequent sections, properties of the semi-discrete scheme (21) and fully-discrete scheme (22) will be examined in more detail. 2.1.1. Monotone schemes. Unfortunately, the numerical flux conditions (19a) and (19b) are insufficient to guarantee convergence to entropy satisfying weak solutions (12) and additional numerical flux restrictions are necessary. Two classes of numerical fluxes that guarantee such convergence for piecewise constant numerical solution data are monotone fluxes and E-fluxes. Specifically, Harten, Hyman and Lax, 1976 provide the following result concerning convergence of the fully-discrete one-dimensional scheme to weak solutions which was later generalized to (22) and irregular grids by Cockburn, Coquel and Lefloch, 1994. Theorem 2.4 (Monotone schemes and weak solutions) Consider a 1-D finite volume discretization of (la-lb) with 2k + 1 stencil on a uniformly spaced mesh in both time and space with corresponding mesh spacing parameters At and Ax u? +1 = H (u i+k ,...,u h ...,Uj-k) = u j - -^(9j+i/2 ~ 93-1/2) (23) and consistent numerical flux of the form Sj+1/2 - 9(uj+k, ■■■, Uj+i,Uj, ..., Uj- k+ i) that is monotone in the sense apr 3 > , V \l\ < k . (24) duj+i Then as At and Ax tend to zero with At/ Ax — constant, u" converges boundedly almost everywhere to u(x,t), an entropy satisfying weak solution of (la-lb). The monotonicity condition (24) motivates the introduction of Lipschitz continuous monotone fluxes satisfying d 9j+l/2 dui dgj+1/2 dui > if l = j (25a) < if l^j (25b) Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 11 together with a CFL (Courant-Priedrichs-Levy) like condition 1 _ At f dg j+1/2 _ dg^ 1/2 \ > Q Ax \ dv,j duj J ~ so that (24) is satisfied. Some examples of monotone fluxes for (la) include • (Godunov flux) min f(u) if Uj < Uj+i 9 j+ i/2 ~ \ • (Lax-Friedrichs flux) 9 i+ 1 /* - \ max /(«) if « 3 - > u i+1 ^ ; g£ l/2 = 5 (/(«,•) + /(«i+i)) - J sup |/'(«)| (u i+1 - «,-) . (27) 2.L2. E-flux schemes. A more general class of numerical fluxes was introduced and analyzed by Osher, 1984 that still guarantees convergence to weak entropy solutions when used in (22) or (23). These fluxes are called E-fluxes, ffj+1/2 = 9 E (uj+k, ■ ■ .,Uj+i,uj, . . .,Uj- k +i), due to the relationship to Olienick's well-known E-condition which characterizes entropy satisfying discontinuities. E-fluxes satisfy the inequality ^ +1/2 " /(M) < , \/u e [uj,u j+1 ] . (28) U j+ i - Uj E-fluxes can be characterized by their relationship to Godunov's flux. Specifically, E-fluxes are precisely those fluxes such that gf+i/2 ^ 9°+i/2 if u i+i < u i ( 29a ) Sj+1/2 > Pj+1/2 if u i+i > u i ■ ( 29b ) Viewed another way, note that any numerical flux can be written in the form 9j+i/2 = j (/(«>) + f( u J+i)) - 2 < 5(%+ 1 _ U A ^ where Q(-) denotes a viscosity for the scheme. When written in this form, E-fluxes are those fluxes that contribute at least as much viscosity as Godunov's flux, i.e. Qf +1/2 < Q j+1/2 . (31) The most prominent E-flux is the Enquist-Osher flux gf+i/2 = \ (/KO + /(«;+i)) ~ \ £ +1 l/'( s )l ds > (32) although other fluxes such as certain forms of Roe's flux with entropy fix fall into this category. From (29a-29b), the monotone fluxes of Godunov gf +1 / 2 and Lax-Friedrichs g^/ 2 are also E-fluxes. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 12 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 2.2. Stability, convergence and error estimates Several stability results are provided here that originate from discrete maximum principle analysis and are straightforwardly stated in multi-dimensions and on general unstructured meshes. In presenting results concerning convergence and error estimates, a notable difference arises between one and several space dimensions. This is due to the lack of a BV bound on the approximate solution in multi-dimensions. Thus, before considering convergence and error estimates for finite volume methods, stability results are presented first together with some a priori bounds on the approximate solution. 2.2.1. Discrete maximum principles and stability. A compelling motivation for the use of monotone and E-fluxes in the finite volume schemes (21) and (22) is the obtention of discrete maximum principles in the resulting numerical solution of nonlinear conservation laws (la). A standard analysis technique is to first construct local discrete maximum principles which can then be applied successively to obtain global maximum principles and stability results. The first result concerns the boundedness of local extrema in time for semi-discrete finite volume schemes that can be written in nonnegative coefficient form. Theorem 2.5 (LED Property) The semi-discrete scheme for each Tj € T ir = W~\ 51 C jk (u h )(u k -Uj), (33) Ve ih £dTj dt |:r,-j is Local Extremum Diminishing (LED), i.e. local maxima are non-increasing and local minima are nondecr easing, if C jk (u h )>0 , Ve jk edTj . (34) Rewriting the semi-discrete finite volume scheme (21) in terms of monotone fluxes or E-fiuxes duj 1 ^ 9jk{uj,u k ) - f{uj) -i>j k at \1A *— ' Uk — Ui = ~W\ ^2 -KTiuJ'UjkXuk-Uj) (35) |J ^Ve ifc e9T, ° Uk for appropriately chosen iij k 6 [uj,u k ] together with the monotone flux conditions (25a-25b) or the E-flux condition (28) reveals that monotone flux and E-flux finite volume schemes (21) are LED. In order to obtain local space-time maximum principle results for the fully- discrete discretization (22) requires the introduction of an additional CFL-like condition for non-negativity of coefficients in space-time. Theorem 2.6 (Local space-time discrete maximum principle) The fully-discrete scheme for the time slab increment [t n ,t n+i ] and each Tj € T u™ +1 = u 1 + Wa E C;*«)(«*-«") (36) exhibits a local space-time discrete maximum principle min ( U J, «?) < U7+ 1 < max («J, «?) (37) Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 13 if C jk (<)>0, Ve jk edTj (38) and satisfies the CFL-like condition 1 -£r-\ E C i*«)>°- (39) Again noting that the flux terms in the fully-discrete finite volume scheme (22) can be written in the form (35) reveals that the monotone flux conditions (25a-25b) or the E-flux condition (28) together with a local CFL-like condition obtained from (39) imply a local space-time discrete maximum principle. By successive application of Theorem 2.6, a global L°°-stability bound is obtained for the scalar initial value problem (la-lb) in terms of initial data uo(x). Theorem 2.7 (L°°-stability) Assume a fully-discrete finite volume scheme (22) for the scalar initial value problem (la-lb) utilizing monotone fluxes or E-fluxes that satisfy a local CFL-like condition as given in Theorem 2.6 for each time slab increment [t n ,t n+1 ]. Under these conditions, the finite volume scheme is L°° -stable and the following estimate holds: inf u {x) < u? < sup u (x), for all (Tj,t n ) G T x [0,r]. (40) Consider now steady-state solutions, u n+1 =u n = u*, using monotone flux or E-flux schemes in the fully-discrete finite volume scheme (22). At steady state, non-negativity of the coefficients C{uh) in (36) implies a discrete maximum principle. Theorem 2.8 (Local discrete maximum principle in space) The fully-discrete scheme (36) exhibits a local discrete maximum principle at steady state, u* h , for each Tj £ T min ut < ul < max ut (41) if C jk (u* h )>0 , \/e jk edTj . Once again by virtue of (25a-25b) and (28), the conditions for a local discrete maximum principle at steady state are fulfilled by monotone flux and E-flux finite volume schemes (22) . Global maximum principles for characteristic boundary valued problems are readily obtained by successive application of the local maximum principle result. The local maximum principles given in (37) and (41) preclude the introduction of spurious extrema and C(l) Gibbs-like oscillations that occur near solution discontinuities computed using many numerical methods (even in the presence of grid refinement). For this reason, discrete maximum principles of this type are a highly sought after design principle in the development of numerical schemes for nonlinear conservation laws. 2.2.2. Convergence results. The Instability bound (40) is an essential ingredient in the proof of convergence of the fully-discrete finite volume scheme (22). This bound permits the subtraction of a subsequence that converges against some limit in the L°° weak-* sense. The primary task that then remains is to identify this limit with the unique solution of the problem. So although Z°°-stabiIity is enough to ascertain convergence of the scheme, stronger estimates are needed in order to derive convergence rates. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 14 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Let BV denote the space of functions with bounded variation, i.e. BV = {g G L 1 (R d ) | \g\BV < °o} with \g\ B v = sup / g V • <p dx . riwiioo<i From the theory of scalar conservation laws, it is known that, provided the initial data is in BV, the solution remains in BV for all times. Therefore, it is desirable to have an analog of this property for the approximate solution as well. Unfortunately, up to now, such an result is only rigorously proved in the one-dimensional case or in the case of tensor product cartesian meshes in multiple space dimensions. In the general multi-dimensional case, the approximate solution can only be shown to fulfill some weaker estimate which is thus called a weak BV estimate (see Vila, 1994; Cockburn, Coquel and Lefioch, 1994; Eymard, Gallouet, Ghilani and Herbin, 1998). Theorem 2.9 (Weak BV estimate) Let T be a regular triangulation, and let J be a uniform partition of [0, r], e.g. At n = At. Assume that there exists some a > such that ah 2 < \Tk\, a \dTk\ < h. For the time step At", assume the following CFL-like condition for a given £ G (0, 1) At" < L g where L g is the Lipschitz constant of the numerical flux function. Furthermore, let uo G L°°(]R d ) n BV(M. d ) l~l L 2 (M d ). Then, the numerical solution of the fully- discrete discretization (22) fulfills the following estimate Y^At^XjihW-uflQuiu^u?) < KVT\B R+h {0)\Vh , (42) n jl where K only depends on a, L g , £ and the initial function uq. In this formula Qji is defined as n , % _ 2 9ji (u,v)~ gji (u,u)- g jt (v, v) Q jt (u,v) = — y and Xji denotes the discrete cutoff function on Br(0) C R d , i.e. _ / 1, if (T j UT i )n5fi(O)^0 Xjl \ 0, else Note that in the case of a strong BV estimate, the right-hand side of (42) would be 0(h) instead of C(Vh). Another important property of monotone finite volume schemes is that they preserve the I/ 1 -contraction property (see Theorem 1.7). Theorem 2.10 (L 1 — contraction property and Lipschitz estimate in time) Letuh,Vh G Vfi be the approximate monotone finite volume solutions corresponding to initial data uo, vq assuming that the CFL-like condition for stability has been fulfilled. Then the following discrete L 1 -contraction property holds \\u h (;t + T) -w /l (-,t-l-r)|| I ,i( R <!) < \\u h (;i) - v h (-, t)||ii( R <«) . Furthermore, a discrete Lipschitz estimate in time is obtained £ |T,|K +1 - U ?| < L g At n £ £ le^lK - u?| . Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 15 The principle ingredients of the convergence theory for scalar nonlinear conservation laws are compactness of the family of approximate solutions and the passage to the limit within the entropy inequality (12). In dealing with nonlinear equations, strong compactness is needed in order to pass to the limit in (12). In one space dimension, due to the BV estimate and the selection principle of Helly, strong compactness is ensured and the passage to the limit is summarized in the well known Lax-Weiidroff theorem (see Lax and Wendroff, 1960). Theorem 2.11 ( Lax-Wendroff theorem) Let (u m ) m€ N be a sequence of discrete solutions defined by the finite volume scheme in one space dimension with respect to initial data uq. Assume that (u m )meN is uniformly bounded with respect to m in L°° and u m converges almost everywhere in E x 1 + against some function u. Then u is the uniquely defined entropy weak solution of (la- lb). With the lack of a BV estimate for the approximate solution in multiple space dimensions, one cannot expect a passage to the limit of the nonlinear terms in the entropy inequality in the classical sense, i.e. the limit of u m will not in general be a weak solution. Nevertheless, the weak compactness obtained by the jL°°-estimate is enough to obtain a measure-valued or entropy process solution in the limit. The key theorem for this convergence result is the following compactness theorem of Tartar (see Tartar, 1983; Eymard, Galluoet and Herbin, 2000). Theorem 2.12 (Tartar's Theorem) Let (u m )meN be a family of bounded functions in Z,°°(K n ). Then, there exists a subsequence (u m ) m€ N, and a function u G L°°(R n x (0,1)) such that for all functions g G C(R) the weak-* limit of g{u m ) exists and lim / g(u m {x))<j>(x)dx = / / g(u(x,a))<f>{x)dx da, for all <f> G L 1 (R n ) . (43) In order to prove the convergence of a finite volume method, it now remains to show that the residual of the entropy inequality (12) for the approximate solution it/, tends to zero if h and At tend to zero. Before presenting this estimate for the finite volume approximation, a general convergence theorem is given which can be viewed as a generalization of the classical Lax-Wendroff result (see Eymard, Galluoet and Herbin, 2000) . Theorem 2.13 (Sufficient condition for convergence) Let uo G L co (R d ) and f G C l (R). Further, let (u m ) m£ N be any family of uniformly bounded functions in L co (R d x R + ) that satisfies the following estimate for the residual of the entropy inequality using the class of Kruzkov entropies n K (see Note 1.5). I j (VK{u m )d t (p + F T)K (u rn ) -V<p) dtdx+ / n K (u )</>(x,0)dx > -R(K,u m ,4>) (44) for all k G R and <f> G Cg (R d x E + , R + ) where the residual R(k, u m , 4>) tends to zero for m — >■ co uniformly in k. Then, u m converges strongly to the unique entropy weak solution of (la-lb) in Lf oc (R d x R+) for all pe [l,oo) . Theorem 2.14 (Estimate on the residual of the entropy inequality) Let (u m ) m eN be a sequence of monotone finite volume approximations satisfying a local CFL-like condition as given in (39) such that h,At tend to zero for m — ¥ oo. Then, there exist measures Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 16 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Hm £ M(R d xR + ) andu m £ M(R d ) such that the residual R(K,u m ,cf>) of the entropy inequality is estimated by R{ K ,u m ,<P)< [ [ (\d t <t>(x,t)\ + \V<P(x,t)\)dn m (x,t)+ [ 4>(x,0)dvm(x) for all k £ R and 4> € Co(R d x R + ,1R + ). The measures /x TO and u m satisfy the following properties: 1. For all compact subsets fi CC R d x E + , lim m ^. 0O /i m (0) = 0. 2. For all g £ Co (R d ) the measure v m is given by {u m ,g) = J^ d g{x) \uq {x) — u m (x, 0) | dx. These theorems are sufficient for establishing convergence of monotone finite volume schemes. Corollary 2.15 (Convergence theorem) Let (u m ) m gN oe a sequence of monotone finite volume approximations satisfying the assumptions of Theorem 2.14- Then, u m converges strongly to the unique entropy weak solution of (la-lb) in Lif oc (M. d x E + ) for all p £ [l,oo). Convergence of higher order finite volume schemes can also be proven within the given framework as long as they are Z^-stable and allow for an estimate on the entropy residual in the sense of Theorem 2.14, for details see Kroner, Noelle and Rokyta, 1995; Chainais-Hillairet, 2000. 2.2.3. Error estimates and convergence rates. There are two primary approaches taken to obtain error estimates for approximations of scalar nonlinear conservation laws. One approach is based on the ideas of Oleinik and is applicable only in one space dimension (see Oleinik, 1963; Tadmor, 1991). The second approach which is widely used in the numerical analysis of conservation laws is based on the doubling of variables technique of Kruzkov (see Kruzkov, 1970; Kuznetsov, 1976). In essence, this technique enables one to estimate the error between the exact and approximate solution of a conservation law in terms of the entropy residual R(k, u m ,$) introduced in (44). Thus, an a posteriori error estimate is obtained. Using a priori estimates of the approximate solution (see Section 2.2.1, and Theorems 2.9, 2.10), a convergence rate or an a priori error estimate is then obtained. The next theorem gives a fundamental error estimate for conservation laws independent of the particular finite volume scheme (see Eymard, Galluoet and Herbin, 2000; Chainais-Hillairet, 1999; Kroner and Ohlberger, 2000). Theorem 2.16 (Fundamental error estimate) Let uq £ BV(R d ) and let u be an entropy weak solution of (la-lb). Furthermore, let v £ L°°(R d x R + ) be a solution of the following entropy inequalities with residual term R: f ! Vn(v)dt4> + F v Jv)-\/^>+ [ n K (u )<j>(;0)>-R( K ,v,4>) (45) JU d JR+ JR d for all k£R and <f> £ C^(R d x E+ , R + ). Suppose that there exist measures fj, v £ M{R d xR + ) and v v £ M(R d ) such that R(k,v,(J>) can be estimated independently of k by R(k, v, 4>) < (\8 t ct>\ + | V0|, im,) + (\4>{; 0)|, u v ). (46) Let K CC R d x !+, u = Lip(f) , and choose T,R and x such that T e]0, ^[ and K lies within its cone of dependence Do, i.e. K C Do where Dg is given as D 5 := |J B R - Ut+S (x ) x {r}. (47) 0<t<T Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 17 Then, there exists a 8 > and positive constants C\ , Ci such that u, v satisfy the following error estimate ll« ~ v\\ LHK) < T{u v {B R+5 {x )) + CiUviPs) + C 2 V»v{D s )). (48) This estimate can be used either as an a posteriori control of the error, as the right-hand side of the estimate (48) only depends on v, or it can be used as an a priori error bound if one is able to estimate further the measures \x v and u v using some a priori bounds on v. Finally, note that comparable estimates to (48) are obtainable in an L°° (0,T; L 1 (R d ))-noim (see Cockburn and Gau, 1995; Bouchut and Perthame, 1998). 2.2.4- A posteriori error estimate. Theorem 2.17 (A posteriori error estimate) Assume the conditions and notations as in Theorem 2.16. Letv = Uh be a numerical approximation to (la-lb) obtained from a monotone finite volume scheme that satisfies a local CFL-like condition as given in (39). Then the following error estimate holds L \u-u h \ < r(|K-u h (-,0)|| L i (BH+h(lo)) +C 1 r ? + C 2 V57), (49) * K where i =£ £ h n+1 - «"iAt n ^ + 2 j] £ At n (At n + h^QiiOi?,*?)!"" - «n , n€l jeM(t") n€/o (j,l)eE(t n ) n / , , A _ 2 9ji (u, v) - 9ji (u, u) - g jt (v, v) Qji{u,v) = u — v with the index sets I ,M(t),E(t) given by I = {n | 0<t n <mm{^^-,T} , M(t) = {j | there exists x £ Tj such that (x,t) & Dr + s} , E(t) = {{j, I) \ there exists xeTjUTi such that (x,t) G D R+S } . Furthermore, the constants C\,C2 only depend on T, u>, \]uo]\bv and \]uo\\l°° (for details see Kroner and Ohlberger, 2000). Note that this o posteriori error estimate is local, since the error on a compact set K is estimated by discrete quantities that are supported in the cone of dependence Dr+s. 2.2.5. A priori error estimate. Using the weak BV estimate (Theorem 2.9) and the Lipschitz estimate in time (Theorem 2.10), the right-hand side of the a posteriori error estimate (49) can be further estimated. This yields an a priori error estimate as stated in the following theorem (for details see Eymard, Galluoet and Herbin, 2000; Chainais-Hillairet, 1999). Theorem 2.18 (A priori error estimate) Assume the conditions and notations as in Theorem 2.16 and let v = Uh be the approximation to (la), (lb) given by a monotone finite volume scheme that satisfies a local CFL-like condition as given in (39). Then there exists a constant C > such that I \u-u h \ < Ch l/i . Moreover, in the one-dimensional case, the optimal convergence rate of h 1 ' 2 is obtained. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 18 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 2.2.6. Convergence proofs via the streamline diffusion discontinuous Galerkin finite element method. It is straightforward to show that the fully-discrete finite volume scheme (22) can be viewed as a specific case of the more general streamline diffusion discontinuous Galerkin (SD-DG) finite element method which utilizes the mesh dependent broken space V* denned as V£ = {v\v\ T eV p (T),VTeT} (50) with V P {T) the space of polynomials of degree < p in the control volume T. By generalizing the notion of gradient and flux to include the time coordinate as well, the discontinuous Galerkin finite element method for a space-time tessellation T n spanning the time slab increment [i n ,i n+1 ] is given compactly by the following variational statement. SD-DG(p) finite element method. Find Uh 6 Vjf such that Viv, € V£ and n = 0, 1, . . . V) ( I (Vh+ S(u h )f'(u h ) ■ Vv h ) V • f(u h ) dx+ I e(u h ) \>u h ■ Vv h dx T€r „ V Jt Jt + v h -(g(u h -,u h +)-f(u h -)-i/)ds)=0 (51) JdT ) where in the integration over dT it is understood for the portion x £ dT tldT' that Uh-,Vh- denotes the trace restriction of uh(T) and Vh(T) onto dT and Uh+ denotes the trace restriction of Uh,(T') onto dT'. Given this space-time formulation, convergence results for a scalar nonlinear conservation law in multi-dimensions and unstructured meshes are given in Jaffre, Johnson and Szepessy, 1995 for specific choices of the stabilization functions 6(uh) : R i-4 R + and e(u/,) : R !->• K + together with a monotone numerical flux function g{uh-,Uh + ). Using their stabilization functions together with a monotone flux function, the following convergence result is obtained: Theorem 2.19 (SD-DG(p) convergence) Suppose that components of f'(u) € C d (M.) are bounded and that uq € L2(M. d ) has compact support. Then the solution Uh of the SD-DG (p) method converges strongly in L l ° c (M. d x R + ), 1 < p < 2, to the unique solution u of the scalar nonlinear conservation law system (la-lb) as H = max(||/i||x, 00 (R<i'), Ai) tends to zero. The proof of convergence to a unique entropy solution on general meshes for p > is based on an extension by Szepessy, 1989 of a uniqueness result by DiPerna, 1985 by providing convergence for a sequence of approximations satisfying: • a uniform Loo bound in time and Li in space, • entropy consistency and inequality for all Kruzkov entropies, • consistency with initial data. By choosing SD-DG(O), the dependence on the as yet unspecified stabilization functions S(u h ) and e(u/,) vanishes identically and the fully-discrete scheme (22) with monotone flux function is exactly reproduced, thus yielding a convergence proof for general scalar conservation laws for the finite volume method as well. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 19 3. Higher order accurate FV generalizations Although an 0(h}/ 2 ) ii-norm error bound for the monotone and E-flux schemes of Sect. 2 is known to be sharp (Peterson, 1991), an 0(h) solution error is routinely observed in numerical experiments with convex flux functions. Even so, first order accurate schemes are generally considered too inaccurate for most quantitative calculations unless the mesh spacing is made excessively small thus rendering the schemes inefficient. Godunov, 1959 has shown that all linear schemes that preserve solution monotonicity are at most first order accurate. The low order accuracy of these monotonicity preserving linear schemes has motivated the development of higher order accurate schemes with the important distinction that these new schemes utilize essential nonlinearity so that monotone resolution of discontinuities and high order accuracy away from discontinuities are simultaneously attained. 3.1. Higher order accurate FV schemes in 1-D A significant step forward in the generalization of Godunov's finite volume method to higher order accuracy is due to van Leer, 1979. In the context of Lagrangian hydrodynamics with Eulerian remapping, van Leer generalized Godunov's method by employing linear solution reconstruction in each cell (see Fig. 3b). Let N denote the number of control volume cells in a. Cell averaging of quartic data b. Linear reconstruction c. Quadratic reconstruction Figure 3. Piecewise polynomial approximation used in the finite volume method: (a) cell averaging of analytic data, (b) piecewise linear reconstruction from cell averages and (c) piecewise quadratic reconstruction from cell averages. space so that the j'-th cell extends over the interval Tj = [xj_i/ 2 ,Xj+i/2] with length Axj such that 1>i<j<nTj = [0, 1] with T{ f)Tj = 0,i ^ j. In a purely Eulerian setting, the higher order accurate schemes of van Leer are of the form duj ~dt + ^(5(Vl/2' U i"+l/2)-^ U j-l/2'Vl/2)) = where g(u,v) is a numerical flux function utilizing states u. ±1 , 2 and ut ±1 , 2 obtained from evaluation of the linear solution reconstructions from the left and right cells surrounding the Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 20 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS interfaces Xj±x/ 2 , By altering the slope of the linear reconstruction in cells, non-oscillatory resolution of discontinuities can be obtained. Note that although obtaining the exact solution of the scalar nonlinear conservation law with linear initial data is a formidable task, the solution at each cell interface location for small enough time is the same as the solution of the Riemann problem with piecewise constant data equal to the linear solution approximation evaluated at the same interface location. Consequently, the numerical flux functions used in Sect. 2 can be once again used in the generalized schemes of van Leer. This single observation greatly simplifies the construction of higher order accurate generalizations of Godunov's method. The ideas of van Leer have been extended to quadratic approximations in each cell (see Fig. 3c) by Colella and Woodward, 1984. Although these generalizations of Godunov's method and further generalizations given later can be interpreted in 1-D as finite difference discretizations, concepts originally developed in 1-D such as solution monotonicity, positive coefficient discretization and discrete maximum principle analysis are often used in the design of finite volume methods in multiple space dimensions and on unstructured meshes where finite difference discretization is problematic. 3.1.1. TVD schemes. In considering the scalar nonlinear conservation law (la-lb), Lax, 1973 made the following basic observation: "the total increasing and decreasing variations of a differentiable solution between any pair of characteristics are conserved". Furthermore, in the presence of shock wave discontinuities, information is lost and the total variation decreases. For the 1-D nonlinear conservation law with compactly supported (or periodic) solution data u(x,t), integrating along the constant time spatial coordinate at times ij and t 2 yields /OO rCO \du(x,t 2 )\< \du(x,ti)\, t 2 >h . (52) -00 J —00 This motivated Harten, 1983 to consider the discrete total variation TV(ufc) = ]P \A j+1/2 u h \ , A j+1/2 Uh. = u j+ i - Uj j and the discrete total variation non-increasing (TVNI) bound counterpart to (52) TV« +1 ) < TV«) (53) in the design of numerical discretizations for nonlinear conservation laws. A number of simple results relating TVNI schemes and monotone schemes follow from simple analysis. Theorem 3.1 (TVNI and monotone scheme properties, Harten, 1983) (i) Monotone schemes are TVNI. (ii) TVNI schemes are monotonicity preserving, i.e. the number of solution extreme is preserved in time. Property (i) follows from the L\ contraction property of monotone schemes. Property (ii) is readily shown using a proof by contradiction by assuming a TVNI scheme with monotone initial data that produces new solution data at a later time with interior solution extrema present. Using the notion of discrete total variation, Harten, 1983 then constructed sufficient algebraic conditions for achieving the TVNI inequality (53). Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 21 Theorem 3.2 (Harten's explicit TVD criteria) The fully discrete explicit 1-D scheme u» +1 = u] + At (C i+1/2 «) A j+1/2 < + D j+1/2 «) A,_ 1/2 <) , j = l,...,N (54) is total variation non-increasing if for each j Cj+i/2 > , (55a) D j+1/ 2 < , (55b) 1 - At (C,_ 1/2 - D J+1/2 ) > . (55c) Note that although the inequality constraints (55a-55c) in Theorem 3.2 insure that the total variation is non-increasing, these conditions are often referred to as total variation diminishing (TVD) conditions. Also note that inequality (55c) implies a CFL-like time step restriction that may be more restrictive than the time step required for stability of the numerical method. The TVD conditions are easily generalized to wider support stencils written in incremental form, see for example Jameson and Lax, 1986 and their corrected result in Jameson and Lax, 1987. While this simple Euler explicit integration scheme may seem too crude for applications requiring true high order space-time accuracy, special attention and analysis is given to this fully-discrete form because it serves as a fundamental building block for an important class of high order accurate Runge-Kutta time integration techniques discussed in Sect. 4.1 that, by construction, inherit TVD (and later maximum principle) properties of the fully-discrete scheme. Theorem 3.3 (Generalized explicit TVD criteria) The fully discrete explicit 1-D scheme Jfc-i u] +1 =u? + At J2 Cfl 1/2 (u n h )A j+l+1/2 u%, j = l,...,N (56) l=~k with stencil width parameter k is total variation non-increasing if for each j cf + ~% > , (57a) <7j;*> 3 < 0, (57b) C i+i/2 - C f-i/2 > , -k + 1 < I < k - 1, / ? , (57c) l-At(cf} 1/2 -C^ /2 ) > 0. (57d) The extension to implicit methods follows immediately upon rewriting the implicit scheme in terms of the solution spatial increments Aj + i +1 / 2 «ft and imposing sufficient algebraic conditions such that the implicit matrix acting on spatial increments has a nonnegative inverse. Theorem 3.4 (Generalized implicit TVD criteria) The fully discrete implicit 1-D scheme u n+i _ At ]T C<; i/2 « +1 )A i+I+1/2 <« = „«, j = l,...,N (58) l=-k Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Ren6 de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 22 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS with stencil width parameter k is total variation non-increasing if for each j C%\% > , (59a) Cj;*> 2 < , (59b) c i+i/2- C i-i/2 ^ °> -k + l<l<k-l, IjLO . (59c) Theorems 3.3 and 3.4 provide sufficient conditions for non-increasing total variation of explicit (56) or implicit (58) numerical schemes written in incremental form. These incremental forms do not imply discrete conservation unless additional constraints are imposed on the discretizations. A sufficient condition for discrete conservation of the discretizations (56) and (58) is that these discretizations can be written in a finite volume flux balance form fc-i 9j+i/2 -flj-i/2 = ^2 C j+1/2 (u h )A j+l+1/2 u h l=-k where gj±i/2 are the usual numerical flux functions. Section 3.1.2 provides an example of how the discrete TVD conditions and discrete conservation can be simultaneously achieved. A more comprehensive overview of finite volume numerical methods based on TVD constructions can be found the books by Godlewski and Raviart, 1991 and LeVeque, 2002. 3.1.2. MUSCL schemes. A general family of TVD discretizations with 5-point stencil is the Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) discretization of van Leer, 1979; van Leer, 1985. MUSCL schemes utilize a re-parameter family of interpolation formulas with limiter function ^{R) : R •-> M u J+i/2 = u i + -^*(l/fli)A i+1/2 U|, + — ^(-R^A^i/su,, u i"+i/2 = ' u i~ ^-^*(^+i) A J+i/2"ft - -^*(l/ J R,-+i)A i+ 3 / au fc (60) where Rj is a ratio of successive solution increments 3 Aj-.j/aUh The technique of incorporating limiter functions to obtain non-oscillatory resolution of discontinuities and steep gradients dates back to Boris and Book, 1973. For convenience, the interpolation formulas (60) have been written for a uniformly spaced mesh although the extension to irregular mesh spacing is straightforward. The unlimited form of this interpolation is obtained by setting $(R) = 1. In this unlimited case, the truncation error for the conservation law divergence in (la) is given by Truncation Error = -- — -^-(Ail^/fu) . 4 ox s This equation reveals that for re = 1/3, the 1-D MUSCL formula yields an overall spatial discretization with 0(Aa; 3 ) truncation error. Using the MUSCL interpolation formulas given in (60), sufficient conditions for the discrete TVD property are easily obtained. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 23 Theorem 3.5 (MUSCL TVD scheme) The fully discrete 1-D scheme yd j +1/2 «J +1 = «? ■tf-1/2) . i = i,-..,iv Ax- ^^^/'^ vi-i/ 2 ' with monotone Lipschitz continuous numerical flux function 9j+i/2=9{u- +1/2 ,u+ +1/2 ) utilizing the n-parameter family of MUSCL interpolation formulas (60) and (61) is total variation non-increasing if there exists a ^(R) such that Vi? € K < V(R) < 1-k K ,, . 1 + K -(1 + a): 1 and R (62a) (62b) with a € [—2, 2 (1 — k)/(1 + k)] under the time step restriction At 2 -(2 + oc)k dg Axj 1-K du >0 where dg du sup dg + dg ~ selu i-l/2- "j + l/2 1 For accuracy considerations away from extrema, it is desirable that the unlimited form of the discretization is obtained. Consequently, the constraint \f (1) = 1 is also imposed upon the limiter function. This constraint together with the algebraic conditions (62a-b) are readily achieved using the well known MinMod limiter, $ MM , with compression parameter /? determined from the TVD analysis * MM (i?)=max(0,min(i?, / 3)) , /? € [1, (3 - «)/(! - k)] . Table I. Members of the MUSCL TVD family of schemes. Unlimited Scheme Pn Truncation Error 1/3 -1 1/2 Third-Order Fully Upwind Fromm's Low Truncation Error 4 2 3 5 I(A*) 2 ^/( U ) ±(Ax)^f(u) -^r(Az) 2 ^/M Table I summarizes the MUSCL scheme and maximum compression parameter for a number of familiar discretizations. Another limiter due to van Leer that meets the technical conditions of Theorem 3.5 and also satisfies \&(1) = 1 is given by R + \R\ * VL (i?) 1 + \R\ Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 24 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS This limiter exhibits differentiability away from R = which improves the iterative convergence to steady state for many algorithms. Numerous other limiter functions are considered and analyzed in Sweby, 1984. Unfortunately, TVD schemes locally degenerate to piecewise constant approximations at smooth extrema which locally degrades the accuracy. This is an unavoidable consequence of the strict TVD condition. Theorem 3.6 (TVD critical point accuracy, Osher, 1984) The TVD discretizations (54), (56) and (58) all reduce to at most first order accuracy at non-sonic critical points, i.e. points u* at which f'(u*) ^ and u* = 0. 3.1.3. ENO/WENO schemes. To circumvent the degradation in accuracy of TVD schemes at critical points, weaker constraints on the solution total variation were devised. To this end, Harten proposed the following abstract framework for generalized Godunov schemes in operator composition form (see Harten et al., 1986; Harten et al., 1987; Harten, 1989) ul+ l =A-E{r).Rl(-,ul) . (63) In this equation, u£ € V£ denotes the global space of piecewise constant cell averages as denned in (20), Rp(x) is a reconstruction operator which produces a cell- wise discontinuous p-th order polynomial approximation from the given solution cell averages, E(r) is the evolution operator for the PDE (including boundary conditions), and A is the cell averaging operator. Since A is a nonnegative operator and E(t) represents exact evolution in the small, the control of solution oscillations and Gibbs-like phenomena is linked directly to oscillation properties of the reconstruction operator, R°{x). One has formally in one space dimension TV« +1 ) = TV(^ ■ E{r) ■ R° p (-; u n h )) < TV(R° p (x; u n h )) so that the total variation depends entirely upon properties of the reconstruction operator Rp(x;u^). The requirements of high order accuracy for smooth solutions and discrete conservation give rise to the following additional design criterion for the reconstruction operator • R° p {x; u h ) = u(x) + e(x) Ax p+l + 0(Ax p+2 ) whenever u is smooth (64a) • AlTjRpix; Uh) = Uh\Tj =uj, j = 1, ... ,N to insure discrete conservation (64b) • TV(R(x; u£)) < TV(u£) + 0(Ax p+1 ) an essentially non-oscillatory reconstruction. (64c) Note that e(x) may not be Lipschitz continuous at certain points so that the cumulative error in the scheme is 0(Ax p ) in a maximum norm but remains 0(Ax p+1 ) in an Li-norm. To achieve the requirements of (64a-64c) , Harten and coworkers considered breaking the task into two parts • Polynomial reconstruction from a given stencil of cell averages • Construction of a "smoothest" polynomial approximation by a solution adaptive stencil selection algorithm. In the next section, a commonly used reconstruction technique from cell averages is considered. This is then followed by a description of the solution adaptive stencil algorithm proposed by Harten et al., 1986. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 25 3.1.4. Reconstruction via primitive function. Given cell averages Uj of a piecewise smooth function u(x), one can inexpensively evaluate pointwise values of the primitive function U(x) c)= u Jxo U(x)= / «(0df by exploiting the relationship ^2 AxjUj = U(x j+1/2 ) 3=30 Let H p (x;u) denote a p-th order piecewise polynomial interpolant of a function u. Since u{x) = faU{x) , an interpolant of the primitive function given pointwise samples U(xj + i/ 2 ) yields a reconstruction operator 3 R° p (x;u h ) = — H p+1 (x;U) . As a polynomial approximation problem, whenever U(x) is smooth one obtains ^H p (x; U) = j^U{x) + 0{Ax p+1 ~ k ) ,0<k<p and consequently <* r, , ^ dl —R p{x ., Uh ) = — l , By virtue of the use of the primitive function U(x), it follows that A\ Tj R p (x;u h ) =Uj and from the polynomial interpolation problem for smooth data R° p (x;u h )=u{x) + 0(Ax p+1 ) as desired. 3.1.5. ENO reconstruction. The reconstruction technique outlined in Section 3.1.4 does not satisfy the oscillation requirement given in (64c). This motivated Harten and coworkers to consider a new algorithm for essentially non-oscillatory (ENO) piecewise polynomial interpolation. When combined with the reconstruction technique of Section 3.1.4, the resulting reconstruction then satisfies (64a-c). Specifically, a new interpolant H p (x; u) is constructed so that when applied to piecewise smooth data v(x) gives high order accuracy r R p {x;u h ) = -zjuix) + 0(Ax? +1 - 1 ) H p (x; v) = -^v{x) + 0(Ax p+l - k ) ,0<k<p —H p (x;v) = — k but avoids having Gibbs oscillations at discontinuities in the sense TV(H p {x;v)) < TV(v) + 0(Ax p+l ) . Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 26 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS The strategy pursued by Harten and coworkers was to construct such an ENO polynomial H p (x;w) using the following steps. Define H™° (x; w) = lf™ /2 (x; w) for Xj < x < x j+1 , j = l,...,N where -PjfS^ / 2 is the p-th degree polynomial which interpolates w(x) at the p + 1 successive points {xi}, i p (j) <i< ip{j) + P that include Xj and Xj+i, i.e. Ifj+i/2^i;w) = w(xi) , i p (j) < i < i p (j) +P , l-p< ip{j) - j < . (65) Equation (65) describes p possible polynomials depending on the choice of i p (j) for an interval (xj,xj + i). The ENO strategy selects the value i p (j) for each interval that produces the "smoothest" polynomial interpolant for a given input data. More precisely, information about smoothness of w(x) is extracted from a table of divided differences of w(x) defined recursively for i — 1,. . .,N by w[xi] = w(xi) w[x i+ i] - w[Xi] w[xi,x i+ i Xi+1 - Xi W[Xi,...,X i+k ] = — ' '- L - : i . ^i+A ~ x i The stencil producing the smoothest interpolant is then chosen hierarchically by setting and for 1 < k < p — 1 iu (i\ - i ik ^ ~ l if H^W-i' • • ■'M\(j)+tll < M*i»Cj)> ■ • • > w N fc O)+fc+i]l (ao^ i k +iK3)-<y h ^ otherwise. w Harten et al., 1986 demonstrate the following properties of this ENO interpolation strategy • The accuracy condition ^ N + °/2(z) = "to + 0{Ax*>+ 1 ) , x € (x jy x j+1 ) . • P™°(x) is monotone in any cell interval containing a discontinuity. • There exists a function z(x) nearby P p NO (x) in the interval (xj,Xj+i) in the sense z(x) = P p E f + ° 1/2 (x) + O(Ax^) , x e {x j>Xj+1 ) that is total variation bounded, i.e. the nearby function z(x) satisfies TV{z) < TV(w) . 3.1.6. WENO reconstruction. The solution adaptive nature of the ENO stencil selection algorithm (66) yields non-differentiable fluxes that impede convergence to steady state. In addition, the stencil selection algorithm chooses only one of p possible stencils and other slightly less smooth stencils may give similar accuracy. When w(x) is smooth, using a linear Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 27 combination of all p stencils with optimized weights yields a more accurate C(Ax 2p_1 ) interpolant. More specifically, let -P., ; +1 / 2 denote the unique polynomial interpolating p + 1 points with stencil {xj + i- p +k,Xj+i+k} then P Pij+l/2 {w{x)) = ^^jW + ©(A* 2 *- 1 ) , J>* = 1 . fc=0 *=0 For example, optimized weights for p = 1,2,3 yielding 0(Ax 2p ~ 1 ) accuracy are readily computed p = 1 : wo = 1, 2 1 p=2: wo = -, wi = -, In the WENO schemes of Jiang and Shu, 1996; Shu, 1999, approximate weights, C, are devised such that for smooth solutions Z k = u k + 0(Ax p - 1 ) so that the 0(Ax 2p ~ 1 ) accuracy is still retained using these approximations P pJ+ i/2(w{x)) =E3»ft /a (x) + OiAx 2 '- 1 ) , £> = 1 . fc=0 fc=0 The approximate weights are constructed using the ad hoc formulas ak = 7 — i oTi ' w * — t-m>-i (e + /3*) 2 E*=S°* where e is an approximation to the square root of the machine precision and /3fc is a smoothness indicator 'd l PHx) fc_l _ / ./ _!., , \ 2 v^ f Xi+U2 oi i / trPSix) \ For a sequence of smooth solutions with decreasing smoothness indicator /?* , these formulas approach the optimized weights, ui k — * oj k . These formulas also yield vanishing weights uik — > for stencils with large values of the smoothness indicator such as those encountered at discontinuities. In this way, the WENO construction retains some of the attributes of the original ENO formulation but with increased accuracy in smooth solution regions and improved differentiability often yielding superior robustness for steady state calculations. 3.2. Higher order accurate FV schemes in multi- dimensions. Although the one-dimensional TVD operators may be readily applied in multi-dimensions on a dimension-by-dimension basis, a result of Goodman and LeVeque, 1985 shows that TVD schemes in two or more space dimensions are only first order accurate. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 28 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Theorem 3.7 (Accuracy of TVD schemes in multi-dimensions) Any two-dimensional finite volume scheme of the form At , „ „ .At u n+l n t „n jn «■> "& ~ |TJ~ ( ^ +1 /2 J ~ Si-1/2,3 ) - jrjTT ( h h+i/2 ~ Kj-1/2) > 1 < * < M, 1 < J < JV wiiA Lipschitz continuous numerical fluxes for integers p, q, r, s 9i+l/2,j ~ 5( u i-p,i-?j---,Ui+r,j+s)> hi,j+l/2 = "( u »— P.j— «> • • • j Ui+r,j+sJ! £/ia£ is toiaZ variation non-increasing in the sense TV(ul +l ) < TV{u n h ) where TV(u) = Y^ [Ay i+ i/ 2 ,j \u i+ ij - Uij\ + Ax i!J+1/ 2 \u itj+1 - u itj \] is at most first-order accurate. Motivated by the negative results of Goodman and LeVeque, weaker conditions yielding solution monotonicity preservation have been developed from discrete maximum principle analysis. These alternative constructions have the positive attribute that they extend to unstructured meshes as well. 3.2.1. Positive coefficient schemes on structured meshes. Theorem 2.6 considers schemes of the form u ; +1 = «? + |^ £ ^*«)« - «?) , VT, e T Ve,- fc €8T,- IT, 1 and provides a local space-time discrete maximum principle min «,<)<u" +1 < max (u?,u?) VTj € T under a CFL-like condition on the time step parameter if all coefficients Cjk are nonnegative. Schemes of this type are often called positive coefficient schemes or more simply positive schemes. To circumvent the negative result of Theorem 3.7, Spekreijse, 1987 developed a family of high order accurate positive coefficient schemes on two-dimensional structured M x N meshes. For purposes of positivity analysis, these schemes are written in incremental form on a M x N logically rectangular 2-D mesh u n+l _ u2j + At ( A^ « +1J . - ulj) + B? J+1 « j+1 - u&) + C?_ ltj «_!,,- - «&) + £>&_! {ul^ - u^) ) , l<i<M,l<j<N (67) where the coefficients are nonlinear functions of the solution ■^i+l,j = -^V— > u i-l,j> u i,ji u i+l,ji •••) Bi,j+1 = B\-—> u i,3-\' U i,ji u i,j+lT--) c?_ hj = C(...,<_ 1J ,< j ,< +1 , i ,...) Dl^ = C(...,< j _ 1 ,< i ,< j+1 ,...) . Once written in incremental form, the following theorem follows from standard positive coefficient maximum principle analysis. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 29 Theorem 3.8 (Positive coefficient schemes in multi-dimensions) The discretization (67) is a positive coefficient scheme if for each 1 < i < M, 1 < j < N and time slab increment [t n ,t n+1 ] A? +1J > 0, Bl j+1 > 0, CVL^ > 0, Dfj^ > 0, (68) and 1 - At (A? +1J + B? J+1 + C^j + Dl^) > (69) with discrete space-time maximum principle and discrete maximum principle at steady state min «-i,jX+i,;Xj-i><;+i) < <j < H«i-ij.«Jnji«y-i.«fj+i) where u* denotes the numerical steady state. Using a procedure similar to that used in the development of MUSCL TVD schemes in 1-D, Spekreijse, 1987 developed a family of monotonicity preserving MUSCL approximations in multi-dimensoins from the positivity conditions of Theorem 3.8. Theorem 3.9 (MUSCL positive coefficient scheme) The fully discrete 2-D finite vol- ume scheme <V = <J - |rj77^i/ 2j - - 9?-i/2j) - J2^r(*&+i/2 - Ki-1/2) , *<i<M,l<j<N utilizing monotone Lipschitz continuous numerical flux functions 9i+l/2,j = ff(«;+l/2,j. U £-l/2j) Ki+1/2 = MV 1/2 ,«tj+l/2) and MUSCL extrapolation formulas U t+l/2,j = u i,i + 2*(V-Rij) Kj - «i+lj) U 7j+l/2 = U iJ + 2*( 5 «J') ( U U ~ U iJ-l) U tj+l/2 = u *J+y iS ( 1 /Si,j)(ll i ,j-Ui lj+1 ) where p _ u i+l,j ~ u ij c _ Ui,j+1 ~ Ui j UiJ 1li—\j UiJ Uij— 1 satisfies the local maximum principle properties of Lemma 3.8 and is second order accurate if the limiter * = $(i?) has the properties that there exist constants /3 € (0, 00), a € [-2, 0] such that Vi? € R a<V(R)<P, -0<^l<2 + a (70) ti Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Ren<5 de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 30 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS with the constraint *(1) = 1 and the smoothness condition $>(R) 6 C 2 near R = 1 together with a time step restriction for stability 1-(1+/?) At dg du n,max + «j <9/i <9u *.j >0 where dg du dh du ',j hi sup ...«. H = 6 I"i"-l/2j'"tl/2,j 1 S€[uT sew: sup ,-l/2'"i,i+l/2 I .,j-l/2'"i,j+l/2 J dh , dh ^S\ „ Many limiter functions satisfy the technical conditions (70) of Theorem 3.9. Some examples include the van Leer limiter • the van Albada limiter W W 1 + \R\ ' W W 1 + R? In addition, Koren, 1988 has constructed the limiter R + 2R 2 * K (i?) = 2-R + 2R 2 which also satisfies the technical conditions (70) and corresponds for smooth solutions in 1-D to the most accurate k = 1/3 MUSCL scheme of van Leer. 3.2.2. FV schemes on unstructured meshes utilizing linear reconstruction. Higher order finite volume extensions of Godunov discretization to unstructured meshes are of the general form (71) (72) -& = - Wl £ 9jk(u jk ,u+ k ) , VT,- € r VejkZdTj with the numerical flux gj k (u,v) given by the quadrature rule Q gjk{uj k ,uj k ) = ^2u q g(vjk(x q );uJ k (x q ),uf k {x q )) , 9=1 where w q £ M. and x q € ejk represent quadrature weights and locations, q = l,...,Q. Given the global space of piecewise constant cell averages, «/, E ^°, the extrapolated states uj k (x) and u\{x) are evaluated using ap-th order polynomial reconstruction operator, R° p : VJ° >->■ V£, u~ k (x) = lim R° (x - e v jk (x) ;u h ) Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 31 Figure 4. Polygonal control volume cell Tj and perimeter quadrature points (solid circles). ut (x) = lim R°(x + e v jk (x) ; u h ) for x € ejk- In addition, it is assumed that the reconstruction satisfies the property YfT! f T . Rp(x~, Uh) dx = Uj stated previously in (64b). In the general finite volume formulation, the control volume shapes need not be convex, see for example Fig. 4. Even so, the solution accuracy and maximum stable time step for explicit schemes may depend strongly on the shape of individual control volumes. In the special case of linear reconstruction, R° (x; u/,), the impact of control volume shape on stability of the scheme can be quantified more precisely. Specifically, the maximum principle analysis presented later for the scheme (71) reveals an explicit dependence on the geometrical shape parameter T geom = gup Q -l^ ( 73 ) O<0<2tt where < a(6) < 1 represents the smallest fractional perpendicular distance from the gravity center to one of two minimally separated parallel hyperplanes with orientation 8 and hyperplane location such that all quadrature points in the control volume lie between or on the hyperplanes as shown in Fig. 5. Table II lists T 9eom values for various control volume shapes in E 1 , R 2 , K 3 , and E d . As might be expected, those geometries that have exact quadrature point symmetry with respect to the control volume gravity center have geometric shape parameters pgeom e q ua i £ 2 regardless of the number of space dimensions involved. Lemma 3.10 (Finite volume interval bounds on unstructured meshes, R.5(x;uh)) The fully discrete finite volume scheme u n + l _ „» _ ^* 1 3i ve^ear, with monotone Lipschitz continuous numerical flux function, nonnegative quadrature weights, and linear reconstructions u~ k (x) = \imR°(x-ei/jk(x);u h ) , x e e jk , u h eV h ° u tk( x ) - limiJ°(x + ei/jjfc(x);u/i) , x e e jk , u h eV h ° , Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 32 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS ^ gravity center • quadrature point Figure 5. Minimally separated hyperplanes h L {9) and h u (0) and the fractional distance ratio a{9) for use in the calculation of r geom . Table II. Reconstruction geometry factors for various control volume shapes utilizing midpoint quadrature rule. control volume shape space dimension pgeom segment 1 2 triangle 2 3 parallelogram 2 2 regular n-gon 2 »/IVl tetrahedron 3 4 parallelepiped 3 2 simplex d d-f-1 hyper-parallelepiped d 2 polytope d Eqn. (73) with extremal trace values at control volume quadrature points £/ min = mm uf k (x q ) , U?** = max uf k (x q ) , x q 6 e jk ' mm - min -^ i<9<<3 exhibits the local interpolated interval bound jmin,n _,_ f-| _ „.\„,n <? „, n +l l<q<Q ajU™'" + (1 - aj)u] < u] +l < (1 - aj) uj + ajUf with the time step proportional interpolation parameter o~j defined by At 3 ~ \Tj pgeom V^ sup v.,-»eer,- se[t/?" •,",•• 1 dg du+ (Vjk(Xq));U,u) (75) (76) l<q<Q = ,_,...min,n ..max.jii that depends on the shape parameter r geom defined in (73). Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 33 Given the two-sided bound of Lemma 3.10, a discrete maximum principle is obtained under a CFL-like time step restriction if the limits Uf 13 * and £f min can be bounded from above and below respectively by the neighboring cell averages. This idea is given more precisely in the following theorem. Theorem 3.11 (Finite volume maximum principle on unstructured meshes, R°) Let u™ ,n and u max denote the minimum and maximum value of solution cell averages for a given cellTj and corresponding adjacent cell neighbors, i.e. uf n = min ( Uj ,u k ) and w raax = max ( Uj ,u k ) . (77) The fully discrete finite volume scheme with monotone Lipschitz continuous numerical flux function, nonnegative quadrature weights, and linear reconstructions u Jk( x ) - limR%(x-ei/j k (x);uh) , x e e jk , u h €V° u tk( x ) - ^^(x + evj k (xy,u h ) ' xGe Jk , u h eV° (79) exhibits the local space-time maximum principle for each Tj 6 T min «,u?)<< +1 < max (u?,u?) as well as the local spatial maximum principle at steady state (u n+1 = u n = u") min ui < u*; < max ut Ve jk edTj " ~ J -v ejk edT i if the linear reconstruction satisfies Vej k € dTj and x q G ej k ,q = 1, . . . , Q ■^.n,-!",,. min.n min,n\ -. — ,n/_ \ ^- _,;„/., max, n max,n\ (qiW max.[u j , u k ) < u jk ' (x q ) < min(u J . ' ,u k ) (80) under the time step restriction x _ A* r g eom £ sup \^S(u jk {x q ));u,H>0 l-i ■» I , -., min.n max, n, I C7U \ !<?<« 5 € [u , 7 ,in '",„ , ?"'' x '"l 1 1 ' 1 ' with r<5 eom defined in (73). Note that a variant of this theorem also holds if the definition of u max and u min are expanded to include more control volume neighbors. Two alternative definitions frequently used when the control volume shape is a simplex are given by uf in = min u k and u? ax = max u k . (81) J T k €T 3 T k €T Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Ren6 de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 34 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS These expanded definitions include adjacent cells whose intersection with Tj in E d need only be a set of measure zero or greater. Slope limiters for linear reconstruction. Given a linear reconstruction R\{x;Uh) that does not necessarily satisfy the requirements of Theorem 3.11, it is straightforward to modify the reconstruction so that the new modified reconstruction does satisfy the requirements of Theorem 3.11. For each control volume Tj € T a modified reconstruction operator R\(x;v,h) of the form _ R\{x\ u h )\ Tj = uj + a Tj (R°(x; u h )\ Tj - uj) is assumed for ax, £ [0, 1]. By construction, this modified reconstruction correctly reproduces the control volume cell average for all values of a^. , i.e. mL** Uh) dx = Uj . (82) The most restrictive value of a^ for each control volume Tj is then computed based on the Theorem 3.11 constraint (80), i.e. ' min(u ■ ,u k )—Vj V MM _ a T = mm \<q<Q S:uV;"' if *?(*.; <* > min K max > ^ ax ) "SiruiC-^ tf^(*.;« h )h<max(«f-,t.r) (83) 1 otherwise where u max and u min are defined in (77). When the resulting modified reconstruction operator is used in the extrapolation formulas (79), the discrete maximum principle of Theorem 3.11 is attained under a CFL-like time step restriction. By utilizing the inequalities max(uj, u fc ) < min(«J iax > u^ ax ) and mm(uj,u k ) > max(uf n ,< in ) it is straightforward to construct a simpler but more restrictive limiter function ' max(uj,Ufc)- Uj l<q<Q %SZ$\£Z, if^(x g ;u,)| Ti >maxK-,u fc ) ^Kofe iim(^nn)\ Ti <min(uj,u k ) (84) 1 otherwise that yields modified reconstructions satisfying the technical conditions of Theorem 3.11. This simplified limiter (84) introduces additional slope reduction when compared to (83). This can be detrimental to the overall accuracy of the discretization. The limiter strategy (84) and other variants for simplicial control volumes are discussed further in Liu, 1993; Wierse, 1994; Batten, Lambert and Causon, 1996. In Barth and Jespersen, 1989, a variant of (83) was proposed R 7 ' a? = mm BUZu^-ui i*fl?(*,;«h)|T,>«r Uj -Uj if i?? (a:,; Ufc) | Tj <uf n ■ ( 85 ) 1 otherwise Although this limiter function does not produce modified reconstructions satisfying the requirements of Theorem 3.11, using Lemma 3.10 it can be shown that the Barth and Jespersen Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 35 limiter yields finite volume schemes (74) possessing a global extremum diminishing property, i.e. that the solution maximum is non-increasing and the solution minimum is nondecreasing between successive time levels. This limiter function produces the least amount of slope reduction when compared to the limiter functions (83) and (84). Note that in practical implementation, all three limiters (83), (84) and (85) require some modification to prevent near zero division for nearly constant solution data. 3.2.3. Linear reconstruction operators on simplicial control volumes. Linear reconstruction operators on simplicial control volumes that satisfy the cell averaging requirement (64b) often exploit the fact that the cell average is also a pointwise value of any valid linear reconstruction evaluated at the gravity center of a simplex. This reduces the reconstruction problem to that of gradient estimation given pointwise samples at the gravity centers. In this case, it is convenient to express the reconstruction in the form #°(x; u h )\ Tj = Uj + (Vu h )T, • (x - x*) (86) where a;® denotes the gravity center for the simplex Tj and (Vu/,)t, is the gradient to be determined. Figure 6 depicts a 2-D simplex A123 and three adjacent neighboring simplices. Also shown are the corresponding four pointwise solution values {A, B, C, 0} located at gravity centers of each simplex. By selecting any three of the four pointwise solution values, a set of four possible gradients are uniquely determined, i.e. {V(ABC), V(ABO), V(BCO), V(CAO)}. Using the example of Fig. 6, a number of slope limited reconstruction techniques are possible Figure 6. Triangle control volume A123 (shaded) with three adjacent cell neighbors. for use in the finite volume scheme (78) that meet the technical conditions of Theorem 3.11. 1. Choose (Vu/ 1 )t 123 = V(ABC) and limit the resulting reconstruction using (83) or (84). This technique is pursued in Barth and Jespersen, 1989 but using the limiter (85) instead. 2. Limit the reconstructions corresponding to gradients V(ABC), V(ABO), V(BCO) and V(CAO) using (83) or (84) and choose the limited reconstruction with largest gradient magnitude. This technique is a generalization of that described in Batten, Lambert and Causon, 1996 wherein limiter (84) is used. 3. Choose the unlimited reconstruction V(ABC), V{ABO), V{BCO) and V(CAO) with largest gradient magnitude that satisfies the maximum principle reconstruction bound inequality (80). If all reconstructions fail the bound inequality, the reconstruction gradient is set equal to zero, see Liu, 1993. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 36 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Figure 7. Triangulation of gravity center locations showing a typical control volume To associated with the triangulation vertex vo with cyclically indexed graph neighbors Tk , k = 1, . . . , iVo . 3.2.4- Linear reconstruction operators on general control volumes shapes. In the case of linear reconstruction on general volume shapes, significant simplification is possible when compared to the general p-exact reconstruction formulation given in Sect. 3.2.5. It is again convenient to express the reconstruction in the form i??(x;u A )| T , = Uj + (Vu/Ot,- • (x - x?) (87) where x® denotes the gravity center for the control volume Tj and (Vtth)^ is the gradient to be determined. Two common techniques for simplified linear reconstruction include a simplified least squares technique and a Green-Gauss integration technique. Simplified least squares linear reconstruction. As was exploited in the linear reconstruction techniques for simplicial control volumes, linear reconstructions satisfying (64b) on general control volume shapes are greatly simplified by exploiting the fact that the cell average value is also a pointwise value of all valid linear reconstructions evaluated at the gravity center of a general control volume shape. This again reduces the linear reconstruction problem to that of gradient estimation given pointwise values. In the simplified least squares reconstruction technique, a triangulation (2D) or tetrahedralization (3D) of gravity centers is first constructed as shown in Fig. 7. Referring to this figure, for each edge of the simplex mesh incident to the vertex v , an edge projected gradient constraint equation is constructed subject to a pre- specified nonzero scaling Wk w k (Vu/,)t • (xl - Xo) = w k(uk ~ u ) ■ The number of edges incident to a simplex mesh vertex in R d is greater than or equal to d thereby producing the following generally non-square matrix of constraint equations w No Ax s Nq toi Ay? WNo^VNa. (Vu fc ) To Wi(u\ -Uq) \ ^jv («JV -«o)/ or in abstract form [U U] vu = / . Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 37 This abstract form can be symbolically solved in a least squares sense using an orthogonalization technique yielding the closed form solution V« = 1 ( l22{Zl ' ^ " /l2 ^ 2 ' ?)\ (&&) I11I22 - l\ 2 \l n (L 2 ■ f) - l 12 (Li ■ f)J K ' with kj = Li ■ Lj . The form of this solution in terms of scalar dot products over incident edges suggests that the least squares linear reconstruction can be efficiently computed via an edge data structure without the need for storing a non-square matrix. Green-Gauss linear reconstruction. This reconstruction technique specific to simplicial meshes assumes nodal solution values at vertices of the mesh which uniquely describes a C° linear interpolant, w/,. Gradients are then computed from the mean value approximation |n | (Vu/On,, * / Vu h dx= / u h dv . (89) JQo JdU For linear interpolants, the right-hand side term can be written in the following equivalent k+l k-i Figure 8. Median dual control volume To demarcated by median segments of triangles incident to the vertex vo with cyclically indexed adjacent vertices Vk, k = 1, ... No. form using the configuration depicted in Fig. 8 Vm dx = ^T l -(u + u k ) vak f JClc No 2 V where u k represents any path integrated normal connecting pairwise adjacent simplex gravity centers, i.e. u ok = du . (90) A particularly convenient path is one that traces out portions of median segments as shown in Fig. 8. These segments demarcate the so-called median dual control volume. By construction, the median dual volume |T | is precisely equal to |n |/3 in 2-D. Consequently, a linear reconstruction operator on non-overlapping median dual control volumes is given by No 1 |To|(Vu h )r « Yl 2( u o + Uk)vok ■ (91) Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Ren6 de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 38 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS The gradient calculation is exact whenever the numerical solution varies linearly over the support of the reconstruction. Since mesh vertices are not located at the gravity centers of median dual control volumes, the cell averaging property (64b) and the bounds of Theorem 3.11 are only approximately satisfied using the Green-Gauss technique. A number of slope limited linear reconstruction strategies for general control volume shapes are possible for use in the finite volume scheme (78) that satisfy the technical conditions of Theorem 3.11. Using the example depicted in Fig. 7, let Vk+i/2 u h denote the unique linear gradient calculated from the cell average set {uo,Uk,Uk+i}- Three slope limiting strategies that are direct counterparts of the simplex control volume case are: 1. Compute (Vu/i)t using the least squares linear reconstruction or any other valid linear reconstruction technique and limit the resulting reconstruction using (83) or (84). 2. Limit the reconstructions corresponding to the gradients Vjfc+i/ 2 u/,,fc = l,...,No and (Vu/i)r using (83) or (84) and choose the limited reconstruction with largest gradient magnitude. 3. Choose the unlimited reconstruction from V* + i/2U/»,fc = l,...,No and (Vuh)r with largest gradient magnitude that satisfies the maximum principle reconstruction bound inequality (80). If all reconstructions fail the bound inequality, the reconstruction gradient is set equal to zero. 3.2.5. General p-exact reconstruction operators on unstructured meshes. Abstractly, the reconstruction operator serves as a finite-dimensional pseudo inverse of the cell averaging operator A whose j'-th component Aj computes the cell average of the solution in Tj iu= mLr The development of a general polynomial reconstruction operator, R®, that reconstructs p- degree polynomials from cell averages on unstructured meshes follows from the application of a small number of simple properties. 1. (Conservation of the mean) Given solution cell averages Uh, the reconstruction RpU^ is required to have the correct cell average, i.e. More concisely, if v = R°Uh then Uh = Av . AR° p = I so that R° is a right inverse of the averaging operator A. 2. (p-exactness) A reconstruction operator i?° is p-exact if RpA reconstructs polynomials of degree p or less exactly. Denoting by V v the space of all polynomials of degree p, if u € V p and v — Au then RpV = u . This can be written succinctly as R° p A\ Vp = I so that R° is a left inverse of the averaging operator A restricted to the space of polynomials of degree at most p. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 39 3. (Compact support) The reconstruction in a control volume Tj should only depend of cell averages in a relatively small neighborhood surrounding Tj. Recall that a polynomial of degree p in R d contains ( p ^ d ) degrees of freedom. The support set for Tj is required to contain at least this number of neighbors. As the support set becomes even larger for fixed p, not only does the computational cost increase, but eventually the accuracy decreases as less valid data from further away is brought into the calculation. Practical implementations of polynomial reconstruction operators fall into two classes: • Fixed support stencil reconstructions. These methods choose a fixed support set as a preprocessing step. Various limiting strategies are then employed to obtain non- oscillatory approximation, see for example Barth and Frederickson, 1990; Delanaye, 1996 for further details. • Adaptive support stencil reconstructions. These ENO-like methods dynamically choose reconstruction stencils based on solution smoothness criteria, see for example Harten and Chakravarthy, 1991; Vankeirsblick, 1993; Abgrall, 1994; Sonar, 1997; Sonar, 1998 for further details. 3.2.6. Positive coefficient schem.es on unstructured meshes Several related positive coefficient schemes have been proposed on multi-dimensional simplicial meshes based on one-dimensional interpolation. The simplest example is the upwind triangle scheme as introduced by Billey et al., 1987; Desideri and Dervieux, 1988; Rostand and StoufHet, 1988 with later improved variants given by Jameson, 1993; Cournede, Debiez and Dervieux, 1998. These schemes are not Godunov methods in the sense that a single multi-dimensional gradient is not obtained in each control volume. The basis for these methods originates from the gradient estimation formula (91) generalized to the calculation of flux divergence on a median dual tessellation. In deriving this flux divergence formula, the assumption has been made that flux components vary linearly within a simplex yielding the discretization formula [ div(f)dx= [ f-dv = Y, \U(uj)+f{u k ))-u jk JT i JdT = Ve jk €dTj where Vjk is computed from a median dual tessellation using (90). This discretization is the unstructured mesh counterpart of central differencing on a structured mesh. Schemes using this discretization of flux divergence lack sufficient stability properties for computing solutions of general nonlinear conservation laws. This lack of stability can be overcome by adding suitable diffusion terms. One of the simplest modifications is motivated by upwind domain of dependence arguments yielding the numerical flux gjk(uj,u k ) = - (f(uj) + f(u k )) - Vj k - -\a\ jk A jk u (92) with a,jk a mean value (a.k.a. Murman-Cole) linearization satisfying Vjk ■ Ajkf = a,jk Aj k u . Away from sonic points where f'(u*) = for u* € [v,j,Uj+i], this numerical flux is formally an E-flux satisfying (28). With suitable modifications of a,j k near sonic points, it is then possible to produce a modified numerical flux that is an E-flux for all data, see Osher, 1984. Theorems Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 40 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 2.6, 2.7 and 2.8 show that schemes such as (22) using E-fluxes exhibit local discrete maximum principles and L^ stability. Unfortunately, schemes based on (92) are too dissipative for most practical calculations. The main idea in the upwind triangle scheme is to add anti-diffusion terms to the numerical flux function (92) such that the sum total of added diffusion and anti-diffusion terms in the numerical flux function vanish entirely whenever the numerical solution varies linearly over the support of the flux function. In all remaining situations, the precise amount of anti-diffusion is determined from maximum principle analysis. Figure 9. Triangle complex used in the upwind triangle schemes showing the lineax extension of e,jt into neighboring triangle for the determination of points x^ and xy . Theorem 3.12 (Maximum Principles for the Upwind Triangle Scheme) Let T de- note the median dual tessellation of an underlying simplicial mesh. Also let Uj denote the nodal solution value at a simplex vertex in one-to-one dual correspondence with the control volume Tj € T such that a C° linear solution interpolant is uniquely specified on the simplicial mesh. Let gjk{uji,Uj,Uk,Uk') denote the numerical flux function with limiter function $(•) : M >-> M. gjkiuj^u^UkiUk^^-ifiu^+fiuk))-^ - -af k (l - ^ / ^ A^'fcW ) AjkU utilizing the mean value speed ajk satisfying Vjk ■ Ajfc/ = ajk Ajfcit and variable spacing parameter hjk = |Ajao;|. The fully discrete finite volume scheme u] +1 = «? - |^ £ 9ik {u} , u?, u n k , ul ) , VT,- € T Veneer, \Tj, with linearly interpolated values Ujt and u^ as depicted in Fig. 9 exhibits the local space-time maximum principle min «X)<u" +1 < max «,<) and the local spatial maximum principle at steady state (u n+1 — u n = u*) min ut < u", < max ut Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 41 if the limiter *(i2) satisfies VR€R < <£>(R)/R , < tf(iZ) < 2 . Some standard limiter functions that satisfy the requirements of Theorem 3.12 include • the MinMod limiter with maximum compression parameter equal to 2 * MM (i?) = max(0,min(fl,2)) • the van Leer limiter 1 ' 1 + \R\ ' Other limiter formulations involving three successive one-dimensional slopes are given in Jameson, 1993; Cournede, Debiez and Dervieux, 1998. 4. Further Advanced Topics The material presented in previous sections gives a brief overview of the derivation and analysis of finite volume methods. For simplicity and brevity of the presentation, exclusive attention has been devoted to scalar nonlinear conservation laws in divergence form. In this overview, special consideration has been given to the formulation and stability analysis of higher order accurate schemes since these developments have had the largest impact on development of industrial computer codes in use at the time of this writing. The remainder of this chapter will consider several extensions of the finite volume method. Section 4.1 considers a class of higher order accurate discretizations in time that still preserve the stability properties of the fully-discrete schemes using Euler time integration. This is followed by a discussion of generalizations of the finite volume method for problems including second order diffusion terms and the extension to systems of nonlinear conservation laws. 4-1. Higher order time integration schemes The derivation of finite volume schemes in Sect. 2 first yielded a semi-discrete formulation (21) that was later extended to a fully-discrete formulation (22) by the introduction of first order accurate forward Euler time integration. These latter schemes where then subsequently extended to higher order accuracy in space using a variety of techniques. For many computing problems of interest, first order accuracy in time is then no longer enough. To overcome this low order accuracy in time, a general class of higher order accurate time integration methods was developed that preserve stability properties of the fully-discrete scheme with forward Euler time integration. Following Gottlieb, Shu and Tadmor, 2001 and Shu, 2001, these methods will be referred to as Strong Stability Preserving (SSP) time integration methods. Explicit SSP Runge-Kutta methods were originally developed by Shu, 1988; Shu and Osher, 1988 and Gottlieb and Shu, 1998 and called TVD Runge-Kutta time discretizations. In a slightly more general approach, total variation bounded (TVB) Runge-Kutta methods were considered by Cockburn and Shu, 1989; Cockburn, Lin and Shu, 1989; Cockburn, Hou and Shu, 1990; Cockburn and Shu, 1998 in combination with the discontinuous Galerkin discretization Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 42 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS in space. Kiither, 2000 later gave error estimates for second order TVD Runge-Kutta finite volume approximations of hyperbolic conservation laws. To present the general framework of SSP Runge-Kutta methods, consider writing the semi- discrete finite volume method in the following form jU{t) = L(U(t)) (93) where U = U(t) denotes the solution vector of the semi-discrete finite volume method. Using this notation together with forward Euler time integration yields the fully-discrete form Jjn+l =U n_ At L (Jjn} ^ (94) where U n is now an approximation of U(t n ). As demonstrated in Section 2.2, the forward Euler time discretization is stable with respect to the Z,°°-norm, i.e. ||^ n+1 |loo < Halloo > (95) subject to a CFL-like time step restriction At <At . (96) With this assumption, a time integration method is said to be SSP (see Gottlieb, Shu and Tadmor, 2001) if it preserves the stability property (95), albeit with perhaps a slightly different restriction on the time step At < c At (97) where c is called the CFL coefficient of the SSP method. In this framework, a general objective is to find SSP methods that are higher order accurate, have low computational cost and storage requirements, and have preferably a large CFL coefficient. Note that the TVB Runge-Kutta methods can be embedded into this class if the following relaxed notion of stability is assumed ll^ n+1 ||co<(l + 0(At))||C/ n |U . (98) 4-1.1. Explicit SSP Runge-Kutta methods. Following Shu and Osher, 1988 and the review articles by Gottlieb, Shu and Tadmor, 2001; Shu, 2001, a general m stage Runge-Kutta method for integrating (93) in time can be algorithmically represented as U° := U n , i-i U l := J2( a "*U k +Pik&tL(U k )), a lk >0, l = l,...,m, (99) fc=0 rrn+l ._ fjm To ensure consistency, the additional constraint X^=o a <-k = 1 is imposed. If, in addition, all Pin are assumed to be non-negative, it is straightforward to see that the method can be written as a convex (positive weighted) combination of simple forward Euler steps with At replaced by ^-At. From this property, Shu and Osher, 1988 concluded the following lemma: Lemma 4.1. If the forward Euler method (94) is L°° -stable subject to the CFL condition (96), then the Rung-Kutta method (99) with pik > is SSP, i.e. the method is L°° -stable under the time step restriction (97) with CFL coefficient c = min— . (100) i,k aik In the case of negative /?u, a similar result can be proven, see Shu and Osher, 1988. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 43 4-1.2. Optimal second and third order nonlinear SSP Runge-Kutta methods. Gottlieb, Shu and Tadmor, 2001 [Proposition 3.1] show that the maximal CFL coefficient for any m-stage, ru- th order accurate SSP Runge-Kutta methods is c = 1. Therefore, SSP Runge-Kutta methods that achieve c = 1 are termed "optimal" . Note that this restriction is not true if the number of stages is higher than the order of accuracy, see Shu, 1988. Optimal second and third order nonlinear SSP Runge-Kutta methods are given in Shu and Osher, 1988. The optimal second order, two-stage non-linear SSP Runge-Kutta method is given by 17° := U n , U 1 := U° + AtL(U°) , £ £ JL This method corresponds to the well known method of Heun. Similarly, the optimal third order, three-stage non-linear SSP Runge-Kutta method is given by U° := U n , U 1 := U° + AtL(U°) ,. U 2 := -V° + -V 1 + ]&t L{U l ) , 4 4 4 Tjn+1 . = 1^0 + 2^2 + 2 Af L{ jj 2) _ 000 Further methods addressing even higher order accuracy or lower storage requirements are given in the review articles of Gottlieb, Shu and Tadmor, 2001 and Shu, 2001 where SSP multi-step methods are also discussed. 4-2. Discretization of elliptic problems Finite volume methods for elliptic boundary value problems have been proposed and analyzed under a variety of different names: box methods, covolume methods, diamond cell methods, integral finite difference methods and finite volume element methods, see Bank and Rose, 1987; Cai, 1991; Siili, 1991; Lazarov, Michev and Vassilevsky, 1996; Viozat et al, 1998; Chatzipantelidis, 1999; Chou and Li, 2000; Hermeline, 2000; Eymard, Galluoet and Herbin, 2000; Ewing, Lin and Lin, 2002. These methods address the discretization of the following standard elliptic problem in a convex polygonal domain ft C M 2 -V-AVu = f in ft , ' (101) u(x) = on dft for A e E 2x2 , a symmetric positive definite matrix (assumed constant). Provided / e H®{Q) then a solution u exists such that u £ Hq +2 (Q), -1 < < 1,/? # ±1/2, where H s (£l) denotes the Sobolev space of order s in ft. Nearly all the above mentioned methods can be recast in Petrov-Galerkin form using a piecewise constant test space together with a conforming trial space. A notable exception is given in Chatzipantelidis, 1999 wherein nonconforming Crouzeix-Raviart elements are utilized and analyzed. To formulate and analyze the Petrov-Galerkin representation, two tessellations Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 44 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS of fi are considered: a triangulation T with simplicial elements K € T and a dual tessellation T* with control volumes T € T*. In the class of conforming trial space methods such as the finite volume element (FVE) method, a globally continuous, piecewise p-th. order polynomial trial space with zero trace value on the physical domain boundary is constructed X h = {ve C°(Q) | v\ K € V V {K), VK€T and v\ dn = 0} using nodal Lagrange elements on the simplicial mesh. A dual tessellation T* of the Lagrange element is then constructed, see for example Fig. 10 which shows a linear Lagrange element with two dual tessellation possibilities. These dual tessellated regions form control volumes for the finite volume method. The tessellation technique extends to higher order Lagrange a. Voronoi dual tessellation b. Median dual tessellation Figure 10. Two control volume valiants used in the finite volume discretization of second order derivative terms: (a) Voronoi dual where edges of the Voronoi dual are perpendicular to edges of the triangulation and (b) median dual formed from median dual segments in each triangle. elements in a straightforward way. A piecewise constant test space is then constructed using r* Yh = {v\v\ T e X (T), VTZT*} where x(T) is a characteristic function in the control volume T. The finite volume element discretization of (101) then yields the following Petrov-Galerkin formulation: Find Uh € Xh such that Y\ I I Wh AVu h -dv+ [ w h fdx) =0 , Vw A £ Y h . (102) The analysis of (102) by Ewing, Lin and Lin, 2002 using linear elements gives an o priori estimate in an L 2 norm that requires the least amount of solution regularity when compared to previous methods of analysis. Theorem 4.2 (FVE a priori error estimate, Ewing, Lin and Lin, 2002) Assume a 2- D quasi-uniform triangulation T with dual tessellation T* such that 3 C > satisfying C^h 2 < \T\<Ch 2 , VTGT* . Assume that u and uh are solutions of (101) and (102) respectively with u £ H 2 (Cl), f € H& , (0 < /? < 1) . Then 3C > such that the a priori estimate holds h-u h \\ L 2 (Q) <C'{h 2 \\u\\ HHn) +h 1+ P\\f\\ HHn) ) . (103) Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 45 Unlike the finite element method, the error estimate (103) reveals that optimal order convergence is obtained only if / £ H® with /? > 1. Moreover, numerical results show that the source term regularity can not be reduced without deteriorating the measured convergence rate. Optimal convergence rates are also shown for the nonconforming Crouzeix-Raviart element based finite volume method analyzed by Chatzipantelidis, 1999 for u € iT 2 (fl) and / G ff 1 (fi). An extensive presentation and analysis of finite volume methods for elliptic equations without utilizing a Petrov-Galerkin formulation is given in Eymard, Galluoet and Herbin, 2000. In this work, general boundary conditions that include non-homogeneous Dirichlet, Neumann and Robin conditions are discussed. In addition, the analysis is extended to general elliptic problems in divergence form including convection, reaction and singular source terms. 4-3. Conservation laws including diffusion terms As demonstrated in Sect. 1, hyperbolic conservation laws are often approximations to physical problems with small or nearly vanishing viscosity. In other problems, the quantitative solution effects of these small viscosity terms are actually sought. Consequently, it is necessary in these problems to include viscosity terms into the conservation law formulation. As a model for these latter problems, a second order Laplacian term with small diffusion parameter is added to the first order Cauchy problem, i.e. d t u + V • f(u) - eAu = in E d x R+ , (104a) u(x, 0) = u in R d . (104b) Here u(x, i) : R d x R+ ->■ K denotes the dependent solution variable, / € C(E) the hyperbolic flux and e > a small diffusion coefficient. Application of the divergence and Gauss theorems to (104a) integrated in a region T yields the following integral conservation law form — udx + f(u) ■ dv — sVu ■ dv = . Ot Jt JdT J8T (105) A first goal is to extend the fully-discrete form (22) of Sect. 2 to the integral conservation law (105) by the introduction of a numerical diffusion flux function dj k (uh) for a control volume Tj e T such that / eVu-di/tt 22 djk(uh) ■ When combined with the general finite volume formulation (22) for hyperbolic conservation laws, the following fully-discrete scheme is produced Uj " +1 = U ?-— £ (g jk (u?,u?)-d jk (uZ)) , VTjeT. (106) 1 jl ejk€dTj In this equation, the index m may be chosen either as n or n + 1, corresponding to an explicit or implicit discretization. 4-3.1. Choices of the numerical diffusion flux dj k . The particular choice of the numerical diffusion flux function dj k depends on the type of control volume that is used. Since the approximate solution Uh is assumed to be a piecewise constant function, the definition of dj k Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 46 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS involves a gradient reconstruction of Uh in the normal direction to each cell interface ejk- The reconstruction using piecewise constant gradients is relatively straightforward if the control volumes are vertex-centered, or if the cell interfaces are perpendicular to the straight lines connecting the storage locations (see Fig. 10). Vertex- centered finite volume schemes. In the case of vertex-centered control volumes such as the median dual control volume, a globally continuous, piecewise linear approximate solution Uh is first reconstructed on the primal mesh. Vu/, is then continuous on the control volume interfaces and the numerical diffusion flux straightforwardly computed as d jk (u^)= I VfiJT ■ */,-* . (107) Cell-centered finite volume schemes. In the case of cell-centered finite volume schemes where an underlying primal-dual mesh relationship may not exist, a simple numerical diffusion flux can be constructed whenever cell interfaces are exactly or approximately perpendicular to the straight lines connecting the storage locations, e.g. Voronoi meshes, quadrilateral meshes, etc. In these cases, the reconstructed gradient of u/, projected normal to the cell interface ejk can be represented by \Xk-Xj\ where n denotes the storage location of cell T< . The numerical diffusion flux for this case is then given by djk (u^) = -^-(u k --uf) . (108) \ x k •Ej\ Further possible constructions and generalizations are given in Eymard, Gallouet and Herbin, 2001; Gallouet, Herbin and Vignal, 2000; Herbin and Ohlberger, 2002. 4-3.2. Note on stability, convergence and error estimates. Stability analysis reveals a CFL- like stability condition for the explicit scheme choice (m = n) in (106) a 3 (h n ■ ) 2 \f n < \"'min) ~ aL 9 h min+ £ where L g denotes the Lipschitz constant of the hyperbolic numerical flux, a is a positive mesh dependent parameter and e is the diffusion coefficient. In constructing this bound, a certain form of shape regularity is assumed such that there exists an a > such that for all j, k with hk = diam (T k ) Qh\ < \T k \, a \dT k \ < h k , ah k < \x k - x t \ . (109) Thus, Ai n is of the order h 2 for large s and of the order h for e < h. In cases where the diffusion coefficient is larger than the mesh size, it is therefore advisable to use an implicit scheme (m = n + 1). In this latter situation, no time step restriction has to be imposed (see Eymard et al., 2002; Ohlberger, 2001b). In order to demonstrate the main difficulties when analyzing convection dominated problems, consider the following result from Feistauer et al., 1999 for a homogeneous diffusive boundary value problem. In this work, a mixed finite volume finite element method sharing similarities Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 47 with the methods described above is used to obtain a numerical approximation u/, of the exact solution u. Using typical energy-based techniques, they prove the following a priori error bound. Theorem 4.3. For initial data u G L°°(R 2 ) D W 1 ' 2 ^ 2 ) and r > there exist constants ci , C2 > independent of e such that \\u{;t n ) - U h {;t n )\\ L 2 (a) < Cl he C * T / e . (110) This estimate is fundamentally different from estimates for the purely hyperbolic problems of Sects. 2 and 3. Specifically, this result shows how the estimate strongly depends on the small parameter e; ultimately becoming unbounded as e tends to zero. In the context of convection dominated or degenerate parabolic equations, Kruzkov- techniques have been recently used by Carrillo, 1999; Karlsen and Risebro, 2000 in proving uniqueness and stability of solutions. Utilizing these techniques, convergence of finite volume schemes (uniform with respect to e -> 0) was proven in Eymard et al., 2002 and a priori error estimates were obtained for viscuous approximations in Jakobsen and Karlsen, 2001 and Eymard, Gallouet and Herbin, 2002. Finally, in Ohlberger, 2001a; Ohlberger, 2001b uniform a posteriori error estimates suitable for adaptive meshing are given. 4-4- Extension to systems of nonlinear conservation laws A positive attribute of finite volume methods is the relative ease in which the numerical discretization schemes of Sects. 2 and 3 can be algorithmically extended to systems of nonlinear conservation laws of the form d t u + V • f{u) = in R d x 1+ , (Ilia) u(x,0) = u (x) inR d (Hlb) where u(x,t) :l <l xi + — > R m denotes the vector of dependent solution variables, f(u) : R m i->- E mxd denotes the flux vector, and u (x) : R d -»■ R m denotes the initial data vector at time t = 0. It is assumed that this system is strictly hyperbolic, i.e. the eigenvalues of the flux jacobian A{v; u) H df/du ■ v are real and distinct for all bounded v G JR d . The main task in extending finite volume methods to systems of nonlinear conservation laws is the construction of a suitable numerical flux function. To gain insight into this task, consider the one-dimensional linear Cauchy problem for u(x, t) : R x R + i-> R m and uq(x) : R i-t R 771 d t u + d x (Au) = inlxR + , u(x,0) = u (x) in R (112) where A £ R mxm is a constant matrix. Assume the matrix A has m real and distinct eigenvalues, Ai < Ai < • • • < A m , with corresponding right and left eigenvectors denoted by r k € R m and l k G R m respectively for k = 1, . . . ,m. Furthermore, let X G R mxm denote the matrix of right eigenvectors, X = [ri,...,r m ], and A G R mxm the diagonal matrix of eigenvalues, A = diag(Ai, . . . , A m ) so that A = XKX~ X . The one-dimensional system (112) is readily decoupled into scalar equations via the transformation into characteristic variables a = X^u d t a + d x {Aa) =0 in R x R + , Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 48 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS a(x,0) = a {x) in E (113) and component-wise solved exactly a {k) (x,t)=a { h) (x-X k t) , k = 1, . . . , m or recombined in terms of the original variables m fc=l Using this solution, it is straightforward to solve exactly the associated Riemann problem for w(£,t) eR m d T w + dz(Aw) = in R x R + with initial data lt ., (u if £ < w «.°) = {t, if£>0 thereby producing the following Godunov-like numerical flux function g(u,v) = Aw(0,R + ) = l(Au + Av)-±\A\(v-v) (114) with | A | = X|A|JY" _1 . When used in one-dimensional discretization together with piecewise constant solution representation, the linear numerical flux (114) produces the well-known Courant-Isaacson-Rees (CIR) upwind scheme for linear systems of hyperbolic equations where A ± = XA ± X _1 . Note that higher order accurate finite volume methods with slope limiting procedures formally extend to this linear system via component wise slope limiting of the characteristic components a^*^ , k = 1, . . . m for use in the numerical flux (114). 4-4-1- Numerical flux functions for systems of conservation laws. In Godunov's original work (see Godunov, 1959), exact solutions of the one-dimensional nonlinear Riemann problem of gas dynamics were used in the construction of a similar numerical flux function g G (u,v) = f(w(0,R + ))-u (115) where w(£, r) £ R m is now a solution of a nonlinear Riemann problem d T w + d(:f {v \w) =0 in K x E + with initial data u fu if £ < ^'°) = {, if£>0 ■ Recall that solutions of the Riemann problem for gas dynamic systems are a composition of shock, contact and rarefaction wave family solutions. For the gas dynamic equations considered by Godunov, a unique solution of the Riemann problem exists for general states u and v Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 49 except those states producing a vacuum. Even so, the solution of the Riemann problem is both mathematically and computationally nontrivial. Consequently, a number of alternative numerical fluxes have been proposed that are more computationally efficient. These alternative numerical fluxes can be sometimes interpreted as approximate Riemann solvers. A partial list of alternative numerical fluxes is given here. A more detailed treatment of this subject is given in Godlewski and Raviart, 1991, Kroner, 1997, and LeVeque, 2002. • Osher-Solomon flux (Osher and Solomon, 1982). This numerical flux is a system generalization of the Enquist-Osher flux of Sect. 2. All wave families are approximated in state space as rarefaction or inverted rarefaction waves with Lipschitz continuous partial derivatives. The Osher-Solomon numerical flux is of the form 9° S (u,v) = \{f{u) + /(«)) v-\[ \A{u;w)\dw where \A\ denotes the usual matrix absolute value. By integrating on m rarefaction wave integral subpaths that are each parallel to a right eigenvector, a system decoupling occurs on each subpath integration. Furthermore, for the gas dynamic equations with ideal gas law, it is straightforward to construct m — 1 Riemann invariants on each subpath thereby eliminating the need for path integration altogether. This reduces the numerical flux calculation to purely algebraic computations with special care taken at sonic points, see Osher and Solomon, 1982. • Roe flux (Roe, 1981). Roe's numerical flux can be a interpreted as approximating all waves families as discontinuities. The numerical flux is of the form 9 R ° e (u,v) = !(/(«) + /(«)) • v - \\A{v;u,v)\ (v - u) where A(u; u, v) is the "Roe matrix" satisfying the matrix mean value identity (f(v) - /(«)) • v = A{v; u, v) (v - u) with A(v; u, u) = A(v; u). For the equations of gas dynamics with ideal gas law, the Roe matrix takes a particularly simple form. Steady discrete mesh-aligned shock profiles are resolved with one intermediate point. The Roe flux does not preclude the formation of entropy violating expansion shocks unless additional steps are taken near sonic points. • Steger- Warming flux vector splitting (Steger and Warming, 1981). Steger and Warming considered a splitting of the flux vector for the gas dynamic equations with ideal gas law that exploited the fact that the flux vector is homogeneous of degree one in the conserved variables. From this homogeneity property, Euler's identity then yields that f(u) - v = A(y\ it) u. Steger and Warming then considered the matrix splitting A = A + + A~ , A ± = XA ± X~ 1 where A 1 * 1 is computed component wise. From this matrix splitting, the final upwind numerical flux function was constructed as g sw (u,v) = A + {y\u)u + A~(u;v)v . Although not part of their explicit construction, for the gas dynamic equations with ideal gas law, the jacobian matrix dg sv/ /du has eigenvalues that are all nonnegative and the Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 50 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS jacobian matrix dg sw /dv has eigenvalues that are all nonpositive whenever the ratio of specific heats 7 lies in the interval [1, 5/3]. The matrix splitting leads to numerical fluxes that do not vary smoothly near sonic and stagnation points. Use of the Steger- Warming flux splitting in the schemes of Sect. 2 and 3 results in rather poor resolution of linearly degenerate contact waves and velocity slip surfaces due to the introduction of excessive artificial diffusion for these wave families. • Van Leer flux vector splitting. Van Leer, 1982 provided an alternative flux splitting for the gas dynamic equations that produces a numerical flux of the form g VL (u,v) = f~(u)+f + (v) using special Mach number polynomials to construct fluxes that remain smooth near sonic and stagnation points. As part of the splitting construction, the jacobian matrix dg sw /du has eigenvalues that are all nonnegative and the matrix dg sw /dv has eigenvalues that are all nonpositive. The resulting expressions for the flux splitting are somewhat simpler when compared to the Steger- Warming splitting. The van Leer splitting also introduces excessive diffusion in the resolution of linearly degenerate contact waves and velocity slip surfaces. • System Lax-Friedrichs flux. This numerical flux is the system equation counterpart of the scalar Lax-Friedrichs flux (27). For systems of conservation laws the Lax-Friedrichs flux is given by 9 LF (u,v) = \{f{u) + /(«)) • v - |a(«/) (v - u) where a(u) is given through the eigenvalues Afc(^; w) of A(u; w) a{v) = max sup |Ajfc(i^; -il»)| . 1 < k < m w€[u,v] The system Lax-Friedrichs flux is usually not applied on the boundary of domains since it generally requires an over specification of boundary data. The system Lax-Friedrichs flux introduces a relatively large amount of artificial diffusion when used in the schemes of Sect. 2. Consequently, this numerical flux is typically only used together with relatively high order reconstruction schemes where the detrimental effects of excessive artificial diffusion are mitigated. • Harten-Lax-van Leer flux (Harten, Lax and van Leer, 1983). The Harten-Lax-van Leer numerical flux originates from a simplified two wave model of more general m wave systems such that waves associated with the smallest and largest characteristic speeds of the m wave system are always accurately represented in the two wave model. The following numerical flux results from this simplified two wave model g nLL (u,v) = -(f( u ) + f(v))-v--- — — (f(v) - /(«)) • v + — - — (v-u) Z Z u m ax ^min C*max "min where amax(") = max (0, sup \ k {v,w)) , a min = min (0, inf \ k (i>;w)) . l<k<m w£[u,v\ \<k<m ui€[«,u] When compared to the Lax-Friedrichs flux, this flux can be considerably more accurate in flow situations where < |(a max + Oj m i n )/(o! ra ax — <2min)| < 1- Using this flux, full upwinding is obtained for supersonic flow. Modifications of this flux are suggested in Einfeldt et al., 1998 to improve the resolution of intermediate waves as well. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley k Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 51 Further examples of numerical fluxes (among others) include the kinetic flux vector splitting due to Deshpande, 1986, the advection upstream splitting flux (AUSM) of Liou and Steffen, 1993, and the convective upwind and split pressure (CUSP) flux of Jameson, 1993 and Tatsumi et al., 1994. 5. Concluding Remarks The literature associated with the foundation and analysis of the finite volume methods is extensive. This article gives a very brief overview of finite volume methods with particular emphasis on theoretical results that have significantly impacted the design of finite volume methods in everyday use at the time of this writing. More extensive presentations and references on various topics in this article can be found in the books by Godlewski and Raviart, 1991, Kroner, 1997, Eymard, Galluoet and Herbin, 2000 and LeVeque, 2002. , REFERENCES Abgrall R. On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation. J. Comp. Phys. 1994; 114:45-58. Angermann L, Knabner P and Thiele K. An error estimate for a finite volume discretization of density driven flow in porous media. Appl. Numer. Math. 1998; 26:179-191. Bank R and Rose DJ. Some error estimates for the box method. SIAM J. Numer. Anal. 1987; 24:777-787. Barth TJ and Jespersen DC. The design and application of upwind schemes on unstructured meshes. American Institite for Aeronautics and Astronautics 1989; Report 89-0366:1-12. Barth TJ and Frederickson PO. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. American Institite for Aeronautics and Astronautics 1990; Report AIAA-90-0013. Batten P, Lambert, C and Causon DM. Positively Conservative high-resolution convection schemes for unstructured elements. Int. J. Numer. Meth. Engrg 1996; 39:1821-1838. Billey V, Periaux J, Perrier P and Stoufflet B. 2-D and 3-D Euler computations with finite element methods in aerodynamics. Lecture Notes in Mathematics (Vol. 1270), Springer- Verlag: Berlin, 1987. Boris JP and Book DL. Flux corrected transport: SHASTA, a fluid transport algorithm that works. J. Comp. Phys. 1973; 11:38-69. Bouchut F and Perthame B. Kruzkov's estimates for scalar conservation laws revisited. Trans. Am. Math. Soc. 1998; 350(7) :2847-2870. Carrillo J. Entropy solutions for nonlinear degenerate problems. Arch. Ration. Mech. Anal. 1999; 147:269-361. Cai Z. On the finite volume element method. Numer. Math. 1991; 58:713-735. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 52 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Chatzipantelidis P. A finite volume method based on the Crouzeix-Raviart element for elliptic problems. Numer. Math. 1999; 82:409-432. Chou SH and Li Q. Error estimates in L 2 , H 1 and L°° in covolume methods for elliptic and parabolic problems: a unified approach. Math. Comp. 2000; 69:103-120. Cockburn B, Coquel F and Lefioch PG. An error estimate for finite volume methods for multidimensional conservation laws. Math. Comput. 1994; 63:77-103. Cockburn B and Gau H. A posteriori error estimates for general numerical methods for scalar conservation laws. Comput. Appl. Math. 1995; 14:37-47. Cockburn B, Hou S and Shu CW. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case. Math. Comp. 1990; 54(190):545-581. Cockburn B, Lin SY and Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. III. One-dimensional systems. J. Comput. Phys. 1989; 84(1):90-113. Cockburn B and Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Math. Comp. 1989; 52:411- 435. Cockburn B and Shu CW. The Runge-Kutta discontinuous Galerkin method for conservation laws. V. Multidimensional systems. J. Comput. Phys. 1998; 141 (2): 199-224. Colella P and Woodward P. The piecewise parabolic methods for gas-dynamical simulations, J. Comp. Phys. 1984; 54:174-201. Chainais-Hillairet C. Finite volume schemes for a nonlinear hyperbolic equation: convergence towards the entropy solution and error estimates. M2AN Math. Model. Numer. Anal. 1999; 33:129-156. Chainais-Hillairet C. Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimates. Math. Methods Appl. Sci. 2000; 23(5):467-490. Cournede P-H and Debiez C and Dervieux A. A positive MUSCL scheme for triangulations. Institut National De Recherche En Informatique Et En Automatique (INRIA), Report 3465, 1998. DiPerna RJ. Measure-valued solutions to conservation laws. Arch. Rational Mech. Anal. 1985; 88(3):223-270. Delanaye M. Polynomial Reconstruction Finite Volume Schemes for the Compressible Euler and Navier-Stokes Equations on Unstructured Adaptive Grids. Ph.D. Thesis 1996, University of Liege, Belgium. Deshpande, SM. On the Maxwellian distribution, symmetric form, and entropy conservation for the Euler equations. NASA Langley, Hampton, Virginia 1986; NASA Report TP-2583. Desideri JA and Dervieux A. Compressible flow solvers using unstructured grids. Von Karman Institute Lecture Notes 1988-05; Von Karman Institute for Fluid Dynamics, Belgium. Einfeldt B, Munz C, Roe P and Sjogreen B. On Godunov-type methods near low densities. J. Comput. Phys. 1992; 92:272-295. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 53 Ewing RE, Lin T and Lin Y. On the accuracy of the finite volume element method bases on piecwise linear polynomials. SIAM J. Numer. Anal. 2002; 39(6):1865-1888. Eymard R, Gallouet T, Ghilani M, and Herbin R. Error estimates for the approximate solution of a nonlinear hyperbolic equation given by finite volume schemes. IMA Journal of Numerical Analysis 1998; 18:563-594. Eymard R, Galluoet T and Herbin R. Finite volume methods. Handbook of Numerical Analysis, North Holland: Amsterdam, 2000; 7:713-1020. Eymard R, Gallouet T and Herbin R. Finite volume approximation of elliptic problems and convergence of an approximate gradient. Appl. Numer. Math. 2001; 37(l-2):31-53. Eymard R, Gallouet T and Herbin R. Error estimates for approximate solutions of a nonlinear convection-diffusion problem. Adv. Differential Equations 2002; 7(4):419-440. Eymard R, Gallouet T, Herbin R. and Michel A. Convergence of a finite volume scheme for nonlinear degenerate parabolic equations. Numer. Math. 2002; 92(l):41-82. Feistauer M, Felcman J, Lukacova-Medvid'ova M, and Warnecke G. Error estimates for a combined finite volume-finite element method for nonlinear convection-diffusion problems. SIAM J. Numer. Anal. 1999; 36(5):1528-1548. Gallouet T, Herbin R and Vignal MH. Error estimates on the approximate finite volume solution of convection diffusion equations with general boundary conditions. SIAM J. Numer. Anal. 2000; 37(6):1935-1972. Godunov SK. A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics. Mat. Sb. 1959; 47:271-290. Godlewski E and Raviart P-A. Hyperbolic systems of conservation laws. Mathematiques & Applications Ellipses, Paris 1991. Goodman JD and LeVeque RJ. On the accuracy of stable schemes for 2D conservation laws. Math. Comp. 1985; 45(171):15-21. Gottlieb S and Shu CW. Total variation diminishing Runge-Kutta schemes. Math. Comput. 1998; 67(221):73-85. Gottlieb S, Shu CW and Tadmor E. Strong stability-preserving high-order time discretization methods. SIAM Rev. 2001; 43(1):89-112. Harten A, Hyman JM, and Lax PD. On finite-difference approximations and entropy conditions for shocks. Comm. Pure and Appl. Math. 1976; 24:297-322. Harten A. High resolution schemes for hyperbolic conservation laws. J. Comp. Phys. 1983; 49:357-393. Harten A, Lax PD and van Leer, B. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 1983; 25:35-61. Harten A, Osher S, Engquist B and Chakravarthy S. Some results on uniformly high order accurate essentially non-oscillatory schemes. Appl. Num. Math. 1986; 2:347-377. Harten A, Osher S, Engquist B and Chakravarthy S. Uniformly high-order accurate essentially nonoscillatory schemes III. J. Comp. Phys. 1987; 71(2):231-303. Harten A. ENO schemes with subcell resolution. J. Comp. Phys. 1989; 83:148-184. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. 54 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS Harten A and Chakravarthy S. Multi-dimensional ENO schemes for general geometries. Institite for Computer Applications in Science and Engineering 1991; Report ICASE- 91-76. Herbin R and Ohlberger M. A posteriori error estimate for finite volume approximations of convection diffusion problems. In proceedings -.Finite volumes for complex applications - problems and perspectives, Porquerolles, 753-760; Hermes Science Publications, Paris, 2002. Hermeline F. A finite volume method for the approximation of diffusion operators on distorted meshes. J. Comput. Phys. 2000; 160(2):481-499. Jaffre J, Johnson C and Szepessy A. Convergence of the discontinuous Galerkin finite element method for hyperbolic conservation laws. Math. Models and Methods in Appl. Sci. 1995; 5(3):367-386. Jakobsen ER and Karlsen KH. Continuous dependence estimates for viscosity solutions of fully nonlinear degenerate parabolic equations. J. Differential Equations 2002; 183(2):497-525. Jameson A and Lax PD. Conditions for the construction of multipoint variation diminishing difference schemes, Appl. Numer. Math. 1986; 2(3-5):335-345. Jameson A and Lax PD. Corrigendum: Conditions for the construction of multipoint variation diminishing difference schemes, Appl. Numer. Math. 1987; 3(3):289. Jameson A. Artificial Diffusion, Upwind biasing, limiters and their effect on accuracy and convergence in transonic and hypersonic flows. American Institite for Aeronautics and Astronautics 1993; Report AIAA-93-3359:l-28. Jiang G and Shu CW. Efficient implementation of weighted ENO schemes. J. Comp. Phys. 1996; 126:202-228. Johnson C and Szepessy A. Adaptive finite element methods for conservation laws based on a posteriori error estimates. Commun. Pure Appl. Math. 1995; 48:199-234. Karlsen KH and Risebro NH. On the uniqueness and stability of entropy solutions of nonlinear degenerate parabolic equations with rough coefficients. Preprint 143, Department of Mathematics, University of Bergen, 2000. Koren, B. Upwind schemes for the Navier-Stokes equations. Proceedings of the Second International Conference on Hyperbolic Problems. Vieweg: Braunschweig, 1988. Kroner D. Numerical Schemes for Conservation Laws. Wiley-Teubner: Stuttgart, 1997. Kroner D, Noelle S and Rokyta M. Convergence of higher order upwind finite volume schemes on unstructured grids for conservation laws in several space dimensions. Numer. Math. 1995; 71:527-560. Kroner D and Ohlberger M. A-posteriori error estimates for upwind finite volume schemes for nonlinear conservation laws in multidimensions. Math. Comput. 2000; 69:25-39. Kruzkov SN. First order quasilinear equations in several independent variables. Math. USSR Sbornik 1970; 10:217-243. Kiither M. Error estimates for second order finite volume schemes using a TVD-Runge- Kutta time discretization for a nonlinear scalar hyperbolic conservation law. East-West J. Numer. Math. 2000; 8(4):299-322. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd. FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 55 Kuznetsov NN. Accuracy of some approximate methods for computing the weak solutions of a first-order quasi-linear equation. USSR, Comput. Math, and Math. Phys. 1976; 16(6):159— 193. Lax PD. Hyperbolic Systems of Conservation Laws. SIAM Pub.rPhiladelphia, 1973. LaxPD and WendroffB. Systems of conservation laws. Comm. Pure Appl. Math. 1960; 13:217- 237. , Lazarov RD, Michev ID and Vassilevsky PS. Finite volume methods for convection-diffusion problems. SIAM J. Numer. Anal. 1996; 33:31-35. LeVeque R. High resolution finite volume methods on arbitrary grids via wave propagation. J. Comp. Phys. 1988; 78:36-83. LeVeque R. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press: Cambridge, 2002. Liou MS and Steffen CJ. A new flux-splitting scheme. J. Comp. Phys. 1993; 107:23-39. Liu X-D. A maximum principle satisfying modification of triangle based adaptive stencils for the solution of scalar hyperbolic conservation laws. SIAM J. Numer. Anal. 1993; 30:701- 716. Malek J, Necas J, Rokyta M, and Ruzicka M. Weak and measure-valued solutions to evolutionary PDEs. Applied Mathematics and Mathematical Computation (Vol 13). Chapman and Hall: London, 1968. Ohlberger M. A posteriori error estimates for finite volume approximations to singularly perturbed nonlinear convection-diffusion equations. Numer. Math. 2001a; 87(4):737-761. Ohlberger M. A posteriori error estimates for vertex centered finite volume approximations of convection-diffusion-reaction equations. M2AN Math. Model. Numer. Anal. 2001b; 35(2):355-387. Oleinik OA. Discontinuous solutions of non-linear differential equations. Amer. Math. Soc. Transl. (2) 1963; 26:95-172. Osher S and Solomon F. Upwind Difference Schemes for Hyperbolic Systems of Conservation Laws. Math. Comp. 1982; 38(158) :339-374. Osher S. Riemann solvers, the entropy condition, and difference approximations. SIAM J. Numer. Anal. 1984; 21(2):217-235. Osher S. Convergence of generalized MUSCL schemes. SIAM J. Numer. Anal. 1985; 22(5):947- 961. Peterson T. A Note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation. SIAM J. Numer. Anal. 1991; 28(1):133-140. Roe PL. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comp. Phys. 1981; 43:357-372. Rostand P and Stoufflet B. TVD schemes to compute compressible viscous flows on unstructured meshes. Proceedings of the Second International Conference on Hyperbolic Problems. Vieweg: Braunschweig, 1988. Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. © 2004 John Wiley & Sons, Ltd.