A Primitive Variable Discrete Exterior Calculus Discretization of Incompressible Navier-Stokes Equations over Surface Simplicial Meshes

A conservative primitive variable discrete exterior calculus (DEC) discretization of the Navier-Stokes equations is performed. An existing DEC method (Mohamed, M. S., Hirani, A. N., Samtaney, R. (2016). Discrete exterior calculus discretization of incompressible Navier-Stokes equations over surface simplicial meshes. Journal of Computational Physics, 312, 175-191) is modified to this end, and is extended to include the energy-preserving time integration and the Coriolis force to enhance its applicability to investigate the late time behavior of flows on rotating surfaces, i.e., that of the planetary flows. The simulation experiments show second order accuracy of the scheme for the structured-triangular meshes, and first order accuracy for the otherwise unstructured meshes. The method exhibits second order kinetic energy relative error convergence rate with mesh size for inviscid flows. The test case of flow on a rotating sphere demonstrates that the method preserves the stationary state, and conserves the inviscid invariants over an extended period of time.


Introduction
Recently, structure preserving numerical methods that satisfy the mathematical and physical properties of the governing equations are increasingly becoming popular [1,2,3]. One of such methods is discrete exterior calculus (DEC). It uses integral quantities / discrete forms as the primary unknowns and transforms the continuous partial differential equation (PDE) system into a discrete algebraic system such that all errors/approximations occur only at the algebraic level and not in the approximation of the differential operators of the PDE [4]. Thus, DEC is known to preserve some of the underlying physical and mathematical properties of the continuous PDE system.
DEC is the discrete counterpart of smooth exterior calculus, which is more natural tool than tensor calculus in numerous situations. The theory of exterior calculus and differential forms dates back to the pioneering work of Poincaré [5] (who presented the notion of line and surface integrals), Cartan [6] ( who formalized the notion of differential forms), and Goursat [7] ( who noted that the three theorems of vector calculus, viz. Gauss, Green, and Stokes are all special cases of a generalized Stokes theorem for differential forms). Subsequently, there were further theoretical developments (for example, see [8,9,10]). However, there were limited applications of this theory to solve physical problems, until the theory of discrete exterior calculus was developed (see [11,12,13] and the references therein). As this theory advanced, DEC made its way to the applications such as computer graphics and computer vision [14,15,16,13,17,18,19] to mitigate the effects of numerical viscosity that causes detrimental visual consequences, and fluid flows [20,21,22,23] to preserve the physical properties of the governing equations, i.e, to conserve secondary quantities, such as kinetic energy and/or discretely preserving Kelvin's circulation theorem. Developing such physically conservative discretizations, for Navier-Stokes equations for example, is favorable for many physical applications such as turbulent flows [24] and shallow-water simulations [25] to avoid undesirable numerical artifacts.
Apart from DEC, some staggered mesh methods have also demonstrated conservation of secondary (kinetic energy and vorticity) physical quantities in addition to the primary ones (mass and momentum) [26]. Nicolaides [27] and Hall et al. [28] developed staggered mesh methods, also known as covolume methods, for unstructured meshes. The covolume discretization commences by taking dot product of the momentum equation with the unit normal to each face of the triangular / tetrahedral elements, reducing the velocity vector to a scalar flux (equal to the mass flux across the face for an incompressible flow with constant density) defined on each face; and the pressure is defined at the circumcenter of each triangular / tetrahedral element. Nicolaides [29] estimated the accuracy of the covolume method to be second order for a mesh with all its triangular elements having the same circumradii (i.e., structured-triangular mesh), and to be first order otherwise. Perot [30] demonstrated the conservation properties of the covolume method. The divergence form of the Navier-Stokes equations conserved the momentum and kinetic energy both locally and globally, and the rotational form conserved the circulation and kinetic energy locally and globally for both 2D [30] and 3D [24] discretizations. A main advantage of DEC discretizations is their applicability to simulate flows over curved surfaces, unlike the covolume approach [21]. Due to its intrinsic coordinate independent nature, DEC based implementations work on both planar domains and curved surfaces without any change. Its independence of the embedding space dimension as well as the details of the embedding characterizes a key difference between DEC and the covolume method. In the latter, face normal vectors, with which dot product of the momentum equation is performed for the discretization, may not be unique for a simplicial mesh approximation of a curved surface.
Here, we extend an existing stream function formulation DEC scheme [21], to the primitive variable formulation. To compute the aerodynamic forces for the flows past bluff / streamlined bodies, the stream function formulation requires a posteriori computation of the pressure field from the stream function field, which may introduce some error. The primitive variable formulation computes the pressure field as a part of the DEC solution and is more relevant for such physical applications. Moreover, the method is extended to include the Coriolis force, in order to enhance its applicability to vortical flows on rotating surfaces which are usually used to model planetary flows. In addition, an energypreserving time integrator [17] is incorporated, to enhance applicability of the method to study the late time behavior of such (vortical flows on rotating surfaces) initial value problems.

Primitive Variable DEC Formulation
A primitive variable DEC formulation is presented here, extending an existing stream function method [21]. First, the Navier-Stokes equations are expressed in exterior calculus notation, and then the DEC discretization is derived. The resultant degrees of freedom are the primitive variables, namely the velocity 1-form and the pressure 0-form. We recommend for readers unfamiliar with DEC to read through the appendix for a primer on discrete exterior calculus before reading through this section.

Exterior Calculus Expression of the Navier-Stokes Equations
For an incompressible flow of a homogeneous fluid with unit density, on a compact smooth Riemannian surface, the vector calculus notation of Navier-Stokes equations in a rotating frame of reference are as follows.
where v is the tangential surface velocity, p is the effective surface pressure (which includes the centrifugal force), µ is the dynamic viscosity, ∆ dR is the surface Laplace-DeRham operator, κ is the surface Gaussian curvature, ∇ v is the covariant directional derivative, grad S is the surface gradient, div S is the surface divergence, f is the Coriolis parameter, andk is the unit vector in the direction of the axis of rotation. The Laplace-DeRham operator and the Gaussian curvature term are a result of the divergence of the deformation tensor and the non-commutativity of the second covariant derivative in curved spaces [22]. Following Mohamed et al. [21], with the additional curvature term −2µκu and Coriolis force fk × v, the 2D coordinate invariant form of equations (1) -(2), for constant κ (for convenience, its variation is addressed later in Sec. 2.2.2), read Here, u is the 1-form velocity, d is the exterior derivative, * is the Hodge star operator, ∧ is the wedge product operator, p d = p + 1 2 (u (v)) is the effective dynamic pressure 0form, N is the space dimension, and the Coriolis force fk × v = * (u ∧ * f dual2 ) [31], where f dual2 is the dual 2-form corresponding to f . Equations (3) -(4) correspond to smooth exterior calculus expressions, and are still continuous.

Discrete Exterior Calculus Expression
The discretization consists of the domain discretization, and discretization of the governing equations, i.e., discrete exterior calculus expression of the governing equations (equations (3) and (4)).

Domain discretization
Two-dimensional (2D) domain discretization is considered assuming a primal simplicial mesh. The objects of this simplicial mesh are the triangles (primal 2-simplices), edges of the triangles (primal 1-simplices), and the vertices / nodes (primal 0-simplices). A circumcentric dual mesh is considered, which consists of dual polygons, the edges of the dual polygons, and the dual vertices / nodes. The dual to a primal triangle is its circumcenter (which is a dual node), that to a primal edge is the perpendicular dual edge connecting the circumcenters of the two triangles sharing this primal edge, and that to a primal node is the polygon formed by the edges dual to the primal edges connected to this node. For a triangular mesh representing a curved surface, the dual edges are kinked lines and the dual cells are non-planar. Counterclockwise orientation of both the primal mesh (triangles/primal 2-simplices) and dual mesh (polygons) is assumed to be positive. The choice of the orientation of the primal edges (primal 1-simplices) is arbitrary, however it induces the dual edges orientations. The orientation of the dual edges is derived simply by rotating the primal edge orientation 90 degrees along the triangles orientation (i.e., counterclockwise). Moreover, there is no requirement such as the triangles need to be Delaunay [32].

Discrete exterior calculus expression of the governing equations
First we need to specify the correspondence between the discrete variables (the forms) and the mesh objects (vertices, edges etc.). This requires that we define the discrete forms, or discretize the smooth forms. We then determine the discrete operators and substitute the smooth exterior calculus operators with their discrete counterparts. For the current discretization, we define the velocity 1-form u in all terms (except one in the nonlinear term which will be discussed later) on dual edges in equation (3), i.e., we chose the velocity 1-form u to be a dual 1-form. The dual 1-form represents the velocity integration along the dual edges and hence corresponds to the mass flux across the primal edges. With this choice, for consistency, each term (time derivative, convection, diffusion, pressure gradient) in equation (3) must be a dual 1-form. This choice of the velocity 1-form requires the dynamic pressure 0-form p d to be defined at the dual nodes or the triangles circumcenters (to be defined as a dual 0-form), so that the pressure gradient term dp d is a dual 1-form, consistent with the choice of u as the dual 1-form. For consistency, since the convective term * (u ∧ * du) is a dual 1-form, the wedge product term (u ∧ * du) is a primal 1-form (defined on primal edges). Since u is a 1-form, and hence * du is a 0-form (in 2D), and therefore the operands of the wedge product are a 1-form and a 0-form, such that the resultant (primal 1-form) must be defined on primal edges. This requires that both wedge product operands are defined on primal mesh objects. Accordingly, the velocity 1-form in the first operand (u) of the wedge product term is defined on the primal edges, whereas the velocity 1-form in the second operand ( * du) is defined on the dual edges, so that the 0-form * du is defined on the primal nodes. Following Mohamed et al. [21], the definitions and expressions for the discrete d, * , and ∧ operators, the DEC discretization of equation (3) now reads (in 2D) Here, U is the vector containing the discrete dual velocity 1-forms for all mesh dual edges, V is the vector containing the discrete primal velocity 1-forms for all mesh primal edges, and P d is the vector containing discrete effective dynamic pressure 0-forms for all mesh dual vertices, ∆t is the discrete time interval between the current time instant n + 1 and the previous time instant n. Note that the time-centering of all terms except the one involving the time derivative is discussed later in sections 2.2.3 and 2.2.4. The matrix W V is defined as W V = 0.5diag(V ) |d 0 | and is a sparse N 1 ×N 0 matrix, where N 1 is the number of primal edges and N 0 is the number of primal nodes. Moreover, the term * −1 0 −d T 0 U represents vorticity and is a vector with N 0 entries. The term W V * −1 0 −d T 0 U is equivalent to taking the average, at each primal edge, of the vorticity * −1 0 −d T 0 U evaluated at the edge's nodes and multiply it by the edge's tangential velocity 1-form v. Thus, this term evaluates the first wedge product. Moreover, since * f dual2 is a primal 0-form, the 1-form velocity u in the second wedge product also has to be the primal 1-form for consistency. Thus the term * 1 W V * −1 0 f dual2 evaluates the second wedge product. The tangential velocity 1-form v at each primal edge is evaluated following Mohamed et al. [21].
The boundary of the dual 2-cells to the primal nodes on the domain boundary includes primal boundary edges. Since −d T 0 matrix is the boundary operator of these dual 2-cells, it is complemented by an additional operator d b to account for the primal boundary edges. The operation −d T 0 U evaluates the circulation of the velocity 1-forms u along the dual 2-cells boundaries. The operation d b V complements this circulation for the dual 2-cells part of whose boundaries consists of primal edges. The matrix d with all entries made positive and 1 is a vector of ones with N 2 (number of primal 2-simplices / triangles) entries, and diag (·) is a diagonal matrix consisting of the enclosed vector entries.
Since the mass flux primal 1-form (the mass flux across the primal edge) is defined as the integral of the normal velocity (normal to the primal edge) along the primal edge, the vector containing mass flux primal 1-form u * for all mesh primal edges can be expressed as U * = − * −1 1 U or U = * 1 U * . After this substitution, Eq. (5) reads * 1 Applying the Hodge star * −1 1 to equation (6), and considering the property * −1 In general, κ may not be uniform. For domains with non-uniform κ, it may be treated as a primal 0-form (defined at the primal vertices). Now the curvature term (2µκU * ) is modified and Eq. (7) reads Here, K is half times a matrix containing values of κ at the nodes of each edge, and it accounts for the wedge product of the primal zero form κ with the primal 1-form u * .
With reference to temporal discretization, we first implemented the Euler first order time integration scheme. Later, in order to enhance the ability to study the late-time behavior of initial value problems, we implemented the energy-preserving time integration [17]. Both time integration schemes are discussed below.

Euler first order time integration
We treat the viscous, curvature and pressure gradient terms implicitly, and the convection and Coriolis force terms explicitly. In order to obtain a linear representation for the convective term, we evaluate the tangential velocity 1-forms V through an interpolation of previously-known normal velocity 1-forms U . Therefore, we have or where I is an identity matrix, and The continuity equation is expressed as Equations (10) -(12) represent a set of linear equations with U * and P d as the solution vectors, which is solved along with usual boundary conditions depending on the considered test case.

Energy-preserving time integration
The midpoint integration is assumed for the update in time of the advection, diffusion, curvature, pressure gradient and Coriolis force terms [17], i.e., these terms are evaluated at time n + 1 2 . Thus we have, Rearranging the terms, we have with Equation (14) along with equation (12), and the relevant boundary conditions (depending on the test case) form a non-linear system of equations. These are solved using Picard's iterative method, with previous time-step value used as an initial guess for each time-step and a convergence criterion of L 2 norm on the residual of R 2 ≤ 10 −8 . It is worth noting that the current discretization is different from that in [17] for the convective term spatial discretization and the inclusion of both surface curvature effect and Coriolis force.

Results and Discussions
In order to assess the accuracy of the scheme, several test cases are performed. All test cases of flow past a cylinder are solved using the Euler time integration scheme. The energy-preserving time integration scheme is then employed for the remaining test cases.

Flow Past a Stationary Circular Cylinder
The computational domain consists of a rectangular plane with a circular cylinder embedded on it. The circular cylinder contour embedded on the surface is generated by intersecting a cylinder of diameter D = 1 with the surface. The domain lengths upstream and downstream of the cylinder are 10D and 60D, respectively, and the domain width on either side of the cylinder is 10D. These domain sizes are usual in computational investigations of flow past a cylinder. Uniform normal velocity of unity is assumed at the inlet. Outflow boundary condition is assumed at the outlet. Free-slip boundary condition is considered on the top and bottom boundaries, whereas no-slip boundary condition is assumed on the cylinder wall. The computational domain is discretized using a triangular mesh (consisting of approximately 160k elements) with locally refined mesh in the vicinity of the cylinder and in the wake. Reynolds numbers (Re = ρu ∞ D/µ) = 40, 100, and 1000 are considered. Here, ρ and u ∞ are the fluid density (taken to be unity) and the free stream velocity at the inlet, respectively.
The flow exerts pressure and viscous forces on the cylinder, and the net force acting on the cylinder is expressed as [33] where p is the static pressure, n is the unit normal facing outward of the cylinder, and ω is the vorticity vector (normal to the embedding surface). The integration is carried out over the total surface area S (perimeter here) of the cylinder. The discrete expression of the force vector reads where the summation is over the cylinder N boundary edges. Here, global (x, y) coordinate system is considered with the cylinder axis intersection with the embedding surface as the origin. The x and y coordinates are along the streamwise and cross-streamwise directions, respectively. The static pressure p k is in the triangles adjacent to the boundary edges, n k is the unit vector normal to the edge, and l k is the edge length. The vorticity vector is ω k = n k,s |ω k |, where n k,s is the unit normal to the embedding surface at the boundary edge midpoint. The vorticity magnitude |ω k | is calculated by averaging the known vorticity on both of the edge nodes. The drag force F d and the lift force F l are respectively the streamwise and cross-streamwise components of the force vector F , i.e., are the x and y components of the force vector F . We calculate the drag coefficient, lift coefficient, and Strouhal number as and St = f D/u ∞ , respectively. Here, ρ is the fluid density (taken to be unity), u ∞ is the flow velocity at the inlet boundary of the domain (also taken to be unity), and f is the vortex shedding frequency. Table 1 reports the drag coefficient mean values. The computed values are 1.605, 1.386, and 1.544 at Re = 40, 100, and 1000, respectively. They agree well with that reported in the literature (e.g., at Re = 40: 1.6 [34], 1.618 [35]; at Re = 100: 1.38 [34], 1.37 [36]; at Re = 1000: 1.5144 [37], 1.5 [38]). Table 1 also reports the RMS values of the lift coefficient. We compute the values to be 0.248, and 1.081 at Re = 100, and 1000, respectively. They are within the range of values reported in the literature (e.g., at Re = 100: 0.233 [39]; at Re = 1000: 1.0494 [37]). In passing, it is worth noting that the predicted values of the mean drag coefficient and RMS lift coefficient using the stream function formulation [21] are 1.337 and 0.171, respectively at Re = 100. The drag coefficient is predicted correctly, but the lift coefficient is under-predicted. The under-prediction of the lift coefficient may be, as already discussed before (see section 1), because of the error induced in posterior computation of the pressure field. Table 1 reports the Strouhal number values. The Strouhal number is computed as 0.171, and 0.239 at Re = 100, and 1000, respectively, which is within the range of values reported in the literature (e.g., at Re = 100: 0.172 [34], 0.166 [39]; at Re = 1000: 0.237 [37], 0.244 [40]). Figure 1 shows the plots of time averaged pressure coefficient C p = (p − p ∞ ) / 1 2 ρu 2 ∞ as a function of the angular position on the cylinder surface (θ) at Re = 100. Here, p ∞ is the free stream pressure. The current results agree well with literature.
There is an attached stationary pair of vortices in the cylinder wake at Re = 40, the length of which is defined as the distance between the cylinder center and the end of the stationary vortex. The vortex pair length is determined to be 2.78, which is within the range reported in the literature (e.g., 3.01 [42], and 2.64 [44]).

Flow Past a Rotating Circular Cylinder
The same flow configuration as described in section in 3.1 is used, except that the cylinder now has a non-zero tangential velocity. Reynolds number of 200 is considered. Non-dimensional tangential velocities α = (D/2) ω/u ∞ = 2.5 and α = 3.0 are employed. Here, ω is the tangential velocity of the cylinder. Since the cylinder is rotating, there is non-zero lift force acting on the cylinder. The lift coefficient (C L ) (as already defined in section 3.1) values are reported in table 3. For α = 2.5 and 3, the values of C L are found to be −7.641 and −10.413, respectively, which agree with the values −7.63 and −10.329 reported in [45]. Figure 2 shows the plots of pressure coefficient as a function of the angular position on the cylinder surface (θ) for α = 2.5 and 3.0 at Re = 200. They agree with that reported in [45].

Flow Past a Stationary Circular Cylinder on a Locally Curved Surface
The same flow configuration as the one described in section 3.1 is employed except that now there is a bump present in the cylinder wake. The bump is described by where the parameters describing the bump location read: x 0 = 6.0, y 0 = 0.5, h = 1.6, 1/ (2r 2 0 ) = 0.4. The global Cartesian coordinate system as described in section 3.1 is employed, with an additional coordinate z to describe the bump. Figure 3 shows an instantaneous cylinder wake at Reynolds number, based on cylinder diameter, Re = 100. A comparison with the case without bump shows some disruption of the wake due to the vorticity generated by the bump. However, the vorticity generated at the cylinder wall dominates and the effect of the presence of bump is insignificant. The mean drag

Flow Past Stationary Square and Triangular Cylinders on a Spherical Surface
The computational domain consists of a square or an equilateral triangular cylinder contour, both with edge length unity, embedded on a spherical surface of radius R = 12. One of the edges of the square cylinder (contour) faces the upstream flow and the opposite edge faces the downstream flow. The vertical edge of the triangular cylinder (contour) faces the downstream flow. The domain size is similar to the one described in section 3.1, except that now the lengths are geodesic distances. The same boundary conditions as described in section 3.1 are employed. Reynolds number, based on the edge length, Re = 100 is considered. Figure 4 shows the instantaneous vortex streets, along with that for the case of flat embedding surface for a reference. Qualitatively, the vortex streets on the spherical embedding surface look similar to that on the flat embedding surface. For the square cylinder case, the mean drag coefficient, RMS lift coefficient and the Strouhal number are computed to be 1.477, 0.2 and 0.148, respectively. We compare them to the corresponding values for the flat embedding surface, and the relative difference is found to be 1.2% and 0.67%, respectively for the mean drag coefficient and the Strouhal number. The RMS lift coefficient is found to be identical. Similarly, for the triangular cylinder, the mean drag coefficient, the RMS lift coefficient and the Strouhal number are computed to be 1.77, 0.305 and 0.2, respectively. The corresponding difference is found to be 2.69%, 4.98% and 1.96% relative to the flat embedding surface case.

Flow Past an Airfoil
Flow over an airfoil, NACA 0012 at different angles of attack α is considered. The computational domain consists of the contour of the airfoil embedded on a rectangular plane. The domain lengths upstream and downstream of the airfoil leading edge are 5C and 15C, respectively, and the domain width on either side of the leading edge is 5C. Here, C stands for the chord length. The flow is set to be uniform at the inlet (the left boundary) and outflow boundary condition is applied at the outlet (the right boundary). The free stream boundary conditions are employed at the top and bottom boundaries of the domain, and a no-slip boundary condition is employed at the airfoil wall boundary. Figure 5 shows the mean pressure distribution over the upper and lower surfaces of the airfoil at Re = 1000. The present result shows good agreement to those obtained by [46].

Taylor-Green Vortices
The Taylor-Green vortices are simulated on a flat square domain of dimension [-0.5,0.5] both in x and y coordinate directions. The non-dimensional analytical solution for the decay of Taylor-Green vortices with time in 2D reads [47] where u x and u y are the x and y components of non-dimensional velocity respectively, x and y are the non-dimensional coordinates, Re is the Reynolds number, and t is the non-dimensional time. The Reynolds number is defined as the ratio of viscous diffusion time to the circulation (turnover) time, i.e., Re = L 2 /ν L/V . Here, L is the characteristic length scale, V is the characteristic velocity scale, and ν is the kinematic viscosity. The   Figure 8 shows the convergence with the mesh element size of the L 2 -norm of the velocity 1- at simulation time t = 10. Here the superscripts analytical and numerical stand for the analytical and numerical solutions, respectively. We observe that the flux error converges with a second order rate for the structured-triangular mesh, and with a rate somewhat better than first order (i.e., 1.45) otherwise, which agrees largely with that reported in [21]. Figure 8 also shows the convergence of the L 2 -norm of the velocity magnitude error. The velocity vector field was interpolated from the flux distribution using Whitney maps [21]. The velocity error converges with a rate slightly better than first order (i.e., 1.2) for both the structured triangular mesh as well as the Delaunay mesh. For the structured triangular meshes, although the flux error converges with a second order rate, the velocity error converges at a rate slightly better than first order, due to the first order interpolation of velocity field from the flux distribution. We note that the stream function formulation [21] also exhibits similar second order flux error convergence rate for the structured-triangular meshes and first order flux error convergence rate for the otherwise unstructured meshes, and first order velocity error convergence rate for both the structured-triangular meshes as well as the unstructured meshes. For the post processing, the vorticity distribution (here, there is only one component of vorticity normal to the domain plane) was determined from the flux distribution as ω = * −1 Figure 9a shows the vorticity distribution at simulation time t = 10. Figure 9b shows, at simulation time t = 10, the distribution of the y-velocity along the horizontal domain center line along with analytical solution, and the two agree well.

Double Periodic Shear Layer
An inviscid flow and a square domain of unit edge length is assumed. Figure 10a shows the initial flow, which is a shear layer of finite thickness with a small magnitude of  vertical velocity perturbation. The expression of the initial velocity field is [48] u x = tanh ((y − 0.25) /ρ) for y ≤ 0.5 tanh ((0.75 − y) /ρ) for y > 0.5 (21) u y = δ sin (2πx)  Figure 10 shows the evolution of the vorticity field with time, using the finest mesh. At time t = 0.8, two well resolved vortices appear. Similar to that reported in Mohamed et al. [21], the shear layer connecting the coherent vortices becomes thinner with time and within a finite time interval reach the resolution of the mesh after which dispersion error manifests as mesh level oscillations. The vorticity distribution plots in figure 10 exhibit similarities with that reported previously by Mohamed et al. [21], and Bell et al. [48]. Figure 10d the convergence of the kinetic energy relative error KE(0)−KE(t) with the mesh size at time t = 2.0. Here, the total kinetic energy (KE) is calculated as 1 2 Ω v ·vdΩ, where the integration is performed over the entire simulation domain, where the velocity vector (v) is calculated in each triangular element via Whitney map interpolation, as discussed before. As expected, the kinetic energy relative error converges in a second order fashion with the mesh size. It is worth noting that the stream function formulation [21] also exhibits similar second order convergence rate of the kinetic energy relative error. Kinetic energy relative error is 0.2% for the coarsest mesh (equivalent to 64×64 Cartesian mesh) and it is only 0.0039% for the finest mesh (equivalent to 512×512 Cartesian mesh). It is 0.06% for the mesh with 32768 triangular elements (equivalent to a 128×128 Cartesian mesh), almost one order of magnitude lower than a second order collocated mesh scheme [48] using almost the same mesh size.

Motion of Harmonic Waves on a Rotating Sphere
There exists a stationary solution to the motion of harmonic waves in the Earth's atmosphere. These waves are the so-called Rossby waves. The stationary solution is [49] where, ψ is the stream function, B = 2Ω l(l+1)−2 , φ is longitude, θ is the colatitude angle, A is an arbitrary constant, P m l is associated Legendre polynomial of degree l and order m, R is earth's radius, and Ω is the rate of earth's rotation. An inviscid flow and f = 2Ω cos θ are considered. The stream function distribution is shown in figure 11 for A = 200, l = 7, m = 6, R = 6.371 × 10 6 , and Ω = 7.2921 × 10 −5 . Figure 11a shows the stationary solution at time t = 0, and 11b shows the one at time t = 3.144 × 10 6 s or 36.39 days. There is no distortion of the stream function distribution and the stationary solution is preserved. However, there is a westward phase displacement of 28.6 • longitude after 36 days. (5.7 • longitude after 7 days). Similar westward phase displacement of 7 • longitude after 8 days has been reported in [50].   Figure 12 shows various inviscid invariants as a function of time, and demonstrates conservation of these invariants for an extended period of time. It is observed that energypreserving time integration has a superior conservation of total kinetic energy and the total enstrophy as compared to the Euler first order time integration. The total kinetic energy is already defined in section 3.7. The total enstrophy is calculated as 1 2 Ω ω 2 dΩ, where the integration is performed over the entire simulation domain, and ω is the vorticity.

Conclusions
We presented a primitive variable DEC formulation of Navier-Stokes equations on surface simplicial meshes. Here, the mass flux 1-forms at all mesh primal edges, and the dynamic pressure 0-forms at the mesh dual nodes (triangle circumcenters) are the degrees of freedom. The discrete operators emerging from the DEC discretization are discussed. Two time integration schemes, i.e., Euler first order, and energy-preserving, are implemented. The Coriolis force is incorporated to facilitate the use of the method to study flows on rotating surfaces. The method exhibits second order flux error convergence rate with the mesh size for the structured triangular meshes, and a rate somewhat better than first order (i.e., 1.45) for unstructured meshes. Moreover, the method exhibits slightly better than first order (i.e., 1.2) error convergence rate with the mesh size for the interpolated velocity vector for both the structured triangular as well as unstructured meshes. The method exhibits second order kinetic energy relative error convergence rate with mesh size for inviscid flows. The method preserves a stationary state for the flow on a rotating sphere over an extended period of time. The energy-preserving time integrator exhibits superior conservation of inviscid invariants such as total kinetic energy and total enstrophy for an extended period of time.
Other key concepts are the exterior derivative (d), the Hodge star operator ( * ), and the exterior product or wedge product (∧). The exterior derivative implies multidimensional differentiation, and is a generalization of the usual curl, divergence and gradient operations of calculus. It operates on a form and results in to the next higher dimensional form. The gradient of a scalar in 3D, i.e., ∇q = p in vector calculus is equivalent to dq 0 = p 1 in exterior calculus. Here, the 1-form p 1 best represents the resulting vector p. The curl of a vector, ∇ × u = a is equivalent to du 1 = a 2 . The resulting vector a (in vector calculus) is represented as a 2-form (in exterior calculus). The divergence of a vector, ∇ · u = b is equivalent to du 2 = b 3 , resulting in to a 3-form which is a scalar. In 2D, gradient ∇, and 2D curl ∇× 2D are the only primary differentiation operations, and the divergence does not add any useful content. The gradient of a scalar looks the same as in 3D, i.e., it is equivalent to the exterior derivative operating on a 0-form, and implies the difference between two node values / discrete 0-forms. The curl of a vector lying in a plane , ∇ × 2D u = ∂uy ∂xx − ∂ux ∂xy , is equivalent to a 2-form resulting from the application of an exterior derivative to a 1-form (on an edge). The discrete exterior derivative d 0 operates on primal 0-forms and results in to primal 1-forms, and the discrete exterior derivative d 1 operates on primal 1-forms and results in to primal 2-forms. Similarly, (in 2D) the discrete exterior derivative −d T 0 (the negative sign is due to the adopted mesh orientation convention) operates on dual 0-forms and results in to dual 1-forms and the discrete exterior derivative d T 1 operates on dual 1-forms and results in to dual 2-forms. The Hodge star operator has a metric term incorporated in it. The Hodge star operating on a form in one dimension (e.g., k) results in to a form in the reflected dimension (N − k), i.e., * a k = b N −k , where N stands for the space dimension. In 3D, Hodge star on a 0-form (scalar value at a point / node) results in to a 3-form (cell average of scalar), and vice-versa. Applying Hodge star to a 1-form (vector component along a line) results in a 2-form (vector component normal to a face) and vice-versa. In practice, this implies that a discrete Hodge star operator transfers the forms on a mesh to a dual mesh (as well as changes their dimension). One need to identify which dual mesh is being used for the definition of the discrete Hodge star matrices, since every well formed mesh has an infinite number of well formed dual meshes. The median dual mesh (barycentric dual mesh), and the voronoi dual mesh (circumcentric dual mesh) are the examples of dual meshes for unstructured (triangular or tetrahedral) primary meshes. The discrete Hodge star ( * 0 ) operating on cell node (vertex) values (0-forms) of a tetrahedral mesh results in cell average values (dual 3-forms) for the surrounding Voronoi cells. Similarly, the inverse operator * −1 1 operates on 1-forms lying on a dual mesh edge (i.e., on the line between circumcenters of two adjacent primal cells), and results in to 2-forms on the primal mesh (normal flux values on the cell faces). In 2D, a discrete * 0 operates on mesh vertex values (primal 0-forms), and results in to Voronoi cell values (dual 2-forms). The operation * −1 0 transfers the cell average values (dual 2-forms) to mesh vertex values (primal 0-forms). The discrete * 1 operator transfers the primal 1-forms (primal edge values) to dual 1-forms (dual edge values). On the other hand, The discrete * −1 1 operator transfers the dual 1forms (dual edge values) to primal 1-forms (primal edge values). The wedge product of two forms implies multiplication of the forms. The dimension of the resultant form is the sum of the dimensions of the two forms in the wedge product, i.e., a j ∧ b k = c j+k for j + k ≤ N , with a property a j ∧ b k = (−1) jk b k ∧ a j . The wedge product is also associative, i.e., a j ∧ b k ∧ c l = a j ∧ b k ∧ c l for j + k + l ≤ N . In 3D, the wedge product implies a regular multiply, a dot-product, or a cross-product depending on the dimension of the two forms involved in the wedge product. The wedge product implies a regular multiplication,