Eulerian and Semi-Lagrangian Methods for Convection-Diffusion for Differential Forms

We consider generalized linear transient convection-diffusion problems for differential forms on bounded domains in $\mathbb{R}^{n}$. These involve Lie derivatives with respect to a prescribed smooth vector field. We construct both new Eulerian and semi-Lagrangian approaches to the discretization of the Lie derivatives in the context of a Galerkin approximation based on discrete differential forms. Details of implementation are discussed as well as an application to the discretization of eddy current equations in moving media.


Introduction
Recall the classical linear transient 2nd-order convection-diffusion problem for an unknown scalar function u = u(x, t) on a bounded domain Ω ⊂ R n : in Ω , u = 0 on ∂Ω , u(0) = u 0 . (1) Here, and in the remainder of the paper, β : Ω → R n stands for a smooth vector field. For 0 < ε ≪ 1 we encounter a singularly perturbed boundary value problem, whose stable discretization has attracted immense attention in numerical analysis, see [32] and the many references cited therein.
The boundary value problem (1) turns out to be a member of a larger family of 2nd-order boundary value problems (BVP), which can conveniently be described using the calculus of differential forms. It permits us to state the generalized transient 2nd-order convection-diffusion problems as * ∂ t ω(t) − ε(−1) l d * dω(t) + * L β ω(t) = ϕ in Ω ⊂ R n , ı * ω = 0 on ∂Ω , ω(0) = ω 0 . (2) These are BVPs for an unknown time dependent l-form ω(t), 0 ≤ l ≤ n, on the domain Ω. The symbol * stands for a Hodge operator mapping an l-form to an n − l-form, and d denotes the exterior derivative. Dirichlet boundary conditions are imposed via the trace ı * ω of the l-form ω. We refer to [20,Sect. 2] and [4,Sect. 2] for more details and an introduction to the calculus of differential forms in the spirit of this article. By L β we denote the Lie derivative of ω for the given velocity field β [1,Ch. 5]. It provides the convective part of the differential operator in (2). Details will be given in Sect. 2.1 below.
In two and three dimensions differential forms can be modelled by means of functions and vector fields through so-called vector proxies, see [7,Sect. 7], [4, 1 Occasionally we will use the operator v. p. to indicate that a form is mapped to its corresponding Euclidean vector proxy 1 the Euclidean metric on R 3 , the following vector analytic incarnations of (2) are obtained. The solution l-form ω is replaced with an unknown vector field u or function u, cf. [19,Sect. 2].
• In the case l = 1 we arrive at a convection-diffusion problem ∂ t u + ε curl curl u − β × curl u + grad(u · β) = f in Ω , u × n = 0 on ∂Ω , u(0) = u 0 . (3) • The corresponding boundary value problems for vector proxies of 2-forms read • For l = 3 the diffusion term vanishes and we end up with pure convection problems for a scalar density u ∂ t u + div(βu) = f in Ω , u(·, 0) = u 0 , where we assume β · n = 0 on ∂Ω Though equivalent to (2) the vector proxy formulations very much conceal the common structure inherent in (3)- (4). Thus, in this article, we consistently adopt the perspective of differential forms. Thinking in terms of co-ordinate free differential forms offers considerable benefits as regards the construction of structure preserving spatial discretizations. By now this is widely appreciated for boundary value problems for d * dω, where the so-called discrete exterior calculus [2,13,22], or, equivalently, the mimetic finite difference approach [9,[23][24][25], or discrete Hodgeoperators [8,19] have shed new light on existing discretizations and paved the way for new numerical methods.
All these methods have in common that l-forms are approximated by l-cochains on generalized triangulations of the computational domain Ω. This preserves deRham co-homology and, thus, plenty of algebraic properties enjoyed by the solutions of the continuous BVP carry over to the discrete setting, see, e.g., [20,Sect. 2.2]. Moreover, simplicial and tensor product triangulations allow the extension of co-chains to discrete differential forms, which furnishes structure preserving finite element methods. Also in this paper the focus will be on Galerkin discretization by means of discrete differential forms and how it can be used to deal with Lie derivatives. In light of the success of discrete differential forms, it seems worthwhile exploring their use for the more general equations (2). This was pioneered in [6,29] and this article continues and extends these considerations. Revealing structure is not our only motivation: as already pointed out the discretization of (1) has been the object of intense research, but (3) and (4) have received scant attention, though, (3) is relevant for numerical modelling, e.g., in magnetohydrodynamics. Apart from standard Galerkin finite element methods and [6,29], we would like to mention [12,30] as attempts to devise a meaningful discretization of the transport term in (3). The first paper relies on a Lagrangian approach along with a divergence preserving projection, which can be interpreted as an interpolation onto co-chains. The second article employs edge degrees of freedom (1-co-chains) combined with a special transport scheme. Both refrain from utilizing exterior calculus.
We stress that this article is confined to the derivation and formulation of fully discrete schemes for (2) on a single fixed triangulation of Ω. Principal attention is paid to the treatment of the convective terms. The discussion centers on structural and algorithmic aspects, whereas rigorous numerical analysis of stability and convergence analysis will be addressed in some forthcoming work, because it relies on certain technical results that must be established for each degree l of the form ω separately.
The plan of the paper is as follows: We first review the definition of Lie derivatives and material derivatives for differential forms and give a short summary on discrete forms. Next, we present a semi-Lagrangian approximation procedure for material derivatives of arbitrary discrete l-forms. Afterwards we present an Eulerian treatment of convection of discrete forms. These two sections focus on the formulation of fully discrete schemes. They are subsequently elaborated for the eddy current model in moving media as an example for non-scalar convection-diffusion, see also [17]. Finally we present numerical experiments for convection of 1-forms in 2 dimensions. They give empiric hints on convergence and stability.

Preliminaries
2.1. Lie derivatives and material derivatives of forms. We introduce a spacetime domain where Ω ⊂ R n is a bounded curvilinear polyhedron and [t 0 , t 1 ] ⊂ R. Given a vector field β : Q T → R n and initial values (x, s) ∈Q T the associated flow X(x, s, t) : For fixed (x, s) the solution X(x, s, ·) is called the characteristic curve through (x, s). A unique solution of problem (6) exists, whenever β is continuous inQ T and Lipschitz continuous inΩ for fixed t ∈ [t 0 , t 1 ]. In this case X(X(x, s, t), t, s) = x, ∀x ∈ Ω (6) hence X(·, t, s) is the inverse of X(·, s, t): X(·, t, s) • X(·, s, t) = id .
For convenience we abbreviate X s,t (·) = X(·, s, t). It follows directly from (6) that the Jacobian DX s,t (x) solves: We write F l (Ω) for the space of differential l-forms on Ω (cf. [26, Chapter V, Section 3]). Differential l-forms can be seen as continuous additive mappings from the set S l (Ω) of compact oriented, piecewise smooth, l-dimensional sub-manifolds of Ω into the real numbers 2 . The spaces F l (Ω) can be equipped with the L 2 -norm ω 2 0 := Ω * ω ∧ ω, see [4,Sect. 2.2]. For Euclidean Hodge operator this agrees with the usual L 2 -norm of the vector proxies. Although it is very common to regard the mapping ω : S l (Ω) → C, ω ∈ F l (Ω) as integration we adopt the notation < ω, M l > of a duality pairing instead of M l ω, M l ∈ S l (Ω). Only for l = n we will use the integral symbol. 0-forms can be identified with real valued functions. Recall the definition of the directional derivative for a smooth function f : Ω → R: 2 This notion of differential forms can be rigorously defined with tools from geometric measure theory [16,31]. Here we forgo these technical aspects and appeal to intuition, cf. [20, Sect. The Lie derivative L β of l-forms ω ∈ F l (Ω) generalizes this. Note that the point evaluation of 0-forms is replaced with evaluation on l-dimensional oriented submanifolds M l ∈ S l (Ω). Thus the Lie derivative of a l-form ω is In terms of the pullback X * t,τ defined by we can also write Following [6] the extrusion Ext t,τ (β, M l ) = {X t,s (x) : t ≤ s ≤ τ, x ∈ M l } is the union of flux lines emerging at M l running from t to τ ( Figure 1). The orientation of M l induces an orientation of the extrusion Ext t,τ (β, M l ) such that the boundary is: Plugging this into the definition of the Lie derivative (11) we get by means of Stokes' theorem (13) The contraction operator is defined as the limit of the dual of the extrusion: and we recover from (13) Cartan's formula [27, page 142, prop. 5.3] for the Lie derivative: For 0-forms the second term vanishes, for n-forms the first one.
For time dependent l-forms ω(t) ∈ F l (Ω) the limit value: corresponding to (11) is the rate of change of the action of l-forms in moving media, hence a material derivative [15, page 62]. We deduce: Remark 2.1. Since the exterior derivative and the Lie derivative commute we have: As a consequence closed forms remain closed when they are advected by the material derivative. If D β ω(t) = 0 and dω(0) = 0 then this property is preserved: dω(t) = 0, ∀t.

Whitney forms.
We rely on a version of discrete exterior calculus spawned by using discrete differential forms. We equip Ω with a simplicial triangulation Ω h . The sets of subsimplices of dimension l are denoted by S l h (Ω h ) and have cardinality N l := |S l h |. In words, S 0 h (Ω h ) is the set of vertices of Ω h , S 1 h (Ω h ) the set of edges, etc. We write W l (Ω h ) ⊂ F l (Ω) for the space of Whitney l-forms defined on Ω h , see [7,Sect. 23] and [4,Sect. 3]. From the finite element point of view a Whitney l-form ω h ∈ W l (Ω h ) is a linear combination The construction of these spaces ensures that the traces ı * ω h on element boundaries are continuous. The evaluation maps (l are the corresponding degrees of freedoms, which give rise to the usual nodal interpolation operators Π l : F l (Ω) → W l (Ω h ), also called DeRham maps [20,Sect. 3.1], Remark 2.3. Discrete differential forms are available for more general triangulations [7,Sect. 25]. The Whitney forms feature piecewise linear vector proxies, but can be extended to schemes with higher polynomial degree [20,Sect. 3.4]. For the sake of simplicity we are not going to discuss those.

Direct and adjoint semi-Lagrangian Galerkin formulation
There are two basically different strategies to discretize material derivatives. We either approximate the difference quotient (16) directly or use the formulation (17) in terms of partial time derivative of forms and Lie derivative. The former policy is referred to as Lagrangian and in our setting of a single fixed mesh it is known as semi-Lagrangian. It is explored in this section. The latter is called Eulerian and reduces to a method of lines approach based on a spatial discretization of the Lie derivatives. Its investigation is postponed to Sect. 4.
We begin the study of the semi-Lagrangian method with a look at the pure transport problem: given some ϕ(t) ∈ F l (Ω) and Here and below we assume in the following that on the boundary of Ω the normal components of the vector field β vanish, hence X(Ω, ·, ·) = Ω. This is a restrictive but very common assumption for semi-Lagrangian methods. Only a few semi-Lagrangian methods like ELLAM [11] can handle a non-vanishing velocity field on the boundary.
Under the assumption X t,τ (Ω) = Ω, ∀t, τ : Replacing the limits in (23) and (24) with a finite difference quotient yields semidiscrete timestepping schemes: Given ϕ(t) ∈ F l (Ω) and ω(t − ∆t) ∈ F l (Ω) find and given ϕ(t) ∈ F l (Ω) and Restricting the semi-discrete formulation to the discrete spaces W l (Ω h ), we end up with the following direct and adjoint schemes: Given ϕ(t) ∈ F l (Ω) and A significant advantage of the semi-Lagrangian approach is the preservation of closeness (see Remark 2.1). This property is important in many physical applications. The adjoint semi-Lagrangian timestepping fulfils this in a weak sense.
For ϕ = 0 the adjoint semi-Lagrangian timestepping boils down to As the exterior derivative and pullback commute we conclude is weakly closed. Remark 3.4. Another advantage of semi-Lagrangian methods is the straightforward stability analysis for the homogeneous transport problem D β ω h = 0. The corresponding direct scheme reads: Testing with ω h we immediately get the stability estimate: For Ω ⊂ R 3 and Euclidean Hodge the concrete constants can be deduced from vector proxy representations of the pullback [20, p. 245]. The stability estimates then read: where · is the Euclidean matrix norm. Stability will depend in general on the distortion effected by X t,t−∆t . For l = 1, 2 even the usual assumption divβ = 0 does not guarantee stability.
It is important to note that both the stability estimates and the preservation of weak closedness as previously stated assume exact evaluation of the bilinear forms Any non-trivial problem will require the use of additional approximations. We propose the following three approximation steps.
(i) First we use some numerical integrator to determine approximationsX t,t ′ (s 0 i ) to the exact evolutions X t,t ′ (s 0 i ) of vertices s 0 i . Note that the direct scheme requires to solve (6) backward in time, while for the adjoint scheme we calculate forward in time. The simplest numerical integrator for this first step is the forward Euler scheme, which givesX t,t ′ (s 0 Then we approximate the evolution X t,t ′ (x) of arbitrary points x as the convex combination of the approximated evolutions of surrounding vertices: This yields a piecewise linear, continuous approximationX t,t ′ (x) of the exact evolution X t,t ′ (x).
Note that for small time steps and Lipschitz continuous β the approximating flowX t,t ′ maps the simplicial triangulation Ω h onto a new simplicial triangulationΩ h , such that each subsimplex s l i of Ω h has an imagē s l i =:X t,t ′ (s l i ) inΩ h . (iii) Finally, we apply the nodal interpolation operators of discrete forms [20, Sect. 3.3] to map the transported discrete formsX * t,t−∆t ω h (t − ∆t) and X * t−∆t,tη h onto the space of discrete forms. Let us illustrate these three steps for the case of 0-forms and 1-forms: To determine now the interpolant of a discrete transported 0-formX * t,t ′ ω h , ω h ∈ W 0 (Ω h ), by linearity it is enough to consider the basis functions. Since < Π 0 ω, s 0 i >=< ω, s 0 i > for all vertices s 0 i and 0-forms ω the matrix operator P 0 t,t ′ with entries maps the expansions coefficients of ω h to the coefficients of Π 0X * t,t ′ ω h . This means that in each time step we not only need to determine the pointsX t,t ′ (s 0 i ) but also the location within the triangulation Ω h . To find the element, in whichX t,t ′ (s 0 i ) is located, we trace the path of the trajectory from one element to the next. Based on this data the matrix entries (30) can be assembled element by element (see fig.  2 for the direct scheme).
The advantage of the second approximation step, the linear interpolation, is elucidated by the treatment of 1-forms. The interpolation of a transported discrete 1-formX * t,t ′ ω h is determined by the interpolation of transported basis forms and This defines a matrix P 1 t,t ′ , mapping the expansion coefficients of ω h to those ofX * t,t ′ ω h . The matrix entries are line integrals along the straight line fromX t,τ (s 0 j ) toX t,τ (s 0 k ), if s 1 i is the edge connecting the vertices s 0 j and s 0 k (see fig. 3 for the direct scheme). To determine the entries of the i-th row, we trace the straight line fromX t,t ′ (s 0 j ) toX t,t ′ (s 0 k ) and To determine the location ofX t,t−∆t (s 0 i ), we move backward along the trajectoryX t,· (s 0 i ) starting from s 0 i and identify the crossed elements s 0 1 , s 0 2 , s 0 3 and s 0 4 . In this case p 0 ik , p 0 il and p 0 im are the only non-zero entries in the i-th row of P 0 t,t−∆t .  4 for the direct scheme), then the element contribution to p 1 ii ′ is: where s 0 j ′ and s 0 k ′ are the terminal vertices of edge Observe that we do not need to calculate additional approximate solutions of (6) to determine the matrix entries (31). Piecewise linear interpolation of the approximate evolution of vertices automatically gives approximate evolutions of all subsimplices.
Although the approximate flowX t,t ′ is not smooth, we can show that the fully discrete adjoint scheme preserve the weak closedness property. Is is easily established that exterior derivative and interpolation of approximate pullbacks commute [20,Sect. 3.2]. dΠ lX * t,τ ω = Π l+1X * t,τ dω , ∀ω ∈ W l (Ω h ) . Finally we state the fully discrete Semi-Lagrangian timestepping schemes for the general convection-diffusion problem (2): Definition 3.6. Given a temporal mesh 0 = t 0 < t 1 < t 2 < . . . , the discrete semi-Lagrangian timestepping scheme for the direct variational formulation of convectiondiffusion of differential forms is: Given ω h (t 0 ) ∈ W l (Ω h ) and ϕ ∈ F l (Ω), find ω h (t k ), k > 0 such that: The discrete Semi-Lagrangian timestepping scheme for the adjoint variational formulation of convection-diffusion of differential forms is: Given ω h (t 0 ) ∈ W l (Ω h ) and ϕ ∈ F l (Ω), find ω h (t k ), k > 0 such that: Remark 3.7. The proposed sequence of approximation steps allows for a comprehensive stability and convergence analysis of fully discrete schemes. Stability of the fully discrete schemes in Ω = R 3 for example follows from remark 3.4. We only need to control perturbations due to some numerical integrator, the linear interpolation and interpolation of forms. We will address this in some future work.
Remark 3.8. The standard finite element approach for approximating bilinear forms like Ω * X * t,t ′ ω h ∧η h would be quadrature. But since the product of X * t,t ′ ω h andη h is smooth only on intersections s n i ∩X t ′ ,t (s n j ), this entails finding all intersections of an element s n i of the triangulation with any other element X t ′ ,t (s n j ) of the transported triangulation, and apply some quadrature there. We conclude that such an approach is even more complex than methods based in interpolation. Further, it is totally unclear how such an approximation affects stability and conservation of closedness.
Remark 3.9. If we use in the scalar case l = 0 a vertex based quadrature rule to approximate the element integrals s n i Π 0X * t,τ ω h (τ ) ∧ η, the resulting scheme corresponds to a Lagrange-Galerkin scheme from [28] using vertex based quadrature.
Remark 3.10. Notice that in both (32) and (33) the unknown ω h (t k+1 ) is obtained by solving a standard discrete diffusion problem arising from a symmetric positive definite bilinear form, which are amenable to fast iterative solvers [3,18,21]. In the words of [33], the semi-Lagrangian method is "solver-friendly".

Eulerian approach: discrete Lie-derivatives
The Eulerian approach treats the discretization in space and the discretization in time separately. We therefore consider the stationary version of the convectiondiffusion equation (2): Since the approximation of the diffusion term is well-known [19], we will discuss only the approximation of the convective terms in (34). We will pursue two different ways to obtain spatial discretization for the Lie-derivative terms. 4.1. Standard scheme. The standard Galerkin technique to approximate Ω * L β ω h ∧ η h , ω h , η h ∈ W l (Ω h ) relies on elementwise quadrature. While evaluation of the part Ω * i β dω h ∧ η h is standard in the finite element context, we encounter some difficulties for the term Ω * di β ω h ∧ η h , which occurs for l > 0. The term di β implicitly requires spatial derivatives that fail to be square integrable for discrete differential forms. For instance, for 1-forms ω ∈ F l (Ω) with smooth vector proxy u := v. p.(ω) we have: while for discrete 1-forms ω h ∈ W l (Ω h ) with vector proxy u h on each element the symmetric part (Du h ) sym := 1 2 (Du h + (Du h ) T ) vanishes. A scheme based on elementwise quadrature will always miss (Du) sym β. Stated differently, the exterior derivative di β ω h of the contraction of a discrete form is not defined in a strong sense. But it is defined in a weak sense: For any smooth n − l form γ ∈ F n−l (Ω), ı * γ = 0, and discrete differential form ω h ∈ W l (Ω h ) elementwise integration by parts yields: where d h is the restriction of d to each element s n i . Since each face inside Ω occurs twice in the right sum we can rewrite this term as sum of face integrals over certain jump terms [i β ω h ] l−1 .
One is now tempted to use the right hand side of (38) to define approximations for the convection terms. The difficulty here are the face integrals, which are not welldefined for γ ∈ W l (Ω h ). Motivated by finite volume approximations, we replace in the case γ ∈ W l (Ω h ) in the face integrals γ with a flux functionγ =γ(γ l , γ r ) depending on the values γ l and γ r on adjacent elements.
We use the characterization (11) of the Lie derivative to determine such consistent flux functions.
We could have used one-sided finite differences as well. Next we show that the limit lim ∆τ →0 b ∆τ (ω h ,η h ), ω h ,η h ∈ W l (Ω h ) exists also for discrete differential forms.
Lemma 4.2. The limit b 0 (ω h ,η h ) := lim ∆τ →0 b ∆τ (ω h ,η h ), ω h ,η h ∈ W l (Ω h ) exists, and we have the representation with upwind and downwind traces ω + h and ω − h and the pointwise limitL β ω h (x) = lim ∆τ →0 L ∆τ β ω h (x). J • ⊂ {1, . . . N n−1 } is the index set of faces, that are not on the boundary of Ω. It's cardinality is N • n−1 . Here we used Proof. Let us first note that for all x ∈ ∂s n i the pointwise limitL β ω h (x) := lim ∆τ →0 L ∆τ β ω h (x) exists. We decompose the central difference into a sum of upwind and downwind finite differences: and look first at the limit for the upwind finite difference. We split the integration over the domain Ω in a sum of integrals over patches P i,j (∆τ ) := s n i ∩ X −1 t,t−∆τ (s n j ) with smooth integrand: by the definition of contraction.
After summation over all elements we obtain: for the upwind finite difference. A similar argument for the downwind finite difference shows: Combining these two results for upwind and downwind finite difference we get the assertion.
Remark 4.3. The careful examination of the proof of Lemma 4.2 shows that the existence of the limit can be shown for quite general approximation spaces for differential forms. We used only the continuity inside each element.  Figure 5. Illustration for proof of Lemma 4.2. In the case of the upwind finite difference we have, that discrete l-forms Ω h ∈ W l , l > 0 are discontinuous on edges of the blue triangulation Ω h . Their pullbacks X * t,t−∆τ ω h are discontinuous on edges of the red dashed triangulation X t−∆τ,t (Ω h ), the image of Ω h .
case the face integrals in (40) vanish and we obtain the standard bilinear form for scalar convection: For discrete 1-forms we have both volume and face contributions: and similar for discrete 2 forms: For discrete 3-forms we get: Remark 4.5. For 3-forms in 3D the limit (40) corresponds to the non-stabilized discontinuous Galerkin method for the advection problem from [10]. In the lowest order case this is identical with the scheme for discrete 3-forms. The limit for the upwind finite difference yields the stabilized discontinuous Galerkin methods with upwind numerical flux.
In practice we will use local quadrature rules to evaluate the volume and face integrals in b 0 (ω h ,η h ). This will not destroy the consistency if the quadrature rules are first order consistent.
In the definition of b ∆τ from Def. 4.1 we may also appeal to Lemma 3.1 and shift the difference quotient onto the test function. This carries over to the limit and yields two different Eulerian timestepping methods for problem (2) that, in analogy to Sect. 3, we may call "direct" and "adjoint".

4.2.
Upwind scheme. Now we use the definition (11) of Lie-derivatives as limit of a difference quotient to define a proper interpolation operator for Lie-derivatives of discrete forms. Let s i l ∈ S l h be a subsimplex of the triangulation. Then both limit values, the limit from the upwind direction , ω h ∈ W l (Ω h ) (54) and the limit from the downwind direction , ω h ∈ W 1 (Ω h ) (55) exist for Lipschitz continuous β. The continuity of traces of discrete differential forms ensures that integrals are well-defined. The existence of upwind and downwind limits then follows from the existence of the limits for locally linear β, which in turn can be calculated explicitly by distinguishing different geometric situations. We skip the tedious details.
Observe that in general e.g. for 0-forms a short calculation shows and with upwind and downwind traces ω − h and ω + h . Only for cases where ω h is continuous in a neighbourhood of subsimplex s l i upwind and downwind limit agree. Since (54) and (55) correspond to integral evaluation of the usual interpolation operators, we could either use the upwind or the downwind limit to define discrete Lie-derivatives. We will use the upwind limit for the direct formulation and the downwind limit for the adjoint formulation. This is motivated by taking the cue from semi-Lagrangian schemes. Definition 4.7. Let β be Lipschitz continuous and ω h ∈ W l (Ω h ). Further let (s l i ) N l i=1 be the subsimplices corresponding to the degrees of freedom of basis functions is the upwind interpolated discrete Lie-derivative. Analogously the downwind interpolated Lie-derivative is It follows from the Cartan formula (15) that Remark 4.8. The exterior derivatives can be evaluated exactly for discrete differential forms. Thus, the contraction operators are the only approximations. But for locally constant discrete forms these approximations are exact and one can prove consistency once one has interpolation estimates for discrete differential forms. This is outside the scope of this paper. Remark 4.9. We finally want to stress that it is important to treat the limits (54) and (55) for l > 0 as limits of co-chains rather than integrals with certain integrands that are defined as pointwise limits. First, these pointwise limits are not welldefined since discrete forms are in general discontinuous across element boundaries. Second, even for special cases where the pointwise limit lim ∆τ →0 To understand this, we may consider the 2D example depicted in Fig. 6 with constant β. We take ω h ∈ W 1 such that < ω h , s 1 i >= 0 for all edges s 1 i ∈ S 1 except the edge s 1 1 between vertices s 0 1 and s 0 2 and calculate the projection of the Lie-derivative onto edge s 1 2 between vertices s 0 2 and s 0 3 . Then it is clear that and therefor < lim On the other hand, one can easily show that: Since an explicit calculation of the upwind limits (54) might be expensive for arbitrary β, we have to find consistent approximations. The preceding discussion shows that this must be done very carefully. Replacing the integration simply with some quadrature will not yield a consistent approximation.  Figure 6. Example that we can not change the order of integration and limit.
We combine the upwind discretization with an implicit Euler method and state the implicit upwind methods for convection-diffusion for discrete differential forms.
Definition 4.10. The implicit direct upwind Eulerian timestepping scheme for the variational formulation of convection-diffusion of differential forms is: Given ω h (t 0 ) ∈ W l (Ω h ) and ϕ ∈ F l (Ω), find ω h (t k ), k > 0 such that: The adjoint upwind Eulerian timestepping scheme for the adjoint variational formulation of convection-diffusion of differential forms is: Given ω h (t 0 ) ∈ W l (Ω h ) and ϕ ∈ F l (Ω), find ω h (t k ), k > 0 such that: In contrast to the semi-Lagrangian schemes (32) and (33), in (69) and (70) we face a non-symmetric algebraic system in each time step, whose iterative solution can be challenging. We therefor consider a second kind of upwind Eulerian scheme, that treats the convection in an explicit fashion.
Definition 4.11. The direct semi-implicit upwind Eulerian timestepping scheme for the variational formulation of convection-diffusion of differential forms is: Given ω h (t 0 ) ∈ W l (Ω h ) and ϕ ∈ F l (Ω), find ω h (t k ), k > 0 such that: The adjoint upwind Eulerian timestepping scheme for the variational formulation of convection-diffusion of differential forms is: Given ω h (t 0 ) ∈ W l (Ω h ) and ϕ ∈ F l (Ω), find ω h (t k ), k > 0 such that: These semi-implicit schemes do not coincide with the usual method of lines approach. Nevertheless there are two heuristic arguments that could justify the semi-implicit schemes. Note first that for discrete 0-forms the Semi-Lagrangian and the semi-implicit Eulerian scheme are identical, if we use a time step length ∆t = O(h) and explicit Euler for solving the characteristic equation (6). The direct Semi-Lagrangian scheme (32) for discrete 0-forms was, written here for vector proxies u and v for ω h and η h ∈ W 0 (Ω h ): Find u such that while the semi-implicit Eulerian scheme was: Find u such that On the other handX t k+1 ,t k (s 0 i ) = s 0 i − (t k+1 − t k )β(s 0 i ) and by the time step condition ∆t = O(h): , hence we have the identity: . The second argument, that supports the definition of semi-implicit schemes applies to forms of any degree and follows from the identity: Hence we can derive the semi-implicit schemes (69) and (70) from the Semi-Lagrangian schemes (32) and (33) in replacing e.g.
. Remark 4.12. The analysis of a semi-implicit discontinuous Galerkin scheme for scalar non-linear convection-diffusion can be found in [14]. For the linear case this scheme is modulo certain quadrature rules identical with the semi-implicit Euler scheme for 3-forms.
Remark 4.13. There is also a very close link between both the semi-Lagrangian and semi-implicit Eulerian methods and arbitrary Lagrangian-Eulerian methods (ALE) [5, p. 322]. ALE methods are operator splitting methods, that treat the diffusion part in a Lagrangian and the convection part in an Eulerian fashion. After a certain number of Lagrangian iterations steps an Eulerian iteration step maps the mesh function on the distorted mesh back to the initial mesh. The semi-Lagrangian and semi-implicit Eulerian methods apply this mapping in each time step. In many problems these remap operations should preserve certain properties of the mesh function, e.g. volume or closedness. It is the discrete differential form approach, that naturally guaranties such properties for the methods presented here and for the ALE method in [30].

Application: Magnetic convection-diffusion
In this section we will derive two different convection-diffusion equations for the electromagnetic part of magnetohydrodynamics (MHD) models. The first one is an equation for the magnetic field h and requires the adjoint methods while the second one for the magnetic vector potential a can be solved via the direct methods.
Commonly in MHD applications, one neglects the displacement current. This reduced model, called eddy current model, is a system of equations for the magnetic field h ∈ F 1 (Ω), the electric field e ∈ F 1 (Ω), the magnetic field density b ∈ F 2 (Ω), the current density j ∈ F 2 (Ω) and external source f ∈ F 2 (Ω): in Ω , in Ω , in Ω .
A uniformly positive conductivity σ will be assumed in the sequel.
To rewrite the eddy current model in terms of a material derivative we substitute e =ẽ+i β b and add −i β db to Faraday's law (73). This leaves the solution unchanged, since db = 0. Hence we end up with the system: in Ω .
We eliminateẽ using Ohm's law and end up with the following variational formulation: Seek h ∈ F 1 (Ω), b ∈ F 2 (Ω) such that: In order to eliminate b as well we use Corollary 3.2 and get This variational formulation can be approximated using a Semi-Lagrangian framework (33) yielding an algebraic system: where M µ and C σ −1 are the mass matrices Ω * µ ω h ∧η h and Ω * σ −1 dω h ∧ dη h , ω h ,η h ∈ W 1 (Ω h ). The evaluation of the right hand side requires the calculation of the matrix P 1 τ,t from (31). Lemma (3.5) shows that the solution of (86) is weakly closed if the initial data is weakly closed and f = 0. Alternatively, we may use one of the implicit Eulerian schemes (53) or (70) where L a is the stiffness matrix for either the standard (53) or the downwind (70) discretization. While these schemes will not preserve the weakly closed property, the semi-implicit upwind Eulerian scheme (72) with algebraic system (M µ + ∆tC σ −1 )h t = (M µ + ∆tL a )h t−∆t + ∆tf (88) does, due to the interpolation based construction.

a-based variational formulation.
Here, for simplicity, we assume homogeneous magnetic boundary conditions on ∂Ω. The ansatzẽ = −D β a and b = da solves Faraday's law (77). We then multiply Ampere's law with a test form a ′ ∈ D 1 and integrate by parts: The material laws (79) and (80) eliminate j and h and we get the a-based variational formulation: Find a ∈ F 1 (Ω) such that: The algebraic system for the direct Semi-Lagrangian scheme (32) then reads as while the systems for the implicit and semi-implicit Eulerian schemes are: and Here M σ and C −1 µ are the stiffness matrices for Ω * σ ω h ∧η h and Ω * µ −1 dω h ∧ dη h , ω h ,η h ∈ W 1 (Ω h ). L d corresponds either to the standard (52) or the upwind Eulerian (69) discretization of the convection part.

Numerical experiments
We finally we present a few numerical examples that illustrate the performance of the methods derived here. We take Ω ⊂ R 2 and look at in Ω for 1-forms ω ∈ F 1 (Ω). In vector proxy notion with u := v. p.(ω) this reads with curl := R grad, curl := divR and u × β := u T Rβ, R = 0 −1 1 0 . We approximate ω by discrete 1-forms ω h on a triangular mesh. We study the adjoint Semi-Lagrangian and Eulerian methods (33), (53), (70) and (72) for this boundary value problem. Experiment I. In the first experiment we take in problem (95) the domain Ω = [−1, 1] 2 , the divergence free velocity and u = cos(2πt) sin(πx) sin(πy) Then with ε = 1 neither the convection nor the diffusion part dominates and we expect for all three schemes first order convergence in the mesh size h, if the time step ∆t := t − τ is of order h. The convergence plots in figures 7 and 8 show the L 2 -error at time t = 0.25 and confirm this for different ratios CF L := ∆t β ∞ h . Experiment II. We consider the non-divergence free velocity field β = sin(πx)(1 − y 2 ) sin(πy)(1 − x 2 ) .
Again, we obtain first order convergence (see figs 9) and 10). Experiment III. Next, we examine the stability properties of the three different schemes for dominating convection, meaning that we choose ε in problem (95) very small. While for ε = 1 all three schemes are stable for small and larger CFLnumbers (see fig. (7) and (8)) we encounter this for small ε only for the semi-Lagrangian and the implicit Eulerian scheme (see fig. (11) and (12)).
As expected in the convection dominated case, one has to use either the semi-Lagrangian or fully implicit Eulerian scheme.
Experiment IV. Another question that we address concerns the stabilization for the Eulerian schemes. For the implicit timestepping schemes we need to solve in each time step a stationary convection-diffusion problem. It is well-known that in the scalar case non-stabilized methods would produce highly oscillating solutions for convection-diffusion problems with dominating convection. To study this issue,     We could either use the standard approximation (40) or the upwind interpolated discrete Lie-derivative (59) in the discretization of (97). Again, we choose u = sin(πx) sin(πy) (1 − x 2 )(1 − y 2 ) . and the divergence free velocity (96). Fig. 13 shows that in the balanced case (ε = 1) both the standard and the stabilizing upwind scheme yield meaningful solutions. For small ε the standard method does not seem to converge while the upwind scheme converges (see Fig. 14).