StructurePreserving Discretization of Incompressible Fluids
Abstract.
The geometric nature of Euler fluids has been clearly identified and extensively studied over the years, culminating with Lagrangian and Hamiltonian descriptions of fluid dynamics where the configuration space is defined as the volumepreserving diffeomorphisms, and Kelvin’s circulation theorem is viewed as a consequence of Noether’s theorem associated with the particle relabeling symmetry of fluid mechanics. However computational approaches to fluid mechanics have been largely derived from a numericalanalytic point of view, and are rarely designed with structure preservation in mind, and often suffer from spurious numerical artifacts such as energy and circulation drift. In contrast, this paper geometrically derives discrete equations of motion for fluid dynamics from first principles in a purely Eulerian form. Our approach approximates the group of volumepreserving diffeomorphisms using a finite dimensional Lie group, and associated discrete Euler equations are derived from a variational principle with nonholonomic constraints. The resulting discrete equations of motion yield a structurepreserving time integrator with good longterm energy behavior and for which an exact discrete Kelvin’s circulation theorem holds.
1. Introduction
The geometric nature of Euler fluids has been extensively studied in the literature in works of Arnold, EbinMarsden and others; however the geometricdifferential standpoint of these studies sharply contrasts with the numerical approaches traditionally used in Computational Fluid Dynamics (CFD). In particular, methods based on particles, vortex particles, staggered Eulerian grids, spectral elements, as well as hybrid LagrangianEulerian formulations were not designed with structure preservation in mind — in fact, recent work pinpoints the loss of Lagrangian structures as a major numerical impediment of current CFD techniques [27]. In contrast, structure preserving methods (socalled geometric integrators) have recently become popular in the context of Lagrangian dynamics in solid mechanics. Based on discrete versions of Hamilton’s principle and its variants, they have been shown to capture the dynamics of the mechanical system they discretize without traditional numerical artifacts such as loss of energy or momenta.
While the variational principles for incompressible fluid mechanics are best expressed in a Lagrangian formalism, computational efficiency often calls for an Eulerian treatment of fluid computations to avoid numerical issues inherent to deforming meshes. In order to circumvent these issues without giving up structure preservation, a new Eulerian formulation of discrete fluid mechanics is thus needed.
Guided by the variational integrators used in the Lagrangian setting, this paper introduces a discrete, structurepreserving theory for incompressible perfect fluids based on Hamiltond’Alembert’s principle. Such a discrete variational approach to fluid dynamics guarantees invariance under the particlerelabeling group action and gives rise to a discrete form of Kelvin’s circulation theorem. Due to their variational character, the resulting numerical schemes also exhibit good longterm energy behavior. In addition, the resulting schemes are not difficult to implement in practice (see Figure 1), and we will derive particular instances of numerical update rules and provide numerical results. We will favor formalism over smoothness in the exposition of our approach in order to better elucidate the correspondences between continuous and discrete expressions.
1.1. Brief Review of the Continuous Case.
Let be an arbitrary compact manifold, possibly with boundary (where denotes the dimension of the domain, typically, 2 or 3), and be the group of smooth volumepreserving diffeomorphisms on . As was shown in [2], the motion of an ideal incompressible fluid in may be described by a geodesic curve in . That is, serves as the configuration space—a particle located at a point at time travels to at time . Being geodesics, the equations of motion naturally derive from Hamilton’s stationary action principle:
(1) 
subject to arbitrary variations vanishing at the endpoints. Here, the Lagrangian is the kinetic energy of the fluid and is the standard volume element on . As this Lagrangian is invariant under particle relabeling—that is, the action of on itself by composition on the right, the principle stated in Eq. (1) can be rewritten in reduced (Eulerian) form in terms of the Eulerian velocity :
(2) 
subject to constrained variations (called Lin constraints), where is an arbitrary divergencefree vector field—an element of the Lie algebra of the group of volume preserving diffeomorphisms—and is the Jacobi Lie bracket (or vector field commutator). There is a complex history behind this reduced variational principle which was first shown for general Lie groups by [34]; see also [2, 22, 5, 32]). As stated above, the reduced Eulerian principle is more attractive in computations because it involves a fixed Eulerian domain (mesh); however, the constrained variations necessary in this context complicates the design of a variational Eulerian algorithm.
1.2. Overview and Contributions.
While time integrators for fluid mechanics are often derived by approximating equations of motion, we instead follow the geometric principles described above and discretize the configuration space of incompressible fluids in order to derive the equations of motion through the principle of stationary action. Our approach uses an Eulerian, finite dimensional representation of volumepreserving diffeomorphisms that encodes the displacement of a fluid from its initial configuration using special orthogonal, signed stochastic matrices. From this particular discretization of the configuration space, which forms a finite dimensional Lie group, one can derive a rightinvariant discrete equivalent to the Eulerian velocity through its Lie algebra, i.e., through antisymmetric matrices whose columns sum to zero. After imposing nonholonomic constraints on the velocity field to allow transfer only between neighboring cells during each time update, we apply the Lagranged’Alembert principle (a variant of Hamilton’s principle applicable to nonholonomic systems) to obtain the discrete equations of motion for our fluid representation. As we will demonstrate, the resulting Eulerian variational Liegroup integrator is structurepreserving, and as such, has numerous numerical properties, from momentum preservation (through a discrete Noether theorem) to good longterm energy behavior.
1.3. Notations.
The spatial discretization (mesh), either simplicial (tetrahedra) or regular (cubes), will be denoted , with being the number of dimensional cells in . The size of a mesh will refer to the maximum diameter of its cells. The Lebesgue measure will be denoted by . Thus, is the volume of cell , is the area of the face common to and , etc (see Figure 2). The dual of is the circumentric dual cell complex [28], formed by connecting the circumcenters of each cell based on the connectivity of . We will further assume that the mesh is Delaunay with well shaped elements [47] to avoid degeneracies of its orthogonal dual as well as to simplify the exposition. We will also use the term regular grid (or Cartesian grid) to designate a mesh that consists of cells that are dimensional cubes of equal size. The notation will denote the set of indices of cells neighboring cell , that is, cell shares a face with cell iff . We will say that a pair of cells is positively oriented around an edge if they share a face containing and they are oriented such that they “turn” clockwise around the edge when viewed along the oriented edge. The same term will be used similarly for triplets of cells where and all three cells contain edge .
The notation and will respectively refer to the inner product of vectors and the pairing of oneforms and vector fields, while their discrete counterparts will be denoted by and . Table 1 summarizes the main variables used in the remainder of this paper, along with their meaning and representation.
Symbol  Meaning  Representation 

Domain of motion  
Dimension of the domain  
Configuration space of ideal fluid  Volumepreserving diffeomorphisms on  
Tangent space of at Id  Divergencefree vector fields on  
Mesh discretizing domain  Simplicial or regular mesh  
Number of cells in  
Cell of  Tetrahedron or cube in  
Discrete analog of volume form  Diagonal matrix of cell volumes,  
Discrete configuration space  orthogonal signed stochastic matrices  
Lie algebra of  antisymmetric nullrow matrices  
Discrete configuration  Matrix  
Discrete Eulerian velocity  Matrix  
Discrete form  dimensional tensor of order  
Space of matrices with sparsity based on cell adjacency  Constrained set of matrices, with  
Space of sparse discrete velocities  Constrained set of velocities, 
Acknowledgments. We thank Daryl Holm and Yann Brénier for helpful early discussions and input, Evan Gawlik for generating the energy plots, and Keenan Crane for generating our 2D tests. This research was partially supported by NSF grants CMMI0757106, CCF0811373, and DMS0453145.
2. Discrete Volume Preserving Diffeomorphisms
We first introduce a finite dimensional approximation to the infinite dimensional Lie group of volume preserving diffeomorphisms that tracks the amount of fluid transfered from one cell to another while preserving two key properties: volume and mass preservation.
2.1. Finite Dimensional Configuration Space.
Suppose that the domain is approximated by a mesh . Our first step in constructing a discrete representation of ideal fluids is to approximate with a finite dimensional Lie group in such a way that the elements of the corresponding Lie algebra can be considered as a discretization of divergencefree vector fields. To achieve this goal, we will not discretize the diffeomorphism itself, but rather the associated operator defined by . Here is the space of square integrable real valued functions on . An important property of is given by the following lemma, which follows from the change of variables formula.
Lemma 1 (Koopman’s lemma^{1} ).
If the diffeomorphism is volumepreserving, then is a unitary operator on .
Another important property of is that it preserves constants, i.e., for every constant function , which can be seen as mass preservation for fluids. Next we present an approach to discretize this operator while respecting its two defining properties.
Discrete Functions. To discretize the operator we first need to discretize the space on which acts. Since the mesh splits the domain of motion into cells of maximum diameter , a function can be approximated by a step function , constant within each cell of the mesh, through a map , which averages per cell:
where is the indicator function for the cell , and is the volume of cell . Since the space of all step functions on is isomorphic to , we can consider the step functions as vectors: using the map defined by
(3) 
we can define a vector of size to represent the step function . To reconstruct a step function from an arbitrary vector we define an operator by
Thus, the operators , and are related through:
The vector will be called a discrete function as it provides an approximation of a continuous function : when ,
We also introduce a discrete approximation of the continuous inner product of functions through:
(4) 
Discrete Diffeomorphisms. Using the fact that a matrix (here is the space of real valued matrices) acts on a vector , we will say that approximates if is close to :
Definition 1.
Consider a family of meshes of size , each consisting of cells . We will say that a family of matrices approximates a diffeomorphism (and denote this property as: ) if the following is true:
In order to better respect the continuous structures at play, we further enforce that our discrete configuration space of diffeomorphisms satisfies two key properties of : volumepreservation, reflecting the fact that is unitary, and total mass preservation, as preserves constants. We will thus only consider matrices that:

preserve the discrete inner product of functions, i.e.,
where the inner product of discrete functions is defined by Eq. (4). Denoting
note that this discrete notion of volume preservation directly implies that for our mesh a volume preserving matrix satisfies
The matrix is thus orthogonal, restricted to matrices of determinant .

preserve constant vectors (i.e., vectors having all coordinates equal) as well:
The matrix must thus be signed stochastic as well.
Consequently, the finite dimensional space of matrices we will use to discretize volumepreserving diffeomorphisms has the following definition:
Definition 2.
Let be a mesh consisting of cells , and be the diagonal matrix consisting of volumes of the cells, i.e., and when (we will abusively use the shorter notation to denote a diagonal element of for simplicity in what follows). We will call a matrix volumepreserving and constantpreserving with respect to the mesh if, for all in
(5) 
and
(6) 
The set of all such orthogonal, signed stochastic matrices of determinant will be denoted , and will be used as a discretization of the configuration space .
Our finite dimensional configuration space for fluid dynamics is thus the intersection of two Lie groups: the orthogonal group, and the group of invertible stochastic matrices; therefore, it is a Lie group. Note that if all cells of have the same volume, i.e., , then a matrix is orthogonal in the usual sense and the equality (5) implies . For such meshes (which include Cartesian grids), the matrix is signed doublystochastic.
Remark. An alternate, arguably more intuitive way to discretize a diffeomorphism on a mesh would be to define a matrix as:
This discretization also satisfies by definition a discrete preservation of mass and a (different) notion of volume preservation. While it has the added benefit of enforcing that has no negative terms (therefore respecting the positivity of ), the class of matrices it generates is, unfortunately, only a semigroup, which would be an impediment for establishing a variational treatment of fluids as an inverse map will be needed in the Eulerian formulation. So instead, we take the orthogonal part of this matrix as our configuration (which can be obtained in practice through the polar decomposition). Notice that polar factorization has often been proposed in the context of fluids (see, e.g., [10]), albeit for more general nonlinear Hodgelike decomposition.
2.2. Discrete Velocity Field.
Now that we have established a finite dimensional configuration space , we describe its associated Lie algebra, and show that elements of this Lie algebra provide a discretization of divergencefree vector fields . We will assume continuous time for simplicity, but a fully discrete treatment of space and time will be introduced in Section 5.
Consider a smooth path in the space of volumepreserving diffeomorphisms with , and let be an approximation of , i.e., for any piecewise constant function approximating a smooth function a discrete version of is given by
Assuming is smooth in time, we define its Eulerian velocity to be
thus yielding
Since , where and is the Lie derivative, the matrix represents an approximation of the Eulerian velocity field , which motivates the following definition:
Definition 3.
Consider a oneparameter family of volumepreserving diffeomorphisms and the associated timedependent vector field . Consider a family of meshes of size consisting of cells and an operator defined by Eq. (3).
We will say that a family of matrices approximates a vector field (denoted by ) if the following statement is true:
Remark. The choice of the minus sign in the definition of stems from the fact that represents (thus, in essence). Since , we picked the sign to make represents , consistent with the continuous case.
If a curve of matrices belongs to the configuration space (i.e., if is orthogonal signed stochastic), then its associated belongs to its Lie algebra that we denote as . Matrices from this Lie algebra inherit the properties that their rows must sum to zero:
and they are antisymmetric:
These two properties can be intuitively understood as discrete statements that represents an advection, and the vector field representing this advection is divergencefree. Lie algebra elements for arbitrary simplicial meshes will be called nullrow antisymmetric matrices. Note that if the mesh is regular (), belongs to the orthogonal group and the matrix has to be antisymmetric with both its rows and columns summing to zero (“doubly null”).
The link between convergence of to and convergence of to is described by the following lemma.
Lemma 2.
Consider the setup of Definition 3 and suppose a family of matrices approximates the Lie derivative (in the sense of Definition 3) uniformly in when for some .
Then there is a family of matrices such that and approximates (in the sense of Definition 1).
Proof.
Consider a family of smooth functions satisfying the advection equation
Suppose that is an approximation to with
and that satisfies the discrete advection equation
Since approximates , given , we can choose such that
Therefore,
and
Thus, we have shown that approximates . However, satisfies
and satisfies
where is the matrix satisfying the equation
Therefore, we see that approximates . Thus, implies that . ∎
2.3. Discrete Commutator.
A spacediscrete flow that approximates a continuous flow is defined to be a smooth path in the space of orthogonal signed stochastic matrices, such that (see Definition 2) and (see Definition 3). It is straightforward to show that the Lie algebra structure of the space of divergencefree vector fields is preserved by our discretization. Indeed, if two matrices and approximate vector fields and then their commutator approximates the commutator of the Lie derivative operators:
Since , we obtain where denotes both the commutator of vector fields and the commutator of matrices. This property will be very useful to deal with Lin constraints later on.
2.4. Nonholonomic Constraints (NHC)
For a smooth path the matrix describes the infinitesimal
exchanges of fluid particles between any pair of cells and . We will thus
assume that is nonzero only if cells and
share a common boundary, i.e., are immediate neighbors. This sparsity will be numerically advantageous later on to reduce the computational complexity of the resulting integration schemes. We thus choose to restrict discrete paths on to those for which satisfies this constraint
More formally, we define the constrained set as the set of matrices corresponding to exchanges between neighboring cells only, i.e., if and only if implies that the cells and are neighbors. In this case the matrix is defined by a set of nonzero values defined on faces between adjacent cells and . As mentioned previously, to indicate their adjacency, we will write that and , where refers to the set of indices of adjacent cells to cell in the mesh . We will say that a matrix belongs to the class if implies . Finally, we will denote by , the constrained set at the identity. Consequently, our treatment of fluid dynamics will only consider matrices in , i.e., matrices in satisfying the sparsity constraints.
Note that if two matrices and both satisfy the constraints, their commutator need not: while the element of the commutator corresponding to any pair of cells which are more than two cells away is zero, the element may be nonzero when cells and are “two cells away” from each other since
Notice that the commutator is zero for neighboring cells since due to their antisymmetry. Writing , one sees that , where is the zero matrix. Therefore, the constraints we just defined are nonholonomic.
Remark. When a discrete vector field is in , the nonzero values of the antisymmetric matrix can be understood as dual chains, i.e., dimensional chains on the dual of [38]. This connection with 1chains will become crucial later when dealing with advection of curves to derive a discrete Kelvin’s theorem in Section 4.3.
2.5. Relation Between Elements of and Fluxes.
Suppose we have a family of discrete flows which approximates a flow such that approximates and satisfies the NHC. Let’s see how individual elements of are related to spatial values of . Recall that
is a discrete version of the advection equation
and in the norm. But it also means that is an approximation to the integral , i.e.,
(7) 
where is the normal vector to the boundary of and denotes the inner product of vectors. However,
where is the face shared by cells and , and the normal vector to oriented from to . By comparing this result to equation (7), it is clear that an element can be considered (up to a constant) as an approximation to the flux of a vector field through :
We know that , where is the barycenter of the boundary and is the area of . Therefore, we obtain that, up to a constant dependent on local mesh measures, approximates the flux through the boundary between and , i.e.,
In the case of a Cartesian grid of size this formula simplifies to:
2.6. Towards Lagrangian Dynamics with Nonholonomic Constraints.
One of the goals of this paper is to approximate geodesic flows on by Lagrangian flows on . To achieve this goal, we first need to define a Lagrangian such that
(8) 
and
(9) 
Such a Lagrangian, depending only on to mimic the continuous case, can then be used to formulate fluid dynamics through a discrete Lagranged’Alembert principle (to account for the nonholonomic constraint we impose on the sparsity of our Eulerian velocity approximation):
Note that the constraint on the variations of will induce a constraint on the variations of , giving rise to a discrete version of the wellknown Lin constraints of the form with (see Section 4.2).
However, we will show in later sections that coming up with a proper Lagrangian will require great care. As is typical with nonholonomic systems, the dynamics on will depend strongly on the values of (i.e., the matrix with as its element) outside of the constraint set because of the commutator present in the Lin constraints. In particular, a conventional discretization of the kinetic energy via the sum of all the squared fluxes on the grid would lead to a matrix with only values on pairs of adjacent cells, resulting in no dynamics. Instead, the Lagrangian must depend on values where .
To satisfy properties (8) and (9), we will look for a Lagrangian of the form
where the discrete innerproduct will be defined to satisfy the following properties (where denotes the continuous inner product of vector fields): for all ,
and for all
(10) 
where is the continuous flat operator (see for instance [1]). These properties will guarantee that conditions (8) and (9) are satisfied, and will lead to the proper dynamics. In the next section we will present a discretization of differential forms and a few operators acting on them to help us construct the discrete inner product (or equivalently, the discrete flat operator ).
3. StructurePreserving Spatial Field Discretization
We now introduce a discrete calculus consistent with our discretization of vector fields. Unlike previous discrete exterior calculus approaches, mostly based on chains and cochains (see [17, 7, 4] and references therein), we clearly distinguish between discrete vector fields and discrete forms acting on them. Moreover, our notion of forms will need to act not only on vector fields satisfying the NHC (being thus very reminiscent of the chain/cochain approach), but also on vector fields resulting from a commutator as imposed by the Lin constraints. We also introduce a discrete contraction operator and a discrete Lie derivative to complete our set of spatial operators—we will later show that the algebraic definition of our Lie derivative matches its dynamic counterpart as expected. We will not make any distinction in symbols between the discrete and continuous exterior calculus operators (, , , , etc) as the context will make their meaning clear.
3.1. Discrete Zeroforms.
In our context, a discrete form is a function that is piecewise constant per cell as previously defined in Section 2. Note that its representation is a vector of cell values,
where represents the value of the function in cell . Also, the volume integral of such a discrete form is obtained by weighting the value of each cell by the Lebesgue measure of this cell, and summing all contributions:
3.2. Discrete Oneforms.
As the space of matrices is used to discretize vector fields, a natural way to discretize oneforms is to also use matrices to respect the duality between these two entities. Moreover, it is in line with the previous definition for forms that were encoded as a 1tensor: forms will now be encoded by a tensor. Notice that this is also reminiscent of the approximation used in discrete mechanics [36].
Discrete Contraction. We define the contraction operator by a discrete vector field , acting on a discrete oneform to return a discrete zeroform, as:
(11) 
Notice the metricindependence of this definition, and that if the discrete vector field contains only nonzero terms for neighboring cells, any term where cell and cell are not neighbors does not contribute to the contraction. In this case, the value of the resulting form for cell is thus: , which is a local sum of the natural pairings of and on each face of cell .
Discrete Total Pairing. With this contraction defined, we derive a total pairing between a discrete form and a discrete vector field as:
This definition satisfies the following connection with the contraction defined in Eq. (11): indeed, for all
Note that the volume form is needed to integrate the piecewiseconstant form over the entire domain as explained in Section 3.1. Finally, since the matrix is antisymmetric, the symmetric component of does not play any role in the pairing.
Therefore, we will assume hereafter that a discrete oneform is defined by an antisymmetric matrix: .
Remark I. When viewed as acting on vector fields in the NHC space , our representation of discrete forms coincides with the use of cochains on the dual of [17]: the value (resp., ) can be understood as the integral of a continuous form on the oriented dual edge going from cell to cell (resp., from to ). However, our use of antisymmetric matrices extends this cochain interpretation. This will become particularly useful when forms need to be paired with vector fields that have the form of the commutator of two vector fields and both in as in Eq. (10).
Remark II. Notice finally that we can also define the notion of contraction of the volume form by a discrete vector field using . The resulting matrix can be thought of as a discrete two form encoding the flux of over each mesh face as derived in Section 2.5. In the notation convention of [28], this would be called a “primal” form, while the forms we will work with in this paper are “dual” forms. We won’t discuss these primal forms further in this paper (as the construction of a consistent discrete calculus of forms and tensors is a subject on its own), but it is clear that they naturally pair with dual forms (forming a discrete wedge product between primal  and dual forms), numerically resulting in the same value as the discrete pairing .
3.3. Discrete Twoforms.
We extend our definition of oneforms to twoforms in a similar fashion: discrete forms will be encoded as tensors that are completely antisymmetric, i.e., antisymmetric with respect to any pair of indices.
Discrete Contraction. Contraction of a form by a vector field is defined as:
Notice again here that the resulting discrete form is indeed an antisymmetric matrix (by construction), and that if , many of the terms in the sum vanish.
Discrete Total Pairing. The total pairing of a discrete form by two discrete vector fields and , the discrete equivalent of , will be defined as:
(12) 
This definition satisfies the expected property linking contraction and pairing: for all
Indeed, using our previous definitions, we have:
3.4. Other Operators on Discrete Forms.
A few more operators acting on , , or forms will be valuable to our discretization of incompressible fluids.
Discrete Exterior Derivative. We can easily define a discrete version of the exterior derivative. For a discrete 0form the oneform is defined as
Similarly, if is a discrete oneform then we can define:
More generally, we define our operator as acting on a form through:
where indicates the omission of a term. This expression respects the antisymmetry of our discrete form representation.
Remark. Notice here again that when the circumcenters of cells form a simplex on the dual of mesh , our definition of simply enforces Stokes’ theorem and thus coincides with the discrete exterior derivative widely used in the literature [17]. Our discrete exterior derivative extends this simple geometric property to arbitrary tuples of cells, while trivially enforcing that on the discrete level as well.
Discrete Lie Derivative. Now that we have defined contraction and derivatives on discrete oneforms we can define the Lie derivative using Cartan’s “magic” formula in the continuous setting