1 Introduction

Structure-Preserving Discretization of Incompressible Fluids


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 volume-preserving 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 numerical-analytic 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 volume-preserving diffeomorphisms using a finite dimensional Lie group, and associated discrete Euler equations are derived from a variational principle with non-holonomic constraints. The resulting discrete equations of motion yield a structure-preserving time integrator with good long-term 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, Ebin-Marsden and others; however the geometric-differential 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 Lagrangian-Eulerian 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 (so-called 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, structure-preserving theory for incompressible perfect fluids based on Hamilton-d’Alembert’s principle. Such a discrete variational approach to fluid dynamics guarantees invariance under the particle-relabeling 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 long-term 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.

Figure 1. Our geometric approach to discretizing the dynamics of incompressible fluids leads to discrete, structure-preserving, Lie group integrators. Here, six frames of an animation simulating heated smoke rising around a round obstacle in a closed box of incompressible fluid.

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 volume-preserving 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:


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 :


subject to constrained variations (called Lin constraints), where is an arbitrary divergence-free 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 volume-preserving 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 right-invariant discrete equivalent to the Eulerian velocity through its Lie algebra, i.e., through antisymmetric matrices whose columns sum to zero. After imposing non-holonomic constraints on the velocity field to allow transfer only between neighboring cells during each time update, we apply the Lagrange-d’Alembert principle (a variant of Hamilton’s principle applicable to non-holonomic systems) to obtain the discrete equations of motion for our fluid representation. As we will demonstrate, the resulting Eulerian variational Lie-group integrator is structure-preserving, and as such, has numerous numerical properties, from momentum preservation (through a discrete Noether theorem) to good long-term energy behavior.

Figure 2. Spatial Discretization: two cells and , with their common face of area and its dual edge of length .

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 one-forms 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 Volume-preserving diffeomorphisms on
Tangent space of at Id Divergence-free 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 null-row 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,
Table 1. Physical/Geometric meaning of the basic (continuous and discrete) variables used throughout this document.

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 CMMI-0757106, CCF-0811373, and DMS-0453145.

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 divergence-free 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 lemma1).

If the diffeomorphism is volume-preserving, 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


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:


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 : volume-preservation, 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 volume-preserving 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 volume-preserving and constant-preserving with respect to the mesh if, for all in




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 doubly-stochastic.

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 semi-group, 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 non-linear Hodge-like 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 divergence-free 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 volume-preserving 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 one-parameter family of volume-preserving diffeomorphisms and the associated time-dependent 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 divergence-free. Lie algebra elements for arbitrary simplicial meshes will be called null-row -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).


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



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 space-discrete 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 divergence-free 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. Non-holonomic 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 non-zero 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 constraint2. In other words, we only consider null-row -antisymmetric matrices satisfying the constraints as valid discrete vector fields. The non-zero elements of these matrices correspond to boundaries between adjacent cells and , and can be interpreted as directional transfer densities (per second) from to —they could abusively be called “fluxes” on regular grids; but we will make the proper link with the integrals of the velocity field over mesh faces in the next section.

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 non-zero 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 non-zero 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 non-holonomic.

Remark. When a discrete vector field is in , the non-zero values of the antisymmetric matrix can be understood as dual -chains, i.e., -dimensional chains on the dual of  [38]. This connection with 1-chains 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.,


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 Non-holonomic 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




Such a Lagrangian, depending only on to mimic the continuous case, can then be used to formulate fluid dynamics through a discrete Lagrange-d’Alembert principle (to account for the non-holonomic 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 well-known 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 -inner-product will be defined to satisfy the following properties (where denotes the continuous inner product of vector fields): for all ,

and for all


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. Structure-Preserving 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 Zero-forms.

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:

Remark. Our definition of -forms is no different from dual -cochains in dimension as used extensively in, e.g., [38, 17]. They naturally pair with dual -chains (i.e., linear combinations of cell circumcenters).

3.2. Discrete One-forms.

As the space of matrices is used to discretize vector fields, a natural way to discretize one-forms 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 1-tensor: -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 one-form to return a discrete zero-form, as:


Notice the metric-independence of this definition, and that if the discrete vector field contains only non-zero 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 piecewise-constant -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 one-form 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 Two-forms.

We extend our definition of one-forms to two-forms 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:


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 0-form the one-form is defined as

Similarly, if is a discrete one-form 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 one-forms we can define the Lie derivative using Cartan’s “magic” formula in the continuous setting