One-dimensional flux-corrected advection algorithms are based on the finite volume approximation equation (7.121)
in which and are again the flux averaged over the adopted time step at the edge of grid cell (j). The presence of spurious oscillations in the solution is avoided by preventing the total variation (TV) in the discrete representation of the solution
from increasing as the integration proceeds. This is accomplished by limiting the amplitude of the upstream and downstream fluxes, so that the following condition:
(7.137)
is fulfilled (total variation diminishing or TVD condition). An increase in the total variation (TV) is a measure of the formation of oscillations in the solution.
Fluxes can be limited in the discrete form of the advection equation by specifying the fluxes at each edge of the finite volume cells [here for (j – 1/2)] as
(7.138)
where represent the fluxes calculated by a low-order and a high-order method, respectively. The flux limiter functions Φ(rj–1/2) are expressed as a function of parameter rj–1/2 defined as
(7.139)
(7.140)
Similar expressions can be established for the flux at the other edge (j + 1/2) of the cell. One can show that, to fulfill the TVD condition, the flux limiter function must be chosen such
In practical terms, the solution can be made monotonic by adopting for the Φ(r) a value close to zero in the vicinity of sharp gradients (low-order method) and a value close to 1–2 (higher-order method) in regions where the solution is expected to be smooth.
To illustrate the flux limiter method, we consider the linear piecewise scheme discussed in Section 7.4. In this particular case, the flux at the left interface of a grid cell (j) can be expressed by (7.128) and (7.129) or
(7.141)
(7.142)
A flux limiter required to preserve monotonicity is introduced by adjusting the second term in expression (7.142). This is accomplished by replacing bj–1 Δx and bj Δx in equations (7.141–7.142) by Φ(rj–1/2) (Ψj – Ψj–1) The concept of flux limiter is therefore similar to the concept of slope limiter, which is sometimes used to characterize the flux-corrected transport. The resulting “corrected flux” is therefore
(7.143)
(7.144)
with rj–1/2 defined by expressions (7.139) or (7.140), depending on the sign of the wind velocity c. By adjusting all indices in the above expression, one finds the value of the corrected flux at the right edge of cell (j), and the solution at time tn+1 is derived by applying (7.121) with an appropriate choice for the limiter Φ(r). Note that the choice of Φ(r) = 0 and Φ(r) = 1 for all values of r, corresponds the donor cell and Lax–Wendroff algorithms, respectively. Similarly, the choice Φ(r) = r and Φ(r) = (1 + r)/2 corresponds respectively to the Beam–Warming and the Fromm schemes discussed in Section 7.4. None of these four limiters satisfy the TVD condition, and as a result the corresponding schemes do not provide monotonic solutions.
Several formulations for flux limiters that satisfy the TVD condition and hence lead to monotonic solutions have been proposed (Roe, 1986). The following expressions
Φ(r) = max [0, min(1, r)]
define the minmod method,
Φ(r) = max [0, min(2r, 1), min(r, 2)]
the superbee method, and
the van Leer algorithm. Figure 7.18 shows as a function of parameter r the values of the three limiters Φ(r), as well as the domain in which the TVD condition is met. Figure 7.19 shows how the application of a minmod and superbee flux correction improves the solution. Several other limiters have been proposed to enforce monotonicity of the solution.
Figure 7.18 Flux/slope limiter functions (blue curve) for the superbee, minmod and van Leer algorithms superimposed on the regions (shaded) in which the TVD condition is met.
Courtesy of Graham W. Griffiths.
Figure 7.19 Advection of a sharp discontinuity (step function) using six different algorithms: the first-order diffusive donor-cell scheme; the second-order non-monotonic Fromm and Beam–Warming algorithms; and the flux-corrected minmod and superbee methods. The results are obtained after 300 time steps with Δt = 0.1 over a grid of 100 points with a spacing Δx = 1.
From C. P. Dullemond with permission.
An interesting numerical method that overcomes the excessive diffusion of the upstream algorithm and provides monotonic solutions is the flux-corrected scheme developed by van Leer (1977, 1979). To describe this algorithm, we start again from the finite volume expression
(7.145)
where Ψj applies to the center of cell j and Fj+1/2 and Fj–1/2 are the time-averaged flux across boundaries j + ½ and j – ½, respectively. We assume here that the wind field is not uniform, and Fj+1/2 is therefore computed as the product of the velocity uj+1/2 by an estimate of Ψ at the grid cell boundary. This estimate is obtained from a Taylor’s series expansion on the gridded Ψ field:
(7.146)
(7.147)
where Ψ/Δx corresponds to the slope bj of the subgrid function ψ(x). Van Leer proposes for Ψ an expression that limits the flux at the edge of the grid cells (see also Allen et al., 1991):
(7.148)
if (Ψj – Ψj–1) (Ψj+1 – Ψj) > 0, and by otherwise. The algorithm is easily extended to two dimensions by using expression (7.135) rather than (7.145) as the initial step. Although considerably less diffusive than the upstream scheme, this algorithm still contains scale-dependent diffusion. Once the spatial distribution of the transported quantity has diffused to a preferred shape, the diffusion decreases considerably.
7.6 Advanced Eulerian Methods
Several advanced schemes for solving the advection equation have been developed with the purpose of avoiding the spurious oscillations found in high-order methods and the excessive numerical diffusion characteristic of low-order methods. They can be viewed as an extension of some of the fundamental methods discussed in the previous sections.
The MPDATA Scheme of Smolarkiewicz
The Multidimensional Definite Advection Transport Algorithm (MPDATA) proposed by Smolarkiewicz (1983, 1984) focuses on compensating the first-order error of the upstream scheme by reducing the implicit numerical diffusion. Starting from the upstream scheme:
(7.149)
which provides a first guess for the solution at time tn+1, the algorithm uses a second step in which the velocity c is replaced by a compensatory “anti-diffusion velocity” uA defined as
(7.150)
with Kvisc = 0.5(cΔx − c2Δt) = 0.5 c Δx(1 − α). Here, the ratio (1/Ψ) ∂Ψ/∂x is calculated iteratively using the latest estimate of the solution (at the first iterative step). The value of the anti-diffusion velocity at half mesh point j + ½ is therefore
(7.151)
where ε is a small value that ensures that uA is equal to zero when and are equal to zero. The “anti-diffusion” step becomes:
(7.152)
Several iterations can be performed to improve accuracy.
This simple and computationally efficient algorithm is positive definite (if the initial condition is positive) with considerably less implicit diffusion than in the upstream method (Figure 7.20). It is stable under the CFL condition. It does not preserve monotonicity of the transported quantities and, in general, the solutions are not free from small oscillations. The algorithm can easily be extended to multiple dimensions, with an anti-diffusion pseudo velocity defined in each direction. Smolarkiewicz (2006) expanded his MPDATA algorithm to arbitrary finite volume frameworks.
Figure 7.20 Comparison between two advection algorithms applied to a cosine-shaped function (resolved with 12 intervals). (a) first-order upwind scheme; (b) second-order accurate Smolarkiewicz scheme in which the numerical diffusion that characterizes the upwind scheme is compensated by the introduction of an “anti-diffusion” velocity. As in Figure 7.9, the adopted grid is uniform with 500 cells. The Courant number is 0.5 and the solution is shown after 1600 time steps.
From Smolarkiewicz (2006).
&nb
sp; The SHASTA Scheme of Boris and Book
The Sharp and Smooth Transport Algorithm (SHASTA) proposed by Boris and Book (1973) is an Eulerian finite difference algorithm that makes use of the flux-corrected transport (FCT) technique described in Section 7.5. It ensures monotonicity of the solution, conserves mass, and handles steep gradients and shocks particularly well. The scheme includes an advection step followed by a corrective step that reduces the effect of the diffusion produced by the first step. We consider here the 1-D case and assume a variable velocity u(x).
Advection step. The SHASTA algorithm first defines fluid elements formed by connecting linearly adjacent values (Ψj and Ψj+1 in Figure 7.21). Each resulting trapezoidal element is displaced by the distance u Δt. Since the wind velocity u(x) is variable in space, the advection is not a simple translation of the initial element; contraction or dilatation along x can take place. We assume a Courant number less than 0.5, so that the function at grid point j can never be advected further than the grid cell boundaries. After the displacement of the function is completed, the displaced elements are interpolated back onto the original Eulerian grid (see Figure 7.21).
Figure 7.21 Advection of a fluid element. (a) Initial condition. (b) Location and shape of the fluid element after the advection step. During a time step Δt the two boundaries of the fluid element at location j and j + 1 are displaced by a distance uj Δt and uj+1 Δt, respectively. Here, the wind velocities are provided at the intermediate time tn+1/2.. At the end of the advection step, the fluid element is deformed if the velocity u(x) is not uniform. (c) Interpolation of the fluid element onto the grid. The light orange fraction remains in cell j while the darker orange fraction goes into cell j + 1.
Adapted from Boris and Book (1973).
In their algorithm, Boris and Book (1973) prescribe the wind velocities at the grid points xj–1, xj, xj+1, and at the intermediate time level tn+1/2. They deduce at each grid point j a first approximate value for the function Ψ and time level tn+1:
(7.153)
where
(7.154)
and
(7.155)
For a uniform velocity c, the displacement of the trapezoid corresponds to a translation without deformation, and expression (7.153) becomes
(7.156)
with α = c Δt/Δx. This expression includes a two-sided differencing expression that approximates the advection, and a diffusion approximation where the diffusion coefficient is the sum of an independent term (1/8) and a velocity-dependent term (α2/2). This second term is smaller since c Δt/Δx is chosen to be less than 0.5. Without the velocity-independent diffusivity, (7.156) is identical to the Lax–Wendroff algorithm. Following the von Neumann analysis, the amplification coefficient associated with the advective step is
(7.157)
and the corresponding amplification factor is
(7.158)
The value of this factor is smaller than one for all wave harmonics if the Courant number is less than 0.5.
Correction (anti-diffusion) step. Assuming that the diffusivity in (7.156) is only weakly velocity-dependent (α ≪ 0.5), a second step is applied to remove the excessive diffusion produced by the advection step. Thus, we write:
(7.159)
where is the approximation for the transported function derived by the first (advective) step. This can be rewritten as
(7.160)
where
(7.161)
represents the amount of material (“flux”) crossing the boundaries of grid cell j during the time step Δt. The amplification coefficient associated with the anti-diffusion step
(7.162)
is real, so that the anti-diffusion step does not affect the phase properties of the solution.
The overall amplification factor for the two consecutive steps is:
(7.163)
in which the velocity-dependent (or α-dependent) term is generally small. Again, the method is stable when the amplification factor is less than or equal to 1.
The anti-diffusion correction step can introduce spurious extrema and negative values, which can again be avoided by applying an FCT constraint. The monotonicity of the solution is indeed preserved if the anti-diffusion flux f never produces values for Ψ at any grid point j that are larger than the values at the neighboring points. This is achieved if, rather than using (7.161), the value of the flux is replaced by the following FCT condition:
(7.164)
where
Δj + 1/2 = Ψj − 1 − Ψj
(7.165)
The FCT step can be improved by replacing the factor 1/8 in (7.164) with a factor that accounts for the velocity dependence of the diffusivity. Further improvements to the SHASTA method that lead to more accurate solutions have been introduced by Boris and Book (1976).
The Piecewise Parabolic Method
In their piecewise parabolic method (PPM), Colella and Woodward (1984) assume that the subgrid distribution ψ(x) of the tracer concentration inside cell j can be represented by a quadratic function:
ψ(x) = ψj − 1/2 + y(x)[ψj + 1/2 − ψj − 1/2 + dj(1 − y(x))]
(7.166)
where y = (x – xj–1/2)/Δx. Coefficients ψj − 1/2 and ψj + 1/2 are the values of ψ(x) at the boundaries xj–1/2 and xj+1/2 of cell j, and
(7.167)
where is the mean concentration in grid cell j at time tn. For constant spacing Δx, it can be shown through interpolation from in nearby zones that the value of ψ(x) at for example xj+1/2 can usually be expressed as:
(7.168)
When the grid cells are unequally spaced, the expressions for ψj + 1/2 and ψj − 1/2 are more complicated. Colella and Woodward (1984) propose a slightly modified interpolation procedure in the presence of sharp discontinuities (shocks) to ensure that these discontinuities remain sharp during the advection step. The effect of this “steepening” process (Carpenter et al., 1990) is shown in Figure 7.22. Large discontinuities can still produce oscillations in the post-shock flow. In this case, it is advised to introduce some dissipation in the neighborhood of the shock. This can be achieved, for example, by flattening the interpolation profile in the vicinity of the discontinuity, which is equivalent to reducing locally the order of the method. In this case, coefficient ψj + 1/2 can be replaced, for example, by
(7.169)
where fj ∈ [0, 1] is an adjustable factor. Far away from discontinuities or if the shock profile is sufficiently broad, coefficient fj should be set to zero.
Figure 7.22 Advection of a square function using the piecewise parabolic method (PPM) (Carpenter et al., 1990) without steepening (a) and with steepening (b). The spatial domain extends over 40 Δx with cyclic boundary conditions. The Courant number α = 0.5. The numerical solution is shown after 1000 time steps (12.5 revolutions). From Müller (1992).
Copyright © American Meteorological Society, used with permission.
The solution at time tn+1 is obtained by:
(7.170)
where, for example,
(7.171)
Other algorithms for rendering the PPM shape-preserving are presented by Colella and Sekora (2008).
The Crowley–Tremback–Bott Scheme
To improve the accuracy of the first-order upstream method, Crowley (1968), Tremback et al. (1987), and Bott (1989a, 1989b) have also proposed to represent the transported quantity within each grid cell j by a polynomial ψ of order (with assumed to be an even integer number). Thus, at time level n, we write
(7.172)
where y = (x – xj)/Δx is a dimensionless variable such that –½ ≤ y ≤ ½. We assume again that the grid spacing Δx is uniform. Coefficients are determined from the requirement that the value of agree with the value of at grid points (i = j – /2, …, j, …j + /2), and that the area covered by in grid cell j equals Δx. Thus, the coefficients are expressed as a function of the values of Ψj, Ψj+1, … at the ( + 1) neighboring points. Table 7.2 provides the values of the coefficients derived for second-order and fourth-order polynomials, while Trem
back et al. (1987) also considers higher order schemes.
Table 7.2 Coefficients aj,k for the = 2 and = 4 versions of the Bott’s area preserving flux form algorithm (after Bott, 1989a and Chlond, 1994)
= 2 = 4
aj, 0
aj, 1
aj, 2
aj, 3 –
aj, 4 –
As before, the solution for grid cell j at time tn+1 is provided by the finite-volume approximation:
In the Crowley–Tremback–Bott scheme, the fluxes Fj+1/2 and Fj–1/2 at the right and left boundaries of grid cell j are estimated from (7.122) and (7.123) in which ψj (y, t) is replaced by its polynomial approximation (7.172) of order .
In the more general case where the velocity u in the x-direction is spatially variable, (7.123) is replaced by (see Bott, 1989a, 1989b; Chlond, 1994)
Modeling of Atmospheric Chemistry Page 39