A unified analysis of algebraic flux correction schemes for convection–diffusion equations

Recent results on the numerical analysis of algebraic flux correction (AFC) finite element schemes for scalar convection–diffusion equations are reviewed and presented in a unified way. A general form of the method is presented using a link between AFC schemes and nonlinear edge-based diffusion schemes. Then, specific versions of the method, that is, different definitions for the flux limiters, are reviewed and their main results stated. Numerical studies compare the different versions of the scheme.


Introduction
Scalar convection-diffusion equations model the convective and molecular transport of a quantity like temperature or concentration. In applications, the convective transport is usually dominant, which is the case of interest in this paper.
Here, we consider the steady-state situation, where the mathematical problem is formulated as follows: Find u : Ω → R such that where Ω ⊂ R d (d = 2, 3) is a bounded polygonal or polyhedral domain with a Lipschitz continuous boundary ∂Ω, ε > 0 is a constant diffusion coefficient, b ∈ W 1,∞ (Ω) d is a solenoidal convection field, c ∈ L ∞ (Ω) is a non-negative reaction coefficient, g ∈ L 2 (Ω) is an outer source of the quantity u, and u D ∈ H 1 2 (∂Ω) ∩ C 0 (∂Ω) is a boundary datum. A characteristic feature of solutions of (1) is the appearance of layers, i.e., of narrow regions where the solution has a large gradient. These regions are usually so narrow that the layers cannot be resolved by affordable grids. It is well known that standard discretizations cannot cope with this situation and they lead to meaningless numerical solutions that are globally polluted with huge spurious oscillations. The remedy consists in using stabilized discretizations. In the context of finite element methods, the proposal of the streamline-upwind Petrov-Galerkin (SUPG) method in [15,22] was the first milestone in this direction. Solutions computed with this method usually have sharp layers at the correct position, but there are still non-negligible spurious oscillations in a vicinity of layers. Since the publication of [15,22] the development and analysis of stabilized discretizations for convection-dominated equations has been an active field of research.
In this research, one can distinguish two directions. The first one is the development of stabilized methods with a provable order of convergence in appropriate norms. Examples of this direction are the continuous interior penalty (CIP) method (see, e.g. [16]) and the local projection stabilization (LPS) method (see [12] for the first application of this method to a convection-dominated equation). The second direction consists in finding stabilized methods that compute solutions without spurious oscillations and still with sharp layers. The property of being free of spurious oscillations can be expressed mathematically with the satisfaction of the discrete maximum principle (DMP). Usually, the satisfaction of the DMP is proved by the sufficient condition that the matrix of a linear discretization 1 is an M-matrix. However, it is well known that, in the limit case ε = 0, there is a barrier of the order of the local discretization error for linear discretizations with M-matrices: these discretizations are at most of first order, e.g., see [42,Thm. 4.2.2].
Since the property of being free of spurious oscillations might be of utmost importance for a method to be applicable in practice, a significant amount of work has been devoted to the development of such methods. Due to the order barrier for linear discretizations, nonlinear discretizations became of interest. One further argument in favor of using nonlinear discretizations for a convection-dominated problem stems from the fact that most of the applications in which convection dominates are modeled by nonlinear partial differential equations. Then, the use of a nonlinear discretization does not constitute a significant overhead. Since the late 1980s, there have been a number of proposals to remove the spurious oscillations of the SUPG method by adding appropriate nonlinear terms. This class of methods is called spurious oscillations at layers diminishing (SOLD) methods, or shock capturing methods. A comprehensive review was carried out in the companion papers [25,26], and the main conclusion of it was that none of the proposed SOLD methods reduced the spurious oscillations sufficiently well.
Algebraic stabilizations, so-called algebraic flux correction (AFC) schemes, became of interest to us as a result of numerical assessments of stabilized discretizations in [5,25,28,29]. The main motivation for the design of AFC methods is the satisfaction of the DMP. In addition, they provide reasonably sharp approximations of the layers. In contrast to SOLD methods, which are based on variational formulations, the main idea of AFC schemes consists in modifying the algebraic system corresponding to a discrete problem, typically the Galerkin discretization, by means of solution-dependent flux corrections. Consequently, AFC schemes are nonlinear. The basic philosophy of flux correction schemes was formulated already in [14,43]. Later, the idea was extended to the finite element context, e.g., in [4,39]. In the last fifteen years, there has been an intensive development of these methods, e.g., see [33][34][35][36][37].
None of the above references deals with the mathematical analysis of the AFC methods. In fact, the first contributions to the numerical analysis of AFC schemes were presented only recently in [8,9,11]. The first paper [8] focuses on the solvability of the nonlinear scheme, while [9] presents the first error analysis of the AFC schemes. Interestingly, the paper [9] also presented negative results, in the sense that it was shown that unless some restrictions are imposed on the mesh, the numerical scheme may not converge. Finally, in the recent paper [11] the role of the linearity preservation was studied. This study is also complemented by the work [10], where a link between the AFC schemes and a nonlinear edge-based diffusion scheme is presented, and the linearity preservation of the scheme is also studied in detail. This latter reformulation offers the applicability of different tools than those used so far for the analysis of AFC schemes. In particular, it facilitated the a posteriori error analysis of the AFC method, presented in [2]. Thus, the present paper aims at providing a review of these works, and performing the analysis in a unified framework.
The rest of the manuscript is organized as follows. After having introduced AFC methods in Sect. 2, a rewriting as an edge diffusion scheme is presented. A unified analysis is given in Sect. 3, covering the existence of a solution, minimal conditions for the validity of the DMP, and finite element error estimates. Three definitions of limiters are provided in Sect. 4. Strategies for the solution of the nonlinear problems are discussed in Sect. 5. Numerical studies for different limiters used in AFC schemes and the edge diffusion scheme proposed in [10] are presented in Sect. 6. Finally, Sect. 7 states the most important open problems in the field of AFC schemes.

The model problem and a unified presentation of AFC schemes
The weak formulation of (1) reads: find u ∈ H 1 (Ω) such that u| ∂Ω = u D , and where (·, ·) Ω denotes the inner product in L 2 (Ω) or L 2 (Ω) d and the bilinear form a(·, ·) is given by Thanks to the Poincaré inequality, and the fact that b is solenoidal and c is non-negative, this problem has a unique solution.
To discretize the problem (1), we introduce the following notation: -{T h } h>0 denotes a family of shape regular simplicial triangulations of Ω.
-For a given triangulation T h , E h denotes the set of its internal edges.
-For every edge E ∈ E h , we denote by h E the length of E and by x E,1 , x E,2 the endpoints of E. Furthermore, for every E ∈ E h , we choose one unit tangent vector t E . Its orientation is of no importance. -For every edge E ∈ E h , we define the neighborhood ω E := ∪{T ∈ T h : T ∩ E = ∅}.
-For a given triangulation T h , {x 1 , . . . , x N } is the set of its nodes. We will assume that the nodes x 1 , . . . , x M are the internal nodes, and x M+1 , . . . , x N are boundary nodes, i.e., the nodes where the Dirichlet boundary condition is imposed. -For a node x i , i = 1, . . . , N , we define -For an interior node x i , i = 1, . . . , M, we define the index set of its neighbors -The finite element spaces used in this work are given by -These spaces have standard nodal basis functions denoted by {ϕ 1 , . . . , ϕ N }, uniquely determined by the conditions In addition, we set The Galerkin scheme associated to (2) is given as follows: This scheme is well known to lead to inaccurate results on affordable grids. The first step towards the building of an AFC scheme is the writing of the Galerkin method (4) in matrix form. For this, we introduce the matrix A = (a i j ) N i, j=1 , where a i j = a(ϕ j , ϕ i ). Then, we represent the discrete solution by a vector U ∈ R N of its coefficients with respect to the basis {ϕ 1 , . . . , ϕ N } of V h . Then U ≡ (u 1 , . . . , u N ) satisfies the following system of linear equations: where g i = (g, ϕ i ) Ω for i = 1, . . . , M. Thanks to the ellipticity of a(·, ·) on V h,0 , the matrix Using the matrix A = (a i j ) N i, j=1 , we introduce a symmetric artificial diffusion matrix The first step of defining an AFC scheme is then to add artificial diffusion to the algebraic system. More precisely, the problem (4) is replaced by In practice, the solution of such a perturbed scheme, which corresponds to simple upwinding, is too diffusive to be of interest. Then, the aim of AFC schemes is to localize this added diffusion in such a way that the DMP is respected, while the internal and boundary layers are not too smeared. This requires a finer analysis of the structure of the product D U. Since the row sums of the matrix D vanish, it follows that . Clearly, f i j = − f ji for all i, j = 1, . . . , N . Then, a further rewriting of (9) reads as follows: The next fundamental step in the building of an AFC scheme is to limit the fluxes f i j . In other words, the idea is to localize the diffusion to the areas surrounding extrema and layers. To this end, we introduce solution-dependent correction factors (or flux limiters) β i j ∈ [0, 1], and replace system (9) For β i j = 0, the original system (5) is recovered. Hence, intuitively, the coefficients β i j should be as close to 0 as possible to limit the modifications of the original problem. So far, these coefficients have been chosen in various ways, and their definition is always based on the fluxes f i j . To guarantee that the resulting scheme is conservative, and to be able to show existence of solutions, one should require that the coefficients β i j are symmetric, i.e., This requirement also has a mathematical justification. As a matter of fact, in [8], the possible non-existence of solutions has been shown if this restriction is ignored. Note that (11) does not involve β i j with i ∈ {M + 1, . . . , N } and hence these values can be chosen arbitrarily. We define them by the above symmetry condition and by the requirement that β i j = 0 if i, j ∈ {M + 1, . . . , N }.

A variational formulation and a rewriting as an edge diffusion scheme
Our starting point is the following variational formulation presented in [9] for problem (11), Here, the nonlinear form D h (·; ·, ·) is given by We now rewrite this nonlinear form using the symmetry of d i j and β i j : where we have denoted β E = β i j = β ji and d E = d i j = d ji for any edge E ∈ E h that has the endpoints x i and x j . Hence, with this rewriting of D h , we can state the following general form of an AFC scheme: Find u h ∈ V h such that u h | ∂Ω = i h u D , and where a h (z; v, w) = a(v, w) + D h (z; v, w) with a(·, ·) being defined in (3) and D h (·; ·, ·) given by Since, for any function from V h , the restriction to any edge E of E h is a linear function, one has From now on we will suppose that, for every edge E ∈ E h , d E is a real number, not necessarily linked to the matrix A. This flexibility will allow us to include a wider class of methods in our presentation.
The solution-dependent limiters β E are still assumed to satisfy β E ∈ [0, 1] and to assure the solvability of (14) (see the next section), we further make the following continuity assumption: It will be shown in Sect. 4 that the limiters defined in [9][10][11] satisfy Assumption (A1).

Remark 1
The fact that the restriction of the functions v and w to the internal edges is a linear function is what makes it possible to obtain the expression (16) for D h . This property also holds for unmapped Q 1 finite elements, and for mapped Q 1 finite elements on parallelepipeds, although in that case this methodology would lead to a completely different method, as the cross-terms would not be included in the method. The implications of this remark are the topic of current investigations and are to be reported elsewhere. On a related note to the previous point, AFC-related schemes using higher order elements combined with Bernstein basis functions have been developed recently in the work [38], but the full stability and error analysis of the methods is lacking.

General properties of the nonlinear scheme
In this section we present the main results associated to the nonlinear scheme (14). More precisely, we present results on its solvability, minimal conditions for the validity of the discrete maximum principle, and a first error estimate for the method. In the following section the conditions imposed herein will be checked for different definitions of the limiters β E .

Existence of solutions
Lemma 1 (Consequence of Brouwer's fixed-point theorem) Let X be a finite-dimensional Hilbert space with inner product (·, ·) X and norm · X . Let T : X → X be a continuous mapping and K > 0 a real number such that (T x, x) X > 0 for any x ∈ X with x X = K . Then there exists x ∈ X such that x X < K and T x = 0.
A proof of Lemma 1 can be found in [40, p. 164, Lemma 1.4]. Now, the existence of solutions for the nonlinear scheme (14) can be proved. (14)) If Assumption (A1) holds, then there exists a solution u h of (14).

Theorem 1 (Existence of a solution of
Proof For this proof only, we will consider constants C > 0 that may depend on the data of (1) and h. In addition, we will make use of a function u h,D ∈ V h , which is an extension of the boundary datum i h u D . Let us first define the nonlinear mapping T : Since a(·, ·) is a continuous bilinear form, Assumption (A1) implies that T is a continuous mapping. Next, from the definition of a(·, ·), it follows that, for any v h ∈ V h,0 , Moreover, (16) and the fact that Then, the definition of the operator T yields The terms involving u h,D are bounded next. The Cauchy-Schwarz and Poincaré inequalities lead to In addition, using the shape regularity of the mesh sequence, β E (·) ≤ 1, and the local trace inequality, one arrives at Finally, the application of the Poincaré and Young inequalities gives (14).

The discrete maximum principle
In this section we shall formulate general properties of the limiters β E under which the AFC scheme (14) satisfies the local and global DMP. The local DMP will be formulated on the patches Δ i defined in Sect. 2.
To prove the DMP, we make the following general assumption, which is a reformulation of an analogous assumption introduced in [32]. (14) with limiters β E satisfying Assumption (A2). Consider any i ∈ {1, . . . , M}. Then

Theorem 2 (Local DMP) Let u h ∈ V h be a solution of
Proof Let u h ∈ V h satisfy (14) and let us denote where Moreover, for i = 1, . . . , M, one derives These properties follow from the fact that To prove (19), let c = 0 in Δ i and assume that max Δ i u h > max ∂Δ i u h . Since the maximum of u h on Δ i is attained at vertices of the elements of T h making up Δ i , this means that u h (x i ) is the strict maximum of u h on Δ i . Then Assumption (A2) implies that the sum in (24) is non-negative. Since A i = 0 (see (23)) and a ii > 0 (see (22)), there is j ∈ S i such that a i j < 0 and hence the left-hand side of (24) is positive, which is a contradiction.
For proving (17), it suffices to consider the case A i > 0. Let us assume that max Δ i u h > max ∂Δ i u + h . Then again max Δ i u h > max ∂Δ i u h and also u i > 0. Like before, the sum in (24) is non-negative and since A i u i > 0, the left-hand side of (24) is positive, which is again a contradiction proving the assertion.
The implications (18) and (20) follow in an analogous way.

Theorem 3 (Global DMP) Let u h ∈ V h be a solution of (14) with limiters β E satisfying Assumptions (A1) and (A2). Then
If, in addition, c = 0 in Ω, then Proof The proof is based on the technique used in [31, Theorems 5.1 and 5.2]. Let u h ∈ V h satisfy (14) and let g ≤ 0 in Ω. Then the nodal values of u h satisfy (21) and, due to (7), one has Note that First, let us show that one has a i j = 0, which completes the proof of (30). Now we want to prove that the relations (21), (23), (29), and (30) imply (25) and (27). If c = 0 in Ω and hence N j=1 a i j = 0 for i = 1, . . . , M (see (23)), then (21) still holds if one adds a constant to all components of the vector (u 1 , . . . , u N ) so that one can assume that s > 0. If N j=1 a i j > 0, then s > 0 can be also assumed since otherwise (25) trivially holds. Thus, let s > 0 and let us assume that (27) does not hold, which implies that J ⊂ {1, . . . , M}. We shall prove that then Assume that (31) does not hold. Then, applying (23) and (30), one derives for any Thus, the matrix ( a i j ) i, j∈J is singular and hence there exist real numbers {v i } i∈J , not all equal to zero, such that i∈J a i j v i = 0 for j = 1, . . . , M. Consequently, the matrix ( a i j ) M i, j=1 is singular, which contradicts (29). Therefore, (31) holds and hence, denoting (21), (30), and (23) (note that the first inequality implies that r > 0). Hence, s ≤ r , which is a contradiction to the definition of J . Therefore (27) and hence also (25) holds. The relations (26) and (28) can be proved analogously.

Remark 2
The global DMP stated in (27) and (28) assures that the global maximum (resp. minimum) is attained on the boundary of Ω but does not exclude that it is also attained at internal nodes. Therefore, it is called the weak DMP. The strong DMP states that the global maximum (resp. minimum) is attained only on the boundary of Ω and easily follows if c = 0 in Ω and g is negative (resp. positive). Then one has . . , M}, then i ∈ J and hence it follows from (21), (23), and (30) that which is a contradiction. The statement (33) follows analogously.

An a priori error estimate
The error estimate will be proven using the following mesh-dependent norm where D h is defined in (15) and c 0 := inf ess Ω c .
Theorem 4 (Error estimate) Let us suppose that the solution of (2) belongs to H 2 (Ω) and that c 0 > 0. Then, there exists C > 0, independent of h and the data of (1), such that Proof We decompose the error in the usual way First, we notice that D h (u h ; ρ h , ρ h ) = 0, and then, standard interpolation estimates lead to To bound the discrete error e h we use the ellipticity of a(·, ·), the properties of D h (·; ·, ·), and the relations (14) and (2) to get Next, the continuity of a gives Moreover, since D h (u h ; ·, ·) is a symmetric positive semi-definite bilinear form, it satisfies the Cauchy-Schwarz inequality, which gives Combining the above relations proves the result.
A simple estimate of the consistency error D h (u h ; i h u, i h u) 1 2 is given in the following lemma.

Lemma 2 (Basic estimate of the consistency error) Denoting
Proof Using β E ≤ 1 and the shape regularity of T h implies that If d E is defined by (8) for an internal edge E with endpoints x i and x j , then which finishes the proof.
Lemma 2 shows that if d E is defined by (8), then the convergence order of u − u h h is reduced to 1/2 in the convection-dominated case and no convergence follows in the diffusiondominated case. It was demonstrated in [9] that these results are sharp. On the other hand, the results of [10,11] indicate that a better convergence behaviour in the diffusion-dominated case may be expected if the AFC scheme is linearity preserving, i.e., if the stabilization originating from the AFC vanishes in regions where the approximate solution is a polynomial of degree 1. This property can be formulated in terms of the limiters β E in the following way.

Assumption (A3)
The limiters β E possess the linearity-preservation property, i.e., The linearity preservation leads to an improved bound of the consistency error provided that the limiters satisfy the following Lipschitz-continuity assumption.

Lemma 3 (Improved estimate of the consistency error) Let the limiters β E satisfy Assumptions (A3) and (A4). Then
Proof The proof is a refinement of the technique used in [10,Theorem 4]. Let us write Then it follows from Assumption (A4) and the shape regularity of T h that, for any Consequently, Like in Lemma 2, one also obtains Using the last two estimates and applying Young's inequality, one obtains To bound the last term, we use the linearity preservation and the Lipschitz continuity of D E . More precisely, for a given E ∈ E h , we introduce the function i E u ∈ P 1 (ω E ) as the unique solution of the problem Using standard finite element approximation results (see [18]), i E u satisfies Outside ω E , the function i E u can be arbitrarily extended to a function from V h . In view of Assumption (A3), one has D E (i E u; i E u, i h u) = 0 and hence, using (34), (35), and the shape regularity of T h , one obtains which completes the proof.

The Kuzmin limiter
In this section we review the results obtained when implementing the method with the definition of the limiters proposed in [34]. In that work, the algorithm, originally proposed by Zalesak in [43] was adapted to the steady-state case, and exploited further. We refer then to this limiter as the Kuzmin limiter. This limiter has been used in numerous works, for example [5,9], where a detailed study of its performance for the convection-diffusion equation is carried out. The numbers d E in the definition of D h are given by (8) in this case. After having presented the Kuzmin limiter, we will show that it satisfies Assumption (A1) and, under an additional assumption on the matrix A, also Assumption (A2). Consequently, the nonlinear problem (14) possesses a solution and satisfies the discrete maximum principle. It was demonstrated in [11,Ex. 7.2] with the help of a numerical example that the AFC scheme with the Kuzmin limiter is not linearity preserving in general.
The definition of the coefficients for the Kuzmin limiter relies on the values where These values can be computed by performing a loop over all internal edges. After this loop, one defines R + i := min 1, If P + i or P − i vanishes, we define R + i := 1 or R − i := 1, respectively. At Dirichlet nodes, these quantities are also set to be 1, i.e., Then, for any i, j ∈ {1, . . . , N } such that a ji ≤ a i j , we set The final step consists in defining β E := 1 − α i j for any internal edge E ∈ E h having the endpoints There is an obvious ambiguity in the definition of β E if a i j = a ji . This ambiguity does not influence the resulting method if min{a i j , a ji } ≤ 0 since then d E = 0 and the respective term with β E does not occur in (15). To fulfill the condition min{a i j , a ji } ≤ 0, which also assures the DMP (cf. Lemma 5), it may help to replace the matrix corresponding to the reaction term by a lumped diagonal matrix, see [9].

Lemma 4 The Kuzmin limiter satisfies Assumption (A1).
Proof Let E be an internal edge that connects the nodes x i and x j . Then it suffices to show we can restrict these considerations to the situation that a ji ≤ a i j . Moreover, it suffices to consider d i j < 0 since otherwise α i j ≡ 1.
As first case,ū h ∈ V h such that f i j (ū h ) > 0 will be considered. Thenū i >ū j and hence f i j (u h ) > 0 in a neighborhood ofū h . Using (39), (37), and (36), we obtain The numerator and the denominator are continuous functions and by assumption, the denominator is positive in a neighborhood ofū h . Hence α i j is a continuous function atū h . In the same way, we get for the case f i j (ū h ) < 0 first thatū i <ū j and second the representation formula Using exactly the same reasoning as above, we conclude that α i j is continuous atū h in this case.
The last case is Remark 3 In [9] it was shown that the terms α i j (u h )(u j − u i ) are even Lipschitz-continuous. The proof of this property is based on the representations (40) and (41) of the coefficients α i j . The sums in these representations are Lipschitz-continuous and then one can show that the function which is obtained by multiplying these representations with (u j − u i ) is Lipschitzcontinuous, too.

Then the Kuzmin limiter satisfies Assumption (A2).
Proof Consider any u h ∈ V h , i ∈ {1, . . . , M}, and j ∈ S i . Let u i := u h (x i ) be a strict local extremum of u h in Δ i . We want to prove that If a i j ≤ 0, then (43) holds since (1 − α i j (u h )) d i j ≤ 0. If a i j > 0, then a ji ≤ 0 due to (42) and hence a ji ≤ a i j and d i j = −a i j < 0. Thus, if u i > u k for any k ∈ S i , then f i j > 0 and f ik ≥ 0 for k ∈ S i , so that α i j = R + i = 0. Similarly, if u i < u k for any k ∈ S i , then f i j < 0 and f ik ≤ 0 for k ∈ S i , so that α i j = R − i = 0. Since a i j + d i j ≤ 0, one concludes that (43) holds.

A limiter leading to linearity preservation and DMP on general meshes (BJK limiter)
Here we present a limiter recently proposed in [11] using some ideas of [35]. This limiter is designed in such a way that the AFC scheme satisfies the discrete maximum principle and linearity-preservation property on arbitrary meshes, which is a substantial improvement in comparison with the Kuzmin limiter. Like in the previous section, the numbers d E used in (15) are given by (8).
The definition of the limiter again relies on local quantities  (37) and (38) and one sets Finally, the limiters are defined by β E := 1 − min{ α i j , α ji } for any internal edge E ∈ E h having the endpoints x i , x j .

Lemma 6 The above limiter satisfies Assumptions (A1) and (A2).
Proof The validity of Assumption (A1) follows analogously as in the proof of Lemma 4. Let us prove Assumption (A2). Consider any u h ∈ V h , i ∈ {1, . . . , M}, and j ∈ S i and assume that u i := u h (x i ) is a strict local extremum of u h in Δ i . Then we want to prove that If d i j = 0, then a i j ≤ 0 and hence (44) holds. Thus, let us assume that d i j < 0. If u i > u k for any k ∈ S i , then f i j > 0 and u max i = u i so that P + i > 0, Q + i = 0 and α i j = R + i = 0. Since a i j + d i j ≤ 0, one obtains (44). If u i < u k for any k ∈ S i , (44) follows analogously.
The constants γ i can be adjusted in such a way that the linearity-preservation assumption (A3) is satisfied. In fact, it suffices to use such constants that It was proved in [11] that (45) holds with γ i = 1 if the patch Δ i is symmetric with respect to the vertex x i , and with

Remark 4
Note that large values of the constants γ i cause that more limiters α i j will be equal to 1 and hence less artificial diffusion is added, which makes it possible to obtain sharp approximations of layers. On the other hand, however, large values of γ i 's also cause that the numerical solution of the nonlinear algebraic problem becomes more involved.

A limiter based on the variation of the discrete solution (BBK limiter)
In this section we review briefly the limiter presented in [10] and its main results. This limiter, also referred to as smoothness-based viscosity, has its origin in the finite volume literature (see, e.g., [24] and [23]), and has also been used (although in a slightly modified way) in the recent work [19]. The numbers d E in the definition of D h are given by where γ 0 is a fixed parameter, dependent on the data of (1). The limiters β E , E ∈ E h , are given by the following algorithm: for w h ∈ V h , one defines ξ w h as the unique element in V h whose nodal values are given by The value for p determines the rate of decay of the numerical diffusion with the distance to the critical points. A value closer to 1 adds more diffusion in the far field, while a larger value makes the diffusion vanish faster, but on the other hand, increasing p may make the nonlinear system more difficult to solve. In our experience, values up to p = 20 are considered safe to use (see [10] for a detailed discussion). In this section we will detail the proof of the results using p = 1, but these extend to p > 1 without difficulty (see [10] for details). Two remarks can be rapidly made about this definition of the limiter. First, if a function v h has en extremum at an internal node x i , and E ∈ E i , then β E (v h ) = 1. This will be of paramount importance for the satisfaction of the DMP. Moreover, for meshes which have a certain structure, the method is linearity preserving, i.e., Assumption (A3) holds. More precisely, we will say that a mesh is symmetric with respect to its inner nodes if, for every node x i , and every j ∈ S i , there exists k ∈ S i such that x k − x i = −(x j − x i ). So, if the mesh is symmetric with respect to its internal nodes, if E ∈ E h has endpoints x i and x j , and v h ∈ P 1 (ω E ), then which gives β E (v h ) = 0. So, the method does not add extra diffusion in smooth regions, whenever the mesh is sufficiently structured.
Remark 5 In [10, Remark 1] a process to generate a method which is linearity preserving on general meshes is described. It involves a minimization process per node to determine a set of weights. The same results that hold for the method presented in this work hold for that variant.

Lemma 8 The limiter defined in this section satisfies Assumptions (A1) and (A4). Moreover, if the triangulation T h is such that
and γ 0 ≥ C 0 b ∞,Ω + C 1 c ∞,Ω h (where C 0 and C 1 are two constants independent of h, but large enough), then Assumption (A2) is fulfilled, too.
Proof To prove (A1) and (A4), let u h , v h ∈ V h . First, in [10, Lemma 1] the following result is proven: for any internal node x i the following holds Let us now suppose that, for E ∈ E h having endpoints In a completely analogous way one obtains which proves Assumption (A4) and hence also (A1).
To prove (A2) let us suppose that u h , solution of (14), has an extremum at the internal node x i . Let j ∈ S i , and let E ∈ E h be the edge with endpoints x i and x j . Then, as was mentioned earlier, β E (u h ) = 1. Thus, using the shape regularity of the mesh, one obtains , and Assumption (A2) follows.
It follows from Lemma 2 that the consistency error D h (u h ; i h u, i h u) can be bounded as follows: 2 1,Ω . Moreover, if the mesh is symmetric with respect to its internal nodes, then Lemma 3 implies that the following bound holds for the consistency error Thus, the method with the definition of the limiters from this section converges for every regular mesh, and, in addition, in the case in which the limiters are linearity preserving, the convergence order increases from O(h

Related recent work
We finish this section by mentioning that, to the research reviewed in this paper, work has been done in parallel, e.g., in [6,7,19,20]. In those references, a stabilizing term similar to the one defined in (15) is added to the formulation, and referred to as the Graph Laplacian.
The stabilizing mechanism of the methods presented in those works and the ones reviewed in this manuscript are very similar. There are, nevertheless, significant differences in the limiters β E , some of the results obtained in terms of the satisfaction of the discrete maximum principle, and the linearity preservation of the final schemes.
For example, in [6] the emphasis is in the regularization of the limiter proposed there (related to the BBK one) in order to make the limiter differentiable, to allow the use of Newton's method to solve the nonlinear system. The regularization proposed in there used regularization parameters that had an impact on the performance of the method. In addition, although the results concerning the discrete maximum principle were not too different from the ones reviewed in this work, the linearity preservation was not guaranteed for the regularized limiters. In [19] the definition of the limiter (non-differentiable this time) is modified using generalized barycentric coordinates in order to make it differentiable on meshes for which the support of basis functions is convex.
A final important difference between the works reviewed in this paper and the abovementioned references consists in the emphasis. While the papers just quoted deal with first order hyperbolic systems, the results reviewed in this paper deal with the convection-diffusion equation. The presence of the Laplacian in the partial differential equation makes the method satisfy very different properties, especially on non-Delaunay meshes, as it can be seen in [9] where an example of non-convergence was given for the Kuzmin limiter.

Iterative schemes for solving the nonlinear problem
Consider the weak formulation (13) and the equivalent formulation (11), (12) in matrix-vector notation. For simplicity, we will restrict the discussion to the case of homogeneous boundary conditions. These formulations represent a nonlinear problem since the coefficients β i j depend on the finite element solution u h . Applying an iterative scheme for solving the nonlinear problem, our experience is that usually damping is necessary to achieve convergence. Let u (m) h , m ≥ 0, be a given approximation of u h . A fixed point iteration can be defined as follows. In a first step, a finite element functioñ u (m+1) h is computed by solving: Findũ The matrix-vector form of (47) is N j=1 a i jũ where β (m) i j = β i j u (m) . In the iterations (47) and (48), the matrix of the problem changes in each iteration.
It is also possible to perform a fixed point iteration in such a way that only the right-hand side changes. Using the relation one can consider instead of (48) the iteration Using a sparse direct solver, then the matrix of (49) has to be factorized only once and in all subsequent iterations, only the solutions of the triangular systems have to be computed. Another approach for solving the nonlinear problem is a (damped) Newton method. Let us consider as starting point for deriving this method the matrix-vector formulation (11), (12). Let the i-th equation be written in the form

then the intermediate solution in Newton's method is computed by solving
where D F u (m) is the Jacobian, which can be computed by applying the product rule and the chain rule, and observing that the derivative of the limiter with respect to the Dirichlet nodes is not needed since these values are fixed Hence, the entries of the matrix that has to be inverted in (50) are given by  the particular limiter that is used in the simulations. Note that the derivation presented above requires smoothness of the limiters. If this property is not given, a possible approach consists in modifying the limiters such that they become sufficiently smooth, as is proposed in [6] for a transport equation. The exploration of strategies for applying Newton-type methods to (modified) AFC schemes for convection-diffusion-reaction equations is currently ongoing and will be reported elsewhere.
Let ω (m+1) ∈ [ω 0 , 1], ω 0 > 0, be a damping factor. The next iterate is given by The choice of appropriate damping parameters is essential for the efficiency of the iteration. In [26], an automatic strategy for adapting the parameter during the iteration is described.
In [5], the use of the so-called Anderson acceleration, proposed in [3,41], is advocated. The Anderson acceleration stores vectors from previous iterations and builds with them second order information.
6 Numerical studies

The Hemker example
We will consider the so-called Hemker example, which was proposed in [21]. It models the convection of temperature from a hot circle (2d cylinder) in a channel. The convection field is constant. There are exponential layers at the circle and interior layers downstream from the circle. The Hemker problem can be considered as a standard benchmark problem for convection-diffusion equations. It was used in [5] for comparing a number of stabilized discretizations. Here, the same setup as in this paper will be considered. This problem is defined in Ω = {(−3, 9) × (−3, 3)} \ {(x, y) : x 2 + y 2 ≤ 1}, the coefficients are b = (1, 0) T , c = 0, g = 0, and the boundary conditions are given by u(x, y) = 0, for x = −3, 1, for x 2 + y 2 = 1, and ε∇u · n = 0, elsewhere on the boundary.
In [5], a reference solution on a fine grid containing 48,252,416 degrees of freedom on a Q 1 mesh (for details, see [5]) was computed for ε = 10 −4 , see Fig. 1. Quantities of interest defined in [5] are the magnitude of the over-and under-shoots, the difference to the reference solution on selected cut lines, and the smearing of the interior layer at a certain cutline downstream from the cylinder. For the concrete definition of these quantities see [5]. The simulations were performed with P 1 finite elements on two types of grids, see Fig. 2. Grid 1 is aligned downstream to the convection and it has edges at the position where the interior layer is expected. The stopping criterion for the nonlinear iteration was based on the Euclidean norm of the residual vector, which should be smaller or equal than 10 −13 (#dof) 1 2 , where #dof is the number of degrees of freedom (including Dirichlet nodes) on the respective grid. As initial iterate, a function that vanishes on all degrees of freedom was used. The linear systems were solved with the sparse direct solver UMFPACK, [17]. The simulations were performed with the code MooNMD [27].
By construction of the Hemker example, the solution takes values in [0, 1]. The first quantity of interest from [5] considers the violation of this range by the numerical solutions. Since the AFC methods satisfy the DMP, it is expected that there are no violations if the nonlinear problems are solved exactly. In fact, we could observe in the numerical results only negligible violations of the order of the stopping criterion for the iteration of the nonlinear problem.
Another quantity of interest studies the smearing of the interior layer at x = 4, see Fig. 3. It can be seen that the smearing introduced by the BJK limiter is always smaller than with the Kuzmin limiter. In particular, on the aligned Grid 1, the results with the BJK limiter are much better. This statement is supported by considering the error to the reference solution at the cutline x = 4, see Fig. 4. To compute the errors, 10001 equidistant points were taken on the cutline and the vector e contains the differences of the reference solution and the numerical solution in these points. The errors in the Euclidean norm e 2 and the maximum norm e ∞ are given.
At the cutline y = 1, the results obtained with both limiters are similar, compare Fig. 5. The negative peak of the error is at the circle in a neighborhood of the point (0, 1). In this neighborhood, there is the transition from the exponential layer to the interior layer.
Finally, the costs for solving the nonlinear problems is studied. In the used code, only the fixed point iteration (49) is implemented. Either, the selection of the damping parameter as described in [26] can be used or the Anderson acceleration with l > 0 vectors and a fixed damping parameter ω. Results are presented for ω = 0.5 and l = 5, 10, 25 vectors in the Anderson acceleration. The numbers of iterations that were necessary for solving the nonlinear problems are illustrated in Fig. 6. It can be seen that generally fewer iterations were needed for the Kuzmin limiter. On the structured grid, the variant with 25 vectors in the Anderson acceleration needed often the smallest number of iterations and the fixed point iteration with an adaptive selection of the damping parameter needed most iterations. But on the unstructured grid, there is no clear picture. Using many vectors in the Anderson iteration did even result in failing to reach the stopping criterion on certain levels. For example, results have been obtained for other values of the damping parameter ω (not reported here due to space restrictions). More precisely, the use of ω = 0.25 gave similar results to the ones shown  Fig. 6, needing a few more iterations in most cases. Using ω = 0.7, we observed that the stopping criterion was not reached in many more cases than for ω = 0.5.
Altogether, these results show a bottleneck of AFC schemes that can hopefully be reduced or cured by using more advanced methods.

Illustration of the smearing of layers
A motivation for studying convection-diffusion equations in channel geometries comes from the simulation of population balance systems in chemical engineering. For experiments, chemical engineers often use long and thin pipes. That means, the diameter of the pipes is of the order of a few millimeters or centimeters and the length of the order of several meters. There are several specific properties when considering convection-diffusion equations in pipes or channels. First, a preferred flow direction exists. Second, the grids are eventually aligned with the flow direction and third, the mesh cells might be anisotropic. For convectiondominated problems there is the experience that it is of advantage to align the grid with the convection. In the literature, one finds already observations that report notable smearing of layers for algebraic stabilizations in examples where the grid is aligned to the convection, e.g., in [13,28].
This example considers a straight 2d channel, where the convection is a constant vector pointing into the direction of the channel. Let Ω = (0, 10) × (0, 1) and let ε = 10 −10 , b = (1, 0) T , c = g = 0 be the coefficients of the problem. The boundary condition is an impulse in a center strip at the inlet of the domain and there is a homogeneous Neumann boundary condition at the outlet: Because of the very small diffusion coefficient, one expects that the initial condition is transported from the inlet to the outlet. The coarsest grid is presented in Fig. 7. There are horizontal lines at both positions where the inlet condition has its jumps. The same stopping criterion for the solution of the nonlinear problem as in the Hemker example was used. Applying the SUPG method, one obtains a solution with sharp layers in the whole channel and with basically no spurious oscillations already on level 0, compare Fig. 8. In contrast, the solutions computed with the AFC schemes showed a notable smearing of the layers, in particular the solutions obtained with the Kuzmin limiter, see Fig. 9. One can see that the layers become sharper when refining the grid. The solutions computed with the BJK limiter are considerably more accurate than those obtained with the Kuzmin limiter, compare Figs. 9 and 10. We could observe that the solution for the Kuzmin limiter on level 3 looks similarly accurate as the solution of the BJK limiter on level 1.
The deeper understanding of the reasons for the smearing effect and the finding of remedies are open problems. So far, the probably best explanation is given in [28]. Algebraic stabilizations are by construction multi-dimensional schemes, i.e., there is no dimensional Fig. 11 The 3d example: the slice at z = 1 of the approximation obtained using the SUPG method on an adaptively refined mesh containing 1,308,237 elements splitting in the construction of the limiters. Such a splitting would be of advantage in this example since it is basically one-dimensional. However, the limiters see the layers of the solution that are vertical to the convection and they do not recognize that it is not necessary to introduce notable diffusion for preventing spurious oscillations.
Slices along the plane z = 1 of the solution obtained with the different methods are shown in Figs. 11, 12, 13, 14. To obtain a reference solution, we pursued this approach further and obtained a sequence of adaptively refined meshes using the same error estimator until a highly refined mesh, containing 135,408,953 elements, was built. A highly accurate (although not fully resolved) SUPG solution was computed in this mesh, and a slice of this solution along the plane z = 1 is presented in Fig. 15. Finally, in Fig. 16, we compare all these approximations and depict the cross section of them on the line y = z = 1.
There is a slight violation of the DMP for the method with the Kuzmin limiter. This violation is due to the fact that the mesh does not respect the hypotheses under which the DMP can be shown, cf. Lemma 5. For this mesh we have found violations of this condition, which explains the numerical results and confirms the sharpness of the analytical results. The boundary and inner layers are significantly sharper for the method with the BJK limiter, although this comes at the price of having to perform significantly more fixed point iterations than with the other methods. In fact, for this example the method using the Kuzmin limiter took 70 iterations to reach convergence, while the method using the BBK limiter took 166 iterations, and the use of the BJK limiter took 1117 iterations to reach convergence. As was mentioned earlier, the BJK and Kuzmin limiters provide sharper profiles than the BBK one. This has been observed not only in this example. This behavior seems to be related to the equal weight given to all fluxes by the construction of the BBK limiter, different from the Kuzmin one (essentially, an upwind limiter), and the BJK limiter, which has the flux associated to  the local extremum as main ingredient, and also includes some explicit mesh information to make it linearity preserving.

Open problems
The improvement of AFC schemes and the further development of their analysis have been listed in [30] among the most important open problems for H 1 -conforming finite elements for convection-diffusion equations. Some concrete issues are the following. It was shown by means of a numerical example that the general a priori estimate given in [9] is sharp. However, one can observe for the Kuzmin limiter and the BJK limiter higher orders of convergence than proved in [9], at least on special grids. So far, there is no concrete characterization of the necessary properties of such grids and no corresponding analysis. A priori analysis of AFC schemes for anisotropic grids remains an open problem. In addition, numerical analysis of AFC schemes for time-dependent equations is not available. Last but not least, efficient numerical methods for solving the nonlinear problems have to be developed.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.