A variational $$\boldsymbol{H}({\rm div})$$ finite-element discretization approach for perfect incompressible fluids

A variational $$\boldsymbol{H}({\rm div})$$ finite-element discretization approach for perfect... Abstract We propose a finite-element discretization approach for the incompressible Euler equations which mimics their geometric structure and their variational derivation. In particular, we derive a finite-element method that arises from a nonholonomic variational principle and an appropriately defined Lagrangian, where finite-element $$\boldsymbol{H}({\rm div})$$ vector fields are identified with advection operators; this is the first successful extension of the structure-preserving discretization of Pavlov et al. (2009) to the finite-element setting. The resulting algorithm coincides with the energy-conserving scheme proposed by Guzmán et al. (2016). Through the variational derivation, we discover that it also satisfies a discrete analogous of Kelvin’s circulation theorem. Further, we propose an upwind-stabilized version of the scheme that dissipates enstrophy while preserving energy conservation and the discrete Kelvin’s theorem. We prove error estimates for this version of the scheme, and we study its behaviour through numerical tests. 1. Introduction The equations of motion for perfect fluids with constant density, namely the incompressible Euler equations, possess a rich geometric structure that is directly responsible for their many properties. More precisely, one can derive these equations from a variational principle using geometric techniques that also apply to classical mechanical systems such as the rigid body. This geometric view of perfect fluids is due to Arnol’d (1966). It has been extended to a wide variety of fluid models including magnetohydrodynamics (MHD) and geophysical fluid dynamics (GFD) models through the incorporation of advected quantities (see Holm et al., 1998). In this article, we derive finite-element numerical methods for perfect incompressible fluids, aiming to preserve as much of this geometric structure as possible. We pursue this objective for two main reasons. On the one hand, by addressing structure preservation in numerical simulations, we expect to preserve features of the solution that are usually lost by standard discretization approaches. On the other hand, preserving the structure of the equations at the discrete level enables us to devise discretizations with similar properties for different fluid models in a unified and systematic way. The idea of structure preservation in numerical simulations of the Euler equations is not new. Probably, the most successful approach in this context in terms of number of conservation laws is the sine truncation due to Zeitlin (1991), further developed in the form of a geometric integrator by McLachlan (1993). This uses a modified spectral truncation of the original system, yielding a Lie–Poisson structure of the equations at the discrete level. Such a truncation, however, is only well defined in the very specific case of a two-dimensional periodic domain, and generalizations are far from trivial. One of the most recent approaches is the one associated to the school of Brenier (1989, 2000), see, e.g., the discretizations proposed in Gallouët & Mérigot (2017) and Mérigot & Mirebeau (2016). Such methods are based on the reformulation of the variational principle underlying the incompressible Euler equations into an optimal transport problem. This leads to the definition of the notion of generalized incompressible flow, in which the motion of the fluid particles is defined in a probabilistic sense. The methods in Gallouët & Mérigot (2017) and Mérigot & Mirebeau (2016) aim at producing approximations of the flow map in this setting. In particular, they construct minimizing geodesics between discrete measure-preserving maps transporting the particle mass distributions. Moreover, they establish convergence to classical solutions of the Euler equations. In this article, we follow a different route, which is inspired by the work of Pavlov et al. (2009). In this case, the discrete flow maps are discretized by considering their action on piecewise constant scalar functions intended as right composition. However, such maps are not constructed explicitly, and the final algorithm is fully Eulerian, i.e., one approximates the evolution of the velocity field at fixed positions in space, without following each fluid particle. One key ingredient to achieve this is the identification of vector fields with discrete Lie derivatives, i.e., advection operators, also acting on piecewise constant scalar functions. Because of this identification, one can reduce the variational principle governing the dynamics from the space of discrete flow maps to the one of discrete velocity fields, i.e., from a Lagrangian (material) description to an Eulerian (spatial) one. The final algorithm can be interpreted as an advection problem for the velocity one-form by means of another discrete Lie derivative operator, defined consistently with the discrete exterior calculus (DEC) formalism (see Desbrun et al., 2005). Desbrun et al. (2014) extended this promising idea to generate a variational numerical scheme for rotating stratified fluid models. Both the algorithms presented by Pavlov et al. (2009) and Desbrun et al. (2014) can be classified as low-order finite difference methods and require a dual grid construction to generate an appropriate discrete Lagrangian. We will be concerned with generalizing the approach presented by Pavlov et al. (2009) using finite-element methods, with the aim of producing higher order variational schemes that require a single computational grid. As we are interested in discretizing incompressible fluids, we need finite-element vector spaces that encode the divergence-free constraint in their construction. This problem is generally solved by means of mixed finite-element methods (see, e.g., Boffi et al., 2013) that effectively enable us to build such spaces as subspaces of $$\boldsymbol{H}(\mathrm{div})$$ (the Hilbert space of square-integrable vector fields with square-integrable divergence on a given domain). The mathematical structure of mixed problems is particularly clear when formulating them in terms of differential forms, as Arnold et al. (2006) did by introducing the finite-element exterior calculus (FEEC) framework. In FEEC, many different mixed finite-element spaces can be seen as specific instances of polynomial differential form spaces. One can analyse different discretizations at once; most importantly, stability can be addressed in a systematic way independently of the domain topology. From our perspective, using differential forms in the discretization will be useful to have a single definition for a discrete Lie derivative, both when acting on scalar functions and differential one-forms (which can be identified with vector fields). Mullen et al. (2011) defined a numerical approach in this direction in the context of DEC using semi-Lagrangian techniques. Heumann & Hiptmair (2011) developed similar techniques formulated as Galerkin methods using the differential form spaces arising from FEEC. In the same work, they defined Eulerian-type discretizations including upwind stabilization, which they applied to linear advection/advection–diffusion problems with Lipschitz-continuous advecting velocity (see also Heumann & Hiptmair, 2013, for an application to the magnetic advection problem). Recently, these techniques were further developed to allow for piecewise Lipschitz-continuous advecting velocities (Heumann et al., 2016). Our strategy will be to adopt this latter approach and apply it to the case where the vector field generating the Lie derivative is in $$\boldsymbol{H}({\rm div})$$ finite-element spaces. In this way, we will be able to reproduce the structure of the variational algorithm in the work of Pavlov et al. (2009), transplanting the main features to the finite-element setting. This will lead us to rediscover the centred flux discretization described in Guzmán et al. (2016), and inspired by the discretization for the Navier–Stokes equations proposed in Cockburn et al. (2005). Our derivation shows that such a discretization, besides conserving kinetic energy, also inherits a discrete version of Kelvin’s circulation theorem. As a matter of fact, variational formulations lead to conservation laws resulting from symmetries in Hamilton’s principle, as a consequence of Noether’s theorem. Unfortunately, just as in Pavlov et al. (2009), this approach does not lead to conservation of enstrophy. In turn, this yields suboptimal convergence rates as observed in Guzmán et al. (2016), because a priori control of enstrophy is required for stability of smooth solutions. To overcome these issues, we introduce an upwind-stabilized version of the algorithm, which is different from the one studied in Guzmán et al. (2016). More specifically, we construct the method in such a way that it still possesses the same properties as the centred flux discretization, even though it cannot be derived from first principles. We prove suboptimal-order $$s$$ convergence when polynomial spaces of order $$s$$ are used, for $$s\,{>}\,1\,{+}\,n/2$$, where $$n$$ is the domain dimension. However, numerical results suggest that, at least in two dimensions, the scheme converges optimally with rate $$s+1$$ for $$s\geq 1$$. This article is organized as follows. In Section 2, we introduce the notions of one-forms and Lie derivatives, and we use these to review the variational structure underlying the incompressible Euler equations. In Section 3, we describe our discretization approach. In particular, we derive the central and upwind discretizations and we describe their properties. In Section 4, we present a convergence proof for the upwind scheme, using as starting point the error estimate for the centred scheme derived in Guzmán et al. (2016). We provide some numerical tests in Section 5. Finally, in Section 6, we conclude by describing further possible developments of our approach. 2. Preliminaries This section provides the basic tools necessary to formulate our numerical approach. We proceed in two steps. First, in Section 2.1, we introduce notation for one-forms and Lie derivatives. Then, in Section 2.2, we review the continuous geometric formulation of perfect incompressible fluids, and show how the equations of motion can be formulated as an advection problem for the momentum one-form. We refer to Abraham et al. (1988) and Marsden & Ratiu (2010) for more details on these topics. 2.1 One-forms and Lie derivatives Let $$V$$ be a vector space. A one-form on $$V$$ is an element of $$V^*$$ the dual of $$V$$, i.e., a linear functional $$a:V\,{\rightarrow}\, \mathbb{R}$$. For any $${\bf v}\,{\in}\, V$$, we denote by $$a(\bf v)\,{\in}\, \mathbb{R}$$ the pairing of $$a$$ with $$\bf v$$. Consider now an $$m$$-dimensional manifold $${\it {\Omega}}\subset \mathbb{R}^n$$. At each point $$p\in {\it {\Omega}}$$, we have a vector space $$T_p{\it {\Omega}}$$, namely the tangent space, whose elements are the tangent vectors to $${\it {\Omega}}$$ at $$p$$. The union of such spaces is the tangent bundle of $${\it {\Omega}}$$, $$T{\it {\Omega}} := \bigcup_{p\in {\it {\Omega}}} T_p{\it {\Omega}}$$. A vector field on $${\it {\Omega}}$$ is a smooth assignment $${\bf v}: {\it {\Omega}} \rightarrow T{\it {\Omega}}$$ such that $${\bf v}_p$$ is tangent to $${\it {\Omega}}$$ at $$p$$, i.e., $${\bf v}_p \in T_p{\it {\Omega}}$$. Differential one-forms generalize the concept of one-form as dual elements of vector fields on $${\it {\Omega}}$$. Specifically, we define the cotangent space $$T^*_p{\it {\Omega}}$$ to be the dual of $$T_p{\it {\Omega}}$$ and introduce the cotangent bundle $$T^*{\it {\Omega}} := \bigcup_{p\in {\it {\Omega}}} T^*_p{\it {\Omega}}$$. Then, a differential one-form is an assignment $$a:{\it {\Omega}} \rightarrow T^*{\it {\Omega}}$$ such that $$p \in {\it {\Omega}} \rightarrow a_p\in T^*_p{\it {\Omega}}$$. For any smooth vector field $$\bf v$$, the pairing $$a({\bf v})$$ is the real-valued function $$a({\bf v}): p \in {\it {\Omega}} \rightarrow a_p({\bf v}_p)\in \mathbb{R}$$. Then, we define the total pairing of a differential one-form $$a$$ with a vector field $$\bf v$$ to be ⟨a,v⟩:=∫Ωa(v)vol, (2.1) with $$\mathrm{vol}$$ being the measure on $${\it {\Omega}}$$ induced by the Euclidean metric. In this article, we restrict our study to the special case when $${\it {\Omega}}$$ is an $$n$$-dimensional domain in $$\mathbb{R}^n$$, with $$n=2,3$$, and admits global Cartesian coordinates $$\{x_i\}_{i=1}^n$$ oriented accordingly to an orthonormal basis $$\{{\bf e}_i\}_{i=1}^n$$ of $$\mathbb{R}^n$$. Then, a vector field $$\bf v$$ on $${\it {\Omega}}$$ has the coordinate representation $${\bf v} = \sum_{i=1}^n v_i {\bf e}_i$$, where $$v_i$$ are real-valued functions on $${\it {\Omega}}$$. Note that $${\bf e}_i$$, for $$i=1,\ldots,n$$, is interpreted as a constant vector field on $${\it {\Omega}}$$. Now, let $$\{{\mathrm{d}} x_i\}_{i=1}^n$$ be the dual basis to $$\{{\bf e}_i\}_{i=1}^n$$, i.e., $${\mathrm{d}} x_i({\bf e}_j) =\delta_{ij}$$. Then, a differential one-form $$a$$ has the coordinate representation $$a = \sum_{i=1}^n a_i\,{\mathrm{d}} x_i$$, where $$a_i$$ are real-valued functions on $${\it {\Omega}}$$, and $$a({\bf v}) = \sum_{i=1}^n a_i v_i$$. Note that $${{\mathrm{d}} x}_i$$, for $$i=1,\ldots,n$$, is interpreted as a constant differential one-form on $${\it {\Omega}}$$. In this simple setting, we have a trivial identification of a differential one-form $$a$$ with a vector field $$\bf a$$ by the equation a=a⋅dx, (2.2) where the right-hand side represents the formal inner product between the vector $${\bf a} = (a_1,\ldots,a_n)$$ and $${\mathrm{d}} {\bf x} = ({\mathrm{d}} x_1, \ldots, {\mathrm{d}} x_n)$$. Let $${\it {\Gamma}}\subset {\it {\Omega}}$$ be a smooth curve, parameterized by the functions $$x_i = x_i(s)$$, for $$i=1,\ldots,n$$ and $$s\in I$$, where $$I$$ is an open interval of $$\mathbb{R}$$ with $$0\in I$$. The quantity $${\bf s} = \sum_{i=1}^n{\mathrm{d}} x_i/{\mathrm{d}} s|_{s=0} \,{\bf e}_i$$ is a vector field on $${\it {\Gamma}}$$, since for each $$p\in {\it {\Gamma}}$$, $${\bf s}_p\in T_p{\it {\Gamma}}$$. Then, the integral of a differential one-form on $${\it {\Gamma}}$$ is defined by ∫Γa:=∫Ia(s)ds, (2.3) which is independent of the parameterization of $${\it {\Gamma}}$$. A smooth vector field $$\bf v$$ on $${\it {\Omega}}$$ can be identified with a one-parameter family of local diffeomorphisms $$\varphi: I \times U \rightarrow {\it {\Omega}}$$, with $$U$$ being an open set of $${\it {\Omega}}$$ and $$I\subset \mathbb{R}$$ an open set around $$0\in \mathbb{R}$$, such that for all $$t\in I$$, $$\varphi_t$$ is a local diffeomorphism defined by ∂φt∂t=v∘φt, (2.4) and $$\varphi_0 = e$$ the identity map on $${\it {\Omega}}$$. The map $$\varphi$$ defines the flow of the vector field $$\bf v$$. The Lie derivative of a differential one-form $$a$$ with respect to the vector field $$\bf v$$ is the one-form $${\sf L}_{\bf v} a$$ satisfying the equation ∫ΓLva=ddt|t=0∫φt∘Γa, (2.5) for any smooth curve $${\it {\Gamma}} \subset {\it {\Omega}}$$. Moreover, if $$n=3$$, we have the identification Lva=(grad(v⋅a)+(curla)×v)⋅dx, (2.6) or if $$n=2$$, Lva=(grad(v⋅a)+(rota)v⊥)⋅dx, (2.7) where at each $$p\in{\it {\Omega}}$$, $${\bf{v}}_p^\perp$$ is the vector $$\bf v_p$$ after a clockwise rotation of $$\pi/2$$, and $$\mathrm{rot}\,{\bf{a}}:=\mathrm{div}\,{\bf{a}}^\perp$$ (see Abraham et al. (1988) for more details). The Lie derivative can be regarded as a generalization of the scalar advection operator. To realize this, we need to find an equivalent of equation (2.5) for scalar functions. For a given scalar function $${b}:{\it {\Omega}}\rightarrow \mathbb{R}$$, we can interpret pointwise evaluation as the natural analogue of the integration of differential one-forms over curves (see, e.g., Heumann & Hiptmair, 2011). With this substitution, equation (2.5) becomes Lvb=ddt|t=0b∘φt=v⋅gradb, (2.8) which is just the directional derivative of $${b}$$ in the direction of $$\bf v$$. We conclude this section by introducing some notation for the function spaces and relative inner products we will use in this article. We denote by $$L^2({\it {\Omega}})$$ (respectively, $$\boldsymbol{L}^2({\it {\Omega}})$$) the space of square-integrable functions (respectively, vector fields) on $${\it {\Omega}}$$. We denote by $$(\cdot, \cdot)_{{\it {\Omega}}}$$ and $$\|\cdot\|_{{\it {\Omega}}}$$ the inner product and associated norm for both $$L^2({\it {\Omega}})$$ and $$\boldsymbol{L}^2({\it {\Omega}})$$. The definitions extend to one-forms, i.e., for any one-form $$a={\bf a}\cdot{{\mathrm{d}}} {\bf x}$$ and $$b={\bf b}\cdot{{\mathrm{d}}} {\bf x}$$, (a,b)Ω:=⟨a,b⟩=(a,b)Ω. (2.9) Moreover, we denote by $$L^p({\it {\Omega}})$$, with $$1\,{\leq}\, p \,{\leq}\, \infty$$, the standard generalization of $$L^2({\it {\Omega}})$$ with norm $$\|\cdot\|_{L^p({\it {\Omega}})}$$. The spaces $$H^k({\it {\Omega}})$$ (respectively, $$W^{k,p}({\it {\Omega}})$$) are the spaces of scalar functions with derivatives up to order $$k>0$$ in $$L^2({\it {\Omega}})$$ (respectively, $$L^p({\it {\Omega}})$$); the relative norms are $$\|\cdot\|_{H^k({\it {\Omega}})}$$ and $$\|\cdot\|_{W^{k,p}({\it {\Omega}})}$$. Similar definitions hold for vector fields; for example, $$\boldsymbol{L}^p({\it {\Omega}})$$ is the space of vector fields with components in $$L^p({\it {\Omega}})$$. Finally, we define H(div,Ω):={v∈L2(Ω):divv∈L2(Ω)}. (2.10) 2.2 Geometric formulation of incompressible perfect fluids In this section, we review the formal variational derivation of the incompressible Euler equations. This will be important in Section 3.2, where we develop our discretization approach by repeating this derivation step by step in a finite-dimensional setting. The diffeomorphisms on $${\it {\Omega}}$$ play a major role in the description of fluids, as they can be used to represent the motion of particles in the absence of discontinuities such as fractures, cavitations or shocks. More specifically, denote by $$\mathrm{Diff}({\it {\Omega}})$$ the group of (boundary-preserving) diffeomorphisms from $${\it {\Omega}}$$ to itself, with group product given by map composition. The subgroup of volume-preserving diffeomorphisms $$\mathrm{Diff}_\mathrm{vol} ({\it {\Omega}})$$ is defined by Diffvol(Ω):={g∈Diff(Ω):det(Dg)=1}, (2.11) where $$Dg$$ denotes the Jacobian matrix of the map $$g$$. This definition ensures that for any $$g \in \mathrm{Diff}_\mathrm{vol} ({\it {\Omega}})$$ and any open set $$U\subset {\it {\Omega}}$$, $$g(U)$$ has the same measure as $$U$$. Then, the motion of an incompressible fluid is represented by a curve on $$\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$, i.e., a one-parameter family of diffeomorphisms $$\{g_t\}_{t\in I}$$, where $$I$$ is an open interval of $$\mathbb{R}$$ containing $$0\in \mathbb{R}$$, so that for all $$t\in I$$, $$g_t \in \mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ and $$g_0 = e$$ the identity map on $${\it {\Omega}}$$. This means that, for each point $$X \in {\it {\Omega}}$$, $$g_t(X)$$ denotes the position occupied at time $$t$$ by the particle that was at $$X$$ when $$t=0$$. We define the material velocity field, ∂gt∂t(X)=U(t,X). (2.12) This is a vector field on $${\it {\Omega}}$$ that gives the velocity at time $$t\in I$$ of the particle that occupied the position $$X\in {\it {\Omega}}$$ at time $$t=0$$, for each $$X\in {\it {\Omega}}$$. The spatial velocity is the vector field u(t,x):=U(t,gt−1(x)). (2.13) This gives the velocity of the particle occupying the position $$x\in {\it {\Omega}}$$ at time $$t$$. Therefore, in the material formulation, $$X\in {\it {\Omega}}$$ is to be considered as a particle label, whereas in the spatial formulation, $$x\in {\it {\Omega}}$$ is a fixed position in the physical space. By the fact that $$g_t\in \mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$, we obtain that any spatial velocity field is tangent to the boundary $$\partial {\it {\Omega}}$$ and has zero divergence. We also see that the spatial velocity is left unchanged when composing on the right $$g_t$$ with any time independent diffeomorphism $$q \in \mathrm{Diff}({\it {\Omega}})$$, i.e., ∂∂t(gt∘q)∘(gt∘q)−1=∂∂tgt∘gt−1. (2.14) The interpretation of this is simple: composing the map $$q$$ with $$g_t$$ on the right is the same as relabelling the particles in the domain $${\it {\Omega}}$$. On the other hand, the spatial velocity field gives the velocity at fixed positions in space and it is not dependent on the specific particle labels. The geometric picture of incompressible perfect fluids is based on interpreting $$G:= \mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ as a configuration space, then applying the reduction techniques of dynamical systems on Lie groups (see, e.g., Marsden & Ratiu, 2010). In particular, $$\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ is not a Lie group, and one should provide a rigorous characterization of $$\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ to proceed in such a direction (see, e.g., Ebin & Marsden, 1970). Here, we are not concerned with these details, and we will proceed formally. A map $$g \in G$$ characterizes the configuration of the system by assigning to each particle label $$X\in {\it {\Omega}}$$ its physical position in the domain $$x = g(X)$$. Let $${\bf U} \in T_g G$$. Then there exists a one-parameter family $$\{g_t\}_{t\in I}$$, defined as above, such that for a fixed $$\bar{t}\in I$$, $$g_{\bar t} = g$$ and $${\frac{\partial}{\partial {t}}}|_{t=\bar{t}}\, g_t = {\bf U}$$. By comparison with equation (2.12), $$\bf U$$ can be interpreted as a material velocity field on $${\it {\Omega}}$$. In the classical Lagrangian setting, the evolution of the system is described by a curve on $$TG$$. In our case, this would be equivalent to the material formulation, where we follow the motion as a time-dependent diffeomorphism $$g_t$$ together with the associated material velocity field as defined in equation (2.12). The tangent space at the identity $$\mathfrak{g}:= T_eG$$ is the Lie algebra of $$G$$. By equation (2.13), we can regard elements of $$\mathfrak{g}$$ as spatial velocity fields on $${\it {\Omega}}$$. Therefore, an element of $$\mathfrak{g}= T_e G$$ can be identified with a vector field on $${\it {\Omega}}$$ tangent to its boundary $$\partial {\it {\Omega}}$$ and with zero divergence. For all $${\bf u},{\bf v}\in \mathfrak{g}$$, we denote by $$({\bf u},{\bf v}) \mapsto \mathrm{ad}_{\bf u} {\bf v}$$ their Lie algebra product. Let $$q_s$$ be a curve on $$G$$ such that $$q_0 = 0$$ and $$\dot{q}_0 = \bf v$$, the map $$\mathrm{ad}:\mathfrak{g} \times \mathfrak{g} \rightarrow \mathfrak{g}$$ is defined by aduv:=ddt|t=0dds|s=0gt∘qs∘gt−1=ddt|t=0(Dgtv)∘gt−1=−[u,v], (2.15) where $$[{\bf u}, {\bf v}]$$ is the Jacobi–Lie bracket, defined in coordinates by [u,v]:=∑i,j=1n(uj∂vi∂xj−vj∂ui∂xj)ei. (2.16) The Lagrangian of the system is a function $$L:TG\rightarrow \mathbb{R}$$ and is given by the kinetic energy of the spatial velocity field, i.e., L(g,g˙)=12∫Ω‖g˙∘g−1‖2vol. (2.17) Such a Lagrangian is right invariant because of equation (2.14). Therefore, the dynamics can be expressed in terms of the reduced Lagrangian $$l({\bf u})\,{:=}\, L(e, {\bf u})$$ only. More precisely, Hamilton’s principle for the Lagrangian in equation (2.17) is equivalent to δ∫t0t1l(u)dt=0, (2.18) with variations $$\delta {\bf u} = \dot{{\bf v}} - \mathrm{ad}_{{\bf u}} {\bf v}$$, where $${\bf v}_t \in \mathfrak{g}$$ is an arbitrary curve and $${\bf v}_{t_0}\,{=}\, {\bf v}_{t_1}\,{=}\,0$$. The variational principle in equation (2.18) produces the Euler–Poincaré equations, ⟨ddtδlδu,v⟩+⟨δlδu,aduv⟩=0, (2.19) for all $${\bf v}\in \mathfrak{g}$$, where we define ⟨δlδu,v⟩:=ddϵ|ϵ=0l(u+ϵv). (2.20) Now $$\delta l/\delta \bf u$$ can be identified with a one-form $$m$$ on $${\it {\Omega}}$$, which we will refer to as momentum, imposing that ⟨δlδu,v⟩=∫Ωm(v)vol, (2.21) for all $${\bf v}\in \mathfrak{g}$$. In our case, we readily see that $${m} = u := {\bf u}\cdot{{\mathrm{d}} {\bf x}}$$; however, we will distinguish between momentum one-forms $${m}$$ and velocity one-forms $$u$$, as this will be useful when describing model problems characterized by more complex Lagrangians. Inserting equations (2.15) and (2.21) into equation (2.19), we obtain ∫Ωm˙(v)vol−∫Ωm([u,v])vol=0. (2.22) If $$n=3$$, for divergence free vector fields $$\bf u$$ and $$\bf v$$ we have the useful identity $$[{\bf u}, {\bf v}] = - {\bf{curl}}({\bf u} \times {\bf v})$$. Hence, integrating by parts, the second integral can be rewritten as follows, ∫Ωm([u,v])vol =−∫Ωm⋅curl(u×v)vol =∫Ω(u×curlm)⋅vvol =−∫ΩLum(v)vol. (2.23) The same result holds for $$n=2$$. Note that the $${\bf{grad}}$$ term in equation (2.6) or (2.7) vanishes upon pairing with $$\bf v$$ because $$\mathrm{div}\,{\bf v} =0$$. Reinserting equation (2.23) into equation (2.22) yields the incompressible Euler equations in variational form ⟨m˙,v⟩+⟨Lum,v⟩=0, (2.24) for all $${\bf v} \in \mathfrak{g}$$, where $$m=u:= \bf u \cdot {{\mathrm{d}} {\bf x}}$$. Remark 2.1 Equation (2.23) establishes a link between the bracket on $$\mathfrak{g}$$ coming from Lagrangian reduction and the advection of the momentum one-form $${m}$$ by the Lie derivative operator. In the following we will develop a discretization approach that maintains this property, and in this way provides a discrete version of both the bracket and the Lie derivative. 3. Finite-element variational discretization In this section, we construct a variational algorithm for the incompressible Euler equations using discrete Lie derivatives intended as finite-element operators as the main ingredient. As starting point for this, we describe the Galerkin discretization of Lie derivatives proposed by Heumann & Hiptmair (2011) (see also Heumann & Hiptmair, 2013), in the case where the advecting velocity is a smooth vector field. Then, in Section 3.1, we consider a generalization for advecting velocities in finite-element spaces (meaning that the velocities have reduced smoothness). Such a generalization coincides with the one proposed in Heumann et al. (2016) for the more general case of piecewise Lipschitz-continuous advecting velocities. Next, in Section 3.2, we construct a variational algorithm using these operators as discrete Lie algebra variables, adapting the approach of Pavlov et al. (2009) to the finite-element setting. Moreover, we show that such an algorithm coincides with the centred flux discretization described in Guzmán et al. (2016). In Section 3.3, we include upwind stabilization to the method derived in the previous section. Finally, in Section 3.4, we show that both the centred flux and upwind schemes conserve energy and possess a discrete version of Kelvin’s circulation theorem. Throughout this section, we assume $$n=3$$, although all the results of this section extend without major modifications to the case $$n=2$$. 3.1 Discrete Lie derivatives Let $${\mathscr T}_h$$ be a triangulation on $${\it {\Omega}}$$, i.e., a decomposition of $${\it {\Omega}}$$ into simplices $$K\,{\subset}\, {\it {\Omega}}$$, with $$h$$ being the maximum simplex diameter. We define $$V_h\subset L^2({\it {\Omega}})$$ to be a chosen scalar finite-element space on $${\mathscr T}_h$$, and $${\bf W}_h \subset\boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ to be a chosen $$\boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$-conforming finite-element space of vector fields on $${\it {\Omega}}$$. We define $$W_h:= {\bf W}_h\cdot{{\mathrm{d}} {\bf x}}$$, i.e., $$W_h$$ is the set of one-forms $$a = {\bf a} \cdot {\mathrm{d}} {\bf x}$$ with $${\bf a} \in {\bf W}_h$$. More specifically, for $$r\geq 0$$, we denote by $$V^r_h\subset L^2({\it {\Omega}})$$ a scalar finite-element space such that, for all elements $$K\in {\mathscr T}_h$$, $$V^r_h|_K = \mathscr P_{r}(K),$$ where $$\mathscr P_r(K)$$ is the set of polynomial functions on $$K$$ of degree up to $$r$$. In other words, $$V_h^r$$ is the standard scalar discontinuous Galerkin space of order $$r$$. Analogously, we denote by $${\bf W}^r_h \subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ either the Raviart–Thomas space of order $$r$$, defined for $$r\geq 0$$ by RTr(Th):={u∈H(div,Ω):u|K∈(Pr(K))n+xPr(K)}, (3.1) or the Brezzi–Douglas–Marini space of order $$r$$, and defined for $$r\geq 1$$ by BDMr(Th):={u∈H(div,Ω):u|K∈(Pr(K))n}. (3.2) Later, we will use the divergence-free subspace of $${\bf W}^r_h$$, which is the same for $${\bf{RT}}_r(\mathscr T_h)$$ and $${\bf{BDM}}_r(\mathscr T_h)$$ (see, e.g., Boffi et al., 2013, Corollary 2.3.1). Let $$\mathscr F$$ be the sets of $$(n-1)$$-dimensional simplices in $${\mathscr T}_h$$, and let $${\mathscr F}^\circ\subset {\mathscr F}$$ be the set of facets $$f\in {\mathscr F}$$ such that $$f \cap \partial {\it {\Omega}} = \emptyset$$. We fix an orientation for all $$f\in {\mathscr F}$$ by specifying the unit normal vector $${\bf n}_f$$ in such a way that on the boundary $${\bf n}_f$$ is equal to the outward pointing normal $${\bf n}_{\partial {\it {\Omega}}}$$. An orientation on the facet $$f\in {\mathscr F}$$ defines a positive and negative side on $$f$$, so that any $$a\in V_h$$ has two possible restrictions on $$f$$, denoted $$a^+$$ and $$a^-$$, with $$a^+:= a|_{K^+}$$ and $$K^+\in \mathscr T_h$$ being the element with outward normal $${\bf{n}}_f$$. Therefore, we can define for each $$a \in V_h$$ and on each $$f\in \mathscr F^\circ$$ the jump operator $$[a] := a^+ - a^-$$ and the average operator $$\{a\} := (a^+ + a^-)/2$$. Such definitions extend without modifications to vector fields $${\bf a} \in {\bf W}_h$$ or to their one-form representation $$a = {\bf a}\cdot {{\mathrm{d}} {\bf x}}\in W_h$$. Lemma 3.1 Let $${\bf W}_h\subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$, then for all $${\bf a}\in {\bf W}_h$$ and for all $$f\in \mathscr F^\circ$$, $${\bf a}\cdot {\bf n}_f$$ is single valued on $$f$$, and in particular $$[{\bf a}]\cdot {\bf n}_f = 0$$. Proof. See Lemma 5.1 in Arnold et al. (2006). □ For all elements $$K\,{\in}\,\mathscr T_h$$ and $$f\,{\in}\,\mathscr F$$, we denote by $$(\cdot,\cdot)_K$$ and $$(\cdot,\cdot)_f$$ the standard $$L^2$$ inner products on $$K$$ and $$f$$. We use the same definition for both scalar and vector quantities, and we extend it to one-forms as in equation (2.9). Furthermore, given a smooth vector field $$\boldsymbol \beta$$, for any $$f\,{\in}\,\mathscr F$$ and for any smooth scalar functions $$a$$ and $$b$$ on $$f$$ we define (a,b)f,β:=(β⋅nfa,b)f. (3.3) Similarly, for any element $$K\in\mathscr T_h$$ and for any smooth scalar functions $$a$$ and $$b$$ on $$\partial K$$ define (a,b)∂K,β:=(β⋅n∂Ka,b)∂K, (3.4) where $${\bf n}_{\partial K}$$ is the unit normal to the boundary of the element $$K$$ pointing outwards. The same definitions hold if $$a$$ and $$b$$ are one-forms or vector fields with smooth components on $$f$$ or $$\partial K$$. A discrete Lie derivative in this article means any finite-element operator from $$V_h$$ to itself, or from $$W_h$$ to itself, which is also a consistent discretization of the Lie derivative. In Heumann & Hiptmair (2011), the authors proposed two types of discrete Lie derivatives1 based on either an Eulerian or a semi-Lagrangian approach, which were used in the context of linear advection and advection–diffusion problems. We will take as starting point their Eulerian discretization (see Section 4 in Heumann & Hiptmair, 2011), where the advecting velocity is a fixed smooth vector field $$\boldsymbol \beta$$. In the following, we will consider extensions of this approach to the case, where $$\boldsymbol \beta$$ is not fully continuous; specifically $$\boldsymbol \beta \in {\bf W}_h\subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$. The resulting operators coincide with the ones proposed in Heumann et al. (2016) for piecewise Lipschitz-continuous advecting velocities. Such a discussion will be useful to discretize the nonlinear advection term in the Euler equations. The Eulerian-type discrete Lie derivative proposed by Heumann & Hiptmair (2011) is given by the following definition. Definition 3.2 Let $$\boldsymbol{\beta}$$ be a smooth vector field on $${\it {\Omega}}$$ tangent to the boundary $$\partial {\it {\Omega}}$$. The finite-element operator $${\sf L}_{\boldsymbol{\beta}}^h: V_h \rightarrow V_h$$ (respectively, $${\sf L}_{\boldsymbol{\beta}}^h: {W}_h \rightarrow {W}_h$$) is defined by (Lβha,b)Ω=∑K(Lβa,b)K+∑f∈F∘(−([a],{b})f,β+([a],[b])f,cfβ), (3.5) where $$c_f:f \rightarrow \mathbb{R}$$ is a scalar function depending on $$\boldsymbol \beta$$, for all $$a,b\in V_h$$ (respectively, $$a,b \in W_h$$). The consistency of such a discretization is proved by observing that if the advected quantity $$a$$ is smooth, $${\sf L}^h_{\boldsymbol{\beta}} a$$ coincides with the Galerkin projection of $${\sf L}_{\boldsymbol{\beta}}a$$ onto $$V_h$$ or $$W_h$$. Note that the Lie derivative discretization on $$V_h$$, as defined in equation (3.5), coincides with the classical discontinuous Galerkin discretization of the advection operator described in Brezzi et al. (2004). This can be verified using integration by parts in equation (3.5), yielding the following equivalent definition for the operator $${\sf L}_{\boldsymbol{\beta}}^h: V_h \rightarrow V_h$$, (Lβha,b)Ω=−∑K(a,div(βb))K+∑f∈F∘(({a},[b])f,β+([a],[b])f,cfβ), (3.6) for all $$a,b\in V_h$$. Similarly, the operator $${\sf L}_{\boldsymbol{\beta}}^h: W_h \rightarrow W_h$$ can be equivalently defined by (Lβha,b)Ω=∑K(a,curl(β×b)−βdivb)K+∑f∈F∘(({a},[b])f,β+([a],[b])f,cfβ), (3.7) for all $$a,b\in W_h$$ with $$a = {\bf a}\cdot {{\mathrm{d}}}{\bf x}$$ and $$b = {\bf b} \cdot {{\mathrm{d}}{\bf x}}$$. Details on the calculations leading to equations (3.6) and (3.7) can be found in Heumann (2011). Remark 3.3 There are two choices for the function $$c_f$$ that are particularly important: $$c_f = 0$$ (centred discretization). In this case, equation (3.6) reduces to (Lβha,b)Ω=−∑K(a,div(βb))K+∑f∈F∘({a},[b])f,β, (3.8) for all $$a,b\in V_h$$. We refer to such discretization as centred, since the facet integrals contain the average of the advected quantity $$a$$; $$c_f = \boldsymbol{\beta} \cdot{\bf n}_f /(2|\boldsymbol{\beta} \cdot{\bf n}_f|)$$ (upwind discretization). In this case, equation (3.6) reduces to (Lβha,b)Ω=−∑K(a,div(βb))K+∑f∈F∘(aup,[b])f,β, (3.9) for all $$a,b\in V_h$$, where $$a^{\rm up} = a^+$$ if $$\boldsymbol{\beta}\cdot {\bf n}_f\geq0$$ and $$a^{\rm up} = a^-$$ if $$\boldsymbol{\beta}\cdot {\bf n}_f<0$$, i.e., in the facet integrals $$a$$ is always evaluated from the upwind side. A similar interpretation holds when considering the action of $${\sf L}_{\boldsymbol{\beta}}^h$$ on $$W_h$$. Lemma 3.4 Let $$n=3$$ and let $$\boldsymbol{\beta}$$ be a smooth vector field on $${\it {\Omega}}$$ tangent to the boundary $$\partial {\it {\Omega}}$$. Then, the discrete Lie derivative $${\sf L}^h_{\boldsymbol{\beta}}:W_h\rightarrow W_h$$ in equation (3.5) can be equivalently defined by (Lβha,b)Ω= ∑K(a,curl(β×b)−βdivb)K +∑f∈F∘((nf×{a},β×[b])f+(cfnf×[a],β×[b])f), (3.10) for all $$a,b\in W_h$$ with $$a = {\bf a}\cdot {{\mathrm{d}}}{\bf x}$$ and $$b = {\bf b} \cdot {{\mathrm{d}}{\bf x}}$$. Proof. The proof relies on the fact that if $${\bf a} \in {\bf W}_h\subset\boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ then, by Lemma 3.1, $${\bf a}\cdot {\bf n}_f$$ is single-valued across elements of $$\mathscr T_h$$, so $$[{\bf a}]\cdot {\bf n}_f = 0$$ on all facets $$f\in \mathscr F^\circ$$. Then, the result follows by applying to the facet integrals in equation (3.7) the vector identity (A×B)⋅(C×D)=(A⋅C)(B⋅D)−(B⋅C)(A⋅D), (3.11) for all $${\bf{A}},{\bf{B}},{\bf{C}},{\bf{D}}\in \mathbb{R}^3$$. See Heumann (2011) for details. □ As we are concerned with discretizing the Euler equations, we want to extend Definition 3.2 to the case where the advecting velocity also belongs to a finite-element space. We will limit our study to the case $$\boldsymbol{\beta} \in {\bf W}_h\subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$. This will prove to be enough for our scope since if $$\boldsymbol{\beta} \in \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$, we can ensure that $$\mathrm{div}\, \boldsymbol{\beta} = 0$$. Remark 3.5 Equation (3.5) is still well defined if $$\boldsymbol{\beta} \in {\bf W}_h$$, since in this case, by Lemma 3.1, $$\boldsymbol{\beta} \cdot {\bf n}_f$$ is single valued on each facet $$f \in {\mathscr F}^\circ$$. In view of Remark 3.5, equation (3.5) can be used without modification to extend the definition of the discrete Lie derivative to the case $$\boldsymbol{\beta} \in {\bf W}_h$$. However, this is not the only possible extension: there exist many other discretization approaches that reduce to Definition 3.2 for continuous advecting velocity. In particular, Lemma 3.4 suggests the following alternative definition, which is also the one adopted in Heumann et al. (2016). Definition 3.6 Let $$n=3$$ and let $$\boldsymbol{\beta}\in{\bf W}_h$$ such that $$\boldsymbol{\beta} \cdot{\bf n}_{\partial {\it {\Omega}}} =0$$ on $$\partial {\it {\Omega}}$$. The finite-element operator $${\sf X}_{\boldsymbol{\beta}}^h:V_h \rightarrow V_h$$ is defined by $${\sf X}_{\boldsymbol{\beta}}^h \,a := {\sf L}_{\boldsymbol{\beta}}^h \, a$$ for all $$a \in V_h$$, whereas the operator $${\sf X}_{\boldsymbol{\beta}}^h:W_h \rightarrow W_h$$ is defined by (Xβha,b)Ω= ∑K(a,curl(β×b)−βdivb)K +∑f∈F∘((nf×{a},[β×b])f+(cfnf×[a],[β×b])f), (3.12) where $$c_f:f \rightarrow \mathbb{R}$$ is a scalar function depending on $$\boldsymbol \beta$$, for all $$a,b\in W_h$$ with $$a ={\bf a} \cdot {\mathrm{d}} {\bf x}$$ and $$b ={\bf b} \cdot {\mathrm{d}} {\bf x}$$. In other words, $${\sf X}_{\boldsymbol{\beta}}^h$$ coincides with $${\sf L}^h_{\boldsymbol{\beta}}$$ when applied to scalar functions, but not when applied to one-forms. Proposition 3.7 For a continuous vector field $$\boldsymbol \beta$$, we have the equivalence $${\sf X}_{\boldsymbol{\beta}}^h={\sf L}_{\boldsymbol{\beta}}^h$$. Therefore, $${\sf X}_{\boldsymbol{\beta}}^h$$ is a discrete Lie derivative, i.e., it is a consistent discretization of the Lie derivative operator. Proof. This is immediate by comparing equation (3.10) with equation (3.12) and observing that if $$\boldsymbol \beta$$ is continuous, for any $${\bf b}\in {\bf W}_h$$, we have $$\boldsymbol{\beta} \times [{\bf b}] = [\boldsymbol{\beta} \times {\bf b}]$$ on all facets $$f\in {\mathscr F}^\circ$$. □ 3.2 Discrete incompressible Euler equations We now use the discrete Lie derivatives defined in the previous section to design a variational discretization approach for the incompressible Euler equations. The derivation presented here closely follows Pavlov et al. (2009), but is adapted to the context of finite-element discretizations and does not require a dual grid. First, we will define a discrete group that approximates the configuration space $$\mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$. Next, we will derive the relative Euler–Poincaré equations for an appropriately chosen Lagrangian. Finally, we will connect the resulting algorithm to the Lie derivative discretizations in Definition 3.6. Such an approach will lead us to rediscover the centred flux discretization proposed in Guzmán et al. (2016). The main issue related to deriving a variational integrator for the incompressible perfect fluids is that the configuration space $$G = \mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$ is infinite dimensional, so one needs to find an appropriate finite-dimensional approximation that converges, in the limit, to the original system. In Pavlov et al. (2009), this issue is solved by applying Koopman’s lemma, i.e., identifying group elements with their action, intended as right composition, on $$L^2$$ functions on $${\it {\Omega}}$$. This is to say that $$G$$ can be equivalently represented as a subgroup $$G(V) \subset GL (V)$$ of the group of invertible linear maps from $$V:= L^2({\it {\Omega}})$$ to itself, by means of a group homomorphism $$\rho:G \rightarrow G(V)$$ defined by ρ(g)⋅a=ρg⋅a:=a∘g−1 (3.13) for any $$g\,{\in}\,G$$ and $$a\,{\in}\,C^\infty({\it {\Omega}})\,{\subset}\,V$$. As the Lie algebra $$\mathfrak g$$ of $$G$$ is the set of divergence-free vector fields $$\bf v$$ tangent to the boundary, the Lie algebra $$\mathfrak g(V)$$ of $$G(V)$$ is the set of Lie derivatives $${\sf L}_{\bf v}: V \rightarrow V$$, interpreted as unbounded operators, with respect to such vector fields. As a matter of fact, we have dds|s=0ρgs⋅a=dds|s=0(a∘gs−1)=−Lva (3.14) for any $$a\in C^\infty({\it {\Omega}}) \subset V$$, and where $$g_s$$ is a curve on $$G$$ such that $$g_0=e$$ and $$\dot{g_s}|_{s=0} = {\bf v}$$. To get a discrete version of this picture, we start by restricting $$V$$ to be the finite-element space $$V_h$$. Therefore, we consider the finite-dimensional Lie group $$GL(V_h)$$ and we seek an appropriate subgroup $$G_h(V_h)$$ that in the limit approximates $$G(V)$$. Clearly, this cannot be accomplished by restricting $$G(V)$$ to $$V_h$$, as we have done for the general linear group $$GL(V)$$, because it would imply losing the group structure. However, one can still construct a subgroup $$G_h(V_h)\subset GL(V_h)$$ that approximates $$G(V)$$ and hence $$G$$. For all $$g^h\in GL(V_h)$$ and $$a\in V_h$$, we denote by $$a \mapsto g^h a$$ the action of $$g^h$$ on $$a$$. Then, $$G_h(V_h)$$ can be defined as follows. Definition 3.8 The Lie group $$G_h(V_h)\subset GL(V_h)$$ is the subgroup of the group of linear invertible operators from $$V_h$$ to itself, defined by the following properties: for all $$g^h\in G_h(V_h)$$,  (gha,ghb)Ω=(a,b)Ω∀a,b∈Vh, (3.15) ghc=c∀c∈R, (3.16) where $$c$$ is considered as a constant function on $${\it {\Omega}}$$. It is trivial to check that both equations (3.15) and (3.16) are verified in the continuous case. In particular for any $$g\in \mathrm{Diff}({\it {\Omega}})$$, we have $$c \circ g^{-1} = c$$. Moreover, as a consequence of volume preservation, for all $$a,b \in V$$ and $$g\in \mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$, ∫Ω(a∘g−1)(b∘g−1)vol=∫Ω(ab)∘g−1vol=∫Ωdet(Dg−1)abvol=∫Ωabvol (3.17) since $$g^{-1}({\it {\Omega}}) = {\it {\Omega}}$$. The Lie algebra of $$GL(V_h)$$ is the set $$\mathfrak{gl}(V_h)$$ of all linear operators from $$V_h$$ to itself. For all $${\sf A}^h\in \mathfrak{gl}(V_h)$$ and $$a\in V_h$$, we denote by $$a \mapsto {\sf A}^h a$$ the action of $${\sf A}^h$$ on $$a$$. Then, the Lie algebra of $$G_h(V_h)$$ can be defined as follows. Definition 3.9 The Lie algebra $$\mathfrak{g}_h(V_h)\subset \mathfrak{gl}(V_h)$$ associated to the Lie group $$G_h(V_h)\subset GL(V_h)$$ is the subalgebra of the linear operators from $$V_h$$ to itself defined by the following properties: for all $${\sf A}^h \in \mathfrak{g}_h(V_h)$$,  (Aha,b)Ω=−(a,Ahb)Ω∀a,b∈Vh, (3.18) Ahc=0∀c∈R. (3.19) Note that equations (3.18) and (3.19) can be derived directly from equations (3.15) and (3.16). In particular, equation (3.18) can be again related to volume preservation. Moreover, we can regard $$\mathfrak{g}_h(V_h)$$ as an approximation $$\mathfrak{g}(V)$$ so its elements are to be thought as discrete Lie derivatives in view of equation (3.14). The collection of all the spaces $$G_h(V_h)$$, for $$h$$ arbitrarily small, is not $$G(V)$$, even if the union of the spaces $$V_h$$ is dense in $$V$$. This is because there exist elements of $$GL(V)$$ that satisfy equations (3.15) and (3.16), but do not belong to $$G(V)$$ (for example, they might be representative of maps that are not even continuous). Therefore, before defining the discrete system we must restrict ourselves to a space smaller than $$G_h(V_h)$$. The same considerations hold at the Lie algebra level, meaning that one needs to constrain the dynamics in a subset $$S_h(V_h) \,{\subset}\, \mathfrak g_h(V_h)$$. Our definition for such a space is based on the discrete Lie derivatives of Definition 3.6. In particular, we fix a finite-element space $${\bf W}_h \subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ and define W∘h:={u∈Wh:divu=0,u⋅n∂Ω=0}. (3.20) We denote by $$\overset{\circ}{W}_{h}:= \overset{\circ}{\mathbf W}_{h}\cdot {{\mathrm{d}} {\bf x}}$$, i.e., $$\overset{\circ}{W}_{h}$$ is the set of one-forms $$u = {\bf u} \cdot {\mathrm{d}} {\bf x}$$ with $${\bf u} \in \overset{\circ}{\mathbf W}_{h}$$. Analogous definitions hold for $$\overset{\circ}{\mathbf W}_{h}^r$$ and $$\overset{\circ}{W}_{h}^r$$. Note that $$\overset{\circ}{\mathbf W}_{h}^r$$ is the same space if constructed using $${\bf W}_h^r = {\bf{RT}}_r(\mathscr T_h)$$ or $${\bf W}_h^r = {\bf{BDM}}_r(\mathscr T_h)$$ (see, e.g., Boffi et al., 2013, Corollary 2.3.1). The set of advection operators associated to $$\overset{\circ}{\mathbf W}_{h}$$ is given by Sh(Vh):={Xuh:Vh→Vhwithcf=0:u∈W∘h}. (3.21) More specifically, we define Shs(Vhr):={Xuh:Vhr→Vhrwithcf=0:u∈W∘hs}. (3.22) Note that setting $$c_f=0,$$ we immediately ensure that $$S_h(V_h)$$ is a linear space; this can be directly verified using Definition 3.6. Lemma 3.10 $$S_h(V_h) \subset \mathfrak{g}_h(V_h)$$. Proof. We need to verify that the elements of $$S_h(V_h)$$ satisfy equations (3.18) and (3.19). By Definition 3.6 and equation (3.5), for all $$a,b\in V_h$$, (Xuha,b)Ω=∑K(Lua,b)K−∑f∈F∘([a],{b})f,u. (3.23) The Lie derivative is just the directional derivative when applied to scalar functions, as shown in equation (2.8). Moreover, since $$\mathrm{div}\, {\bf u}=0$$, integration by parts yields ∑K(Lua,b)K =−∑K(a,Lub)K+∑K(a,b)∂K,u =−∑K(a,Lub)K+∑f∈F∘((a+,b+)f,u−(a−,b−)f,u). (3.24) Note that the decomposition of the integrals on $$\partial K$$ is justified by the fact that each facet is shared by the boundary of two adjacent elements, therefore it appears twice in the sum, but with two opposite normals (since $${\bf n}_{\partial K}$$ is always oriented outwards). Moreover, we have the identity (a+,b+)f,u−(a−,b−)f,u=({a},[b])f,u+([a],{b})f,u. (3.25) Combining equations (3.23), (3.24) and (3.25) gives (Xuha,b)Ω=−(a,Xuhb)Ω, (3.26) which is the same as in equation (3.18). Equation (3.19) is immediately verified by observing that if we set $$a=c$$ in equation (3.23), with $$c \in \mathbb{R}$$ being a constant function on $${\it {\Omega}}$$, we obtain (Xuhc,b)Ω=0, (3.27) for all $$b\in V_h$$, and therefore $${\sf X}_{\bf u}^h c =0$$. □ Lemma 3.11 Let $${\it {\Omega}}$$ be a domain in $$\mathbb{R}^n$$ and let $$\mathscr T_h$$ be a triangulation on $${\it {\Omega}}$$. For any $$r\geq s$$ and $$r\geq1$$, the spaces $$S_h^s(V^r_h)$$ and $$\overset{\circ}{\mathbf W}_{h}^s$$ are isomorphic. Proof. We assumed that $$\{{\bf e}_i\}_{i=1}^n$$ is a global orthonormal reference frame on $$\mathbb{R}^n$$ with coordinates $$\{x_i\}_{i=1}^n$$. Then, since $$r\geq 1$$, $$x_i\in V_h^r$$ and we can define the map $$\hat{\cdot}: {\mathfrak g}_h(V_h^r) \rightarrow (V_h^r)^n$$ by A^h:=∑i=1n(Ahxi)ei, (3.28) for all $${\sf A}^h \in {\mathfrak g}_h$$. By a standard result of mixed finite-element theory (see, e.g., Boffi et al., 2013, Corollary 2.3.1), for any $${\bf u} \in \overset{\circ}{\mathbf W}_{h}^s$$ we have that $${\bf u}|_K \in ({\mathscr P}_s(K))^n$$ for any $$K\in {\mathscr T}_h$$. Then, since $$r \geq s$$, for all $${\bf u} \in \overset{\circ}{\mathbf W}_{h}^s$$, $$u_i = {\bf u} \cdot {\bf e}_i \in V_h^r$$. Moreover, we can verify by direct calculation, i.e., applying Definition 3.6, that for all $${\sf X}^h_{\bf u}\in S_h^s(V_h^r)$$, $${\sf X}^h_{\bf u}\, x_i = u_i$$. Hence the restriction of the $$\hat{\cdot}$$ map to $$S_h^s(V^r_h)$$ is the map $${\sf X}^h_{\bf u} \mapsto {\bf u}$$ which is a surjection from $$S_h^s(V_h^r)$$ to $$\overset{\circ}{\mathbf W}_{h}^s$$. On the other hand, the map $${\bf u} \mapsto {\sf X}^h_{\bf u}$$ from $$\overset{\circ}{\mathbf W}_{h}^s$$ to $$S_h^s(V_h^r)$$ is also a surjection by definition of $$S_h^s(V_h^r)$$, hence the result. □ The different spaces introduced so far can be collected in the diagrams represented in Fig. 1. The definition of the $$\hat{\cdot}$$ map suggests the choice of a Lagrangian $$l:\mathfrak{g}_h(V_h)\rightarrow \mathbb{R}$$ given by Fig. 1. View largeDownload slide Diagrams representing the discretization of the group $$G:=\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ and its Lie algebra. The dashed arrows represent relations that hold only in an approximate sense. Fig. 1. View largeDownload slide Diagrams representing the discretization of the group $$G:=\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ and its Lie algebra. The dashed arrows represent relations that hold only in an approximate sense. l(Ah):=12(A^h,A^h)Ω, (3.29) so that the restriction of $$l$$ to $$S_h(V_h)$$ is just the total kinetic energy, i.e., $$l({\sf X}_{\bf u}^h) = \|{\bf u}\|^2_{{\it {\Omega}}}/2$$. Such a choice for the Lagrangian was already suggested in Pavlov et al. (2009) as an alternative to the Lagrangian construction proposed therein based on dual grids. To get a geometric picture analogous to the continuous case, we need a Lagrangian defined on the whole tangent space $$TG_h(V_h)$$ and approximating equation (2.17). We achieve this by extending $$l$$ by right translations to get $$L:TG_h(V_h)\rightarrow \mathbb{R}$$ defined by L(gh,g˙h):=l(g˙h(gh)−1)=12∫Ω‖g˙h(gh)−1^‖2vol, (3.30) which is right invariant by construction. Analogously, the space $$S_h(V_h)$$ can be extended by right translation to give a right invariant sub-bundle of $$TG_h(V_h)$$. To see how this is done, consider the right translation map $$R_{g^h} : q^h \in G_h(V_h) \rightarrow q^hg^h \in G_h(V_h)$$ and denote by $$T_eR_{g^h} : {\sf A}^h \in {\mathfrak g}_h(V_h) \rightarrow {\sf A}^h g^h \in T_{g^h}G_h(V_h)$$ its tangent map at the identity (see Marsden & Ratiu, 2010, for more details). Then, the right invariant extension of $$S_h(V_h)$$ is obtained by collecting the spaces $$T_eR_{g^h} S_h(V_h)$$ for all $$g^h\,{\in}\, G_h(V_h)$$, where $$T_eR_{g^h} S_h(V_h) \,{\subset} T_{g^h} G_h(V_h)$$ is the space of linear operators which can be written as the composition $${\sf A}^h g^h$$, for a given $${\sf A}^h\in {\mathfrak g}_h(V_h)$$. Lemma 3.12 The Euler–Poincaré–d’Alembert equations relative to the right invariant Lagrangian $$L$$ with constraint $$\dot{g}^h\in T_eR_{g^h} S_h(V_h)$$ are given by: Find $${\sf A}^h\in S_h(V_h)$$ such that (ddtA^h,B^h)Ω+(A^h,[Ah,Bh^])Ω=0, (3.31) for all $${\sf B}^h\in S_h(V_h)$$, where $$[{{\sf A}^h, {\sf B}^h}]:= {\sf A}^h {\sf B}^h - {{\sf B}^h {\sf A}^h}$$ is the commutator of the linear operators $${\sf A}^h$$ and $${\sf B}^h$$. Proof. The proof is analogous to that of Theorem 1 in Pavlov et al. (2009). In essence, it is just an application of the reduction theorem (Marsden & Ratiu, 2010, Theorem 13.5.3). In particular, note that proceeding as equation (2.15), we get $$\mathrm{ad}_{{\sf A}^h} {\sf B}^h = [{{\sf A}^h, {\sf B}^h}]$$. Then, equation (3.31) is the finite-dimensional analogue of equation (2.19). The only difference is that we added a nonholonomic constraint on the velocity so that both the solution and the variations are constrained in the space $$S_h(V_h)$$. □ Theorem 3.13 Under the hypotheses of Lemma 3.11, the Euler–Poincaré–d’Alembert equations (3.31) relative to the right invariant Lagrangian in equation (3.30) with $$S_h(V_h) = S^s_h(V^r_h)$$ are equivalent to: Find $$u \in \overset{\circ}{W}_{h}^s$$ such that (u˙,v)Ω+(Xuhu,v)Ω=0, (3.32) with $$u = {\bf u} \cdot {{\mathrm{d}} {\bf x}}$$, for all $${v} \in \overset{\circ}{W}_{h}^s$$. Proof. First, by Lemma 3.11 for each $${\sf A}^h\in S_h^s(V_h^r)$$ and $${\sf B}^h\in S_h^s(V_h^r),$$ there is a unique $${\bf u}\in \overset{\circ}{\mathbf W}_{h}^s$$ and $${\bf v}\in \overset{\circ}{\mathbf W}_{h}^s$$ such that $${\sf A}^h=-{\sf X}^h_{\bf u}$$ and $${\sf B}^h=-{\sf X}^h_{\bf v}$$. We now examine the two terms in equation (3.31) separately. For the first term, we have (ddtA^h,B^h)Ω=(u˙,v)Ω, (3.33) by definition of the one-form pairing in equation (2.9). As for the second term, we have (A^h,[Ah,Bh^])Ω =−(X^uh,[Xuh,Xvh^])Ω =−∑i=1n(ui,Xuhvi−Xvhui)Ω =−∑i=1n∑K(ui,Luvi−Lvui)K+∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v) =−∑K(u,[u,v])K+∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v), (3.34) where we used the identity $$\sum_{i=1}^n({\sf L}_{\bf u}v_i - {\sf L}_{\bf v} u_i) {\bf e}_i = [{\bf u}, {\bf v}]$$, which is a consequence of equations (2.8) and (2.16). As for the facet integrals, we first notice that ∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v)=∑f∈F∘(({u},[v])f,u−({u},[u])f,v). (3.35) Moreover, since $$\bf u$$ has continuous normal component on the mesh facets, we can substitute ({u},[v])f,u=({u},[v])f,{u}=(nf×{u},{u}×[v])f, (3.36) where for the last equality we have used the same reasoning as in the proof of Lemma 3.4. Similarly, for the second term in equation (3.35), we have ({u},[u])f,v=(nf×{u},{v}×[u])f=−(nf×{u},[u]×{v})f. (3.37) Therefore, equation (3.35) becomes ∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v) =∑f∈F∘(nf×{u},{u}×[v]+[u]×{v})f =∑f∈F∘(nf×{u},[u×v])f. (3.38) The last equality in equation (3.38) follows from the linearity of the cross product, i.e., {u}×[v]+[u]×{v} =(u++u−)/2×(v+−v−)+(u+−u−)×(v++v−)/2 =u+×v+−u−×v−=[u×v]. (3.39) Inserting equation (3.38) into equation (3.34) yields (A^h,[Ah,Bh^])Ω=∑K(u,curl(u×v))K+∑f∈F∘(nf×{u},[u×v])f, (3.40) where we used the identity $$[{\bf u},{\bf v}] = -{\bf{curl}}({\bf u}\times {\bf v})$$ valid for any divergence-free vector fields $${\bf u}$$ and $${\bf v}$$. Recalling that $$c_f=0$$, comparison with equation (3.12) concludes the proof. □ Proposition 3.14 The centred discretization in equation (3.32), i.e., the case $$c_f\,{=}\,0$$, on a simplicial triangulation coincides with the discontinuous Galerkin discretization with centred fluxes proposed in Guzmán et al. (2016), which can be written as follows: Find $$({\bf u},p) \,{\in}\, {\bf W}_h^s \times V_h^k$$, where $$k=s$$ if $${\bf W}_h^s = {\bf{RT}}_s(\mathscr T_h)$$, or $$k=s-1$$ if $${\bf W}_h^s = {\bf{BDM}}_s(\mathscr T_h)$$, such that {(u˙,v)Ω−∑K(u,u⋅∇v)K+∑f∈F∘(u⋅nf{u},[v])f−∑K(p,divv)K=0∑K(divu,q)K=0  (3.41) for all $$({\bf v}, q) \in {\bf W}_h^s \times V_h^k$$, with $${\bf u} \cdot {\bf n}_{\partial {\it {\Omega}}} = 0$$ on $$\partial {\it {\Omega}}$$, and $$\int_{{\it {\Omega}}} p \, {\mathrm{d}} x = 0$$. Proof. First, we observe that restricting equation (3.41) on the subspace $$\overset{\circ}{\mathbf W}_{h}^s$$ yields the equivalent problem: Find $${\bf u} \in \overset{\circ}{\mathbf W}_{h}^s$$ such that (u˙,v)Ω−∑K(u,u⋅∇v)K+∑f∈F∘(u⋅nf{u},[v])f=0, (3.42) for all $${\bf v} \in \overset{\circ}{\mathbf W}_{h}^s$$. On the other hand, we can write equation (3.32) in the following form (u˙,v)Ω+∑K(u,curl(u×v))K+∑f∈F∘(nf×{u},[u×v])f=0. (3.43) We rewrite the volume integrals in equation (3.43) as follows, (u,curl(u×v))K=(u,−u⋅∇v+v⋅∇u)K=−(u,u⋅∇v)K+(∇u2/2,v)K, (3.44) for any $$K \in \mathscr {T}_h$$. We rewrite the facet integrals in equation (3.43) as follows, ∑f∈F∘(nf×{u},[u×v])f =∑f∈F∘((u⋅nf{u},[v])f−(v⋅nf{u},[u])f) =∑f∈F∘((u⋅nf{u},[v])f−(v⋅nf,[u2/2])f) =∑f∈F∘(u⋅nf{u},[v])f−∑K∫Kdiv(v(u2/2))dx =∑f∈F∘(u⋅nf{u},[v])f−∑K(v,∇(u2/2))K, (3.45) where we used equation (3.38) for the equality in the first line. Inserting equations (3.44) and (3.45) into equation (3.43) gives the equivalence of the two algorithms. □ 3.3 Upwind scheme In this section, we study the upwind version of the scheme in equation (3.32). In particular, we extend the discussion of the previous section to the upwind case, with the aim of preserving the main features of the variational derivation. Including upwinding in the variational framework described above is not straightforward. This is because, if we define $$S_h(V_h)$$ by Sh(Vh):={Xuh:Vh→Vhwithcf=c¯f(u):u∈W∘h}, (3.46) where $$\bar{c}_f$$ is a given function of $$\bf u$$, then $$S_h(V_h)$$ ceases to be a linear space, and furthermore $$S_h(V_h)\nsubseteq \mathfrak{g}_h(V_h)$$. This last issue can be solved easily by choosing a group larger than $$G_h(V_h)$$ as configuration space. More specifically, we need to allow for discrete diffeomorphisms, which do not satisfy equation (3.15). We call the new configuration space $$\tilde{G}_h(V_h)$$ and we define it as follows. Definition 3.15 The Lie group $$\tilde{G}_h(V_h)\subset GL(V_h)$$ is the subgroup of the group of linear invertible operators from $$V_h$$ to itself, defined by the following property: for all $$g^h\in \tilde{G}_h(V_h)$$, ghc=c∀c∈R, (3.47) where $$c$$ is considered as a constant function on $${\it {\Omega}}$$. Denote with $$\tilde{\mathfrak{g}}_h(V_h)$$ the Lie algebra of $$\tilde{G}_h(V_h)$$. Then, by the same arguments as in the last section, it is easy to verify that, for any choice of $$c_f$$ in the definition of $$S_h(V_h)$$, $$S_h(V_h)\subset \tilde{\mathfrak{g}}_h(V_h)$$. Next, we note that also Lemma 3.11 holds independently of the choice of $$c_f$$ in the definition of $$S_h(V_h)$$, and in particular, for fixed polynomial orders $$r$$ and $$s$$ with $$r\,{\geq}\,1$$ and $$r\,{\geq}\,s$$, if $${\sf A}^h \in S_h^s(V_h^r)$$ then $${\bf u}:= \hat{\sf A}^h \,{\in}\, \overset{\circ}{\mathbf W}_{h}^s$$. This observation leads to the following result. Proposition 3.16 (Upwinding). Under the hypotheses of Lemma 3.11, the problem described by equation (3.31), where $$S_h(V_h) = S_h^s(V_h^r)$$ is defined to be Shs(Vhr):={Xvh:Vhr→Vhrwithcf=c¯f(u):v∈W∘hs}, (3.48) with $${\bf u}\,{:=}\, \hat{\sf A}^h\,{\in}\, \overset{\circ}{\mathbf W}_{h}^s$$ is equivalent to equation (3.32), where $${\sf X}_{\bf u}^h$$ is defined as in Definition 3.6 with $$c_f=-{\bar{c}_f({\bf u})}$$. {In particular, if we choose $$\bar{c}_f({\bf u})=-{{\bf u} \cdot {\bf n}_f}/({2 |{\bf u} \cdot {\bf n}_f|})$$, we obtain the upwind version of equation (3.31). Proof. The proof is analogous to Theorem 3.13. □ Remark 3.17 Note that the upwind version of the scheme in equation (3.32) cannot be derived directly from Hamilton’s principle because the constraint space in equation (3.48) is dependent on the solution of the problem. However, in the next section, we will see that the equivalence between equations (3.32) and (3.31), ensured by Proposition 3.16, is enough to maintain the main properties of the centred discretization (i.e., energy conservation and a discrete version of Kelvin’s circulation theorem) also when upwinding is introduced. 3.4 Energy conservation and discrete Kelvin’s circulation theorem A consequence of using a nonholonomic constraint in our approach is that our discretization does not preserve all the Casimirs, i.e., the conserved quantities, of the original system. This is because the discrete bracket induced by the advection term in equation (3.32) is only a quasi-Poisson bracket (it is only skew symmetric and does not satisfy the Jacobi identity). Nonetheless, because of the variational derivation of the equations of motion, we are still able to ensure energy conservation and derive a discrete version of Kelvin’s circulation theorem, with or without upwinding. Energy conservation is established in the following proposition. Proposition 3.18 (Energy conservation/$$L^2$$ stability) The discrete system in equation (3.32), with or without upwinding, satisfies conservation of the total kinetic energy, i.e., $${\mathrm{d}}{l}/{\mathrm{d}} t=0$$. Proof. This can be verified directly by setting $$v=u$$ in equation (3.32) and exploiting the antisymmetric structure of the equations. □ We now turn to derive a discrete Kelvin’s circulation theorem. At the continuous level, the theorem states that for any closed loop $${\it {\Gamma}}$$, ddt∫gt∘Γu=0, (3.49) where $$g_t$$ is the flow of $${\bf u}$$, with $$u= {\bf u}\cdot {{\mathrm{d}} {\bf x}}$$. To show in what sense equation (3.49) is verified at the discrete level, we will proceed as in Pavlov et al. (2009) by replacing the integral in equation (3.49) by a pairing with currents, i.e., vector fields supported in a small region around a curve. A discrete current is a vector field $${\bf c}\in \overset{\circ}{\mathbf W}_{h}$$ that is nonzero only on a closed loop of elements $$K\,{\in}\,\mathscr T_h$$ and such that the flux of $${\bf c}$$ through adjacent elements of the loop is equal to 1. Clearly, with mesh refinement, a discrete current approaches the notion of closed loop. In other words, if $${\it {\Gamma}}$$ is a closed loop such that there exists a discrete current $${\bf c}$$ with $${\it {\Gamma}}\in \overline{\mathrm{supp}({\bf c})}$$, then ⟨u,c⟩≈∫Γu. (3.50) Let $$g_t$$ be the flow of $$\bf u$$ as defined by equations (2.12) and (2.13). If both $$g_t$$ and a given discrete current $$\bf c$$ are smooth, then $$\bf c$$ is advected by $$g_t$$ via the map $${\bf c} \mapsto (D g_t \, {\bf c}) \circ g_t^{-1}$$, as it can be derived from equation (2.15). Therefore, Kelvin’s theorem can be reformulated using currents, by stating that ddt⟨u,(Dgtc)∘gt−1⟩=0, (3.51) for any smooth current $$\bf c$$. In the discrete setting, $$g_t$$ can only be defined in weak sense, and in general, it is not differentiable, because we only have $${\bf u} \in \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$. Nonetheless, we can still use the discrete flow map $$g_t^h$$ to reproduce an analogous of equation (3.51), as shown in the following proposition. Proposition 3.19 Let $${\bf u} \in \overset{\circ}{\mathbf W}_{h}$$, with $$u={\bf u}\cdot {{\mathrm{d}} {\bf x}}$$, be a solution of the discrete Euler equations (3.32). Then, for all discrete currents $${\bf c}\in \overset{\circ}{\mathbf W}_{h}$$ ddt|t=0⟨u,(gth)−1Xchgth^⟩=0, (3.52) where $$g_t^h\in G_h(V_h)$$ for the centred discretization, or $$g_t^h\in \tilde{G}_h(V_h)$$ for the upwind discretization, is the discrete flow of $$-{\sf X}_{\bf u}^h$$, i.e., it satisfies $$\dot{g}^h_t =- {\sf X}_{\bf u}^h g^h_t$$ and $$g^h_0 = e$$. The result extends to any time $$t$$ by appropriately translating the definition of $$g_t^h$$ in time. Proof. For the centred discretization, let $$q_t^h\in G_h(V_h)$$ be the discrete flow of $$-{\sf X} ^h_{\bf c}$$ (or $$q_t^h\in \tilde{G}_h(V_h)$$, for the upwind discretization), then [Xuh,Xch]=adXuhXch=−ddt|t=0dds|s=0(gth)−1qshgth=ddt|t=0(gth)−1Xchgth. (3.53) Finally, by equation (3.31), we have ddt|t=0⟨u,(gth)−1Xchgth^⟩=⟨u˙,c⟩|t=0+⟨u,[Xuh,Xch^]⟩|t=0=0. (3.54) □ 4. Convergence analysis In this section, we analyse the convergence properties of the semidiscrete system in equation (3.32). In particular, as the centred discretization, i.e., the case $$c_f =0$$, was studied in Guzmán et al. (2016), we will concentrate on the upwind formulation, i.e., the case $$c_f = {\bf u} \cdot {\bf n}_f/ (2 |{\bf u} \cdot {\bf n}_f|)$$. For simplicity, we continue to restrict to the case $$n=3$$, although the results of this section extend to $$n=2$$. We keep the dependence on $$n$$ explicit when needed. We start with characterizing the space $${\bf W}_h^s$$ and its approximation properties. We will assume that $$\mathscr T_h$$ is quasi-uniform and shape regular and that $${\bf W}_h^s$$ is either $${\bf{RT}}_{s}(\mathscr T_h)$$ or $${\bf{BDM}}_{s}(\mathscr T_h)$$, i.e., the Raviart–Thomas or Brezzi–Douglas–Marini finite-element spaces of order $$s$$ on $${\mathscr T}_h$$. Then, as shown in Boffi et al. (2013) (or in Arnold et al., 2006, in the context of FEEC), there exists a projection operator $${\bf P}_h$$ of smooth vector fields onto $${\bf W}^s_h$$ such that ‖u−Phu‖Ω≤Chr|u|Hr(Ω), (4.1) for any vector field $${\bf u}\in \boldsymbol{H}^{r}({\it {\Omega}})$$ and $$1\leq r \leq s+1$$. From now on, we will denote by $${\bf u}$$ the exact solution of the Euler equation, and by $${\bf u}_h\,{\in}\,\overset{\circ}{\mathbf W}_{h}^s$$ the discrete solution obtained by solving equation (3.32). The estimate for the centred flux scheme in Guzmán et al. (2016) was derived as follows. First, the error is decomposed using the triangle inequality ‖u−uh‖Ω≤‖u−Phu‖Ω+‖Phu−uh‖Ω. (4.2) The first term on the right-hand side is the approximation error, therefore for the convergence estimate is sufficient to obtain a bound on $$\boldsymbol{\gamma}_h:= {\bf P}_h {\bf u} - {\bf u}_h$$. We collect the results of Guzmán et al. (2016) in the following lemma. Lemma 4.1 (Theorem 2.1 in Guzmán et al. (2016)). The error $$\boldsymbol{\gamma}_h$$ satisfies the following inequality ddt‖γh‖Ω2≤C0h2s+C1h2s+2+C2‖γh‖Ω2, (4.3) where $$C_0,C_1,C_2>0$$ are constants depending on $$\|{\bf u}\|_{\boldsymbol{W}^{1,\infty}({\it {\Omega}})}$$, $$\|{\bf u}\|_{\boldsymbol{H}^{r+1}({\it {\Omega}})}$$, $$\|\dot{\bf u}\|_{\boldsymbol{H}^{r+1}({\it {\Omega}})}$$, but independent of $$h$$. Then, if the norms above are uniformly bounded for $$t\in [0,T]$$, the following estimate holds ‖u−uh‖Ω≤Chs, (4.4) for an appropriate constant $$C>0$$ independent of $$h$$, and for all $$t\in [0,T]$$. When upwinding is introduced, i.e., for the case $$c_f ={\bf u} \cdot {\bf n}_f/ (2 |{\bf u} \cdot {\bf n}_f|)$$, we just need to add the upwind contribution to the left-hand side of equation (4.3) and proceed with the estimate. In particular, from the proof of Theorem 2.1 in Guzmán et al. (2016), it can be easily verified that, assuming the exact solution to be continuous, equation (4.3) needs to be replaced by ddt‖γh‖Ω2−∑f∈F∘(cfnf×[uh],[uh×γh])f≤C0h2s+C1h2s+2+C2‖γh‖Ω2. (4.5) Proceeding as in equation (3.34), we rewrite the upwind term as follows, ∑f∈F∘(cfnf×[uh],[uh×γh])f =∑f∈F∘(cfnf×[uh],{uh}×[γh]+[uh]×{γh})f =∑f∈F∘(([uh],[γh])f,cfuh−([uh],[uh])f,cfγh) =∑f∈F∘(([Phu],[γh])f,cfuh−([γh],[γh])f,cfuh−([uh],[uh])f,cfγh). (4.6) Reinserting equation (4.6) into equation (4.5) gives ddt‖γh‖Ω2+∑f∈F∘([γh],[γh])f,cfuh ≤C0h2s+C1h2s+2+C2‖γh‖Ω2 +∑f∈F∘([Phu],[γh])f,cfuh−∑f∈F∘([uh],[uh])f,cfγh, (4.7) where the facet integrals on the left-hand side are always non-negative due to the definition of $$c_f$$. Define I1up:=|∑f∈F∘([Phu],[γh])f,cfuh|,I2up:=|∑f∈F∘([uh],[uh])f,cfγh|, (4.8) then ddt‖γh‖Ω2+∑f∈F∘([γh],[γh])f,cfuh≤C0h2s+C1h2s+2+C2‖γh‖Ω2+I1up+I2up. (4.9) Estimate for$$I_1^{\rm up}$$ Since $$c_f {\bf u}_h\cdot {\bf n}_f = |{\bf u}_h \cdot {\bf n}_f|/2$$, the following estimate holds, I1up=|∑f∈F∘([Phu],[γh])f,cfuh| ≤12|∑f∈F∘([Phu],[Phu])f,cfuh+∑f∈F∘([γh],[γh])f,cfuh| ≤12∑f∈F∘([Phu−u],[Phu−u])f,cfuh+12∑f∈F∘([γh],[γh])f,cfuh ≤Ctrh−1‖uh‖L∞(Ω)‖Phu−u‖Ω2+12∑f∈F∘([γh],[γh])f,cfuh, (4.10) where we used the trace inequality $$\|\cdot\|^2_{\partial K} \leq C_{\rm tr} h^{-1} \|\cdot\|^2_{K}$$, with $$C_{\rm tr}>0$$ and independent of $$h$$. We estimate the $$\boldsymbol{L}^\infty$$ norm of $${\bf u}_h$$ as follows, ‖uh‖L∞(Ω)≤‖u‖L∞(Ω)+Cinvh−n/2(‖u−Phu‖Ω+‖γh‖Ω), (4.11) where we used the inverse inequality $$\|\cdot\|_{\boldsymbol{L}^{\infty}({\it {\Omega}})}\,{\leq}\,C_{\rm inv} h^{-n/2}\|\cdot\|_{\it {\Omega}}$$, with $$C_{\rm inv}\,{>}\,0$$ and independent of $$h$$. Inserting equation (4.11) into equation (4.10) and applying Young’s inequality to the term involving $$\|\boldsymbol{\gamma}_h\|_{\it {\Omega}}$$, we obtain I1up≤C3(h2s+1+h4s+2−n+h3s+2−n/2)+12‖γh‖Ω2+12∑f∈F∘([γh],[γh])f,cfuh, (4.12) for an appropriate constant $$C_3>0$$ independent of $$h$$. Estimate for$$I_2^{\rm up}$$ Proceeding as for $$I_1^{\rm up}$$, we obtain I2up=|∑f∈F∘([uh],[uh])f,cfγh| ≤2|∑f∈F∘([γh],[γh])f,cfγh|+2|∑f∈F∘([Phu−u],[Phu−u])f,cfγh| ≤4CtrCinvh−1−n/2(‖γh‖Ω3+‖γh‖Ω‖Phu−u‖Ω2) ≤C4h−1−n/2‖γh‖Ω3+C4h4s+2−n+12‖γh‖Ω2, (4.13) for an appropriate constant $$C_4>0$$ independent of $$h$$. Final estimate Assuming $$s\,{\geq}\,1$$ we can neglect higher order terms in $$h$$ in the estimates above. In particular, inserting equations (4.12) and (4.13) into equation (4.9), we obtain ddt‖γh‖Ω2−‖γh‖Ω2+∑f∈F∘([γh],[γh])f,cfuh≤C5h2s+C5h−1−n/2‖γh‖Ω3, (4.14) for a constant $$C_5>0$$ independent of $$h$$. We multiply both sides of equation (4.14) by $$e^{-t}$$, and we let $$q(t):= e^{-t}\|\boldsymbol{\gamma}_h\|_{{\it {\Omega}}}^2$$, so we obtain ddtq(t)+e−t∑f∈F∘([γh],[γh])f,cfuh≤e−tC5h2s+et/2C5h−1−n/2q(t)3/2. (4.15) Integrating from $$0$$ to $$T$$ gives q(T) ≤q(0)+(1−e−T)C5h2s+C5h−1−n/2∫0Tes/2q(s)3/2ds ≤C6h2s+C6h−1−n/2∫0Tes/2q(s)3/2ds, (4.16) where $$C_6$$ depends on $$\|{\bf u}|_{t=0}\|_{\boldsymbol{H}^{r+1}({\it {\Omega}})}$$, but not on $$h$$. Then, the nonlinear generalization of the Gronwall’s inequality in Butler & Rogers (1971) yields q(T)≤h2s(1C61/2−2C6hs−1−n/2(eT/2−1))−2, (4.17) with the condition hs−1−n/2<12C63/2(eT/2−1). (4.18) If $$s>1+n/2,$$ such a condition can be verified for any $$T$$, if $$h$$ is sufficiently small. Moreover, for hs−1−n/2<14C63/2(eT/2−1), (4.19) we get q(T)≤4C6h2s, (4.20) which proves the following theorem. Theorem 4.2 Under the same assumptions on $$\bf u$$ as in Lemma 4.1, the discrete solution $${\bf u}_h\in {\bf W}_h^s$$, with $$s>1+n/2$$, obtained by solving the system in equation (3.32) with $$c_f = {\bf u} \cdot {\bf n}_f/ (2 |{\bf u} \cdot {\bf n}_f|)$$, satisfies ‖u−uh‖Ω≤Chs, (4.21) with $$C>0$$ dependent $$T$$, $$\|{\bf u}|_{t=0}\|_{\boldsymbol{H}^{s+1}({\it {\Omega}})}$$, $$\|{\bf u}\|_{H^1([0,T],\boldsymbol{H}^{s+1}({\it {\Omega}}))}$$, $$\|{\bf u}\|_{\boldsymbol{W}^{1,\infty}([0,T]\times {\it {\Omega}})}$$, but not on $$h$$. Remark 4.3 Theorem 4.2 gives a suboptimal rate of convergence for the upwind scheme, i.e., $$c_f\,{=}{{\bf u}_h \cdot {\bf n}_f}/({2|{\bf u}_h \cdot {\bf n}_f|})$$, only for the case $$s>1+n/2$$. This result is rather unsatisfactory, since the numerical tests in the next section suggest that the scheme converges with the optimal rate $$s\,{+}\,1$$ for $$s\,{\geq}\,1$$. Unfortunately, we were not able to prove this result analytically. 5. Numerical tests We now show some numerical results illustrating our analytical results. All the tests presented here are computed in the domain $${\it {\Omega}} \,{=}\, [0,2\pi]\,{\times}\,[0,2\pi]$$ with spatial dimension $$n\,{=}\,2$$ and global Cartesian coordinates $$(x_1,x_2)$$. We combine the semidiscrete scheme introduced in the previous section with the implicit midpoint time discretization, with fixed time step $${\it {\Delta}} t$$. More specifically, let $$t^n := n {\it {\Delta}} t$$ for $$n \geq 0$$, and $$u^n \,{=}\, {\bf u}^n \cdot {\mathrm{d}} {\bf x} \,{:=}\, u|_{t=t^n}$$, then the fully discrete scheme is defined by Given $$u^n \,{\in}\, \overset{\circ}{W}_{h}$$, find $$u^{n+1}\,{\in}\, \overset{\circ}{W}_{h}$$ such that 1Δt(un+1−un,v)Ω+(Xun+un+12hun+un+12,v)Ω=0, (5.1) for all $$v \in \overset{\circ}{W}_{h}$$. It is well known that the implicit midpoint rule preserves the quadratic invariants of the semidiscrete system exactly. In our case, kinetic energy is conserved via Proposition 3.18. The tests of this section were performed using the Firedrake software suite (see Rathgeber et al., 2016), which allows for symbolic implementation of finite-element problems of mixed type. The nonlinear system was solved using Newton’s method, with LU factorization for the inner linear system implemented using the PETSc library (see Hendrickson & Leland, 1993; Balay et al., 1997, 2015; Dalcin et al., 2011). We start by examining the convergence rates with respect to a given manufactured solution. The manufactured solution solves a modified system of equations, where an appropriate forcing is introduced so that we can study the approximation properties of our discretization (although we cannot make conclusive remarks on stability). We pick as manufactured solution the Taylor–Green vortex (see Guzmán et al., 2016), u(t,x1,x2)=sin⁡(x1)cos⁡(x2)e−2t/σe1−cos⁡(x1)sin⁡(x2)e−2t/σe2, (5.2) with $${\it {\Omega}} = [0,2\pi]\times [0,2\pi]$$, $$\sigma=100$$ and $$t\in [0,1]$$. Table 1 lists the $$L^2$$-error and the convergence rates for both the centred and upwind scheme, for the case $${\bf W}_ h = {\bf{RT}}_s(\mathscr T_h)$$ and time step $${\it {\Delta}} t = 10^{-2}$$. The centred scheme was already analysed in Guzmán et al. (2016), where it was pointed out that without upwinding the rate of convergence $$s$$ is sharp for $$s=1$$, whereas convergence rates higher than $$s$$ can be observed for $$s=0$$ and $$s=2$$. The form of upwinding proposed in this article is different from the one of Guzmán et al. (2016). However, it gives the same behaviour in terms of order of convergence. In particular, we notice that introducing upwinding yields the optimal convergence rate $$s+1$$ for $$s\geq 1$$. Table 1 Comparison between centred ($$c_f=0$$) and upwind $$\left(c_f = \frac{{\mathbf u}\cdot{\mathbf n}_f}{2|{\mathbf u}\cdot{\mathbf n}_f|}\right)$$ scheme in terms of the error $$\|{\mathbf{u}}-\mathbf{u}_h\|_{{\it {\Omega}}}$$ and the order of convergence, for the Taylor–Green vortex test case and $${\mathbf W}_h={\mathbf{RT}}_s(\mathscr T_h)$$ Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 Table 1 Comparison between centred ($$c_f=0$$) and upwind $$\left(c_f = \frac{{\mathbf u}\cdot{\mathbf n}_f}{2|{\mathbf u}\cdot{\mathbf n}_f|}\right)$$ scheme in terms of the error $$\|{\mathbf{u}}-\mathbf{u}_h\|_{{\it {\Omega}}}$$ and the order of convergence, for the Taylor–Green vortex test case and $${\mathbf W}_h={\mathbf{RT}}_s(\mathscr T_h)$$ Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 In the following, we will refer to the scheme of Guzmán et al. (2016) as the GSS scheme (by the initial of the authors), whereas we refer to the scheme proposed in this article as the Lie derivative scheme. We now consider a double shear problem (see, e.g., Liu & Shu, 2000; Guzmán et al., 2016) obtained by setting the initial velocity $${\bf u} = u_1{\bf e}_1 + u_2 {\bf e}_2$$ as follows u1|t=0:={tanh⁡((x2−π/2)/ρ)x2≤πtanh⁡((3π/2−x2)/ρ)x2>π ,u2|t=0:=δsin⁡(x1), (5.3) where $$\rho = \pi/15$$, $$\delta = 0.05$$ and with periodic boundary conditions on $${\it {\Omega}}$$. We use $${\it {\Delta}} t = 8/200$$ as in the references above to integrate the solution in time. In Figs 2–4, we plot the vorticity $$\mathrm{rot}\, {\bf u}$$ at $$t=8$$, for the centred scheme and the upwind GSS and Lie derivative schemes. The centred scheme produces oscillations in the vorticity field, as it was already observed in Guzmán et al. (2016). The presence of oscillations is reflected in the growth in enstrophy $$Z(t)\,{:=}\,\int_{{\it {\Omega}}} (\mathrm{rot}\, {\bf u})^2\, {\mathrm{d}} {\bf x}$$, which accumulates at small scales, see Fig. 5. As a matter of fact, unfortunately, our variational derivation does not imply enstrophy conservation, even though this is an invariant of the continuous system. On the other hand, the energy $$K(t) \,{:=}\, \int_{{\it {\Omega}}} \|{\bf u}\|^2\, {\mathrm{d}} {\bf x}$$ is exactly conserved by the scheme. Furthermore, we observe from Fig. 5 that with $$h$$ refinement, the enstrophy at the final time does not seem to converge to its initial value, when $$s=1$$, whereas we do observe convergence for $$s\,{=}\,2$$. This is expected since the suboptimal error estimates provided in Guzmán et al. (2016) are sharp for $$s\,{=}\,1$$, so we can only guarantee convergence of the vorticity in $$L^2$$ for $$s\geq 2$$ via an inverse estimate. Fig. 2. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, centred scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. White colour is used for values outside of range. Fig. 2. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, centred scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. White colour is used for values outside of range. Fig. 3. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind Lie derivative scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 3. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind Lie derivative scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 4. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind GSS scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 4. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind GSS scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 5. View largeDownload slide Enstrophy history $$Z(t)$$ for the centred scheme with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 5. View largeDownload slide Enstrophy history $$Z(t)$$ for the centred scheme with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. The vorticity fields for the upwind GSS and Lie derivative schemes are similar and provide much better agreement to the reference solution in Liu & Shu (2000). The enstrophy time evolution in Fig. 6 indicates that the introduction of upwinding induces dissipation in enstrophy. However, the upwind Lie derivative schemes still preserves energy exactly, whereas the upwind GSS scheme also dissipates energy. Fig. 6. View largeDownload slide Kinetic energy, $$K(t)$$, and enstrophy history, $$Z(t)$$, for the upwind Lie derivative scheme (a) and the upwind GSS scheme (b), with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 6. View largeDownload slide Kinetic energy, $$K(t)$$, and enstrophy history, $$Z(t)$$, for the upwind Lie derivative scheme (a) and the upwind GSS scheme (b), with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. 6. Summary and outlook In this article, we extended the variational discretization proposed in Pavlov et al. (2009) for the incompressible Euler equations to the finite-element setting. We largely based our efforts on the work of Heumann & Hiptmair (2011, 2013) and Heumann et al. (2016), where the authors introduced Galerkin discretizations for the Lie derivative operator. Specifically, we used such operators as discrete Lie algebra variables, so that the discrete equations of motion could be derived directly from an appropriately defined Lagrangian and the Hamilton–d’Alembert’s principle. We found that the discretization obtained using this strategy coincides with the centred flux scheme proposed in Guzmán et al. (2016). Moreover, it has a built-in energy conservation property and satisfies a discrete Kelvin’s circulation theorem. We also introduced an upwind-stabilized version of the algorithm, which preserves both properties. Finally, we provided a convergence analysis proving (suboptimal) convergence rates for the upwind scheme for sufficiently high-order finite elements. Numerical tests suggest that this result might not be sharp, as we obtain optimal convergence rates even using low-order spaces. The methodology that we used in this article to discretize the incompressible Euler equations is general and can be applied to different fluid models that share a similar variational structure. The simplest example of such possible extensions is given by the Euler-alpha model (see, e.g., Holm et al., 1998). In this case, we keep the setting of the Euler equations described in Section 2.2, whereas the Lagrangian is given by l(u)=12∫Ω(‖u‖2+α2‖gradu‖2)vol, and $$\delta{l}/\delta{\bf u}$$ can be identified with an opportunely defined momentum $$m = A {\bf u} \cdot {\mathrm{d}}{\bf x}$$, with $$A = (I-\alpha^2{\it {\Delta}})$$ and $$I$$ being the identity map. Suppose that we have a discretization for $$A$$ given by $$A_h$$. Then, in the notation of Section 3, we can take $$A_h:{\bf W}_h \rightarrow {\bf W}_h$$ to be a linear invertible operator acting on the velocity finite-element space $${\bf W}_h$$. If we assume $$A_h$$ to be symmetric, the variational derivation we discussed in Section 3 can be applied without changes so that eventually we get an advection equation for the momentum $$m$$ expressed in terms of our discrete Lie derivative operator. Many MHD and GFD models can also be derived from a variational principle similar to the one that applies for perfect incompressible fluids (see, e.g., Holm et al., 1998). In these cases, one needs to add advected quantities to the system, which again require an appropriate definition for a discrete Lie derivative. Moreover, the variational derivation of Section 3 needs to be appropriately modified, although we expect to maintain its general features. This provides a direction for our future research. Footnotes 1The discrete Lie derivatives proposed in Heumann & Hiptmair (2011) are more general than considered in this article, since they account for a larger class of finite-element spaces, either conforming or not. In the language of differential forms, we use spaces that are conforming with respect to the $$L^2$$-adjoint of the exterior derivative; this case is considered explicitly in Section 4.1.1 of Heumann (2011). References Abraham R. , Marsden J. E. & Ratiu T. ( 1988 ) Manifolds, Tensor Analysis, and Applications , vol. 75 , 2nd edn . New York : Springer . Google Scholar CrossRef Search ADS Arnold D. N. , Falk R. S. & Winther R. ( 2006 ) Finite element exterior calculus, homological techniques, and applications. Acta Numer. , 15 1 – 155 . Google Scholar CrossRef Search ADS Arnol’d V. ( 1966 ) Sur la g éom étrie diff érentielle des groupes de Lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits. Ann. Inst. Fourier , 16 319 – 361 . Google Scholar CrossRef Search ADS Balay S. , Gropp W. D. , McInnes L. C. & Smith B. F. ( 1997 ) Efficient management of parallelism in object oriented numerical software libraries. Modern Software Tools in Scientific Computing ( Arge E. Bruaset A. M. & Langtangen H. P. eds). Boston, MA : Birkhäuser Press , pp. 163 – 202 . Google Scholar CrossRef Search ADS Balay S. , Abhyankar S. , Adams M. F. , Brown J. , Brune P. , Buschelman K. , Dalcin L. , Eijkhout V. , Gropp W. D. , Kaushik D. , Knepley M. G. , McInnes L. C. , Rupp K. , Smith B. F. , Zampini S. & Zhang H. ( 2015 ) PETSc users manual. Technical Report ANL-95/11—Revision 3.6. Lemont, IL : Argonne National Laboratory . Boffi D. , Fortin M. & Brezzi F. ( 2013 ) Mixed Finite Element Methods and Applications . Springer Series in Computational Mathematics. Berlin, Heidelberg : Springer . Google Scholar CrossRef Search ADS Brenier Y. ( 1989 ) The least action principle and the related concept of generalized flows for incompressible perfect fluids. J. Amer. Math. Soc. , 2 225 – 255 . Google Scholar CrossRef Search ADS Brenier Y. ( 2000 ) Derivation of the Euler equations from a caricature of coulomb interaction. Commun. Math. Phys. , 212 93 – 104 . Google Scholar CrossRef Search ADS Brezzi F. , Marini L. D. & Süli E. ( 2004 ) Discontinuous Galerkin methods for first-order hyperbolic problems. Math. Models Methods Appl. Sci. , 14 1893 – 1903 . Google Scholar CrossRef Search ADS Butler G. & Rogers T. ( 1971 ) A generalization of a lemma of Bihari and applications to pointwise estimates for integral equations. J. Math. Anal. Appl. , 33 77 – 81 . Google Scholar CrossRef Search ADS Cockburn B. , Kanschat G. & Schötzau D. ( 2005 ) A locally conservative LDG method for the incompressible Navier–Stokes equations. Math. Comp. , 74 1067 – 1095 . Google Scholar CrossRef Search ADS Dalcin L. D. , Paz R. R. , Kler P. A. & Cosimo A. ( 2011 ) Parallel distributed computing using Python. Adv. Water Resour. , 34 1124 – 1139 . Google Scholar CrossRef Search ADS Desbrun M. , Hirani A. N. , Leok M. & Marsden J. E. ( 2005 ) Discrete exterior calculus. arXiv preprint math/0508341 . Desbrun M. , Gawlik E. S. , Gay-Balmaz F. & Zeitlin V. ( 2014 ) Variational discretization for rotating stratified fluids. Discrete Contin. Dyn. Syst. , 34 477 – 509 . Ebin D. G. & Marsden J. E. ( 1970 ) Groups of diffeomorphisms and the motion of an incompressible fluid. Ann. Math. , 92 102 – 163 . Google Scholar CrossRef Search ADS Gallouët T. & Mérigot Q. ( 2017 ) A Lagrangian scheme à la Brenier for the incompressible Euler equations . Accepted for publication in Foundations of Computational Mathematics . Guzmán J. , Shu C.-W. & Sequeira F. A. ( 2016 ) H (div) conforming and DG methods for incompressible Euler’s equations. IMA J. Numer. Anal. Available at https://doi.org/10.1093/imanum/drw054. Hendrickson B. & Leland R. ( 1993 ) A Multilevel Algorithm for Partitioning Graphs . Tech. Rep. SAND93-1301. Albuquerque, NM : Sandia National Laboratories . Heumann H. ( 2011 ) Eulerian and semi-Lagrangian methods for advection-diffusion of differential forms. Ph.D. Thesis , ETH Zürich , Switzerland . Heumann H. & Hiptmair R. ( 2011 ) Eulerian and semi-Lagrangian methods for advection-diffusion for differential forms. Discrete Contin. Dyn. Syst. , 29 1471 – 1495 . Google Scholar CrossRef Search ADS Heumann H. & Hiptmair R. ( 2013 ) Stabilized Galerkin methods for magnetic advection. ESAIM Math. Model. Numer. Anal. , 47 1713 – 1732 . Google Scholar CrossRef Search ADS Heumann H. , Hiptmair R. & Pagliantini C. ( 2016 ) Stabilized Galerkin for transient advection of differential forms. Discrete Contin. Dyn. Syst. Ser. S , 9 185 – 214 . Holm D. D. , Marsden J. E. & Ratiu T. S. ( 1998 ) The Euler-Poincaré equations and semidirect products with applications to continuum theories. Adv. Math. , 137 1 – 81 . Google Scholar CrossRef Search ADS Liu J.-G. & Shu C.-W. ( 2000 ) A high-order discontinuous Galerkin method for 2D incompressible flows. J. Comput. Phys. , 160 577 – 596 . Google Scholar CrossRef Search ADS Marsden J. E. & Ratiu T. S. ( 2010 ) Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems . New York : Springer Science & Business Media . McLachlan R. I. ( 1993 ) Explicit Lie–Poisson integration and the Euler equations. Phys. Rev. Lett. , 71 3043 – 3046 . Google Scholar CrossRef Search ADS PubMed Mérigot Q. & Mirebeau J.-M. ( 2016 ) Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport. SIAM J. Numer. Anal. , 54 3465 – 3492 . Google Scholar CrossRef Search ADS Mullen P. , McKenzie A. , Pavlov D. , Durant L. , Tong Y. , Kanso E. , Marsden J. E. & Desbrun M. ( 2011 ) Discrete Lie advection of differential forms. Found. Comput. Math. , 11 131 – 149 . Google Scholar CrossRef Search ADS Pavlov D. , Mullen P. , Tong Y. , Kanso E. , Marsden J. E. & Desbrun M. ( 2009 ) Structure-preserving discretization of incompressible fluids. Phys. D , 240 443 – 458 . Google Scholar CrossRef Search ADS Rathgeber F. , Ham D. A. , Mitchell L. , Lange M. , Luporini F. , Mcrae A. T. T. , Bercea G.-T. , Markall G. R. & Kelly P. H. J. ( 2016 ) Firedrake: automating the finite element method by composing abstractions. ACM Trans. Math. Software , 43 , 24:1 – 24:27 . Google Scholar CrossRef Search ADS Zeitlin V. ( 1991 ) Finite-mode analogs of 2D ideal hydrodynamics: Coadjoint orbits and local canonical structure. Phys. D: Non. Phen. , 49.3 , 353 – 362 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

A variational $$\boldsymbol{H}({\rm div})$$ finite-element discretization approach for perfect incompressible fluids

Loading next page...
 
/lp/ou_press/a-variational-boldsymbol-h-rm-div-finite-element-discretization-bdm65VtREr
Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx033
Publisher site
See Article on Publisher Site

Abstract

Abstract We propose a finite-element discretization approach for the incompressible Euler equations which mimics their geometric structure and their variational derivation. In particular, we derive a finite-element method that arises from a nonholonomic variational principle and an appropriately defined Lagrangian, where finite-element $$\boldsymbol{H}({\rm div})$$ vector fields are identified with advection operators; this is the first successful extension of the structure-preserving discretization of Pavlov et al. (2009) to the finite-element setting. The resulting algorithm coincides with the energy-conserving scheme proposed by Guzmán et al. (2016). Through the variational derivation, we discover that it also satisfies a discrete analogous of Kelvin’s circulation theorem. Further, we propose an upwind-stabilized version of the scheme that dissipates enstrophy while preserving energy conservation and the discrete Kelvin’s theorem. We prove error estimates for this version of the scheme, and we study its behaviour through numerical tests. 1. Introduction The equations of motion for perfect fluids with constant density, namely the incompressible Euler equations, possess a rich geometric structure that is directly responsible for their many properties. More precisely, one can derive these equations from a variational principle using geometric techniques that also apply to classical mechanical systems such as the rigid body. This geometric view of perfect fluids is due to Arnol’d (1966). It has been extended to a wide variety of fluid models including magnetohydrodynamics (MHD) and geophysical fluid dynamics (GFD) models through the incorporation of advected quantities (see Holm et al., 1998). In this article, we derive finite-element numerical methods for perfect incompressible fluids, aiming to preserve as much of this geometric structure as possible. We pursue this objective for two main reasons. On the one hand, by addressing structure preservation in numerical simulations, we expect to preserve features of the solution that are usually lost by standard discretization approaches. On the other hand, preserving the structure of the equations at the discrete level enables us to devise discretizations with similar properties for different fluid models in a unified and systematic way. The idea of structure preservation in numerical simulations of the Euler equations is not new. Probably, the most successful approach in this context in terms of number of conservation laws is the sine truncation due to Zeitlin (1991), further developed in the form of a geometric integrator by McLachlan (1993). This uses a modified spectral truncation of the original system, yielding a Lie–Poisson structure of the equations at the discrete level. Such a truncation, however, is only well defined in the very specific case of a two-dimensional periodic domain, and generalizations are far from trivial. One of the most recent approaches is the one associated to the school of Brenier (1989, 2000), see, e.g., the discretizations proposed in Gallouët & Mérigot (2017) and Mérigot & Mirebeau (2016). Such methods are based on the reformulation of the variational principle underlying the incompressible Euler equations into an optimal transport problem. This leads to the definition of the notion of generalized incompressible flow, in which the motion of the fluid particles is defined in a probabilistic sense. The methods in Gallouët & Mérigot (2017) and Mérigot & Mirebeau (2016) aim at producing approximations of the flow map in this setting. In particular, they construct minimizing geodesics between discrete measure-preserving maps transporting the particle mass distributions. Moreover, they establish convergence to classical solutions of the Euler equations. In this article, we follow a different route, which is inspired by the work of Pavlov et al. (2009). In this case, the discrete flow maps are discretized by considering their action on piecewise constant scalar functions intended as right composition. However, such maps are not constructed explicitly, and the final algorithm is fully Eulerian, i.e., one approximates the evolution of the velocity field at fixed positions in space, without following each fluid particle. One key ingredient to achieve this is the identification of vector fields with discrete Lie derivatives, i.e., advection operators, also acting on piecewise constant scalar functions. Because of this identification, one can reduce the variational principle governing the dynamics from the space of discrete flow maps to the one of discrete velocity fields, i.e., from a Lagrangian (material) description to an Eulerian (spatial) one. The final algorithm can be interpreted as an advection problem for the velocity one-form by means of another discrete Lie derivative operator, defined consistently with the discrete exterior calculus (DEC) formalism (see Desbrun et al., 2005). Desbrun et al. (2014) extended this promising idea to generate a variational numerical scheme for rotating stratified fluid models. Both the algorithms presented by Pavlov et al. (2009) and Desbrun et al. (2014) can be classified as low-order finite difference methods and require a dual grid construction to generate an appropriate discrete Lagrangian. We will be concerned with generalizing the approach presented by Pavlov et al. (2009) using finite-element methods, with the aim of producing higher order variational schemes that require a single computational grid. As we are interested in discretizing incompressible fluids, we need finite-element vector spaces that encode the divergence-free constraint in their construction. This problem is generally solved by means of mixed finite-element methods (see, e.g., Boffi et al., 2013) that effectively enable us to build such spaces as subspaces of $$\boldsymbol{H}(\mathrm{div})$$ (the Hilbert space of square-integrable vector fields with square-integrable divergence on a given domain). The mathematical structure of mixed problems is particularly clear when formulating them in terms of differential forms, as Arnold et al. (2006) did by introducing the finite-element exterior calculus (FEEC) framework. In FEEC, many different mixed finite-element spaces can be seen as specific instances of polynomial differential form spaces. One can analyse different discretizations at once; most importantly, stability can be addressed in a systematic way independently of the domain topology. From our perspective, using differential forms in the discretization will be useful to have a single definition for a discrete Lie derivative, both when acting on scalar functions and differential one-forms (which can be identified with vector fields). Mullen et al. (2011) defined a numerical approach in this direction in the context of DEC using semi-Lagrangian techniques. Heumann & Hiptmair (2011) developed similar techniques formulated as Galerkin methods using the differential form spaces arising from FEEC. In the same work, they defined Eulerian-type discretizations including upwind stabilization, which they applied to linear advection/advection–diffusion problems with Lipschitz-continuous advecting velocity (see also Heumann & Hiptmair, 2013, for an application to the magnetic advection problem). Recently, these techniques were further developed to allow for piecewise Lipschitz-continuous advecting velocities (Heumann et al., 2016). Our strategy will be to adopt this latter approach and apply it to the case where the vector field generating the Lie derivative is in $$\boldsymbol{H}({\rm div})$$ finite-element spaces. In this way, we will be able to reproduce the structure of the variational algorithm in the work of Pavlov et al. (2009), transplanting the main features to the finite-element setting. This will lead us to rediscover the centred flux discretization described in Guzmán et al. (2016), and inspired by the discretization for the Navier–Stokes equations proposed in Cockburn et al. (2005). Our derivation shows that such a discretization, besides conserving kinetic energy, also inherits a discrete version of Kelvin’s circulation theorem. As a matter of fact, variational formulations lead to conservation laws resulting from symmetries in Hamilton’s principle, as a consequence of Noether’s theorem. Unfortunately, just as in Pavlov et al. (2009), this approach does not lead to conservation of enstrophy. In turn, this yields suboptimal convergence rates as observed in Guzmán et al. (2016), because a priori control of enstrophy is required for stability of smooth solutions. To overcome these issues, we introduce an upwind-stabilized version of the algorithm, which is different from the one studied in Guzmán et al. (2016). More specifically, we construct the method in such a way that it still possesses the same properties as the centred flux discretization, even though it cannot be derived from first principles. We prove suboptimal-order $$s$$ convergence when polynomial spaces of order $$s$$ are used, for $$s\,{>}\,1\,{+}\,n/2$$, where $$n$$ is the domain dimension. However, numerical results suggest that, at least in two dimensions, the scheme converges optimally with rate $$s+1$$ for $$s\geq 1$$. This article is organized as follows. In Section 2, we introduce the notions of one-forms and Lie derivatives, and we use these to review the variational structure underlying the incompressible Euler equations. In Section 3, we describe our discretization approach. In particular, we derive the central and upwind discretizations and we describe their properties. In Section 4, we present a convergence proof for the upwind scheme, using as starting point the error estimate for the centred scheme derived in Guzmán et al. (2016). We provide some numerical tests in Section 5. Finally, in Section 6, we conclude by describing further possible developments of our approach. 2. Preliminaries This section provides the basic tools necessary to formulate our numerical approach. We proceed in two steps. First, in Section 2.1, we introduce notation for one-forms and Lie derivatives. Then, in Section 2.2, we review the continuous geometric formulation of perfect incompressible fluids, and show how the equations of motion can be formulated as an advection problem for the momentum one-form. We refer to Abraham et al. (1988) and Marsden & Ratiu (2010) for more details on these topics. 2.1 One-forms and Lie derivatives Let $$V$$ be a vector space. A one-form on $$V$$ is an element of $$V^*$$ the dual of $$V$$, i.e., a linear functional $$a:V\,{\rightarrow}\, \mathbb{R}$$. For any $${\bf v}\,{\in}\, V$$, we denote by $$a(\bf v)\,{\in}\, \mathbb{R}$$ the pairing of $$a$$ with $$\bf v$$. Consider now an $$m$$-dimensional manifold $${\it {\Omega}}\subset \mathbb{R}^n$$. At each point $$p\in {\it {\Omega}}$$, we have a vector space $$T_p{\it {\Omega}}$$, namely the tangent space, whose elements are the tangent vectors to $${\it {\Omega}}$$ at $$p$$. The union of such spaces is the tangent bundle of $${\it {\Omega}}$$, $$T{\it {\Omega}} := \bigcup_{p\in {\it {\Omega}}} T_p{\it {\Omega}}$$. A vector field on $${\it {\Omega}}$$ is a smooth assignment $${\bf v}: {\it {\Omega}} \rightarrow T{\it {\Omega}}$$ such that $${\bf v}_p$$ is tangent to $${\it {\Omega}}$$ at $$p$$, i.e., $${\bf v}_p \in T_p{\it {\Omega}}$$. Differential one-forms generalize the concept of one-form as dual elements of vector fields on $${\it {\Omega}}$$. Specifically, we define the cotangent space $$T^*_p{\it {\Omega}}$$ to be the dual of $$T_p{\it {\Omega}}$$ and introduce the cotangent bundle $$T^*{\it {\Omega}} := \bigcup_{p\in {\it {\Omega}}} T^*_p{\it {\Omega}}$$. Then, a differential one-form is an assignment $$a:{\it {\Omega}} \rightarrow T^*{\it {\Omega}}$$ such that $$p \in {\it {\Omega}} \rightarrow a_p\in T^*_p{\it {\Omega}}$$. For any smooth vector field $$\bf v$$, the pairing $$a({\bf v})$$ is the real-valued function $$a({\bf v}): p \in {\it {\Omega}} \rightarrow a_p({\bf v}_p)\in \mathbb{R}$$. Then, we define the total pairing of a differential one-form $$a$$ with a vector field $$\bf v$$ to be ⟨a,v⟩:=∫Ωa(v)vol, (2.1) with $$\mathrm{vol}$$ being the measure on $${\it {\Omega}}$$ induced by the Euclidean metric. In this article, we restrict our study to the special case when $${\it {\Omega}}$$ is an $$n$$-dimensional domain in $$\mathbb{R}^n$$, with $$n=2,3$$, and admits global Cartesian coordinates $$\{x_i\}_{i=1}^n$$ oriented accordingly to an orthonormal basis $$\{{\bf e}_i\}_{i=1}^n$$ of $$\mathbb{R}^n$$. Then, a vector field $$\bf v$$ on $${\it {\Omega}}$$ has the coordinate representation $${\bf v} = \sum_{i=1}^n v_i {\bf e}_i$$, where $$v_i$$ are real-valued functions on $${\it {\Omega}}$$. Note that $${\bf e}_i$$, for $$i=1,\ldots,n$$, is interpreted as a constant vector field on $${\it {\Omega}}$$. Now, let $$\{{\mathrm{d}} x_i\}_{i=1}^n$$ be the dual basis to $$\{{\bf e}_i\}_{i=1}^n$$, i.e., $${\mathrm{d}} x_i({\bf e}_j) =\delta_{ij}$$. Then, a differential one-form $$a$$ has the coordinate representation $$a = \sum_{i=1}^n a_i\,{\mathrm{d}} x_i$$, where $$a_i$$ are real-valued functions on $${\it {\Omega}}$$, and $$a({\bf v}) = \sum_{i=1}^n a_i v_i$$. Note that $${{\mathrm{d}} x}_i$$, for $$i=1,\ldots,n$$, is interpreted as a constant differential one-form on $${\it {\Omega}}$$. In this simple setting, we have a trivial identification of a differential one-form $$a$$ with a vector field $$\bf a$$ by the equation a=a⋅dx, (2.2) where the right-hand side represents the formal inner product between the vector $${\bf a} = (a_1,\ldots,a_n)$$ and $${\mathrm{d}} {\bf x} = ({\mathrm{d}} x_1, \ldots, {\mathrm{d}} x_n)$$. Let $${\it {\Gamma}}\subset {\it {\Omega}}$$ be a smooth curve, parameterized by the functions $$x_i = x_i(s)$$, for $$i=1,\ldots,n$$ and $$s\in I$$, where $$I$$ is an open interval of $$\mathbb{R}$$ with $$0\in I$$. The quantity $${\bf s} = \sum_{i=1}^n{\mathrm{d}} x_i/{\mathrm{d}} s|_{s=0} \,{\bf e}_i$$ is a vector field on $${\it {\Gamma}}$$, since for each $$p\in {\it {\Gamma}}$$, $${\bf s}_p\in T_p{\it {\Gamma}}$$. Then, the integral of a differential one-form on $${\it {\Gamma}}$$ is defined by ∫Γa:=∫Ia(s)ds, (2.3) which is independent of the parameterization of $${\it {\Gamma}}$$. A smooth vector field $$\bf v$$ on $${\it {\Omega}}$$ can be identified with a one-parameter family of local diffeomorphisms $$\varphi: I \times U \rightarrow {\it {\Omega}}$$, with $$U$$ being an open set of $${\it {\Omega}}$$ and $$I\subset \mathbb{R}$$ an open set around $$0\in \mathbb{R}$$, such that for all $$t\in I$$, $$\varphi_t$$ is a local diffeomorphism defined by ∂φt∂t=v∘φt, (2.4) and $$\varphi_0 = e$$ the identity map on $${\it {\Omega}}$$. The map $$\varphi$$ defines the flow of the vector field $$\bf v$$. The Lie derivative of a differential one-form $$a$$ with respect to the vector field $$\bf v$$ is the one-form $${\sf L}_{\bf v} a$$ satisfying the equation ∫ΓLva=ddt|t=0∫φt∘Γa, (2.5) for any smooth curve $${\it {\Gamma}} \subset {\it {\Omega}}$$. Moreover, if $$n=3$$, we have the identification Lva=(grad(v⋅a)+(curla)×v)⋅dx, (2.6) or if $$n=2$$, Lva=(grad(v⋅a)+(rota)v⊥)⋅dx, (2.7) where at each $$p\in{\it {\Omega}}$$, $${\bf{v}}_p^\perp$$ is the vector $$\bf v_p$$ after a clockwise rotation of $$\pi/2$$, and $$\mathrm{rot}\,{\bf{a}}:=\mathrm{div}\,{\bf{a}}^\perp$$ (see Abraham et al. (1988) for more details). The Lie derivative can be regarded as a generalization of the scalar advection operator. To realize this, we need to find an equivalent of equation (2.5) for scalar functions. For a given scalar function $${b}:{\it {\Omega}}\rightarrow \mathbb{R}$$, we can interpret pointwise evaluation as the natural analogue of the integration of differential one-forms over curves (see, e.g., Heumann & Hiptmair, 2011). With this substitution, equation (2.5) becomes Lvb=ddt|t=0b∘φt=v⋅gradb, (2.8) which is just the directional derivative of $${b}$$ in the direction of $$\bf v$$. We conclude this section by introducing some notation for the function spaces and relative inner products we will use in this article. We denote by $$L^2({\it {\Omega}})$$ (respectively, $$\boldsymbol{L}^2({\it {\Omega}})$$) the space of square-integrable functions (respectively, vector fields) on $${\it {\Omega}}$$. We denote by $$(\cdot, \cdot)_{{\it {\Omega}}}$$ and $$\|\cdot\|_{{\it {\Omega}}}$$ the inner product and associated norm for both $$L^2({\it {\Omega}})$$ and $$\boldsymbol{L}^2({\it {\Omega}})$$. The definitions extend to one-forms, i.e., for any one-form $$a={\bf a}\cdot{{\mathrm{d}}} {\bf x}$$ and $$b={\bf b}\cdot{{\mathrm{d}}} {\bf x}$$, (a,b)Ω:=⟨a,b⟩=(a,b)Ω. (2.9) Moreover, we denote by $$L^p({\it {\Omega}})$$, with $$1\,{\leq}\, p \,{\leq}\, \infty$$, the standard generalization of $$L^2({\it {\Omega}})$$ with norm $$\|\cdot\|_{L^p({\it {\Omega}})}$$. The spaces $$H^k({\it {\Omega}})$$ (respectively, $$W^{k,p}({\it {\Omega}})$$) are the spaces of scalar functions with derivatives up to order $$k>0$$ in $$L^2({\it {\Omega}})$$ (respectively, $$L^p({\it {\Omega}})$$); the relative norms are $$\|\cdot\|_{H^k({\it {\Omega}})}$$ and $$\|\cdot\|_{W^{k,p}({\it {\Omega}})}$$. Similar definitions hold for vector fields; for example, $$\boldsymbol{L}^p({\it {\Omega}})$$ is the space of vector fields with components in $$L^p({\it {\Omega}})$$. Finally, we define H(div,Ω):={v∈L2(Ω):divv∈L2(Ω)}. (2.10) 2.2 Geometric formulation of incompressible perfect fluids In this section, we review the formal variational derivation of the incompressible Euler equations. This will be important in Section 3.2, where we develop our discretization approach by repeating this derivation step by step in a finite-dimensional setting. The diffeomorphisms on $${\it {\Omega}}$$ play a major role in the description of fluids, as they can be used to represent the motion of particles in the absence of discontinuities such as fractures, cavitations or shocks. More specifically, denote by $$\mathrm{Diff}({\it {\Omega}})$$ the group of (boundary-preserving) diffeomorphisms from $${\it {\Omega}}$$ to itself, with group product given by map composition. The subgroup of volume-preserving diffeomorphisms $$\mathrm{Diff}_\mathrm{vol} ({\it {\Omega}})$$ is defined by Diffvol(Ω):={g∈Diff(Ω):det(Dg)=1}, (2.11) where $$Dg$$ denotes the Jacobian matrix of the map $$g$$. This definition ensures that for any $$g \in \mathrm{Diff}_\mathrm{vol} ({\it {\Omega}})$$ and any open set $$U\subset {\it {\Omega}}$$, $$g(U)$$ has the same measure as $$U$$. Then, the motion of an incompressible fluid is represented by a curve on $$\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$, i.e., a one-parameter family of diffeomorphisms $$\{g_t\}_{t\in I}$$, where $$I$$ is an open interval of $$\mathbb{R}$$ containing $$0\in \mathbb{R}$$, so that for all $$t\in I$$, $$g_t \in \mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ and $$g_0 = e$$ the identity map on $${\it {\Omega}}$$. This means that, for each point $$X \in {\it {\Omega}}$$, $$g_t(X)$$ denotes the position occupied at time $$t$$ by the particle that was at $$X$$ when $$t=0$$. We define the material velocity field, ∂gt∂t(X)=U(t,X). (2.12) This is a vector field on $${\it {\Omega}}$$ that gives the velocity at time $$t\in I$$ of the particle that occupied the position $$X\in {\it {\Omega}}$$ at time $$t=0$$, for each $$X\in {\it {\Omega}}$$. The spatial velocity is the vector field u(t,x):=U(t,gt−1(x)). (2.13) This gives the velocity of the particle occupying the position $$x\in {\it {\Omega}}$$ at time $$t$$. Therefore, in the material formulation, $$X\in {\it {\Omega}}$$ is to be considered as a particle label, whereas in the spatial formulation, $$x\in {\it {\Omega}}$$ is a fixed position in the physical space. By the fact that $$g_t\in \mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$, we obtain that any spatial velocity field is tangent to the boundary $$\partial {\it {\Omega}}$$ and has zero divergence. We also see that the spatial velocity is left unchanged when composing on the right $$g_t$$ with any time independent diffeomorphism $$q \in \mathrm{Diff}({\it {\Omega}})$$, i.e., ∂∂t(gt∘q)∘(gt∘q)−1=∂∂tgt∘gt−1. (2.14) The interpretation of this is simple: composing the map $$q$$ with $$g_t$$ on the right is the same as relabelling the particles in the domain $${\it {\Omega}}$$. On the other hand, the spatial velocity field gives the velocity at fixed positions in space and it is not dependent on the specific particle labels. The geometric picture of incompressible perfect fluids is based on interpreting $$G:= \mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ as a configuration space, then applying the reduction techniques of dynamical systems on Lie groups (see, e.g., Marsden & Ratiu, 2010). In particular, $$\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ is not a Lie group, and one should provide a rigorous characterization of $$\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ to proceed in such a direction (see, e.g., Ebin & Marsden, 1970). Here, we are not concerned with these details, and we will proceed formally. A map $$g \in G$$ characterizes the configuration of the system by assigning to each particle label $$X\in {\it {\Omega}}$$ its physical position in the domain $$x = g(X)$$. Let $${\bf U} \in T_g G$$. Then there exists a one-parameter family $$\{g_t\}_{t\in I}$$, defined as above, such that for a fixed $$\bar{t}\in I$$, $$g_{\bar t} = g$$ and $${\frac{\partial}{\partial {t}}}|_{t=\bar{t}}\, g_t = {\bf U}$$. By comparison with equation (2.12), $$\bf U$$ can be interpreted as a material velocity field on $${\it {\Omega}}$$. In the classical Lagrangian setting, the evolution of the system is described by a curve on $$TG$$. In our case, this would be equivalent to the material formulation, where we follow the motion as a time-dependent diffeomorphism $$g_t$$ together with the associated material velocity field as defined in equation (2.12). The tangent space at the identity $$\mathfrak{g}:= T_eG$$ is the Lie algebra of $$G$$. By equation (2.13), we can regard elements of $$\mathfrak{g}$$ as spatial velocity fields on $${\it {\Omega}}$$. Therefore, an element of $$\mathfrak{g}= T_e G$$ can be identified with a vector field on $${\it {\Omega}}$$ tangent to its boundary $$\partial {\it {\Omega}}$$ and with zero divergence. For all $${\bf u},{\bf v}\in \mathfrak{g}$$, we denote by $$({\bf u},{\bf v}) \mapsto \mathrm{ad}_{\bf u} {\bf v}$$ their Lie algebra product. Let $$q_s$$ be a curve on $$G$$ such that $$q_0 = 0$$ and $$\dot{q}_0 = \bf v$$, the map $$\mathrm{ad}:\mathfrak{g} \times \mathfrak{g} \rightarrow \mathfrak{g}$$ is defined by aduv:=ddt|t=0dds|s=0gt∘qs∘gt−1=ddt|t=0(Dgtv)∘gt−1=−[u,v], (2.15) where $$[{\bf u}, {\bf v}]$$ is the Jacobi–Lie bracket, defined in coordinates by [u,v]:=∑i,j=1n(uj∂vi∂xj−vj∂ui∂xj)ei. (2.16) The Lagrangian of the system is a function $$L:TG\rightarrow \mathbb{R}$$ and is given by the kinetic energy of the spatial velocity field, i.e., L(g,g˙)=12∫Ω‖g˙∘g−1‖2vol. (2.17) Such a Lagrangian is right invariant because of equation (2.14). Therefore, the dynamics can be expressed in terms of the reduced Lagrangian $$l({\bf u})\,{:=}\, L(e, {\bf u})$$ only. More precisely, Hamilton’s principle for the Lagrangian in equation (2.17) is equivalent to δ∫t0t1l(u)dt=0, (2.18) with variations $$\delta {\bf u} = \dot{{\bf v}} - \mathrm{ad}_{{\bf u}} {\bf v}$$, where $${\bf v}_t \in \mathfrak{g}$$ is an arbitrary curve and $${\bf v}_{t_0}\,{=}\, {\bf v}_{t_1}\,{=}\,0$$. The variational principle in equation (2.18) produces the Euler–Poincaré equations, ⟨ddtδlδu,v⟩+⟨δlδu,aduv⟩=0, (2.19) for all $${\bf v}\in \mathfrak{g}$$, where we define ⟨δlδu,v⟩:=ddϵ|ϵ=0l(u+ϵv). (2.20) Now $$\delta l/\delta \bf u$$ can be identified with a one-form $$m$$ on $${\it {\Omega}}$$, which we will refer to as momentum, imposing that ⟨δlδu,v⟩=∫Ωm(v)vol, (2.21) for all $${\bf v}\in \mathfrak{g}$$. In our case, we readily see that $${m} = u := {\bf u}\cdot{{\mathrm{d}} {\bf x}}$$; however, we will distinguish between momentum one-forms $${m}$$ and velocity one-forms $$u$$, as this will be useful when describing model problems characterized by more complex Lagrangians. Inserting equations (2.15) and (2.21) into equation (2.19), we obtain ∫Ωm˙(v)vol−∫Ωm([u,v])vol=0. (2.22) If $$n=3$$, for divergence free vector fields $$\bf u$$ and $$\bf v$$ we have the useful identity $$[{\bf u}, {\bf v}] = - {\bf{curl}}({\bf u} \times {\bf v})$$. Hence, integrating by parts, the second integral can be rewritten as follows, ∫Ωm([u,v])vol =−∫Ωm⋅curl(u×v)vol =∫Ω(u×curlm)⋅vvol =−∫ΩLum(v)vol. (2.23) The same result holds for $$n=2$$. Note that the $${\bf{grad}}$$ term in equation (2.6) or (2.7) vanishes upon pairing with $$\bf v$$ because $$\mathrm{div}\,{\bf v} =0$$. Reinserting equation (2.23) into equation (2.22) yields the incompressible Euler equations in variational form ⟨m˙,v⟩+⟨Lum,v⟩=0, (2.24) for all $${\bf v} \in \mathfrak{g}$$, where $$m=u:= \bf u \cdot {{\mathrm{d}} {\bf x}}$$. Remark 2.1 Equation (2.23) establishes a link between the bracket on $$\mathfrak{g}$$ coming from Lagrangian reduction and the advection of the momentum one-form $${m}$$ by the Lie derivative operator. In the following we will develop a discretization approach that maintains this property, and in this way provides a discrete version of both the bracket and the Lie derivative. 3. Finite-element variational discretization In this section, we construct a variational algorithm for the incompressible Euler equations using discrete Lie derivatives intended as finite-element operators as the main ingredient. As starting point for this, we describe the Galerkin discretization of Lie derivatives proposed by Heumann & Hiptmair (2011) (see also Heumann & Hiptmair, 2013), in the case where the advecting velocity is a smooth vector field. Then, in Section 3.1, we consider a generalization for advecting velocities in finite-element spaces (meaning that the velocities have reduced smoothness). Such a generalization coincides with the one proposed in Heumann et al. (2016) for the more general case of piecewise Lipschitz-continuous advecting velocities. Next, in Section 3.2, we construct a variational algorithm using these operators as discrete Lie algebra variables, adapting the approach of Pavlov et al. (2009) to the finite-element setting. Moreover, we show that such an algorithm coincides with the centred flux discretization described in Guzmán et al. (2016). In Section 3.3, we include upwind stabilization to the method derived in the previous section. Finally, in Section 3.4, we show that both the centred flux and upwind schemes conserve energy and possess a discrete version of Kelvin’s circulation theorem. Throughout this section, we assume $$n=3$$, although all the results of this section extend without major modifications to the case $$n=2$$. 3.1 Discrete Lie derivatives Let $${\mathscr T}_h$$ be a triangulation on $${\it {\Omega}}$$, i.e., a decomposition of $${\it {\Omega}}$$ into simplices $$K\,{\subset}\, {\it {\Omega}}$$, with $$h$$ being the maximum simplex diameter. We define $$V_h\subset L^2({\it {\Omega}})$$ to be a chosen scalar finite-element space on $${\mathscr T}_h$$, and $${\bf W}_h \subset\boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ to be a chosen $$\boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$-conforming finite-element space of vector fields on $${\it {\Omega}}$$. We define $$W_h:= {\bf W}_h\cdot{{\mathrm{d}} {\bf x}}$$, i.e., $$W_h$$ is the set of one-forms $$a = {\bf a} \cdot {\mathrm{d}} {\bf x}$$ with $${\bf a} \in {\bf W}_h$$. More specifically, for $$r\geq 0$$, we denote by $$V^r_h\subset L^2({\it {\Omega}})$$ a scalar finite-element space such that, for all elements $$K\in {\mathscr T}_h$$, $$V^r_h|_K = \mathscr P_{r}(K),$$ where $$\mathscr P_r(K)$$ is the set of polynomial functions on $$K$$ of degree up to $$r$$. In other words, $$V_h^r$$ is the standard scalar discontinuous Galerkin space of order $$r$$. Analogously, we denote by $${\bf W}^r_h \subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ either the Raviart–Thomas space of order $$r$$, defined for $$r\geq 0$$ by RTr(Th):={u∈H(div,Ω):u|K∈(Pr(K))n+xPr(K)}, (3.1) or the Brezzi–Douglas–Marini space of order $$r$$, and defined for $$r\geq 1$$ by BDMr(Th):={u∈H(div,Ω):u|K∈(Pr(K))n}. (3.2) Later, we will use the divergence-free subspace of $${\bf W}^r_h$$, which is the same for $${\bf{RT}}_r(\mathscr T_h)$$ and $${\bf{BDM}}_r(\mathscr T_h)$$ (see, e.g., Boffi et al., 2013, Corollary 2.3.1). Let $$\mathscr F$$ be the sets of $$(n-1)$$-dimensional simplices in $${\mathscr T}_h$$, and let $${\mathscr F}^\circ\subset {\mathscr F}$$ be the set of facets $$f\in {\mathscr F}$$ such that $$f \cap \partial {\it {\Omega}} = \emptyset$$. We fix an orientation for all $$f\in {\mathscr F}$$ by specifying the unit normal vector $${\bf n}_f$$ in such a way that on the boundary $${\bf n}_f$$ is equal to the outward pointing normal $${\bf n}_{\partial {\it {\Omega}}}$$. An orientation on the facet $$f\in {\mathscr F}$$ defines a positive and negative side on $$f$$, so that any $$a\in V_h$$ has two possible restrictions on $$f$$, denoted $$a^+$$ and $$a^-$$, with $$a^+:= a|_{K^+}$$ and $$K^+\in \mathscr T_h$$ being the element with outward normal $${\bf{n}}_f$$. Therefore, we can define for each $$a \in V_h$$ and on each $$f\in \mathscr F^\circ$$ the jump operator $$[a] := a^+ - a^-$$ and the average operator $$\{a\} := (a^+ + a^-)/2$$. Such definitions extend without modifications to vector fields $${\bf a} \in {\bf W}_h$$ or to their one-form representation $$a = {\bf a}\cdot {{\mathrm{d}} {\bf x}}\in W_h$$. Lemma 3.1 Let $${\bf W}_h\subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$, then for all $${\bf a}\in {\bf W}_h$$ and for all $$f\in \mathscr F^\circ$$, $${\bf a}\cdot {\bf n}_f$$ is single valued on $$f$$, and in particular $$[{\bf a}]\cdot {\bf n}_f = 0$$. Proof. See Lemma 5.1 in Arnold et al. (2006). □ For all elements $$K\,{\in}\,\mathscr T_h$$ and $$f\,{\in}\,\mathscr F$$, we denote by $$(\cdot,\cdot)_K$$ and $$(\cdot,\cdot)_f$$ the standard $$L^2$$ inner products on $$K$$ and $$f$$. We use the same definition for both scalar and vector quantities, and we extend it to one-forms as in equation (2.9). Furthermore, given a smooth vector field $$\boldsymbol \beta$$, for any $$f\,{\in}\,\mathscr F$$ and for any smooth scalar functions $$a$$ and $$b$$ on $$f$$ we define (a,b)f,β:=(β⋅nfa,b)f. (3.3) Similarly, for any element $$K\in\mathscr T_h$$ and for any smooth scalar functions $$a$$ and $$b$$ on $$\partial K$$ define (a,b)∂K,β:=(β⋅n∂Ka,b)∂K, (3.4) where $${\bf n}_{\partial K}$$ is the unit normal to the boundary of the element $$K$$ pointing outwards. The same definitions hold if $$a$$ and $$b$$ are one-forms or vector fields with smooth components on $$f$$ or $$\partial K$$. A discrete Lie derivative in this article means any finite-element operator from $$V_h$$ to itself, or from $$W_h$$ to itself, which is also a consistent discretization of the Lie derivative. In Heumann & Hiptmair (2011), the authors proposed two types of discrete Lie derivatives1 based on either an Eulerian or a semi-Lagrangian approach, which were used in the context of linear advection and advection–diffusion problems. We will take as starting point their Eulerian discretization (see Section 4 in Heumann & Hiptmair, 2011), where the advecting velocity is a fixed smooth vector field $$\boldsymbol \beta$$. In the following, we will consider extensions of this approach to the case, where $$\boldsymbol \beta$$ is not fully continuous; specifically $$\boldsymbol \beta \in {\bf W}_h\subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$. The resulting operators coincide with the ones proposed in Heumann et al. (2016) for piecewise Lipschitz-continuous advecting velocities. Such a discussion will be useful to discretize the nonlinear advection term in the Euler equations. The Eulerian-type discrete Lie derivative proposed by Heumann & Hiptmair (2011) is given by the following definition. Definition 3.2 Let $$\boldsymbol{\beta}$$ be a smooth vector field on $${\it {\Omega}}$$ tangent to the boundary $$\partial {\it {\Omega}}$$. The finite-element operator $${\sf L}_{\boldsymbol{\beta}}^h: V_h \rightarrow V_h$$ (respectively, $${\sf L}_{\boldsymbol{\beta}}^h: {W}_h \rightarrow {W}_h$$) is defined by (Lβha,b)Ω=∑K(Lβa,b)K+∑f∈F∘(−([a],{b})f,β+([a],[b])f,cfβ), (3.5) where $$c_f:f \rightarrow \mathbb{R}$$ is a scalar function depending on $$\boldsymbol \beta$$, for all $$a,b\in V_h$$ (respectively, $$a,b \in W_h$$). The consistency of such a discretization is proved by observing that if the advected quantity $$a$$ is smooth, $${\sf L}^h_{\boldsymbol{\beta}} a$$ coincides with the Galerkin projection of $${\sf L}_{\boldsymbol{\beta}}a$$ onto $$V_h$$ or $$W_h$$. Note that the Lie derivative discretization on $$V_h$$, as defined in equation (3.5), coincides with the classical discontinuous Galerkin discretization of the advection operator described in Brezzi et al. (2004). This can be verified using integration by parts in equation (3.5), yielding the following equivalent definition for the operator $${\sf L}_{\boldsymbol{\beta}}^h: V_h \rightarrow V_h$$, (Lβha,b)Ω=−∑K(a,div(βb))K+∑f∈F∘(({a},[b])f,β+([a],[b])f,cfβ), (3.6) for all $$a,b\in V_h$$. Similarly, the operator $${\sf L}_{\boldsymbol{\beta}}^h: W_h \rightarrow W_h$$ can be equivalently defined by (Lβha,b)Ω=∑K(a,curl(β×b)−βdivb)K+∑f∈F∘(({a},[b])f,β+([a],[b])f,cfβ), (3.7) for all $$a,b\in W_h$$ with $$a = {\bf a}\cdot {{\mathrm{d}}}{\bf x}$$ and $$b = {\bf b} \cdot {{\mathrm{d}}{\bf x}}$$. Details on the calculations leading to equations (3.6) and (3.7) can be found in Heumann (2011). Remark 3.3 There are two choices for the function $$c_f$$ that are particularly important: $$c_f = 0$$ (centred discretization). In this case, equation (3.6) reduces to (Lβha,b)Ω=−∑K(a,div(βb))K+∑f∈F∘({a},[b])f,β, (3.8) for all $$a,b\in V_h$$. We refer to such discretization as centred, since the facet integrals contain the average of the advected quantity $$a$$; $$c_f = \boldsymbol{\beta} \cdot{\bf n}_f /(2|\boldsymbol{\beta} \cdot{\bf n}_f|)$$ (upwind discretization). In this case, equation (3.6) reduces to (Lβha,b)Ω=−∑K(a,div(βb))K+∑f∈F∘(aup,[b])f,β, (3.9) for all $$a,b\in V_h$$, where $$a^{\rm up} = a^+$$ if $$\boldsymbol{\beta}\cdot {\bf n}_f\geq0$$ and $$a^{\rm up} = a^-$$ if $$\boldsymbol{\beta}\cdot {\bf n}_f<0$$, i.e., in the facet integrals $$a$$ is always evaluated from the upwind side. A similar interpretation holds when considering the action of $${\sf L}_{\boldsymbol{\beta}}^h$$ on $$W_h$$. Lemma 3.4 Let $$n=3$$ and let $$\boldsymbol{\beta}$$ be a smooth vector field on $${\it {\Omega}}$$ tangent to the boundary $$\partial {\it {\Omega}}$$. Then, the discrete Lie derivative $${\sf L}^h_{\boldsymbol{\beta}}:W_h\rightarrow W_h$$ in equation (3.5) can be equivalently defined by (Lβha,b)Ω= ∑K(a,curl(β×b)−βdivb)K +∑f∈F∘((nf×{a},β×[b])f+(cfnf×[a],β×[b])f), (3.10) for all $$a,b\in W_h$$ with $$a = {\bf a}\cdot {{\mathrm{d}}}{\bf x}$$ and $$b = {\bf b} \cdot {{\mathrm{d}}{\bf x}}$$. Proof. The proof relies on the fact that if $${\bf a} \in {\bf W}_h\subset\boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ then, by Lemma 3.1, $${\bf a}\cdot {\bf n}_f$$ is single-valued across elements of $$\mathscr T_h$$, so $$[{\bf a}]\cdot {\bf n}_f = 0$$ on all facets $$f\in \mathscr F^\circ$$. Then, the result follows by applying to the facet integrals in equation (3.7) the vector identity (A×B)⋅(C×D)=(A⋅C)(B⋅D)−(B⋅C)(A⋅D), (3.11) for all $${\bf{A}},{\bf{B}},{\bf{C}},{\bf{D}}\in \mathbb{R}^3$$. See Heumann (2011) for details. □ As we are concerned with discretizing the Euler equations, we want to extend Definition 3.2 to the case where the advecting velocity also belongs to a finite-element space. We will limit our study to the case $$\boldsymbol{\beta} \in {\bf W}_h\subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$. This will prove to be enough for our scope since if $$\boldsymbol{\beta} \in \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$, we can ensure that $$\mathrm{div}\, \boldsymbol{\beta} = 0$$. Remark 3.5 Equation (3.5) is still well defined if $$\boldsymbol{\beta} \in {\bf W}_h$$, since in this case, by Lemma 3.1, $$\boldsymbol{\beta} \cdot {\bf n}_f$$ is single valued on each facet $$f \in {\mathscr F}^\circ$$. In view of Remark 3.5, equation (3.5) can be used without modification to extend the definition of the discrete Lie derivative to the case $$\boldsymbol{\beta} \in {\bf W}_h$$. However, this is not the only possible extension: there exist many other discretization approaches that reduce to Definition 3.2 for continuous advecting velocity. In particular, Lemma 3.4 suggests the following alternative definition, which is also the one adopted in Heumann et al. (2016). Definition 3.6 Let $$n=3$$ and let $$\boldsymbol{\beta}\in{\bf W}_h$$ such that $$\boldsymbol{\beta} \cdot{\bf n}_{\partial {\it {\Omega}}} =0$$ on $$\partial {\it {\Omega}}$$. The finite-element operator $${\sf X}_{\boldsymbol{\beta}}^h:V_h \rightarrow V_h$$ is defined by $${\sf X}_{\boldsymbol{\beta}}^h \,a := {\sf L}_{\boldsymbol{\beta}}^h \, a$$ for all $$a \in V_h$$, whereas the operator $${\sf X}_{\boldsymbol{\beta}}^h:W_h \rightarrow W_h$$ is defined by (Xβha,b)Ω= ∑K(a,curl(β×b)−βdivb)K +∑f∈F∘((nf×{a},[β×b])f+(cfnf×[a],[β×b])f), (3.12) where $$c_f:f \rightarrow \mathbb{R}$$ is a scalar function depending on $$\boldsymbol \beta$$, for all $$a,b\in W_h$$ with $$a ={\bf a} \cdot {\mathrm{d}} {\bf x}$$ and $$b ={\bf b} \cdot {\mathrm{d}} {\bf x}$$. In other words, $${\sf X}_{\boldsymbol{\beta}}^h$$ coincides with $${\sf L}^h_{\boldsymbol{\beta}}$$ when applied to scalar functions, but not when applied to one-forms. Proposition 3.7 For a continuous vector field $$\boldsymbol \beta$$, we have the equivalence $${\sf X}_{\boldsymbol{\beta}}^h={\sf L}_{\boldsymbol{\beta}}^h$$. Therefore, $${\sf X}_{\boldsymbol{\beta}}^h$$ is a discrete Lie derivative, i.e., it is a consistent discretization of the Lie derivative operator. Proof. This is immediate by comparing equation (3.10) with equation (3.12) and observing that if $$\boldsymbol \beta$$ is continuous, for any $${\bf b}\in {\bf W}_h$$, we have $$\boldsymbol{\beta} \times [{\bf b}] = [\boldsymbol{\beta} \times {\bf b}]$$ on all facets $$f\in {\mathscr F}^\circ$$. □ 3.2 Discrete incompressible Euler equations We now use the discrete Lie derivatives defined in the previous section to design a variational discretization approach for the incompressible Euler equations. The derivation presented here closely follows Pavlov et al. (2009), but is adapted to the context of finite-element discretizations and does not require a dual grid. First, we will define a discrete group that approximates the configuration space $$\mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$. Next, we will derive the relative Euler–Poincaré equations for an appropriately chosen Lagrangian. Finally, we will connect the resulting algorithm to the Lie derivative discretizations in Definition 3.6. Such an approach will lead us to rediscover the centred flux discretization proposed in Guzmán et al. (2016). The main issue related to deriving a variational integrator for the incompressible perfect fluids is that the configuration space $$G = \mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$ is infinite dimensional, so one needs to find an appropriate finite-dimensional approximation that converges, in the limit, to the original system. In Pavlov et al. (2009), this issue is solved by applying Koopman’s lemma, i.e., identifying group elements with their action, intended as right composition, on $$L^2$$ functions on $${\it {\Omega}}$$. This is to say that $$G$$ can be equivalently represented as a subgroup $$G(V) \subset GL (V)$$ of the group of invertible linear maps from $$V:= L^2({\it {\Omega}})$$ to itself, by means of a group homomorphism $$\rho:G \rightarrow G(V)$$ defined by ρ(g)⋅a=ρg⋅a:=a∘g−1 (3.13) for any $$g\,{\in}\,G$$ and $$a\,{\in}\,C^\infty({\it {\Omega}})\,{\subset}\,V$$. As the Lie algebra $$\mathfrak g$$ of $$G$$ is the set of divergence-free vector fields $$\bf v$$ tangent to the boundary, the Lie algebra $$\mathfrak g(V)$$ of $$G(V)$$ is the set of Lie derivatives $${\sf L}_{\bf v}: V \rightarrow V$$, interpreted as unbounded operators, with respect to such vector fields. As a matter of fact, we have dds|s=0ρgs⋅a=dds|s=0(a∘gs−1)=−Lva (3.14) for any $$a\in C^\infty({\it {\Omega}}) \subset V$$, and where $$g_s$$ is a curve on $$G$$ such that $$g_0=e$$ and $$\dot{g_s}|_{s=0} = {\bf v}$$. To get a discrete version of this picture, we start by restricting $$V$$ to be the finite-element space $$V_h$$. Therefore, we consider the finite-dimensional Lie group $$GL(V_h)$$ and we seek an appropriate subgroup $$G_h(V_h)$$ that in the limit approximates $$G(V)$$. Clearly, this cannot be accomplished by restricting $$G(V)$$ to $$V_h$$, as we have done for the general linear group $$GL(V)$$, because it would imply losing the group structure. However, one can still construct a subgroup $$G_h(V_h)\subset GL(V_h)$$ that approximates $$G(V)$$ and hence $$G$$. For all $$g^h\in GL(V_h)$$ and $$a\in V_h$$, we denote by $$a \mapsto g^h a$$ the action of $$g^h$$ on $$a$$. Then, $$G_h(V_h)$$ can be defined as follows. Definition 3.8 The Lie group $$G_h(V_h)\subset GL(V_h)$$ is the subgroup of the group of linear invertible operators from $$V_h$$ to itself, defined by the following properties: for all $$g^h\in G_h(V_h)$$,  (gha,ghb)Ω=(a,b)Ω∀a,b∈Vh, (3.15) ghc=c∀c∈R, (3.16) where $$c$$ is considered as a constant function on $${\it {\Omega}}$$. It is trivial to check that both equations (3.15) and (3.16) are verified in the continuous case. In particular for any $$g\in \mathrm{Diff}({\it {\Omega}})$$, we have $$c \circ g^{-1} = c$$. Moreover, as a consequence of volume preservation, for all $$a,b \in V$$ and $$g\in \mathrm{Diff}_{\mathrm{vol}}({\it {\Omega}})$$, ∫Ω(a∘g−1)(b∘g−1)vol=∫Ω(ab)∘g−1vol=∫Ωdet(Dg−1)abvol=∫Ωabvol (3.17) since $$g^{-1}({\it {\Omega}}) = {\it {\Omega}}$$. The Lie algebra of $$GL(V_h)$$ is the set $$\mathfrak{gl}(V_h)$$ of all linear operators from $$V_h$$ to itself. For all $${\sf A}^h\in \mathfrak{gl}(V_h)$$ and $$a\in V_h$$, we denote by $$a \mapsto {\sf A}^h a$$ the action of $${\sf A}^h$$ on $$a$$. Then, the Lie algebra of $$G_h(V_h)$$ can be defined as follows. Definition 3.9 The Lie algebra $$\mathfrak{g}_h(V_h)\subset \mathfrak{gl}(V_h)$$ associated to the Lie group $$G_h(V_h)\subset GL(V_h)$$ is the subalgebra of the linear operators from $$V_h$$ to itself defined by the following properties: for all $${\sf A}^h \in \mathfrak{g}_h(V_h)$$,  (Aha,b)Ω=−(a,Ahb)Ω∀a,b∈Vh, (3.18) Ahc=0∀c∈R. (3.19) Note that equations (3.18) and (3.19) can be derived directly from equations (3.15) and (3.16). In particular, equation (3.18) can be again related to volume preservation. Moreover, we can regard $$\mathfrak{g}_h(V_h)$$ as an approximation $$\mathfrak{g}(V)$$ so its elements are to be thought as discrete Lie derivatives in view of equation (3.14). The collection of all the spaces $$G_h(V_h)$$, for $$h$$ arbitrarily small, is not $$G(V)$$, even if the union of the spaces $$V_h$$ is dense in $$V$$. This is because there exist elements of $$GL(V)$$ that satisfy equations (3.15) and (3.16), but do not belong to $$G(V)$$ (for example, they might be representative of maps that are not even continuous). Therefore, before defining the discrete system we must restrict ourselves to a space smaller than $$G_h(V_h)$$. The same considerations hold at the Lie algebra level, meaning that one needs to constrain the dynamics in a subset $$S_h(V_h) \,{\subset}\, \mathfrak g_h(V_h)$$. Our definition for such a space is based on the discrete Lie derivatives of Definition 3.6. In particular, we fix a finite-element space $${\bf W}_h \subset \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$ and define W∘h:={u∈Wh:divu=0,u⋅n∂Ω=0}. (3.20) We denote by $$\overset{\circ}{W}_{h}:= \overset{\circ}{\mathbf W}_{h}\cdot {{\mathrm{d}} {\bf x}}$$, i.e., $$\overset{\circ}{W}_{h}$$ is the set of one-forms $$u = {\bf u} \cdot {\mathrm{d}} {\bf x}$$ with $${\bf u} \in \overset{\circ}{\mathbf W}_{h}$$. Analogous definitions hold for $$\overset{\circ}{\mathbf W}_{h}^r$$ and $$\overset{\circ}{W}_{h}^r$$. Note that $$\overset{\circ}{\mathbf W}_{h}^r$$ is the same space if constructed using $${\bf W}_h^r = {\bf{RT}}_r(\mathscr T_h)$$ or $${\bf W}_h^r = {\bf{BDM}}_r(\mathscr T_h)$$ (see, e.g., Boffi et al., 2013, Corollary 2.3.1). The set of advection operators associated to $$\overset{\circ}{\mathbf W}_{h}$$ is given by Sh(Vh):={Xuh:Vh→Vhwithcf=0:u∈W∘h}. (3.21) More specifically, we define Shs(Vhr):={Xuh:Vhr→Vhrwithcf=0:u∈W∘hs}. (3.22) Note that setting $$c_f=0,$$ we immediately ensure that $$S_h(V_h)$$ is a linear space; this can be directly verified using Definition 3.6. Lemma 3.10 $$S_h(V_h) \subset \mathfrak{g}_h(V_h)$$. Proof. We need to verify that the elements of $$S_h(V_h)$$ satisfy equations (3.18) and (3.19). By Definition 3.6 and equation (3.5), for all $$a,b\in V_h$$, (Xuha,b)Ω=∑K(Lua,b)K−∑f∈F∘([a],{b})f,u. (3.23) The Lie derivative is just the directional derivative when applied to scalar functions, as shown in equation (2.8). Moreover, since $$\mathrm{div}\, {\bf u}=0$$, integration by parts yields ∑K(Lua,b)K =−∑K(a,Lub)K+∑K(a,b)∂K,u =−∑K(a,Lub)K+∑f∈F∘((a+,b+)f,u−(a−,b−)f,u). (3.24) Note that the decomposition of the integrals on $$\partial K$$ is justified by the fact that each facet is shared by the boundary of two adjacent elements, therefore it appears twice in the sum, but with two opposite normals (since $${\bf n}_{\partial K}$$ is always oriented outwards). Moreover, we have the identity (a+,b+)f,u−(a−,b−)f,u=({a},[b])f,u+([a],{b})f,u. (3.25) Combining equations (3.23), (3.24) and (3.25) gives (Xuha,b)Ω=−(a,Xuhb)Ω, (3.26) which is the same as in equation (3.18). Equation (3.19) is immediately verified by observing that if we set $$a=c$$ in equation (3.23), with $$c \in \mathbb{R}$$ being a constant function on $${\it {\Omega}}$$, we obtain (Xuhc,b)Ω=0, (3.27) for all $$b\in V_h$$, and therefore $${\sf X}_{\bf u}^h c =0$$. □ Lemma 3.11 Let $${\it {\Omega}}$$ be a domain in $$\mathbb{R}^n$$ and let $$\mathscr T_h$$ be a triangulation on $${\it {\Omega}}$$. For any $$r\geq s$$ and $$r\geq1$$, the spaces $$S_h^s(V^r_h)$$ and $$\overset{\circ}{\mathbf W}_{h}^s$$ are isomorphic. Proof. We assumed that $$\{{\bf e}_i\}_{i=1}^n$$ is a global orthonormal reference frame on $$\mathbb{R}^n$$ with coordinates $$\{x_i\}_{i=1}^n$$. Then, since $$r\geq 1$$, $$x_i\in V_h^r$$ and we can define the map $$\hat{\cdot}: {\mathfrak g}_h(V_h^r) \rightarrow (V_h^r)^n$$ by A^h:=∑i=1n(Ahxi)ei, (3.28) for all $${\sf A}^h \in {\mathfrak g}_h$$. By a standard result of mixed finite-element theory (see, e.g., Boffi et al., 2013, Corollary 2.3.1), for any $${\bf u} \in \overset{\circ}{\mathbf W}_{h}^s$$ we have that $${\bf u}|_K \in ({\mathscr P}_s(K))^n$$ for any $$K\in {\mathscr T}_h$$. Then, since $$r \geq s$$, for all $${\bf u} \in \overset{\circ}{\mathbf W}_{h}^s$$, $$u_i = {\bf u} \cdot {\bf e}_i \in V_h^r$$. Moreover, we can verify by direct calculation, i.e., applying Definition 3.6, that for all $${\sf X}^h_{\bf u}\in S_h^s(V_h^r)$$, $${\sf X}^h_{\bf u}\, x_i = u_i$$. Hence the restriction of the $$\hat{\cdot}$$ map to $$S_h^s(V^r_h)$$ is the map $${\sf X}^h_{\bf u} \mapsto {\bf u}$$ which is a surjection from $$S_h^s(V_h^r)$$ to $$\overset{\circ}{\mathbf W}_{h}^s$$. On the other hand, the map $${\bf u} \mapsto {\sf X}^h_{\bf u}$$ from $$\overset{\circ}{\mathbf W}_{h}^s$$ to $$S_h^s(V_h^r)$$ is also a surjection by definition of $$S_h^s(V_h^r)$$, hence the result. □ The different spaces introduced so far can be collected in the diagrams represented in Fig. 1. The definition of the $$\hat{\cdot}$$ map suggests the choice of a Lagrangian $$l:\mathfrak{g}_h(V_h)\rightarrow \mathbb{R}$$ given by Fig. 1. View largeDownload slide Diagrams representing the discretization of the group $$G:=\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ and its Lie algebra. The dashed arrows represent relations that hold only in an approximate sense. Fig. 1. View largeDownload slide Diagrams representing the discretization of the group $$G:=\mathrm{Diff}_\mathrm{vol}({\it {\Omega}})$$ and its Lie algebra. The dashed arrows represent relations that hold only in an approximate sense. l(Ah):=12(A^h,A^h)Ω, (3.29) so that the restriction of $$l$$ to $$S_h(V_h)$$ is just the total kinetic energy, i.e., $$l({\sf X}_{\bf u}^h) = \|{\bf u}\|^2_{{\it {\Omega}}}/2$$. Such a choice for the Lagrangian was already suggested in Pavlov et al. (2009) as an alternative to the Lagrangian construction proposed therein based on dual grids. To get a geometric picture analogous to the continuous case, we need a Lagrangian defined on the whole tangent space $$TG_h(V_h)$$ and approximating equation (2.17). We achieve this by extending $$l$$ by right translations to get $$L:TG_h(V_h)\rightarrow \mathbb{R}$$ defined by L(gh,g˙h):=l(g˙h(gh)−1)=12∫Ω‖g˙h(gh)−1^‖2vol, (3.30) which is right invariant by construction. Analogously, the space $$S_h(V_h)$$ can be extended by right translation to give a right invariant sub-bundle of $$TG_h(V_h)$$. To see how this is done, consider the right translation map $$R_{g^h} : q^h \in G_h(V_h) \rightarrow q^hg^h \in G_h(V_h)$$ and denote by $$T_eR_{g^h} : {\sf A}^h \in {\mathfrak g}_h(V_h) \rightarrow {\sf A}^h g^h \in T_{g^h}G_h(V_h)$$ its tangent map at the identity (see Marsden & Ratiu, 2010, for more details). Then, the right invariant extension of $$S_h(V_h)$$ is obtained by collecting the spaces $$T_eR_{g^h} S_h(V_h)$$ for all $$g^h\,{\in}\, G_h(V_h)$$, where $$T_eR_{g^h} S_h(V_h) \,{\subset} T_{g^h} G_h(V_h)$$ is the space of linear operators which can be written as the composition $${\sf A}^h g^h$$, for a given $${\sf A}^h\in {\mathfrak g}_h(V_h)$$. Lemma 3.12 The Euler–Poincaré–d’Alembert equations relative to the right invariant Lagrangian $$L$$ with constraint $$\dot{g}^h\in T_eR_{g^h} S_h(V_h)$$ are given by: Find $${\sf A}^h\in S_h(V_h)$$ such that (ddtA^h,B^h)Ω+(A^h,[Ah,Bh^])Ω=0, (3.31) for all $${\sf B}^h\in S_h(V_h)$$, where $$[{{\sf A}^h, {\sf B}^h}]:= {\sf A}^h {\sf B}^h - {{\sf B}^h {\sf A}^h}$$ is the commutator of the linear operators $${\sf A}^h$$ and $${\sf B}^h$$. Proof. The proof is analogous to that of Theorem 1 in Pavlov et al. (2009). In essence, it is just an application of the reduction theorem (Marsden & Ratiu, 2010, Theorem 13.5.3). In particular, note that proceeding as equation (2.15), we get $$\mathrm{ad}_{{\sf A}^h} {\sf B}^h = [{{\sf A}^h, {\sf B}^h}]$$. Then, equation (3.31) is the finite-dimensional analogue of equation (2.19). The only difference is that we added a nonholonomic constraint on the velocity so that both the solution and the variations are constrained in the space $$S_h(V_h)$$. □ Theorem 3.13 Under the hypotheses of Lemma 3.11, the Euler–Poincaré–d’Alembert equations (3.31) relative to the right invariant Lagrangian in equation (3.30) with $$S_h(V_h) = S^s_h(V^r_h)$$ are equivalent to: Find $$u \in \overset{\circ}{W}_{h}^s$$ such that (u˙,v)Ω+(Xuhu,v)Ω=0, (3.32) with $$u = {\bf u} \cdot {{\mathrm{d}} {\bf x}}$$, for all $${v} \in \overset{\circ}{W}_{h}^s$$. Proof. First, by Lemma 3.11 for each $${\sf A}^h\in S_h^s(V_h^r)$$ and $${\sf B}^h\in S_h^s(V_h^r),$$ there is a unique $${\bf u}\in \overset{\circ}{\mathbf W}_{h}^s$$ and $${\bf v}\in \overset{\circ}{\mathbf W}_{h}^s$$ such that $${\sf A}^h=-{\sf X}^h_{\bf u}$$ and $${\sf B}^h=-{\sf X}^h_{\bf v}$$. We now examine the two terms in equation (3.31) separately. For the first term, we have (ddtA^h,B^h)Ω=(u˙,v)Ω, (3.33) by definition of the one-form pairing in equation (2.9). As for the second term, we have (A^h,[Ah,Bh^])Ω =−(X^uh,[Xuh,Xvh^])Ω =−∑i=1n(ui,Xuhvi−Xvhui)Ω =−∑i=1n∑K(ui,Luvi−Lvui)K+∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v) =−∑K(u,[u,v])K+∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v), (3.34) where we used the identity $$\sum_{i=1}^n({\sf L}_{\bf u}v_i - {\sf L}_{\bf v} u_i) {\bf e}_i = [{\bf u}, {\bf v}]$$, which is a consequence of equations (2.8) and (2.16). As for the facet integrals, we first notice that ∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v)=∑f∈F∘(({u},[v])f,u−({u},[u])f,v). (3.35) Moreover, since $$\bf u$$ has continuous normal component on the mesh facets, we can substitute ({u},[v])f,u=({u},[v])f,{u}=(nf×{u},{u}×[v])f, (3.36) where for the last equality we have used the same reasoning as in the proof of Lemma 3.4. Similarly, for the second term in equation (3.35), we have ({u},[u])f,v=(nf×{u},{v}×[u])f=−(nf×{u},[u]×{v})f. (3.37) Therefore, equation (3.35) becomes ∑i=1n∑f∈F∘(({ui},[vi])f,u−({ui},[ui])f,v) =∑f∈F∘(nf×{u},{u}×[v]+[u]×{v})f =∑f∈F∘(nf×{u},[u×v])f. (3.38) The last equality in equation (3.38) follows from the linearity of the cross product, i.e., {u}×[v]+[u]×{v} =(u++u−)/2×(v+−v−)+(u+−u−)×(v++v−)/2 =u+×v+−u−×v−=[u×v]. (3.39) Inserting equation (3.38) into equation (3.34) yields (A^h,[Ah,Bh^])Ω=∑K(u,curl(u×v))K+∑f∈F∘(nf×{u},[u×v])f, (3.40) where we used the identity $$[{\bf u},{\bf v}] = -{\bf{curl}}({\bf u}\times {\bf v})$$ valid for any divergence-free vector fields $${\bf u}$$ and $${\bf v}$$. Recalling that $$c_f=0$$, comparison with equation (3.12) concludes the proof. □ Proposition 3.14 The centred discretization in equation (3.32), i.e., the case $$c_f\,{=}\,0$$, on a simplicial triangulation coincides with the discontinuous Galerkin discretization with centred fluxes proposed in Guzmán et al. (2016), which can be written as follows: Find $$({\bf u},p) \,{\in}\, {\bf W}_h^s \times V_h^k$$, where $$k=s$$ if $${\bf W}_h^s = {\bf{RT}}_s(\mathscr T_h)$$, or $$k=s-1$$ if $${\bf W}_h^s = {\bf{BDM}}_s(\mathscr T_h)$$, such that {(u˙,v)Ω−∑K(u,u⋅∇v)K+∑f∈F∘(u⋅nf{u},[v])f−∑K(p,divv)K=0∑K(divu,q)K=0  (3.41) for all $$({\bf v}, q) \in {\bf W}_h^s \times V_h^k$$, with $${\bf u} \cdot {\bf n}_{\partial {\it {\Omega}}} = 0$$ on $$\partial {\it {\Omega}}$$, and $$\int_{{\it {\Omega}}} p \, {\mathrm{d}} x = 0$$. Proof. First, we observe that restricting equation (3.41) on the subspace $$\overset{\circ}{\mathbf W}_{h}^s$$ yields the equivalent problem: Find $${\bf u} \in \overset{\circ}{\mathbf W}_{h}^s$$ such that (u˙,v)Ω−∑K(u,u⋅∇v)K+∑f∈F∘(u⋅nf{u},[v])f=0, (3.42) for all $${\bf v} \in \overset{\circ}{\mathbf W}_{h}^s$$. On the other hand, we can write equation (3.32) in the following form (u˙,v)Ω+∑K(u,curl(u×v))K+∑f∈F∘(nf×{u},[u×v])f=0. (3.43) We rewrite the volume integrals in equation (3.43) as follows, (u,curl(u×v))K=(u,−u⋅∇v+v⋅∇u)K=−(u,u⋅∇v)K+(∇u2/2,v)K, (3.44) for any $$K \in \mathscr {T}_h$$. We rewrite the facet integrals in equation (3.43) as follows, ∑f∈F∘(nf×{u},[u×v])f =∑f∈F∘((u⋅nf{u},[v])f−(v⋅nf{u},[u])f) =∑f∈F∘((u⋅nf{u},[v])f−(v⋅nf,[u2/2])f) =∑f∈F∘(u⋅nf{u},[v])f−∑K∫Kdiv(v(u2/2))dx =∑f∈F∘(u⋅nf{u},[v])f−∑K(v,∇(u2/2))K, (3.45) where we used equation (3.38) for the equality in the first line. Inserting equations (3.44) and (3.45) into equation (3.43) gives the equivalence of the two algorithms. □ 3.3 Upwind scheme In this section, we study the upwind version of the scheme in equation (3.32). In particular, we extend the discussion of the previous section to the upwind case, with the aim of preserving the main features of the variational derivation. Including upwinding in the variational framework described above is not straightforward. This is because, if we define $$S_h(V_h)$$ by Sh(Vh):={Xuh:Vh→Vhwithcf=c¯f(u):u∈W∘h}, (3.46) where $$\bar{c}_f$$ is a given function of $$\bf u$$, then $$S_h(V_h)$$ ceases to be a linear space, and furthermore $$S_h(V_h)\nsubseteq \mathfrak{g}_h(V_h)$$. This last issue can be solved easily by choosing a group larger than $$G_h(V_h)$$ as configuration space. More specifically, we need to allow for discrete diffeomorphisms, which do not satisfy equation (3.15). We call the new configuration space $$\tilde{G}_h(V_h)$$ and we define it as follows. Definition 3.15 The Lie group $$\tilde{G}_h(V_h)\subset GL(V_h)$$ is the subgroup of the group of linear invertible operators from $$V_h$$ to itself, defined by the following property: for all $$g^h\in \tilde{G}_h(V_h)$$, ghc=c∀c∈R, (3.47) where $$c$$ is considered as a constant function on $${\it {\Omega}}$$. Denote with $$\tilde{\mathfrak{g}}_h(V_h)$$ the Lie algebra of $$\tilde{G}_h(V_h)$$. Then, by the same arguments as in the last section, it is easy to verify that, for any choice of $$c_f$$ in the definition of $$S_h(V_h)$$, $$S_h(V_h)\subset \tilde{\mathfrak{g}}_h(V_h)$$. Next, we note that also Lemma 3.11 holds independently of the choice of $$c_f$$ in the definition of $$S_h(V_h)$$, and in particular, for fixed polynomial orders $$r$$ and $$s$$ with $$r\,{\geq}\,1$$ and $$r\,{\geq}\,s$$, if $${\sf A}^h \in S_h^s(V_h^r)$$ then $${\bf u}:= \hat{\sf A}^h \,{\in}\, \overset{\circ}{\mathbf W}_{h}^s$$. This observation leads to the following result. Proposition 3.16 (Upwinding). Under the hypotheses of Lemma 3.11, the problem described by equation (3.31), where $$S_h(V_h) = S_h^s(V_h^r)$$ is defined to be Shs(Vhr):={Xvh:Vhr→Vhrwithcf=c¯f(u):v∈W∘hs}, (3.48) with $${\bf u}\,{:=}\, \hat{\sf A}^h\,{\in}\, \overset{\circ}{\mathbf W}_{h}^s$$ is equivalent to equation (3.32), where $${\sf X}_{\bf u}^h$$ is defined as in Definition 3.6 with $$c_f=-{\bar{c}_f({\bf u})}$$. {In particular, if we choose $$\bar{c}_f({\bf u})=-{{\bf u} \cdot {\bf n}_f}/({2 |{\bf u} \cdot {\bf n}_f|})$$, we obtain the upwind version of equation (3.31). Proof. The proof is analogous to Theorem 3.13. □ Remark 3.17 Note that the upwind version of the scheme in equation (3.32) cannot be derived directly from Hamilton’s principle because the constraint space in equation (3.48) is dependent on the solution of the problem. However, in the next section, we will see that the equivalence between equations (3.32) and (3.31), ensured by Proposition 3.16, is enough to maintain the main properties of the centred discretization (i.e., energy conservation and a discrete version of Kelvin’s circulation theorem) also when upwinding is introduced. 3.4 Energy conservation and discrete Kelvin’s circulation theorem A consequence of using a nonholonomic constraint in our approach is that our discretization does not preserve all the Casimirs, i.e., the conserved quantities, of the original system. This is because the discrete bracket induced by the advection term in equation (3.32) is only a quasi-Poisson bracket (it is only skew symmetric and does not satisfy the Jacobi identity). Nonetheless, because of the variational derivation of the equations of motion, we are still able to ensure energy conservation and derive a discrete version of Kelvin’s circulation theorem, with or without upwinding. Energy conservation is established in the following proposition. Proposition 3.18 (Energy conservation/$$L^2$$ stability) The discrete system in equation (3.32), with or without upwinding, satisfies conservation of the total kinetic energy, i.e., $${\mathrm{d}}{l}/{\mathrm{d}} t=0$$. Proof. This can be verified directly by setting $$v=u$$ in equation (3.32) and exploiting the antisymmetric structure of the equations. □ We now turn to derive a discrete Kelvin’s circulation theorem. At the continuous level, the theorem states that for any closed loop $${\it {\Gamma}}$$, ddt∫gt∘Γu=0, (3.49) where $$g_t$$ is the flow of $${\bf u}$$, with $$u= {\bf u}\cdot {{\mathrm{d}} {\bf x}}$$. To show in what sense equation (3.49) is verified at the discrete level, we will proceed as in Pavlov et al. (2009) by replacing the integral in equation (3.49) by a pairing with currents, i.e., vector fields supported in a small region around a curve. A discrete current is a vector field $${\bf c}\in \overset{\circ}{\mathbf W}_{h}$$ that is nonzero only on a closed loop of elements $$K\,{\in}\,\mathscr T_h$$ and such that the flux of $${\bf c}$$ through adjacent elements of the loop is equal to 1. Clearly, with mesh refinement, a discrete current approaches the notion of closed loop. In other words, if $${\it {\Gamma}}$$ is a closed loop such that there exists a discrete current $${\bf c}$$ with $${\it {\Gamma}}\in \overline{\mathrm{supp}({\bf c})}$$, then ⟨u,c⟩≈∫Γu. (3.50) Let $$g_t$$ be the flow of $$\bf u$$ as defined by equations (2.12) and (2.13). If both $$g_t$$ and a given discrete current $$\bf c$$ are smooth, then $$\bf c$$ is advected by $$g_t$$ via the map $${\bf c} \mapsto (D g_t \, {\bf c}) \circ g_t^{-1}$$, as it can be derived from equation (2.15). Therefore, Kelvin’s theorem can be reformulated using currents, by stating that ddt⟨u,(Dgtc)∘gt−1⟩=0, (3.51) for any smooth current $$\bf c$$. In the discrete setting, $$g_t$$ can only be defined in weak sense, and in general, it is not differentiable, because we only have $${\bf u} \in \boldsymbol{H}(\mathrm{div},{\it {\Omega}})$$. Nonetheless, we can still use the discrete flow map $$g_t^h$$ to reproduce an analogous of equation (3.51), as shown in the following proposition. Proposition 3.19 Let $${\bf u} \in \overset{\circ}{\mathbf W}_{h}$$, with $$u={\bf u}\cdot {{\mathrm{d}} {\bf x}}$$, be a solution of the discrete Euler equations (3.32). Then, for all discrete currents $${\bf c}\in \overset{\circ}{\mathbf W}_{h}$$ ddt|t=0⟨u,(gth)−1Xchgth^⟩=0, (3.52) where $$g_t^h\in G_h(V_h)$$ for the centred discretization, or $$g_t^h\in \tilde{G}_h(V_h)$$ for the upwind discretization, is the discrete flow of $$-{\sf X}_{\bf u}^h$$, i.e., it satisfies $$\dot{g}^h_t =- {\sf X}_{\bf u}^h g^h_t$$ and $$g^h_0 = e$$. The result extends to any time $$t$$ by appropriately translating the definition of $$g_t^h$$ in time. Proof. For the centred discretization, let $$q_t^h\in G_h(V_h)$$ be the discrete flow of $$-{\sf X} ^h_{\bf c}$$ (or $$q_t^h\in \tilde{G}_h(V_h)$$, for the upwind discretization), then [Xuh,Xch]=adXuhXch=−ddt|t=0dds|s=0(gth)−1qshgth=ddt|t=0(gth)−1Xchgth. (3.53) Finally, by equation (3.31), we have ddt|t=0⟨u,(gth)−1Xchgth^⟩=⟨u˙,c⟩|t=0+⟨u,[Xuh,Xch^]⟩|t=0=0. (3.54) □ 4. Convergence analysis In this section, we analyse the convergence properties of the semidiscrete system in equation (3.32). In particular, as the centred discretization, i.e., the case $$c_f =0$$, was studied in Guzmán et al. (2016), we will concentrate on the upwind formulation, i.e., the case $$c_f = {\bf u} \cdot {\bf n}_f/ (2 |{\bf u} \cdot {\bf n}_f|)$$. For simplicity, we continue to restrict to the case $$n=3$$, although the results of this section extend to $$n=2$$. We keep the dependence on $$n$$ explicit when needed. We start with characterizing the space $${\bf W}_h^s$$ and its approximation properties. We will assume that $$\mathscr T_h$$ is quasi-uniform and shape regular and that $${\bf W}_h^s$$ is either $${\bf{RT}}_{s}(\mathscr T_h)$$ or $${\bf{BDM}}_{s}(\mathscr T_h)$$, i.e., the Raviart–Thomas or Brezzi–Douglas–Marini finite-element spaces of order $$s$$ on $${\mathscr T}_h$$. Then, as shown in Boffi et al. (2013) (or in Arnold et al., 2006, in the context of FEEC), there exists a projection operator $${\bf P}_h$$ of smooth vector fields onto $${\bf W}^s_h$$ such that ‖u−Phu‖Ω≤Chr|u|Hr(Ω), (4.1) for any vector field $${\bf u}\in \boldsymbol{H}^{r}({\it {\Omega}})$$ and $$1\leq r \leq s+1$$. From now on, we will denote by $${\bf u}$$ the exact solution of the Euler equation, and by $${\bf u}_h\,{\in}\,\overset{\circ}{\mathbf W}_{h}^s$$ the discrete solution obtained by solving equation (3.32). The estimate for the centred flux scheme in Guzmán et al. (2016) was derived as follows. First, the error is decomposed using the triangle inequality ‖u−uh‖Ω≤‖u−Phu‖Ω+‖Phu−uh‖Ω. (4.2) The first term on the right-hand side is the approximation error, therefore for the convergence estimate is sufficient to obtain a bound on $$\boldsymbol{\gamma}_h:= {\bf P}_h {\bf u} - {\bf u}_h$$. We collect the results of Guzmán et al. (2016) in the following lemma. Lemma 4.1 (Theorem 2.1 in Guzmán et al. (2016)). The error $$\boldsymbol{\gamma}_h$$ satisfies the following inequality ddt‖γh‖Ω2≤C0h2s+C1h2s+2+C2‖γh‖Ω2, (4.3) where $$C_0,C_1,C_2>0$$ are constants depending on $$\|{\bf u}\|_{\boldsymbol{W}^{1,\infty}({\it {\Omega}})}$$, $$\|{\bf u}\|_{\boldsymbol{H}^{r+1}({\it {\Omega}})}$$, $$\|\dot{\bf u}\|_{\boldsymbol{H}^{r+1}({\it {\Omega}})}$$, but independent of $$h$$. Then, if the norms above are uniformly bounded for $$t\in [0,T]$$, the following estimate holds ‖u−uh‖Ω≤Chs, (4.4) for an appropriate constant $$C>0$$ independent of $$h$$, and for all $$t\in [0,T]$$. When upwinding is introduced, i.e., for the case $$c_f ={\bf u} \cdot {\bf n}_f/ (2 |{\bf u} \cdot {\bf n}_f|)$$, we just need to add the upwind contribution to the left-hand side of equation (4.3) and proceed with the estimate. In particular, from the proof of Theorem 2.1 in Guzmán et al. (2016), it can be easily verified that, assuming the exact solution to be continuous, equation (4.3) needs to be replaced by ddt‖γh‖Ω2−∑f∈F∘(cfnf×[uh],[uh×γh])f≤C0h2s+C1h2s+2+C2‖γh‖Ω2. (4.5) Proceeding as in equation (3.34), we rewrite the upwind term as follows, ∑f∈F∘(cfnf×[uh],[uh×γh])f =∑f∈F∘(cfnf×[uh],{uh}×[γh]+[uh]×{γh})f =∑f∈F∘(([uh],[γh])f,cfuh−([uh],[uh])f,cfγh) =∑f∈F∘(([Phu],[γh])f,cfuh−([γh],[γh])f,cfuh−([uh],[uh])f,cfγh). (4.6) Reinserting equation (4.6) into equation (4.5) gives ddt‖γh‖Ω2+∑f∈F∘([γh],[γh])f,cfuh ≤C0h2s+C1h2s+2+C2‖γh‖Ω2 +∑f∈F∘([Phu],[γh])f,cfuh−∑f∈F∘([uh],[uh])f,cfγh, (4.7) where the facet integrals on the left-hand side are always non-negative due to the definition of $$c_f$$. Define I1up:=|∑f∈F∘([Phu],[γh])f,cfuh|,I2up:=|∑f∈F∘([uh],[uh])f,cfγh|, (4.8) then ddt‖γh‖Ω2+∑f∈F∘([γh],[γh])f,cfuh≤C0h2s+C1h2s+2+C2‖γh‖Ω2+I1up+I2up. (4.9) Estimate for$$I_1^{\rm up}$$ Since $$c_f {\bf u}_h\cdot {\bf n}_f = |{\bf u}_h \cdot {\bf n}_f|/2$$, the following estimate holds, I1up=|∑f∈F∘([Phu],[γh])f,cfuh| ≤12|∑f∈F∘([Phu],[Phu])f,cfuh+∑f∈F∘([γh],[γh])f,cfuh| ≤12∑f∈F∘([Phu−u],[Phu−u])f,cfuh+12∑f∈F∘([γh],[γh])f,cfuh ≤Ctrh−1‖uh‖L∞(Ω)‖Phu−u‖Ω2+12∑f∈F∘([γh],[γh])f,cfuh, (4.10) where we used the trace inequality $$\|\cdot\|^2_{\partial K} \leq C_{\rm tr} h^{-1} \|\cdot\|^2_{K}$$, with $$C_{\rm tr}>0$$ and independent of $$h$$. We estimate the $$\boldsymbol{L}^\infty$$ norm of $${\bf u}_h$$ as follows, ‖uh‖L∞(Ω)≤‖u‖L∞(Ω)+Cinvh−n/2(‖u−Phu‖Ω+‖γh‖Ω), (4.11) where we used the inverse inequality $$\|\cdot\|_{\boldsymbol{L}^{\infty}({\it {\Omega}})}\,{\leq}\,C_{\rm inv} h^{-n/2}\|\cdot\|_{\it {\Omega}}$$, with $$C_{\rm inv}\,{>}\,0$$ and independent of $$h$$. Inserting equation (4.11) into equation (4.10) and applying Young’s inequality to the term involving $$\|\boldsymbol{\gamma}_h\|_{\it {\Omega}}$$, we obtain I1up≤C3(h2s+1+h4s+2−n+h3s+2−n/2)+12‖γh‖Ω2+12∑f∈F∘([γh],[γh])f,cfuh, (4.12) for an appropriate constant $$C_3>0$$ independent of $$h$$. Estimate for$$I_2^{\rm up}$$ Proceeding as for $$I_1^{\rm up}$$, we obtain I2up=|∑f∈F∘([uh],[uh])f,cfγh| ≤2|∑f∈F∘([γh],[γh])f,cfγh|+2|∑f∈F∘([Phu−u],[Phu−u])f,cfγh| ≤4CtrCinvh−1−n/2(‖γh‖Ω3+‖γh‖Ω‖Phu−u‖Ω2) ≤C4h−1−n/2‖γh‖Ω3+C4h4s+2−n+12‖γh‖Ω2, (4.13) for an appropriate constant $$C_4>0$$ independent of $$h$$. Final estimate Assuming $$s\,{\geq}\,1$$ we can neglect higher order terms in $$h$$ in the estimates above. In particular, inserting equations (4.12) and (4.13) into equation (4.9), we obtain ddt‖γh‖Ω2−‖γh‖Ω2+∑f∈F∘([γh],[γh])f,cfuh≤C5h2s+C5h−1−n/2‖γh‖Ω3, (4.14) for a constant $$C_5>0$$ independent of $$h$$. We multiply both sides of equation (4.14) by $$e^{-t}$$, and we let $$q(t):= e^{-t}\|\boldsymbol{\gamma}_h\|_{{\it {\Omega}}}^2$$, so we obtain ddtq(t)+e−t∑f∈F∘([γh],[γh])f,cfuh≤e−tC5h2s+et/2C5h−1−n/2q(t)3/2. (4.15) Integrating from $$0$$ to $$T$$ gives q(T) ≤q(0)+(1−e−T)C5h2s+C5h−1−n/2∫0Tes/2q(s)3/2ds ≤C6h2s+C6h−1−n/2∫0Tes/2q(s)3/2ds, (4.16) where $$C_6$$ depends on $$\|{\bf u}|_{t=0}\|_{\boldsymbol{H}^{r+1}({\it {\Omega}})}$$, but not on $$h$$. Then, the nonlinear generalization of the Gronwall’s inequality in Butler & Rogers (1971) yields q(T)≤h2s(1C61/2−2C6hs−1−n/2(eT/2−1))−2, (4.17) with the condition hs−1−n/2<12C63/2(eT/2−1). (4.18) If $$s>1+n/2,$$ such a condition can be verified for any $$T$$, if $$h$$ is sufficiently small. Moreover, for hs−1−n/2<14C63/2(eT/2−1), (4.19) we get q(T)≤4C6h2s, (4.20) which proves the following theorem. Theorem 4.2 Under the same assumptions on $$\bf u$$ as in Lemma 4.1, the discrete solution $${\bf u}_h\in {\bf W}_h^s$$, with $$s>1+n/2$$, obtained by solving the system in equation (3.32) with $$c_f = {\bf u} \cdot {\bf n}_f/ (2 |{\bf u} \cdot {\bf n}_f|)$$, satisfies ‖u−uh‖Ω≤Chs, (4.21) with $$C>0$$ dependent $$T$$, $$\|{\bf u}|_{t=0}\|_{\boldsymbol{H}^{s+1}({\it {\Omega}})}$$, $$\|{\bf u}\|_{H^1([0,T],\boldsymbol{H}^{s+1}({\it {\Omega}}))}$$, $$\|{\bf u}\|_{\boldsymbol{W}^{1,\infty}([0,T]\times {\it {\Omega}})}$$, but not on $$h$$. Remark 4.3 Theorem 4.2 gives a suboptimal rate of convergence for the upwind scheme, i.e., $$c_f\,{=}{{\bf u}_h \cdot {\bf n}_f}/({2|{\bf u}_h \cdot {\bf n}_f|})$$, only for the case $$s>1+n/2$$. This result is rather unsatisfactory, since the numerical tests in the next section suggest that the scheme converges with the optimal rate $$s\,{+}\,1$$ for $$s\,{\geq}\,1$$. Unfortunately, we were not able to prove this result analytically. 5. Numerical tests We now show some numerical results illustrating our analytical results. All the tests presented here are computed in the domain $${\it {\Omega}} \,{=}\, [0,2\pi]\,{\times}\,[0,2\pi]$$ with spatial dimension $$n\,{=}\,2$$ and global Cartesian coordinates $$(x_1,x_2)$$. We combine the semidiscrete scheme introduced in the previous section with the implicit midpoint time discretization, with fixed time step $${\it {\Delta}} t$$. More specifically, let $$t^n := n {\it {\Delta}} t$$ for $$n \geq 0$$, and $$u^n \,{=}\, {\bf u}^n \cdot {\mathrm{d}} {\bf x} \,{:=}\, u|_{t=t^n}$$, then the fully discrete scheme is defined by Given $$u^n \,{\in}\, \overset{\circ}{W}_{h}$$, find $$u^{n+1}\,{\in}\, \overset{\circ}{W}_{h}$$ such that 1Δt(un+1−un,v)Ω+(Xun+un+12hun+un+12,v)Ω=0, (5.1) for all $$v \in \overset{\circ}{W}_{h}$$. It is well known that the implicit midpoint rule preserves the quadratic invariants of the semidiscrete system exactly. In our case, kinetic energy is conserved via Proposition 3.18. The tests of this section were performed using the Firedrake software suite (see Rathgeber et al., 2016), which allows for symbolic implementation of finite-element problems of mixed type. The nonlinear system was solved using Newton’s method, with LU factorization for the inner linear system implemented using the PETSc library (see Hendrickson & Leland, 1993; Balay et al., 1997, 2015; Dalcin et al., 2011). We start by examining the convergence rates with respect to a given manufactured solution. The manufactured solution solves a modified system of equations, where an appropriate forcing is introduced so that we can study the approximation properties of our discretization (although we cannot make conclusive remarks on stability). We pick as manufactured solution the Taylor–Green vortex (see Guzmán et al., 2016), u(t,x1,x2)=sin⁡(x1)cos⁡(x2)e−2t/σe1−cos⁡(x1)sin⁡(x2)e−2t/σe2, (5.2) with $${\it {\Omega}} = [0,2\pi]\times [0,2\pi]$$, $$\sigma=100$$ and $$t\in [0,1]$$. Table 1 lists the $$L^2$$-error and the convergence rates for both the centred and upwind scheme, for the case $${\bf W}_ h = {\bf{RT}}_s(\mathscr T_h)$$ and time step $${\it {\Delta}} t = 10^{-2}$$. The centred scheme was already analysed in Guzmán et al. (2016), where it was pointed out that without upwinding the rate of convergence $$s$$ is sharp for $$s=1$$, whereas convergence rates higher than $$s$$ can be observed for $$s=0$$ and $$s=2$$. The form of upwinding proposed in this article is different from the one of Guzmán et al. (2016). However, it gives the same behaviour in terms of order of convergence. In particular, we notice that introducing upwinding yields the optimal convergence rate $$s+1$$ for $$s\geq 1$$. Table 1 Comparison between centred ($$c_f=0$$) and upwind $$\left(c_f = \frac{{\mathbf u}\cdot{\mathbf n}_f}{2|{\mathbf u}\cdot{\mathbf n}_f|}\right)$$ scheme in terms of the error $$\|{\mathbf{u}}-\mathbf{u}_h\|_{{\it {\Omega}}}$$ and the order of convergence, for the Taylor–Green vortex test case and $${\mathbf W}_h={\mathbf{RT}}_s(\mathscr T_h)$$ Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 Table 1 Comparison between centred ($$c_f=0$$) and upwind $$\left(c_f = \frac{{\mathbf u}\cdot{\mathbf n}_f}{2|{\mathbf u}\cdot{\mathbf n}_f|}\right)$$ scheme in terms of the error $$\|{\mathbf{u}}-\mathbf{u}_h\|_{{\it {\Omega}}}$$ and the order of convergence, for the Taylor–Green vortex test case and $${\mathbf W}_h={\mathbf{RT}}_s(\mathscr T_h)$$ Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 Centred Upwind $$s$$ $$h$$ Error Order Error Order 0 7.40e-1 2.84e-1 — 4.01e-1 — 3.70e-1 1.42e-1 1.00 2.24e-1 0.84 2.47e-1 9.50e-2 1.00 1.58e-1 0.87 1.85e-1 7.12e-2 1.00 1.22e-1 0.89 1 7.40e-1 1.42e-1 — 2.15e-2 — 3.70e-1 7.13e-2 0.99 5.38e-3 1.99 2.47-1 4.76e-2 1.00 2.39e-3 2.00 1.85e-1 3.57e-2 1.00 1.35e-3 2.00 2 7.40e-1 1.81e-3 — 7.61e-4 — 3.70e-1 2.09e-4 3.11 9.02e-5 3.08 2.47e-1 6.28e-5 2.97 2.59e-5 3.08 1.85e-1 2.69e-5 2.96 1.07e-5 3.07 In the following, we will refer to the scheme of Guzmán et al. (2016) as the GSS scheme (by the initial of the authors), whereas we refer to the scheme proposed in this article as the Lie derivative scheme. We now consider a double shear problem (see, e.g., Liu & Shu, 2000; Guzmán et al., 2016) obtained by setting the initial velocity $${\bf u} = u_1{\bf e}_1 + u_2 {\bf e}_2$$ as follows u1|t=0:={tanh⁡((x2−π/2)/ρ)x2≤πtanh⁡((3π/2−x2)/ρ)x2>π ,u2|t=0:=δsin⁡(x1), (5.3) where $$\rho = \pi/15$$, $$\delta = 0.05$$ and with periodic boundary conditions on $${\it {\Omega}}$$. We use $${\it {\Delta}} t = 8/200$$ as in the references above to integrate the solution in time. In Figs 2–4, we plot the vorticity $$\mathrm{rot}\, {\bf u}$$ at $$t=8$$, for the centred scheme and the upwind GSS and Lie derivative schemes. The centred scheme produces oscillations in the vorticity field, as it was already observed in Guzmán et al. (2016). The presence of oscillations is reflected in the growth in enstrophy $$Z(t)\,{:=}\,\int_{{\it {\Omega}}} (\mathrm{rot}\, {\bf u})^2\, {\mathrm{d}} {\bf x}$$, which accumulates at small scales, see Fig. 5. As a matter of fact, unfortunately, our variational derivation does not imply enstrophy conservation, even though this is an invariant of the continuous system. On the other hand, the energy $$K(t) \,{:=}\, \int_{{\it {\Omega}}} \|{\bf u}\|^2\, {\mathrm{d}} {\bf x}$$ is exactly conserved by the scheme. Furthermore, we observe from Fig. 5 that with $$h$$ refinement, the enstrophy at the final time does not seem to converge to its initial value, when $$s=1$$, whereas we do observe convergence for $$s\,{=}\,2$$. This is expected since the suboptimal error estimates provided in Guzmán et al. (2016) are sharp for $$s\,{=}\,1$$, so we can only guarantee convergence of the vorticity in $$L^2$$ for $$s\geq 2$$ via an inverse estimate. Fig. 2. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, centred scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. White colour is used for values outside of range. Fig. 2. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, centred scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. White colour is used for values outside of range. Fig. 3. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind Lie derivative scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 3. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind Lie derivative scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 4. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind GSS scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 4. View largeDownload slide Vorticity function, $$\mathrm{rot}\, {\bf u}$$, at $$t=8$$, upwind GSS scheme, $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 5. View largeDownload slide Enstrophy history $$Z(t)$$ for the centred scheme with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 5. View largeDownload slide Enstrophy history $$Z(t)$$ for the centred scheme with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. The vorticity fields for the upwind GSS and Lie derivative schemes are similar and provide much better agreement to the reference solution in Liu & Shu (2000). The enstrophy time evolution in Fig. 6 indicates that the introduction of upwinding induces dissipation in enstrophy. However, the upwind Lie derivative schemes still preserves energy exactly, whereas the upwind GSS scheme also dissipates energy. Fig. 6. View largeDownload slide Kinetic energy, $$K(t)$$, and enstrophy history, $$Z(t)$$, for the upwind Lie derivative scheme (a) and the upwind GSS scheme (b), with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. Fig. 6. View largeDownload slide Kinetic energy, $$K(t)$$, and enstrophy history, $$Z(t)$$, for the upwind Lie derivative scheme (a) and the upwind GSS scheme (b), with $${\bf W}_h = {\bf{BDM}}_s(\mathscr T_h)$$. 6. Summary and outlook In this article, we extended the variational discretization proposed in Pavlov et al. (2009) for the incompressible Euler equations to the finite-element setting. We largely based our efforts on the work of Heumann & Hiptmair (2011, 2013) and Heumann et al. (2016), where the authors introduced Galerkin discretizations for the Lie derivative operator. Specifically, we used such operators as discrete Lie algebra variables, so that the discrete equations of motion could be derived directly from an appropriately defined Lagrangian and the Hamilton–d’Alembert’s principle. We found that the discretization obtained using this strategy coincides with the centred flux scheme proposed in Guzmán et al. (2016). Moreover, it has a built-in energy conservation property and satisfies a discrete Kelvin’s circulation theorem. We also introduced an upwind-stabilized version of the algorithm, which preserves both properties. Finally, we provided a convergence analysis proving (suboptimal) convergence rates for the upwind scheme for sufficiently high-order finite elements. Numerical tests suggest that this result might not be sharp, as we obtain optimal convergence rates even using low-order spaces. The methodology that we used in this article to discretize the incompressible Euler equations is general and can be applied to different fluid models that share a similar variational structure. The simplest example of such possible extensions is given by the Euler-alpha model (see, e.g., Holm et al., 1998). In this case, we keep the setting of the Euler equations described in Section 2.2, whereas the Lagrangian is given by l(u)=12∫Ω(‖u‖2+α2‖gradu‖2)vol, and $$\delta{l}/\delta{\bf u}$$ can be identified with an opportunely defined momentum $$m = A {\bf u} \cdot {\mathrm{d}}{\bf x}$$, with $$A = (I-\alpha^2{\it {\Delta}})$$ and $$I$$ being the identity map. Suppose that we have a discretization for $$A$$ given by $$A_h$$. Then, in the notation of Section 3, we can take $$A_h:{\bf W}_h \rightarrow {\bf W}_h$$ to be a linear invertible operator acting on the velocity finite-element space $${\bf W}_h$$. If we assume $$A_h$$ to be symmetric, the variational derivation we discussed in Section 3 can be applied without changes so that eventually we get an advection equation for the momentum $$m$$ expressed in terms of our discrete Lie derivative operator. Many MHD and GFD models can also be derived from a variational principle similar to the one that applies for perfect incompressible fluids (see, e.g., Holm et al., 1998). In these cases, one needs to add advected quantities to the system, which again require an appropriate definition for a discrete Lie derivative. Moreover, the variational derivation of Section 3 needs to be appropriately modified, although we expect to maintain its general features. This provides a direction for our future research. Footnotes 1The discrete Lie derivatives proposed in Heumann & Hiptmair (2011) are more general than considered in this article, since they account for a larger class of finite-element spaces, either conforming or not. In the language of differential forms, we use spaces that are conforming with respect to the $$L^2$$-adjoint of the exterior derivative; this case is considered explicitly in Section 4.1.1 of Heumann (2011). References Abraham R. , Marsden J. E. & Ratiu T. ( 1988 ) Manifolds, Tensor Analysis, and Applications , vol. 75 , 2nd edn . New York : Springer . Google Scholar CrossRef Search ADS Arnold D. N. , Falk R. S. & Winther R. ( 2006 ) Finite element exterior calculus, homological techniques, and applications. Acta Numer. , 15 1 – 155 . Google Scholar CrossRef Search ADS Arnol’d V. ( 1966 ) Sur la g éom étrie diff érentielle des groupes de Lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits. Ann. Inst. Fourier , 16 319 – 361 . Google Scholar CrossRef Search ADS Balay S. , Gropp W. D. , McInnes L. C. & Smith B. F. ( 1997 ) Efficient management of parallelism in object oriented numerical software libraries. Modern Software Tools in Scientific Computing ( Arge E. Bruaset A. M. & Langtangen H. P. eds). Boston, MA : Birkhäuser Press , pp. 163 – 202 . Google Scholar CrossRef Search ADS Balay S. , Abhyankar S. , Adams M. F. , Brown J. , Brune P. , Buschelman K. , Dalcin L. , Eijkhout V. , Gropp W. D. , Kaushik D. , Knepley M. G. , McInnes L. C. , Rupp K. , Smith B. F. , Zampini S. & Zhang H. ( 2015 ) PETSc users manual. Technical Report ANL-95/11—Revision 3.6. Lemont, IL : Argonne National Laboratory . Boffi D. , Fortin M. & Brezzi F. ( 2013 ) Mixed Finite Element Methods and Applications . Springer Series in Computational Mathematics. Berlin, Heidelberg : Springer . Google Scholar CrossRef Search ADS Brenier Y. ( 1989 ) The least action principle and the related concept of generalized flows for incompressible perfect fluids. J. Amer. Math. Soc. , 2 225 – 255 . Google Scholar CrossRef Search ADS Brenier Y. ( 2000 ) Derivation of the Euler equations from a caricature of coulomb interaction. Commun. Math. Phys. , 212 93 – 104 . Google Scholar CrossRef Search ADS Brezzi F. , Marini L. D. & Süli E. ( 2004 ) Discontinuous Galerkin methods for first-order hyperbolic problems. Math. Models Methods Appl. Sci. , 14 1893 – 1903 . Google Scholar CrossRef Search ADS Butler G. & Rogers T. ( 1971 ) A generalization of a lemma of Bihari and applications to pointwise estimates for integral equations. J. Math. Anal. Appl. , 33 77 – 81 . Google Scholar CrossRef Search ADS Cockburn B. , Kanschat G. & Schötzau D. ( 2005 ) A locally conservative LDG method for the incompressible Navier–Stokes equations. Math. Comp. , 74 1067 – 1095 . Google Scholar CrossRef Search ADS Dalcin L. D. , Paz R. R. , Kler P. A. & Cosimo A. ( 2011 ) Parallel distributed computing using Python. Adv. Water Resour. , 34 1124 – 1139 . Google Scholar CrossRef Search ADS Desbrun M. , Hirani A. N. , Leok M. & Marsden J. E. ( 2005 ) Discrete exterior calculus. arXiv preprint math/0508341 . Desbrun M. , Gawlik E. S. , Gay-Balmaz F. & Zeitlin V. ( 2014 ) Variational discretization for rotating stratified fluids. Discrete Contin. Dyn. Syst. , 34 477 – 509 . Ebin D. G. & Marsden J. E. ( 1970 ) Groups of diffeomorphisms and the motion of an incompressible fluid. Ann. Math. , 92 102 – 163 . Google Scholar CrossRef Search ADS Gallouët T. & Mérigot Q. ( 2017 ) A Lagrangian scheme à la Brenier for the incompressible Euler equations . Accepted for publication in Foundations of Computational Mathematics . Guzmán J. , Shu C.-W. & Sequeira F. A. ( 2016 ) H (div) conforming and DG methods for incompressible Euler’s equations. IMA J. Numer. Anal. Available at https://doi.org/10.1093/imanum/drw054. Hendrickson B. & Leland R. ( 1993 ) A Multilevel Algorithm for Partitioning Graphs . Tech. Rep. SAND93-1301. Albuquerque, NM : Sandia National Laboratories . Heumann H. ( 2011 ) Eulerian and semi-Lagrangian methods for advection-diffusion of differential forms. Ph.D. Thesis , ETH Zürich , Switzerland . Heumann H. & Hiptmair R. ( 2011 ) Eulerian and semi-Lagrangian methods for advection-diffusion for differential forms. Discrete Contin. Dyn. Syst. , 29 1471 – 1495 . Google Scholar CrossRef Search ADS Heumann H. & Hiptmair R. ( 2013 ) Stabilized Galerkin methods for magnetic advection. ESAIM Math. Model. Numer. Anal. , 47 1713 – 1732 . Google Scholar CrossRef Search ADS Heumann H. , Hiptmair R. & Pagliantini C. ( 2016 ) Stabilized Galerkin for transient advection of differential forms. Discrete Contin. Dyn. Syst. Ser. S , 9 185 – 214 . Holm D. D. , Marsden J. E. & Ratiu T. S. ( 1998 ) The Euler-Poincaré equations and semidirect products with applications to continuum theories. Adv. Math. , 137 1 – 81 . Google Scholar CrossRef Search ADS Liu J.-G. & Shu C.-W. ( 2000 ) A high-order discontinuous Galerkin method for 2D incompressible flows. J. Comput. Phys. , 160 577 – 596 . Google Scholar CrossRef Search ADS Marsden J. E. & Ratiu T. S. ( 2010 ) Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems . New York : Springer Science & Business Media . McLachlan R. I. ( 1993 ) Explicit Lie–Poisson integration and the Euler equations. Phys. Rev. Lett. , 71 3043 – 3046 . Google Scholar CrossRef Search ADS PubMed Mérigot Q. & Mirebeau J.-M. ( 2016 ) Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport. SIAM J. Numer. Anal. , 54 3465 – 3492 . Google Scholar CrossRef Search ADS Mullen P. , McKenzie A. , Pavlov D. , Durant L. , Tong Y. , Kanso E. , Marsden J. E. & Desbrun M. ( 2011 ) Discrete Lie advection of differential forms. Found. Comput. Math. , 11 131 – 149 . Google Scholar CrossRef Search ADS Pavlov D. , Mullen P. , Tong Y. , Kanso E. , Marsden J. E. & Desbrun M. ( 2009 ) Structure-preserving discretization of incompressible fluids. Phys. D , 240 443 – 458 . Google Scholar CrossRef Search ADS Rathgeber F. , Ham D. A. , Mitchell L. , Lange M. , Luporini F. , Mcrae A. T. T. , Bercea G.-T. , Markall G. R. & Kelly P. H. J. ( 2016 ) Firedrake: automating the finite element method by composing abstractions. ACM Trans. Math. Software , 43 , 24:1 – 24:27 . Google Scholar CrossRef Search ADS Zeitlin V. ( 1991 ) Finite-mode analogs of 2D ideal hydrodynamics: Coadjoint orbits and local canonical structure. Phys. D: Non. Phen. , 49.3 , 353 – 362 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jun 29, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off