Comparison of several high-order advection schemes for vertex-based triangular discretization

Abstract

Possible extensions of the standard 3rd (upwind) and 4th (centered) order vertex-centered finite-volume advection schemes on unstructured triangular meshes are explored. A finite-volume analog of the P₁ finite-element advection scheme with a consistent mass matrix is proposed, referred to as the compact scheme. In a simple two-dimensional test of advection by a shearing rotating flow, the compact scheme leads to a substantial error reduction compared with the traditionally used algorithms of the 3rd or 4th order, and can be nearly as accurate as the extension of the standard schemes to the 5th and 6th order. In full three-dimensional simulations of a turbulent baroclinically unstable flow, its eddy kinetic energy (EKE) only weakly increases for more accurate tracer advection. In terms of wall-clock CPU time, the compact scheme is the fastest despite the associated approximate inversion of mass matrix, and hence can be recommended.

Appendices

A Appendix: Polynomial reconstruction versus gradient estimate

In FV schemes based on field reconstruction, one writes a polynomial reconstruction around vertex i

$$ {\mathcal T}_{i}=a_{0}+a_{1}x+a_{2}y+a_{3}x^{2}+a_{4}y^{2}+a_{5}xy+..., $$

where x,y are the horizontal coordinates in the local reference frame associated with vertex i, the coefficients a_n, n = 0, 1, 2, ... are found by imposing a strong constraint ${\int \limits }_{S_{i}} {\mathcal T} dS=T_{i}S_{i}$ at vertex i and similar constraints at a sufficient number of neighboring locations (given by index j), but in a weak sense. The resultant least square problem

$$ {\mathcal L}={\sum}_{j\in N(i)} {w_{j}^{2}}|S_{j}^{-1}{\int}_{S_{j}}{\mathcal T} dS-T_{j}|^{2}+\lambda(S_{i}^{-1}{\int}_{S_{i}}{\mathcal T} dS-T_{i})=\min $$

is solved for the coefficients a_n, and their values are used to estimate fields and thus flux leaving the control volume. Here N(i) is a set containing a sufficient number of neighbor vertices; the weights w_j are (commonly) the inverse distances from vertex i to a neighbor vertex j, and λ is the Lagrange multiplier. In this polynomial expansion, a₀ is the tracer value at location i which generally differs from the area averaged value T_i. The estimates of the fluxes leaving control volumes should be carried out at properly selected Gaussian quadrature points at the boundary segments to ensure accuracy. All matrices needed to compute the coefficients a_n in terms of T_i and surrounding T_j (for j ∈ N(i)) are computed in advance and stored. This is the essence of schemes based on high-order reconstructions as proposed by Barth and Frederickson (1990) and developed further in many works that followed. Instead of neighboring locations, one may select any additional locations as, for example, in Chen et al. (2012) or Miura and Skamarock (2013). Generalizations may rely on several reconstruction stencils and WENO procedure (see, e.g., Dumbser and Käser2007).

On good-quality triangular meshes, a vertex has 6 nearest neighbors, which is sufficient for a quadratic reconstruction (QR). A scheme based on quadratic reconstruction was implemented in prototype FESOM2 (Danilov 2012) for median-dual control volumes. In its standard version, a reconstruction from the upwind vertex is used, but any combination of upwind and centered versions is possible. On 3D z-coordinate meshes, the number of neighbors might vary with depth, and quadratic reconstruction is replaced by the linear one when bottom topography is encountered. In practice, it turns out that QR schemes are nearly same accurate as the third- and fourth-order schemes based on the gradient estimate (GE) to be discussed below. This statement is similar to the conclusion in Skamarock and Gassmann (2011), who proposed a scheme that can be reformulated in terms of gradient estimates (see Miura 2013). In early implementation in FESOM, the QR schemes were about twice as expensive as the GE schemes explained further, which was an argument in favor of the latter. They are even better now.

A wider stencil is needed for higher order polynomial reconstructions, and, although such reconstructions open unlimited possibilities even on strongly distorted meshes (see examples in Dumbser and Käser 2007), their computational cost increases, as well as the halo size in parallel implementations. According to Skamarock and Menchaca (2010), the QR is optimal judged by accuracy per computational cost in standard tests where the convergence rate is about the second order.

The gradient estimate schemes achieve high order only on uniform meshes, and stay second order otherwise. We begin with a 1D simplification to explain them. Let the velocity u be positive and uniform. Traditionally, to estimate fluxes leaving a control volume around vertex i, one needs to estimate the tracer at i + 1/2 based on T_i and T_j of neighboring control volumes. The procedure can be cast in terms of gradients.

One first introduces the upwind and downwind estimates

$$ T_{i+1/2}^{-} = T_{i}+(h/2)G_{i+1/2}^{-},\quad T_{i+1/2}^{+} = T_{i+1}-(h/2)G_{i+1/2}^{+}, $$

where h is the uniform grid spacing. The simplest high-order choice is to estimate the gradients $G_{i+1/2}^{\pm }$ as

$$ G_{i+1/2}^{-}=(2/3)G_{i+1/2}^{c}+(1/3)G_{i+1/2}^{u},\quad G_{i+1/2}^{+}=(2/3)G_{i+1/2}^{c}+(1/3)G_{i+1/2}^{d}. $$

Here $G_{i+1/2}^{c}=(T_{i+1}-T_{i})/h$, $G_{i+1/2}^{u}=(T_{i}-T_{i-1})/h$, and $G_{i+1/2}^{d}=(T_{i+2}-T_{i+1})/h$ are the centered, upwind, and downwind estimates of gradient, respectively. Writing the flux

$$ F_{i+1/2}=(1/2)(F_{i+1/2}^{-}+F_{i+1/2}^{+})+(\lambda/2)(F_{i+1/2}^{-}-F_{i+1/2}^{+})\text{sign}(u), $$

one obtains the standard third-order upwind method for λ = 1 (GE3) or the forth-order-centered method for λ = 0 (GE4). Intermediate λ can be used to reduce dissipation and increase accuracy (λ = 0.15–0.25 works well in practice). One readily sees that the estimates $T_{i+1/2}^{\pm }$ can be rewritten as

$$ T_{i+1/2}^{-}=(1/2)(T_{i}+T_{i+1})-(1/6)h^{2}\delta^{2}T|_{i}, $$

(5)

and

$$ T_{i+1/2}^{+}=(1/2)(T_{i}+T_{i+1})-(1/6)h^{2}\delta^{2}T|_{i+1}, $$

(6)

where δ² stands for the operator of second derivative. They are the familiar estimates used to construct the standard third-order finite-difference method.

The advantage of the GE is that it can be generalized to arbitrary triangular meshes and also to higher orders (see, for example, Abalakin et al. 2002) by extending estimates to wider stencils. Consider the arrangement shown in Figs. 1 and 11.

Fig. 11

Schematic of the arrangement. Edge e with vertices i and j is characterized by the edge vector x_ij. u and d are up- and down-edge triangles, black circles are points where the continuation of edge e intersects the sides of u and d triangles, u₁,u₂ and d₁,d₂ are the vertices related to these sides. The gradients on triangles are computed based on three vertex values, and the gradients at vertices are obtained by averaging over neighboring triangles

We denote the vector connecting vertices i and j of edge e in Fig. 11 as x_ij. For each edge e, we store the indices to up-edge (u) and down-edge triangles (d) that contain the continuation of the edge line. Such storage is sufficient for the third- and fourth-order schemes. For the fifth and sixth order, one additionally stores the indices of vertices forming the edges intersected by the continuation of edge e in u and d triangles (u₁,u₂ and d₁,d₂) and the coefficients to interpolate the vertex gradients at these vertices to the intersection points. Similarly to the 1D case, one writes

$$ T_{ij}=T_{i}+(1/2)\textbf{x}_{ij}\textbf{G}_{ij},\quad T_{ji}=T_{j}-(1/2)\textbf{x}_{ji}\textbf{G}_{ji}. $$

Here we use pairs ij or ji to indicate up-edge or down-edge reconstructions. The reconstruction is done to the edge mid-point, because the boundary of control volume passes through it for median-dual control volumes on triangular meshes. As an alternative to median-dual control volumes, one can use Voronoi polygons of dual mesh. According to Abalakin et al. (2002), the following estimate for the gradients can be used

$$ \textbf{G}_{ij}=(1-\upbeta)\textbf{G}^{c}+\upbeta\textbf{G}^{u}+ \delta_{c}(\textbf{G}^{u}+\textbf{G}^{d}-2\textbf{G}^{c})+\delta_{d}(\textbf{G}^{j}+\textbf{G}^{u*}-2\textbf{G}_{i}) $$

and

$$ \textbf{G}_{ji}=(1-\upbeta)\textbf{G}^{c}+\upbeta\textbf{G}^{d}+ \delta_{c}(\textbf{G}^{u}+\textbf{G}^{d}-2\textbf{G}^{c})+\delta_{d}(\textbf{G}^{i}+\textbf{G}^{d*}-2\textbf{G}_{j}). $$

In these expressions, G^u and G^d are the gradients on triangles u and d, $\textbf {x}_{ij}\textbf {G}^{c}=T_{j}-T_{i}$ is the centered estimate, Gⁱ and G^j are the estimate at vertices i and j, and G^u∗ and G^d∗ are the vertex gradients interpolated to the edge continuation intersection points (see Fig. 11). The gradients on triangles are computed assuming linear interpolation. The gradients on vertices are computed as area-weighted averages of the gradients over neighboring triangles, or equivalently, by applying the divergence theorem. The selection β = 1/3 and δ_c = δ_d = 0 leads to a third/forth-order methods depending on the upwind parameter λ in the expression for fluxes. Keeping β = 1/3 but taking δ_c =-1/30 and δ_d =-2/15 leads to a fifth/sixth-order method, once again, depending on the value of the upwind parameter. On 3D z-coordinate meshes, any of or both u and d triangles can be absent for edges touching boundaries or bottom topography. We replace u and d gradients by vertex gradients at i and j in this case. Technically, the fifth/sixth-order method adds computations of vertex gradients and increases the amount of computations needed to estimate fluxes. However, logistics related to the z-coordinate bottom is expensive for the third/fourth order, and even more so for the fifth/sixth-order method. The third/fourth-order method requires an extended halo for triangles, and the fifth/sixth-order method needs additionally an extended halo for vertices including neighbors of neighbors. We note that if we were only interested in the methods of third or fourth order, the implementation of Skamarock and Gassmann (2011), based on computing second-order derivatives, as in Eqs. 5 and 6, is the same convenient as operations on gradients. The results are not equivalent, but similar. The coefficients and indices to vertices needed to compute the second-order derivatives can be computed before, but the z-coordinate bottom introduces complications as in the other case.

B Appendix

Table 4 L₂ errors in the 2D test case