Extremum-Preserving Limiters for MUSCL and PPM

# Extremum-Preserving Limiters for MUSCL and PPM

Limiters are nonlinear hybridization techniques that are used to preserve positivity and monotonicity when numerically solving hyperbolic conservation laws. Unfortunately, the original methods suffer from the truncation-error being order accurate at all extrema despite the accuracy of the higher-order method [1, 2, 3, 4]. To remedy this problem, higher-order extensions were proposed that relied on elaborate analytic and geometric constructions [5, 6, 7, 8]. Since extremum-preserving limiters are applied only at extrema, additional computational cost is negligible. Therefore, extremum-preserving limiters ensure higher-order spatial accuracy while maintaining simplicity. This report presents higher-order limiting for computing van Leer slopes and adjusting parabolic profiles. This limiting preserves monotonicity and accuracy at smooth extrema, maintains stability in the presence of discontinuities and under-resolved gradients, and is based on constraining the interpolated values at extrema (and only at extrema) by using nonlinear combinations of derivatives. The van Leer limiting can be done separately and implemented in MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) [2] or done in concert with the parabolic profile limiting and implemented in PPM (Piecewise Parabolic Method) [9, 10]. The extremum-preserving limiters elegantly fit into any algorithm which uses conventional limiting techniques. Limiters are outlined for scalar advection and nonlinear systems of conservation laws. This report also discusses the order correction to the point-valued, cell-centered initial conditions that is necessary for implementing higher-order limiting. The material herein complements Colella and Sekora [11]. Lastly, there is no guarantee that extremum-preserving limiters preserve positivity. To ensure this property, one should combine the limiting with FCT (Flux-Corrected Transport) [3].

## 1 Algorithms for Scalar Advection

Consider the following scalar equation in one spatial dimension:

 ∂a∂t+λ∂a∂x=0,  λ=constant. (1)

At time-step , the average-valued, cell-centered quantity over a finite volume of length is:

 ani≈1h∫(i+12)h(i−12)ha(x,nΔt)dx. (2)

MUSCL/PPM are conservative finite volume methods that are used to compute :

 an+1i=ani−λΔtΔx(an+12i+12−an+12i−12), (3)

where is the average of a linear/parabolic interpolant over the interval swept out by the characteristics crossing the cell face at and is given by:

 an+12i+12=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩a+i+12=1σh∫(i+12)h(i+12−σ)haIi(x)dxλ≥0a−i+12=1σh∫(i+12+σ)h(i+12)haIi+1(x)dxλ≤0, (4)

where is the CFL number and is the linear/parabolic interpolant, such that .

There are three variations of Godunov-type methods for which limiters can be implemented:

1. van Leer Limiter in MUSCL

2. van Leer Limiter + Parabolic Profile Limiter in PPM

3. Parabolic Profile Limiter in PPM

Each of these algorithm variations are discussed below.

### 1.1 Muscl

1. van Leer limit the differences , giving order results. Apply the corresponding boundary conditions.

2. Use the van Leer limited differences to compute order differences:

 Δ4ai=23((ai+1−14Δai+1)−(ai−1+14Δai−1)). (5)
3. Employ piecewise linear reconstruction by computing spatially extrapolated face-centered values at the low and high (left and right) edges of cells:

 ai,±=ai+12(±1−σ)Δ4ai. (6)

### 1.2 Ppm

1. van Leer limit the differences , giving order results. Apply the corresponding boundary conditions.

2. Employ either or order piecewise parabolic reconstruction:

• Method 1: use the van Leer limited differences and compute spatially extrapolated face-centered values at the low and high (left and right) edges of cells:

 4ai,+ = 12(ai+1+ai)−16(Δai+1−Δai), (7) 4ai,− = 12(ai+ai−1)−16(Δai−Δai−1), (8) 6ai,+ = 4ai,+−130(3(Δai+1−Δai)−(Δai+2−Δai−1)), (9) 6ai,− = 4ai,−−130(3(Δai−Δai−1)−(Δai+1−Δai−2)), (10) α± = a±−ai. (11)
• Method 2: employ either or order piecewise parabolic reconstruction without using the van Leer limited differences in Step 1:

 4ai+12 = 712(ai+1+ai)−112(ai+2+ai−1), (12) 6ai+12 = 3760(ai+1+ai)−860(ai+2+ai−1)+160(ai+3+ai−2), (13) α± = ai±12−ai. (14)

It is important to note that when (centered difference) the two formulations for the piecewise parabolic reconstruction are identical.

3. Limit the parabolic profile .

4. Use the PPM predictor values to reconstruct the parabolic profile:

 a±=ai+α±+σ2(±(α−−α+))−(α−+α+)(3−2σ). (15)

### 1.3 Update Solution for MUSCL or PPM

1. Compute fluxes:

 Fi+12={λa+  λ≥0,λa−  λ≤0.. (16)
2. Use the divergence of the fluxes to update the solution:

 an+1i = ani−ΔtΔx(Fi+12−Fi−12), (17) = ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ani−σ(a+i+12−a+i−12)  λ≥0,ani−σ(a−i+12−a−i−12)  λ≤0.. (18)

### 1.4 Conventional Limiters

To complete the specification of the MUSCL and PPM schemes, one defines the conventional limiters used to constrain the interpolated profiles within each cell.

#### Conventional van Leer Limiter

Given a sequence of average-valued, cell-centered quantities, the conventional van Leer limiter proceeds with the following steps [2, 9, 10]:

1. Compute one-sided and centered differences:

 Δ−ai = ai−ai−1, (19) Δcai = 12(ai+1−ai−1), (20) Δ+ai = ai+1−ai. (21)
2. Apply the conventional van Leer limiter:

 Δlimai = 2min(|Δ−ai|,|Δ+ai|), (22) S = sign(Δcai), (23) Δai = {min(|Δcai|,|Δlimai|)SΔ−ai Δ+ai>00Δ−ai Δ+ai≤0. (24)

One significant defect of this method is the clipping of extremum when . This clipping sets as a precautionary measure for suppressing spurious oscillations.

How one arrives at the formula for the conventional van Leer limiter can be understood by considering the following example. Assume that one is not at an extremum, , and . The value of at face-center is approximated by:

Clearly, and . Therefore:

and one arrives at Eq. 22. However, at extrema , which leads one to Eq. 24.

#### Conventional Parabolic Profile Limiter

Given the higher-order reconstruction of such that , the conventional parabolic profile limiter proceeds with the following steps [9, 10]:

1. Adjust according to the following cases:

 α±→0α+α−≥0,α+→−2α−α2+>4α2−,α−→−2α+α2−>4α2+. (27)
2. Reconstruct given the adjusted values for .

One significant defect of this method is the adjustments to make the reconstruction a monotone profile as a precautionary measure for suppressing spurious oscillations. However, this constraint is more restrictive than is required to preserve monotonicity [13].

The adjustments to are derived from the interpolation polynomial described in the original Piecewise Parabolic Method:

 ai = a−+σ(δ±a+a6(1−σ)), (28) δ±a = a+−a−=α+−α−, (29) a6 = 6(ai−12(a++a−))=−3(α++α−). (30)

where is the CFL number, which also corresponds to a dimensionless length scale. This scale is associated with each grid cell such that the left side of a cell is designated and the right side of a cell is designated . Differentiating with respect to gives:

 daidσ=δ±a+a6(1−2σ)=(α+−α−)−(α++α−)(1−2σ). (31)

By evaluating at the left and right sides of the cell:

 daidσ∣∣−σ→0 = δ±a+a6=−2(α++2α−), (32) daidσ∣∣+σ→1 = δ±a−a6=2(2α++α−). (33)

Maximize with respect to by setting the derivatives equal to zero and solving the resulting equations. One arrives at the following result:

 α+→−2α−σ→0,α−→−2α+σ→1. (34)

## 2 Extremum-Preserving Limiters

For smooth solutions away from extrema, the MUSCL scheme is order accurate for linear advection whereas the PPM is order accurate for linear advection and order accurate in the limit of vanishing CFL number. However, the monotonicity constraints at extrema reduce the truncation error to even at smooth extrema. This reduction in the overall accuracy of the method also introduces a non-smooth component to the error. To eliminate this problem, one constructs a new limiting scheme at extrema.

The defect of the standard approach to limiting is most easily seen in the MUSCL limiter. Away from extrema, the magnitude of the slope is computed as the minimum of three undivided differences: the centered difference and twice the one-sided differences. In smooth regions away from extrema, the centered difference and the one-sided differences all approximate such that the minimum is always defined by the centered difference. At discontinuities, one of the one-sided differences is typically much smaller than the other two differences. Therefore, this one-sided difference is chosen because it leads to a reduction in the slope and suppresses oscillations. However, the idea behind this method fails at extrema because the derivative vanishes. Furthermore, the one-sided differences have opposite signs and this non-constant multiple bounds the centered difference. In the original van Leer and PPM limiters, the solution is to simply drop the order of the method to order.

In the approach used in this report, one changes the limiters at extrema, and only at extrema, by using comparisons of different estimates of the derivatives as a basis for whether to limit the underlying linear scheme. If the solution is smooth at the extremum, then all of the estimates of the derivative are comparable and the limiter leaves the underlying linear scheme unchanged. Discontinuities, underresolved gradients, and high-wavenumber oscillations are detected either by one of the estimates of the derivative being much smaller than the others or by the various estimates of the derivatives changing sign. Either effect triggers a nontrivial limiting of the interpolating function in the cell and a resulting suppression of oscillations.

### 2.1 Extremum-Preserving van Leer Limiter

The extremum-preserving van Leer limiter parallels the conventional van Leer limiter:

1. Compute one-sided and centered differences as well as one-sided differences that are an additional spatial step-size away:

 Δ−−ai = ai−1−ai−2, (35) Δ−ai = ai−ai−1, (36) Δcai = 12(ai+1−ai−1), (37) Δ+ai = ai+1−ai, (38) Δ++ai = ai+2−ai+1. (39)
2. Test for extrema. An extremum is defined if the following condition is satisfied:

 min(Δ− aiΔ+ai, Δ−−ai Δ++ai)<0. (40)
3. If the above extremum condition is not satisfied, then limiting follows the conventional van Leer method:

 Δlimai = 2min(|Δ−ai|,|Δ+ai|), (41) S = sign(Δcai), (42) Δai = min(|Δcai|,|Δlimai|)S. (43)
4. If the above extremum condition is satisfied, then:

• Compute one-sided and centered derivatives:

 D2−ai = 1h2(a(ih)−2a((i−1)h)+a((i−2)h)), (44) D2cai = 1h2(a((i+1)h)−2a(ih)+a((i−1)h)), (45) D2+ai = 1h2(a((i+2)h)−2a((i+1)h)+a(ih)). (46)
• Find the minimum derivative over the five-cells in question:

 S2 = sign(D2cai), (47) D2limai = min(|D2cai|,max(S2D2−ai,0),max(S2D2+ai,0)). (48)
• Apply the modified van Leer limiter at the extremum:

 Δlimai = ⎧⎪⎨⎪⎩min(CVL3h22D2limai,2|Δ−ai|)S2Δcai<0min(CVL3h22D2limai,2|Δ+ai|)else (49) S = sign(Δcai), (50) Δai = min(|Δcai|,|Δlimai|)S. (51)

designates a derivative while designates a difference. Derivatives and differences are similar operators that can be transformed back-and-forth when one considers the relevant Taylor expansion and multiplies/divides each operator by factors of . is a constant that is independent of the mesh spacing and as , the extremum-preserving van Leer limiter reduces to the conventional van Leer limiter. For most calculations, .

How one arrives at the tighter bound of for peak height in the extremum-preserving van Leer limiter can be understood by considering the following example. Assume that one is at a local maximum such that and . The van Leer limiting condition bounds the derivative on the high (right) edge of the cell:

 Δai≥2(ai+1−ai). (52)

When looking for a bound on from the low (left) edge of the cell, an extremum near implies that:

Therefore, one arrives at the following bound on :

These bounds are summarized in the following inequality:

where the comes from the condition that is the nearest face-centered point to the cell-centered point for which one can unequivocally assert that .

### 2.2 Extremum-Preserving Parabolic Profile Limiter

If the piecewise parabolic reconstruction is done without using van Leer limited differences, then one has to include an additional step that limits extrema at cell faces such that lies between adjacent cell averages. Van Leer limiting automatically enforces this constraint.

• Test for extremum at cell faces:

 (ai+12−ai)(ai+1−ai+12)<0. (56)

If no extremum is found, then proceed without any adjustment to .

• If an extremum is found, compute one-sided and centered derivatives:

 D2−ai+12 = 1h2(a((i+1)h)−2a(ih)+a((i−1)h)), (57) D2cai+12 = 3h2(a((i+1)h)−2a((i+12)h)+a(ih)), (58) D2+ai+12 = 1h2(a((i+2)h)−2a((i+1)h)+a(ih)). (59)
• Find the minimum difference with respect to that is greater than zero:

 S2i+12 = sign(D2cai+12), (60) D2limai+12 = max(min(CPPMS2i+12D2−ai+12,S2i+12D2cai+12, (61) CPPMS2i+12D2+ai+12),0).
• Adjust according to:

 ai+12=12(ai+1+ai)−16D2limai+12. (62)

The above formulas are arrived at by considering the centered derivative:

 D2cai+12=1(h/2)2(a((i+1)h)−2a((i+12)h)+a(ih)). (63)

The spatial resolution of enters because one is considering differences between . By substituting in the following expressions:

 ai = a(ih)−h224D2cai+O(h4), (64) ai+1 = a((i+1)h)−h224D2cai+1+O(h4), (65)

which convert point-valued to average-valued quantities, one arrives at:

 h24D2cai+12=ai+1−2a((i+12)h)+ai−h212D2cai+12+O(h4). (66)

Rearranging the above equation gives the expression:

 D2cai+12=3h2(ai+1−2a((i+12)h)+ai). (67)

After finding the minimum derivative , one can again use the above equation to assign a value to :

 ai+12=12(ai+1+ai)−16D2limai+12. (68)

Lastly, and one proceeds with limiting .

Now given the higher-order reconstruction of such that , adjust according to the following cases:

• Compute the derivative with respect to from the the interpolation polynomial described in the original Piecewise Parabolic Method:

 D2PPMai=d2aid2σ=−2a6=6(α++α−). (69)

It is important to note that this derivative with respect to is also a order difference for . Furthermore, this difference represents the maximum difference that can occur across a cell given a parabolic profile.

• Compute one-sided and centered order derivatives:

 D2−ai = 1h2(a(ih)−2a((i−1)h)+a((i−2)h)), (70) D2cai = 1h2(a((i+1)h)−2a(ih)+a((i−1)h)), (71) D2+ai = 1h2(a((i+2)h)−2a((i+1)h)+a(ih)). (72)
• Find the minimum difference that is greater than zero, given the derivatives over the five-cells in question as well as the above difference that was derived from the interpolation polynomial:

 S2PPM = sign(D2PPMai), (73) D2limai = max(min(S2PPMD2PPMai,CPPMS2PPMD2−ai, (74) CPPMS2PPMD2cai,CPPMS2PPMD2+ai),0).
• Adjust according to and :

 α±→α±|D2limai||D2PPMai|. (75)
• Compute the maximum value of over a given cell:

 αmax−=−α2+4(α++α−). (76)

This formula is arrived at by averaging the interpolation polynomial from the original Piecewise Parabolic Method on the left side of cell :

 f− = 1σ∫σ0a(ξ)ξ (77) = a−+σ2(δ±a+(1−2σ3a6)) (78) = a−+σ2(α+−α−)−3σ2(1−2σ3)(α++α−). (79)

Maximize with respect to :

 df−dσ=0⇒σmax=α++2α−2(α++α−). (80)

Therefore:

 αmax−=f−(σmax)−ai=−α2+4(α++α−). (81)
• Preserve monotonicity by ensuring that the following condition is just satisfied:

 I−≤ai−1−ai⇒I−=ai−1−ai=αmax−. (82)
• Solve the above equation for :

 S− = sign(α−), (83) α+ = −2I−−2S−(I2−−α−I−)1/2. (84)
• Compute the maximum value of over a given cell:

 αmax+=−α2−4(α++α−). (85)

This formula is arrived at by averaging the interpolation polynomial from the original Piecewise Parabolic Method on the right side of cell :

 f+ = 1σ∫11−σa(ξ)ξ (86) = a+−σ2(δ±a−(1−2σ3a6)) (87) = a+−σ2(α+−α−)−3σ2(1−2σ3)(α++α−). (88)

Maximize with respect to :

 df+dσ=0⇒σmax=2α++α−2(α++α−). (89)

Therefore:

 αmax+=f+(σmax)−ai=−α2−4(α++α−). (90)
• Preserve monotonicity by ensuring that the following condition is just satisfied:

 I+≥ai+1−ai⇒I+=ai+1−ai=αmax+. (91)
• Solve the above equation for :

 S+ = sign(α+), (92) α− = −2I+−2S+(I2+−α+I+)1/2. (93)
1. Reconstruct given the adjusted values for .

## 3 Numerical Results

Results are presented for 1D scalar advection. The standard test problems were employed to demonstrate improvements in accuracy [14]. The following parameters were used in all test problems:

 σ=0.2,  x∈[0,1],  t=10,λ=1. (94)

Periodic boundary conditions were used with the following number of ghost cells:

 4th Order Reconstruction ⇒ 3 ghost cells per side (95) 6th Order Reconstruction ⇒ 4 ghost cells per side. (96)

The test problems were defined by the following initial conditions:

 Gaussian Wave: a = exp(−256(x−0.5)2) (97) Semi-Circle Wave: a = {(0.252−(x−0.5)2)0.25

The conventional initialization for order accurate methods is obtained by approximating the average over a cell by the value at the center of the cell, i.e., the midpoint rule for integrals. However, that initialization is only order accurate, and less accurate than what one expects for PPM when applied to the advection equation. For that reason, one uses a order accurate approximation to the cell average, following [12]:

 ai=a(ih)+h224D2cai,  D2cai=1h2(a(i+1)h)−2a(ih)+a((i−1)h)). (100)

This order correction is only administered when defining the initial conditions. The corresponding boundary conditions are applied and one begins advancing the temporal loop. For smooth solutions away from extrema, PPM is order accurate for linear advection and order accurate in the limit of vanishing CFL number . The following definitions for the n-norms and convergence rates are used throughout this note. Given the numerical solution at resolution and the analytic solution , the error at a given point is:

 ϵri=ari−ui. (101)

The 1-norm and -norm of the error are:

 L1=∥ϵ∥1=∑i|ϵi|Δx,    L∞=∥ϵ∥∞=maxi|ϵi|. (102)

The convergence rate is measured using Richardson extrapolation:

 Rn=ln(Ln(ϵr)/Ln(ϵr+1))ln(Δxr/Δxr+1). (103)