Understand Slope Limiter - Graphically

Ling Zou

Nuclear Science and Engineering Division, Argonne National Laboratory

9700 S. Cass Ave, Lemont, IL 60439 ([email protected])
arXiv:2102.04435v1 [math.NA] 8 Feb 2021

Abstract. In this article, we illustrate how the concept of slope limiter can be interpreted
graphically, i.e., how the slope of reconstructed piecewise linear function is limited by four bounding
lines that connect cell-averaged data and its neighboring cell-averaged data. It is then conjectured
that the same graphical rule can be generalized from uniform mesh to non-uniform mesh, such
that the high-resolution total variance diminishing (TVD) region of slope limiter for non-uniform
meshes can be obtained.
Keywords. Slope limiter, high-resolution TVD, irregular mesh

1 Introduction
Slope limiter is central to high-resolution1 total variance diminishing (TVD) finite volume meth-
ods, which are widely used in the numerical solving of hyperbolic partial different equations. These
methods can produce high-order spatial accuracy in regions with smooth solutions, and can avoid
spurious oscillations near discontinuities, such as shock waves. In one dimension, this is realized by
performing piecewise linear data reconstruction in each finite volume cell, such that better-quality
edge 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 these
earlier 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 are
not physically straightforward or intuitive, especially to people with engineering background (like
me) instead of mathematics background. Inspired by the work of Berger et al. [3], who provided
a different perspective to examine the properties of slope limiters, we found that the concept can
be further interpreted in a more straightforward and physically intuitive way using graphics - the
main idea discussed in this article.
In the following sections, we will at first have a very short introduction to the finite volume
method 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. In
the last section, we propose a conjecture that will extend high-resolution TVD slope limiters from
uniform mesh to non-uniform mesh, based on which the second-order TVD region of slope limiter
for non-uniform meshes can be obtained.
1
These methods are commonly referred to as high-resolution instead of second-order methods, as they drop to
first-order accuracy at local extrema.

1
2 Finite Volume Method
Considering the following simple hyperbolic partial differential equation
ut + fx = 0, x ∈ (−∞, ∞)
(1)
f = au,
where ut ≡ ∂u(x, t)/∂t, fx ≡ ∂f (x, t)/∂x, and a is a constant non-zero real number. This equation
is also called the linear advection equation, and can be rewritten in a simpler form as
ut + aux = 0, (2)
which could be interpolated as the advection of a tracer material, with concentration u, in a fluid
field with flow speed a.
Let’s solve the above hyperbolic equation (1) using the finite volume method on a mesh with a
uniform cell size, as illustrated in figure 1. The i-th grid cell is denoted by
Ci = (xi−1/2 , xi+1/2 ),
and we use Uin to approximate the average value of u over the i-th interval at time tn
Z xi+1/2 Z
n 1 1
Ui ≈ u(x, tn )dx ≡ u(x, tn )dx (3)
∆x xi−1/2 ∆x Ci

in which, ∆x = xi+1/2 − xi−1/2 , which is constant for the uniform-sized mesh.

𝑈"#&%
tn+1
𝐹"#$%/) 𝐹"#&%/)

𝑈"#$% 𝑈"# 𝑈"#&%

tn
xi-1/2 xi+1/2

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 of
the conservation law, equation (1), on the i-th interval,
Z
d
u(x, t)dx = f (u(xi−1/2 , t)) − f (u(xi+1/2 , t)) (4)
dt Ci
which is obtained by applying Leibniz integral rule and the divergence theorem (Gauss’s theorem).
Then, performing a time integral from tn to tn+1 on the above equation, it is easy to see that
Z Z
1 1
u(x, tn+1 )dx = u(x, tn )dx
∆x Ci ∆x Ci
Z tn+1 Z tn+1  (5)
1
− f (u(xi+1/2 , t))dt − f (u(xi−1/2 , t))dt .
∆x tn tn

Now let’s introduce a numerical flux, F , to approximate the average flux along the cell edge, i.e.,
Z tn+1
n 1
Fi−1/2 ≈ f (u(xi−1/2 , t))dt. (6)
∆t tn

2
 "#$ 𝑈 "#$ 𝑈 "#$
𝑈!%$ ! !#$
"
𝑈!#$

𝑈!"

aΔt
"
𝑈!%$

(a) (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 now
have a fully discrete method to solve for U n+1 , i.e.,
∆t  n 
Uin+1 = Uin − n
Fi+1/2 − Fi−1/2 . (7)
∆x

2.1 REA Algorithm

For this advection equation, the simplest way to compute the numerical flux is the so-called
donor-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 is
entirely determined by its upwind neighboring cell. For example, on the left edge of the i-th cell,
the flux is determined as
n n
Fi−1/2 = aUi−1 , (8)
and similarly on the right edge,
n
Fi+1/2 = aUin . (9)
Using this upwind method, the discrete equation becomes
a∆t
Uin+1 = Uin − Uin − Ui−1
n


, (10)
∆x
from which Uin+1 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 tn , and then perform a simple
piecewise constant reconstruction of the solution in each cell, shown as figure 2a. We then evolve
the solution using the upwind method within the time step, ∆t, i.e., the entire solution profile
moves to the right (as a > 0) for a distance of a∆t, figure 2b. Finally, perform local average in
each cell to obtain the solution U n+1 at tn+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:
1. 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 function
for a time step ∆t.
3. Average the evolved function over each cell to obtain the new cell-averaged value for the
new time step.

3
3 Slope Limiter
In the previous section, in the context of finite volume method/REA algorithm, we discussed
the simplest reconstruction method, namely the piecewise constant reconstruction in each cell. It
is 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 linear
reconstruction 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 (xi , Ui ) in the i-th
cell, marked as ‘x’ in figure 3a. So this really leaves only one freedom we can play with, the slope
of 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 or
counterclockwise (the ‘hand’ symbol) about this pinned center.

"
𝑈!#$

𝑈!"
E

📌
x

"
𝑈!%$

(a) (b)

Figure 3: (a) Illustration of the piecewise linear reconstruction within each cell; (b) Illustration of
how 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 can
be found in textbooks on finite volume method, such as LeVeque [4]. In general, TVD method
can be interpreted as a monotonicity-preserving method, or more intuitively, a method that avoids
non-physical overshoot/undershoot solutions.

𝑈" '$
∆)
∆' 2𝑠( ∆( 𝑠& =
𝑠( = 2∆𝑥
∆( ∆𝑥
∆#
𝑈" 𝑠# =
∆# ∆𝑥
𝑈" #$ 2𝑠#

𝑥 "#$ 𝑥" 𝑥 "'$ 𝑥 "#$ 𝑥" 𝑥 "($ 𝑥 "#$ 𝑥" 𝑥 "+$

(a) (b) (c)

Figure 4: (a) The differences between neighboring cell-averaged values; (b) Slopes relevant to
construct slope limiter; (c) The reference slope;

4
 Before we discuss the slope limiter, several important quantities ought to be introduced and
discussed. First, the differences between neighboring cell-averaged values are defined as:

∆− := Ui − Ui−1
, (11)
∆+ := Ui+1 − Ui

and
∆t = ∆− + ∆+ , (12)
which are illustrated in figure 4a. In above equations, we use := for ‘equal by definition’. It
is noted that their values can be positive or negative. We then introduce the non-dimensional
location indicator, f , defined as
∆−
f := (13)
∆t
which indicates the relative location of Ui in between Ui−1 and Ui+1 . It is noted that, if Ui−1 ,
Ui , and Ui+1 are monotonic, either increasing or decreasing, we have 0 ≤ f ≤ 1. For all other
conditions, Ui is a local extremum, and we have −∞ < f < 0 or 1 < f < +∞. It is also noted
that, traditionally, the indicator is expressed as r := ∆− /∆+ (e.g., Sweby [2]), which makes the
concept of slope limiter less physically intuitive.
The two slopes associated with the ‘-’ (left) and ‘+’ (right) sides are defined as:

Ui − Ui−1 ∆−
s− := =
∆x ∆x , (14)
Ui+1 − Ui ∆+
s+ := =
∆x ∆x
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 2s− , and the orange line has a slope of 2s+ . We will see that, these four
bounding slopes (lines) play the utmost important roles to construct high-resolution
TVD slope limiters.
In figure 4c, a reference slope is also defined, which is the slope of the line that connects the
two neighboring data, i.e., (xi−1 , Ui−1 ) and (xi+1 , Ui+1 ),

Ui+1 − Ui−1 ∆t
sR := = . (15)
xi+1 − xi−1 2∆x

Using sR , we can non-dimensionalize the slopes discussed above. It is easy to see:

s−
φ− := = 2f ,
sR
s+ (16)
φ+ := = 2(1 − f ) ,
sR
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 4f ; and the orange line has a non-dimensional slope of
4(1 − f ).
There are several special cases when two of the four slopes are equal, which are illustrated in
figure 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. We
would like to learn how s should behave when Ui is at different relative locations in between Ui−1

5
 𝑠% = 2𝑠# 𝑠% = 𝑠#

𝑠# = 2𝑠%

𝑥 "#$ 𝑥" 𝑥 "%$ 𝑥 "#$ 𝑥" 𝑥 "%$ 𝑥 "#$ 𝑥" 𝑥 "%$

(a) s+ = 2s− , when f = 1/3 (b) s+ = s− , when f = 1/2 (c) s− = 2s+ , when f = 2/3

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 Ui+1 . Let’s define the non-dimensional slope of the reconstructed piecewise linear function, φ,
as
s
φ := . (17)
sR
It is now equivalent to ask: how φ should behave given different f . We will learn that, to achieve
high-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.

3.1 TVD Rule

The rule for slope limiters to achieve TVD is very simple. The mathematical theory behind the
rules is well explained in Sweby [2], which we do not intend to repeat. Instead, we will illustrate
the two important rules in a graphical way.
TVD rules:

1. If Ui is a local extremum, −∞ < f < 0 and 1 < f < +∞, the piecewise linear reconstruction
must have a zero-slope, i.e., piecewise constant, as illustrated in figure 6a.

2. For all other conditions (0 ≤ f ≤ 1), Ui is in between Ui−1 and Ui+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 Ui is closer to Ui−1 , the reconstructed
linear function must not go beyond the red line to allow its left-edge value to pass Ui−1 . This
is illustrated by the ‘!’ symbol and the ‘lock’-B symbol in figure 6b. Similarly, figure 6c
illustrates the case when Ui is closer to Ui+1 .
In both figures 6b and 6c, the TVD region is simply bounded by the zero-slope blue line and
the red/orange line. Anything between them, marked by the green arrow, is acceptable as a
TVD scheme.

To summarize, for TVD schemes, the non-dimensional slope of the reconstructed piecewise

6
 🔒B ❗

𝑈!"# 𝑈! 𝑈!$# 𝑈!"# 𝑈! 𝑈!$# 𝑈!"# 𝑈!$#

📌
🔒A
🔒B

📌
📌
🔒 ❗ 🔒A

𝑥 !"# 𝑥! 𝑥 !$# 𝑥 !"# 𝑥! 𝑥 !$# 𝑥 !"# 𝑥! 𝑥 !$#

(a) (b) (c)

Figure 6: The limiting (a) Ui is a local extremum; (b) Ui is in between Ui−1 and Ui+1 , closer to
Ui−1 ; (c) Ui is in between Ui−1 and Ui+1 , closer to Ui+1 .

linear function will respect the following conditions:


0 ≤ φ(f ) ≤ 4f ;
 for 0 ≤ f ≤ 0.5
0 ≤ φ(f ) ≤ 4(1 − f ); for 0.5 < f ≤ 1 ; (18)

φ(f ) = 0; otherwise


or in a more compact form:

(
0 ≤ φ(f ) ≤ min (4f, 4(1 − f )) ; for 0≤f ≤1
. (19)
φ(f ) = 0; otherwise

On the φ-f plot, the TVD region is illustrated as the shaded area in figure 8a.

3.2 High-resolution TVD

In the previous section, we have discussed the slope limiter region for being TVD. If we would
like to further improve the slope limiter to achieve higher order spatial algorithm, clearly, it will
only be a subset of the TVD scheme. The mathematical derivation of high-resolution TVD schemes
can also be found in Sweby [2]. Again, we will proceed our discussion on a set of monotonically
increasing data, such that the four slopes, s− , 2s− , s+ , and 2s+ , 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 four
bounding slopes.

Figure 7 illustrates the two smallest slopes for various f values. Especially, figure 7c shows a
special case where the two smallest slopes are identical, and also equal to the reference slope. On
the φ-f plot, the high-resolution TVD region can be easily found in between the two lines having
smallest values among the four. This is illustrated as the shaded area in figure 8b. It is noted that
any high-resolution TVD schemes must pass the point:
 
1
φ = 1, (20)
2

7
 2𝑠% 𝑠%
𝑠%
𝑠%

𝑠#
𝑠# 𝑠' = 𝑠% = 𝑠#
𝑠#
2𝑠#

𝑥 "#$ 𝑥" 𝑥 "%$ 𝑥 "#$ 𝑥" 𝑥 "%$ 𝑥 "#$ 𝑥" 𝑥 "%$ 𝑥 "#$ 𝑥" 𝑥 "%$ 𝑥 "#$ 𝑥" 𝑥 "%$
(a) 0 < f < 1/3 (b) 1/3 < f < 1/2 (c) f = 1/2 (d) 1/2 < f < 2/3 (e) 2/3 < f < 1

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 colored
lines, which have the smallest slopes among the four lines. For condition (c), the reconstructed
linear function must have a slope of sR .

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.

4 4

3 3

4(1-f) 4f 4(1-f) 4f
ɸ(f)

ɸ(f)

2 2
2(1-f) 2f

1 1
TVD
1 2
3 3
0 1/2 1 0 1/2 1
f f

(a) (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 TVD
region.

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/2,
and it seems natural to construct both linear functions with the same slope. This requires the
following condition:
φ(1 − f ) = φ(f ). (21)
Many well-known high-resolution TVD schemes do possess this symmetric feature. Figure 10 shows
several such examples, which also include the sin slope limiter proposed by Berger et al. [3] to

8
demonstrate how a new symmetric slope limiter could be easily constructed. However, we note
that 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

0 ≤ f ≤ 1/3

 4f

2(1 − f ) 1/3 < f ≤ 1/2

φ(f ) = (22)


 2f 1/2 < f ≤ 2/3

4(1 − f ) 2/3 < f ≤ 1

• minmod, the lower bound of the high-resolution TVD region

φ(f ) = min[2f, 2(1 − f )] (23)

• Barth-Jespersen, also known as MC

φ(f ) = min[1, 4f, 4(1 − f )] (24)

• van Leer
φ(f ) = 4f (1 − f ) (25)

• van Albada
2f (1 − f )
φ(f ) = (26)
f 2 + (1 − f )2

• sin
φ(f ) = sin(πf ) (27)

9
 1.5

1.0

ɸ(f)
superbee
0.5 minmod
Barth-Jespersen
van Leer
van Albada
sin
0.0
0 0.25 0.5 0.75 1
f

Figure 10: Illustration of several symmetric high-resolution TVD slope limiters.

4 Non-uniform Mesh
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 theoretical
derivation to define the upper bound of TVD region and the generalized min slope limiter for
high-resolution schemes on non-uniform meshes. They also specified the necessary condition when
dealing with linear data (see f2 in equation (31)).
In the previous section, we have interpreted the high-resolution TVD region on uniform meshes
graphically using the bounding slopes. It seems natural that the same rule can be applied to
irregular 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-resolution
TVD 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:

∆xi−1
a :=
∆xi
(28)
∆xi+1
b :=
∆xi

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 sR is defined as the slope of the line that connects
(xi−1 , Ui−1 ) and (xi+1 , Ui+1 ):

Ui+1 − Ui−1 2∆t

sR := = , (29)
xi+1 − xi−1 (2 + a + b)∆xi

as illustrated in figure 11a. The four bounding lines are shown in figure 11b, and their non-
dimensional slopes are:
s− 2+a+b
φ− := = f
sR 1+a
slef t
φlef t := = (2 + a + b)f
sR
(30)
s+ 2+a+b
φ+ := = (1 − f )
sR 1+b
sright
φright := = (2 + a + b)(1 − f )
sR

in which, f is the same as defined for uniform-sized mesh, and f := ∆− /∆t . Similar to the uniform
mesh case, at several special f values, two of these four slopes are equal. These special f values
are given as follows:
1
f1 = when φlef t = φ+
2+b
1+a
f2 = when φ− = φ+ (31)
2+a+b
1+a
f3 = when φ− = φright
2+a
The direct outcome of Conjecture 1 is that, on non-uniform meshes, the high-resolution TVD
slope is restricted by the two smallest slopes among the four colored lines. On the φ-f plot, this is
the shaded region illustrated in figure 12b. In the same figure, the high-resolution TVD region on
uniform mesh is also shown for comparison. It is easy to see that figure 12b reduces to figure 12a
on a uniform mesh, where a = b = 1.
At this stage, it is tempting to obtain high-resolution TVD scheme on non-uniform meshes by
projecting existing high-resolution TVD schemes on uniform meshes. This can be done by finding
the transform matrix that will homogeneously transform the two triangles of figure 12a into the
two in 12b. The transformed high-resolution TVD schemes can then be obtained by using this

11
transform matrix. This however does not seem to be the most efficient way, and there does not
seem 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 et
al. used a different reference slope, so in our plot, the y-axis is scaled by a factor of (2+a+b)/2; and
2) Berger et al. only marked the TVD region, not the high-resolution TVD region as we illustrate
in figure 12b. Berger et al. [3] also proposed two functions that lay in the TVD region, which can
be seen as generalized van Leer slope limiters, as they recover the usual van Leer slope limiter on
uniform meshes. It can also be easily proved that, with proper scaling of (2 + a + b)/2, they lay in
the high-resolution TVD region proposed in this article. To clarify the discussion, one of the slope
limiter provided by Berger et al. [3], i.e., their equation (38), is scaled and provided as:
 "  1/a #
f 1 − a f
f ≤ f2


2+a+b

 1 + a f 2
φ(f ) = " 1/b # (32)
2 
 b

1 − f
(1 − f ) 1 − f > f2


1+b f2

As discovered by Berger et al. [3], this slope limiter does a very good job on various non-uniform
meshes. The details on the performance of this slope limiter on non-uniform meshes are referred to
Berger et al. [3]. The same slope limiter has also been used by the author in a previous numerical
simulation study of the Welander natural circulation problem, which showed very good results [5].

𝑎 𝑏
1+ +
2 2

2
2 +𝑎 + 𝑏 𝑓 2 + 𝑎 +𝑏 1 −𝑓

4f 4(1-f)
ɸ(f)

1
ɸ(f)

2f 2(1-f)

2 +𝑎 + 𝑏 2 + 𝑎 +𝑏
𝑓 1 −𝑓
1+ 𝑎 1+𝑏

0 1/3 1/2 2/3 1 0 f1 f2 1/2 f3 1

f f

(a) (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 TVD
region.

12
5 Summary
In this paper, we have illustrated how the concept of slope limiter can be interpreted graphically
for 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 for
non-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, Office
of Nuclear Energy, under contract DE-AC02-06CH11357.

References
[1] van Leer, B., 1974. Towards the ultimate conservative difference scheme. II. Monotonicity
and 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 conservation
laws. 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. Progress in Nuclear Energy,
94, pp.162-172.

13