A group of new conservative remapping schemes based on nonpolynomial approximations is proposed. The remapping schemes rely on the conservative cascade scheme (CCS), which employs an efficient sequence of 1D remapping operations to solve a multidimensional problem. The present study adapts three new nonpolynomial-based reconstructions of subgrid variation to the CCS: the Piecewise Hyperbolic Method (PHM), the Piecewise Double Hyperbolic Method (PDHM), and the Piecewise Rational Method (PRM) for comparison with the baseline method: the Piecewise Parabolic Method (PPM). Additionally, an adaptive hybrid approximation scheme, PPM-Hybrid (PPM-H), is constructed using monotonic PPM for smooth data and local extrema and using PHM for steep jumps where PPM typically suffers large accuracy degradation because of its original monotonic filter. Smooth and nonsmooth data profiles are transported in 1D, 2D Cartesian, and 2D spherical frameworks under uniform advection, solid-body rotation, and deformational flow. Accuracy is compared via the L1 global error norm. In general, PPM outperformed PHM, but when the majority of the error came from PPM degradation at sharp derivative changes (e.g., the vicinity near sine wave extrema), PHM was more accurate. PRM performed very similarly to PPM for nonsmooth functions, but the order of convergence was worse than PPM for smoother data. PDHM performed the worst of all of the nonpolynomial methods for nearly every test case. PPM-H outperformed PPM and all of the nonpolynomial methods for all test cases in all geometries, offering a robust advantage in the CCS scheme with a negligible increase in computational time.
Conservative remapping is essentially the process of transferring data from one grid to another while conserving the global and local integrals and it has several very promising applications in atmospheric sciences. Perhaps the most immediate application would be conservative interpolation for pre- or postprocessing of numerical model data (see Lauritzen and Nair 2008, and the references therein). Another potential application is in developing conservative baroclinic three-dimensional (3D) models, which exploit Lagrangian coordinates in the vertical based on an “evolve and remap” approach. Since the Lagrangian surfaces are subject to deformation, the data must be remapped onto a reference Eulerian grid at a regular interval of time (Lin 2004; Nair and Tufo 2007; Machenhaur et al. 2008). At the top and bottom boundary line, the remapping scheme needs to be as local as possible requiring only minimal neighboring cell information. Finally, one could apply conservative remapping to semi-Lagrangian (SL) transport since it essentially involves the transferring of mass from one grid to another at each time step. In a backward trajectory mode, data on an upstream departure grid is remapped onto a static arrival grid (Rancic 1995; Nair and Machenhauer 2002; Nair et al. 2002, hereinafter NSS02; Zerroukat et al. 2004). This is the application focused upon in the present study.
The SL transport methods have proved to be an effective means of overcoming the well-known Courant–Friedrichs–Lewy (CFL) limitation and stably increasing the temporal truncation error to the same order of magnitude as the spatial truncation error (Staniforth and Cote 1991). The problem with SL transport schemes in the past is that they typically did not conserve the global and local integrals of the transported quantity. Inherent conservation is important because without it, the global integral drifts requiring a fairly arbitrary postprocessing redistribution of mass globally leading to the artificial global transport of local mass losses and the slight artificial diffusion of gradients (Gates et al. 1971; Priestley 1993; Gravel and Staniforth 1994).
Near the mid-1990s, mass-conserving interpolation schemes began to surface via differentiation of a mass interpolator in order to gain a conservative density estimate (Leslie and Purser 1995; Rancic 1995). In recent years, SL transport schemes have been developed for spherical application using a finite-volume paradigm, integrating over control volumes (or cells) in order to not only conserve mass but also preserve monotonicity (Nair and Machenhauer 2002; NSS02; Zerroukat et al. 2002). The method adopted in this study is called the spherical Conservative Cascade Scheme (CCS) from NSS02, which is based on a finite-volume cascade splitting in which the Cell-Integrated Semi-Lagrangian (CISL; Nair and Machenhauer 2002) scheme is employed in each 1D sweep. This cascade splitting is more efficient than prior Cartesian splittings in that it is dictated by the flow and requires only two sets of sweeps (Nair et al. 1999). In Lauritzen (2007), this CCS method has been shown to be unconditionally stable for transport problems.
The CISL (and thus the CCS) is currently based on the Piecewise Parabolic Method (PPM) of Colella and Woodward (1984, hereinafter CW84) and Carpenter et al. (1990) in which parabolas are fit to the cells based on each cell mean and a fourth-order conservative reconstruction of the left and right interface values. Because the integration of parabolas renders a cubic representation of mass and the cell boundaries are calculated to fourth-order accuracy, PPM is formally fourth-order accurate for a regular mesh and smooth data when no limiting is applied. However, to ensure a monotonic solution, the preprocessing monotonic limiter of CW84 is used in this study. This limiter constructs a piecewise constant representation of the data in the presence of local extrema (only a first-order accurate solution) and suffers arbitrary accuracy degradation at steep jumps as the previously fourth-order accurate interface boundaries are altered to keep the parabola monotonic within the cell.
This study introduces three nonpolynomial approximations to the CCS framework for comparison against the monotonic PPM in Cartesian space and on the sphere. The first of the schemes is the Piecewise Hyperbolic Method (PHM) of Marquina (1994), which was further modified in Serna (2006, hereinafter S06). This scheme is especially promising because of its three-cell stencil (as compared with the four-cell stencil of PPM) and the flexibility in directly controlling the inner-cell variation via the exponent of the power limiter introduced in S06. Additionally, S06 showed mathematically that hyperbolas are by nature less oscillatory than parabolas. The other two methods tested are the Piecewise Double Hyperbolic Method (PDHM) of Artebrant and Schroll (2006, hereinafter AS06) and the Piecewise Rational Method (PRM) of Xiao et al. (2002, hereinafter X02). All three methods are adapted to the semi-Lagrangian context from Godunov-type Eulerian frameworks. Thus, each method has some measure of implicit or explicit limiting such that in the Godunov framework, they are Local Total Variation Bounded (LTVB) within a given cell. The LTVB constraint in the Godunov context turns out to provide bounded monotonicity violation in the semi-Lagrangian context. Additionally, this study constructs an adaptive hybrid approximation scheme called PPM-Hybrid (PPM-H) in which the generally most accurate method, PPM, is used for smooth data and for extrema and PHM is used for steep jumps. This is because PHM suffers no formal accuracy degradation in the presence of steep jumps as PPM does. Only local (or narrow stencil) methods are of interest in the present study. Therefore, such methods as spline approximations (Zerroukat et al. 2006), which require a global cell reconstruction, are not considered here.
In section 2, the conservative cascade scheme will be reviewed in detail including the formulation of the integral continuity equation (section 2a), the 1D CISL remapping procedure (section 2b), the 2D conservative cascade splitting (section 2c), and the extension of the CCS to spherical coordinates (section 2d). Section 3 will present the new subgrid reconstructions in detail including the hybrid scheme, PPM-H, and the respective methods of limiting local total variation. Section 4 presents the numerical experiments in 1D–2D Cartesian and spherical geometries and examine the comparative accuracy of each simulation in terms of global error norms. A summary of the study and conclusions will be given in section 5.
2. Conservative cascade scheme
a. Finite-volume semi-Lagrangian formulation
Consider the 1D continuity equation or mass conservation law where ρ and u represents density and fluid velocity, respectively:
As shown in Laprise and Plante (1995), integrating both sides of (1) over a time-dependent control volume and applying the Leibniz rule to the left-hand side renders strict conservation of the volumetric integral of ρ following fluid motion. This integral conservation is also known as the integral form of the continuity equation, where xi−1/2 and xi+1/2 are the left and right boundaries respectively of cell i:
b. 1D CISL scheme
The CISL scheme follows directly from (2) as mass is the quantity transported from departure cells to arrival cells. This process is schematically illustrated in Fig. 1. Note that throughout this sections x*i always denotes the departure location of xi, and xi±1/2 will always refer to the right and left boundary locations of arrival cell i. First, approximations are fit to the cell means in order to describe the variation within each cell. To insure inherent mass conservation, the foremost constraint on all approximations is that the integral across the cell match the cell’s mean density. To formulate this constraint mathematically, the following must hold true:
where Mi is the mass, Δxi is the grid spacing, ρi is the average density, and ρ̃i is the density approximation function for an arbitrary cell, i.
Next, the cell boundaries are traced upstream from their static locations to find their departure locations. Each arrival cell boundary must have a corresponding departure cell boundary, and the departure cells are bounded by the departure boundaries. In Fig. 1, this is illustrated for cell i where boundaries, xi±1/2 are traced upstream to xi±1/2*. In this study, to isolate the errors due to the CISL scheme and subgrid approximations only, these trajectories are calculated analytically.
After the departure cells are defined, the mass within each departure cell is calculated by integrating between each departure cell’s boundaries over the approximations formed earlier. Generally, the departure cells span more than one arrival cell. Thus, in order to compute the mass within each departure cell efficiently, mass is first accumulated up to each arrival boundary (Nair and Machenhauer 2002). Then, the integrated mass from the leftmost arrival boundary up to the right departure boundary is added to the accumulation up to the leftmost arrival boundary (mathematically expressed in NSS02). Finally, the departure cell mass simply replaces the arrival cell value, and the arrival cell grid spacing is divided out to obtain the new arrival cell density in a conservative manner. Mathematically, the process could be expressed as follows, where n represents the current time index:
c. 2D CCS scheme
The CISL can be extended into two dimensions by use of a finite-volume cascade-splitting technique (i.e., the CCS). In the CCS, only two sweeps of remapping are performed: one from the arrival grid cells to the intermediate grid cells and another from the intermediate grid cells to the departure grid cells. Refer to Fig. 1 of NSS02 for a schematic of this process using a rectangular Cartesian grid in (λ, μ) coordinates.
The departure grid is strictly rectangular being bounded by arrival longitudes (straight lines bounding a column of arrival cells to the east and west) and by arrival latitudes (straight lines bounding a row of arrival cells to the north and south). The intermediate grid remains bounded to the east and west by the arrival longitudes, but the north and south intermediate boundaries move with the fluid for each column of arrival cells. Therefore, for the first CCS sweep, mass is remapped northward and southward parallel to arrival longitudes for each column of cells via the CISL method. Before the CISL remapping is performed, however, each column (i) of scalars is multiplied by the zonal grid spacing, ΔλI, so that when CISL integration in the μ direction is performed, total mass is calculated for the intermediate cells. Mass in each intermediate cell is then divided by the intermediate cell area, Aij = ΔλiΔμ*j, to obtain the intermediate cell densities.
Each straight row of cells on the arrival grid corresponds to a curved row of cells on the intermediate grid, which are bounded by intermediate latitudes (curved lines that bound intermediate rows to the north and south). The departure grid cells are bounded to the north and south by intermediate latitudes and to the east and west by tracing the arrival east–west boundaries along intermediate rows following fluid motion upstream. Therefore, the second CCS sweep oriented parallel to intermediate latitudes remaps the intermediate grid cell masses to the departure grid cells to obtain the mass within each departure cell, again using the CISL method. As in the first sweep, before the CISL remapping is performed, each row of scalars is multiplied by the intermediate meridional grid spacing, Δμ*j, so that when CISL integration along intermediate rows is performed, total mass is calculated for the departure cells. Mass in each departure cell then replaces arrival cell mass to perform transport, and the mass is divided by the arrival cell area, Aij = ΔλiΔμj, to conservatively obtain the new density. See NSS02 for additional details on the CCS.
A positive definite filter is employed in 2D simulations in which one sweep is made across the grid to add mass to cells with negative values. Then, the amount of mass added is taken away from the rest of the grid weighted linearly by the mass already in the grid cell.
d. Extension to spherical geometry
To extend the CCS into spherical geometry, an area-preserving transformation is first applied in which the spherical coordinates (λ, θ) are changed into (λ, μ), where μ = sinθ (Nair and Machenhauer 2002). Where θ was previously in the domain θ ∈ [−π/2, π/2], μ is now in the domain μ ∈ [−1, 1], and (λ, μ) represents a rectangular coordinate set with a meridional stretch which can be manipulated in an ordinary Cartesian manner as before.
Also, because of the deformed mesh of the (λ, μ) rectangular grid, a linear intersection calculation to obtain the intermediate grid is not accurate enough. Intersections of intermediate latitudes and arrival longitudes must be calculated with a cubic intersection calculate, and in this implementation, a standard cubic Lagrange interpolant is used. Additionally, because of the excessive grid deformation near the Poles, the constant approximation of the east–west boundaries of departure cells causes undesirable accuracy loss in the rows just outside the polar caps. For this reason, the cells are refined such that the departure grid east–west boundaries may be multiple constant values within a computational cell to allow for a more accurate description of the zonal remapping near the Poles from the intermediate grid to the departure grid. In this study, the intermediate row adjacent to each Pole is refined to give three constant zonal boundaries within the cell, and the intermediate row two cells equatorward is refined to give two constant zonal boundaries.
The Poles of a spherical coordinate set are a well-known problem because the meridians converge to a singular point. For this reason, any row of intermediate cells that contains the arrival Pole point cannot be integrated as the other cells are integrated because the rectangular Cartesian integration procedure breaks down at this singularity. In the present implementation of spherical transport, the meridional Courant number is restricted to less than unity such that the departure polar cap always includes the arrival polar singularity; that is, the departure mesh’s polar point never advects out of the arrival grid’s polar cap. This is not a formal restriction of the CCS method, however (see Nair 2004). The extension to large Courant numbers does not provide additional insight into the comparative performance of these subgrid reconstructions and is therefore not implemented in the present study.
Therefore, in the polar cap region for the second (zonal) CCS sweep, a traditional pointwise cubic SL method is employed since it does not suffer the singular restriction and the polar cap mass is redistributed based on weights formed by the traditionally interpolated SL values. The mass is known in the polar caps from the first sweep, which calculates mass in the meridional direction. Integration is not necessary over the singularity in the first sweep because it can be calculated via subtraction using the accumulated mass in the arrival cells. Suppose each point in a polar cap is given an interpolated SL departure value of ψi. The weights wi are formed by
and the mass is redistributed zonally across the polar cap such that a cell in the polar cap has a mass m, of
where MPC denotes the total mass in the polar cap calculated in the first sweep.
3. Subgrid approximations
Because of the extensive descriptions of PPM in the literature and its common use in modern applications, this method will only be described briefly here. PPM was originally introduced in CW84 for gas-dynamical simulations in an Eulerian context based on fluxes through cell interfaces. It was later applied to meteorological modeling by Carpenter et al. (1990). PPM fits parabolas to each cell to describe the inner-cell variation using a cell’s mean density (ρi) and the left and right interface values (ρi−1/2 and ρi+1/2, respectively). Interface values are reconstructed before parabolas are fit using a conservative, monotonic, cubic interpolation for fourth-order accurate values, which requires four cells of information rendering a PPM requirement of a four-cell stencil. The parabolas are constructed on a normalized domain, ξ ∈ [−1/2, 1/2] with the following form:
where ρS is defined as
To limit PPM such that all parabolas are monotonic within the local cell domain, Carpenter et al. (1990) gives a detailed explanation of how the interface values are altered. It should be noted that much of the cause of PPM error comes from the damping effects of this filter. For smooth functions with extrema (e.g., the sine wave), nearly all error can be attributed to the effects of this filter. Less-damping filters exist at the cost of wider stencils (Zerroukat et al. 2005), but the focus of this study is on local stencil methods.
PHM was originally introduced in a Godunov finite-volume context in Marquina (1994) where hyperbolas are fit to the cells to describe the subgrid variation. The hyperbola is constructed based on second-order approximation to the interface derivatives rather than using interface values. The hyperbola takes on the following form:
where hi is the grid spacing and x0 is the cell center location for cell i, and αi and di are the derived parameters used to define a unique hyperbola within a cell. The value di is the approximation to the derivative at the cell center, and αi is tweaked in order to interpolate one of the lateral derivatives. The limiting of PHM is performed in the evaluation of di, and a recent advance in PHM limiting has been given in S06, in which a power limiter is used to calculate di rather than the Marquina (1994) harmonic limiter. The S06 definition of di is
where di−1/2 and di+1/2 are centered second-order approximations to the interface derivatives given by
the power limiter is expressed as
and the minsign function is given by
There exist two situations in which the constructed hyperbola would become ill-defined across the interval. First, if |αi| ≥ 2, the natural logarithm term would no longer be defined. Therefore, the definition of αi must restrict it to the interval, αi ∈ (−2, 2) which, according to Marquina (1994) and S06, also serves as a sufficient condition for PHM to be Local Total Variation Bounded (LTVB) within a given cell. Second, if αi = 0, a divide-by-zero singularity occurs. In the Godunov context in which only the interface hyperbola derivatives are of interest, this is a removable singularity, which can be computed. However, in the CISL context, a continuously integrable function is necessary. It turns out that αi = 0 simply signifies that di−1/2 and di+1/2 are equal, which in turn means that the function may be approximated accurately by a line. Therefore, if |αi| is close to zero, the hyperbola is replaced by a line of the following form:
The parameter αi is defined based on the value di such that
where the ratio ηi is defined as
where ɛ is a small value to avoid a floating-point divide-by-zero error and may be set to a value as small as machine precision (see X02). S06 formulates the ratio ηi such that no division is necessary based on a fixed power limiter exponent value of p = 3, but this study adopts (17) such that the exponent may be varied within the set of real numbers rather than be restricted to the set of integers. The S06 improvement was attributed to allowing more variation within the cells by increasing the αi range from αi ∈ [−2(2 − 1), −2(2 − 1)] formed by harmonic limiting to αi ∈ [−2(3 − 1), 2(3 − 1)] formed by powerp limiting with p = 3. It was stated in S06 that maxa,bηi(a, b) = 3 when p = 3. Because ηi increases as the difference in magnitude of a and b increases, Fig. 2a shows numerically that for any real power limiter exponent p,
Since the LTVB constraint restricting α to the interval, α ∈ (−2, 2), implies that ηi ∈ (−4, 4), it is clear that a LTVB hyperbola may be obtained with a power limiter exponent in the range p ∈ (0, 4), and thus, this range would also render the CCS solution essentially monotonic with a bounded magnitude of monotonicity violation. Therefore, p = 3.999 is used in this study for PHM-based simulations.
AS06 developed a logarithmic approximation to computational cell variation in the Godunov context called the Piecewise Logarithmic Method (PLM), but this method is not trivially transferred to the CCS context because the logarithm functions are not naturally integrable across the cell. Later, AS06 came out with the third-order Piecewise Double Logarithmic Method (PDLM) with the stated advantage of constructing extrema at full accuracy. This method was implemented in the CCS context because it could be limited and was naturally integrable across a given cell domain. AS06 also developed a Piecewise Double Hyperbolic Method (PDHM) with very similar parameters as the PDLM. Cancellation errors caused accuracy degradation for PDLM on fine meshes; and for this reason, PDHM is the method shown in this paper.
It should be stated that PDHM has very little in common with PHM except that they are both in rectangular hyperbolic form. AS06 defines a double rectangular hyperbola within an arbitrary cell i:
where h is the cell grid spacing; xi is the center location of cell i; and a, b, c, and d̂ are parameters defined by the cell interface derivatives as follows:
where ɛ = 0.1hq and h is the average grid spacing. The parameter d̂i should not be confused with the interface derivatives di±1/2. The limiting of PDHM is performed via the parameter q, which as suggested by AS06 is set to unity. There exists a possibility of singularity if a is not confined to the interval a ∈ (0, 1). Since the ɛ limiting of (20) confines the parameter a to the interval a ∈ [ɛ(1 − ɛ), 1 − ɛ2) (shown in AS06), this case of singularity never occurs as ɛ > 0. Thus, the function is integrable across an arbitrary cell domain.
The double hyperbolic function in (19) does not inherently conserve mass. Therefore, to obtain a unique conservative double hyperbola, the mean integrated value of (19) is subtracted off to render an integrated mean of zero, and the cell mean density is added on as follows:
This method is developed in X02 wherein unique rational functions (ratio of two parabolas) are fit to cells based on the cell mean and the interface values. Unlike with PPM, and unlike the typical Godunov schemes, the PRM was developed in a context in which the interface values are advected in a SL manner, and the cell means are transported via flux form transport. The rational functions are defined as
where ai, bi, and βi, respectively, are given by
As seen before with PHM, ɛ is a small number used to avoid a floating-point divide by zero and may be set as small as machine precision. No singularities exist in the actual function itself as shown in X02 because the definition of βi does not allow the expression 1 + βi (x − xi−1/2) to reach the value zero within a given cell. One notable feature of this method is that in the CISL context, the integration from the closest arrival boundary, xî−1/2, to a given departure boundary, x*i−1/2, up to that departure boundary reduces to a relatively simple expression:
where δx = x*i−1/2 − xî−1/2. Therefore, construction and integration using PRM is very efficient. However, because this method requires values at the cell interfaces, the conservative cubic approximations of PPM are used, which requires additional computation.
In light of the degradation of PPM accuracy in the presence of steep jumps and local extrema in the transported data, a hybrid method was devised in which PPM is used for smooth data and local extrema and PHM is used for steep jumps. It was experimentally found that no replacement scheme yielded consistent gain in accuracy for local extrema, and PHM provided the best solution when used at steep jumps. The reasoning for this adaptive replacement is that PHM retains full order of convergence at steep jumps while remaining essentially monotonic where PPM suffers arbitrary degradation of accuracy at steep jumps due to altered interface values. The hybrid scheme is automatically adaptive because PPM already tests for cases in which parabolas are nonmonotonic within the cell domain due to steep jumps in the preprocessing limiter in CW84.
It was observed experimentally that PPM-H violations of monotonicity were unacceptably large: O(10−2) for a C0 discontinuous jump of magnitude unity. Therefore, two methods were employed to reduce the violation. First, the larger the value for the PHM power limiter exponent, p, the more the variation allowed within a cell. Therefore, a measure of jump severity is employed based up on the relative change in magnitude of the left and right interface derivatives. The larger the difference in interface derivatives, the greater the overshoot or undershoot generally. Therefore, the severity, S, can be calculated as a normalized quantity such that S ∈ [0, 1]:
Therefore, if S is greater than some cutoff value, S*, a PHM power limiter exponent of p = 3 is used and below this cutoff value, p = 3.999 is used. This was observed to successfully reduce overshoots and undershoots by restricting the PHM variation for severe jumps and retain optimal accuracy with a cutoff value of S* = 0.8.
Second, a natural overshoot damping mechanism was observed and verified in PHM. Therefore, PHM is adaptively used for local extrema only when those extrema are caused by overshoots. Overshoot-induced extrema are detected utilizing the same parameter, S, because these particular extrema naturally have one derivative much larger in magnitude than the other. A cutoff value of S* = 0.95 is found to effectively reduce monotonicity violation without altering accuracy to a large extent. Therefore, only when S > 0.95, is PHM used at local extrema.
4. Numerical experiments
In each of these experiments, a certain initial spatial distribution of a scalar quantity (it will be referred to as density) will be transported by a time constant flow (though not generally uniform in space) and compared to the known analytical result for evaluation. It should be noted that performance of the nonpolynomial methods in this CCS context does not necessarily imply similar comparative performance in another transport context such as the Godunov-type transport.
a. 1D uniform advection
1) Sine wave
This test case involves transporting a sine wave uniformly and cyclically through a domain of 100 cells Ω ∈ [0, 1] for 10 revolutions (2000 time steps) with a wind of u = 1 m s−1 time step of Δt = 0.005 yielding a Courant number of C = 0.5. The initial data are given by
such that the profile remains above zero. Figure 2b shows a plot of the new methods versus PPM for the sine wave test with the displayed domain constrained to the sine wave maximum to better visualize the damping errors that cause the bulk of global error. Table 1 gives the errors norms of the methods for quantitative comparison. These L1, L2, L∞, Lmin, and Lmax error measures are suggested and defined by Williamson et al. (1992) and are used throughout the present study for a comparison of methods. Following the findings of Willmott and Matsuura (2005), the L1 error is given preference over the L2 error as most of the comments about root-mean-squared error also apply to the L2 norm. Basically, to square the errors themselves as in L2 gives greater weight to the larger errors, which is not as straightforward as the L1, which exhibits a linear relationship with the true mean error.
Visually analyzing Fig. 2b, it can be immediately seen that two methods perform the best in this test case: PHM and PPM-H. There is a well-known phenomenon with PPM strictly due to its original filter in which it clips the extrema of a sine wave causing an artificial flattening that extends typically a couple of grid points out from the extremum center. Note that without filtering, PPM performs much better than both PHM and PPM-H for the sine wave, but oscillates wildly for highly discontinuous profiles. The same phenomenon is seen with PPM-H, but the steepness is more accurately described outside the flattening. Likewise, two methods clearly perform much worse than PPM: namely PRM and PDHM. PRM shares the PPM extrema clipping problem that is attributable to the fact that in this implementation, PRM shares the monotonic cell interface values calculated for PPM. Clearly, PDHM does not exhibit such a clipping phenomenon, but rather simply diffuses out the sine wave excessively leading to a poor preservation of the original data profile. Table 1 reveals numerically what is easy to realize graphically, and PDHM suffers the worst accuracy of all the methods. Impressively, PPM-H performs about 3 times better than PPM, and PHM has about 38% of PPM L1 error showing a substantial improvement in accuracy for this smooth function. PRM and PDHM have 150%–170% the L1 error of PPM. It should be noted that PHM, as made manifest by the L∞ norm, gives the best approximation to the sine wave peak amplitude.
Table 2 gives L1 error norms and convergence information for the sine wave problem. It is clear that the PHM advantage over PPM due to the PPM monotonic filter increases with increasing resolution. On the whole, PPM-H shows the highest numerical convergence. PRM shows the lowest convergence rate giving rise to relatively poor accuracy at high resolutions.
2) Rectangle wave
In this test case, the data profile consists of a rectangular profile with two C0 discontinuous jumps between zero and unity. All transport parameters are the same as the 1D sine wave experiment, but the data is transported for only one cyclic rotation (200 time steps). Figure 2c shows a plot of the different solutions to this experiment on a truncated domain. Table 1 also gives the error norms for this experiment. This test case represents the worst case scenario for any oscillatory scheme and would immediately reveal if a scheme is nonmonotonic in implementation. Clearly, PPM-H gives the best steepness, and all of the strictly nonpolynomial methods render a worse steepness than PPM. A pattern is shown here that tends to hold throughout the study, and that is that PRM tends to have similar accuracy to PPM in discontinuous profiles but is almost always slightly worse. Again, PDHM diffuses the discontinuity the most and renders the worst accuracy of this study’s methods. The PPM-H improvement in the L1 error norm is about 8% relative to PPM error. Also, PHM is about 40% less accurate than PPM.
It may seem contradictory that PHM replaces PPM for steep jumps in the hybrid method yet PHM performs worse than PPM in the rectangle wave test case. In any total variation bounded scheme, steep slopes will inevitably relax to a shallower slope. When this happens with PHM, the center segment of the slope is steep but the derivatives are not changing sharply. Thus, PHM is only advantageous at the edges of the slope where the derivatives are varying more, and PPM performs better along the center segment of the slope. This notion is confirmed in the convergence analysis, which shows that PHM performance degrades relative to PPM with increasing resolution (meaning more cells in the center segment of the relaxing discontinuity).
Error norms for the rectangle wave test case and convergence information is given in Table 2. The relative ratios of error tend not to change much for this particular test case. One exception would be that PHM converges relatively poorly with refinement in comparison with PPM, which shows poorer PHM accuracy at higher resolutions. The low convergences in general are expected because the data themselves are not even C0 continuous (i.e., not differentiable).
3) Irregular signal
This test case involves the transport of an irregular signal (Zerroukat et al. 2006) with the exact same transport parameters as given in the 1D rectangle wave experiment. Figure 2d shows the numerical solutions to this experiment and Table 3 also gives the various error norms. This test case was designed to be difficult to transport because of the multiple changing inflections and steep jump near x = 0.1. What is most impressive in this figure visually is how much more PPM-H reacts to the sharply changing inflections than the other methods. It is clear that, again, PRM performs similarly but less accurately than PPM and that PDHM diffuses the data the most once again. The PPM-H accuracy improvement in preserving the alternating inflections in the data do not show up much in the overall error norm, which is dominated by the steep discontinuities near x = 0.1, x = 0.3, and x = 0.55. In the L1 norm, PPM-H constitutes an 11% gain in accuracy over PPM for this experiment.
Because computational optimization is fairly trivial in 1D and the only difference in run time is the construction and integration of the subgrid approximations, Table 4 gives a summary of the average and standard deviation of run times for the different methods on an Intel Core2 Duo 6400 2.12-GHz processor with simulation parameters provided in the caption. Generally, it is found that the there is little separation in run time for a realistic computing environment between the methods. One thing to note is that for a regular mesh, methods that use interface values for reconstruction of cell variation may calculate those values with fourth-order accuracy much more efficiently (see CW84). In a realistic spherical context, this occurs only in the zonal sweep accounting for roughly half of a given time step.
b. 2D Cartesian tests
In the 2D experiments (Cartesian and spherical), focus is given to the best candidate schemes only which determined by the error norms would be PHM and PPM-H. Clearly, in 1D, PPM-H gave robust improvements in all cases over the PPM standard. Also, PHM though it can be as much as 2 times less accurate than PPM can also be as much as 3 times more accurate for very smooth functions. Therefore, PPM, PHM, and PPM-H are the only methods shown and analyzed in the following sections of this paper.
1) 1) Solid-body rotation of a cosine hill
Solid-body rotation in Cartesian geometry consists of rotating the domain with uniform angular velocity about the domain center point. This experiment uses a domain of 33 × 33 cells Ω ∈ [0, 32 × 105 m]2 for one revolution (71 time steps) with an angular velocity of ω = 10−5 s−1, and a time step of Δt = 8849.56 s yielding a maximum Courant number of C = 2 within the domain. The definition of the cosine function is given in Zerroukat et al. (2002). Figure 3 shows surface plots of the analytical and numerical solutions and the error norms are given in Table 5. For external comparison, the error measures of the Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme of Zerroukat et al. (2002) and the monotonic SLICE (SLICE-M) scheme of Zerroukat et al. (2005) are included.
This test case could only be thought of as a smooth test case on a finer mesh; however in this context, it resembles a cone more than a cosine hill making it fairly nonsmooth. Clearly, the contact discontinuity at the top is diffused by all of the methods, but PPM-H keeps the profile the most accurately. PHM also clearly diffuses more than PPM leading to the expected result that the methods tested in 1D perform similarly on a 2D regular mesh with a cascade splitting. The error norms, of course, reflect this in which PPM-H gives a 31% increase in L1 accuracy. The L1 error shows that PPM-H performs better than SLICE and SLICE-M in an average sense, but the larger errors that dominate L2 and L∞ (occurring at the cosine hill center) are much better for SLICE and SLICE-M. Clearly, the piecewise constant limiting for the maximum is what causes the excessive diffusion in PPM-H and PPM, which is why SLICE is performing better. SLICE-M on the other hand uses a six-cell stencil postprocessing monotonic filter to enforce monotonicity while keeping extrema well resolved. The performance of this filter keeps it from experiencing the accuracy degradation of PPM-based methods. However, even with the deficiency at the local maximum relative to SLICE-M, PPM-H still gives overall better accuracy with only a four-cell stencil.
2) Solid-body rotation of a slotted cylinder
Here, a Cartesian solid-body rotation is used to transport a slotted cylinder (Zalesak 1979). This experiment uses a domain of 101 × 101 cells Ω ∈ [0, 100 m]2 for one revolution (96 time steps) with an angular velocity of ω = 3.635 × 10−5 s−1, and a time step of Δt = 1800 s yielding a maximum Courant number C = 3.27 within the domain. The specific parameters for the initial distribution are given in (Zerroukat et al. 2002). Surface plots of the results are shown in Fig. 4, and the errors are given in Table 5. Clearly, the steepness of the slot and sides are best preserved by PPM-H and most poorly preserved by PHM. The thin strip along the back of the slotted cylinder is especially sensitive to this diffusion and is clearly preserved well by PPM-H. The error norms show a more modest increase in PPM-H L1 accuracy, about a 9% improvement relative to PPM. PHM, on the other hand, is about 35% worse than PPM in L1 error. Comparing against SLICE and SLICE-M, it is clear that the cubics must have some difficulty when steep jumps are present in the data. It is clear that PPM-H gives a better representation of the slopes of the steep jumps (manifested in the L∞ error) than all methods except SLICE, which is not monotonically filtered.
c. 2D spherical tests
1) Polar solid-body rotation of a cosine hill
Solid-body rotation over the sphere induces a degree of freedom regarding the axis about which the flow rotates. In this test case, the flow is poleward such that it rotates about an equatorial axis. Given the α parameter described in (Nair and Machenhauer 2002), this experiment is run with α = π/2. The initial distribution is a cosine hill formation in which the radius is formed using arc length from the cosine hill center. (Nair and Machenhauer 2002; Nair et al. (2002, hereafter NSS02) also describe this initial distribution mathematically along with the wind components induced from constant angular rotation about an equatorial axis. As before, the trajectories are calculated analytically before the simulation is run to obtain errors due to the spatial scheme and CCS only. The grid is composed of 128 zonal and 65 meridional grid points with a grid spacing of about 2.8° in both directions. One rotation is performed over 12 days with 256 time steps and a maximum meridional Courant number of 0.5.
Figure 5 displays the orthographic projection of the analytical and numerical solutions, and Table 6 gives the error norms for this experiment. Clearly, PHM performs more poorly than the other schemes in preserving the shape and magnitude of the original hill. This flow is not deformational and thus should preserve the shape perfectly. PPM-H did not elongate the cosine hill as severely as PPM, and the L∞ norm shows that PPM-H kept the center magnitude 24% better than PPM. The overall PPM-H solution was 34% better than PPM, and the PHM solution was 140% worse than PPM in the L1 norm. The PHM degradation may be due to the irregular grid spacing in the meridional direction on the (λ, μ) grid. The left and right interface derivatives are approximated with centered measures, but the interface location is not geometrically centered between the cell means.
Additionally, to trace the errors as the flow proceeds over the Pole of the sphere, Fig. 6 gives the L1 and L∞ error norm over time for PPM, PHM, and PPM-H. There L∞ oscillation at the polar crossing is due to the grid singularity. It is evident that there is definitely some irregularity at the Poles as manifested by the L∞ norm at days 3 and 9 (corresponding to the days when the cosine hill center crosses the North and South Poles, respectively). However, this is expected given the handling of the Poles, and the error does not seem to jump too steeply for any of the given methods.
2) Smooth deformational flow
This experiment introduces a flow that deforms the original data into two opposing vortices centered at the Poles of a rotated sphere. In this experiment, the Pole of the rotated sphere is located at λp = π + 0.025 and θp = π /2.2 such that the vortices are formed about 81°N–S. (Nair and Machenhauer 2002; NSS02) give the initial data and flow specifications for this experiment. The grid is the same as with experiment 3-A, and the maximum Courant numbers in the meridional and zonal directions are Cθ = 0.82 and Cλ = 19.1, respectively.
Figure 7 shows the orthographic projection of the analytical solution and absolute errors of this experiment, and Table 6 gives the error norms. It is interesting that the spatial pattern of PHM error differs substantially from the other two methods with a more uniform distribution. Also, it is clear that PPM-H has a very similar error distribution as PPM, but the error magnitude is smaller. PPM-H gives about a 7% increase in L1 accuracy relative to PPM, and PHM performs about 3% worse in the same error norm. Not much changes for the infinity norm, and the minimum and maximum norms reveal that there is a slight violation of monotonicity for this flow for all methods. Part of this (namely, that revealed in the PPM norms) is attributed to numerical errors the cascading approach associated with the deformational character of the flow, and is known to occur (see NSS02). Additionally, since PHM is not strictly monotonic, it undershoots by about 7 times more than PPM yet still only about 0.01% of the interval in the data. PPM-H shows the same overshoot character as does PPM with identical Lmin and Lmax norms demonstrating no addition to the errors introduced in the cascading.
5. Summary and conclusions
Four nonpolynomial methods and one hybrid method based on adaptive use of PPM and PHM have been tested in the Cell-Integrated Semi-Lagrangian context to perform conservative semi-Lagrangian transport. For extension into two dimensions, the Conservative Cascade Scheme of NSS02 was employed. Further extension into spherical geometry was achieved by converting the spherical coordinate system (λ, θ) into a stretched rectangular coordinate system (λ, μ) such that μ = sinθ. Also, special treatment was used for the polar caps such that a typical pointwise semi-Lagrangian interpolating method was used to generate weights for redistributing the mass.
The new nonpolynomial methods did not give any robust or substantial improvement to the existing baseline of comparison, PPM, in terms of accuracy alone in the general case. It is important to note, however, that PHM involves only a three-cell stencil as compared to the four-cell stencil of PPM, which would render PHM more efficient in a parallel architecture. Additionally, near a rigid boundary, a local stencil is advantageous because a higher-order approximation can be used closer to the boundary. PHM only improved upon PPM accuracy when the majority of PPM error was caused by degradation near sharp derivative changes due to the PPM monotonic filter. PPM-H was found to robustly improve upon PPM’s performance in the 1D, 2D Cartesian, and 2D spherical geometries for both smooth and nonsmooth data. Particularly, the nondeformational transport of smooth data proved to show the best PPM-H improvement over PPM. The PPM-H improvement comes at no measurable increase in computational cost rendering superior PPM-H efficiency.
This conservative remapping process has other applications such as conservative interpolation and remapping of vertical coordinates to a reference grid in which the relative performance of these new methods may differ. For this reason, in the future it would be important to investigate how these subgrid approximations perform in the context of grid-to-grid remapping. PHM stands to improve upon PPM especially in this application since its three-cell stencil would allow full-order interpolation one cell closer to a rigid boundary as in the remapping of Lagrangian vertical coordinates to a reference grid near the land surface. Also, it would be important to perform a comparative study of PPM and the new nonpolynomial and hybrid methods in an Eulerian Godunov framework including simple transport and the shallow-water equations since the SL is not appropriate for all applications, especially considering the growing need for massively parallel model architectures.
The first author extends much gratitude to the Summer Internships in Parallel Computational Science (SIParCS) headed by Dr. Rich Loft in association with the Institute for Mathematics Applied to Geosciences (IMAGe) and the National Center for Atmospheric Research (NCAR) for funding this research. The authors are also thankful to Dr. Peter Lauritzen for revising the manuscript. Many thanks are extended to Dr. Fredrick Semazzi of the North Carolina State University for coordinating the opportunity for this research to occur. The authors also extend their appreciation to the reviewers for their insightful comments.
Corresponding author address: Matthew Norman, North Carolina State University, NCSU Campus, Box 8208, Raleigh, NC 27695. Email: firstname.lastname@example.org
* The National Center for Atmospheric Research is sponsored by the National Science Foundation.