A review of some geometric integrators

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. A review of some geometric integrators Dina Razafindralandy, Aziz Hamdouni, Marx Chhay


Introduction
With the increasing performance of computers (speed, parallel processing, storage capacity, …), one might think that it is not worth to design completely new algorithms to solve numerically basic problems in mechanics. It is tempting to think that to simulate physics with a fair precision, one just has to take small enough time and space steps. Yet, experiments show that, even with an academic problem such as an harmonic oscillator, classical numerical schemes which rely on a direct discretization of the equations are unable to predict correctly the solution over a long time period, even with a fine time and space grids. In fact, the equation of a system may hide physically very important properties, such as conservation laws, which are destroyed by these schemes, leading to a meaningless prediction or a blow up.
Many physical properties of a system are encoded within a geometric structure of the equation. A natural way to correctly predict these properties is then to preserve exactly the geometric structure during the simulation. This is the foundation principle of geometric integrators. This approach of discretization generally enables a better representation of the physics of the system and a particular robustness for a long time or/and large space simulation.
The aim of this article is to draw the attention of computational mechanics specialists to geometric integrators. However, since many geometric structures, and consequently many geometric integrators, exist, it is not conceivable to describe all of them. Instead, we will focus on the most widely used. The main properties of these integrators will be listed. Examples of construction will be presented and illustration on mechanical systems will be given. Many references will be provided throughout the article for each part of the article for more complete details.
One of the oldest geometric structures used in physics is the symplecticity of Hamiltonian systems. When it can be written in a symplectic framework, an ordinary differential equation (ODE) exhibits a preservation of a differential form, the symplectic form, over time. The symplectic formulation also permits to obtain conservation laws via Noether's theorem [1][2][3]. Symplectic integrators are built to preserve the symplecticity of the flow at the discrete level. One of the first works on symplectic integrators is that of Vogeleare in 1956 (see [4]), followed by many papers and books in the 1980's [5][6][7][8][9][10][11]. Since then, an abundant literature can be found on the subject.
Many symplectic systems come from a variational problem. In this case, the equation represents a trajectory which minimizes a Lagrangian action. A way of building a symplectic integrator for such a problem is to discretize the Lagrangian function and solve the corresponding discrete variational problem. This approach yields variational integrators [4,[12][13][14]. Variational integrators are automatically symplectic and momentum preserving. They have been developed since the 1960's in optimal control theory and the 1970's in mechanics ( [15][16][17][18][19][20][21][22], see [13] for a more complete historical overview). Note that variational integrators can be extended to solve efficiently non-variational problems by embeeding the latter into a larger Lagrangian system [23,24].
Other equations of mechanics cannot straightforwardly formulated with a Lagrangian functional nor a (multi-)symplectic structure. It is the case of many mechanical dissipative systems. Yet, many of them have important invariance properties, called symmetries [54], under some transformations. Symmetries play a fundamental role since they may encode exact model solutions (self-similar, vortex, shock solutions, …), conservation laws via Nother's theorem or physical principles (Galilean invariance, scale invariance, …) [55][56][57][58][59][60][61][62][63][64][65][66][67][68]. They have also been used for modelling purposes such as the establishment of wall laws in turbulent flows or the development of turbulence models [69][70][71][72][73][74][75][76]. It is then important that numerical schemes do not break symmetries if one wishes to reproduce numerically the cited model solutions and properties. Preserving the symmetries of the equations at the discrete level is the founding key of invariant integrators [77][78][79][80][81][82][83], as will be seen in "Invariant integrators" section. Note that the symmetry structure does not exclude simplecticity structure. So, combining symplectic/variational algorithm and invariant scheme is feasible, but not done in this article.
In the next section, we recall the basic principle of symplectic geometry. We then show why and how to build symplectic and multisymplectic integrators. In "Variational integrators" section, variational systems and variational integrators are presented. A total variation approach, which extends the usual presentation of Lagrangian mechanics, is used. Indeed, with a total variation consideration, an energy equation (and momentum equation in case of PDE) naturally arises. This equation can be used to build, for instance, energy-preserving schemes. Invariant integrators are dealt with in "Invariant integrators" section. The approach used is that of Kim [79,133] and Chhay and Hamdouni [80] which consists in modifying classical schemes to make them invariant. At last, discrete exterior calculus is described in "Discrete exterior calculus" section. Simulations of fluid flows convecting a passive scalar are presented.
Some theoretical tools of differential geometry are used in this article to highlight the geometry structures. Their definition can be found, for example, in [1,54,[134][135][136][137][138][139][140][141][142][143]. However, their use has been limited (at the risk of using formal definitions) and the examples are simple enough for readers more interested in the computational aspect.
the symbol ¬ standing for the interior product 1 and dH is the differential or exterior derivative of H. A Hamiltonian system on (S, ω), associated to the Hamiltonian function H, is a dynamical system s(t) governed by the equation Let t be the Hamiltonian flow, i.e. the map t : which, to any initial condition s 0 , associates the solution at time t ≥ 0. The Hamiltonian flow t is a symplectomorphism, meaning that along trajectories, * t ω = ω (4) where * t is the pull-back 2 of t . In other words, the Hamiltonian flow preserves the symplectic structure of the equation.
The symplectic manifold S has necessarily an even dimension, say 2n q , and there exists canonical coordinates s = (q, p) T belonging to an open set of R n q × R n q in which Eqs. (1) and (2) can be written as follows: 1 The interior product or contraction of a vector field X and a differential form μ is often noted if μ is a 2-form and Y any vector field.
the Hamiltonian H being a function (t, q, p) ∈ R × R n q × R n q → H(t, q, p) ∈ R. Equation (5) can be written in a more compact way as follows: where ∇H = (∇ q H, ∇ p H) T and J is the skew-symmetric matrix J = 0 −I n q I n q 0 , I n q being the identity matrix of R n q . The matrix J is the matricial representation of the symplectic form ω. In canonical coordinates, ω writes: where ∧ is the exterior product symbol. 3 In canonical coordinates, the symplecticity property (4) reduces to in matricial form, and to dp(t) ∧ dq(t) = dp(0) ∧ dq(0) (9) in exterior calculus notation.

Flow of a numerical scheme
The flow of a numerical integration scheme of Eq. (2) is defined as n : where s n is the numerical approximation, at a discrete time t = t n , of the solution s(t) corresponding to an initial condition s 0 . A numerical integrator is said symplectic if, like the exact flow, it preserves the symplectic structure of the equation. More rigorously, a numerical scheme is symplectic if its flow verifies * n ω = ω.
The numerical pull-back * n is an approximation of the exact pull-back * t at t = t n .
In practice, we only consider one-step incremental integrators. For such integrators, it is more convenient to work with the one-step flow t : (12) instead of n . Note that t may depend on n. As examples, the flow of the explicit Euler integration scheme of (2) is The flow of a second order Runge-Kutta scheme is where And the flow of a fourth order Runge-Kutta scheme is where A sufficient condition for a one-step incremental integrator to be symplectic is that t is symplectic (for any n), that is: * In matricial and in exterior calculus forms in canonical coordinates, this condition writes: and dp n+1 ∧ dq n+1 = dp n ∧ dq n ∀n ∈ N.

Failure of classical integrators
In this section, we show, through elementary examples, that the use of non-structurepreserving algorithms may have dramatic consequences when simulating a long-time evolution problem.

Harmonic oscillator
Consider the 1D harmonic oscillator, described by the autonomous Hamiltonian This Hamiltonian is conserved along exact trajectories. Since n q = 1, the symplectic form is a volume (or area) form. The one-step flow of the first order explicit Euler scheme on the harmonic oscillator is The explicit Euler scheme is not symplectic. Indeed, for this scheme, It can easily be checked that condition (17) is not fulfilled. The symplectic form evolves between two time steps as follows: dp n+1 ∧ dq n+1 = (1 + t 2 )(dp n ∧ dq n ).
This relation shows that an infinitesimal volume in the phase space grows each time step with a rate 1 + t 2 . As a consequence, when used to simulate a long-time evolution problem, the scheme fails to represent correctly the physics of the equation. The implicit Euler scheme has a similar behaviour. Its one-step flow can be written as Like the explicit version, this scheme does not preserve the symplectic form since dp n+1 ∧ dq n+1 = 1 1 + t 2 (dp n ∧ dq n ).
An infinitesimal volume in S is multiplied by (1 + t 2 ) −1 at each iteration.
To illustrate the above properties of the explicit and implicit Euler schemes, we take a set of initial conditions forming a circle centered at (p = 1, q = 0) and with radius 0.2. The equation is solved with the two Euler schemes. The solution is plotted on Fig. 1 each π/6 s. As anticipated, the area delimited by the circle is not preserved over time. It increases drastically with the explicit scheme and decreases with the implicit one.
The growth of the initial volume can be interpreted as a numerical production of energy, whereas the diminution of the volume is associated to a numerical dissipation. Hence, through this elementary example, it was shown that the non preservation of an intrinsic property of the system yields an accumulation of numerical errors, modfying the initial area and contradicting Liouville's theorem.  In fact, the discrete energy H n is not constant and evolves numerically as with the explicit scheme and as with the implicit scheme. Figure 2 shows the evolution of the absolute relative error with an initial condition (q = 0, p = 1). As can be observed, the absolute error increases for both schemes. A refinement of the time step reduces the growth rate but does not avoid the problem. Even with t = 2.10 −2 , the relative error is 10% at t = 5 s. The use of such schemes on long-time evolution problem may be catastrophic.

Kepler's problem
As a second example, consider Kepler's problem [146,148,149], described by the Hamiltonian function This Hamiltonian being autonomous, it is an invariant of the system. The (third component of the) angular momentum is also a first integral. The equation, with initial conditions (q 1 = 1/2, q 2 = 0, p 1 = 0, p 2 = √ 3) is solved with the second order Runge-Kutta scheme. The flow is described by Eq. (14). It can be shown that this scheme is not symplectic.
The absolute relative error on the Hamiltonian and the angular momentum are presented in Fig. 3. As can be observed, the numerical error on both invariants has an increasing trend over time. A time-step refinement reduces the error but the growth rate is essentially the same whatever the time step. Moreover, increasing the order of the scheme is not the solution. Indeed, with a 4th-order Runge-Kutta integrator, defined in Eq. (15), the trend stays similar, as will be exhibited in "Kepler's problem" section.
Other examples, including molecular dynamics and population evolution problems, which show the failure of classical methods are presented in [150]. In the next section, we present some symplectic schemes and show how they are suitable for the resolution of long-time evolution problem having a Hamiltonian structure.

Some symplectic integrators
There are some ways to build symplectic integrators [151][152][153][154][155]. One of them is through a generating function. Another way is to adapt existing integrators. In the present paper, we adopt the latter method. Two types of integrators are investigated, those based on Euler schemes and those belonging to the Runge-Kutta method family.

Euler-type symplectic integrators
The classical explicit and implicit schemes on a general Hamiltonian Eq. (5) can be written, respectively, as follows: By mixing these schemes, a symplectic integrator can be deduced as follows [151]: Indeed, using the relations dp n+1 = dp n − t ∂∇ q H ∂p n+1 dp n+1 + ∂∇ q H ∂q n dp n , a straightforward computation shows that scheme (24), called symplectic Euler scheme, is symplectic. Another symplectic integrator is the centered Euler scheme obtained with the mid-point method: with p n+1/2 = p n+1 + p n 2 , q n+1/2 = q n+1 + q n 2 .

Symplectic Runge-Kutta methods
An s-stage Runge-Kutta integrator of Eq. (5) can be formulated as follows [9,151]: with . . , s, for some real numbers α ij , β i , i, j = 1, . . . , s. It was proven that the s-stage Runge-Kutta (26) is symplectic if and only if the coefficients verify the relation [10,11]: For example, when s = 1, α 1 = 1 2 and β 1 = 1, scheme (26) reduces to the centered symplectic Euler scheme. With s = 2, we have a 4th-order Gauss-Legendre scheme with Another 4th-order scheme, but with a lower triangular matrix α, is the 3-stage, diagonally implicit, 4th-order Runge-Kutta, defined by Many other variants of symplectic Runge-Kutta methods can be found in the literature (see for example [152]). Some of them are suited only to particular Hamiltonian systems. Generally, symplectic integrators do not preserve exactly the Hamiltonian [156] (when the latter is autonomous). However, the symplecticity condition seems to be strong enough that, experimentally, symplectic integrators exhibit a good behaviour toward the preservation property. In fact, we have the following error estimation on the Hamiltonian [157,158]: for some constant γ > 0, r being the order of the symplectic scheme. This relation states that the error is bounded over an exponentially long discrete time. Moreover, it was shown in [159] that symplectic Runge-Kutta methods preserve exactly quadratic invariants.
In fact, a backward analysis shows that symplectic integrators solve exactly a Hamiltonian system, with a Hamiltonian functionH which is a perturbation of H [160]. This explains the bounded error on H.

Harmonic oscillator
As a first test, we solve the harmonic oscillator problem with the symplectic Euler scheme and with its centered version. By construction, both schemes preserve the symplectic form up to machine precision. So, any area in the phase space is preserved. This property can be observed in Fig. 4 which shows the evolution of the disc delimited by the same circle as in "Harmonic oscillator" section. Indeed, even if, contrarily to the exact evolution, the disc is deformed by the Euler symplectic scheme, its area is invariant. With the centered symplectic scheme which is second order accurate, the form of the disc is unmodified and its area is invariant. Figure 5 presents the evolution of the relative error on the Hamiltonian H. As can be observed on the left of this figure, H is not preserved by the non-centered symplectic Euler scheme. However, the error is bounded, almost periodic, contrarily to the result of the classical, non-symplectic Euler schemes in "Harmonic oscillator" section. Moreover, a closer comparison with Fig. 2 shows that the error is much smaller with the symplectic scheme.
With the centered symplectic scheme, the error on H is below the machine precision (Fig. 5b). In fact, it can be shown that, since this scheme is a Runge-Kutta one, it preserves H exactly.

Kepler's problem
As a second test, we consider again Kepler's problem. The equation is solved with the classical 4th-order Runge-Kutta scheme (RK4) and its symplectic version (RK4sym). The time evolution of the error on H and are plotted in Fig. 6. As seen in subfigures (a) and (c), the error of RK4 on H and is smaller than that of the second order scheme in "Kepler's problem", but both schemes have similar trend when time grows. The error increases without upper limit.
As for the symplectic scheme, the error of RK4sym on H stays almost the same over a long time period, as could be predicted with estimation (29).
The above results show that symplectic integrators preserve many first integrals (in the sense that the error is bounded). However, not necessarily all first integrals are preserved. For example, standard symplectic integrators do not necessarily preserve the Laplace-Runge-Lenz vector. Recall that this vector is an invariant of Kepler's problem which makes the equation super-integrable (see [161]). It is possible to build specific schemes which preserve this invariant. It is done in [162] using the canonical Levi-Civita transformation [163].

Vortex dynamic
To end this section, a qualitative study of the convection of passive point vortices in an ideal incompressible fluid is presented [164][165][166][167]. A system of N point vortices or tracers is described by a singular vorticity field with the conjugate variables p i = γ i x i and q i = y i . We take N = 4 identical vortices. In this case, some typical behaviours have to be observed according to the initial configuration of the vorticies. When the vortices are initially positionned with a square symmetry as in Fig. 7, left, then the trajectory of each particle is stable and periodic (Fig. 7, right). In fact, the system is integrable and the configuration of the vortices remains symmetric at any time.
If the initial position of two opposite vortices are shifted from the square configuration with a small angle α as in Fig. 8 then the system is quasi-integrable. Each particle has a quasi-periodic trajectory. But if only one vortex (instead of two opposite ones) is shifted from the square configuration, symmetries are lost and the trajectory becomes chaotic (see [168][169][170]).
We now solve the system numerically, with RK4 and its symplectic version RK4sym. When the initial configuration is square symmetric, both integrators reproduce the periodic behaviour of the solution, even over a long time period, as can be noticed in Fig. 9. In this figure, the approximate solution is plotted over 5.10 5 iterations with a time step t = 5.10 −3 . When the initial position of two opposite vortices are perturbed as in Fig. 8, left, with an angle α = π/5, the violation of the symplectic structure of the equations makes the classical Runge-Kutta scheme unable to reproduce the physics. The quasi-periodicity is indeed lost and the trajectory seems chaotic, as shown in Fig. 10. As for it, RK4sym reproduces the quasi-periodicity of the solution thanks to its symplecticity property.
These results show the importance of the preservation of the equation's structure when simulating Hamiltonian ODEs. In order to extend symplectic integrators to PDEs, multisymplectic geometry will be introduced in the next subsection.

Multisymplectic integrators
Mulstisymplectic integrators have been intensively developped for PDE, specifically in the framework of water waves [47,171]. Indeed, a large class of models for water waves inherits a Hamiltonian structure of infinite dimension. Thus a multisymplectic structure can be exhibited, for example, for Serre-type equations in deep water configuration [172], and for the Serre-Green-Naghdi equations in shallow water configuration [173]. Therefore multisymplectic schemes appear as natural structure preserving integrators applied to these models. The efficiency of such geometric numerical methods is presented in [174] where a long-time simulation of the Korteweg-de Vries dynamics is performed with more robustness and more accuracy using a multisymplectic scheme than a Fourier-type pseudo-spectral method.
A precise presentation of the multisymplectic structure on a manifold necessitates the introduction of many mathematical tools. As our article aims to be a review paper, we do not wish to do it here. The reader can refer to [175, §1.1-1.3] for a good mathematical approach. Instead, we propose a simplified framework, in a vector space.
Consider be an open subset U of R 2 and a regular function s : with n s ≥ 2. Let H : U × R n s −→ R be a regular Hamiltonian function.
In most of applications, a multisymplectic system is a system of partial differential equations of the following form: In (32), W and K are skew-symmetric matrices in R n s and the gradient of H is relative to s. Equation (32) is a generalisation of the (one-)symplectic Eq. (6).
To W and K can be associated two skew-symmetric two-forms ω and κ defined by 4

Conservation laws
There are conservation laws of Eq. (32) that should be considered when designing a numerical integrator for the equation. First, in the (one-)symplectic case, the two-form is conserved along trajectories. In the 2-symplectic case, it is a sum of the variation in each direction which is conserved. More precisely, it can be shown that [171]: This is a consequence of a linearization of (32) and the symmetry of the Hessian of H. Furthermore, in the (one-)symplectic case, the Hamiltonian is a first integral, provided it does not depend explicitely on time. In multisymplectic case, this conservation law is split and we have one conservation law in each direction. More precisely, if H does not depend explicitly on t then one can deduce an energy conservation law: 4 More precisely, if the W ij are the components of W. Moreover, ω(X 1 , X 2 ) = (WX 1 )·X 2 for two arbitrary vectors X 1 and X 2 of R ns , the dot symbol · being the Euclidean inner product of R ns . Similar remarks apply to κ.
where F 1 and F 2 are the functions Conservation law (35)

Examples
As an example, consider the following Hamiltonian function with s = (u, v, w): Let ω = dv ∧ du and κ = du ∧ dw. The corresponding matrix expressions are where V is a source function. This system is equivalent to the wave equation The quantities involved in conservation laws (35) and (36) are Another multisymplectic equation is the Korteweg de-Vries equation then Eq. (32) is equivalent to Korteweg-de Vries Eq. (38). The matrix representations of ω and κ are The conserved quantities are defined by Lastly, the non-linear Schrödinger equation and ω = dp ∧ dq, κ = dv ∧ dp + dw ∧ dq.
A multisymplectic integrator is a discretization scheme of (32) which satisfies the corresponding discrete form of the multisymplectic conservation law (34) [176]. As for ODEs, many classical PDE integrators have their multisymplectic version. In this subsection, some examples of multisymplectic integrators are presented. To this aim, the expression of the Hamilton equation in local coordinates (32) is used. As for the symplectic case, Euler-type and Runge-Kutta-type schemes are considered.

Euler-type multisymplectic schemes
To build a multisymplectic Euler scheme, the skew-symmetric matrices W and K are split into [177]: For example, the splitting may be in upper and lower triangular parts. Next, define a grid formed with points (t n , x i ). Recall the forward and backward Euler discretization schemes in each direction: One multisymplectic Euler discretization of (32) is defined by It is first order in time and second in x. Any solution of this equation verifies the following discrete multisymplectic conservation law This scheme also verifies semi-discrete forms of energy and momentum conservation laws (35) and (36) [170]. Another multisymplectic discretization scheme of (32) is the Preissman box scheme This scheme is second order in both t and x, and verifies the discrete multisymplicity conservation + t ω n Lastly, the multisymplectic midpoint scheme is defined by and satisfies the discrete multisymplicticity property

Multisymplectic Runge-Kutta schemes
It was shown in "Symplectic Runge-Kutta methods" section that a symplectic Runge-Kutta scheme can be obtained with a rather simple condition of the coefficients in the Butcher tableau which guarantees the symplecticity. However, no extension has been established yet in the generic multisymplectic case. Multisymplectic RK schemes were presented and studied in [178][179][180] for the partitioned case. Another particular RK scheme is the implicit Gauss-Legendre integrator [181]. It will be illustrated hereafter on the Sine-Gordon equation. Other multisymplectic schemes, based on the Runge-Kutta-Nyström method, can be found in [182][183][184].
Consider the multisymplectic formulation of Eq. (37), discretized with some discrete time and space derivative operators t and x :

. A multisymplectic Gauss-Legendre discretization scheme of this equation is a combination of space and time Runge-Kutta integrators with (s
..,s w as follows: The multisymplecticity condition on the coefficients of this RK scheme is similar to the symplectic case for our scheme: The discrete conservation law reads:

Numerical test
Consider the Sine-Gordon equation Equation (39) is discretized with the following three different schemes. • A leapfrog (LF) scheme: This scheme is in fact a symplectic (not multisymplectic) scheme in the sense that it preserves a spatial symplectic two-form over time. • An energy conserving but not multisymplectic scheme developped in [185] that we call EC: • A nine-point box multisymplectic (MS) scheme, which is a Runge-Kutta scheme, simplified by variable substitution [186]: where the time and space discretization operators are defined by When the space step is small enough and the time step is smaller, the three schemes reproduce the exact solution with a fair precision. However, when the time step grows, the Leapfrog integrator looses stability, despite its symplecticity in one direction. This behaviour can be seen in Fig. 11a, where t = x. Figures 11b and 12 show that, with the same time and space steps, the energy conserving and the multisymplectic integrators remain stable, even over a long time period.
When the time step is bigger than the space step, the scheme EC based on energy conservation blows up from the first iterations (the corresponding approximate solution is not plotted here). As for it, the multisymplectic scheme gives a relatively well behaved solution, even with a coarse time grid, as can be seen in Fig. 13.
Some results on the quality of the multisymplectic scheme regarding conservation laws are given in [170]. It can be concluded that symplectic and multisymplectic schemes are particularly suitable to the numerical resolution of long time evolution problems.  In the next section, geometric integrators for variational ODEs and PDEs are presented.

Variational integrators
In this section, we deal with Lagrangian systems, that are systems coming from a calculus of variation.

Reminders on Lagrangian mechanics
In Lagrangian mechanics, a Lagrangian system is described by a Lagrangian density. Taking the variation of the corresponding Lagrangian action over the configuration variable q, the Euler-Lagrange equation of the system is deduced from Hamilton's principle of least action. This is the traditional approach of presenting the Euler-Lagrange equation [1,187,188].
In a similar way, variational integrators are obtained by taking the variation of a discrete Lagrangian action over the discrete variable q n . Abundant literature on the traditional presentation of variational integrators exists [13,14,150].
Generally, the numerical solution of the discrete Euler-Lagrange equation preserves the evolution law of energy of the system with a good precision, but not exactly. In [189], Kane et al. proposed to associate an ad hoc equation to the discrete equations in order to satisfy a discrete energy conservation.
In fact, the algorithm obtained in [189] can be viewed as a generalisation of variational calculus, where variation along time is permitted [170,190]. In this new approach, the energy equation is not defined as in [189] but results naturally from the variation of the action in time direction. In the present paper, this second approach is adopted.

Total variation calculus
Consider a dynamic system, which configuration over time is defined by q(t). Denote Q the configuration manifold, TQ its tangent bundle and M = R × Q the extended configuration space including time. A Lagrangian density on Q is a C 2 function defined on R × TQ: The associated Lagrangian action over a time interval [t 0 , t 1 ] is the functional L : The variation of L in the direction of a tangent vector δm = (δt, δq) ∈ T m M is: where D is the functional derivative operator and The total infinitesmal variation of L is Liouville's 1-form below emerges from the last term of (41): A system is said Lagrangian if it evolves under Hamilton's principle of stationary action, stating that the trajectory corresponds to an extremum of L among q(t) ∈ C 2 ([t 0 , t 1 ], M) with fixed endpoints. In other words, the trajectory is a solution of This condition leads to the Euler-Lagrange equation and to the energy evolution equation Equation (45), which reduces to the energy conservation law when L is autonomous, is automatically verified by any solution of Eq. (44). But numerically, this is not always the case.
In the standard way of deducing the Euler-Lagrange equation, only tangent vectors δq ∈ T q Q are considered. However, by considering the time as a variable, and thus the tangent vectors (δt, δq), conservation laws, as defined in Noether's theorem [191][192][193][194], appear naturally.

Noether's theorem
Consider a curve in M defined by the one-parameter transformation: passing at (t, q) when a = 0 (see Fig. 14). Take (δt, δq) as the tangent vector to this curve at (t, q): First Noether's theorem states that if g a preserves the Lagrangian action L, that is for any a belonging to a vicinity of 0 then the following conservation law holds along trajectories of the system: This equation is obtained from δL = 0 in (41) along a solution trajectory without imposing the fix endpoint condition (43). The flow t L : where ω L = d L is the Lagrangian symplectic form. The form ω L endows T M with a symplectic structure.

Variational integrators
A variational integrator is a scheme which verifies a discrete Hamiltonian's variational principle. The most popular way of obtaining a variational integrator is from a discretization of the Lagrangian density and a derivation of a discrete version of Eqs. (44), (45) by mimicking the variational procedure used in "Total variation calculus" section. This is the approach of Marsden and West [13,14], based on a finite difference discretization of the Lagrangian density. Other approaches can be found, for instance, in [22,[195][196][197][198][199][200]. As examples, a higher-order, spectral method is presented in [22]. A finite element method is described in [195]. Another approach, where at each time iteration an exact variational problem is solved, is developped in [198,201]. Fig. 15).
With a suitable discretization ofq we can get a discrete Lagrangian function L(m). And with a suitable quadrature method, a discrete Lagrangian action L over [t 0 , t N ] is deduced: L is a discrete function M N +1 → R. Its variation reads for some functions E n t and EL n and some 1-forms − L and + L which depend on the discretization scheme of the Lagrangian density and the quadrature. Imposing that (D L, δm) = 0 (51) for any δm such that δm 0 = δm N = 0 generates the discrete Euler-Lagrange equations EL n (m) = 0, n = 1, . . . , N − 1 (52) and the discrete energy equations Let us illustrate the above discrete variational equations through explicit examples.

Rectangle rule
A simple example of a variational integrator is obtained with a first order finite difference discretization ofq and a rectangle rule quadrature method of the action. In this case, the discrete Lagrangian density is The discrete action is: (t n+1 − t n )L (t n , q n ,q n ) .

Midpoint rule
If a midpoint rule is used as quadrature method then the discrete Lagrangian action can be written: In this case, the discrete variational equations are In Eqs. (55a), (55b) and (56a), (56b), the unknowns are t n+1 and q n+1 . The time step adapts itself to satisfy both the Euler-Lagrange and the energy equations. Some other variational integrators, such as Newmark scheme or symplectic Runge-Kutta written for separated Lagrangian, can be found in [13,150].

Symplecticity of a variational integrator
A basic stencil of a variational integrator like in the previous example involves three points: m n−1 , m n and m n+1 . The numerical one step flow is the map As seen, two different discrete forms − L and + L appear from the discretization of a variational problem. However, their exterior derivatives are equal up to a sign. If we call then it can be shown that the numerical flow is symplectic and preserves the 2-form ω L . Note that a discrete Noether's theorem can also be derived.

Numerical test
Consider the 1D harmonic oscillator, described by the Lagrangian Hamilton's principle yields the Euler-Lagrange equation Since the Lagrangian is autonomous, Eq. (45) reduces to the energy conservation equation: which is automatically verified by any solution of (57). The Lagrangian density and the action are discretized as in section . In this case, Eqs. (55a), (55b) become In contrast to the continuous case, these two equations are independent. The first one is the discrete Euler-Lagrange equation. The second expresses the conservation of the discrete energy Observing the recurrence relation, (58b) can be written: Knowing (t n−1 , q n−1 ) and (t n , q n ), Eqs. (58a) and (59) are solved for (t n+1 , q n+1 ). For the test, the initial conditions are choosen such that the analytic solution is a sine function. The approximate solution is very close to the exact one, as presented in Fig. 16. Figure 17 shows that the error on the energy is below the machine precision. This is due to Eq. (59) which ensures that H n = H 0 at any discrete time.

Variational integrator on PDE
So far, we only considered variational integrator for ODEs. As we did for Hamiltonian mechanics, the Lagrangian approach can be extended to PDEs.

Lagrangian approach of PDE
Consider a Lagrangian function L(z, u, u z )  where z belongs to an orientable base manifold Z with boundary ∂Z, u lies in a configuration manifold U and u z designates the partial derivatives of u. Denote M = Z × U. 5 In a simplified way, the Lagrangian action can be viewed as L : where C 2 (Z, M) is the space of possible C 2 trajectories defined from Z to M.
To simplify, we assume a two-dimensional base space, that is z = (t, x), and a onedimensional configuration space (dim M = 3). The components of u z are denoted u t and u x . As done previously, the total variation of the action is (see [170,202] Hence, when we impose the Hamilton's principle of least action, under the condition that (δt, δx, δu) = (0, 0, 0) on ∂Z, we get the so-called variational equations Equation (61a) is the Euler-Lagrange equation. Equations (61b) and (61c) are respectively called the energy and the momentum evolution equations. When L does not explicitly depend on t and x, Eqs. (61b) and (61c) reduce respectively to an energy and to a momentum conservation laws. As in the case of ODE, these two equations are automatically verified by any solution of the Euler-Lagrange equation. But, this is generally not true when passing to discrete scale with a classical numerical scheme. Noether's theorem can also be extended to Lagrangian PDEs as follows. If the Lagrangian action is invariant under a transformation t, x, u, a),x(t, x, u, a),û(t, x, u, a)) then we have a conservation law: Indeed, the corresponding Euler-Lagrange equation is the wave equation The energy and momentum equations read Since the Lagrangian depends explicitly neither on t nor on x, Eqs. (64) and (65) could also be obtained from Noether's theorem.

Example of variational integrator
The domain Z is discretized into grid points (t n , x n j ), n = 0, . . . , N and j = 0, . . . , J (see Fig. 18). The maximum space index J may depend on n.
The dependent variable is discretized into (u n i ). As in the case of ODE, a variational integrator results from a choice of a discretization of the derivatives in the Lagrangian density and a quadrature of the integral in the action. One of the simpliest choices is rectangle rules for quadrature: The variation of the action is After expliciting δL( n j ), shifting summation indices and imposing Hamilton's principle, we get the following discrete variational equations: In these equations, the shortcuts In this section, we presented geometric integrators for ODEs and PDEs having a symplectic structure, coming from a variational problem. In the next section, we consider more general equations which may not have any symplectic structure but a Lie symmetry group. We then show how to construct geometric integrator for such equations.

Invariant integrators
We first set some theoretical background, by defining symmetry of an equation and precising the notion of invariance.

Symmetry group
Consider again the manifold M introduced in section and a partial differential equation is called a symmetry of Eq. (68) if it leaves the solution manifold E unchanged, that is

When (71) holds, Eq. (68) is said invariant under g.
A set G of transformations, which has a group (respectively a Lie-group [54,203]) structure, is called a transformation (resp. Lie transformation) group. And if each element of G is a symmetry of an equation, then G is called a symmetry (resp. Lie symmetry) group of that equation. We are particularly interested in one-parameter Lie groups, i.e. Lie groups of transformations g a : (z, u) → ẑ (z, u, a),û(z, u, a) , which depend continuously on a parameter a.
As an example, consider the heat equation with z = (t, x) and u = u, and the scaling transformations where a is a real parameter. It is straightforward to show that transformations (74) form a Lie symmetry group G of the heat Eq. (73). Indeed, So, if u(t, x) is a solution of (73) then u(e 2a t, e a x) is also a solution for any a ∈ R. Moreover, G has a manifold structure of dimension 1, with local coordinate a. It is also easy to show that the Lie symmetry group G of scaling transformations (74) acts freely and regularly on the solution manifold M of Eq. (73). Besides the notion of symmetry of an equation, we will also need the concept of symmetry of a function. A function F (z, u) is said invariant under a transformation g if F (g(z, u)) = F (z, u). (76) Note that if the function E in (68) is invariant under a transformation g, then g is a symmetry of equation (68), but the converse does not hold. Indeed, a symmetry of Eq. (68) modifies the function E in a multiplicative way:

E(g(z, u)) = e(z, u) E(z, u)
for some function e depending on g.
Computing the Lie symmetry groups of an equation is often a tremendous task. Fortunately, it can be made algorithmic with use of infinitesimal generators [54], such that many computer algebra packages can be used [204][205][206][207]. Some of them even compute conservation laws.
The knowledge of symmetries of an equation gives precious information on the equation and on the physical phenomenon it modelises. For example, from one known solution, symmetries enable to find other solutions. Symmetries may also be used to lower the order of the equation and to compute self-similar solutions [54]. More fundamentaly, as stated by Noether's theorem [191][192][193], symmetries are linked to conservation laws.
As introduced, preserving the symmetry group through the discretization process is necessary if one wishes not to loose self-similar solutions and conservation laws during simulations. A way of building integrators that are compatible with, or invariant under, the symmetry group is to make use of independent differential invariants of the equation [77]. However, combining these invariants into a numerically stable scheme is rather complicated. Instead, we follow the idea in [79][80][81], which consists in modifying classical schemes so as to make them invariant under the symmetry group. To show how to do this, we need to formalize the notion of a discretization.

Numerical scheme
A finite difference discretization of the domain Z is a network z = (z 1 , ..., z J ) ⊂ Z of discrete points z j linked by a relation The mesh is defined by the function . Similarly, the unknown variable u is discretized into a sequence where u l is to be understood as an approximate value of u at z l . We denote m j = (z j , u j ) and m = (m 1 , ..., m J ) ⊂ M. A discretizetion scheme of equation (68), with an accuracy order (q 1 , ..., q n z ), is a couple of functions (N, ) such that as soon as that is, as soon as each m j belongs to the solution manifold. In Eq. (78), x i is the step size in the i-th direction, i = 1, ..., n z . In (79), has been extended to m such that N and has the same arguments. This extention is also necessary when the mesh changes with the values of u (adaptative mesh, ...).
As an example, consider again the heat Eq. (73). The domain is descretized into z = (t n j , x n j ) n,j as shown in Fig. 19. An orthogonal (cartesian) mesh on Z is defined by This mesh is regular if Relations (81) and (82) define the function corresponding to an orthogonal and regular mesh. Note that when the mesh is orthogonal, we can denote t n j = t n and x n j = x j thanks to (81).
The unknown u is approximated by a discrete sequence (u n j ) n,j . The Euler explicit descretization scheme of the heat Eq. (73) on an orthogonal and regular mesh is then defined by: Note that conditions (84) but conditions (85) are not required. Starting from any existing scheme (N, ), our aim is to derive a new scheme (Ñ ,˜ ) which is invariant under the symmetry group of the equation. This will be done using the concept of moving frame.

Invariantization by moving frame
Consider a Lie group G acting on M. In applications, G is a Lie symmetry group of the equation. We indicate the group action with a centered dot (·) so that g·m = g(m) for any (g, m) ∈ G × M. Assume that the action is regular and free. A (right) moving frame on M relative to G is a map ρ : M → G verifying the equivariance condition ( [135,208]): This is a generalization of Cartan's moving frame (repère mobile) definition [135,209]. If ρ(m) is seen as the frame at m then condition (86) means that the frame at m and at g·m are linked together with a right translation R g −1 (see Fig. 20). A direct and important consequence of the equivariance condition (86) is that if ρ is a moving frame relative to G then Relation (87) is illustrated in Fig. 21.  (87) is a moving frame relative to G. Indeed, for any m ∈ M and any t a ∈ G t , Note that the moving frame is not unique. Indeed, ρ m 0 [m] = t −m+m 0 is a moving frame for any m 0 ∈ M. Olver proposed a general way of constructing moving frames via cross sections [135,208]. In its approach, fixing the arbitrary constant is equivalent to choosing a cross section. Another example is the scaling symmetry group G s = {s a : (t, x, u) → (a 2 t, ax, u/a), a > 0} acting on the solution manifold of Burgers'equation For any a > 0, the map: is a moving frame relative to G s at any (t, x, u) ∈ M where x = 0. With a moving frame, an invariant integrator can simply be derived as follows. Consider a numerical discretization scheme (N, ) of an equation and a transformation group G. A fundamental theorem [208] shows that if ρ is a moving frame relative to G then the discretization scheme (Ñ ,˜ ) defined bỹ is an invariant (under G) numerical scheme of the same order for the same equation. Let us illustrate the invariantization process (91) on Burgers' equation.

Application
Consider again Burgers'equation (89). The symmetry transformations of (89) are: • time translations: • space translations: • scaling transformations: • projections: • and Galilean boosts: Consider the usual Euler forward-time centered-space (FTCS) scheme: on the orthogonal and regular mesh (81), (82). The scheme (N, ) defined by relations (97), (81) and (82) is invariant under time translation, space translation and scaling transformation groups, like most of standard schemes on Burgers' equation. But (N, ) is not invariant under the Galilean boost and projection groups. We need to apply the invariantization process only under these two groups. However, for convenience, we take also into account the time and space translation groups. These groups can be gathered into the four-parameter group G 0 which generic element is defined by: 1 − a 4 (t + a 1 ) u = (1 − a 4 (t + a 1 ))u + (x + a 2 )a 4 + a 5 .
Each transformation g i is obtained from g 0 by setting a j = 0 for all j = i. A moving frame ρ associated to g 0 is an element of G 0 . This means that ρ[m]·m is of the form (98), with particular values of the parameters a i , depending on m. Determining ρ[m] is then equivalent to deciding the values of the a i 's.

Transformation of the grid
A basic stencil for the FTCS scheme is represented by: where z n j = (t n j , x n j ). At each point m n j , we choose the moving frame such that In this way, the transformed scheme does not depend explicitly on x m i 's and t m i 's but on the step sizes: Indeed, if we denoteẑ m i = ρ[m n j ]·z m i the transformed stencil, the choice (100) of a 1 and a 2 leads to: , , , With an orthogonal and regular mesh, defined by (81), (82), we have: The transformed grid stays regular in space becausê In fact,ĥ = h. Moreover, time layers stay flat as for the original mesh grid sincê t n j+1 −t n j = 0. But the time step size is modified: This induces a translationσ n j =x n+1 j −x n j of the spatial layers at each time increment, with:σ as seen in Fig. 22. The new grid˜ is defined by (105)- (108). The next step is to defineÑ .

Invariantization of the scheme
To apply the invariantization process (91), we need to calculate the imageû m k = ρ[m n j ]·u m k , for allû m k appearing in the stencil. We have: u = (1 − a 4 (t + a 1 ))u + (x + a 2 )a 4 + a 5 .
With the previous choice of a 1 and a 2 and the regularity of the mesh, it follows: The new schemeÑ is then defined by We now have to specify the values of a 4 and a 5 to make the transformed scheme invariant under the group G 0 .

Determination of a 4 and a 5
With the previous choices a 1 = −t n j and a 2 = −x n j , the group element g 0 reduces to: , The moving frame ρ[m n j ] has the same form as g 0 in (112) but with different values a 4 and a 5 of a 4 and a 5 . The scheme is invariant only if a 4 and a 5 are such that the equivariance condition (87) (113) for any m m i belonging to the stencil. Condition (113) is illustrated in Fig. 23. This figure is an application of Fig. 21 to the stencil. u is not represented.
Since g, ρ[m n j ] and ρ[g·m n j ] all have the same form (112) but with different values of the parameters a 4 and a 5 , we have to call the parameters of each of them differently in order to avoid confusion. To this aim, we keep a 4 and a 5 for ρ[m n j ]; we denote μ and η the parameters of g, and a 4 and a 5 those of ρ[g·m n j ]. We have, on the one hand, from (104) and (110): and, on the other hand: ρ[g·m n j ]g·m n j = 0, 0, u n j + (η +ā 5 ) , ρ[g·m n j±1 ]g·m n j = 0, h, u n j±1 ± (μ +ā 4 )h + (η +ā 5 ) .
The discrete invariance condition (113) is verified when To simplify, an algebraic expression is assumed for a 4 . And in order to preserve the explicit form of the numerical scheme, we assume that a 4 does not depend on u n+1 . Moreover, to keep the order of the convective term as in the initial scheme, a 4 must be at most of degree one. Next, the invariance of the scheme under g 4 requires that there are no term in h and k alone. Finally, as a 4 is homogeneous to the inverse of a time scale, we choose: where a, b, and c are real constants. Similar arguments for a 5 , which is homogeneous to a velocity, allow to take a general form: where d, e and f are real constants. Since a 4 and a 5 must satisfy (116), these constants verify: More information on these constants can be obtained by optimizing the order of accuracy.

Order of accuracy
A fundamental result guarantees that the invariantization of a numerical scheme using the moving frame technique preserves the consistency. However, the order of consistency may change.
As the invariantization process does not affect the diffusion term, we only need to consider the unsteady term T t and the convective T c term: corresponding to the non-viscous Burgers' equation: The consistency condition for the symmetrized scheme is: whereT t andT c are the invariantized unsteady and convection terms. A Taylor development gives: The invariantized scheme has the same order of consistency if This condition determines a 4 , which takes the expression: but gives no additional constraint on a 5 . In fact, with (123), the Taylor expansion ofT c is: This means that the invariantized convective term vanishes, and, therefore, the convective phenomenon are produced by the symmetrized unsteady term. Notice that a 5 does not appear when the invariantized numerical scheme is expressed in the regular and orthogonal original mesh. It is no longer the case if the mesh grid is not orthogonal, nor if the numerical solution is expressed in the transformed frame of reference.
We end up with some numerical results.

Numerical tests
The with ν = 5.10 −3 . We consider the standard Euler FTCS, the Lax-Wendroff and the Crank-Nicolson schemes, and their invariantized versions. The space step size is h = 2.10 −2 and the CFL number is 1/2 in the original referential. At each time step, the frame of reference is shifted by a Galilean translation It is expected that u follows this shifting according to (96). Figure 24 presents, in the original frame of reference, how the standard schemes react to this shifting. It clearly shows that when λ is high, the standard schemes produce locally important errors. In and Crank-Nicholson schemes, oscillations appear just before the sharp slope. On the contrary, the invariantized version of these schemes present no oscillation, as can be observed on Fig. 25. Consistency analysis shows that the original schemes are no longer consistent with the equation when λ = 0. This inconsistency introduces a numerical error which grows with λ, independently of the step sizes. As for them, the invariantized schemes respect the physical property of the equation and provide quasi-identical solutions when λ changes.
Another numerical test was carried out with the FTCS scheme. The exact solution is a self-similar solution under projections (95): It corresponds to a viscous shock. The shock tends to be entropic when ν becomes closer to 0. Figure 26 shows the numerical solutions obtained with the standard and the invariantized FTCS schemes at t = 2 s, with ν = 10 −2 , k = 5.10 −2 and h = 5.10 −2 . It can be observed on it that the solution given by the invariantized scheme remains close to the exact solution whereas the original FTCS scheme presents a poor performance close to the shock location.
Since the invariant scheme has the same invariance property as the analytical solution under projections, it does not produce an artificial error like the non-invariant FTCS. This shows the ability of invariantized scheme to respect the physics of the equation and the importance of preserving symmetries at discrete scales. Other numerical tests are presented in [74,81,170].
The last geometric integrator that we shall present is the discrete exterior calculus.

Discrete exterior calculus
Discrete exterior calculus (DEC) can be seen as a differential geometry and exterior calculus theory upon a discrete manifold. The primary calculus tools of DEC are (discrete) differential forms. As we shall see, they are built by duality with the grid elements. Discrete operators on differential forms (exterior derivative, Hodge star, …) are then defined in a way which mimics their continuous counterpart.
The first step to discretizing an equation with DEC is to define a grid, which plays the role of a discrete manifold. Then, each piece of the equation is replaced by its discrete versions. But since equations of mechanics are usually formulated with vector and tensor fields, they must beforehand rewritten in exterior calculus language, that is with differential forms.
One advantage of discretizing with DEC is that, by construction, Stokes'theorem is always verified exactly (up to machine precision). Another advantage is the following. Vector calculus tells us that, in a 3-dimensional manifold, curl grad f = 0 and div curl u = 0 (128) for any vector field u and any function f . However, with classical discretization schemes, there is no guarantee that these relations are verified numerically, since each differential operator is generally discretized individually. This may lead to numerical artefacts. In exterior calculus, relations such as (128) are simply particular cases of the relation where d is the exterior derivative. And, by construction of DEC, relation (129) is always verified exactly. Hence, DEC ensures that relations (128) hold exactly after discretization, even on an arbitrary coarse mesh. Some other advantages of using differential forms over tensors are discussed in [89][90][91][92][99][100][101][102][103][210][211][212]. Note that property (129) is much more general than (128). Indeed, equations (128) are meaningful only in a three-dimensional space. Moreover, equations (128) are metric dependent whereas (129) is not.
In fluid mechanics, Elcott et al. [110] used DEC to solve Euler equation on a flat and a curved surface, with exact verification of Kelvin's theorem on the preservation of circulation along a closed curve. Works on Darcy and Navier-Stokes equations were carried out in [111] and [112,213].
In the next subsections, we recall the basis of DEC and present some test cases.

Domain discretization and elements of algebraic topology
Consider a differential equation defined on a n x -dimensional spatial domain M in R n , n ≥ n x , written within the exterior calculus framework. In order to solve this equation numerically, the domain is meshed into an oriented simplicial complex K that we recall hereafter the definition and the associated algebraic topology elements (chain, cochain and boundary operator). Next, we shall see how to discretize differential forms and the exterior calculus operators.
is the equivalence class of (v 1 , v 2 , v 3 ) with respect to relation (130), the oriented face f 2 is the equivalence class of (v 2 , v 4 , v 3 ). The oriented edge e 1 is the equivalence class of (v 1 , v 2 ), etc.

Oriented simplicial simplex
A k-simplex of R n is the convex hull of k + 1 affinely independent points v 0 , . . . , v k of R n , often denoted (v 0 , . . . , v k ). It is the subset of R n determined by The domain is discretized into a simplicial complex K , that is a finite set composed of 0-, 1-, …and n x -simplices such that • for any simplex σ ∈ K , any face of σ belongs to K , • and the intersection of two simplices σ 1 , σ 2 ∈ K is either empty or a common face of σ 1 and σ 2 .
Moreover, K is assumed to be manifold-like. This means that the union of the simplices of K , treated as a subset of R n , is a C 0 manifold, with or without boundary.
Each element of K is assigned an orientation, defined as a choice of an equivalence class under the relation For k ≥ 1, this definition confers two possible orientations to any k−simplex, which can be understood as a choice of ordering of its vertices. So, to each geometric simplex (v 0 , . . . , v k ) corresponds two oriented simplices For 0-simplices (vertices), definition (130) confers only one trivial orientation. However, it is common to give also two possible orientations to vertices by assigning to each of them a sign (+ or −). Using the trivial orientation from definition (130) corresponds then to assigning the same sign + to every vertex. For simplicity, we use the trivial orientation for vertices. Figure 27 represents a sample oriented 2D mesh.

Chain and cochain
For k = 0, . . . , n x , we denote K k the set of oriented k-simplices of K and k (K ) = spanK k the vector space 6 of formal linear combinations of elements of K k , with coefficients in Z. An element of k (K ) is called a k-chain. For example, the 1-chain e 1 + e 4 + e 5 − e 3 ∈ 1 (K ) (131) in Fig. 27 constitutes a closed loop. For k = 0, . . . , n x , we denote k (K ) the algebraic dual of k (K ), that is, the space of linear maps from k (K ) to R. An element of k (K ) is called a k-cochain. Since the ksimplices of K k form a basis of k (K ), a k-cochain can simply be viewed as a map which, to each oriented k-simplex, assigns a real number: It can be represented by an array of dimension dim( k (K )) = card K k .

Boundary operator
The boundary of an oriented k-simplex is the chain of k−1 defined as the sum of the oriented (k −1)-simplices directly surrounding it. In this sum, each (k −1)-simplex is given a sign, depending on whether its orientation is consistent with that of the considered ksimplex. For example, the boundary of the oriented triangle f 1 in the example represented in Fig. 27 is More formally, the boundary of an oriented k-simplex The boundary operator ∂ defines a linear map from the vector space of the chains of K into itself. We denote ∂ k its restriction to k (K ), that is ∂ k : k (K ) → k−1 (K ). As a linear map, it can be represented by a matrix. For instance, the non-zero boundary operators on the mesh in Fig. 27 are The first column of ∂ 1 contains the components of the boundary ∂ 1 e 1 = v 2 − v 1 of the oriented edge e 1 in the basis (v 1 , v 2 , v 3 , v 4 ) of 0 (K ). The i-th column contains the components of ∂ 1 e i . In the same way, the entry of ∂ 2 at row i and column j is 0 if the edge e i does not belong to the boundary of the face f j . It is 1 (respectively −1) when e i belongs to the boundary of f j and their orientations are compatible (respectively, opposite). An important homological property of the boundary operator is: [(check on ∂ 2 and ∂ 1 in example (134)]. Indeed, the boundary of a boundary is the empty set.

Discrete differential form
We will denote k (M) the space of differential k−forms of M. A differential k-form on M defines naturally a k-cochain on K through integration over the k-chains of K . This enables to define the discretization (or reduction) of differential forms via the de Rham map: R : where Rω is the k-cochain defined on the basis K k of k (K ) by Definition (137) extends by linearity to k (K ).

Discrete exterior derivative
The discrete exterior derivative operator (that we denote d like its continuous counterpart) is defined as the coboundary operator, that is the adjoint of the boundary operator ∂ relatively to the pairing ·, · : dω, c = ω, ∂c , for any ω ∈ k (K ) and any c ∈ k+1 (K ).
As can be remarked, relation (138) is simply the expression of Stokes' theorem if the pairing ·, · is an integration. The discrete d is a linear map from k (K ) to k+1 (K ). Computationally, the restriction d k of d on the space k (K ) of discrete k-forms is represented by the transpose of the matrix of ∂ k+1 . For instance, if a discrete 0-form ω ∈ 0 (K ) takes values ω 1 and ω 2 at vertices v 1 and v 2 in Fig. 27, then And if ω is a 1-form, then Figure 28 presents an arbitrary discrete 0-form ω, its exterior derivatives dω and d 2 ω. On this example, In fact, this relation holds for any discrete k-form. This is due to the homological property (135) of the boundary operator which induces property (141) by duality.

Wedge product
The above defined discrete operators and the discrete Hodge star which will be dealt with in section . already enable to solve some interesting linear problems such as Poisson, Helmholtz or Stokes equations. However, other tools such as the discrete wedge product are also necessary in general. The definition of the discrete wedge product depends on the ranks of the involved forms and their locations (primal or dual mesh) (see [214]). As an example, the wedge product of a primal 1-form ω and a primal 0-form θ may be defined as follows at an oriented edge [v i , v j ]: If the duality ·, · is though as an integration, then the right-hand side of Eq. (142) is a linear approximation of the integral of ω ∧ θ over the edge [v i , v j ]. Another possibility of defining a discrete wedge product is to interpolate the cochain into piecewise continuous differential form (see next section), apply the continuous definition of the wedge product, and finally, discretize the result.

From cochain to differential form
We saw in "Discrete differential form" section how to discretize a differential form into a cochain. Conversely, one can reconstruct a cochain into a piecewise continuous differential form using an interpolation technique. The lowest interpolation commonly used in DEC is the Whitney map W : The Whitney k-form ϕ c is defined on a k-simplex c as follows. For a vertex v i , ϕ v i is the barycentric coordinate function. It verifies ϕ v i (v j ) = δ i j for any vertex v j . For k ≥ 1, the Whitney k-form corresponding to a k-simplex [v 0 . . . , v k ] is (see [226][227][228]): One has ϕ c i (c j ) = δ i j for two arbitrary k-simplices c i and c j ∈ K k . On a 3D mesh, the Whitney 0-, 1-and 2-forms are piecewise linear, whereas the Whitney 3-forms are constant on each tetrahedron.
The discretization and the interpolation operators verify RWω = ω for any cochain ω.
The de Rham and the Whitney maps are essential to convergence analysis (see for example [101,120,229]). Higher order Whitney forms are proposed, for example, in [230,231]. As said, Whitney maps provide a single interpolation within each tetrahedron (if the mesh is three-dimensional). One can also interpolate a cochain differently, for example by dividing each tetrahedron in many subregions and piecewisely construct an interpolated form in each subregion. This approach is used in covolume method and mimetic reconstruction [116].

Discrete Hodge star and codifferential operators
In exterior calculus framework, the expression of constitutive laws generally makes use of the Hodge star operator . Recall that the Hodge star maps a continuous k−form into its "volume complement" (n x − k)-form. It is metric dependant and verifies, at each point of M, for any orthonormal frame (u 1 , . . . , u n x ) of the tangent space. The Hodge operator is an isomorphism between k (M) and n x −k (M). These two vector spaces have the same dimension The discrete Hodge star operator has to realise an isomorphism between k (K ) and n x −k (K ). This is not possible on the same grid because generally dim k (K ) = dim n x −k (K ). For example, with the mesh on Fig. 27, n x = 2 and dim 0 (K ) = 4, dim 2 (K ) = 2.
In order to define a discrete Hodge star operator, a "dual mesh" K is needed [97]. This dual mesh is such that there is a one-to-one correspondance between k-simplices of K and (n x −k)-simplices of K . One of the simplest dual meshes used in DEC is the circumcentric dual. In a circumcentric dual of a 2D mesh, the dual of a triangle is its circumcenter; the dual of a primal edge is the orthogonal edge connecting the two triangles which share this primal edge; and the dual of a vertex is the 2-cell formed by connecting the circumcenters of the primal triangles which share this vertex (see Fig. 29 for an example). The orientation on the primal mesh induces an orientation on its dual [220,221]. The boundary operator and discrete exterior derivative can also be transposed to the dual mesh. As pointed out in [112], some care has to be taken when computing the exterior derivative of a form on the dual mesh. Indeed, the dual cells situated at the boundary of the domain are not closed. So, a boundary complement to the discrete exterior derivative has to be added. This complement can be obtained from boundary conditions for many problems.
There are many ways to design a discrete Hodge star operator, each of them leading to a different scheme. One definition, coming from the finite-difference time-domain method [232], is In Eq. (146), c ∈ K n x −k is the dual of the k-simplex c. For a k-dimensional cell σ , |σ | is the Euclidean measure of σ if k ≥ 1 and |σ | = 1 if k = 0. Definition (146) can be seen as a consequence of relation (144) with a low order approximation, knowing that, with a circumcentric dual mesh, a simplex c and its dual c are orthogonal. The matrix of the discrete resulting from (146) is diagonal and positive defined. As said, relation (146) is not the unique way to define the Hodge operator. In particular, it has the drawback that, like in covolume method [117], it does not handle non-acute triangulation. Amelioration and alternatives can be found in literature [116,221,227,[233][234][235][236]. For example, as presented in [221], on can choose a dual mesh based on the barycenter or on the incenter, instead of the circumcenter, to remove the angle condition. In [116], a discrete Hodge star operator is built from relation (148). A discrete Hodge, seen as a mass matrix of a Galerkin method, was also introduced in [103,237]. The construction of efficient discrete Hodge star operators is still an open question.
In the following, as we work exclusively on acute triangulated domains, we use definition (146) of the discrete Hodge operator for our numerical experiments.
The continuous codifferential operator is defined as for any differential k-form ω ∈ k (M). Relatively to the inner product the codifferential can be viewed as the adjoint of the differential operator d when M has an empty boundary, due to Stokes' theorem, since When M has a non empty boundary, Eq. (149) holds provided that ω or θ has a compact support in the interior of M. More generaly, Eq. (149) holds if ω or θ belongs to the Sobolev space H 1 0 ( k (M)) of forms having finite norm [(relatively to the inner product (148)] and vanishing at the boundary of M (see [124,125]).
With the previously built discrete exterior derivative and discrete Hodge star, Eq. (147) enables to define the discrete codifferential operator. In fact, applying δ on the primal mesh boils down to applying d on the dual mesh, and vice-versa, up to a sign and multiplications by measures of simplices.

Interior product and Lie derivative
To deal with some problems such as the Navier-Stokes equations, a discrete interior product or a discrete Lie derivative is needed. One possibility to define a discrete interior product is to use the discrete Hodge star operator and apply the following property of the interior product (see [214]): for any vector field u and any k-form ω. And Lie derivative L can be broken down to wedge products and interior products using Cartan's magic formula. A remark has to be done on formula (150). It necessitates a metric to define the discrete Hodge star and then the interior product. In many applications, a metric is available (for instance, a metric is required to build constitutive laws). So, formula (150) can be used in these applications. However, it has the serious drawback that it defines the discrete interior product and Lie derivative, which are metric independent in the continuous case, from the metric dependant Hodge operator. Alternatives can be found in literature. For instance, a discrete Lie derivative can be built from extrusion, as developed in [238,239].

Applications to fluid mechanics
Consider an incompressible fluid flow, along with a mass transfer, governed by the Navier-Stokes equations In these equations, u and p are respectively the fluid velocity and pressure, ρ and ν are the density and the kinematic viscosity. θ is a passive scalar, and κ is its diffusivity coefficient. Let us reformulate equations (151) in exterior calculus language by applying the flat operator on them.

Numerical scheme
For Eqs. (152a), (152b), we use the same scheme as in [112]. To remove the pressure term, we apply d to Eq. (152a): In time, a fractional step method is used.
Boundary conditions have to be provided for ψ. Whenever needed, the dynamic pressure p = p + 1 2 ρ||u|| 2 can be computed by applying the operator δ on Eq. (152a): Note that applying δ on (152a) is the counterpart, in tensor calculus, of the application of the divergence operator to the momentum equation in (151). So, if div curl were not exactly zero, then the resolution of (157) would introduce an additional spurious pressure. In DEC, since d 2 is exactly zero, there is no such spurious pressure. In space, a discretization with DEC is used. ψ and θ are placed on primal vertices (the velocity form ω is then on dual edges). The pressure is placed on the dual vertices.

Numerical tests
The first numerical example is a 2D channel flow. The length of the channel is L = 2π and its height is H = 1. The viscosity is set to ν = 1 and the diffusivity to κ = 1.
A zero flow is given as initial condition. The value of θ at the top wall is set to 20 and at the bottom wall and the inlet to 0. The grid is composed of 2744 triangles. The mean length of the edges is x mean = 2.989·10 −2 . The time step is t = 10 −2 .
The computed velocity at t = 5s is presented in Fig. 30, left. As can be observed, the computed profile is matches the exact one. The average absolute error is about 2.4·10 −4 . Figure 31 presents the magnitude of θ and some iso-contours. It shows the evolution of θ along the flow. In particular, it can be observed that θ has a linear profile, like the exact solution, in the fully developped region downstream the flow. It is confirmed by Fig. 30, right.
For the numerical analysis of convergence, Couette and Poiseuille flows were simulated, with x mean ·10 2 = 8.408, 4.917, 2.091, 1.393, 1.044, 0.8354 and 0.6961, in a square domain. For primal 0-cochain ω (stream function and temperature), the discrete norm is defined as It corresponds to the usual L 2 norm for a 0-form which is constant inside each dual cell. This definition is used instead of the usual L 2 norm of the reconstructed formω because it is easier to compute. For dual 0-forms (such as the pressure), the following norm is used for the error analysis: For dual 1-forms (like the flux v), we use the edge-based norm: ω, e 2 |e|. Figure 32 shows the evolution of the error with the grid size. As can be observed, the rate of convergence of the stream function ψ is about 1.84. The temperature converges with a second order rate. The flux and the velocity are first order while the pressure convergence rate is 1.65.
For the Poiseuille flow, the convergence rate of the stream function is 1.37, as shown on Fig. 33. For the temperature, it is about 1.95. The flux and the velocity are first order as in [112]. Figure 34 shows the evolution of the discrete norm of δv with time. As can be seen, δv remains below the machine precision (about 2·10 −16 ), as well for the Couette flow as for the Poiseuille flow. This is due to the fact that, with DEC, d 2 = 0 at machine precision.

Relative error
Δx mean As a consequence, the incompressibility condition (152b) is verified exactly, as well localy as globally.
To show the potential of DEC in predicting more realistic problems, we deal, in what follows, with the convection of a polluant in a ventilated room. The room geometry is presented in Fig. 35. The aspect ratio is set to L/H = 1. The height of the inlet and outlet is h = H/5.
A polluted air flow, with a parabolic profile, is injected at the inlet. The Reynolds number based on centerline velocity and h is 100. The polluant is carbon dioxyde, with θ = 1000 at the inlet. At walls, θ = 0 (polluant particles which stick to the wall are extracted).
The grid, composed of 3436 vertices and 6670 acute triangles, is presented in Fig. 36. The return period τ , that is the time needed to cross a distance equal to the perimeter of the room with the inlet centerline velocity, is taken as time unit. The time step is set to 2.10 −3 τ . Some velocity streamlines and the repartition of the passive scalar are plotted in Fig. 37. This figure shows the evolution of the velocity and θ at t = 0.5, 1 and 2τ .
DEC, as presented here, is used only as a spatial discretization method. The canonical metric of the space is used to define the discrete Hodge operator. This can easily be extended to a discretization on a Riemanian surface [213]. Extension to space-time is possible. For example, for the exterior calculus discretization of problems in gauge theory, the Hodge operator has been defined in [244] relatively to the Minkowski metric. An application to general relativity is presented in [245].

Conclusion
With this article, we aimed at raising the awarness among the readers, and particularly among high-performance computing specialists, about the importance of the geometric structure of their equations. This structure traduces fundamental physical properties and its preservation at the discrete level leads to robust numerical schemes. Without being exhaustive, we tried to give an overview of the most popular geometric structures met in mechanics. For each of them, we presented a way to build a corresponding structure-preserving integrator.
Firstly, for Hamiltonian problems, we showed that symplectic integrators are particularly robust for long-time evolution problems. Their error on the energy of the system is bounded over an exponentially large time interval. More generally, symplectic integrators present generally good properties towards the preservation of conservation laws.
In the case of Hamiltonian PDEs, we presented a way of constructing multisymplectic integrators. It was shown that these multisymplectic schemes are more robust regarding the grid, in the sense that they allow more freedom in the choice of time and space steps, compared to classical schemes.
When the equations derive from a variational principle, we showed in section that with total variation approach, discrete energy and momentum evolution equations are obtained naturally, along with the Euler-Lagrange equation. We then presented a way to build variational integrators for ODEs and PDEs. We observed that, when the Lagrangian is time independent, these schemes preserve exactly the energy.
For more general equations, with or without symplectic or variational structure, we presented invariant schemes which preserve the Lie symmetry group of the equation. We saw that, contrarily to classical schemes, they do not introduce unphysical oscillations in presence of a pseudo-shock solution. They also do not give rise to spurious solution when the grid undergoes a Galilean transformation.
At last, we presented DEC which is a space integrator based on exterior calculus and which reproduces exactly Stokes theorem and the relation d 2 = 0 at the discrete scale. The second property enables, for example, to verify exactly the incompressibility condition. This also ensures the exact preservation of circulation in an ideal incompressible fluid flow. The use of DEC in mechanics being very recent, we presented some applications in fluid mechanics with passive scalar convection.
Other geometry-based integrators exist. We can cite, for instance, algorithms built to preserve a Lie-Poisson structure on a discrete manifold [246,247]. Other examples are schemes based on a port-Hamiltonian, or more generally, on a Dirac structure, which aim at correctly handling the interconnection between subsystems [248]. Note that a combination of port Hamiltonian structure and DEC has been considered in [249]. Liegroup integrators can also be cited [250].