Edge-based nonlinear diffusion for finite element approximations of convection-diffusion equations and its relation to algebraic flux-correction schemes

For the case of approximation of convection--diffusion equations using piecewise affine continuous finite elements a new edge-based nonlinear diffusion operator is proposed that makes the scheme satisfy a discrete maximum principle. The diffusion operator is shown to be Lipschitz continuous and linearity preserving. Using these properties we provide a full stability and error analysis, which, in the diffusion dominated regime, shows existence, uniqueness and optimal convergence. Then the algebraic flux correction method is recalled and we show that the present method can be interpreted as an algebraic flux correction method for a particular definition of the flux limiters. The performance of the method is illustrated on some numerical test cases in two space dimensions.

The weak problem (1.2) has a unique solution u ∈ H 1 (Ω) and its solution satisfies the following maximum principle (see [10]).
Definition 1 (Maximum Principle) Assume that f ≥ 0, g ≥ 0 (resp.≤ 0) and the solution u of (1.2) is smooth enough.Then, if σ = 0 and u attains a strict minimum (resp.maximum) at an interior point x ∈ Ω, then u is constant in Ω.If σ > 0, then the same conclusion remains valid if we suppose in addition that u(x) < 0 (resp.u(x) > 0).This work deals with the development of a method that satisfies the discrete analogous of the last definition.The quest for such a method has been a constant for the last couple of decades.Several methods have been proposed over the years, both in the finite element and finite volume contexts (see [21] for a review).Overall, the common point of all discretisations that satisfy a discrete maximum principle (DMP) is that they add some diffusion to the equations.This extra diffusion can lead to a linear method, but in that case it is a well-known fact that a linear method will provide very diffused numerical solutions, which will converge suboptimally.Due to the previous fact, several methods that add nonlinear diffusion have been proposed.
One approach has been to add a so-called shock-capturing term to the finite element formulation.This typically amounts to a nonlinear diffusion term where the diffusion coefficient depends nonlinearly on the finite element residual, making it large in the zones where the solution is underresolved, but vanish in smooth regions.An analysis showing that nonlinear shock capturing methods may lead to a DMP was first proposed in [5], and then developed further for the Laplace operator in [6], and for the convection-diffusion equation in [7].For a review of shock capturing methods, designed to reduce spurious oscillations, without necessarily satisfying a DMP, see [14].More recent nonlinear discretisations, these ones based on the idea of blending in order to satisfy the DMP, are the works [9,1], where the emphasis has been given to prove the convergence to an entropy solution.Most shock capturing techniques suffer from the strong nonlinearity introduced when the diffusion coefficient is made to depend on the finite element residual (and therefore the gradient of the approximation function).Because of this the analysis of such methods is incomplete even when linear model problems with constant coefficients are considered.In particular, in most cases uniqueness of solutions can not be proved, and the convergence theory is incomplete.
On the other hand, driven initially by the design of explicit time stepping schemes for compressible flows, so called Flux Corrected Transport (FCT) schemes and the related algebraic flux correction (AFC) schemes were introduced [20,19,15].These schemes act on the algebraic level by first modifying the system matrix so that it has suitable properties to make the system monotonous, while perturbing the method as little as possible.This crude strategy, however, necessarily results in a first order scheme.Then a nonlinear switch, or flux limiter, is introduced making the low order monotone scheme active only in the zones where the DMP may be violated.These schemes have also resisted mathematical analysis for a long time, but a number of results have been proved recently in [3,2].Indeed, in these references, existence of solutions and positivity have been proved, and a first error analysis has been performed.Nevertheless, it was shown that the DMP, and even the convergence of the discrete solution to the continuous one, depend on the geometry of the mesh.
Another approach to combine monotone (low order) finite element methods with linear diffusion and high order FEM using flux-limiters was proposed very recently in [13].It then appears that a cross pollination between the idea of AFC and shock-capturing could be fruitful.
The objective of the present paper is to further bridge the gap between the shock capturing approach and the algebraic flux correction.Indeed we will consider a generalisation of the shockcapturing term first introduced in [4] to several dimensions, using an anisotropic diffusion operator along element edges similar to that introduced in [7].We show that the resulting scheme satisfies the DMP and give an analysis of the method.In particular we show that the new shock capturing term is Lipschitz continuous, and, if the mesh is sufficiently regular, linearity preserving (see § 2.1), which allows us to improve greatly on previous results.In § 2.2 we prove existence of solutions, the discrete maximum principle, and noticeably, uniqueness in the diffusion dominated regime.We then show error estimates, which, thanks to the combined use of linearity preservation and Lipschitz continuity, turn out to be optimal in the diffusion dominated regime, for a special class of meshes (see § 3).In § 4, we revisit the design principles of AFC and show that the proposed shock-capturing term can be interpreted as an AFC scheme using a special flux, allowing both for a DMP and Lipschitz continuity.Some numerical results are finally shown in § 5.

Notations
We now introduce some notation that will be needed for the discrete setting.We consider a family {T h } h>0 of shape-regular triangulations of Ω consisting of disjoint d-simplices K.We define h K := diam(K), and h = max{h K : K ∈ T h }.We associate with the triangulation T h the finite element spaces where P (D) is the space of polynomials of degree at most on D. The nodes of T h are denoted by {x i } N i=1 , and the usual associated basis functions of V h are denoted by {ψ i } N i=1 .We let E h be the set of the interior edges of T h .For every edge E ∈ E h , we define h E := |E| and ω E := {K ∈ T h : K ∩ E = ∅}, and fix one unit tangent vector, denoted by t.
For an interior node x i , we define the neighborhoods E i := {E ∈ E h : x i ∈ E}, Ω i := {K ∈ T h : x i ∈ K}, and the set S i := {j ∈ {1, . . ., N } \ {i} : x j shares an internal edge with x i } . (1.4) Finally, we will say that the triangulation T h is symmetric with respect to its internal nodes if for every internal node x i the following holds: for all j ∈ S i there exists k ∈ S i such that x j − x i = −(x k − x i ) (see Figure 1 for examples in two space dimensions).

The nonlinear discretisation
The standard finite element method for the problem (1.2) takes the form: Find Here, u bh ∈ V h is introduced to approximate the boundary condition g.Then, we propose the following stabilised method to discretise (1.2): Find 2) The stabilisation term d h (• ; •, •) is defined by where γ 0 > 0, and α E : with ξ w h ∈ V h given by (2.5) The value for p will determine the rate of decay of the numerical diffusion with the distance to the critical points.A value closer to 1 will add more diffusion in the far field, while a larger value will make the diffusion vanish faster, but on the other hand, increasing p may make the nonlinear system more difficult to solve.In principle, as p goes to infinity the method will add the perturbations only in points with local extrema.In our calculations we have tested several different values for p, and have presented those for p = 1 and p = 4, the latter giving better numerical results, while keeping the nonlinear solver converging within a reasonable number of iterations.
Remark 1 All the proofs presented in this work are valid in both the two and three-dimensional cases.Just to simplify the presentation we will restrict ourselves from now on to the two-dimensional case d = 2.

Properties of d h (•; •, •)
We start noticing that This prevents the method from adding artificial diffusion to the equations in regions in which the solution is constant.Moreover, the method is as well linearity preserving if the mesh is symmetric with respect to its interior nodes.In fact, if E ∈ E h has endpoints x i and x j , and which gives α E (v h ) = 0.Then, the method does not add extra diffusion in smooth regions, whenever the mesh is sufficiently structured.The next step is to show that d h (•; •, •) is continuous.More precisely, it is Lipschitz continuous, and the next result is the first step towards this.
Lemma 1 For any v h , w h ∈ V h , and any given internal node x i , the following holds otherwise the claim is obvious.A quick calculation gives The following estimate can be proved in an analogous way Then, which gives the desired result upon applying the estimate min{a −1 , b −1 } ≤ 2 a+b , for two positive numbers a and b.
The Lipschitz continuity of d h (•; •, •) appears then as a consequence of the previous result.

Lemma 2
The nonlinear form d h (•; •, •) is Lipschitz continuous.More precisely, there exists C lip > 0, independent of h, such that, for all v h , w h , z h ∈ V h , the following holds (2.9) Proof We have The first term in the above estimate is bounded using the fact that |α E (v h )| ≤ 1, the Cauchy-Schwarz inequality, a local trace inequality, and the mesh regularity, to give The second term is bounded next.For this, a general edge E ∈ E h will be considered as having x i and x j as endpoints, where x i is chosen to be the vertex such that α where and the second term in (2.10) reduces to We now remark that for two numbers a, b ∈ [0, 1] we have and the term in E 1 is bounded using Lemma 1.In fact, from the mesh regularity there exists C > 0, independent of h, such that for all E, F ∈ E i , h F ≤ Ch E .Moreover, the number of edges in E i is uniformly bounded, independently of h.Then, using Cauchy-Schwarz's inequality and a local trace inequality we arrive at The sum over E 2 is bounded next.First, using (2.12) we get In an analogous way we obtain since the last term in the middle inequality is always non-positive, since by construction, for The result then follows collecting (2.10)-(2.13).

Solvability of the discrete problem
This section is devoted to analyse the existence of solutions for (2.2).It is interesting to remark that, thanks to the Lipschitz continuity of d h (•; •, •), the solution can be proved to be unique in the diffusion-dominated regime.
Lemma 3 Let T h : V 0 h → V 0 h be the operator defined by Then, where c 1 , c 2 are positive constants independent of z h , f , and g.
Proof For this proof only, we will consider constants C > 0 that may depend on the the physical coefficients.From the definition of a it follows that Moreover, the definition of d h (•; •, •) and the fact that 0 ≤ α E (z h + u bh ) give Then, the definition of the operator T h gives Next, the Cauchy-Schwarz and Poincaré inequalities lead to the following bound In addition, using the mesh regularity, α E (•) ≤ 1, and the local trace inequality, we arrive at We can thus conclude that The claimed result arises by applying the Poincaré and Young inequalities to the last relation.
The solvability of the nonlinear problem (2.2) appears as a consequence of the above result and Brower's fixed point Theorem.
Theorem 1 The discrete problem (2.2) has at least one solution.Moreover, if C lip γ 0 h < ε, where C lip is the constant from Lemma 2, then the solution is unique.

The discrete maximum principle
This section is devoted to prove that Method (2.2) preserves positivity.For this, we will impose the following geometric hypothesis on the mesh.This hypothesis can be tracked back to [22], and in two space dimensions it reduces to impose that the mesh is Delaunay.
Assumption 1 [Hypothesis of Xu and Zikatanov, cf.[22]] For every internal edge E ∈ E h with end points x i and x j the following inequality holds where θ K ij is the angle between the two facets in K opposite to x i and x j (denoted by F i,K and F j,K , respectively), and ω K ij is the (d − 2)-dimensional simplex F i,K ∩ F j,K opposite to the edge E.
We now introduce the discrete analogous of the maximum principle.This definition is very close to the one from [7], and it leads to results which are, essentially, identical to those from that reference.
Definition 2 (DMP) The semilinear form ã(•; •) is said to satisfy the strong DMP property if the following holds: For all u h ∈ V h and for all interior vertices x i , if u h is locally minimal (resp.maximal) on the vertex x i over the macro-element Ω i , then there exist negative quantites (c . Furthermore, we will say that the semilinear form satisfies the weak DMP property, if the local minimum above is supposed to be negative (resp. the local maximum is supposed positive).
A direct consequence of this definition is the following result, whose proof is completely analogous to the one from [7, Proposition 2.5].
Lemma 4 Assume that the semilinear form ã(•; •) satisfies the DMP property.Assume that u h ∈ V h solves (2.2) and that f ≥ 0.Then, u h reaches its minimum on the boundary ∂Ω.

The following result states the DMP for (2.2).
Theorem 2 Let us suppose that the mesh T h satisfies Assumption 1, and that the parameter γ 0 is large enough.Then, the semilinar form ã(•; •) satisfies the DMP property.
Proof Let us suppose that u h has a negative local minimum at an interior node x i .Then, α E (u h ) = 1 for all E ∈ E i , which gives (2.26) We will analyse the expression above term-by-term.First, if u h ≤ 0 in the support of ψ i , then (σu h , ψ i ) Ω ≤ 0. Let us suppose now that u h changes sign in the support of ψ i , and let K ∈ Ω i be an element in which u h changes sign.Let x j and x k be the other nodes of K, E ij and E ik the correspoding edges, and xij ∈ E ij , xik ∈ E ik the points in those edges in which u h vanishes (if u h vanishes in only one edge, then the proof is analogous).Then, using a quadrature formula and a Taylor expansion we obtain and, adding up over all K ∈ Ω i and using the mesh regularity we obtain Also, as in [7] (see also [21]), Assumption 1 on the mesh leads to ε(∇u h , ∇ψ i ) Ω ≤ 0 .
(2.28) Moreover N j=1 ψ j = 1 gives j∈Si (b • ∇ψ j , ψ i ) Ω = 0, and then which, using the mesh regularity gives Finally, since u h (x i ) is a local minimum, then in every E ∈ E i , ∂ t u h and ∂ t ψ i have different signs (independently of the orientation of the tangential vector in E), which gives Hence, gathering all the above computations, we arrive at and the result follows assuming that γ 0 > C 0 σh E + C 1 b ∞,E .Finally, we notice that if σ = 0 then the sign of the strict minimum is irrelevant.
Remark 2 It is interesting to remark that the hypothesis on the meshes of the triangulation can be avoided if the problem is supposed to be strongly convection-dominated.In fact, following analogous steps to those used to prove (2.29) we can arrive at Replacing this into the steps leading to (2.32) gives and the proof follows by assuming that E stays bounded, which means this is applicable only in the case the problem is highly convection-dominated.In this sense, the method proposed in this work can be applied to scalar conservation laws, regardless of the geometrical impositions on the mesh.Similar results have been obtained recently in [12,13].

Convergence
The error will be analysed using the following mesh-dependent norm: As usual, the error e := u − u h is split as follows where i h : C 0 (Ω) ∩ H 1 0 (Ω) → V 0 h stands for the Clément interpolation operator.Using standard interpolation estimates (see [8]), the fact that α E (•) ≤ 1, and the mesh regularity, the following bound for ρ h follows: 3) It only remains to estimate the discrete error.This is done in the next result.
Lemma 5 Let us suppose u ∈ H 2 (Ω) ∩ H 1 0 (Ω).Then, there exists C > 0, independent of h and ε, such that Proof First, from the definition of a and d h we get Next, the continuity of a gives Moreover, since d h (u h ; •, •) is a symmetric positive semi-definite bilinear form it satisfies Cauchy-Schwarz's inequality, which gives Then, inserting (3.6) and (3.7) into (3.5), and using Young's inequality, we arrive at It only remains to bound the consistency error d h (u h ; i h u, i h u) in (3.8).The definition of d h (•; •, •), α E (u h ) ≤ 1, a local trace inequality, the mesh regularity, and the H 1 (Ω)-stability of i h , give Then, the result arises inserting (3.9) into (3.8).
Collecting (3.3) and Lemma 5 we then obtain the following error estimate for (2.2).
Theorem 3 Let us suppose u ∈ H 2 (Ω) ∩ H 1 0 (Ω).Then, there exists C > 0, independent of h and ε, such that The following result states that for meshes which are symmetric with respect to their interior nodes, the method converges with a higher order.This result's main interest lies in the diffusion dominated regime, due to the factor ε − 1 2 present in the estimate.The combination of Lipschitz continuity and linearity preservation seems to be novel, and that is why we do detail it now.
and that the mesh is symmetric with respect to its internal nodes.Then, there exists C > 0, independent of h and ε, such that Proof It is enough to bound the consistency error d(u h ; i h u, i h u).We have The first term is bounded as in the proof of Lemma 2. In fact, in that proof, the bound for the second term in (2.10) leads to the following where we have also used the H 1 (Ω)-stability of i h .To bound II we use the linearity preservation and the Lipschitz continuity of d h (•; •, •).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 [8]), i E u satisfies Since the mesh is symmetric with respect to its internal nodes, α E (i E u) = 0.Then, proceeding as in the bound for I we obtain Then, inserting (3.13) and (3.16) into (3.12)we obtain and the result follows by rearranging terms.
4 A link to algebraic flux correction schemes Method (2.2) has been presented having as motivation the study of the effect of adding edge-based diffusion into the equations to impose the discrete maximum principle.Another family of methods that are built with the same purpose is the AFC schemes.This section is devoted to study the relationship between the two approaches, and that is why we now summarise the main building principles of AFC schemes.
The starting point of an algebraic flux-correction scheme is a discretisation of the convectiondiffusion-reaction equation which leads to the linear system where . The first step of these schemes is to identify which parts of the system matrix A are responsible for the violation of the discrete maximum principle.To achieve this, the diffusion matrix D = (d ij ) N i,j=1 is built, where Adding DU both sides of (4.1) we obtain where Ã := A + D. Since the matrix Ã fullfils the hypothesis to guarantee the discrete maximum principle, then the oscillations that appear in a non-stabilised discretisation (4.1) are due to the right-hand side.This is why the right-hand side is now rewritten.Using that the row-sums of D are zero, then The quantities f ij are called fluxes.Then, the AFC schemes are based on introducing limiters α ij (u h ) such that α ij ∈ [0, 1], α ij = α ji , and α ij = 1 if x i and x j are both Dirichlet nodes.Then, after introducing these limiters, the method reads as follows: The most popular limiters in practice are Zalesak's limiters (see, [23,15,16,17], and the recent review [18] for examples).The analysis of this class of methods for a class of limiters that includes the Zalesak one has been carried out recently in [3,2].In particular, in [2] an O(h 1 2 ) convergence rate was proved for the case in which the mesh used satisfies Assumption 1.In the case of meshes that do not satisfy this assumption, then no convergence can be proved.This result is optimal, as the numerical results in [2] show.
Following [2], (4.3) can be written as the following weak problem: Find u h ∈ V h such that u h − u bh ∈ V 0 h , and where the nonlinear form dh (•; •, •) is given by Next, to link this to the method analysed in the last sections, we use the symmetry of D, and of the limiters α ij = α ji , and a simple calculation gives: Then, since where we have adopted the convention that an edge E ∈ E h has endpoints x i and x j , and used that α ij = 1 for edges included in the Dirichlet boundary.Method (2.2) then appears as an algebraic flux-correction scheme, with a different definition of the limiters.Indeed comparing (2.2) with (4.7) we get the equivalent AFC scheme if we choose α ij (u h ) such that The new definition of the limiters made it possible to write some convergence and existence results, also present in [2], in a more precise way, and improve in some of them.In particular, the new limiters make it possible to prove convergence for general meshes, as well as to prove uniqueness of solutions and optimal convergence in the diffusion dominated regime.

Numerical Results
In this section we present three sets of numerical results for bi-dimensional problems.The nonlinear system (2.2) has been solved using the following fixed-point algorithm with damping: Starting with the Galerkin solution u 0 h , then compute a sequence {u k h } defined by where ω ∈ (0, 1) is a damping parameter, and ũk+1 In all our calculations we have used ω = 0.1, and stopped the iterations when the residual R k := ã(u k+1 h ; ψ i ) − (f, ψ i ) Ω i=1,...,dim(V 0 h ) has an euclidean norm smaller than, or equal to, 10 −8 .

Convergence for a smooth solution
We take b = (2, 1), σ = 1, and different values for ε.We have selected the right-hand-side and boundary conditions in such a way that the solution is given by u(x, y) = sin(2πx) sin(2πy).The meshes used were the three-directional mesh (c) and the non-Delaunay mesh (d) in Figure 1.In these calculations we have used γ 0 = 3 and p = 4.
The results in Tables 1-4 match the theoretical results.In particular we observe a first order convergence in the diffusion-dominated regime for the Mesh (c), as predicted by Theorem 4, and a second order convergence in the L 2 norm of the error for both the convection and diffusiondominated regimes.The latter is in accordance with the empirical observations that linearity preservation implies such a convergence.For Mesh (d), which is non-symmetric, and hence the method is no longer linearity preserving, we can observe a first order convergence in both regimes.This convergence is not affected by the non-Delaunay character of the mesh.as Dirichlet condition at entry.We have solved this problem on a uniform refinement of the threedirectional from Mesh (c) in Figure 1.The parameter γ 0 has been set to 1, and the results show no violation of the DMP.The results for this case are depicted in Figure 2. We can observe that the increase in the value of p provides a solution whose inner layer is much sharper than the choice p = 1.We have solved this problem on a criss-cross mesh as shown in Mesh (a) in Figure 1.We have used the parameter γ 0 = 0.75, and, again, no violations of the DMP have been observed.The results are depicted in Figure 3, where we can observe much sharper layers (especially the internal one) when p = 4 is used.

1 2 (
∂Ω), are given data.In this work we adopt the standard notation for Sobolev spaces.In particular, for D ⊂ R d we denote (•, •) D the L 2 (D) (or L 2 (D) d ) inner product, and by • l,D (| • | l,D ) the norm (seminorm) in H l (D) (with the usual convention that H 0 (D) = L 2 (D)).