Analysis of algebraic flux correction schemes for semi-discrete advection problems

Based on recent developments regarding the analysis of algebraic flux correction schemes, we consider a locally bound-preserving discretization of the time-dependent advection equation. Specifically, we analyze a monolithic convex limiting scheme based on piecewise (multi-)linear continuous finite elements in the semi-discrete formulation. To stabilize the discretization, we use low order time derivatives in the definition of raw antidiffusive fluxes. Our analytical investigation reveals that their limited counterparts should satisfy a certain compatibility condition. The conducted numerical experiments suggest that this prerequisite is satisfied unless the size of mesh elements is vastly different. We prove global-in-time existence of semi-discrete approximations and derive an a priori error estimate for finite time intervals with a worst-case convergence rate of 12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document} w. r. t. the L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{L}}^2$$\end{document} error. This rate is optimal in the setting under consideration because we allow all correction factors of the flux-corrected scheme to become zero. In this case, the algorithm reduces to the bound-preserving discrete upwinding method but the limited counterpart of this scheme converges much faster, in practice. Additional numerical experiments are performed to verify the provable convergence rate for a few variants of the scheme.


Introduction
Algebraic flux correction (AFC) schemes were proposed in [25] and have since then become an active research area [3,5,18,21,28]. These methods provide a robust framework guaranteeing that discrete maximum principles hold [6,27] and/or that entropy conditions are satisfied in the case of nonlinear equations [22,23]. Nonlinear AFC approaches combine a high order baseline scheme, such as the Galerkin finite element discretization, with a provably bound-preserving low order approximation. In this manner, global and/or local constraints can be imposed on the values of AFC solutions resulting from discretizations of various partial differential equations.
The focus of most efforts dealing with AFC schemes was on the development of numerical algorithms, while theoretical aspects have only recently started to attract significant interest. Barrenechea et al. [5] were the first to show solvability of a nonlinear AFC system arising from discretization of stationary convection-diffusion-reaction equations. Moreover, they prove that the scheme convergences with a rate of at least 1 2 in the AFC energy norm. In their subsequent work [6], they derived a sharper, first order error estimate under the assumption that the limiter is linearity preserving. Unfortunately, the proof technique that was used to obtain this superconvergence result relies on the presence of diffusive terms. Lohmann [27] extended the analysis of AFC schemes to the case without diffusive terms and obtained similar theoretical results for linear hyperbolic problems, again with a provable rate of 1 2 . Other theoretically investigated aspects of AFC procedures include their connection to edge-based diffusion [4], proofs of invariant domain preservation for the low order method [13], and a study of a posteriori error estimators [17]. The recent work of Jha and Ahmed [18] presents the first theoretical foundation of AFC schemes for parabolic convection-diffusionreaction equations. The AFC schemes analyzed therein are based on flux-corrected transport algorithms that are fully discrete and employ implicit time stepping.
In contrast to [18], this manuscript presents semi-discrete stability and a priori error analysis of AFC schemes for finite element discretizations of the time-dependent linear advection equation. To cure the oscillatory behavior of continuous Galerkin methods, we stabilize the antidiffusive fluxes using low order time derivatives (defined by (2.14) below). Flux correction is performed using Kuzmin's [21] monolithic convex limiting (MCL) scheme. For analytical purposes, we make an assumption regarding compatibility of the semi-discrete approximations and corresponding time derivatives. As shown in [15,Sec. 3.3], it is possible to enforce this condition by adapting the limiting procedure of the standard MCL approach. The results obtained in this manner are slightly more diffusive but exhibit the same second-order convergence rates in practice. However, based on our experience, the additional fix only rarely needs to be activated because compatibility seems to be automatically satisfied for the standard MCL scheme in most cases. Evidence for this claim is provided in Sect. 5.4. Therefore, we do not discuss enforcement of the compatibility condition in this work and instead refer the interested reader to [15,Sec. 3.3].
We prove that the nonlinear semi-discrete scheme is stable and that for finite times its spatial accuracy w. r. t. the L 2 norm is at least of order 1 2 for linear finite elements and generally unstructured meshes. In practice, second order superconvergence can be expected for smooth solutions and uniform meshes, as evidenced by the numerical examples of this paper and, for instance [3,18,28].
The structure of this article is as follows. We begin with the formulation of the continuous problem, and review the construction of the monolithic AFC schemes under discussion. Subsequently, in individual sections, we present the main theoretical outcomes of our work, an energy estimate and an a priori error analysis. Finally, we report the results of numerical experiments and draw conclusions. The contents of this paper are to a large degree based on [14,Ch. 5], which improves upon the analysis presented in our preprint [15]. In particular, we improved both the theoretical parts and the numerical examples by properly addressing the treatment of boundary conditions and presenting numerical results not just for simple 1D problems but also in the 2D case.

Discretization of the advection equation
In this section, we summarize the MCL strategy for linear transport problems [21]. Our presentation includes a brief discussion of the continuous model problem, a summary of the design principles of the low order method, as well as the formulation of the corresponding monolithic flux-correction schemes. Algebraic limiters of this kind have only recently been applied to different target discretizations such as high-order discontinuous Galerkin methods, e. g., [30]. These algorithms exploit ideas originally proposed in the context of continuous finite elements. It is therefore natural to perform our analysis in this framework as well.

Continuous model problem
Let Ω ⊂ R d , d ∈ {1, 2, 3} be a polyhedral domain, v ∈ C(Ω ×R + ) d a known velocity field, and n ∈ R d the unit outward normal to ∂Ω. We define the time-dependent inand outflow boundaries of Ω as In what follows, we suppress the dependence of Γ ± (t) on time t. The initial-boundary value problem for the linear advection equation reads which allows us to interpret (2.1a) as a hyperbolic conservation law with flux function f (u, x, t) = v(x, t)u. Let us remark that the flux correction tools discussed in this section can also be applied to problem (2.1) in the case of more general velocities.
To derive the weak formulation of (2.1), we multiply (2.1a) by a test function w and perform integration by parts. Replacing the consistent flux v · n u appearing in the resulting boundary integral with the upwind flux to incorporate the boundary dataû, we obtain With regard to the continuous weak formulation of (2. and define a weak solution to (2.1) as follows. Since our analysis is based on the same variational expression, we admit that the theory dictates this assumption on the model. However, we remark that some classical benchmarks for advection problems, such as LeVeque's solid body rotation [26] (see Sect. 5.3) do not satisfy (2.7).

Finite element discretization and low order method
The low order method that is employed in this work is the algebraic Lax-Friedrichs scheme [1,13,29,33] adapted to linear advection problems. In the AFC literature, this linear version is called the discrete upwinding method because of its equivalence to the node-centered upwind finite volume scheme [25,Sec. 6]. Let us now review the main steps of deriving this low order method. First, we discretize (2.4) in space using continuous linear finite elements.
be the vertices of the mesh and ϕ 1 , . . . , ϕ N be the corresponding piecewise linear Lagrange basis polynomials, satisfying ϕ i (x j ) = δ i j . For simplicity, we assume that the mesh has no hanging nodes. The corresponding finite element space shall be denoted as V h := {w h ∈ C(Ω) : w h | K ∈ P 1 (K ) ∀K ∈ K h } and the semi-discrete numerical solution is expanded as follows .
where m i j are scalar-valued entries of the consistent mass matrix Let us now briefly summarize the steps to construct the low order method used in [13,21], among others. We perform row sum mass lumping, i. e., replace the entries of M in the left hand side of (2.8) with those of In addition, we use the partition of unity property of basis functions, i. e., the fact that they sum to one everywhere in Ω to rewrite Here N i is the nodal stencil defined by where int(·) denotes the interior of a set and supp is the support of a function. Moreover, we add diffusive fluxes of the form As a final modification to (2.8), we employ a lumped approximation of boundary terms. This step involves a localization of boundary integrals to individual faces on the domain boundary.
Definition 2 (Nodal boundary faces, [14]) Let F ∂Ω denote the set of (d − 1)dimensional boundary faces of K h . Then the set F i contains all boundary faces that meet at node x i ∈ ∂Ω, i ∈ {1, . . . , N }. For Γ k ∈ F i , we defineû k i :=û(x i ) as the value ofû corresponding to Γ k ⊆ Γ − .
The above modifications made to (2.8) yield the low order method Note that b k i ≥ 0. We may also write (2.10) in the bar state form [13,16,21] where the bar statesū i j are defined bȳ (2.12) Remark 3 Definition (2.9) ensures thatū i j is a convex combination of u i and u j because Instead of (2.9), the classical version of the discrete upwinding method uses [25]

Monolithic convex limiting
The low order method (2.10) produces very diffusive approximations. To recover the accuracy of the standard finite element discretization (2.8), we perform algebraic flux correction. First, we define raw antidiffusive fluxes and their limited counterparts f * i j , which are specified below. Hereu h = N i=1u i ϕ i is a suitable approximation to the time derivative du h dt . Following [21], we employ the low order nodal valueṡ 14) to computeu h in practice. This approach can be interpreted as a modification of the target scheme corresponding to the standard continuous Galerkin discretization that otherwise exhibits a suboptimal first order convergence rate [31, Sec. 14.3.1]. As illustrated in Sect. 5.2, the use of low order time derivativesu i (instead of their consistent Galerkin counterparts defined by (5.2) below) also has a stabilizing effect on the overall approximation. This approach was proposed in the original publication on MCL schemes [21] and corresponds to a cheap and effective target scheme. Advanced stabilization techniques for higher-order methods applied to linear and nonlinear hyperbolic problems can be found in [28] and [23], respectively. By the above definition of raw antidiffusive fluxes, we have f i j = − f ji . To preserve the conservation property of the limited scheme, we enforce the corresponding constraint f * i j = − f * ji for the limited antidiffusive fluxes as is common in the AFC methodology [25]. Using the equivalence of formulations (2.10) and (2.11), we obtain a similar bar state form for the semi-discrete flux correction scheme [21] where the limited bar states are defined by [21] u Thus, a forward Euler update for (2.15) reads where Δt is the time step. Hence, the updated solutionũ i is a convex combination of u i , theū * i j , and theû k i , provided that the Courant-Friedrichs-Lewy (CFL) condition is satisfied. In other words, if (2.16) holds, the forward-Euler updated solutionũ i preserves all local bounds that these states are constrained by. This argument made here for a forward Euler step directly carries over to p-stage, pth-order accurate strongstability-preserving Runge-Kutta (SSP p-RK) methods [12,34], where p ∈ {1, 2, 3}.
In the process of flux correction, we enforce the local maximum principles in addition to skew symmetry of antidiffusive fluxes. Rearranging these constraints, we obtain Kuzmin's formula for the limited antidiffusive fluxes of his monolithic convex limiter [21] Lemma 1 (Conservation property of the MCL scheme, [21,25,33]) The semi-discrete scheme (2.15) Proof Summing over all degrees of freedom, we exploit the symmetry of diffusion coefficients d i j , skew symmetry of antidiffusive fluxes f * i j , and the zero row sum property of matrix A.

Remark 4 Continuous weak solutions u defined by (2.4) satisfy the conservation relation
Thus, (2.19) is a semi-discrete counterpart of (2.20) that accounts for the flux-lumped quadrature rule used in the AFC setting.
Let us now rewrite the bar state form (2.15b) of the semi-discrete MCL scheme in a formulation that is more amenable to theoretical investigations. Despite the fact that using MCL, the fluxes f * i j can be calculated directly via (2.18), we introduce correction factors The dependence of correction factors on the discrete solution makes AFC schemes nonlinear. Using the above definition of f i j , the semidiscrete MCL scheme (2.15) reads and can equivalently be written as

The bilinear and linear forms
are associated with the (stabilized) Galerkin finite element discretization corresponding to The nonlinear forms [5,18,27] in (2.22) are due to algebraic flux correction.

Energy estimate
Let us now derive an energy estimate for approximations obtained via (2.22). In the proof of this stability result, we rely on the assumption that the following requirement is satisfied.  This argument is our main motivation for relying on (3.1) for theoretical purposes. Clearly, the pair (u h , 0) satisfies (3.1). Thus if we do not compensate the mass lumping error in the process of limiting, the scheme automatically satisfies (3.1). As illustrated in Sect. 5.4, the standard MCL scheme using low order time derivatives (2.14) for stabilization purposes is also prone to producing compatible pairs. In [15,Sec. 3.3] we present a modified MCL procedure with which (3.1) can be guaranteed. Due to the complicated nature of this approach we chose not to discuss it any further in this work. Before presenting our energy estimate, we need to prove the following technical result.

Lemma 3 Any function
Proof Invoking the partition of unity property of basis functions, we obtain , the divergence theorem and Young's inequality to show that Multiplying by factor 2 and combining the integrals over Γ − , we write this inequality as Employing Lemma 3 and integrating in time produces (3.2).
Note that as a consequence of Lemma 2 and of the nonnegativity of basis functions, all terms appearing on the left hand side of inequality (3.2) are nonnegative. To guarantee that the assumptions of Proposition 1 are satisfied in practice, one can use the scheme proposed in [15,Sec. 3.3], which enforces (3.1) for user defined values of γ . In our experience, failure to apply this limiter has no negative practical effects, however.

Remark 5
The reader may wonder what significance is attached to Proposition 1. Since the fully discrete MCL scheme produces locally bound-preserving approximations, it is stable by design. Preservation of global bounds in the semi-discrete setting can be shown as in [24] under the assumption that a solution exists. The semi-discrete MCL scheme represents a nonlinear system of ordinary differential equations. Wellposedness of such initial value problems can be shown by invoking the Picard-Lindelöf theorem, which guarantees the existence of solutions on finite time intervals. Once local existence is established, we exploit a global existence and uniqueness result for ordinary differential equations [2,Thm. 7.6]. According to this theorem, solutions that cannot be extended to arbitrary times must in fact blow up, which, in our case, is prevented by Proposition 1. It follows that the semi-discrete MCL scheme (2.22) possesses a unique solution that exists for all times t ≥ 0.

Error analysis
Compared to the energy estimate derived in the previous section, our error analysis is rather involved. In particular, we need to make additional assumptions on the data of the continuous problem (2.4) as well as on the mesh sequences. These aspects are discussed in Sect. 4.1. Subsequently, in Sect. 4.2, we recall some auxiliary results from the literature on numerical analysis of finite element methods including AFC schemes. Finally, in Sect. 4.3, we state, prove, and discuss the main result of this work, which is a semi-discrete a priori error estimate for MCL approximations. Throughout this section, the letter C (possibly with a subscript) denotes a generic positive constant that is independent of the mesh size h. Moreover, we assume that h ≤ 1 and therefore h p ≤ h q for p ≥ q.

Preliminaries
Recall that we only consider meshes that are affine and geometrically conforming triangulations of Ω ⊂ R d , d ∈ {1, 2, 3}. Additionally, we restricted ourselves to simplicial meshes, which allows us to exploit the linearity of finite element approximations inside mesh cells.
The a priori error estimate that we present in Sect. 4.3 is valid only for quasi-uniform families of meshes, i. e., there has to exist C > 0 such that [8,Sec. 3.1.2] h := max where h K = diam(K ). As is standard in finite element analysis, we also assume shape-regularity of (K h ) h>0 . For this requirement to be satisfied, there has to exist C > 0 such that Ch K ≤ r K , where r K is the radius of the largest open ball that fits into K [8, Sec. 1.4.1]. Additionally, we assume that the mesh faces, which are simplices in R d−1 , are also shape regular in this sense. Our final assumption regarding the mesh sequence is that there exists C > 0 such that h ≤ Ch, whereh = min Γ ∈F ∂Ω diam(Γ ) and F ∂Ω is the set of boundary faces (cf. Definition 2). We do not need to assume contact regularity of the mesh sequence [8, Def. 1.38] as Di Pietro and Ern do because this requirement is automatically satisfied for simplicial triangulations.
Following [5,27], we assume H 2 (Ω) regularity of the exact solution u(·, t) for all t ≥ 0. We also require the time derivative ∂ t u to have this regularity. Specifically, we restrict our investigations to exact solutions of (2.4) that satisfy For simplicity, we set u h (·, 0) equal to the continuous interpolant Also for simplicity, we assume that the boundary dataû is linear on every boundary This assumption corresponds to a particular choice of the quadrature rule for boundary integrals.

Auxiliary statements
To prepare the ground for the derivation of our error estimate, we first summarize a few important ingredients of its proof, beginning with some standard inequalities.
Then we discuss aspects that are peculiar to algebraic flux correction schemes. Most of the AFC results were originally proven by Barrenechea et al. [5].

A priori error estimate
To state our main result, we need to define some auxiliary quantities. For t ≥ 0, let

Remark 6
We do not see any practical problems if (3.1) is invalid, only the theoretical results would no longer apply.

Corollary 1 (Convergence order of the semi-discrete MCL scheme) Under the assumptions of Proposition 2, the a priori error estimate
holds with a constant C 4 = C 4 (C 1 , C 2 , C 3 , T , u,û) > 0, which behaves as e C 1 T /2 .
Proof (of Corollary 1) Since q and z are nonnegative functions, we may use the estimate y(T ) ≤ C 2 h z(T ) in (4.1). The claim follows by calculating the integral of the exponential function.
Proof (of Proposition 2) This proof combines recent results on AFC schemes [5,27] with a new way of proving a priori error estimates for nonconforming discretizations of the advection equation [32]. A particular similarity of the approach developed in [32] to our theory is that both apply to semi-discrete formulations. We introduce the interpolation error Θ(t) = Θ(u, h; t) := u(t) − I h u(t) and subtract (2.22) from (2.4). Setting w = w h = ϑ h , we obtain the error equation Recall that the identity m i = N j=1 m i j holds for row sum mass lumping. Using this decomposition of m i and the identities u Arguing as in the proof of Proposition 1, we invoke the divergence theorem, Lemma 3 as well as the identities u = Θ + I h u and u h = I h u − ϑ h , which yields As in [19,Thm. 3.43], we exploit transformation to the reference element, shaperegularity, and the equivalence of norms in finite dimensional spaces to show that To derive an estimate for Ξ 3 , we rewrite the boundary integrals as a sum of integrals over faces. On each face Γ ∈ F − , we use the estimate |û k i −û| ≤ Ch Γ |∇û|, where h Γ = diam(Γ ). In addition, we invoke Young's inequality, estimate (4.3), Lemma 6, and incorporate λ = v L ∞ (Ω×R + ) d into the constant C, which yields For the nonlinear terms in Ξ 4 , we use Lemma 2, Young's inequality, the compatibility condition (3.1) with constant γ ∈ (0, 1), and Lemma 10 to deduce where the factor 1/γ was incorporated into the constant C. Combining the above identities for Ξ 1 and Ξ 2 with the inequalities for Ξ 3 and Ξ 4 produces the estimate (4.5) The terms on the right hand side of inequality (4.5) are now bounded using standard arguments. Specifically, we make use of the assumptions on the mesh and of Young's inequality, apply Lemma 4 to Θ = u − I h u and ∂ t Θ, invoke Lemma 5 through 9 and argue as in the derivation of (4.4) to obtain We now integrate in time observing that, by our definition of the discrete initial data, we have ϑ h (0) ≡ 0. At this stage, we recall the previously given definitions of q and z, which enables us to write the resulting inequality as Using Grönwall's Lemma as in [9, Lem. 1.9], we obtain The triangle inequality applied to u − u h = Θ + ϑ h then yields the error estimate (4.1) by Lemma 4.
We conclude the theoretical discussion with a few remarks regarding the derived error estimate. Let us first point out that in the general setting with unspecified correction factors α i j our result is indeed optimal (cf. Sect. 5.1). If we set all correction factors equal to zero, we obtain the low order method, which cannot be expected to be more than 1 2 order accurate in general. A drawback of our current approach is that the constant on the right hand side of the a priori error estimate (4.2) depends exponentially on the time T . Kučera and Shu [20] demonstrate that exponentially increasing constants can be avoided in some situations. They discretize the advection equation using discontinuous Galerkin methods and derive an error estimate without invoking Grönwall's inequality. It would be interesting to investigate the merit of their approach for the purposes of our analysis.
Let us briefly remark that we assumed all integrals appearing in the bilinear and linear forms a h (·, ·) and b h (·) to be evaluated exactly. In fact, even the energy estimate stated in Proposition 1 was derived under this assumption. For polynomial velocities one can indeed employ a quadrature rule of sufficiently high order to accurately compute all integrals. For general velocities, the theory we present needs to be adapted to include quadrature errors. As is common for linear finite elements [10, Thm. 8.5], we recommend to employ quadrature rules that are exact for polynomials in P 2 and P 3 for volume and boundary integrals, respectively.
Admittedly, a major limitation of Proposition 2 is the fact that the estimate is valid only for problems with exact solutions of very high regularity. In particular, the assumption that ∂ t u is H 2 in space is restrictive. In our opinion, the adaptation of the proofs in [5,27] to the time-dependent setting necessitates this regularity. One can argue that if the exact solution is smooth enough for Proposition 2 to be applicable, a limiter may not even be needed and we could instead employ a stabilized Galerkin method. Since this strategy does not guarantee the validity of discrete maximum principles, AFC schemes provide an appealing alternative. Therefore, theoretical investigations of these methods should be undertaken. It is hoped that our results may serve as a stepping stone for further efforts in this direction.

Numerical examples
Let us now corroborate the theoretical results of this work with numerical experiments. The following acronyms are used to distinguish the numerical methods under investigation • LOW: low order method (2.10), • MCL: monolithic convex limiting scheme (2.15), • target: target scheme, i. e., (2.15) with f * i j = f i j for all i ∈ {1, . . . , N }, j ∈ N i \{i}. Other methods used for comparative purposes are specified below. Recall that our stability and convergence proofs rely on the compatibility condition (3.1). In general, this condition is not fulfilled by the standard MCL approach. However, ifu h is set to zero, i. e., if the mass lumping error is not compensated, (3.1) holds due to Lemma 2. To distinguish between the standard MCL scheme employing low order time derivativeṡ u i given by (2.14) and the lumped-mass version, in whichu h ≡ 0, we employ the Table 1 Convergence history for the one-dimensional advection equation on a sequence of uniform periodic meshes. The · L 2 (Ω) errors at T = 1 and the corresponding EOC for u 0 (x) = exp(−100(x − 0.5) 2 ) acronyms MCL-L and MCL-0, respectively. Here the letter L stands for low order time derivatives, while 0 stands for zero time derivatives.
In all simulations, we choose the time step according to (2.16) by setting Δt = νΔt max , where ν ∈ (0, 1] is a user-defined value and Δt max is the right hand side of inequality (2.16). Discrete initial conditions are obtained by interpolating the continuous initial datum in the discrete space.
In the following sections, we verify that approximations converge at least as fast as the provable rate of 1 2 . Moreover, we stress the need for stabilization by the use of low order time derivatives (2.14) and present a comparison of results obtained with various definitions of antidiffusive fluxes in the MCL scheme. Finally, we perform an a posteriori check to see for which values of the parameter γ the compatibility condition (3.1) is satisfied by the MCL-L scheme.

Experimental orders of convergence
In this section, we solve the one-dimensional advection equation with constant velocity v = 1. The spatial domain Ω = (0, 1) has periodic boundaries. Thus, at each time instant T ∈ N 0 , the exact solution coincides with the initial condition. In this example, we use u 0 (x) = exp(−100(x − 0.5) 2 ).
We study the experimental orders of convergence for discrete upwinding (LOW), MCL-L, and MCL-0 schemes using SSP2-RK time stepping and CFL parameter ν = 0.5. While values as large as ν = 1 can safely be employed without causing violations of maximum principles, smaller values may be necessary to observe certain rates of convergence. Alternatively, SSP3-RK time stepping can be used to improve the temporal accuracy. In this study, we employ sequences of nested meshes with generally nonuniform mesh size h obtained by randomly perturbing the positions of the interior mesh vertices of the coarsest grid. The relative mesh sizes min K ∈K h h K /h of the three sequences are 1 (uniform), ≈ 0.69 (mildly perturbed), and ≈ 0.087 (severely perturbed), respectively. We present the L 2 (Ω) errors at the final time T = 1 and the corresponding experimental orders of convergence (EOC) in Tables 1, 2, 3.
The observed rates are in accordance with our expectations. As suggested by Corollary 1, discrete upwinding converges at least with the rate of 1 2 . Actually, the low order method becomes first order accurate on very fine uniform meshes. Our preferred MCL-L scheme produces second order accurate results in this test. If no correction of the mass lumping error is performed, the order of accuracy deteriorates, while still exceed- ing the provable rate of 1 2 . In this example, the influence of mesh perturbations on the results is insignificant. A decay in the convergence rate of the standard MCL scheme for a steady problem was observed on perturbed 2D meshes in [21, Sec. 6.1].

On the stabilizing effect of low order time derivatives
Let us now compare the standard Galerkin approach to methods that are stabilized by incorporating low order time derivatives (2.14) via antidiffusive fluxes. We consider the same setup as in the previous section with the exception that the initial condition u 0 is replaced by [15] otherwise.
This profile features discontinuities as well as a C ∞ region. In Fig. 1 we display standard continuous Galerkin approximations obtained with four different combinations of time stepping schemes and CFL parameters on a uniform, a mildly perturbed and a severely perturbed mesh with 128 elements in each case. Spurious ripples that are not local to the vicinity of the discontinuities can be observed in all profiles. Although limiters can remove these oscillations, the quality of approximations obtained in this We observe significant improvements in the solution quality for the unlimited target scheme compared to the consistent Galerkin approximations displayed in Fig. 1. Numerical results obtained with LOW, MCL-L and MCL-0 exhibit behavior similar to that observed in Sect. 5.1 In particular, the MCL-0 scheme produces a nonsymmetric profile in the left part of the domain, which can be attributed to dispersive errors that occur commonly if the mass lumping error is not compensated [35].

Solid body rotation
Let us now apply the MCL scheme to a 2D solid body rotation benchmark [26] in which Ω = (0, 1) 2 , v(x, y) = 2π (0.5 − y, x − 0.5) T ,û = 0 and In this example, a cone, a smooth bump and a slotted cylinder rotate around the domain center. At each time instant T ∈ N 0 , the exact solution is equal to the initial condition, which is shown in Fig. 3. The numerical results displayed in this section are visualized with the open source C++ software GLVis [11]. We solve this problem numerically using triangular meshes and h = c/128, where c = √ 2 for uniform grids and c = 1 for unstructured ones. For time stepping we employ the SSP2-RK method with constant time steps Δt = 5 · 10 −4 and Δt = 3.125 · 10 −4 , respectively. In addition to MCL-L and MCL-0 schemes, we test the MCL-G approach, in which the consistent Galerkin time derivativeu is employed to correct the mass lumping error through the antidiffusive fluxes The numerical results of this study are displayed in Although all obtained profiles appear to be similar to each other, we can make out some differences by closely examining the results. First, we observe that on the uniform mesh MCL-0 is noticeably more diffusive than MCL-L and MCL-G. In particular, the smooth hump is not well resolved in this approach due to dispersive errors that arise in lumped Galerkin approximations [35]. The unstabilized MCL-G scheme does not suffer from this deficiency and, in fact, produces the smallest approximate error values among the three approaches. However, the lack of stabilization leads to a distortion in the shape of the sharp cone and, on the uniform mesh, a similar feature can be spotted in the region of the slotted cylinder. In the case of the advection equation, this issue does not seem to have a dominating effect on the obtained profiles but, in our experience [14,Sec. 3.4.3.2], the MCL-G scheme produces more pronounced spurious oscillations if applied to more involved problems like the Euler equations of gas dynamics. The quality of approximations obtained in this manner can be improved by employing smaller time steps or higher order SSP-RK methods [21, Fig. 2 (e)]. Nevertheless, some form of stabilization should be used in combination with continuous Galerkin discretizations of hyperbolic problems. Therefore, MCL-L is the preferable option among the three schemes under investigation. Alternative stabilization techniques for standard finite elements can be found in [23,27] for linear and nonlinear problems, respectively.

A posteriori compatibility check
The compatibility condition (3.1) turned out to be an invaluable tool for our theoretical investigations. Unfortunately, we are unable to prove that the MCL-L scheme automatically produces compatible pairs (u h ,u h ) under suitable assumptions on the mesh. However, it is easy to check for which values of γ ∈ (0, 1) condition (3.1) is fulfilled a posteriori. Indeed, (  if the denominator in the right hand side of (5.3) is nonzero (it is nonnegative due to Lemma 2). If the numerator in (5.3) is also nonnegative, this criterion yields an upper bound on γ . Having calculated these a posteriori bounds via (5.3), one can check how they behave upon mesh refinement. Two issues that lead to a violation of compatibility can occur in practice. First, (5.3) can produce a negative upper bound on γ , a case that is not covered by our theory. Secondly, γ may approach zero upon mesh refinement, which would cause the constant in the leading order term of our error estimate to approach infinity. We found the former concern to be valid on perturbed one-dimensional meshes. Using smaller time steps does not resolve incompatibility issues, which seem to be caused by triangulations of bad quality. In such instances, our stability and error estimates are not applicable to MCL-L but remain valid for the LOW and MCL-0 schemes, as well as for the method proposed in [15,Sec. 3.3] that enforces compatibility.
Having performed the described a posteriori check for various test problems, we conjecture that compatibility of (u h ,u h ) holds for the MCL-L scheme on uniform meshes. Below we report the results of our experiments in which we compute the values of the right hand side of (5.3). First, we consider four one-dimensional test problems. In each case, the spatial domain Ω = (0, 1) is equipped with periodic boundaries and the velocity is v = 1. The first and second tests are the same as in Sects. 5.1 and 5.2. In the third and fourth tests, the final time is T = 0.5 and the initial conditions read u 0 (x) = 0.5 1 + cos π 0.15 (x − 0.25) if |x − 0.25| < 0.15, 0 otherwise, and u 0 (x) = max{0, 1 − 10|x − 0.2|}, respectively. We solve each of these problems on a hierarchy of uniform meshes using SSP2-RK time stepping with CFL parameters ν ∈ {1, 0.1}. The largest value of γ for which (3.1) is satisfied during the whole simulation is presented in Table 4.
Additionally, we repeat the solid body rotation test [26] from Sect. 5.3 on sequences of uniform (c = √ 2) and unstructured (c = 1) triangular meshes and compute a pos-  Table 5, where #TS refers to the total number of employed time steps. We observe slightly larger maximum values of γ on coarse girds than on fine meshes. The use of smaller time steps seems to have marginal influence on the results. In all cases, we have γ > 0.4, which is consistent to the value that we used in [15] to enforce (3.1). Contrary to the 1D case, (5.3) does not produce negative values for γ even on unstructured meshes in 2D. This observation leads us to believe that the low order time derivativeu h given by (2.14) is compatible to u h even for a certain class of unstructured meshes. Further theoretical and numerical studies are required to pinpoint, exactly which conditions a sequence of unstructured meshes has to satisfy in order to produce compatible pairs (u h ,u h ) via the standard MCL approach. The opposite point of view is that the compatibility condition (3.1) can actually be used to determine the mesh quality for mesh optimization purposes. Feasibility and benefits of this approach are yet to be determined. aspects of inexact numerical integration may need to be taken into account. We invite the interested reader to participate in these research endeavors.