UUnderstand Slope Limiter - Graphically
Ling ZouNuclear Science and Engineering Division, Argonne National Laboratory9700 S. Cass Ave, Lemont, IL 60439 ([email protected])
Abstract.
In this article, we illustrate how the concept of slope limiter can be interpretedgraphically, i.e., how the slope of reconstructed piecewise linear function is limited by four boundinglines that connect cell-averaged data and its neighboring cell-averaged data. It is then conjecturedthat the same graphical rule can be generalized from uniform mesh to non-uniform mesh, suchthat the high-resolution total variance diminishing (TVD) region of slope limiter for non-uniformmeshes can be obtained.
Keywords.
Slope limiter, high-resolution TVD, irregular mesh
Slope limiter is central to high-resolution total variance diminishing (TVD) finite volume meth-ods, which are widely used in the numerical solving of hyperbolic partial different equations. Thesemethods can produce high-order spatial accuracy in regions with smooth solutions, and can avoidspurious oscillations near discontinuities, such as shock waves. In one dimension, this is realized byperforming piecewise linear data reconstruction in each finite volume cell, such that better-qualityedge values can be obtained to compute numerical fluxes on the interface between neighboring cells.Owing to the early contributions from van Leer [1], Sweby [2], and many other researchers,mathematical theories have been developed to form the concept of flux/slope limiter. Among theseearlier groundbreaking researches, Sweby’s diagram has been widely used to illustrate the second-order TVD region of slope limiters. These theories are mathematically rigorous, however they arenot physically straightforward or intuitive, especially to people with engineering background (likeme) instead of mathematics background. Inspired by the work of Berger et al. [3], who provideda different perspective to examine the properties of slope limiters, we found that the concept canbe further interpreted in a more straightforward and physically intuitive way using graphics - themain idea discussed in this article.In the following sections, we will at first have a very short introduction to the finite volumemethod to layout the context where the slope limiter concept will be discussed. Following that,we will then discuss the concept of slope limiter, and how it can be interpreted graphically. Inthe last section, we propose a conjecture that will extend high-resolution TVD slope limiters fromuniform mesh to non-uniform mesh, based on which the second-order TVD region of slope limiterfor non-uniform meshes can be obtained. These methods are commonly referred to as high-resolution instead of second-order methods, as they drop tofirst-order accuracy at local extrema. a r X i v : . [ m a t h . NA ] F e b Finite Volume Method
Considering the following simple hyperbolic partial differential equation u t + f x = 0 , x ∈ ( −∞ , ∞ ) f = au, (1)where u t ≡ ∂u ( x, t ) /∂t , f x ≡ ∂f ( x, t ) /∂x , and a is a constant non-zero real number. This equationis also called the linear advection equation, and can be rewritten in a simpler form as u t + au x = 0 , (2)which could be interpolated as the advection of a tracer material, with concentration u , in a fluidfield with flow speed a .Let’s solve the above hyperbolic equation (1) using the finite volume method on a mesh with auniform cell size, as illustrated in figure 1. The i-th grid cell is denoted by C i = ( x i − / , x i +1 / ) , and we use U ni to approximate the average value of u over the i-th interval at time t n U ni ≈ x (cid:90) x i +1 / x i − / u ( x, t n ) dx ≡ x (cid:90) C i u ( x, t n ) dx (3)in which, ∆ x = x i +1 / − x i − / , which is constant for the uniform-sized mesh. x i-1/2 x i+1/2 t n t n+1 𝑈 " 𝑈 " $% 𝑈 " &% 𝑈 " 𝐹 " $%/) 𝐹 " &%/) Figure 1: Illustration of a finite volume method in the x - t space.Assuming that U n is already known, we wish to solve for U n +1 . Consider the integral form ofthe conservation law, equation (1), on the i-th interval, ddt (cid:90) C i u ( x, t ) dx = f ( u ( x i − / , t )) − f ( u ( x i +1 / , t )) (4)which is obtained by applying Leibniz integral rule and the divergence theorem (Gauss’s theorem).Then, performing a time integral from t n to t n +1 on the above equation, it is easy to see that1∆ x (cid:90) C i u ( x, t n +1 ) dx = 1∆ x (cid:90) C i u ( x, t n ) dx − x (cid:20)(cid:90) t n +1 t n f ( u ( x i +1 / , t )) dt − (cid:90) t n +1 t n f ( u ( x i − / , t )) dt (cid:21) . (5)Now let’s introduce a numerical flux, F , to approximate the average flux along the cell edge, i.e., F ni − / ≈ t (cid:90) t n +1 t n f ( u ( x i − / , t )) dt. (6)2 !" 𝑈 ! 𝑈 !%$" (a) a Δ t (b) 𝑈 !" 𝑈 !%$" 𝑈 ! (c) Figure 2: (a) Reconstruct; (b) Evolve; (c) Average.If we can approximate this numerical flux based on the values U n (for explicit method), we nowhave a fully discrete method to solve for U n +1 , i.e., U n +1 i = U ni − ∆ t ∆ x (cid:16) F ni +1 / − F ni − / (cid:17) . (7) For this advection equation, the simplest way to compute the numerical flux is the so-calleddonor-cell upwind method. Without loss of generality, assuming the flow is in the positive x direction, i.e., a >
0. For the donor-cell upwind method, the flux on the edges of each cell isentirely determined by its upwind neighboring cell. For example, on the left edge of the i-th cell,the flux is determined as F ni − / = aU ni − , (8)and similarly on the right edge, F ni +1 / = aU ni . (9)Using this upwind method, the discrete equation becomes U n +1 i = U ni − a ∆ t ∆ x (cid:0) U ni − U ni − (cid:1) , (10)from which U n +1 i can be updated. This process can be illustrated as the process shown in figure 2.We start from the known cell-averaged solutions, U n , from time step t n , and then perform a simplepiecewise constant reconstruction of the solution in each cell, shown as figure 2a. We then evolvethe solution using the upwind method within the time step, ∆ t , i.e., the entire solution profilemoves to the right (as a >
0) for a distance of a ∆ t , figure 2b. Finally, perform local average ineach cell to obtain the solution U n +1 at t n +1 , shown as figure 2c.This process can be seen as a simplified version of the reconstruct-evolve-average (REA) algo-rithm, (see section 4.10 of LeVeque [4]). We quote the REA algorithm and rewrite it as: REA Algorithm: Reconstruct a piecewise polynomial function defined for all x from the cell-averaged values.In the case demonstrated above, a picecwise constant function was reconstructed.2. Evolve the hyperbolic equation exactly (or approximately) with this reconstructed functionfor a time step ∆ t .3. Average the evolved function over each cell to obtain the new cell-averaged value for thenew time step. 3
Slope Limiter
In the previous section, in the context of finite volume method/REA algorithm, we discussedthe simplest reconstruction method, namely the piecewise constant reconstruction in each cell. Itis however well-understood that such an upwind scheme is only first-order. To achieve higher-order(e.g. second-order) spatial accuracy, instead of piecewise constant reconstruction, piecewise linearreconstruction is needed.To maintain the total “mass” in each cell, e.g. the total amount of tracer material, the recon-structed piecewise linear function must pass through, for example, the point ( x i , U i ) in the i-thcell, marked as ‘x’ in figure 3a. So this really leaves only one freedom we can play with, the slopeof the piecewise linear function. As illustrated in figure 3b, to perform linear function reconstruc-tion, it is like that the center of the bar is pinned, while the bar could be rotated clockwise orcounterclockwise (the ‘hand’ symbol) about this pinned center. 𝑈 !" 𝑈 ! 𝑈 !%$" x (a) 📌 E (b) Figure 3: (a) Illustration of the piecewise linear reconstruction within each cell; (b) Illustration ofhow piecewise linear reconstruction is performed.If the local “mass” conservation is the only constraint, the slope can take any values. However,we will see that the slope cannot take arbitrarily values, if we would like to construct a (high-resolution) TVD scheme. For TVD schemes, the slope can only be chosen from a limited range,that is how slope limiter got its name. We do not repeat the concept of TVD here, which canbe found in textbooks on finite volume method, such as LeVeque [4]. In general, TVD methodcan be interpreted as a monotonicity-preserving method, or more intuitively, a method that avoidsnon-physical overshoot/undershoot solutions. 𝑈 " 𝑈 " 𝑥 " ∆ ∆ ’ ∆ ( 𝑥 " 𝑥 "’$ 𝑈 " ’$ (a) 𝑥 " 𝑠 = ∆ ∆𝑥 𝑥 " 𝑥 "($ 𝑠 ( = ∆ ( ∆𝑥2𝑠 ( (b) 𝑥 " 𝑠 & = ∆ ) " 𝑥 "+$ (c) Figure 4: (a) The differences between neighboring cell-averaged values; (b) Slopes relevant toconstruct slope limiter; (c) The reference slope; 4efore we discuss the slope limiter, several important quantities ought to be introduced anddiscussed. First, the differences between neighboring cell-averaged values are defined as:∆ − := U i − U i − ∆ + := U i +1 − U i , (11)and ∆ t = ∆ − + ∆ + , (12)which are illustrated in figure 4a. In above equations, we use := for ‘equal by definition’. Itis noted that their values can be positive or negative. We then introduce the non-dimensionallocation indicator, f , defined as f := ∆ − ∆ t (13)which indicates the relative location of U i in between U i − and U i +1 . It is noted that, if U i − , U i , and U i +1 are monotonic, either increasing or decreasing, we have 0 ≤ f ≤
1. For all otherconditions, U i is a local extremum, and we have −∞ < f < < f < + ∞ . It is also notedthat, traditionally, the indicator is expressed as r := ∆ − / ∆ + (e.g., Sweby [2]), which makes theconcept of slope limiter less physically intuitive.The two slopes associated with the ‘-’ (left) and ‘+’ (right) sides are defined as: s − := U i − U i − ∆ x = ∆ − ∆ xs + := U i +1 − U i ∆ x = ∆ + ∆ x , (14)which are illustrated in figure 4b as the blue and the green lines, respectively. In the same figure,the red line has a slope of 2 s − , and the orange line has a slope of 2 s + . We will see that, these fourbounding slopes (lines) play the utmost important roles to construct high-resolutionTVD slope limiters .In figure 4c, a reference slope is also defined, which is the slope of the line that connects thetwo neighboring data, i.e., ( x i − , U i − ) and ( x i +1 , U i +1 ), s R := U i +1 − U i − x i +1 − x i − = ∆ t x . (15)Using s R , we can non-dimensionalize the slopes discussed above. It is easy to see: φ − := s − s R = 2 f ,φ + := s + s R = 2(1 − f ) , (16)in which φ − and φ + are the non-dimensional slopes of the blue and the green lines, respectively.The red line has a non-dimensional slope of 4 f ; and the orange line has a non-dimensional slope of4(1 − f ).There are several special cases when two of the four slopes are equal, which are illustrated infigure 5. The special values of f will define the intervals where the slope limiter behaves differently.Finally, let s be the slope of the reconstructed piecewise linear function in the i-th cell. Wewould like to learn how s should behave when U i is at different relative locations in between U i − " 𝑥 " 𝑥 "%$ 𝑠 % = 2𝑠 (a) s + = 2 s − , when f = 1 / 𝑥 " 𝑥 " 𝑥 "%$ 𝑠 % = 𝑠 (b) s + = s − , when f = 1 / 𝑥 " 𝑥 " 𝑥 "%$ 𝑠 = 2𝑠 % (c) s − = 2 s + , when f = 2 / Figure 5: Several special cases when two of the four slopes are equal. It is noted that cases (a) and(c) are symmetric about (b).and U i +1 . Let’s define the non-dimensional slope of the reconstructed piecewise linear function, φ ,as φ := ss R . (17)It is now equivalent to ask: how φ should behave given different f . We will learn that, to achievehigh-resolution TVD scheme, φ can only be expressed as a nonlinear function of f , i.e., φ = φ ( f ),which will be dictated by the four colored lines shown in figure 4b. The rule for slope limiters to achieve TVD is very simple. The mathematical theory behind therules is well explained in Sweby [2], which we do not intend to repeat. Instead, we will illustratethe two important rules in a graphical way.
TVD rules:
1. If U i is a local extremum, −∞ < f < < f < + ∞ , the piecewise linear reconstructionmust have a zero-slope, i.e., piecewise constant, as illustrated in figure 6a.2. For all other conditions (0 ≤ f ≤ U i is in between U i − and U i +1 .First, the reconstructed slope must not be negative (positive) for monotonically increasing(decreasing) data, as illustrated by the ‘lock’-A symbols in figures 6b and 6c.Second, the edge values of the linear function must not go beyond the neighboring cell-averaged values. This is illustrated in figure 6b, when U i is closer to U i − , the reconstructedlinear function must not go beyond the red line to allow its left-edge value to pass U i − . Thisis illustrated by the ‘!’ symbol and the ‘lock’-B symbol in figure 6b. Similarly, figure 6cillustrates the case when U i is closer to U i +1 .In both figures 6b and 6c, the TVD region is simply bounded by the zero-slope blue line andthe red/orange line. Anything between them, marked by the green arrow, is acceptable as aTVD scheme.To summarize, for TVD schemes, the non-dimensional slope of the reconstructed piecewise6 !" 𝑈 ! 𝑥 !" 𝑥 ! 𝑥 !$ 𝑈 !$ 📌 🔒 (a) 𝑥 !" 𝑥 ! 𝑥 !$ 📌 ❗ 🔒 A 🔒 B 𝑈 !" 𝑈 ! 𝑈 !$ (b) 𝑥 !" 𝑥 ! 𝑥 !$ 📌 ❗ 🔒 A 🔒 B 𝑈 !" 𝑈 !$ (c) Figure 6: The limiting (a) U i is a local extremum; (b) U i is in between U i − and U i +1 , closer to U i − ; (c) U i is in between U i − and U i +1 , closer to U i +1 .linear function will respect the following conditions: ≤ φ ( f ) ≤ f ; for 0 ≤ f ≤ . ≤ φ ( f ) ≤ − f ); for 0 . < f ≤ φ ( f ) = 0; otherwise ; (18)or in a more compact form: (cid:40) ≤ φ ( f ) ≤ min (4 f, − f )) ; for 0 ≤ f ≤ φ ( f ) = 0; otherwise . (19)On the φ - f plot, the TVD region is illustrated as the shaded area in figure 8a. In the previous section, we have discussed the slope limiter region for being TVD. If we wouldlike to further improve the slope limiter to achieve higher order spatial algorithm, clearly, it willonly be a subset of the TVD scheme. The mathematical derivation of high-resolution TVD schemescan also be found in Sweby [2]. Again, we will proceed our discussion on a set of monotonicallyincreasing data, such that the four slopes, s − , 2 s − , s + , and 2 s + , are all positive.The rule for high-resolution TVD is fascinatingly simple: High-resolution TVD slopes must lay in between the two smallest among the fourbounding slopes.
Figure 7 illustrates the two smallest slopes for various f values. Especially, figure 7c shows aspecial case where the two smallest slopes are identical, and also equal to the reference slope. Onthe φ - f plot, the high-resolution TVD region can be easily found in between the two lines havingsmallest values among the four. This is illustrated as the shaded area in figure 8b. It is noted thatany high-resolution TVD schemes must pass the point: φ (cid:18) (cid:19) = 1 , (20)7 " 𝑥 " 𝑥 "%$ 𝑠 (a) 0 < f < / 𝑥 " 𝑥 " 𝑥 "%$ 𝑠 𝑠 % (b) 1 / < f < / 𝑥 " 𝑥 " 𝑥 "%$ 𝑠 ’ = 𝑠 % = 𝑠 (c) f = 1 / 𝑥 " 𝑥 " 𝑥 "%$ 𝑠 % 𝑠 (d) 1 / < f < / 𝑥 " 𝑥 " 𝑥 "%$ % 𝑠 % (e) 2 / < f < Figure 7: Illustration of how slope should be limited to achieve high-order TVD for different condi-tions. For all cases, the acceptable reconstructed linear function must lay between the two coloredlines, which have the smallest slopes among the four lines. For condition (c), the reconstructedlinear function must have a slope of s R .i.e., the special case shown in figure 7c, which means that, if the solution is already a linear function,the reconstructed linear function must follow the same linear function to maintain the linearity.This is the special case as illustrated in figure 7c. f0 1/2 11234 ɸ ( f ) TVD (a) f0 1/2 11234 ɸ ( f )
13 23 (b)
Figure 8: Illustration of (a) the TVD region (shaded); and (b) the high-resolution TVD region(shaded) on the φ - f plot. It is clear that the high-resolution TVD region is a subset of the TVDregion. Symmetry:
It is noted that some high-resolution TVD scheme are so-called being symmetric.This is can be illustrated in figure 9, in which conditions (a) and (b) are symmetric about f = 1 / φ (1 − f ) = φ ( f ) . (21)Many well-known high-resolution TVD schemes do possess this symmetric feature. Figure 10 showsseveral such examples, which also include the sin slope limiter proposed by Berger et al. [3] to8emonstrate how a new symmetric slope limiter could be easily constructed. However, we notethat being symmetric is NOT required to be a high-resolution TVD scheme. 𝑥 " 𝑥 " 𝑥 "%$ ∆ (a) 𝑥 " 𝑥 " 𝑥 "%$ ∆ % (b) Figure 9: Illustration of a symmetric condition about f = 1 /
2, ∆ − ( a ) = ∆ + ( b ) and ∆ t ( a ) = ∆ t ( b ).The mathematical form of these example limiters are given as follows. • superbee , the upper bound of the high-resolution TVD region φ ( f ) = f ≤ f ≤ / − f ) 1 / < f ≤ / f / < f ≤ / − f ) 2 / < f ≤ • minmod , the lower bound of the high-resolution TVD region φ ( f ) = min[2 f, − f )] (23) • Barth-Jespersen , also known as MC φ ( f ) = min[1 , f, − f )] (24) • van Leer φ ( f ) = 4 f (1 − f ) (25) • van Albada φ ( f ) = 2 f (1 − f ) f + (1 − f ) (26) • sin φ ( f ) = sin( πf ) (27)9 .00.51.01.5 0 0.25 0.5 0.75 1 ɸ ( f ) f superbeeminmodBarth-Jespersenvan Leervan Albadasin Figure 10: Illustration of several symmetric high-resolution TVD slope limiters.
To the author’s best knowledge, there do not exist a complete theoretical derivation to de-fine high-resolution TVD region on non-uniform meshes. Berger et al. [3] provided a theoreticalderivation to define the upper bound of TVD region and the generalized min slope limiter forhigh-resolution schemes on non-uniform meshes. They also specified the necessary condition whendealing with linear data (see f in equation (31)).In the previous section, we have interpreted the high-resolution TVD region on uniform meshesgraphically using the bounding slopes. It seems natural that the same rule can be applied toirregular meshes, which is proposed as the following conjecture: Conjecture 1:Uniform mesh is a special case of the non-uniform mesh, and thus the high-resolutionTVD rule for uniform meshes can be generalized and applied to non-uniform meshes.
To illustrate this, the reference slope and non-dimensional slopes need to be defined first. Fol-lowing Berger et al. [3], we first define the mesh stretching ratios: a := ∆ x i − ∆ x i b := ∆ x i +1 ∆ x i (28)10 𝑥 𝑠 ’ Δ𝑥 Δ𝑥 (a) Δ𝑥 𝜙 $ Δ𝑥 Δ𝑥 𝜙 ’ 𝜙 ()*+ 𝜙 , (b) Figure 11: (a) The reference slope on a irregular mesh; (b) Differences slopes on a irregular mesh.Same as the uniform mesh, the reference slope s R is defined as the slope of the line that connects( x i − , U i − ) and ( x i +1 , U i +1 ): s R := U i +1 − U i − x i +1 − x i − = 2∆ t (2 + a + b )∆ x i , (29)as illustrated in figure 11a. The four bounding lines are shown in figure 11b, and their non-dimensional slopes are: φ − := s − s R = 2 + a + b a fφ left := s left s R = (2 + a + b ) fφ + := s + s R = 2 + a + b b (1 − f ) φ right := s right s R = (2 + a + b )(1 − f ) (30)in which, f is the same as defined for uniform-sized mesh, and f := ∆ − / ∆ t . Similar to the uniformmesh case, at several special f values, two of these four slopes are equal. These special f valuesare given as follows: f = 12 + b when φ left = φ + f = 1 + a a + b when φ − = φ + f = 1 + a a when φ − = φ right (31)The direct outcome of Conjecture 1 is that, on non-uniform meshes, the high-resolution TVDslope is restricted by the two smallest slopes among the four colored lines. On the φ - f plot, this isthe shaded region illustrated in figure 12b. In the same figure, the high-resolution TVD region onuniform mesh is also shown for comparison. It is easy to see that figure 12b reduces to figure 12aon a uniform mesh, where a = b = 1.At this stage, it is tempting to obtain high-resolution TVD scheme on non-uniform meshes byprojecting existing high-resolution TVD schemes on uniform meshes. This can be done by findingthe transform matrix that will homogeneously transform the two triangles of figure 12a into thetwo in 12b. The transformed high-resolution TVD schemes can then be obtained by using this11ransform matrix. This however does not seem to be the most efficient way, and there does notseem to have direct benefit to do so, e.g., it is difficult to define symmetry on the irregular mesh,and the transformation of symmetric slope limiter does not seem to be very useful.It is noted that figure 12b is equivalent to figure 4 of Berger et al. [3], except that 1) Berger etal. used a different reference slope, so in our plot, the y-axis is scaled by a factor of (2+ a + b ) /
2; and2) Berger et al. only marked the TVD region, not the high-resolution TVD region as we illustratein figure 12b. Berger et al. [3] also proposed two functions that lay in the TVD region, which canbe seen as generalized van Leer slope limiters, as they recover the usual van Leer slope limiter onuniform meshes. It can also be easily proved that, with proper scaling of (2 + a + b ) /
2, they lay inthe high-resolution TVD region proposed in this article. To clarify the discussion, one of the slopelimiter provided by Berger et al. [3], i.e., their equation (38), is scaled and provided as: φ ( f ) = 2 + a + b f (cid:34) − a a (cid:18) ff (cid:19) /a (cid:35) f ≤ f (1 − f ) (cid:34) − b b (cid:18) − ff (cid:19) /b (cid:35) f > f (32)As discovered by Berger et al. [3], this slope limiter does a very good job on various non-uniformmeshes. The details on the performance of this slope limiter on non-uniform meshes are referred toBerger et al. [3]. The same slope limiter has also been used by the author in a previous numericalsimulation study of the Welander natural circulation problem, which showed very good results [5]. f0 1/2 112 4f 4(1-f)2f 2(1-f)1/3 2/3 ɸ ( f ) (a) f0 1/2 11 f f ɸ ( f ) f (b) Figure 12: Illustration of (a) the TVD region (shaded); and (b) the high-resolution TVD region(shaded) on the φ - f plot. It is clear that the high-resolution TVD region is a subset of the TVDregion. 12 Summary
In this paper, we have illustrated how the concept of slope limiter can be interpreted graphicallyfor uniform meshes. We then conjecture that the same graphical rule can be generalized for non-uniform meshes. The high-resolution total variance diminishing (TVD) region of slope limiter fornon-uniform meshes can then be obtained.Future work is needed to prove the conjecture, which seems possible if following Sweby’s deriva-tions for uniform meshes [2].
Acknowledgment
Argonne National Laboratory’s work was supported by the U.S. Department of Energy, Officeof Nuclear Energy, under contract DE-AC02-06CH11357.
References [1] van Leer, B., 1974. Towards the ultimate conservative difference scheme. II. Monotonicityand conservation combined in a second-order scheme.
Journal of computational physics ,14(4), pp.361-370.[2] Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservationlaws.
SIAM journal on numerical analysis , 21(5), pp.995-1011.[3] Berger M., Aftosmis M. J., Murman, 2005. Analysis of slope limiters on irregular grids.43rd AIAA Aerospace Sciences Meeting and Exhibit.[4] LeVeque, R. J., 2002. Finite-volume methods for hyperbolic problems. Cambridge Univer-sity Press.[5] Zou, L., Zhao, H., and Kim, S. J., 2017. Numerical study on the Welander oscillatory nat-ural circulation problem using high-order numerical methods.