Advanced Search
Article Contents

Implementation of a Conservative Two-step Shape-Preserving Advection Scheme on a Spherical Icosahedral Hexagonal Geodesic Grid

doi: 10.1007/s00376-016-6097-8

Get Citation+


Share Article

Manuscript History

Manuscript received: 13 April 2016
Manuscript revised: 06 August 2016
Manuscript accepted: 24 August 2016
通讯作者: 陈斌,
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Implementation of a Conservative Two-step Shape-Preserving Advection Scheme on a Spherical Icosahedral Hexagonal Geodesic Grid

  • 1. State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, China Meteorological Administration, Beijing 100081, China

Abstract: An Eulerian flux-form advection scheme, called the Two-step Shape-Preserving Advection Scheme (TSPAS), was generalized and implemented on a spherical icosahedral hexagonal grid (also referred to as a geodesic grid) to solve the transport equation. The C grid discretization was used for the spatial discretization. To implement TSPAS on an unstructured grid, the original finite-difference scheme was further generalized. The two-step integration utilizes a combination of two separate schemes (a low-order monotone scheme and a high-order scheme that typically cannot ensure monotonicity) to calculate the fluxes at the cell walls (one scheme corresponds to one cell wall). The choice between these two schemes for each edge depends on a pre-updated scalar value using slightly increased fluxes. After the determination of an appropriate scheme, the final integration at a target cell is achieved by summing the fluxes that are computed by the different schemes. The conservative and shape-preserving properties of the generalized scheme are demonstrated. Numerical experiments are conducted at several horizontal resolutions. TSPAS is compared with the Flux Corrected Transport (FCT) approach to demonstrate the differences between the two methods, and several transport tests are performed to examine the accuracy, efficiency and robustness of the two schemes.

1. Introduction
2. Implementation of TSPAS on a spherical IHP grid
  • The construction of an icosahedral grid starts from a raw icosahedron (inscribed inside a sphere), which has 12 vertices, 30 edges and 20 triangular faces. One-level finer grids are then generated by bisecting the edges of the former coarser grids (e.g., Baumgardner and Frederickson, 1985; Heikes and Randall, 1995a; Satoh et al., 2014). The nth bisection of the icosahedron is typically called grid level (Glevel) n. The division generates more spherical triangular faces that are formed by vertices. The vertices are then projected onto the sphere.

    The spherical hexagons/pentagons are constructed by connecting the circumcenters (Fig. 1) of six/five neighboring spherical triangles (that share a common vertex). Pentagons only occur at the 12 vertices of the raw icosahedron. A hexagon that is formed by the circumcenters of spherical triangles is also referred to as a Voronoi (e.g., Heikes and Randall, 1995a). This differs from the barycenter-type hexagon, which is formed by connecting the barycenters of neighboring spherical triangles (e.g., Tomita et al., 2001).

    Because spherical hexagons are not regular hexagons and are sort of distorted, grid optimizations are typically used for IHP grids to reduce truncation errors (Miura and Kimoto, 2005; Heikes et al., 2013). Many methods have been proposed for optimizing IHP grids. (Heikes and Randall, 1995b) slightly adjusted the intersection points of triangular and hexagonal/pentagonal grids to improve the accuracy. (Tomita et al., 2001) designed a spring grid system to reduce the grid noise that originates from the use of the Arakawa-A (Arakawa and Lamb, 1977) grid. Based on methods proposed by (Du et al., 1999) and (Du et al., 2003), (Ringler et al., 2008) introduced the Spherical Centroid Voronoi Tessellation (SCVT) for use in earth system modeling. The SCVT method is a more generalized surface tessellation approach that does not necessarily require starting from a subdivided icosahedron. When this method is used to optimize an IHP grid, the resultant grid can achieve the collocation and orthogonality properties, while the bisection property is not constrained (Miura and Kimoto, 2005). The SCVT optimization is used in this study. Open source icosahedral mesh tools are adopted to construct the mesh infrastructure (Renka, 1997; Peixoto and Barros, 2013), including the mesh generation and optimization.

    Figure 1.  Illustration of the C grid discretization with a spherical hexagonal grid. The red triangles are the vertices of the spherical triangles (formed by the dashed lines), and the solid circles denote the circumcenters of spherical triangles; the solid lines that connect the solid circles form a hexagon. The red triangles are also the centroids of the hexagons given the SCVT optimization. See section 2.2 for additional details.

  • For simplicity, only a two-dimensional case is considered in this study. The flux-form advection equation may be written in the following form: \begin{equation} \dfrac{\partial\rho q}{\partial t}=-\nabla\cdot(\rho q{V}) , (1) \end{equation} where ρ is a density-like quantity, q is a mixing-ratio-like quantity, and V is the velocity vector. The continuity equation can be achieved by simply setting q to 1, which gives \begin{equation} \dfrac{\partial\rho}{\partial t}=-\nabla\cdot(\rho {V}) . (2)\end{equation}

    Figure 1 shows the Arakawa-C control volume discretization in the context of an IHP mesh. The scalar values (e.g., ρ and q) are placed at the centroid of a control volume cell, which is also the vertex of a spherical triangle under the SCVT optimization. The velocity fields are defined at the midpoint of each hexagonal cell wall. Let Q=ρ q; then, the right side of Eq. (1) can be transformed into a line integral using Gauss's theorem: \begin{equation} \dfrac{\partial Q_0}{\partial t}=-\dfrac{1}{S_0}\sum_{i=1}^{i=N}{F}_i\cdot {n}_il_i , (3)\end{equation} where Q0 is the scalar value at the centroid of a hexagonal/pentagonal cell 0, S0 is the area of cell 0. A variable with a subscript denotes the target variable at a particular location (similar hereafter). ni is a unit outward normal vector of the cell wall i, N is 6 for a hexagonal cell and 5 for a pentagonal cell, li is the spherical length of the cell wall i, and is located between the vertex 0 and i. Fi has a form of \(\tilde{Q}_iV_i\), which denotes the mass flux across the cell wall i (\(\tilde{Q}_i\) is an imaginary quantity at the cell wall). Fi is assumed to be located at the intersecting point between the Voronoi's and triangle's edges (Fig. 1).

    With these configurations, the next step is to build the divergence operator with an appropriate spatial scheme to calculate the flux at the cell wall and an appropriate time scheme for the time integration. The transport equation can then be numerically solved. Combining TSPAS with the configuration of the IHP grid, it may be formulated as follows.

    The first step is to calculate an intermediate flux at each cell wall using the high-order scheme but with a slightly increased flux (* denotes an intermediate stage): \begin{equation} F_i^*={F}_i\cdot {n}_i=\beta_0\left[0.5U_i(Q_0+Q_i)-0.5|U_i|(Q_i-Q_0)\dfrac{|U_i|\Delta t}{\Delta m_{i0}}\right] ,(4) \end{equation} where β0 is a scalar that is placed at vertex 0 and is used to add a slight increment to the actual flux. The formula for calculating β0 will be given in section 2.3, and β0≥ 1 is a prerequisite. Ui=Vi· ni is the wind speed normal to the cell wall i, (the brown arrow in Fig. 1; the sign is outward positive). Q0 and Qi (Fig. 1) are scalar values defined at vertices 0 and i. ∆ t is the time step, and ∆ mi0 is the spherical length between vertex 0 and i. Next, we need to pre-update the scalar field at vertex 0 using Eqs. (3) and (4) and evaluate the updated \(Q_0^*[Q_0^*-Q_0-(∆ t/S)\sum_{i=1}^{i=N}F_i^*l_i\). An A value is then introduced and defined at each vertex as follows: \begin{eqnarray} A_0&=&(Q_0^*-Q_{\max})(Q_0^*-Q_{\min}) ,(5)\\ Q_{\max}&=&{\rm MAX}(Q_0,Q_1,Q_2,\ldots,Q_N) ,(6)\\ Q_{\min}&=&{\rm MIN}(Q_0,Q_1,Q_2,\ldots,Q_N) , (7)\end{eqnarray} where MAX and MIN denote taking the maximum and minimum values of vertex 0 and all of its neighboring vertices. By applying this operation globally, the entire A field will be evaluated.

    The second step is to choose an appropriate scheme to perform the real/final integration. If A values at both vertex 0 and i are less than zero, the flux at the cell wall i will be calculated using the high-order LW scheme as follows: \begin{eqnarray} F_{{\rm LW},i}&=&{F}_{{\rm LW},i}\cdot {n}_i\nonumber\\ &=&0.5U_i(Q_0+Q_i)-0.5|U_i|(Q_i-Q_0)\dfrac{|U_i|\Delta t}{\Delta m_{i0}} ;(8) \end{eqnarray} otherwise, the flux will be calculated using the UP scheme as follows: \begin{equation} F_{{UP},i}={F}_{{\rm UP},i}\cdot {n}_i=0.5U_i(Q_0+Q_i)-0.5|U_i|(Q_i-Q_0) . (9)\end{equation} The choice between the two schemes can be achieved by using an "if statement" or Eqs. (9)-(12) in (Yu, 1994). This is the entire integration procedure of TSPAS on the IHP grid. Note that by using Eqs. (8) or (9) for each cell wall, we are assuming a one-dimensional finite-difference operator in each direction. The flux in each direction of a given cell wall is calculated using the same equation [Eq. (8) or Eq. (9)] but with opposite signs. Thus, the global sum of the fluxes will vanish, and the conservation is guaranteed (it has been checked that the scheme conserves the tracer mass to machine precision). A concise flowchart of the computational procedure is presented in Appendix B.

  • Up to this point, only β is unknown. Starting from Eq. (3), the discretized two-step integration with a forward-in-time scheme from time level n to n+1 may be written as \begin{eqnarray} Q_0^*&=&Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_i^*l_i ,(10)\\ Q_{n+1,0}&=&Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_il_i . (11)\end{eqnarray} Meanwhile, considering Eqs. (4), (8) and (9), we may obtain \begin{eqnarray} F_{{\rm UP},i}&=&F_{{\rm LW},i}+0.5|U_i|(Q_i-Q_0)\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right) ,(12)\\ F_i^*&=&\beta_0F_{{\rm LW},i} . (13)\end{eqnarray} Note that Q0 and Qn,0 are the same unless otherwise explicitly mentioned. As stated previously, TSPAS combines the fluxes that are computed by the LW and UP schemes. Hence, only three possible situations may occur:

    (1) All of the fluxes are calculated by the LW scheme. Thus, \begin{equation} Q_{n+1,0}=Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_{{\rm LW},i}l_i .(14) \end{equation} Using Eqs. (10) and (13), we obtain \begin{equation} -\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_{{\rm LW},i}l_i=\dfrac{Q_0^*-Q_{n,0}}{\beta_0} . (15)\end{equation} Thus, \begin{equation} Q_{n+1,0}=Q_{n,0}+\dfrac{Q_0^*-Q_{n,0}}{\beta_0}=\dfrac{Q_0^*+(\beta_0-1)Q_{n,0}}{\beta_0} . (16)\end{equation}

    Note that as long as one cell wall around vertex 0 chooses the LW scheme (A0<0), there must be \(Q_\min\leq Q_0^*\leq Q_\max\). Because β0≥ 1, we obtain \begin{eqnarray} \dfrac{Q_{\min}+(\beta_0-1)Q_{\min}}{\beta_0}&\leq& \dfrac{Q_0^*+(\beta_0-1)Q_{n,0}}{\beta_0}\nonumber\\ &\leq& \dfrac{Q_{\max}+(\beta_0-1)Q_{\max}}{\beta_0} . (17)\end{eqnarray} Therefore, \(Q_\min\leq Q_n+1,0\leq Q_\max\). Thus, the shape-preserving property is guaranteed when β0≥ 1.

    (2) All of the fluxes are calculated by the UP scheme. The shape-preserving property is guaranteed when the stability condition ([(|Ui|∆ t)/∆ mi0]-1≤ 0) is satisfied.

    (3) A total of k fluxes are calculated by the UP scheme, and N-k fluxes are calculated by the LW scheme. Using Eqs. (11) and (12), we obtain \begin{eqnarray} Q_{n+1,0}&\!=\!&Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_{{\rm LW},i}l_i-\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|(Q_i-Q_0)\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] ,(18)\\ Q_{n+1,0}&\!=\!&\dfrac{Q_0^*\!+\!(\beta_0\!-\!1)Q_{n,0}}{\beta_0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\! \left[0.5|U_i|Q_i\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\!+\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i| Q_0\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] ,(19)\\ Q_{n+1,0}&\!=\!&\dfrac{Q_0^*}{\beta_0}+\!\left\{1\!-\!\dfrac{1}{\beta_0}\!+\!\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\right\}Q_0\!-\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|Q_i\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] . (20)\end{eqnarray} Because 1/β0≥ 0 and [(|Ui|∆ t)/∆ mi0]-1≤ 0 when the stability condition is satisfied, it only requires \begin{equation} 1-\dfrac{1}{\beta_0}+\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\geq 0 \ \ (21)\end{equation} to achieve \begin{eqnarray} &&\dfrac{Q_{\min}}{\beta_0}+\left\{1-\dfrac{1}{\beta_0}+\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)\right]l_i\right\}Q_{\min}-\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|Q_{\min}\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\leq Q_{n+1,0}\leq\nonumber\\ &&\dfrac{Q_{\max}}{\beta_0}+\left\{1-\dfrac{1}{\beta_0}+\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)\right]l_i\right\}Q_{\max}-\quad\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|Q_{\max}\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] . (22)\end{eqnarray} That is, \(Q_\min\leq Q_n+1,0\leq Q_\max\).

    Therefore, the shape-preserving property is guaranteed when Eq. (21) is satisfied. Thus, we require \begin{equation} \beta_0\geq \dfrac{1}{1-\frac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|\left(1-\frac{|U_i|\Delta t}{\Delta m_{i0}}\right)l_i\right]} \ \ (23)\end{equation} under the condition \begin{equation} 1-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|\left(1-\dfrac{|U_i|\Delta t}{\Delta m_{i0}}\right)l_i\right]>0 . (24)\end{equation} Let γi=|Ui|[1-(|Ui|∆ t)/∆ mi0]li and note that γi≥ 0 under the stability condition. We can find a value of γmax that satisfies \begin{equation} \gamma_{\max}={\rm MAX}(\gamma_1,\gamma_2,\gamma_3\cdots \gamma_N) . (25)\end{equation} Then, \begin{equation} 1-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}0.5\gamma_i>1-\dfrac{k\Delta t}{S_0}0.5\gamma_{\max} ;(26) \end{equation} namely, \begin{equation} \dfrac{1}{1-\frac{k\Delta t}{S_0}0.5\gamma_{\max}}>\dfrac{1}{1-\frac{\Delta t}{S_0}\sum_{i=1}^{i=k}0.5\gamma_i} . (27)\end{equation} Thus, we let \begin{equation} \beta_0={\rm MAX}\left(1,\dfrac{2}{2-\frac{k\Delta t}{S_0}\gamma_{\max}}\right) . (28)\end{equation}

    Eq. (28) is a generalized form of β that was used by (Yu, 1994), and it is a formula that deals with all possible situations. Ideally, the value of k may be determined for each vertex at each time step using iterative routines. However, this will substantially increase the computational burden. For simplicity, we use k=3 globally throughout the integration. Normally, a larger (smaller) k value leads to a larger (smaller) β value, and the scheme may choose the monotone scheme more (less).

3. Numerical tests
  • The numerical tests that are conducted in this study are given in (Williamson et al., 1992) and (Nair and Lauritzen, 2010). We use the solid body rotation test of (Williamson et al., 1992) and the deformational flow cases (1-4) of (Nair and Lauritzen, 2010). The numerical errors are calculated following (Williamson et al., 1992), including \begin{eqnarray} L_1&=&\dfrac{I(|Q_{\rm C}-Q_{\rm T}|)}{I(|Q_{\rm T}|)} ,(29)\\[1mm] L_2&=&\sqrt{\dfrac{I[(Q_{\rm C}-Q_{\rm T})^2]}{I[(Q_{\rm T})^2]}} ,(30)\\[1mm] L_{\infty}&=&\dfrac{\max(|Q_{\rm C}-Q_{\rm T}|)}{\max(|Q_{\rm T}|)} ,(31)\\[1mm] h_{\max}&=&\dfrac{\max(Q_{\rm C})-\max(Q_{\rm T})}{\Delta h} , (32)\end{eqnarray} \begin{eqnarray} h_{\min}&=&\dfrac{\min(Q_{\rm C})-\min(Q_{\rm T})}{\Delta h} ,(33)\\[1mm] \Delta h&=&\max(Q_{\rm T})-\min(Q_{\rm T}) , (34)\end{eqnarray} where I denotes the global integration, Q C denotes the computational solution, and Q T denotes the true/analytical solution, \(h_\max,h_\min\) and ∆ h represent deviations. The operators max and \(\min\) denote taking the maximum and minimum values on the globe, respectively. All of the numerical tests are performed on a unit sphere (radius = 1) for the normalization. The duration of integration is T=5, as suggested by (Nair and Lauritzen, 2010). We solve Eqs. (1) and (2) with a constant initial density field (ρ=1). Results are presented and diagnosed for Q in this paper. The initial spatial distributions of the tracer mixing ratio q are given in each subsection below. The resolutions of the icosahedral mesh are Glevel 4, 5 and 6, which have 2562, 10242 and 40962 points, respectively. The resolutions correspond to approximately 4°, 2° and 1° radians. The integration is divided into 600, 1200 and 2400 steps for Glevel 4, 5 and 6, such that the maximum Courant number has an approximately constant value at all resolutions. All of the results are post-interpolated to a uniform rectangular grid. We run each test with TSPAS and FCT. The error statistics from the two schemes are given in Tables 1-7. Figures that show the FCT's results are provided in a supplementary document.

  • In the solid body rotation test, a cosine bell is advected by prescribed ambient winds on the sphere. The wind field and initial scalar values are as follows: \begin{eqnarray} u&=&u_0(\cos\theta\cos\alpha+\sin\theta \cos\lambda\sin\alpha) ,(35)\\ v&=&-u_0\sin\lambda\sin\alpha ,(36)\\ u_0&=&\dfrac{2\pi}{T} ,(37)\\ q&=&\left\{ \begin{array}{l@{\quad}l} 0.5\left[1+\cos\left(\dfrac{\pi r}{R}\right)\right] & if\ r<R\\ 0 & if\ r\geq R \end{array} \right.,\quad (38)\end{eqnarray} where θ is the latitude, Λ is the longitude, α is the angle between the axis of the solid body rotation and the polar axis of the spherical coordinate system, and r is the great circle distance between vertex i and the center of the initial cosine bell. This can be calculated as follows: \begin{equation} r=\arccos(\sin\theta_i\sin\theta_0+\cos\theta_i\cos\theta_0\cos(\lambda_0-\lambda_i)) ,(39)\end{equation} where i and 0 denote the positions of the vertex and the center of the initial cosine bell, respectively. R=1/3 is the radius of the cosine bell. After the duration T, the cosine bell returns to its initial position. Tests are performed for α=0 and π/2.

    Figure 2 shows the results from the test with α=0, in which the cosine bell is located at the equator (in the sense of earth). The results are shown at t=0, T/2, T for three resolutions. The error measures for these tests are given in Table 1. The numerical precision increases with increasing resolution. No overshooting or undershooting values occur in this case for TSPAS. A comparison of the results of TSPAS and FCT shows that TSPAS generally has slightly smaller L1,L2, and L\infty errors at the different resolutions. FCT maintains larger maximum values, while TSPAS more effectively prevents the undershooting values.

    Figure 2.  Solid body rotation tests with $\alpha=0$ for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

    Figure 3 shows the results from a test with α=π/2, in which the cosine bell starts moving from the South Pole (in the sense of earth). The numerical precision (Table 2) is quite similar to those of the test with α=0 (Table 1), owing to the quasi-uniform mesh in the polar and equator regions. Note that although most of the results in Tables 1 and 2 appear to be the same, they are actually different if more decimal places are considered. The shape-preserving feature of TSPAS is also guaranteed. The differences between TSPAS and FCT are similar to those in Table 1.

    Figure 3.  Solid body rotation tests with $\alpha=\pi/2$ for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

    Figure 4.  Deformational flow test Case 1 for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T/2$ (bottom).

  • (Nair and Lauritzen, 2010) proposed a class of benchmark deformational flow test cases for the two-dimensional horizontal linear transport problems on a sphere. The scalar field undergoes severe deformation during the simulation, and the flow reverses its course at the half time and the scalar field returns to its initial position and shape at the final time step. This process makes the analytical solution available at the end of the simulation, and facilitates an assessment of the accuracy of the underlying transport scheme. We conducted four cases, as described in (Nair and Lauritzen, 2010). The selected initial scalar fields are two symmetrically located cosine bells that are defined as follows: \begin{equation} h_i(\lambda,\theta)=\dfrac{1}{2}\left[1+\cos\left(\dfrac{\pi r_i}{R}\right)\right]if\ r_i<R , (40)\end{equation} where R=0.5 is the base radius of the bells, and ri=ri(Λ,θ) is the great circle distance between (Λ,θ) and a specified center (Λii) of the initial cosine bell. The initial conditions consist of a background value b and two cosine bells with centers at position 1 and 2: \begin{equation} q(\lambda,\theta)=\left\{ \begin{array}{l@{\quad}l} b+ch_1(\lambda,\theta) & if\ r_1<R\\ b+ch_2(\lambda,\theta) & if\ r_2<R\\ b & {\rm otherwise} \end{array} \right., (41)\end{equation} where b=0.1 and c=0.9. The centers of the cosine bells depend on the choice of different test cases as follows.

    3.2.1. Case 1

    This is a non-divergent flow test case. The wind field is given as follows: \begin{eqnarray} u(\lambda,\theta,t)&=&k\sin^2\left(\dfrac{\lambda}{2}\right)\sin(2\theta)\cos\left(\dfrac{\pi t}{T}\right) ,(42)\\ v(\lambda,\theta,t)&=&\dfrac{k}{2}\sin\lambda\cos\theta\cos\left(\dfrac{\pi t}{T}\right) . (43)\end{eqnarray}

    Figure 5.  Deformational flow test Case 2 for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

    The flow parameter k is 2.4, and the centers of the initial cosine bells are kept at 11)=(π,π/3) and 22)=(π,-π/3). Figure 4 shows the results from this test case, and Table 3 shows the results from the error statistics. The numerical errors generally decrease as resolution increases. TSPAS shows slight undershooting values, which are generally smaller (absolute values, similar hereafter) than those produced by FCT. The L1, L2, and L\infty errors from FCT are generally smaller than those from TSPAS at different resolutions, and FCT maintains larger maximum values.

    Figure 6.  Deformational flow test Case 3 for Glevel 4 (left column), 5 (center column), and 6 (right column) at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

    3.2.2. Case 2

    This is also a non-divergent case. The wind field is given as follows: \begin{eqnarray} u(\lambda,\theta,t)=k\sin^2\lambda\sin(2\theta)\cos\left(\dfrac{\pi t}{T}\right) ,(44)\\ v(\lambda,\theta,t)=k\sin(2\lambda)\cos\theta\cos\left(\dfrac{\pi t}{T}\right) . (45)\end{eqnarray} The flow parameter k is 2. The centers of the initial cosine bells are kept at (5π/6,0) and 22)=(7π/6,0). Figure 5 shows the results for this test case, and Table 4 shows the error statistics. For both schemes, the numerical errors for this case are overall larger than those of Case 1, but the undershooting values are generally smaller than those of Case 1. The differences between TSPAS and FCT are similar to those in Case 1. TSPAS generally produces smaller undershooting values than FCT, while FCT has relatively smaller errors and maintains larger maximum values.

    3.2.3. Case 3

    This is a divergent case. The wind field is given as follows: \begin{eqnarray} u(\lambda,\theta,t)&=&-k\sin^2(\lambda/2)\sin(2\theta)\cos^2\theta\cos\left(\dfrac{\pi t}{T}\right),(46)\\ v(\lambda,\theta,t)&=&\dfrac{k}{2}\sin\lambda\cos^3\theta\cos\left(\dfrac{\pi t}{T}\right) . (47)\end{eqnarray} The flow parameter k is 1. The centers of the initial cosine bells are kept at 11)=(3π/4,0) and 22)=(5π/4,0). Figure 6 shows the results from this test case, and Table 5 shows the error statistics. This case produces the smallest numerical errors among all four cases with cosine bells, but generates the largest undershooting values. The differences between the two schemes are similar to those in Case 1 and Case 2.

    Figure 7.  Deformational flow test Case 4 for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

    3.2.4. Case 4

    \begin{eqnarray} u(\lambda,\theta,t)&=&k\sin^2\lambda'\sin 2\theta\cos\left(\dfrac{\pi t}{T}\right)+2\pi\cos\theta/T ,(48)\\[-0.5mm] v(\lambda,\theta,t)&=&k\sin(2\lambda')\cos\theta\cos\left(\dfrac{\pi t}{T}\right) ,(49) \end{eqnarray} \begin{eqnarray} \lambda'&=&\lambda-\dfrac{2\pi t}{T}. (50)\end{eqnarray}

    Figure 8. 

    This is a non-divergent case. The flow parameter k is 2. The centers of the initial cosine bells are kept at 11)=(5π/6,0) and 22)=(7π/6,0). Figure 7 shows the results from this test case, and Table 6 shows the error statistics. This case generates the largest errors of the four cases. The undershooting values for this case are quite small. TSPAS prevents undershooting values more effectively, while FCT still maintains larger maximum values. In contrast to the previous three deformational tests, TSPAS produces smaller L1, L2, and L\infty errors in this case.

    3.2.5. Case 3 with slotted cylinders

    To test the shape-preserving properties of the schemes in a more challenging way, we rerun Case 3 but with discontinuous slotted cylinders (Fig. 8). The errors statistics are presented in Table 7. Both schemes produce overshooting and undershooting values in this case. Similar to Case 3 with cosine bells, TSPAS generally produces smaller overshooting and undershooting values at various resolutions but with slightly larger errors.

    3.2.6. Rates of convergence

    To better check the order of accuracy, rates of convergence were estimated for each case, as shown in Table 8. The rate of convergence is computed by regressing the logarithms of error norms (L1 and L2) against the logarithms of the maximum angular distance between vertices at different resolutions. The regression coefficient is taken as the convergence rate. As shown in the table, the computed rates of convergence are around first-order accuracy, which is generally consistent with the fact that two schemes are composed by first-order and second-order schemes in practice.

  • In the above text, the generalized TSPAS is described. Several numerical tests were conducted to examine the simulation results, and a comparison between TSPAS and FCT was made. The table constituting Appendix C summarizes the overall performance of TSPAS and FCT. One may notice that these two schemes share a similar philosophy to satisfy the shape-preserving property; namely, find an optimal combination of schemes to maintain accuracy under the shape-preserving condition. The similar error norms and order of accuracy are consistent with the fact that both methods are composed of two same schemes. Nevertheless, they generate combinations in quite different ways.

    The FCT approach blends fluxes computed by two schemes at each edge. The blending coefficient is computed by explicitly prohibiting those conditions when maximum net anti-diffusive fluxes are higher than total inflow/outflow anti-diffusive fluxes. This sets the range of the blending coefficient from zero to one. In the TSPAS approach, the value of β is used to modify the LW-updated scalar field. Only when the pre-updated scalar values at both sides of a cell edge satisfy the shape-preserving rule will the scheme choose the LW scheme; otherwise, it will choose the monotone scheme. Therefore, the combination generated by TSPAS may contain more monotone and diffused solutions than FCT. This explains why TSPAS generally produces smaller undershooting/overshooting values, while FCT maintains larger maximum values of the scalar fields. Meanwhile, the procedures needed by TSPAS to generate a combination are more simplified, which helps TSPAS to be less time-consuming (FCT costs 1.3-1.4 times more than TSPAS).

    One may also notice that both schemes sometimes generate undershooting and overshooting. Reasons for this problem are complex, including the challenging flows and distributions of scalar values. From the aspect of the algorithm, the multidimensional application is a key reason. In the case of FCT, (Zalesak, 1979) discussed the impact of multidimensional applications on the monotonicity of the algorithm. We report its impact on TSPAS in the following text.

    As formulated in section 2.3, the value of k is a consequence of multi-dimensional applications. This parameter will not be introduced in a one-dimensional case. k is used to quantify how many edges will be computed by UP. This is a precondition because the scheme will always choose a combination of UP and LW. The value of k does not set a limitation on the number of fluxes that are computed by the UP scheme. Once k is set, it determines how likely a cell wall will choose UP (or not choose). But this possibility is not entirely associated with k, because β is the parameter that balances the options, and β is only partially contributed by k. The problem of using a fixed k value is that the options between two schemes, sometimes, may not rigidly satisfy the shape-preserving rule. Hence, the scheme may generate under/overshooting values, especially under the environment of severe deformational flows and sharp gradients. Tests on different k values suggest that the sensitivity is not significant. We use k=3 because it is the median of the number of edges of a hexagonal mesh.

4. Concluding remarks
  • In this study, a finite-difference advection scheme, TSPAS, was generalized and tested on an icosahedral hexagonal-pentagonal geodesic grid. The generalized scheme conserves the tracer mass and preserves monotonicity in a simple and efficient approach. Several transport tests were used to check the accuracy, efficiency and robustness of the generalized scheme when it is applied on an optimized IHP grid. The results demonstrate the feasibility of the generalized two-step scheme on an unstructured geodesic grid. Tests at various resolutions show that the numerical precision increases as the resolution increases.

    A comparison between TSPAS and FCT suggests that, generally, TSPAS produces smaller undershooting values, whereas FCT maintains larger maximum values. The errors derived from the L1, L2 and L\(\infty\) metrics show that the two schemes generate comparable errors. TSPAS generates smaller errors in two solid-body rotation tests and in the deformational test Case 4, whereas FCT produces smaller errors in the other four cases. Both schemes show an approximate first-order accuracy in idealized tests. Meanwhile, TSPAS is relatively more efficient in terms of the computational expense due to its simpler formulations, and it can be conveniently implemented on other unstructured grids. Moreover, for a practical global model application, the scheme can be easily adapted for a massively parallel environment given that it only requires a very limited node communication.




    DownLoad:  Full-Size Img  PowerPoint