Semi-Lagrangian Finite-Element Exterior Calculus for Incompressible Flows

We develop a mesh-based semi-Lagrangian discretization of the time-dependent incompressible Navier-Stokes equations with free boundary conditions recast as a non-linear transport problem for a momentum 1-form. A linearly implicit fully discrete version of the scheme enjoys excellent stability properties in the vanishing viscosity limit and is applicable to inviscid incompressible Euler flows. Conservation of energy and helicity are enforced separately.


Incompressible Navier-Stokes Equations
We consider the incompressible Navier-Stokes equations-a standard hydrodynamic model for the motion of an incompressible, potentially-viscous fluid-in a container with rigid walls, where we impose so-called "free boundary conditions" in the parlance of (Mitrea and Monniaux, 2009, p. 346) and (Temam and Ziane, 1996, p. 502), see the latter article for further references.We search the fluid velocity field u(t, x) and the pressure p(t, x) as functions of time t and space x on a bounded, Lipschitz domain 1 arXiv:2301.04923v2[math.NA] 2 Feb 2024 Ω ⊂ R d such that they solve the evolution boundary-value problem ∇ • u = 0, on ]0, T [×Ω, (1b) ϵn × ∇ × u = 0, on ]0, T [×∂Ω, (1d) where ϵ ≥ 0 denotes a (non-dimensional) viscosity, f a given source term, T > 0 the final time, ∂Ω the boundary of Ω, and n(x) the outward normal vector at x ∈ ∂Ω.
The initial condition u 0 is to satisfy ∇•u 0 = 0 in Ω and u 0 •n = 0, ϵn×∇×u 0 = 0 on ∂Ω.Based on the variational description of the Navier-Stokes equations as described in (Arnold, 1966), u can be interpreted as a differential 1-form (Natale and Cotter, 2018) and we can recast system (1) in the following way.Let Λ k (Ω) for k ∈ N denote the space of differential k-forms on Ω.Then we search ω ∈ Λ 1 (Ω) and p ∈ Λ 0 (Ω) such that tr ⋆ω = 0, on ]0, T [×∂Ω, ϵ tr ⋆dω = 0, on ]0, T [×∂Ω, (2d) where D u ω denotes the material derivative of ω with respect to u, d : Λ k−1 (Ω) → Λ k (Ω) the exterior derivative, δ : Λ k (Ω) → Λ k−1 (Ω) the exterior coderivative, and the trace tr is the pullback under the embedding ∂Ω ⊂ Ω.Here, u is related to ω through ω := u Z , i.e. u is the vector proxy of ω w.r.t. the Euclidean metric.Similarly, we have that Λ 1 (Ω) ∋ ω 0 := u 0 Z and Λ 1 (Ω) ∋ f := f Z .We would like to emphasize that ω does not represent the vorticity, but the 1-form representation of the velocity in this work.Note that (2) can be derived through classical vector calculus for vector proxies as shown in section A for d = 3.As shown in (Chorin and Marsden, 1993;De Rosa, 2020), sufficiently-smooth solutions ω :]0, T [ → Λ 1 (Ω) of the incompressible Navier-Stokes equations as given in system (2) satisfy an energy relation, that is, This relation implies energy conservation for ϵ = 0 and f = 0. Note that the Onsager conjecture tells us that in the case ϵ = 0 the solutions need to be at least Hölder regular with exponent 1 3 for energy conservation to hold (Isett, 2018).Remark 1.We acknowledge that the boundary condition (1d) is non-standard.This boundary condition was chosen because it is the natural boundary condition associated to system (2).At first glance, an intuitive method to enforce the standard no-slip boundary conditions is to replace (1d) by ϵu × n = 0 on ]0, T [×∂Ω.To obtain wellposedness, it is then also required to impose p = 0 on ]0, T [×∂Ω.However, in that case, we lose the natural boundary condition u • n = 0. Instead, ϵu × n = 0 can be enforced using boundary penalty methods as will be shown in future work.In the case ϵ = 0, the only imposed boundary condition (1c) is standard.Remark 2. Boundary conditions (1c),(1d) can be interpreted as slip boundary conditions.However, on smooth domains Ω, they are only equivalent to Navier's slip boundary conditions if the Weingarten map related to ∂Ω vanishes (Mitrea and Monniaux, 2009, section 2).

Outline and Related Work
We propose a semi-Lagrangian approach to the discretization of the reformulated Navier-Stokes boundary value problem (2).This method revolves around the discretization of the material derivative D u ω in the framework of a finite-element Galerkin discretization on a fixed spatial mesh.The main idea is to approximate D u ω by backward difference quotients involving transported snapshots of the 1-form ω, which can be computed via the pullback induced by the flow of the velocity vector field u.
Semi-Lagrangian methods for transient transport equations like (2) are wellestablished for the linear case when u is a given Lipschitz-continuous velocity field.In particular, for ω a 0-form, that is, a plain scalar-valued function, plenty of semi-Lagrangian approaches have been proposed and investigated, see, for instance, (Bercovier and Pironneau, 1982;Bercovier et al, 1983;Douglas and Russell, 1982;Ewing et al, 1984;Hasbani et al, 1983;Pironneau, 1982;Russell, 1985;Süli, 1988;Bause and Knabner, 2002;Bermejo and Saavedra, 2012;Wang and Wang, 2010).We refer to (Heumann and Hiptmair, 2013, Chapter 5) for a comprehensive pre-2013 literature review on the analysis of general semi-Lagrangian schemes.Most of these methods focus on mapping point values under the flow, with the exception of a particularly interesting class of semi-Lagrangian methods known as Lagrange-Galerkin methods.Lagrange-Galerkin methods do not transport point values, but rather triangles (in 2D) or tetrahedra (in 3D).Refer to (Bermejo and Saavedra, 2016) for a review of those methods.
Meanwhile semi-Lagrangian methods for transport problems for differential forms of any order have been developed (Heumann et al, 2012;Heumann andHiptmair, 2011, 2013).The next section will review these semi-Lagrangian methods for linear transport problems with emphasis on 1-forms.We will also introduce a new scheme which is second-order in space and time based on so-called "small edges", see section 3.1.2for details.
It is important to note that these methods require the evaluation of integrals of transported quantities and, in case these integrals cannot be computed exactly, instabilities can occur (Bermejo and Saavedra, 2012;Morton et al, 1988).A possible remedy is to add an additional stabilization term that includes artificial diffusion (Bermejo and Saavedra, 2016).Other semi-Lagrangian methods for incompressible Navier-Stokes equations directly transport point values with the nodes of a mesh instead of evaluating integrals of transported quantities, see (Patera, 1984;Karniadakis and Sherwin, 2005;Xiu and Karniadakis, 2001;Xiu et al, 2005;Bonaventura et al, 2018) and (Celledoni et al, 2016), where the last work makes use of exponential integrators (Celledoni et al, 2003).Most authors employ spectral elements for the discretization in space (Patera, 1984;Karniadakis and Sherwin, 2005), but any type of finite-element space with degrees-of-freedom relying on point evaluations can be used.The methods proposed in (Xiu and Karniadakis, 2001;Xiu et al, 2005) are also based on finite-element spaces with degrees-of-freedom on nodes, but employ backward-difference approximations for the material derivative.The work (Bonaventura et al, 2018) proposes an explicit semi-Lagrangian method still built around the transport of point values in the nodes of the mesh.The diffusion term is also taken into account in a semi-Lagrangian fashion and the incompressibility constraint is enforced by means of a Chorin projection.Also (Bonaventura et al, 2018) proposes an explicit semi-Lagrangian scheme using the same principles, but based on the vorticity-streamfunction form of the incompressible Navier-Stokes equations.
All the mentioned semi-Lagrangian schemes rely on the transport of point values of continuous vector fields, which is the perspective embraced in formulation (1).However, we believe that, in particular in the case of free boundary conditions (1c) and (1d), the semi-Lagrangian method based on (2) offers benefits similar to the benefits bestowed by the use of discrete differential forms (finite-element exterior calculus, FEEC (Arnold, 2018;Arnold et al, 2006)) for the discretization of electromagnetic fields.Section 4 will convey that the boundary conditions (2c), (2d), and the incompressiblity constraint can very naturally be incorporated into a variational formulation of (2) posed in spaces of 1-forms.This has been the main motivation for pursuing the new idea of a semi-Lagrangian method for (2) that employs discrete 1-forms.Another motivation has been the expected excellent robustness of the semi-Lagrangian discretization in the limit ϵ → 0. Numerical tests reported in section 5 will confirm this.
Two more aspects of our method are worth noting: Firstly, a discrete 1-form ω h will not immediately spawn a continuous velocity field u h = ω Z h , However, continuity is essential for defining a meaningful flow.We need an additional averaging step, which we present in section 4.1.Secondly, since semi-Lagrangian methods fail to respect the decay/conservation law (3) exactly, we present a way how to enforce them in section 4.3.(a) 9 small edges of a second-order element in 2D.All the edges between the different connection points are small edges.In 3D, we simply have all these small edges on the faces of the simplex.edge no.l.s.f.edge no.l.s.f.
x 2 (b) Local shape functions (l.s.f.) for the unit triangle associated with second-order, discrete differential forms in 2D as proposed in (Rapetti and Bossavit, 2009).Each shape function corresponds to the small edge in (a) with the same numbering.3 Semi-Lagrangian Advection of differential forms

Discrete differential forms
We start from a simplicial triangulation K h (Ω) of Ω, which may rely on a piecewise linear approximation of ∂Ω so that it covers a slightly perturbed domain.

Lowest-order case: Whitney forms
For Λ 0 (Ω)-the space of 0-forms on Ω, which is just a space of real-valued functions-the usual (Lagrange) finite-element space of continuous, piecewise-linear, polynomial functions provides the space Λ 0 h,1 (Ω) of lowest-order discrete 0-forms.Let d ∈ {2, 3}, K a d-simplex with edges {e 1 , .., e 3(d−1) }.To construct lowest-order discrete 1-forms on K, we associate to every edge e i a local shape function w ei .Let the edge e i be directed from vertex v 1 i to v 2 i , then the local shape function w ei ∈ Λ 1 (K) associated with edge e i is where λ v represents the barycentric coordinate associated with vertex v.We define the lowest-order, local space of discrete 1-forms Λ 1 h,1 (K) := span{w e ; e an edge of K}. (5) Using these local spaces, we can define the global space of lowest-order, discrete 1-forms where Λ 1 (Ω) again denotes the space of differential 1-forms on Ω.We demand that for every ω ∈ Λ 1 (Ω) integration along any smooth oriented path yields a unique value.Thus, the requirement ω ∈ Λ 1 (Ω) imposes tangential continuity on the vector proxy of ω.

Second-order discrete forms
Similar to the lowest-order case, the space Λ 0 h,2 (Ω) of second-order discrete 0-forms is spawned by the usual (Lagrange) finite-element space of continuous, piecewisequadratic, polynomial functions.
Let d ∈ {2, 3}, K a d-simplex with edges {e 1 , .., e 3(d−1) } and vertices {v 1 , .., v d+1 }.To construct second-order discrete 1-forms, we associate local shape functions to "small edges".We can construct 3(d + 1)(d − 1) small edges (Rapetti and Bossavit, 2009, Definition 3.2) by defining ∀i ∈ {1, .., d + 1} and ∀j ∈ {1, .., where {v i , e j } denotes the small edge.In Figure 1a we illustrate the 9 small edges of a 2-simplex.For example, we see that small edge 9 can be written as {(0, 0), [(1, 0), (0, 1)]}.To make the difference between small edges and edges of the mesh explicit, we will sometimes refer to the latter as "big edges".The local shape function (Rapetti and Bossavit, 2009, Definition 3.3) associated with {v i , e j } is given by w {vi,ej } := λ vi w ej , where w ej denotes the Whitney 1-form associated with the big edge e j as defined in (4).In Figure 1b we give explicit expressions for the shape functions associated with the small edges in Figure 1a.Note that the local shape functions of the form w {v,e} associated with small edges in the interior (d = 2) or on the same face (d = 3) of the form {v, e} such that v / ∈ ∂e (example: small edge 7, 8, and 9 in Figure 1a) are linearly dependent.We define the second-order, local space of discrete 1-forms (Rapetti and Bossavit, 2009, Definition 3.3)   Λ 1 h,2 (K) := span{w {v,e} ; v a vertex of K, e a (big) edge of K}.
Using these local spaces, we can define the global space of second-order, discrete 1-forms where again we have tangential continuity by a similar argument as in section 3.1.1.

Projection operators
We denote by E h,p (Ω) the global set of big edges (p = 1) or small edges (p = 2) associated with K h (Ω).We will define the projection operator is minimized.Note that for p = 1, this mismatch can be made to vanish.In this case, I h,1 agrees with the usual edge-based nodal projection operator (Hiptmair, 2002, Eq. (3.11)).
In practice, we can compute the projection locally as follows.Let K ∈ K h (Ω) be a dsimplex, d ∈ {2, 3}, and let {s 1 , .., s N p,d } and {w s1 , .., w s N p,d } denote the corresponding big (p = 1) or small (p = 2) edges and corresponding shape functions as introduced above.Specifically, we have N 1,2 = 3, N 1,3 = 6, N 2,2 = 9, and N 2,3 = 24.We can define the matrix We will say that there is an interaction from edge s j to s i if (M) i,j ̸ = 0. Note that for p = 1, M is the identity matrix.For p = 2 the local shape functions are linearly dependent and, thus, the above matrix is not invertible.However, we can decompose M into invertible and singular sub-matrices.For illustrative purposes we display for p = 2 and d = 2 the decomposition of M in Figure 2b.The three top-left sub-matrices in Figure 2b are invertible 2×2 matrices that describe the interaction between the two small edges that lie on the same big edge, that is, the blue, red, and green submatrix in Figure 2b correspond to the blue, red, and green small edges in Figure 2a, respectively.The orange submatrix in Figure 2b is a 3 × 3 matrix with rank 2 that describes the interaction between the three small edges that lie in the interior of the simplex Fig. 2: For p = 2 and d = 2 the matrix M corresponding to the 2-simplex K in (a) has the form given in (b).Each row and column in M is associated to a small edge in (a).Each submatrix in (b) describes the interactions between edges with the same color in (a).The gray submatrix is an exception as it describes the one-directional interaction between the small edges that lie on a big edge and the small edges that lie in the interior.For d = 3, M has the structure as shown in Figure 2c, where the purple submatrices are 2 × 2 invertible matrices and the orange submatrices are 3 × 3 matrices of rank 2.
in Figure 2a, that is, the orange small edges.The gray submatrix encodes the onedirectional interaction from the the small edges that lie on a big edge to the small edges in the interior.Note that the decomposition of M as given in Figure 2b is not limited to d = 2.The idea can be extended to d = 3 by considering each face of a 3simplex as a 2-simplex.This is sufficient, since for d = 3 we have no small edges in the interior and there is no interaction between small edges that do not lie on the same face.We give the general structure of M in Figure 2c.Note that the small, purple sub-matrices represent invertible 2 × 2 matrices and the bigger, orange sub-matrices represent 3 × 3 matrices with rank 2. In order to find We can then compute ⃗ η K as a least-squares solution of Without loss of generality we assume that M has the form as given in Figure 2c.Then, we solve (12) as follows: 1.The local shape functions related to small edges that lie on a big edge of the simplex are linearly independent.We solve for their coefficients first, that is, we solve the system corresponding to the invertible blue sub-matrices in Figure 2c first.2. Using the results from step 1, we can move the gray submatrix in Figure 2c to the right-hand side.Then, we solve the matrix-system corresponding to the orange sub-matrices in Figure 2c in a least-squares sense.
If we perform the above steps for all K ∈ K h (Ω), we find ω h = I h,p ω ∈ Λ 1 h,p (Ω).Note that only the shape functions associated to small edges on a face contribute to the tangential fields on that face.Therefore, the above procedure yields tangential continuity.Remark 3.For p = 1, (12) reduces to with s i a big edge of the d-simplex K for all i ∈ {1, .., 3(d−1)}.This yields the standard nodal interpolation operator of (Hiptmair, 2002, Eq. (3.11)).

Semi-Lagrangian material derivative
The method described in this section is largely based on (Heumann and Hiptmair, 2013;Heumann et al, 2012).Throughout this section, unless stated otherwise, we fix the stationary, Lipschitz-continuous velocity field This means that we consider a linear transport problem and our main concern will be the discretization of the material derivative D u ω for a 1-form ω.We can define the flow ]0, Given that flow we can define the material derivative for a time-dependent differential 1-form ω We employ a first-or second-order, backward-difference method to approximate the derivative.Writing X * t,t−τ for the pullback of forms under the flow, we obtain for sufficiently-smooth t → ω(t) and a timestep 0 < τ → 0 or respectively.Note that both backward-difference methods are A-stable (Süli and Mayers, 2003).In the remainder of this section we restrict ourselves to ( 16), but exactly the same considerations apply to (17).Given a temporal mesh .. < t n < t n+1 < .., we approximate ω(t n , •) ∈ Λ 1 (Ω) by a discrete differential form ω n h ∈ Λ 1 h,p (Ω) with p ∈ {1, 2}.Using the backward-difference quotient (16), we can define the discrete material derivative for fixed timestep τ > 0 where we need to use the projection operator Ω) in general.Recall from section 3.1 that the degrees of freedom for discrete 1-forms are associated to small (p = 2) or big (p = 1) edges.As discussed in section 3.1.3,evaluating the interpolation operator entails integrating X * t,t−τ ω n−1 h over small (p = 2) or big (p = 1) edges.We can approximate these integrals as follows where e is a small or big edge and with v 1 , v 2 the vertices of e. Instead of transporting the edge e using the exact flow X t,t−τ , we follow (Bermejo and Saavedra, 2012;Heumann and Hiptmair, 2013;Heumann et al, 2012) and only transport the vertices of the small edges (p = 2) or big edges (p = 1) and obtain a piecewise linear (second-order) approximation Xt,t−τ (e) of the transported edge X t,t−τ (e) as illustrated in Figure 3.We can approximate the movement of the endpoints of e under the flow as defined by ( 14) by solving ( 14) using the explicit Euler method or Heun's method for the first-and second-order case, respectively.We will elaborate on this further in section 4.1.
In Figure 3, we can also see that the approximate transported edge may intersect several different elements of the mesh.When we evaluate the integral in (19), it can happen that there are discontinuities of ω n−1 h along Xt,t−τ (e).Therefore, we cannot apply a global quadrature rule to the entire integral.Instead, we split Xt,t−τ (e) into several segments defined by the intersection of Xt,t−τ (e) with cells of the mesh.In our implementation, for the sake of stability, we find the intersection points by transforming back to a reference element as visualised in Figure 4. Algorithm 1 gives all details.Note that we can forgo the treament of any special cases (e.g.intersection with vertices) without jeopardizing stability.After we split the transported edge into segments, we can evaluate the integrals over these individual pieces exactly, because we know that ω n−1 h is of polynomial form when restricted to individual elements of the mesh (see section 3.1).
When simulating the fluid model (2), we will not have access to an exact velocity field.Instead we only have access to an approximation of the velocity field.This Fig. 3: Edge e (in red) is transported using the flow β (in blue).The exact transported edge X τ (e) and the approximate transported edge Xτ (e) are given in orange and green.
Algorithm 1 Splitting 1-simplex over mesh elements (see Figure 4 for illustration).Here, K ref denotes the reference simplex.
Require: x 0 ∈ K 0 ∈ K h (Ω) and x 1 vertices of a 1-simplex e. Ensure: Find the isoparametric mapping ϕ K : E ← E ∪ {K} 13: end while approximation may not satisfy exact vanishing normal boundary conditions.Therefore, a part of Xt,t−τ (e) may end up outside the domain.This can also happen due to an approximation of the flow by explicit timestepping.Since ω n−1 h is not defined Fig. 4: The red line indicates the line that spans multiple elements.On the left we see the reference triangle associated with the yellow element in the mesh on the right.outside the domain, we set where Xt,t−τ (e) R d \Ω is the part of Xt,t−τ (e) that lies outside the domain Ω and len •) gives the arclength of the argument.This is motivated by the situation displayed in Figure 5-a case where an edge gets transported out of the domain due to the use of approximate flow maps despite vanishing normal components of the velocity.If we set the value defined in (21) to zero in this case, it would be equivalent to applying vanishing tangential boundary conditions, which is inconsistent with (1).Instead, (21) just takes the tangential components from the previous timestep.
We arrive at the following approximation of the material derivative where the only difference between ( 18) and ( 22) is that X t,t−τ was replaced by Xt,t−τ and I h,p is implemented based on (21).Note that in our scheme X * t,t−τ is always evaluated in conjunction with I h,p , which means that we need define Xt,t−τ only on small (p = 2) or big (p = 1) edges.In fact, Xt,t−τ is defined through (20) for all points that lie on small (p = 2) or big (p = 1) edges.
Given a velocity field u ∈ W 1,∞ (Ω) with u • n = 0, it was shown in (Heumann et al, 2012, section 4) that using a first-order backward difference scheme and lowestorder elements for the spatial discretization, we can approximate a smooth solution Fig. 5: A coarse triangulation of Ω = {x ∈ R 2 ; ||x|| < 1} with the velocity field u = [−y, x] T satisfying u • n = 0. Despite the vanishing normal components of the velocity, the blue edge gets transported out of the domain to the green edge.
where h is the spatial meshwidth and τ > 0 is the timestep size.However, numerical experiments (Heumann et al, 2012, section 6) performed with τ = O(h) show an error of O(h)-a slight improvement over the a-priori estimates.This motivates us to link the timestep to the mesh width as τ = O(h).

Semi-Lagrangian Advection applied to the Incompressible Navier-Stokes Equations
Given a temporal mesh t 0 < t 1 < ... < t N −1 < t N , we elaborate a single timestep

Approximation of the flow map
In the Navier-Stokes equations, the flow is induced by the unknown, time-dependent velocity field u(t, x).Therefore, (14) becomes The discretization of the material derivative requires us to approximate the flow map X t,t−τ in order to evaluate (20).

A first-order scheme
We use the explicit Euler method to approximate the (backward) flow according to where t n is a node in the temporal mesh and τ denotes the timestep size.We only have access to an approximation Note that a direct application of the explicit Euler method would require an evaluation of the velocity field at t n .Instead, we perform a constant extrapolation and evaluate the velocity field at t n−1 , that is, we use u n−1 h in ( 26).The approximation u n−1 h resides in the space of vector proxies for discrete differential 1-forms as discussed in section 3.1.This means that only tangential continuity of u n−1 h across faces of elements of the mesh is guaranteed, while discontinuities may appear in the normal direction of the faces.Therefore, u n−1 h is not defined pointwise-even though (26) requires point-wise evaluation.For that reason, we will replace u n−1 h by a globally-continuous, smoothened velocity field ūn−1 h that approximates u n−1 h (see section 4.1.3for the construction).We then have which yields a first-order-in-time approximation of X tn,tn−τ (x), provided that ūn−1 h is a first-order approximation of u(t n−1 , •).

A second-order scheme
A second-order approximation can be achieved by using Heun's method (Vuik et al, 2015) instead of explicit Euler.We find the following second-order in time approximations where we approximate the velocity field at t n by the linear extrapolation u * h = 2u n−1 h − u n−2 h .As described in section 4.1.1,we replace the velocity fields by suitable smooth approximations.We obtain where ū• h with • = * , n − 1, n − 2 denotes the smoothened version of u • h as it will be constructed in the next section.

Smooth reconstruction of the velocity field
Given a discrete velocity field u Z h ∈ Λ 1 h,p (Ω), we can define a smoothened version ūh of u h that is • Lipschitz continuous to ensure stable evaluation of ( 27), • well-defined on every point of the meshed domain, • practically computable, and • second-order accurate.
We introduce ūh as follows.Let h min denote the length of the shortest edge of the mesh and (u i h ) i=1,..,d the components of u h .Then, provides a second-order, Lipschitz-continuous approximation of u h .In the above definition, we can also replace h min by a localized parameter that scales as O(h) with h the length of the edges "close" to x.Note that the above integral can be evaluated up to machine precision using the algorithm as described in section 3.2 for ( 19).The averaging (32) provides a second-order approximation of u h on every point in the mesh.

A first-and second-order SL scheme
We are now ready to turn the ideas of section 3 into a concrete numerical scheme for the incompressible Navier-Stokes equations as given in (2).We cast (2a) and ( 2b) into weak form and, subsequently, do Galerkin discretization in space relying on those spaces of discrete differential forms introduced in section 3.1.For the first-order scheme, we have the following discrete variational formulation.Given for all η h ∈ Λ 1 h,1 (Ω) and ψ h ∈ Λ 0 h,1 (Ω).I h,p denotes the projection operator as defined in section 3.1.For the second-order scheme, we use second-order timestepping and second-order discrete differential forms.Given for all η h ∈ Λ 1 h,2 (Ω) and ψ h ∈ Λ 0 h,2 (Ω).Numerical experiments reported in section 5 give evidence that these schemes indeed do provide first-and second-order convergence for smooth solutions.Note that the schemes presented in this section only require solving symmetric, linear systems of equations at every time-step.

Conservative SL schemes
In order to enforce energy-tracking-the correct behavior of the total energy E(t) over time as expressed in (3)-we add a suitable constraint plus a Lagrange multiplier to the discrete variational problems proposed in section 4.2.Given for all η h ∈ Λ 1 h,1 (Ω) and ψ h ∈ Λ 0 h,1 (Ω).Note that the last scalar equation enforces energy conservation for ϵ = 0 and f = 0. To solve the nonlinear system (35) for ω n h , p n h , µ n , we propose the following iterative scheme.Assume that we have a sequence Then we can employ the Newton-like linearization We use the above expansion to replace the quadratic terms (ω n , ω n ) Ω and (dω n , dω n ) Ω and arrive at the following linear variational problem to be solved in every step of the inner iteration.Given . This is a symmetric, linear system that is equivalent to the original system in the limit (ω In numerical experiments we observe that it takes around 2-3 steps of the inner iteration to converge to machine precision using an initial guess ω n h,0 = ω n−1 h .We can apply the same idea for energy-tracking to our second-order scheme as proposed in section 4.2.

Numerical Results
In this section, we present multiple numerical experiments to validate the new scheme.In the following, we will always consider schemes that include energy-tracking as introduced in section 4.3 unless explicitly stated otherwise.The experiments are based on a C++ code that heavily relies on MFEM (Anderson et al, 2021).The source code is published under the GNU General Public License in the online code repository https://gitlab.com/WouterTonnon/semi-lagrangian-tools.Unless specified otherwise, we use uniform meshes for the experiments.

Experiment 1: Decaying Taylor-Green Vortex
We consider the incompressible Navier-Stokes equations with Ω = [− 1 2 , 1 2 ] 2 , T = 1, varying ϵ ≥ 0, f = 0, and vanishing boundary conditions.An exact, classical solution is the following Taylor-Green vortex (Taylor and Green, 1937) We ran a h-convergence analysis for different values of ϵ ≥ 0 and summarize the results in Figure 6.We also track the energy for different values of ϵ and compare the energy to the exact solution in Figure 7.

Experiment 2: Taylor-Green Vortex
We consider the incompressible Navier-Stokes equations with Ω = [−1, 1] 2 , T = 1, varying ϵ ≥ 0, f and the boundary conditions chosen such that u(t, x) = cos(πx 1 ) sin(πx 2 ) − sin(πx 1 ) cos(πx 2 ) (39) is an exact, classical solution.We ran a h-convergence analysis for all parameters and summarize the results in Figure 8.We observe first-and second-order algebraic convergence for the corresponding schemes.Note that the error of the scheme is stable as ϵ → 0. This is in agreement with the analysis performed on the vectorial advection equations presented in (Heumann and Hiptmair, 2013).This experiment thus suggests that this analysis can be extended to the scheme presented in this work.
The exact solution to this problem is unknown, so we compare the solution computed by our scheme to the solution produced by the incompressible Euler solver Gerris (Popinet, 2007).The algorithm used in this solver is described in (Popinet, 2003).We computed the solution to this problem using the second-order, energy-tracking scheme presented in this work.Then, we plotted the magnitude of the computed velocity vector field for different mesh-sizes and time-steps at different time instances in Figures 10 to 13.Note that different visualisation tools were used to visualize the fields computed using the different solvers, but we observe that the solution computed by the semi-Lagrangian scheme comes visually closer to the solution computed by Gerris as we decrease the mesh width and time step.This is confirmed by Figure 9, where we display the L2 error between the solution computed using the semi-Lagrangian scheme L2 Error u second-order, cons Fig. 9: Convergence results for Experiment 3 using the second-order, conservative scheme on simplicial meshes with mesh width h, timestep τ = 0.06580h, and final time T = 1.The reference solution is a solution computed by Gerris (Popinet, 2003) to achieve an order of convergence that is between first-and second-order, but this may be pre-asymptotic behaviour.

Experiment 5: Lid-driven cavity with slippery walls
In this section, we simulate a situation that resembles a lid-driven cavity problem.Consider the incompressible Navier-Stokes with Ω = [− 1 2 , 1 2 ] 2 , T = 7.93, ϵ = 0, vanishing normal boundary conditions and the initial velocity field is set equal to zero.Then, to simulate a moving lid at the top, we apply the force-field f This force field gives a strong force in the x 1 -direction close to the top lid, but quickly tapers off to zero as we go further from the top lid.In Figure 18, we display the computed velocity field.Note that, because we apply slip boundary conditions, we do not expect to observe vortices.The numerical solution reproduces this expectation.

Experiment 6: A more complicated domain
The numerical experiments given above, show the convergence and conservative properties of the introduced schemes.However, these experiments are all performed on very simple, rectangular domains.In this experiment, we consider a more complicated domain and mesh (generated using (Geuzaine and Remacle, 2009)) as shown in Figure 19.
We consider the case of the incompressible Navier-Stokes equations on the domain as given in Figure 19, T = 100, ϵ = 0, f = 0 and vanishing normal boundary conditions.We need to construct an initial condition that is divergence-free with vanishing normal boundary conditions.Following an approach close to a Chorin projection, we start with We use this definition to define a scalar function, ϕ, as  16: The L2 norm of the computed solutions for Experiment 3 using different variants of the semi-Lagrangian scheme on a simplicial mesh with mesh width h = 0.04748975 and time-step τ = 0.003125.In the legend, 'cons' is short for 'conservative' and refers to energy-tracking schemes.We use ϵ = 0 and f = 0.In this figure, the time-step is too small to indicate samples by stars as done in similar figures.
We can define our initial condition, u 0 , as Note that u 0 is divergence-free and has vanishing normal boundary conditions.
The above system of equations can be solved using an appropriate finite-element implementation.
Note that in this experiment, the field outside the domain is unknown.This is well-defined on a continuous level, since vanishing boundary conditions imply that no particle will flow in from outside the domain.However, on the discrete level we cannot guarantee that the same will happen.It could happen that a part of a transported edge (as discussed in section 3.2) ends up outside the domain.In this case, we will assume that the average of the vector field along the part of the edge that lies outside the domain, will have the same value as the average of the corresponding edge in its original location (before transport) at the previous timestep.
The first ten seconds were simulated and a video of the results can be found at https://youtu.be/Eica8XHLtxY.For the different schemes, we also tracked the energy in Figure 20.

Conclusion
We have developed a mesh-based semi-Lagrangian discretization of the time-dependent incompressible Navier-Stokes equations with free boundary conditions recast as a nonlinear transport problem for a momentum 1-form.A linearly implicit fully discrete version of the scheme enjoys excellent stability properties in the vanishing-viscosity limit and is applicable to inviscid incompressible Euler flows.However, in this case conservation of energy has to be enforced separately.Making the reasonable choice of a time-step size proportional to the mesh width, the algorithm involves only local computations.Yet, these are significantly more expensive compared to those required for purely Eulerian finite-element and finite-volume methods.At this point the verdict on the competitiveness of our semi-Lagrangian scheme is still open.0 0.5 1 1.5 2 2.5 Fig. 18: Velocity field at T = 7.93s of Experiment 5 computed using the secondorder, non-conservative semi-Lagrangian scheme on a simplicial mesh with mesh width h = 0.189959 and τ = 0.01.Using the gradient of the dot-product, we find This identity allows us to rewrite the momentum equation to From (Hiptmair, 2002;Heumann and Hiptmair, 2013), we obtain the identity Fig. 20: The L2 norm of the computed solutions for Experiment 6 using different variants of the semi-Lagrangian scheme on a simplicial mesh as given in Figure 19 and time-step τ = 0.01.In the legend, 'cons' is short for 'conservative' and refers to energy-tracking schemes.In this figure, the time-step is too small to indicate samples by stars as done in similar figures.
where ω is the differential 1-form such that u = ω \ .Since the material derivative for this 1-form is we find that the momentum equation can be written as D u ω + ϵδdω + dp = 0, where p = − 1 2 u • u + p.

Fig. 1 :
Fig. 1: Illustration of small edges (a) and corresponding local shape functions (b) for the unit triangle.
s.t.F ⊂ ∂K and K ̸ = K old (K on the other side of face F )

Fig. 7 :
Fig.7: Energy of the discrete and exact solution for Experiment 1 using the firstand second-order, energy-tracking schemes on a simplicial mesh with mesh width h = 0.0949795, timestep τ = 0.00625.
Velocity field for Experiment 3 computed using the second-order, conservative semi-Lagrangian scheme on a simplicial mesh with mesh width h = 0.189959 and time-step τ = 0.0125.The colors indicate the magnitude of the vector.

Fig. 17 :
Fig. 17: Convergence results for Experiment 4 using the first-and second-order schemes without energy-tracking on simplicial meshes with mesh width h, timestep τ = 1 √ 2 h.