Stabilized Galerkin for Transient Advection of Differential Forms

We deal with the discretization of generalized transient advection problems for diﬀerential forms on bounded spatial domains. We pursue an Eulerian method of lines approach with explicit time-stepping. Concerning spatial discretization we extend the jump stabilized Galerkin discretization proposed in [ H. Heumann and R. Hiptmair , Stabilized Galerkin methods for magnetic advection , Math. Modelling Numer. Analysis, 47 (2013), pp. 1713–1732] to forms of any degree and, in particular, advection velocities that may have discontinuities resolved by the mesh. A rigorous a priori convergence theory is established for Lipschitz continuous velocities, conforming meshes and standard ﬁnite element spaces of discrete diﬀerential forms. However, numerical experiments furnish evidence of the good performance of the new method also in the presence of jumps of the advection velocity.


Introduction
The equations of magneto-hydrodynamics (MHD) [16,Section 3.8], [22,Section 4.1] provide a consistent description of the interaction of electromagnetic fields and conducting non-magnetic fluids like plasmas. The standard model for resistive MHD under a quasi-neutrality assumption comprises balance equations for mass, momentum and energy together with material laws and Maxwell's equations in their magneto-quasistatic reduction (eddy current model) for the electromagnetic fields.
The traditional formulation of the linear eddy current model in the presence of a conducting fluid moving with velocity β = β(x, t) boils down to the evolution PDE governing the evolution of the unknown magnetic vector potential u and with ε being the magnetic diffusion coefficient. Alternatively, one may rely on the magnetic induction field as unknown, again denoted by u, which yields ∂ t u − gradεdivu + αu + βdivu + curl(u × β) = f .
where ω(t) is a time-dependent differential k-form on the bounded domain Ω ⊂ R n , β : Ω × I → R n is a given velocity field and f (t) ∈ Λ n−k (Ω) a source term. The scalar diffusivity parameter ε and the reaction coefficient α are non-negative and bounded functions Ω → R, and the boundary conditions are imposed at the inflow boundary Γ in := {x ∈ ∂Ω : β · n(x) < 0} and at the "elliptic boundary" (where the diffusion parameter ε > 0). All other notations are borrowed from [24, Section 2] and summarized in Table 1.1.
The so-called vector proxies 1 establish the connection between (1.1), (1.2) and (1.3). Indeed, in R 3 endowed with the Euclidean inner product natural isomorphisms between Λ k (R 3 ) and R or R 3 can be defined. The fields associated to differential forms are called proxy fields for the forms and exterior calculus operations on forms correspond to operations on scalar functions and vector fields, see [24,Section 2], [5,  For k = 0 the evolution operator in (1.3) written in vector proxies becomes the familiar and widely studied second order advection-diffusion equation for the unknown scalar function u ∂ t u − divεgradu + αu + β · gradu = f. (1.4) ω ∈ Λ k (Ω) H 1 (Ω) H(curl, Ω) H(div, Ω) L 2 (Ω) Exterior calculus Proxy calculus By analogy we conclude that in the generalized advection-diffusion problem (1.3), the diffusion operator is d⋆d, the zero-th order term amounts to a reaction term and the advection operator is the Lie derivative L β associated to the velocity field β.
It is well known that for scalar advection-diffusion equation (1.4) straightforward Galerkin discretization with Lagrangian finite elements will break down in the singular perturbation limit of vanishing diffusion. Thus, robustness for ε ց 0 will also be a key issue for the spatial discretization of (1.1) and (1.2). In this article we tackle the challenge of robust Eulerian spatial finite element discretization for the general advection-diffusion problem (1.3). In fact, we will focus on the pure advection problem obtained from (1.3) for ε = 0; if a scheme performs well in this case, it will also be suitable for (1.3) when augmented with a standard Galerkin discretization of the diffusion term.
These results require β to be Lipschitz continuous in space. However, MHD solutions feature shocks that give rise to discontinuous velocities; discontinuous transport velocities are relevant in the context of magneto-quasistatic Maxwell's equations, also in the limit of small diffusion.
A well-posedness theory for velocity fields with less regularity is available only for scalar advection. In [18] DiPerna and Lions showed well-posedness of the scalar advection problem for velocity fields β ∈ L 1 loc (0, T ; W 1,1 (R n )) with divβ ∈ L 1 (0, T ; L ∞ (R n )) through the concept of renormalized solutions. More recently, Ambrosio in [1] provided an extension of this breakthrough to transport velocity fields in L 1 loc (0, T ; BV loc (R n )) and divβ ∈ L 1 (0, T ; L 1 loc (R n )). Moreover, a notion of generalized flow associated with low regular velocity fields (the regular Lagrangian flow) and an extension of the characteristics theory to beyond the smooth context have been subject of investigation of several authors, see [15], [2], [8] and the references therein. To the best of our knowledge, beside the case of scalar transport, a well-posedness theory for the generalized transport problem (1.5) with low regular advection velocities has not been developed.
Even though the above mentioned results have been established for nearly incompressible velocity fields (see [17] for a detailed overview), the assumption on the boundedness of the divergence of the velocity (absolute continuity with respect to the Lebesgue measure in the BV case) is of crucial importance for the well-posedness of the scalar advection problem. In the context of the generalized transport problem for a differential k-form, this corresponds to require the operator L β + L β to be bounded in space, which conceals a rather strong assumption on the regularity of the velocity itself, when k = 1, 2.

Novelty and outline
A full discretization of (1.5) was already presented in [24]. There, the authors introduced a semi-Lagrangian approach. Conversely, in the present paper we pursue a mesh based Eulerian method of lines approach to (1.5), employing a (jump) stabilized Galerkin discretization and piecewise polynomial discrete differential forms for spatial discretization. Our new methods will be constructed to accommodate discontinuous velocities aligned with the mesh.
A jump-stabilized discontinuous Galerkin method for the stationary advection problem for 0-forms in R 3 and Lipschitz continuous velocities β ∈ W 1,∞ (Ω), was introduced and theoretically analyzed in [12]. An extension of these results to the magnetic advection problem (1-forms in R 3 , cf. (1.1)) was proposed in [25], where a priori convergence rates were derived for both fully discontinuous piecewise polynomial functions and H(curl, Ω)-conforming finite elements. Discontinuous velocity fields were not taken into account. We remark that for discontinuous velocities, even the spatial discretization of the scalar transport problem (1.4), for which existence and uniqueness of weak solutions are known, is discussed only rarely ( [9], [38]).
The remainder of the paper is organized as follows. In Section 2, we devise a stabilized Galerkin spatial semi-discretization for the generalized stationary advection problem (1.5) for merely piecewise smooth velocity β. It is an extension of the method introduced in [12] for 0-forms and in [25] for 1-forms. Trial and test functions are polynomial discrete differential forms, which will be introduced in Section 2.5. The new method is a substantial extension of the scheme presented in [25] to forms of arbitrary degree, any spatial dimension and velocities with jumps.
Next, Section 3 establishes stability a priori convergence estimates for the stabilized Galerkin discretization in the stationary setting. For want of well-posedness results for the generalized advection problem in case of discontinuous β, these investigations are confined to Lipschitz continuous velocities β ∈ W 1,∞ (Ω). The stability and consistency results obtained in that Section are instrumental for the convergence analysis of the fully discrete scheme in Section 4. We study explicit time-stepping following the approach of [31] and [13].
Finally, in Section 5 and Section 6 the performance of the new method is tested in various numerical experiments for both the stationary and transient generalized advection problem (1.5) in 2D. The tests cover both continuous and discontinuous velocities and employ tensor product grids and triangular meshes.

Stationary generalized advection problem
The Eulerian method of lines policy applies time-stepping after discretization in space. Therefore, we will first address the spatial discretization of (1.5) and we start from the stationary generalized advection boundary value problem for a k-form ω on the bounded computational domain Ω ⊂ R n :

Transmission conditions
We aim for stabilized Galerkin methods that, crudely speaking, involve a penalization of suitable jumps across interfaces inside Ω. In order to select the right jump terms, we have to understand the natural transmission conditions across an internal interface f ⊂ Ω satisfied by a solution ω of (2.1). For smooth velocity β ∈ W 1,∞ (Ω), the requirement L β ω ∈ L 2 Λ k (Ω) read in distributional sense, involves the transmission condition for any oriented (piecewise) smooth n−1-dimensional surface f ⊂ Ω. Here, we wrote [·] f for the jump of a function across the surface f . This formula is a consequence of the integration by parts formula for the Lie derivative 3) The transmission conditions (2.2) carry over to Lipschitz continuous velocity β ∈ W 1,∞ (Ω). Clearly, no transmission conditions are imposed across surfaces tangential to β (characteristic surfaces).

Stabilized Galerkin variational formulation
In the following, let T h = {T } be a cellular partition (generalized triangulation) of Ω ⊂ R n into (curved) polyhedra T . Denote by F • and F ∂ the set of interior and boundary n−1-faces of T h (named facets) and F = F • ∪ F ∂ . The set of facets at the inflow boundary is defined as is the set of facets at the outflow boundary. An oriented facet f has a distinguished normal n f . Any facet f , as part of the boundary of some element T ∈ T h , has either n f = n T | f or n f = −n T | f . Then, given ω ∈ Λ k (Ω), its two different restrictions to f are denoted by ω + and ω − , e.g. ω + := ω | T + where element T + has outward normal n f . Hence, we can introduce the notion of jump and average across a facet f ∈ F • as For f ⊂ ∂Ω we assume f to be oriented such that n f points outwards and [ω] f = {ω} f := ω. We also write h T := diam T and h := max T ∈T h h T . Further, let Λ k h (T h ) denote some piecewise polynomial approximation space for differential k-forms. Here Λ k h (T h ) could be either a HΛ k (Ω)-conforming space Λ k h (T h ) ⊂ HΛ k (Ω) or a non-conforming space Λ k h (T h ) ⊂ L 2 Λ k (Ω) for which Λ k h (T h ) ⊂ HΛ k (Ω) . The method is formulated in the general framework of time-dependent velocity fields β = β(x, t)and relies on the assumption that the possible (space) discontinuities of the velocity are resolved by the mesh: Assumption 2.1. For every t ∈ I we have β(·, t) | T ∈ W 1,∞ (T ) for each T ∈ T h , that is the velocity field is assumed to be T h -piecewise Lipschitz continuous.
This may seem to be a severe limitation but for our purposes it represents a reasonable condition in view of the fact that the velocity field is obtained from numerically solving the MHD system.
Next, multiplying equation (2.1a) by a test form η h ∈ Λ k h (T h ) and applying the integration by parts rule (2.3), results in Let j β be the formal adjoint of the contraction operator i β as in Table 1.1. Applying the following product rule to the boundary terms, results in, for all Moreover, it can be easily verified that, for all For ω ∈ W solution of problem (2.1), the transmission conditions (2.4) at the mesh facets As it is well-known, classical Galerkin finite element discretization of advection problems suffer from instabilities. Therefore, devising stabilization techniques to counteract this limitation has been investigated widely. We consider the following stabilization operator (2.8) where the stabilization scalar functions c f (x) andc f (x) may depend on the velocity field and on the facets diameter h f . Throughout, the stabilization parameters are assumed to satisfy the following: In particular, by considering the direction of the numerical fluxes as given by the average of the velocity field, the choice gives a scheme with upwind fluxes (see [23,Remark 4.1.2] in the case β ∈ W 1,∞ (Ω)). Indeed, from (2.7) together with (2.8) the facets contribution, for Note that, since the velocity field is discontinuous, the upwind direction at the mesh facets may not be well defined. Here we consider the direction of the stream as the one given by the average of the velocity. However, other possibilities are feasible: an upwind direction given locally by the velocity field can be used, even if this choice will lead to non-unique numerical fluxes at mesh facets. The evaluation of the terms in (2.7) involving the Lie derivative L β η h requires the knowledge of the first order derivatives of the velocity field β. Note that since the velocity is assumed to be a smooth function in all elements T ∈ T h , the quantity (ω h , L β η h ) T is well defined for all T ∈ T h . However, as suggested in [25], a different equivalent formulation of the bilinear form a h (·, ·) is convenient for implementation purposes.

Proposition 2.3. The following equality holds for all
(2.10) Proof. By using a Leibniz rule for the exterior derivative with respect to the wedge product and Stokes' theorem [37, Theorem 1.2.7], it easily follows that Hence, using (2.11) together with Cartan's homotopy formula for the adjoint of the Lie derivative L β results in Furthermore, exploiting (2.5), yields where the outflow boundary terms can be recast as Therefore, substituting (2.12) and (2.13) into the bilinear form (2.7) yields the conclusion.
Let us consider the particular case of velocity fields that feature Lipschitz continuity in space, that is β ∈ W 1,∞ (Ω). An easy computation allows to write, for all and similarly for the stabilization terms in (2.8). Since trivially [β] f ≡ 0 for all f ∈ F • , all the terms involving the jump of the velocity can be dropped and the variational problem reduces to: If c f =c f , using (2.5) the bilinear form above can be recast as (2.14) and the formulation in [23,Equation 4.8, p. 61] is recovered. Note that, owing to the fact that {β} f = β| f , the choice of stabilization given in (2.9) yields a scheme with genuine upwind fluxes. Moreover, since the stabilization terms vanish for ω ∈ W solution of (2.1), the variational formulation with stabilized bilinear form given by (2.7) and (2.8) is consistent with (2.1), namely Galerkin orthogonality Observe that the stabilized Galerkin formulation (2.6), (2.14) for Lipschitz continuous velocities β ∈ W 1,∞ (Ω) can be equivalently derived by imposing on the mesh facets the transmission conditions (2.2).
Analogously, the bilinear form corresponding to the reformulated variational problem (2.16)

Stabilized Galerkin formulation in terms of vector proxies
For the sake of completeness, we present the vector proxy representation of the stabilized reformulated bilinear form (2.10), (2.8) corresponding to the variational formulation associated with the transport problem of the corresponding k-form.

Trial and test spaces of discrete differential forms
From now, we restrict ourselves to special types of meshes: • T h is either a simplicial decomposition of Ω ⊂ R n as defined in [5, Section 4.1 and Section 5.3], • or a tensor product mesh, namely a compatible, locally quasi-uniform, affine mesh partition of Ω into non-degenerate axiparallel parallelotopes.
On such meshes various families of piecewise polynomial discrete differential forms of any degree have been constructed, see [5], [3], [26, Section 3] and [4] for a detailed overview. For a simplicial decomposition T h , the spaces of polynomial totally discontinuous discrete differential k-forms on T h are defined as is the space of differential k-forms with polynomial coefficients of degree at most r on the n-cell T ∈ T h , obtained as the restriction of P r Λ k (R n ) to T . The corresponding space of HΛ k (Ω)-conforming discrete differential forms , T ∈ T h } allows to introduce another family of polynomial differential k-forms on T h , namely is the space of homogeneous polynomial differential k-forms of degree r and κ : denotes the Koszul differential [5, Section 3.2]. The so-called "first family" of finite element differential k-forms is hence defined as P − r Λ k (T h ) = {ω ∈ HΛ k (Ω) : ω | T ∈ P − r Λ k (T ), T ∈ T h } whose functions satisfy the continuity requirement that the trace tr ω is single-valued on all n−1-cells which in turn ensures inclusion in HΛ k (Ω).
The family Q − r Λ k (T h ) of finite element differential forms on a tensor product mesh can be constructed by iteratively applying a tensor product strategy from the 1-dimensional interval.
This tensor product construction allows to build a subcomplex of the de Rham complex. We refer the interested reader to [4,Section 5] for details on such constructions.
The HΛ k (Ω)-conforming finite element spaces presented above over a cell complex T h form a discrete de Rham sequence as a cochain projection from the de Rham complex through projection operators Π k h , namely the following diagram commutes [5,Section 5.5] and [3]).
In the subsequent analysis we will make use of pairs of HΛ k (Ω)-conforming spaces Λ k h (T h ) and non-conforming spaces Λ d,k h (T h ) as in the following: Remark 2.5. We pay particular attention to HΛ k (Ω)-conforming trial/test spaces because they allow the straightforward Galerkin discretization of the diffusion form d ⋆ d present in (1.3).

Stationary transport: Estimates for continuous velocity fields
As explained in Section 1.2, a rigorous convergence analysis is only possible in the case β ∈ W 1,∞ (Ω), for want of a well-posedness result for (2.1) with discontinuous velocity fields. Hence, all theoretical results will rely on the assumption β ∈ W 1,∞ (Ω). Moreover, without loss of generality we can assume the scaling β L ∞ (Ω) = 1.
Also, in this section, we admit only the same special classes of meshes as in Section 2.5. We recall the notion of shape regularity of a mesh as presented, e.g., in [26, Section 3.6] following [14, Section 3.1]. Throughout, let Λ k h (T h ) denote some piecewise polynomial approximation space for differential k-forms. If not otherwise specified, Λ k h (T h ) could be either a HΛ k (Ω)conforming approximation space P r Λ k (T h ) or P − r+1 Λ k (T h ), but also the totally discontinuous space P d r Λ k (T h ) on a simplicial mesh and Q − r+1 Λ k (T h ) on a tensor product mesh. Note that we will always assume that Λ n with a h (·, ·) and s h (·, ·) as in (2.14). Note that the bilinear form s h (·, ·) associated to the stabilization operator is symmetric and nonnegative on where Λ := −(2α id +L β + L β ). Let us introduce the following norms on V (h), and Note that the above norms are well-defined in view of the definition of inflow and outflow boundary and the fact that tr i β (ω ∧ ⋆ω) = (β · n f ) tr i n f (ω ∧ ⋆ω) together with Assumption 2.2 on the stabilization coefficients c f , f ∈ F • . The above norms are combined into where |·| H 1 Λ k (T h ) stands for the broken HΛ k -seminorm on T h . Convergence estimates for the spatial discretization of the stationary boundary value problem are key to analyzing the convergence of the fully discrete scheme. They hinge on stability results for the differential operator L h := A h + S h . In order for these results to hold for both non-conforming and HΛ k (Ω)-conforming space discretization, we approximate discontinuous differential forms by differential forms in HΛ k (T h ) as in the following: with constant C > 0 depending only on the polynomial degree and the shape regularity of the mesh.
The proof of this result is constructive and relegated to Appendix A. Note that the construction of the HΛ k (Ω)-conforming approximation is based on an averaging interpolation which is the extension to discrete differential k-forms of the operator introduced and studied for scalar functions in R 3 in [28] and [27,Appendix].
Lemma 3.2. There exists a constant C S depending on the stabilization coefficients |c f | 1/2 , the polynomial degree r and the shape regularity of the mesh, such that where the constant C L depends on α and |β| W 1,∞ (Ω) , and the constant C ′ L depends on the stabilization coefficients |c f | 1/2 , |c f | −1/2 , the polynomial degree r and the shape regularity of the mesh.
In order to show (3.5), let η h ∈ Λ k h (T h ). Integration by parts yields Inserting into the expression of a h (ω, η h ) in (2.16) (given in Proposition 2.3) results in (3.6) Therefore, using Cauchy-Schwarz inequality gives The last term can be bounded using inverse trace inequalities as follows (3.8) Combining (3.7) with (3.8) and using inverse trace inequalities on the boundary terms yields with C ′ L > 0 depending on |c f | 1/2 and |c f | −1/2 , the polynomial degree and the shape regularity of the mesh. Note that the above result applies to both HΛ k (Ω)-conforming and non-conforming space discretization.
For the sake of conciseness, the following bounds on orthogonal subscales are presented in the case of full polynomial spaces of discrete differential k-forms on simplices, case (I), namely for the spaces Λ d,k h (T h ) = P d r Λ k (T h ) and Λ k h (T h ) = P r Λ k (T h ). The proof for the cases (II) and on tensor product meshes (III), follows mutatis mutandis. Lemma 3.3. The following statements hold true: where the constants C π and C ′ π depend on α, |β| W 1,∞ (Ω) , the stabilization coefficients |c f | 1/2 , |c f | −1/2 , the polynomial degree r and the shape regularity of the mesh.
Proof. In order to show (i), we proceed as in [25,Theorem 3.1], namely one can add to the bilinear form (L h (ω − Π h ω), η h ) Ω given in (2.14) the zero term T ∈T h ω − Π h ω, L β h η h T where β h is the L 2 -projection of β ∈ W 1,∞ (Ω) onto piecewise constant vector fields. Hence, using estimates on the projection error for β in the L ∞ -norm, Cauchy-Schwarz inequality and inverse inequalities, results in where the interior facet terms have been bounded has in (3.8).
In the case (ii) of HΛ k (Ω)-conforming discretization, we can proceed as above, use estimate (3.9), but the non-zero term T ∈T h ω − Π h ω, −L β h η h T has to be bounded. We show that for all ω ∈ V (h), η h ∈ P r Λ k (T h ), with the constant C > 0 depending only on the polynomial degree and the shape regularity of the mesh. In order to do that, we exploit the fact that since β h is piecewise constant, for all T ∈ T h , L β h = −L β h and we build HΛ k (Ω)-conforming approximations for each of the two terms appearing in Cartan's formula L β h = i β h d + di β h . In particular, in view of Proposition 3.1, let γ c,k ∈ P r Λ k (T h ) be the HΛ k (Ω)-conforming approximation of i β h dη h ∈ P d r Λ k (T h ) and let γ c,k−1 ∈ P r+1 Λ k−1 (T h ) be the HΛ k (Ω)-conforming approximation of i β h η h ∈ P d r+1 Λ k−1 (T h ). Since γ c,k , dγ c,k−1 ∈ P r Λ k (T h ),

and by Cauchy-Schwarz inequality and inverse inequality, one has
.

By the approximation results in Proposition 3.1, the projection errors can be bounded as
Upper bounds for the facet terms can be derived as follows: Let us decompose the velocity field β into its normal component β n := (β · n)n and its tangential component β t := (n × β) × n. Then we can write If f = ∂T + ∩ ∂T − , using estimates on the projection error for β, trace and inverse inequalities together with (3.11) and the fact that dη h ∈ HΛ k+1 (Ω) due to (2.17), results in which leads to the desired estimate (3.10). Finally, combining the estimates (3.9) and (3.10) yields

Note that (3.5) together with the definition of · * in (3.3), inverse and trace inequalities gives
and C depends only on the polynomial degree and the shape regularity of the mesh. As it will be shown in the following, stability of second order Runge-Kutta schemes can be achieved with the standard CFL condition if the space discretization is performed with piecewise linear finite elements. Therefore, this case is tackled separately. In particular, we can establish the following estimate.
there exists a constant C π which depends on α, |β| W 1,∞ (Ω) , the stabilization coefficients |c f | 1/2 , |c f | −1/2 and the shape regularity of the mesh, such that for all Proof. The proof we propose is very similar to the proof in [13, Lemma 2.1]. Let us first consider the case of non-conforming truly discontinuous discretization, namely Λ k h (T h ) = P d 1 Λ k (T h ). Using the formulation in (3.6) and proceeding similarly as in the proof of (i) from

Lemma 3.3, one has
where we have used the fact that, for where we have used the estimate derived in (3.13) and again the fact that L β h ω h , y h T = 0 for all T ∈ T h . Moreover, using the bound (3.10), there exists C > 0 depending on the shape regularity of the mesh such that which concludes the proof.
As a consequence of the estimates shown in Lemma 3.2, Lemma 3.3 and Lemma 3.4, we present a convergence result for the stationary advection problem with Lipschitz continuous velocity fields. Analogous estimates were proposed in [23, Theorem 4.1.8] for non-conforming differential forms in R n . The present result extends to HΛ k (Ω)-conforming discrete differential forms in R n the a priori convergence estimates derived in [12, Section 5] for 0-forms in R 3 and in [23, Theorem 4.1.13 and 4.1.14] for 1-and 2-forms in R 3 . Note that the numerical experiments presented in Section 6 for non-Lipschitz velocities indicate that the following convergence result might hold in a more general setting. Theorem 3.5. Let α ∈ L ∞ (Ω) and β ∈ W 1,∞ (Ω) in (1.5) satisfy the monotonicity condition (1.6). Furthermore, let the stabilization parameters c f fulfill the non-negativity Assumption 2.2. Then Moreover, if ω ∈ H r+1 Λ k (Ω) is solution of the advection problem (1.5) and ω h ∈ Λ k h (T h ) is solution of the discrete variational formulation with bilinear form given in (2.14), then with the constant C > 0 depending on |c f |, |c f | −1 , α, β, the polynomial degree r and the shape regularity of the mesh.
Proof. The proof of stability (3.14) immediately follows by (3.1), the positivity condition (1.6) and the definition of the h-norm.
Let Π h denote the L 2 -projection into Λ k h (T h ). By stability and consistency (2.15), one has We proceed as in the proof of Lemma 3.3 (i) to get (3.9) and use a multiplicative trace inequality (see [10, Theorem 1.6.6]) for the interior facets terms, i.e.
with C depending only on the shape of T . Moreover, in the case of HΛ k (Ω)-conforming discrete differential forms the extra non-zero terms are bounded as in (3.10). The approximation estimates [5,Theorem 5.3] inf for C > 0 independent of h, yield the conclusion.
Note that in the non-stabilized case (c f = 0), the bilinear form in (2.14) is coercive in the L 2 Λ k -norm but only sub-optimal convergence is attained, namely ω − ω h L 2 Λ k (Ω) ≤ Ch r ω H r+1 Λ k (Ω) holds with C > 0 independent of the mesh width h.

Fully discrete problem
In the present section, we formulate the fully discrete advection problem for a differential k-form by coupling the stabilized Galerkin spatial discretization introduced in Section 2 with explicit time-stepping schemes. In particular, the forward Euler method and explicit secondorder and third-order Runge-Kutta (RK) schemes are investigated.
On the time interval I = [0, T ], we consider a uniform partition N −1 n=0 [t n , t n+1 ] for a given positive integer N and t n = nτ with uniform time step τ such that T = N τ . The semi-discrete problem reads: find ω h (t) ∈ Λ k h (T h ) such that where the bilinear forms a h (·, ·), s h (·, ·) and the linear functional l(·) are obtained at each time step through spatial discretization as in (2.7), (2.8) and (2.6) with forcing term f (x, t) and velocity field β(x, t). The semi-discrete problem (4.1) can equivalently be recast as the finite dimensional operator evolution equation . In light of the results established in Lemma 3.2, Lemma 3.3 and Lemma 3.4, quasi-optimal convergence rates for the L 2 Λ k -error in space L ∞ -error in time can be proven for smooth solutions of the problem (4.1), along the lines of the analysis proposed by Burman et al. in [13]. In particular, under CFL-type conditions, the efficacy of the proposed space-time discretization lies in the fact that the anti-diffusive nature of explicit RK schemes is compensated by the artificial dissipation introduced through the stabilized spatial discretization.
In the following paragraphs, we introduce the fully discrete problem for (4.2) by explicitly stating the stages corresponding to the time-stepping. Moreover, we present the convergence results corresponding to each fully discrete scheme.

Explicit Euler Scheme
The first order explicit Euler scheme for the problem (4.2) reads where ω n h = ω h (t n ), L n h := L h (t n ) and F n h := F h (t n ).
Then there exist constants C, C CFL > 0 depending only on the constants in Lemma 3.

2, Lemma 3.3 and Lemma 3.4 and the trial/test spaces
Proof. For the sake of clarity we briefly sketch here the underlying idea proposed in [31] and [13], valid also for the higher order time schemes introduced below. For ω n = ω(t n ) exact solution at the n-th time step, the proof starts with writing the error generated at each stage (here one stage) of the scheme as and bounding the error e n h by the approximation error e n π . In order to do that, by deriving the equations governing the time evolution of the error e n h , an energy identity associated with the time scheme can be identified. Starting from this identity, the desired estimate is obtained by deriving upper bounds on each terms via the estimates established in Lemma 3.2, Lemma 3.3 and Lemma 3.4. Note that (3.12) is crucial in order to achieve stability under the CFL-condition in (ii), while the bound on the orthogonal subscales inferred in Lemma 3.3 and Lemma 3.4 is instrumental in the derivation of quasi-optimal error estimates. The conclusion follows by a Gronwall type argument and standard estimates on the projection error e n π .
Remark 4.2. Note that the mild time step constraint in Theorem 4.1 (i) valid for spatial approximations based on piecewise constants discontinuous elements (finite volume) Λ k h (T h ) = P d 0 Λ k (T h ) hinges on the trivial observation that the bound (3.5) in Lemma 3.2 reduces to . This is the standard CFLcondition for upwind finite volume or finite difference schemes for scalar advection problems.

Explicit RK2 Schemes
We consider, as in [13, Section 3.1], explicit Runge-Kutta scheme of order two (RK2) for the problem (4.2) of the form Similarly to the explicit Euler scheme, convergence of the fully discrete problem with second order two-stage Runge-Kutta schemes of the form (4.4), (4.5) can be established. A detailed proof of the following theorem can be found in [13, Theorem 3.1 and Theorem 3.2].
Then there exist constants C, C CFL > 0 depending only on the constants in Lemma 3.2, Lemma 3. 3

Numerical experiments in 2D: continuous velocity
The two dimensional transport problem written for a 1-form amounts of solving (as in [25,Section 5]) where R corresponds to a clockwise rotation of π/2. We consider the pure transport (α       Remark 5.1. Though we observed an increased convergence rate of the spatial discretization (also in the stationary case) compared to the predictions of the theory, our results are probably sharp. In the case of scalar transport it is well-known [35] that on very special meshes, sometimes called Peterson meshes, the L 2 -norm estimates are sharp.
6 Numerical experiments in 2D: discontinuous velocity Lacking a sound convergence theory for the generalized advection problem in the presence of discontinuous velocity fields, the first set of experiments aims at testing the numerical performances of the stabilized Galerkin spatial discretization (proposed in Section 2) for the stationary advection problem.

Stationary problem: Test of convergence
Let us consider the stationary pure transport (α = 0) problem corresponding to (5.1) in the unit square Ω = [0, 1] 2 . We perform a set of numerical simulations on unstructured meshes {T h } h obtained by uniform refinement of an initial partition T 0 which resolves the jump discontinuity in the velocity field β = (β 1 , β 2 ) ⊤ . The velocity is assumed to be piecewise polynomial with respect to the open subdomain partition Ω 1 = (0, 0.5) × (0, 1) and Ω 2 = (0.5, 1) × (0, 1). Namely, The data f and g in (5.1) are chosen such that the strong solution of the BVP is given by the discontinuous vector field u with components Note that the exact solution is tangentially continuous and such that the forcing term f belongs to L 2 (Ω). We perform a numerical discretization based on:  Figure 6.3 show the behavior of the L 2 -error as the mesh is refined for the non-stabilized and stabilized Galerkin spatial scheme introduced in Section 2. As with Lipschitz continuous velocity fields, convergence rate r+1 of the L 2 -error is attained by the stabilized scheme with edge elements of polynomial degree r. For lowest order conforming elements, the rate deteriorates by a factor of 1 when the non-stabilized scheme is applied, whereas higher order polynomial discretization yields numerical solutions which suffer of large oscillations.
Note that a discretization based on the full polynomial space (case (ii), Figure 6.2) the error behaves as in the case of Lipschitz continuous velocity fields when the polynomial degree r is odd. Even polynomial degrees lead to a deteriorated convergence rate of the error in the L 2 -norm.

Stationary problem: Velocity with non-resolved discontinuities
The derivation of the method, as of Section 2, relies on the assumption that the mesh resolves the possible discontinuities of the velocity field. In the following experiment we investigate an example of normal jump discontinuity in the velocity field not resolved by the mesh or any of its refinements and observe that the instabilities arising downstream of the discontinuity irremediably compromise the accuracy of the numerical solution and wreck the performance of the method. The failure of the numerical scheme in this test case may be ascribed to the fact that, since across the mesh facets the jump of the velocity vanishes, the scheme itself does not capture the discontinuity of the velocity and hence of the solution. Jump discontinuities are only taken into account through numerical quadrature. In more details, the pure magnetic transport problem is solved in the domain Ω = [0, 1] 2 with a tensor product mesh and velocity field β = (β 1 , β 2 ) defined component-wise as The data f and g are such that the strong solution of the stationary problem corresponding to (5.1) is given by u = (u 1 , u 2 ) with as shown in Figure 6.4 (first column). The numerical discretization is performed with lowest order edge elements i.e. Λ 1 h (T h ) = Q − 1 Λ 1 (T h ) and upwind stabilization. On the basis of Figure 6.4, it can be inferred that the numerical solution obtained with the upwind stabilized scheme is not affected by spurious oscillations but fails to reproduce the exact solution. An error analysis provides evidence of large errors along the discontinuity and no convergence is achieved.

Transient problem: Test cases, shear and collisional velocities
Encouraged by the promising performances of the scheme presented in Section 6.1 for the stationary advection problem in the presence of resolved discontinuous velocities, we tackle the full discretization of the transient problem.
Let us consider the two dimensional transient magnetic advection problem (5.1) in the space domain Ω = [0, 2] 2 with periodic boundary conditions at the boundary ∂Ω and time domain I = [0, 2]. The numerical discretization is performed on a tensor product mesh with lowest order rotated Raviart-Thomas elements Λ 1 h (T h ) = Q − 1 Λ 1 (T h ) and the second order Heun method as time integrator (in order to avoid restrictive time steps). Let 1 S be the characteristic function on the subset S ⊂ Ω, we consider different velocity fields whose discontinuities are resolved by the domain partition Ω = Ω 1 ∪ Ω 2 with Ω 1 = (0, 1) × (0, 2) and Ω 2 = (1, 2) × (0, 2). Numerical simulations have been conducted with the following velocity fields and initial conditions: Figure 6.5 (a)) and u 0 (x, y) = (x(2 − x)y(2 − y), sin(2πx)) ⊤ as initial condition: The initial condition is in H(curl, Ω), and its component in the direction of the velocity field vanishes along the discontinuity; (ii) β = (0, 1 Ω 1 − 1 Ω 2 ) ⊤ (see Figure 6.5 (a)) and initial condition u 0 (x, y) = (sin(2πx), 4) ⊤ : The initial condition is curl-free, while its contraction with the velocity field, namely the vector component in the velocity direction is not in H 1 (Ω); (iii) β = (1 Ω 1 − 1 Ω 2 , 0) ⊤ (see Figure 6.5 (b)) and initial condition u 0 (x, y) = (sin(2πx), 4) ⊤ : The initial condition is curl-free, and its component in the direction of the velocity field vanishes along the discontinuity; (iv) β = (1 Ω 1 − 1 Ω 2 , 0) ⊤ (see Figure 6.5 (b)) and u 0 (x, y) = (x(2 − x)y(2 − y), sin(2πx)) ⊤ as initial condition: u 0 is in H(curl, Ω), and its component in the direction of the velocity field is not in H 1 (Ω). For the case (i), we run the simulation for an entire period, namely until T = 2 and compare the solution with the initial condition. Figure 6.6 shows that the initial datum is transported in the two different domains and the L 2 -error computed at the final time reaches the expected first order convergence (see Figure 6.7).  In case (iv), even if the initial condition is smooth, the magnetic advection with normally discontinuous collisional velocity yields the formation of a shock along the discontinuity ( Figure 6.8 (b)) until complete blow up. An analogous behavior of the numerical magnetic potential obtained from the stabilized scheme can be reported in the case (ii), where instantaneous blow up of the solution along the discontinuity is observed. However, we expect a blow-up of the solution in these situations: the observed behavior of the numerical solution is not engendered by instabilities produced within the numerical scheme, but accurately reflects "physical reality".

Orszag-Tang benchmark
As last numerical experiment, we rely on a widely used two dimensional benchmark problem, the so-called Orszag-Tang vortex system [33], which describes the transition to supersonic turbulence in the MHD equations. The development of shock waves and the complex interaction between various shocks with different speed which characterized the solution, makes the Orszag-Tang benchmark a challenging test for numerical methods.
Let us consider the two dimensional magnetic advection problem in Ω = [0, 2] 2 with periodic boundary conditions at the boundary ∂Ω. The time interval is I = [0, 1] with uniform time step τ = 5 · 10 −3 . The initial condition is the smooth vector field u 0 (x, y) = (sin(2πx), sin(πy)) ⊤ and the velocity field is piecewise constant with respect to the mesh. In particular, it is given at each time step as the outcome of a second order Finite Volume simulation of the full MHD system (from [32] and [21]). Note that even if the initial velocity is smooth, complex structures such as shocks and shock interactions develop in time. Concerning the discretization, the numerical scheme has been implemented on a tensor product mesh with 200 × 200 elements and H(curl, Ω)-conforming lowest order rotated Raviart-Thomas elements Λ 1 h (T h ) = Q − 1 Λ 1 (T h ) have been used for the spatial discretization, while a second order two-stage Runge-Kutta time-stepping is deployed in order to exploit the mild CFL-condition of the scheme. The stabilization parameter is as in (2.9), i.e. the upwind direction is assumed to be given by the average of the velocity field.
The rotated magnetic potential obtained using the scheme (4.1) has been compared with the magnetic induction derived from the full MHD system and a second order Finite Volume discretization (from [32] and [21]). The numerical method we have proposed well resolves shocks, is stable and no spurious oscillations occur. As can be inferred from Figure 6.10, the current sheet characterizing the second component of the magnetic induction (black box in Figure 6.10 (d)) is not captured by the H(curl, Ω)-conforming scheme due to the low order polynomial space discretization which is highly diffusive.  A Conforming Approximation of Non-Conforming Discrete Differential Forms T ∈ T h , the degrees of freedom of P r Λ k (T ) are defined on every j-simplex f j ∈ ∆ j (T ), with k ≤ j ≤ min{n, r + k − 1}, as ω ∈ P r Λ k (T ) −→ W ℓ f j (ω) := f j tr ω ∧ η ℓ j ∀ℓ = 1, . . . , N j , (A.1) where {η ℓ j } N j ℓ=1 is a basis of P − r−j+k Λ j−k (f j ) (see for example [26,Sections 3.2 and 3.4]). Similarly, the degrees of freedom of P − r Λ k (T ) are defined on every j-simplex f j ∈ ∆ j (T ), with k ≤ j ≤ min{n, r + k − 1}, as ℓ=1 is a basis of P r−j+k−1 Λ j−k (f j ). Finally, on a n-dimensional hypercube T ∈ T h the degrees of freedom of Q − r Λ k (T ) are defined on every j-face f j ∈ ∆ j (T ), with k ≤ j ≤ min{n, r + k − 1}, as . Furthermore, if f n−1 ∈ ∆ n−1 (T h ) then we rely on a trace property of the shape function space to infer tr ω ∈ P r Λ k (f n−1 ) which is hence determined by the degrees of freedom corresponding to the j-faces f j ⊂ f n−1 . Note that the interelement continuity required for the finite element space Λ k h (T h ) to belong to HΛ k (Ω) is directly imposed through the choice of degrees of freedom associated to the j-dimensional faces in T h , with k ≤ j ≤ n − 1.
Without loss of generality, in what follows we can assume that n < r + k − 1. Moreover, for the sake of conciseness, we restrict the analysis to the case Λ k h (T h ) = P r Λ k (T h ). The extension to the cases Λ k h (T h ) = P − r Λ k (T h ) and Λ k h (T h ) = Q − r Λ k (T h ) follows alike.
Lemma A.1. Let T be the reference n-simplex. Let φ : T → T be an affine isomorphism.
Then there exist C 1 , C 2 > 0 such that, for allω ∈ P r Λ k (T ), it holds Proof. The proof relies on standard transformation techniques. Details can be found in [26,Section 3.6].
Definition A.2. We define the averaging interpolation operator P h : P d r Λ k (T h ) → P r Λ k (T h ) through the degrees of freedom of ψ := P h ω for ω ∈ P d r Λ k (T h ) with degrees of freedom defined locally as in (A.1). In particular, • Degrees of freedom on n-simplices: Moreover, by definition of trace and jump across a face f n−1 = ∆ • n−1 (T + ) ∩ ∆ • n−1 (T − ), the degrees of freedom of tr [ω] f n−1 are given by the set {W ℓ f j ,T + − W ℓ f j ,T − } for j = k, . . . , n − 1 and ℓ = 1, . . . , j. Hence, applying Lemma A.1 yields .
Using a scaling argument yields the conclusion.

B A coordinate based representation of Lie derivatives
Let {e k I } ( n k ) I=1 be the orthonormal basis of alternating k-linear forms, e.g. ω = ( n k ) I=1 ω I e k I with 0-forms ω I for arbitrary ω ∈ Λ k (Ω) and e k I ∧ ⋆e k J = δ I,J . Then the projection of L β ω onto e J yields Hence, if u = (ω 1 , . . . , ω I , . . . ) ⊤ is a vector proxy of ω ∈ Λ k (Ω) then n i=1 β i ∂ x i u + Cu with C JI := δ J,I ( n i=1 ∂ x i β i ) + e k I ∧ ⋆L β e k J is a vector proxy of L β ω.