A mass-, kinetic energy- and helicity-conserving mimetic dual-field discretization for three-dimensional incompressible Navier-Stokes equations, part I: Periodic domains

We introduce a mimetic dual-field discretization which conserves mass, kinetic energy and helicity for three-dimensional incompressible Navier-Stokes equations. The discretization makes use of a conservative dual-field mixed weak formulation where two evolution equations of velocity are employed and dual representations of the solution are sought for each variable. A temporal discretization, which staggers the evolution equations and handles the nonlinearity such that the resulting discrete algebraic systems are linear and decoupled, is constructed. The spatial discretization is mimetic in the sense that the finite dimensional function spaces form a discrete de Rham complex. Conservation of mass, kinetic energy and helicity in the absence of dissipative terms is proven at the discrete level. Proper dissipation rates of kinetic energy and helicity in the viscous case is also proven. Numerical tests supporting the method are provided.


Relevance of structure preserving methods with focus on kinetic energy and helicity conservation
In this work we address the discretization of the incompressible Navier-Stokes equations, defined on a periodic domain Ω ⊂ R 3 and time interval (0, t F ]. These well known equations govern the dynamics of an incompressible fluid's velocity, u : Ω × (0, t F ] → R 3 , and pressure, p : Ω × (0, t F ] → R, subject to a body force, f : Ω × (0, t F ] → R 3 , and an initial condition, u 0 : Ω → R 3 . A general dimensionless form of these equations is ∂u ∂t in Ω × (0, t F ] , (1b) where C(u) and D(u) represent the nonlinear convective term and the linear dissipative term, respectively, and Re is the Reynolds number. The operators C and D can take different forms, all analytically equivalent at the continuous level, see for example [1,2,3,4,5,6,7].
The four most common forms of the nonlinear convective term C(u) present in the literature, e.g., [3,5,8], are where ω := ∇ × u is the vorticity field. Besides these most common forms, it is also possible to construct a wide range of nonlinear convective terms as linear combinations of the above mentioned ones Preprint submitted to Journal of Computational Physics October 25, 2021 and/or employing vector calculus identities. For example, one such choice with interesting properties is the EMAC scheme [9]. Following similar ideas, it is possible to construct analytically equivalent representations for the dissipative term D, for example, where the latter representation can be derived from the former by using the identity ∆u = ∇ (∇ · u) − ∇ × (∇ × u) and the divergence free condition (1b). As mentioned before, these different forms are equivalent at the continuous level and, therefore, may be used interchangeably. At the discrete level, see for example [3,5,8], a particular choice of convective term used as the starting point of the discretization process leads to numerical schemes with substantially different properties.
One interesting aspect of the incompressible Navier-Stokes equations (1) is the fact that, in the inviscid limit (Re → ∞) and when the external body force is conservative (there exists a scalar field ϕ such that f = ∇ϕ), its dynamics conserves several invariants. Some of these invariants are the total kinetic energy K (in 2D and 3D), total enstrophy E (in 2D), and the total helicity H (in 3D), (4) K := 1 2 Ω u · u , E := 1 2 Ω ω · ω, and H := provided there is no net in-or out-flow of kinetic energy, enstrophy or helicity over the domain boundary. Note that, in 2D, vorticity can be regarded as a vector field constrained to the direction orthogonal to the planar 2D domain and velocity can be regarded as a vector field whose component along the direction orthogonal to the planar 2D domain is zero, i.e, ω = [0, 0, ω] T and u = [u, v, 0] T . Thus helicity is trivially zero in 2D flows.
The proofs for these conservation laws are straightforward. For illustration purposes and as an introduction to some of the ideas discussed later in this work, we present these proofs here for the case of no external force, i.e., f = 0, and periodic boundary condition. For simplicity, and without loss of generality, we use the rotational (or Lamb) form for the nonlinear convective term, (2d). The total (or Bernoulli) pressure is defined as Kinetic energy conservation (in 2D and 3D) corresponds to dK dt = 0. Differentiating K as defined in (4) with respect to time and taking (1) in the inviscid limit, Re → ∞, leads to where we have used (i) the vector calculus relation that the cross product of two vectors is perpendicular to either vector, i.e.,  Enstrophy conservation (in 2D) equates to dE dt = 0. As done above for kinetic energy, time differentiation of E as defined in (4) gives Computing the curl of the momentum equation in (1) with Re → ∞ and substituting ω = ∇ × u into (6) results in where we first used the vector calculus identity ∇ × (ω × u) = 1 2 (u · ∇ω) + 1 2 ∇ · (uω) , 2 followed by integration by parts on the first term of the second equality.
If we now use the momentum equation in (1) and its curl, (7) may be rewritten as where we have used (i) the definition of vorticity ω := ∇ × u, (ii) integration by parts on the second and fourth terms in the right side of the first identity, (iii) the vector calculus relation (5), and (iv) the identities ∇ × ∇ (·) ≡ 0 and ∇ · ∇ × (·) ≡ 0. These conservation laws for kinetic energy (in 2D and 3D), enstrophy (in 2D), and helicity (in 3D), are the expression of a more general structure underlying the incompressible Euler equations: the Hamiltonian structure, [10,11,12,13,14,15,16]. A system of partial differential equations (PDEs) is Hamiltonian if it can be cast in the general form, see for example [12,17], where S is a skew-adjoint operator, such that the induced bilinear form must also be a derivation and satisfy the Jacobi-identity, and H is the Hamiltonian functional. The system of equations (1), in the inviscid limit, is not in Hamiltonian form but may be rewritten in this form if pressure is eliminated and the Hamiltonian functional is set to the kinetic energy K. This can be achieved by either: (i) restricting the momentum equation to divergence free velocity fields (e.g., by making use of the stream function (in 2D) or the stream vector field (in 3D) ψ such that u = ∇ × ψ) [18], or (ii) taking the curl of the momentum equation (transforming it into the vorticity equation) [12].
Noether's theorem establishes a connection between conservation laws of a Hamiltonian system and its underlying symmetries, [19,20,21,22], thus highlighting the strong connection between the (geometric) structure of a system of PDEs and its dynamics. For example, spatial translation symmetry gives rise to conservation of linear momentum, and temporal translation symmetry results in energy conservation.
Helicity, on the other hand, is a more subtle quantity. Helicity as introduced in (4) is a particular case of the general concept of helicity of a divergence-free vector field, v, tangent to the boundary ∂Ω of a simply connected domain Ω ⊂ R 3 , see for example [20,23], which measures the average linking of its field lines.
The 19th century works of Helmholtz [24] and Kelvin [25] contain the seminal ideas for the modern concept of helicity, [26]. A renewed interest in these ideas appeared only later in the mid 20th century, first in the context of magnetohydrodynamics (MHD), [27], and then for hydrodynamics, [28,29]. Moreau, [28], discovered the law of conservation of helicity, and the term helicity appeared first in the work by Moffatt, [29], where the topological nature of this quantity was highlighted. For a detailed historical discussion of helicity see the very informative works by Moffatt [30,31].
It is possible to show, see for example [23], that helicity of any divergence-free vector field is preserved under the action of any volume preserving diffeomorphism. This property shows that helicity is not a dynamical invariant but a topological invariant, since its conservation is independent of the specific diffeomorphism, [20]. In fact, helicity is associated to the nontrivial kernel of the operator S, and is a Casimir for the Hamiltonian formulation of the inviscid Navier-Stokes equations, [20]. In the same way, in 2D, enstrophy is also a Casimir of the inviscid Navier-Stokes equations (as are all integral powers of vorticity). This very brief digression into Hamiltonian formalism intends to show the connection between the physical properties of a system of PDEs and its underlying geometrical structure. Invariants are not mere incidental features of the dynamics of a system, they are expressions of the underlying structure of the equations.
Helicity plays an important role in the generation and evolution of turbulence, [32,33,34]. The joint cascade of energy and helicity, [35], is an active field of research, [36,37,38,39,40]. Particularly important is the interaction between the two and how helicity impacts the energy cascade and, therefore, turbulence, [29,32,34,38,41,42,43]. This complex interaction between the energy cascade and the helicity cascade and especially the suppressive role of helicity motivates the focus on the development of discretization schemes that, besides conserving energy, conserve helicity. In the same way as energy conserving schemes have shown to substantially contribute to a higher fidelity in simulations, see for example [3,44,45,46,47,48,49], due to the connection between the cascades of energy and helicity, helicity conserving schemes should also present a positive impact towards improving the simulation accuracy.

Overview of structure-preserving methods for fluid flows
As highlighted above, the solutions to systems of PDEs (of which the Navier-Stokes are a particular example) satisfy strong constraints [50,51]. These constraints reflect the underlying mathematical structure of the equations (e.g., Hamiltonian structure, Poisson structure, de Rham sequence). These fundamental mathematical structures have long played an essential role in modern physics and pure mathematics. Owing to the fundamental nature of these structures and their impact on the dynamics of the systems under study, in recent years there has been an increasing interest in the various aspects of structure preservation at the discrete level [50,52,53]. This interest is rooted in three important points. First, there are well known connections between discrete structure preservation and standard properties of numerical methods [50,54,55]. Second, standard properties only guarantee physical fidelity in the limit of fully (at least highly) resolved discretizations. Reaching this limit requires infeasible computational resources (e.g., [56]). In contrast, structure preserving discretizations, by construction, generate solutions that satisfy the underlying physics even in highly under-resolved simulations. This is extremely relevant since most (if not all) simulations are inherently under-resolved. Third, physics preservation is fundamental when coupling systems in multiphysics problems [57].
The underlying principle behind structure preserving discretizations is to construct discrete approximations that retain as much as possible the structure of the original system of PDEs. A departure from this principle introduces spurious unphysical modes that pollute the physics of the system being modeled [54,58,59]. For example, as seen before, turbulence plays a fundamental role in the dynamics of the flow. A correct representation of the turbulent dynamics of a fluid is paramount in order to achieve accurate simulations. For this reason, if a numerical discretization introduces spurious unphysical energy dissipation into the system, it will fail to accurately capture the energy cascade and consequently the turbulent dynamics, [60,61,62]. The main focus of structure preserving discretizations for flow problems has been on energy conservation, e.g. [3,44,45,61]. As noted in the previous section, there is a growing knowledge on the role played by helicity and its impact on the energy cascade. For this reason, more recently, helicity conservation at the discrete level has been addressed in the literature, see for example [5,62,63].
Most standard structure preserving discretizations can be seen as variations of staggered grid methods which date back to the pioneering works of Harlow and Welch [64], and Arakawa and colleagues [65,66]. These methods employ a discretization that distributes the different physical quantities (pressure, velocity, vorticity, etc) at different locations in the mesh (vertices, faces, cell centres). It can be shown that, by doing so, important conservation properties can be maintained. Since then, much work has been produced and a rich variety of different flavours of structure preserving discretizations have been presented: finite differences/finite volumes [67,68,69,70,71], discrete exterior calculus (DEC) [72], finite element exterior calculus (FEEC) [55,73,74] and the works by the authors [8,75,76,77,78].
More recently, another approach develops a discretization of the physical field laws based on discrete variational principles. This approach has been used in the past to construct variational integrators for Lagrangian systems, e.g. [79,80]. These ideas have been extended to magneto-hydrodynamics [81,82,83], incompressible flows [84], and geophysical flow [85,86].

Objective
In this work, extending the initial ideas introduced for the 2D case, see [8], we combine (i) a particular choice for the formulation of the Navier-Stokes equations with (ii) a structure preserving discretization. Specifically, we will present two velocity evolution equations (dual-field) in a rotational form, discretized by the mimetic spectral element method (MSEM) [75,87,88].
This formulation attempts to address the dual character of the velocity field in the incompressible Navier-Stokes equations. This dual character implies that it is natural to look for a solution for the velocity field in H(div; Ω) ∩ H(curl; Ω). At the continuous level this is easily achievable, but that is not true at the discrete level since the space H(div; Ω) ∩ H(curl; Ω) is hard to discretize. The use of two velocity field evolution equations enables the representation of this dual character. It is shown that in this way the resulting discretization conserves mass, kinetic energy, and helicity in 3D.
The vorticity fields in the rotational form of the nonlinear convective term, see (2d), serve as a means of exchanging information between the two evolution equations. Additionally, this leads to a leap-frog like scheme that handles the nonlinear rotational term by staggering in time the velocity and vorticity such that the resulting discrete algebraic systems are linearized and decoupled.
Overall, the objective of this novel approach is the construction of a discretization which conserves mass, kinetic energy and helicity for the incompressible Navier-Stokes equations in the absence of dissipative terms and predicts the proper decay rate of kinetic energy and helicity based on the global enstrophy and an integral quantity of vorticity, respectively.

Outline of paper
The outline of the paper is as follows: In Section 2, we introduce a dual-field mixed weak formulation and prove that it preserves the desired conservation properties. In Section 3, a conservative staggered temporal discretization scheme is applied to the formulation, which is followed by a mimetic spatial discretization in Section 4. Numerical results that support the method are presented in Section 5. Finally, a summary is given and potential future work is listed in Section 6.

A mass-, kinetic energy-and helicity-conserving formulation
In this section, we propose a new conservative formulation for the Navier-Stokes equations in periodic domains. As we will only consider periodic domains in this paper, from now on, Ω represents a 3D periodic domain. The function spaces are the classic Hilbert spaces which form an exact complex, namely, the well-known de Rham (or Hilbert) complex [8,55,58,59]: This complex plays a fundamental role in the proofs and analysis of the presented work.

The rotational form of the incompressible Navier-Stokes equations
If in (1) we use the rotational (or Lamb) form for the nonlinear convective term, (2d), and use the representaion D(u) = −∇ × ω for the linear dissipative term, (3), we obtain the rotational form of the incompressible Navier-Stokes equations, We have proven that, in 3D and in the inviscid limit (Re → ∞), these equations preserve total kinetic energy and total helicity over time for the case of no external body force, f = 0, in Section 1.1. For non-zero conservative external body force, f = ∇ϕ = 0, we can include it by replacing the total pressure by an extended total pressure (10) P := P − ϕ .
All analysis and proofs remain valid. Without loss of generality, in this paper we will only use zero external body force for the analysis and proofs. 5 When the flow is viscous, Re < ∞, the viscosity dissipates kinetic energy of the incompressible Navier-Stokes equations at rate (11) dK while it dissipates or generates helicity at rate where ·, · Ω denotes the inner product, i.e., if a, b are vectors and c, d are scalars. The viscosity always dissipates kinetic energy because the total enstrophy cannot be negative, E ≥ 0, see the definition of the total enstrophy in (4). It either dissipates or generates helicity because the term ω, ∇ × ω Ω generally can be either positive or negative (or zero). This means the dissipation rate of helicity can be negative.

A conservative dual-field mixed weak formulation
We propose the following dual-field mixed weak formulation for the rotational form of the incompressible Navier-Stokes equations: Remark 1. In this formulation, the terms (ω i × u i ) · i are not known to be L 2 -integrable for the vector fields that belong to the infinite dimensional function spaces H(curl; Ω) (i = 1) and H(div; Ω) (i = 2). Showing this integrability requires proving additional regularity of the velocity and vorticity variables, which we currently are unable to do. However, in the finite dimensional case, the known regularity is sufficient, see Section 4. Thus, despite the potential mathematical issue, we still write this formulation above for its clear interpretation and to motivate the discrete scheme.
Therefore, in practice, ω 2 may be dropped from (13) if we replace it by ∇ × u 1 . We leave in ω 2 above to maintain the clearness of the formulation.

Properties of the formulation
We now show that the proposed dual-field formulation (13) conserves (i) the mass in terms of u 2 and, in the case of conservative external body force and zero viscosity, (ii) the kinetic energy in the formats and (iii) the helicity in the formats We will also analyze the dissipation rate of kinetic energy and helicity in the viscous case for the proposed formulation.
Note that, in this subsection, everything is still at the continuous level. The purpose is to show that the proposed weak formulation possesses the same properties as the strong formulation does.

Mass conservation
For the mass conservation, since we have restricted u 2 to space H(div; Ω), the de Rham complex (8) and the constraint (13f) ensure that the relation is strongly satisfied; no integration by parts is required. Therefore, the mass conservation is satisfied for velocity u 2 . Such an approach is widely used to construct mass conserving discretizations. While for u 1 ∈ H(curl; Ω), the mass conservation is only weakly satisfied, see (13c).

Time rate of change of kinetic energy
In the inviscid limit (Re → ∞) and when f = 0, the kinetic energy conservation is equivalent to Because (13a) is valid for all 1 ∈ H(curl; Ω), we can select 1 to be u 1 ∈ H(curl; Ω). As a result, we get The second term vanishes because of (5). Meanwhile, from (13c), we know that u 1 , ∇P 0 Ω = 0 because P 0 ∈ H 1 (Ω). Therefore, the third term also vanishes, which accomplishes the proof of kinetic energy conservation for K 1 . Similarly, by selecting 2 of (13d) to be u 2 , we can get where the second and third terms vanish because (5) and (13f), respectively. Thus we can conclude that K 2 is also preserved over time.
In the viscous case, Re < ∞, if we repeat the above analysis, the viscous terms will remain. We will eventually obtain the following kinetic energy dissipation rates, where the total enstrophy E 1 and E 2 are defined as This is in agreement with the kinetic energy dissipation rate of the strong formulation, see (11).
We now reuse (13e) and select 1 to be u 1 ∈ H(curl; Ω). As a result, we get Thus both H 1 and H 2 are preserved over time.
In the viscous case, Re < ∞, if we repeat above analysis, the viscous contribution will not cancel and we will obtain the following helicity dissipation rate, which is consistent with that of the strong formulation, see (12).

Temporal discretization
Inspired by a mass, energy, enstrophy and vorticity conserving (MEEVC) [8] scheme for the 2D incompressible Navier-Stokes equations, we construct a staggered temporal discretization for the two evolution equations in the dual-field formulation (13). The MEEVC scheme, as well as the presented method, starts with a formulation of two evolution equations. The two evolution equations are discretized temporally at two sequences of time steps respectively using a Gauss integrator. The two sequences of time steps are staggered such that the endpoints of time steps in one sequence are exactly the midpoints of time steps in the other sequence. Thus at each time step either discrete evolution equation can use the solution from the other one as known variable at the midpoint, see Fig. 1.
We use a lowest order Gauss integrator as the time integrator [89,90,91]. For example, if we apply the integrator to an ordinary differential equation (ODE) of the form at a time step from time instant t k−1 to time instant t k , we obtain . Additionally, we will use the midpoint rule, namely, We further introduce two time sequences, the integer time steps and the half-integer time steps. The integer time steps use time instants indicated with integer superscripts. For example, kth (k = 1, 2, · · · ) integer time step (denoted by S k ) is from t k−1 to t k . The half-integer time steps use time instants indicated with half-integer superscripts. For example, kth (k = 1, 2, · · · ) half-integer time step (denoted byŜ k ) is from t k− 1 2 to t k+ 1 2 . These time steps satisfy In other words, we restrict ourselves to constant time intervals equal for both time sequences.

Temporal discretizations at staggered time steps
We now apply the time integrator (24) to evolution equations (13d) and (13a) at integer and halfinteger time steps, respectively.

Temporal discretization at integer time steps
If we apply the time integrator (24) to the evolution equation for u 2 (13d) at integer time steps, with the midpoint rule, see (25), and constraints (13e) and (13f), we can obtain a semi-discrete weak formulation at, for example, kth integer time step S k : Given ω k−1 where ω k− 1 2 2 is borrowed from the other time sequence, in particular, is the solution of ω 2 at (k − 1)th half-integer time step and, therefore, is known.

Temporal discretization at half-integer time steps
Similarly, we apply the time integrator (24) to the evolution equation for u 1 (13a) at half-integer time steps. With the midpoint rule, see (25), and constraints (13b) and (13c), we can get a second semidiscrete weak formulation at, for example, kth half-integer time stepŜ k : Given u where ω k 1 is borrowed from the other time sequence and, more specifically, is the solution of ω 1 at kth integer time step, see (26). Thus it is known. The solution ω can be sequentially used for the next, the (k + 1)st, integer time step. Thus iterations can proceed.

Overall temporal discretization
One may notice that to start the iterations we need to know u 2 . The simplest approach for the 0th time step is applying the explicit Euler method to evolution equation (13a) which, together with constraints (13b) and (13c) at t 1 2 , leads to a semi-discrete system similar to (27). More accurate approaches, like directly applying the Gauss integrator (24) or other (higher order) integrators to formulation (13), could also be used. These methods will eventually lead to nonlinear discrete algebraic systems for which more expensive iterative methods like the Newton-Raphson method are needed. After the 0th time step,ŝ 0 , standard iterations, S andŜ, can proceed. 1 The overall temporal scheme is illustrated in Fig. 1.
The kth integer time step also computes P It is easy to see that, instead of applying a standard temporal discretization directly to the dual-field mixed weak formulation (13), using the presented staggered temporal discretization can greatly reduce the computational cost. Although the dual-field formulation doubles the variables, we will only solve for half of them at each time step as the staggered temporal discretization decouples the dual-field formulation. Meanwhile, since each semi-discrete formulation borrows the solution from the other for the nonlinear terms, see the second terms of (26a) and (27a), the semi-discrete formulations will lead to linearized discrete algebraic systems.

Properties after temporal discretization
In this part, we check whether the conservation (in the inviscid case) and dissipation (in the viscous case) properties proven at the continuous level, see Section 2.3, are preserved after the proposed staggered temporal discretization. Note that we have not yet applied a spatial discretization; the function spaces are still the infinite dimensional Hilbert spaces in the de Rham complex, (8).

Mass conservation after temporal discretization
The mass conservation is not influenced by the temporal discretization. Due to the same proof as given in Section 2.3.1, for u k 2 ∈ H(div; Ω), ∇ · u k 2 = 0 is exactly satisfied at all integer time instants, and for u k+ 1 2 1 ∈ H(curl; Ω) the mass conservation is only weakly imposed through integration by parts, see (27c). , u

Time rate of change of kinetic energy after temporal discretization
Thus the kinetic energy is preserved at both integer and half-integer time steps. If Re < ∞, with the same analysis, we will obtain , ∇ × u With (26b), (27b), we can conclude that This shows that the boundedness of the kinetic energy (15) and (16) is preserved by this staggered temporal discretization, which can contribute to the stability of the scheme.

Time rate of change of helicity after temporal discretization
Let Re → ∞ and f = 0. We select 1 in (27a) to be ω k 1 and perform the same process for proving (17). We will get (32) u Analogously, by repeating the proof for (23), we can obtain Equations (32) and (33) together imply At the half-integer time stepŜ k−1 (assume k ≥ 2), (32) reads With this relation, we can extend (34) to If we further apply the midpoint rule, (25), to the first two terms and the last two terms of above equation, we obtain In addition, since (26b) holds for all 1 ∈ H(curl; Ω), we can fill u k Again, as ω 2 is only solved at half-integer time instants, see (27), we use the midpoint rule, (25), to bring it to the integer time instants, namely, . As a result, (37) implies And with (36), we can finally conclude that In the viscous case, Re < ∞, repeating above analysis at the half-integer time stepŜ k and at the integer time step S k (see (32) and (33)) leads to (39) u If we combine the above two equations, i.e., (39) + (40), and use the midpoint rule, (25), we obtain Again, (39) is still valid at (k − 1)st half-integer time step,Ŝ k−1 , (assume k ≥ 2) where it reads If we now combine (40) and (42), and use the midpoint rule, (25), we get We now can combine (41) and (43) and obtain Finally, because (38) still holds, the viscosity dissipates H k 1 and H k 2 at the same rate, denoted by

Mimetic spatial discretization
It has been shown that the de Rham complex plays an essential role in the proofs and analysis of the conservation properties and the dissipation rates for the proposed dual-field formulation at both continuous and semi-discrete levels. For example, (19) is valid because we have chosen P 0 ∈ H 1 (Ω) such that ∇P 0 ∈ H(curl; Ω) is guaranteed. Choosing u 1 ∈ H(curl; Ω) and ω 2 ∈ H(div; Ω) ensures that the relation ω 2 = ∇ × u 1 is satisfied exactly. In addition, as shown in Section 2.3.1, the de Rham complex is essential for the mass conservation ∇ · u 2 = 0 where u 2 ∈ H(div; Ω).
In this work we consider a set of discrete function spaces, i.e., they constitute a discrete de Rham complex. In order to enable the validity of the proofs and analysis at the fully discrete level, we need to employ such a set of discrete spaces for the spatial discretization. Any sequence of discrete function spaces that satisfies (45)  of degree N , see [92], RT N are the Raviart-Thomas spaces of degree N , see [92,93], and DG N −1 are the discontinuous Lagrange spaces of degree (N − 1). Another possible exact sequence of discrete function spaces employing b-splines is employed in the works by Hiemstra et al. [94], Buffa et al. [95], and Ratnani and Sonnendrücker [96]. We call these spaces structure-preserving or mimetic spaces, see another example, the mimetic polynomial spaces [8,75,87,88,97]. Note that variables in the finite dimensional spaces C(Ω) and D(Ω) possess the regularity that ensures the L 2 -integrability of the convective terms in the weak formulation (13), see Remark 1.

Fully discrete systems
Applying a particular set of mimetic spaces to the semi-discrete problems (26) and (27) leads to two local fully discrete linear algebraic systems, one for the kth integer time step S k , i.e., and one for the kth half-integer time stepŜ k , namely, If we rearrange the systems (46) and (47) and write them in linear algebra format, we can obtain following linear systems,    14 A similar spatial discretization can be applied to the semi-discrete system for the 0th time stepŝ 0 , see Fig. 1. Suppose a mesh has been generated in the computational domain Ω. We can perform such discretizations in all elements. After applying the initial condition and assembling the local systems, we will eventually obtain global linear systems ready to be solved in the sequence shown in Fig. 1.

Properties of the fully discrete systems
Since we have used a sequence of function spaces which form a discrete de Rham complex, the proofs for the conservation properties and the analysis for the dissipation rates of kinetic energy and helicity at the semi-discrete level, see Section 3.2, remain valid at the fully discrete level.

Numerical experiments
We now test the proposed mimetic dual-field method with two manufactured solutions and a more general flow, the well-known Taylor-Green vortex.
For all tests, we use the mimetic polynomial spaces as our mimetic spaces and do the spatial discretization under the framework of the MSEM. Meshes are uniform orthogonal structured hexahedral meshes. The mesh size, namely, the edge length of the cubic element cell, is denoted by h. The degree of the mimetic polynomials is denoted by N . And we use the explicit Euler method for the temporal discretization of the 0th time step, i.e.,ŝ 0 in Fig. 1. The implementation is conducted in Python.

Manufactured solution tests
Two manufactured solutions are taken from [62]; one for testing the conservation properties and one for investigating the convergence rate of the method. The domain is selected to be the periodic unit cube Ω := [0, 1] 3 .
Such an initial condition possesses kinetic energy K| t=0 = 0.75 and helicity H| t=0 = −6.283. The problem is solved until t = 10 on an extremely coarse mesh of h = 1/3 and N = 2.
We first try to verify that the proposed method does preserve mass, kinetic energy and helicity if in the inviscid limit Re → ∞ and f = 0. In Fig. 2 some results are presented. The results of ∇ · u k 2 L ∞ in the bottom-right diagram imply that the pointwise mass conservation is always satisfied. In the bottom-left diagram, the results show that both H k 1 and H k 2 are preserved. The fact that the two lines coincide with each other up to O(10 −10 ) verifies (38). As for kinetic energy, the results are present in the top diagrams where the discrete conservation for both K k− 1 2 1 and K k 2 at their corresponding time steps are shown. We then keep f = 0 and use a Re < ∞; we let the viscosity dissipate kinetic energy and helicity. Some results for Re = 100 are presented in Fig. 3 where the results shown in the top diagrams verify the dissipation rate of kinetic energy derived in (30) and (31) and the results in the bottom-left diagram are in agreement with the dissipation rate of helicity, see (44). The pointwise conservation of mass is still satisfied at all time steps as shown in the bottom-right diagram of Fig. 3.
In Fig. 4, some results of the magnitude of ∇ · u . It is not surprising that the error is large especially for the inviscid case as we have used an extremely coarse mesh. This is consistent with the analysis that the constraint of mass conservation is only weakly imposed for u Note that these tests are also valid for the non-zero conservative external body force. If ϕ is known and f = ∇ϕ = 0, we can still first conduct the test with f = 0 and get the same results. The only difference is that we now obtain the solution for the extended total pressure P , see (10). We can post-process P with the known ϕ to retrieve the solution for total pressure P .

Convergence tests
We now investigate whether the proposed method produces converging solutions and, if yes, what is the convergence rate of the proposed method with a manufactured solution. Assume solve the Navier-Stokes equations with Reynolds number Re = 1 and the body force f which can be calculated from u, p and Re using the Navier-Stokes equations. The exact solutions of vorticity ω and total pressure P can also be calculated. We use u| t=0 as initial condition and let the flow evolve for different mesh element sizes and polynomial space degrees. Errors are then measured at t = 2.
Results are presented in Fig. 5 where the optimal convergence rates are observed for all variables of the dual-field formulation when the mesh is h-refined under different polynomial degrees. The plot of ∇ · u h 2 L ∞ shows that the pointwise conservation of mass is satisfied up to the machine precision in all cases.
We can also measure the difference between the two dual solutions of one physical variable. The results of u h 2 − u h 1 L 2 and ω h 2 − ω h 1 L 2 at t = 2 are shown in Fig. 6. Note that, since u h 1 and u h 2 (ω h 1 and ω h 2 ) are staggered in time, we have used the midpoint rule, (25), to u h 1 (ω h 1 ) such that it can compared to u h 2 at t = 2, an integer time instant. It is not surprising that they converge under p-or h-refinement. This suggests that we can use them for accuracy indicators, for example, which can be very helpful for general (non-manufactured) simulations. More interestingly, one can measure the local difference of the dual solutions and use it as an indicator for mesh adaptivity, which is outside of the scope of the current paper. From this aspect, the existence of dual representations of the solution for one variable can be regarded as an advantage for the proposed method. Despite the existence of the difference between the dual representations, both of them should be considered as equally important solutions of the variable. Recall the dual character of the velocity field which is hard to capture in one discrete space, see Section 1.3. The dual representations together can be regarded as a discretization of its dual character.

Taylor-Green vortex
We now test the method with a more general flow, the Taylor-Green vortex (TGV) flow. The domain is given as Ω := [−π, π] 3 and is periodic. V = 8π 3 denotes the volume of the domain. The body force is  Such an initial condition possesses kinetic energy K| t=0 = 0.125 and zero helicity. We solve the flow using the proposed mimetic dual-field method at Re = 500. Iso-surfaces of ω x 1 = −3 ω h 1 = (ω x 1 , ω y 1 , ω z 1 ) at some time instances are shown in Fig. 7. It is seen that the flow initially induces vortices of clear structures which then break down and finally are dissipated by the viscosity.
In Fig. 8 and Fig. 9, results of total kinetic energy and total enstrophy are presented. These results are compared to benchmarks taken from [98]. In Fig. 8 we can see that the proposed mimetic dual-field method, compared to a discontinuous Galerkin (DG) method of the same order (p = 2) and in the same mesh (32 × 32 × 32 elements), produces better results in terms of the error to the results produced by a reference, a very high (128th) order spectral method. This is mostly clear in the enstrophy results near t = 9 when the total enstrophy reaches its peak; the DG method is not able to capture the peak of the total enstrophy while the mimetic dual-field method captures it well for both of the dual solutions. Similar comparisons are made for more resolved simulations in Fig. 9, where improved results are seen especially near the peak of total enstrophy; the DG method now is able to capture the peak and the mimetic dual-field method captures the shape of the peak better.
In Fig. 10, some results of the total helicity versus time for the TGV flow is shown. It is seen that, as the flow evolves, the total helicity remains zero (to the machine precision). Such a phenomenon is consistent with the fact that the dissipation rate of helicity, see (44), is constantly zero (to the machine precision) as shown in the same diagram.
In Fig. 11, the results of kinetic energy spectra at t = 9.1 are presented. In the left diagram, it is seen that, in terms of kinetic energy, the mimetic dual-field method has similar accuracy as the DG method for large scales (k ≤ 10). For medium scales (10 < k ≤ 35), both methods start to deviate from the high order spectral reference results with the proposed dual-field method showing less overdissipation. For small scales (k > 35), both methods show large deviations from the reference results. The interesting aspect is that, for small scales, the DG method and the proposed dual-field method present different behaviors: the DG method over dissipates the energy and the dual-field method accumulates energy. The accumulation of energy at small scales is expected due to the energy conservation properties of the dual-field method. The energy cascade occurs up to the resolved scales and then it is stored (and accumulates at the smaller scales). It is the authors opinion that this can be an advantage of this method since subscale grid methods can specifically target these small scales and introduce the required dissipation that is not resolved. In opposition, the DG method already over dissipates the energy, therefore it is challenging for a dissipation based sub-grid scale model to improve the results for these smaller scales. This is a topic of interest for the authors and will be further researched in the future. A partial support for this claim is the results presented in the right diagram of Fig. 11 where it is seen that the value of k where energy accumulation starts decreases when a less resolved discretization is employed.   Figure 8: A comparison of kinetic energy and enstrophy results of the TGV test. The reference and the DG-32p2 results are taken from [98]. The reference method is a 128th order spectral method. DG stands for a discontinuous Galerkin method. MDF stands for the mimetic dual-field method. 32p2 represents h = 1/32 and the degree of the polynomial spaces is 2. The time step interval is selected to be ∆t = 1/50 for MDF-32p2. MDF-24p3, E k 1 Figure 9: A comparison of kinetic energy and enstrophy results of the TGV test. The reference and the DG-24p3 results are taken from [98]. The reference method is a 128th order spectral method. DG stands for a discontinuous Galerkin method. MDF stands for the mimetic dual-field method. 24p3 represents h = 1/24 and the degree of the polynomial spaces is 3. The time step interval is selected to be ∆t = 1/50 for MDF-24p3.   Figure 11: Kinetic energy spectra of the TGV test. The reference and the DG-32p2 results are taken from [98]. The reference method is a 128th order spectral method. DG stands for a discontinuous Galerkin method. MDF stands for the mimetic dual-field method. 32p2, 24p2 and 16p2 represents that the mesh size h = 1/32, 1/24 and 1/16, respectively, and the degree of the polynomial spaces is 2. The time step interval is selected to be ∆t = 1/50, 1/40 and 1/30 for MDF-32p2, MDF-24p2 and MDF-16p2, respectively.
6. Summary and future work

Summary
In this paper, we introduce a discretization which satisfies pointwise mass conservation and, if in the absence of dissipative terms, conserves total kinetic energy and total helicity and, otherwise, properly captures the dissipation rates of total kinetic energy and total helicity for the 3D incompressible Navier-Stokes equations. The discretization is based on a novel dual-field mixed weak formulation where two evolution equations are employed. A staggered temporal discretization linearizes the convective terms and reduces the size of the discrete systems, which can be regarded as a big advantage of the proposed method in terms of the computational efficiency. A mimetic spatial discretization enables the validity of the conservation properties and the dissipation rates at the fully discrete level.

Future work
In this paper, u 1 and u 2 (ω 1 and ω 2 ) are approximated in different function spaces, and, therefore, equality between them will not hold unless the flow is fully resolved. In other words, we are not able to construct a square, time-independent and explicit discrete Hodge operator. Instead, this method implicitly defines a time-dependent discrete Hodge operator. By allowing the time evolution of the discrete Hodge operator we can construct a helicity conserving scheme. Based on the promising results reported in this paper, we want to apply the algebraic dual polynomial spaces, [99], such that solutions u h 1 22 and u h 2 (ω h 1 and ω h 2 ) are two representations in a pair of algebraic dual polynomial spaces. As a result, we expect the difference between u h 1 and u h 2 (ω h 1 and ω h 2 ) to be smaller and using the vorticity from the other subset of equations, see (26) and (27), to be more consistent.
In the kinetic energy spectra of the dual field formulation, Fig. 11, we see that for high wave numbers the energy decay is insufficient. This is attributed to the fact that the scheme is non-dissipative and the grids are too coarse for energy at the small scales to dissipate. In future work we want to add a sub-grid scale model on the momentum equations for u 1 and u 2 of the form ∆(u 1 − u 2 ) to the u 1 equation, (13a) and ∆(u 2 − u 1 ) to the u 2 equation, (13d). If we definē u = 1 2 (u 1 + u 2 ) , this sub-grid scale diffusion cancels from the average, while it only acts on the difference between the two fields u = 1 2 (u 1 − u 2 ) .
So the numerical dissipation only acts on the difference of the two fields. This implies that the added diffusion is only active for the large wave numbers, where the difference between u 1 and u 2 is significant, while for the small wave numbers where u 1 and u 2 are almost the same, no dissipation takes place of u . In this sense the dual field formulation can be used as a turbulence model. Future work needs to establish how the parameter should be chosen. Other steps we want to report in the future includes, for example, error analysis, mesh adaptivity based on the local difference between u h 1 and u h 2 (ω h 1 and ω h 2 ) and the extension from periodic boundary conditions to general boundary conditions.