On algebraically stabilized schemes for convection-diffusion-reaction problems

An abstract framework is developed that enables the analysis of algebraically stabilized discretizations in a unified way. This framework is applied to a discretization of this kind for convection-diffusion-reaction equations. The definition of this scheme contains a new limiter that improves a standard one in such a way that local and global discrete maximum principles are satisfied on arbitrary simplicial meshes.



There are no comments yet.


page 1

page 2

page 3

page 4


An ADI Scheme for Two-sided Fractional Reaction-Diffusion Equations and Applications to an Epidemic Model

Reaction-diffusion equations are often used in epidemiological models. I...

Analysis of Flux Corrected Transport Schemes for Evolutionary Convection-Diffusion-Reaction Equations

We report in this paper the analysis for the linear and nonlinear versio...

FCNN: Five-point stencil CNN for solving reaction-diffusion equations

In this paper, we propose Five-point stencil CNN (FCNN) containing a fiv...

Efficient discretization and preconditioning of the singularly perturbed Reaction Diffusion problem

We consider the reaction diffusion problem and present efficient ways to...

A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations

In this paper, we propose to improve the stabilized POD-ROM introduced b...

Convergence analysis of the variational operator splitting scheme for a reaction-diffusion system with detailed balance

We present a detailed convergence analysis for an operator splitting sch...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The modeling of physical processes is usually performed on the basis of physical laws, like conservation laws. The derived model is physically consistent if its solutions satisfy the respective laws and, in addition, other important physical properties. Convection–diffusion–reaction equations, which will be considered in this paper, are the result of modeling conservation of scalar quantities, like temperature (energy balance) or concentrations (mass balance). Besides conservation, bounds for the solutions of such equations can be proved (if the data satisfy certain conditions) that reflect physical properties, like non-negativity of concentrations or that the temperature is maximal on the boundary of the body if there are no heat sources and chemical processes inside the body. Such bounds are called maximum principles, e.g., see GT01 . A serious difficulty for computing numerical approximations of solutions of convection–diffusion–reaction equations is that most proposed discretizations do not satisfy the discrete counterpart of the maximum principles, so-called discrete maximum principles (DMPs), and thus they are not physically consistent in this respect. One of the exceptions are algebraically stabilized finite element schemes, e.g., algebraic flux correction (AFC) schemes, where DMPs have been proved rigorously. Methods of this type will be studied in this paper.

The theory developed in this paper is motivated by the numerical solution of the scalar steady-state convection–diffusion–reaction problem


where , , is a bounded domain with a Lipschitz-continuous boundary that is assumed to be polyhedral (if ). Furthermore, the diffusion coefficient is a constant and the convection field , the reaction field , the right-hand side , and the Dirichlet boundary data are given functions satisfying


where is a constant.

In applications, one encounters typically the convection-dominated regime, i.e., it is , where is a characteristic length scale of the problem and denotes the norm in . Then, a characteristic feature of (weak) solutions of (1) is the appearance of layers, which are thin regions where the solution possesses a steep gradient. The thickness of layers is usually (much) below the resolution of the mesh. It is well known that the standard Galerkin finite element method cannot cope with this situation and one has to utilize a so-called stabilized discretization, e.g., see RST08 .

Linear stabilized finite element methods that satisfy DMPs, usually with restrictions to the type of mesh, like the upwind method from BT81 , compute in general very inaccurate results with strongly smeared layers. In order to compute accurate solutions, a nonlinear method has to be applied, typically with parameters that depend on the concrete numerical solution. A nonlinear upwind method was proposed in MH85 and improved in Knobloch06 . In BE05 , a nonlinear edge stabilization method was proposed, see BE02 ; BH04 for related methods, for which a DMP was proved providing that a certain discretization parameter is chosen to be sufficiently large and the mesh is of a certain type. However, already the numerical results presented in BE05 show spurious oscillations. Our own experience from JK08 is that the nonlinear problems for sufficiently large parameters often cannot be solved numerically.

A class of methods that has been developed intensively in recent years is the class of algebraically stabilized schemes, e.g., see BB17 ; BJK17 ; GNPY14 ; Kuzmin06 ; Kuzmin07 ; Kuzmin09 ; Kuzmin12 ; KuzminMoeller05 ; KS17 ; KT04 ; LKSM17 . The origins of this approach can be tracked back to BorisBook73 ; Zalesak79 . In these schemes, the stabilization is performed on the basis of the algebraic system of equations obtained with the Galerkin finite element method. Then, so-called limiters are computed, which maintain the conservation property and which restrict the stabilized discretization mainly to a vicinity of layers to ensure the satisfaction of DMPs without compromising the accuracy. There are several limiters proposed in the literature, like the so-called Kuzmin Kuzmin07 , BJK BJK17 , or BBK BBK17 limiters. Both, the Kuzmin and the BBK limiters were utilized in BBK17a for defining a scheme that blends a standard linear stabilized scheme in smooth regions and a nonlinear stabilized method in a vicinity of layers.

An advantage of algebraically stabilized schemes is that they satisfy a DMP by construction, often under some assumptions on the mesh, and they usually provide sharp approximations of layers, cf. the numerical results in, e.g., ACF+11 ; GKT12 ; JS08 ; Kuzmin12 . In numerical studies presented in JJ19 , it turned out that the results with the BJK limiter were usually more accurate than with the Kuzmin limiter, if the nonlinear problems for the BJK limiter could be solved. However, solving these problems was often not possible for strongly convection-dominated cases. Numerical studies in BJK18 show that using the Kuzmin limiter leads to solutions with sharper layers compared with the solutions obtained with the BBK limiter. As a consequence of these experiences, it seems to be advisable from the point of view of applications to use algebraically stabilized schemes on the basis of the Kuzmin limiter. The AFC scheme with the Kuzmin limiter was analyzed in BJK16

, thereby proving the existence of a solution, the satisfaction of a local DMP, and an error estimate. The local DMP requires lumping the reaction term and using certain types of meshes, e.g., Delaunay meshes in two dimensions, analogously as for the methods from

BBK17 ; BBK17a .

The conservation and stability properties of algebraically stabilized schemes are given if the added stabilization is a symmetric term. For many schemes, this term consists of two factors, an artificial diffusion matrix and the matrix of the limiters, and usually the methods are constructed in such a way that both factors are symmetric. Only recently, motivated by BB17 , a more general approach where only the product is symmetric but not the individual factors was considered in Kno21 .

The first main goal of this paper is the development of an abstract framework that allows to analyze algebraically stabilized discretizations in a unified way. Although our main interest is the numerical solution of problem (1), many considerations will be more general and then problem (1) and its discretizations will only serve as a motivation for our assumptions. Hence, this framework covers a larger class of algebraically stabilized discretizations than the available analysis.

The second main goal consists in proposing and analyzing a modification of the Kuzmin limiter such that, if applied in the framework of the algebraic stabilization of Kno21 , the positive features of the AFC method with the Kuzmin limiter are preserved on meshes where it works well and, in addition, local and global DMPs can be proved on arbitrary simplicial meshes. In particular, our intention was to preserve the upwind character of the AFC method with the Kuzmin limiter. There are already proposals in this direction in the framework of AFC methods. In Kno17 , the Kuzmin limiter is replaced in cases where it does not lead to the validity of the local DMP in a somewhat ad hoc way by a value that introduces more artificial diffusion. The satisfaction of the local DMP on arbitrary simplicial meshes could be proved for this approach. Whether or not the assumption for the existence of a solution of the nonlinear problem is satisfied with this limiter is not discussed. A combination of the Kuzmin and the BJK limiters to obtain a limiter of upwind type for which the AFC scheme satisfies a local DMP on arbitrary simplicial meshes and is linearity preserving is proposed in Kno19 . The definition of this limiter is closer to the BJK than to the Kuzmin limiter. As already mentioned, in Kno21 , a new algebraically stabilized method was proposed that does not require the symmetry of the limiter. Initial numerical results for a nonsymmetric modification of the Kuzmin limiter are presented in Kno21 , but a numerical analysis is missing. The abstract framework mentioned in the previous paragraph covers in particular the method from Kno21 .

In the present paper, the limiter from Kno21 is written in a simpler form, without using internodal fluxes typical for AFC methods. Moreover, a novel modification is performed that improves the accuracy in some computations using non-Delaunay meshes. Of course, this modification is performed in such a way that the resulting method still fits in the abstract analytic framework. The definition of the new method does not contain any ambiguity, in contrast to the AFC method with Kuzmin limiter, which is not uniquely defined in some cases (cf. Remark 8 in BJK16 ). A further advantage of the considered approach is that, in contrast to the AFC method with Kuzmin limiter, lumping of the reaction term is no longer necessary for the satisfaction of the DMP, which enables to obtain sharper layers as we will demonstrate by numerical results.

This paper is organized as follows. Sect. 2 introduces the basic discretization of (1) and its algebraic form. An abstract framework for an algebraic stabilization is presented in Sect. 3. The following section studies the solvability and the satisfaction of local and global DMPs for the abstract algebraic stabilization and Sect. 5 provides an error analysis. In Sect. 6, the AFC scheme with Kuzmin limiter as an example of algebraic stabilization from Sect. 3 is presented, its properties are discussed for the discretizations from Sect. 2 and the definition of the limiter is reformulated. The reformulation is utilized in Sect. 7 for proposing a new limiter such that the resulting algebraically stabilized scheme is of upwind type and satisfies DMPs on arbitrary simplicial meshes. Sect. 8 presents numerical examples which show that the algebraically stabilized scheme with the new limiter in fact cures the deficiencies of the AFC scheme with Kuzmin limiter.

2 The convection–diffusion–reaction problem and its finite element discretization

The weak solution of the convection–diffusion–reaction problem (1) is a function satisfying the boundary condition on and the variational equation



As usual, denotes the inner product in or . It is well known that the weak solution of (1) exists and is unique (cf. Evans ).

An important property of problem (1) is that, for in , its solutions satisfy the maximum principle. The classical maximum principle (cf. Evans ) states the following: if solves (1) and the functions and are bounded in , then, for any set , one has the implications


where , . If, in addition, in , then


Analogous statements also hold for the weak solutions, cf. GT01 .

To define a finite element discretization of problem (1), we consider a simplicial triangulation of which is assumed to belong to a regular family of triangulations in the sense of Ciarlet . Furthermore, we introduce finite element spaces

consisting of continuous piecewise linear functions. The vertices of the triangulation will be denoted by and we assume that and . Then the usual basis functions of are defined by the conditions , , where is the Kronecker symbol. Obviously, the functions form a basis of . Any function can be written in a unique way in the form


and hence it can be identified with the coefficient vector


Now an approximate solution of problem (1) can be introduced as the solution of the following finite-dimensional problem:

Find such that , , and


where is a bilinear form approximating the bilinear form . In particular, one can use . Another possibility is to set


for any and , i.e., to consider a lumping of the reaction term in . This may help to satisfy the DMP for problem (9), cf. Sect. 6. We assume that is elliptic on the space , i.e., there is a constant such that


where is a norm on the space but generally only a seminorm on the space . This guarantees that the discrete problem (9) has a unique solution. In view of (2), the ellipticity condition (11) holds for both and defined by (10) with and


We denote


Then is a solution of the finite-dimensional problem (9) if and only if the coefficient vector corresponding to satisfies the algebraic problem

As discussed in the introduction, the above discretizations are not appropriate in the convection-dominated regime and a stabilization has to be applied. In the next sections, algebraic stabilization techniques will be studied. As already mentioned, a general framework will be presented and the numerical solution of convection–diffusion–reaction equations serves just as a motivation for the assumptions.

3 An abstract framework

In this section we assume that we are given a system of linear algebraic equations of the form


(with ) corresponding to a discretization of a linear boundary value problem for which the maximum principle holds. An example is the algebraic problem derived in the preceding section.

We assume that the row sums of the system matrix are nonnegative, i.e.,


and that the submatrix is positive definite, i.e.,


For the discretizations from the previous section, the latter property follows from (11), whereas (18) is a consequence of the nonnegativity of and the fact that .

Since the algebraic problem (16), (17) is assumed to approximate a problem satisfying the maximum principle, it is natural to require that an analog of this property also holds in the discrete case, at least locally. Then an important physical property of the original problem will be preserved and spurious oscillations of the approximate solution will be excluded. To formulate a local DMP, we have to specify a neighborhood

of any (i.e., of any interior vertex if the geometric interpretation from the previous section is considered). For example, one can set


Then, under the assumptions (18) and (19), the solution of (16), (17) satisfies the local DMP


(with any ) if and only if (cf. (BJK16, , Lemma 21))


Moreover, the stronger local DMP


holds (again with any ) if and only if the conditions (22) and


are satisfied (cf. (BJK16, , Lemma 22)). For the discretizations from the previous section, (24) holds if in , which is a condition used for proving (6) and (7), i.e., a counterpart of (23).

In many cases, the condition (22) is violated (like for the discretizations from the previous section in the convection-dominated regime) and hence the local DMPs (21) and (23) do not hold. To enforce the DMP, one can add a sufficient amount of artificial diffusion to (16), e.g., in the following way. First, the system matrix is extended to a matrix , typically using the matrix corresponding to the underlying discretization in the case when homogeneous natural boundary conditions are used instead of the Dirichlet ones (i.e., using (13) if the setting of the previous section is considered). Then one can define a symmetric artificial diffusion matrix possessing the entries


The matrix has zero row and column sums and is positive semidefinite (cf. (BJK16, , Lemma 1)), the matrix has nonpositive off-diagonal entries by construction and the submatrix is positive definite. Consequently, the stabilized algebraic problem


is uniquely solvable and its solution satisfies the local DMP (21) with


If the condition (24) holds, then the solution of (26), (27) also satisfies the stronger local DMP (23), again with defined by (28). Moreover, if the above stabilization is applied to the discretizations from the previous section, then, for weakly acute triangulations, the approximate solutions converge to the solution of (1), see BJK16 .

However, the amount of artificial diffusion added in (26) is usually too large and leads to an excessive smearing of layers if it is applied to stabilize discretizations of (1) in the convection-dominated regime. To suppress the smearing, the artificial diffusion should be added mainly in regions where the solution changes abruptly and hence it should depend on the unknown approximate solution . This motivates us to introduce a general artificial diffusion matrix having analogous properties as the matrix , i.e., for any , we assume that


Like above, we introduce local index sets such that


and, for any ,


Let us mention that if the algebraic problem (16), (17) corresponds to a finite element discretization based on piecewise linear functions as in the preceding section, one can usually use index sets


, where are the vertices of the underlying simplicial triangulation, numbered as in the preceding section.

Now, we consider the nonlinear algebraic problem


Note that, in view of (31) and (33), system (35) can be written in the form


In view of (29) and (30), one obtains the important property (cf. (BJK16, , Lemma 1))


Thus, due to (31), the matrix is positive semidefinite for any .

4 Analysis of the abstract nonlinear algebraic problem

The aim of this section is to investigate the solvability and the validity of the DMP for the nonlinear algebraic problem (35), (36). These investigations will generalize the results obtained in BJK16 ; Kno17 ; BJK18 .

To prove the solvability of the system (35), (36), we make the following assumption.

Assumption (A1): For any and any , the function is a continuous function of and, for any and any , the function is a bounded function of .

Theorem 1.

Let (19) and (29)–(31) hold and let Assumption (A1) be satisfied. Then there exists a solution of the nonlinear problem (35), (36).


The proof follows the lines of the proof of Theorem 3 in BJK16 . We denote by the elements of the space and, if with occurs, we assume that . To any , we assign . Let us define the operator by

Then is a solution of the nonlinear problem (35), (36) if and only if . The operator is continuous and, in view of (19) and (38), there exist constants , such that (cf. (BJK16, , Theorem 3) for details)

where is the usual inner product in and the corresponding (Euclidean) norm. Then, for any satisfying , one has and hence it follows from Brouwer’s fixed-point theorem (see (Temam77, , p. 164, Lemma 1.4)) that there exists such that . ∎

Remark 1.

For proving the solvability of (35), (36), it would be sufficient to assume that the functions are continuous. However, since should depend on local variations of with respect to , the assumed continuity of is more useful. The functions themselves are often not continuous, cf. Remark 7.

Remark 2.

The solution of (35), (36) is unique if is Lipschitz–continuous with a sufficiently small constant. As pointed out in Lohman19 , this condition can be further refined by introducing a positive semidefinite matrix , e.g., the one defined in (25), and investigating the Lipschitz continuity of . Since, in view of (19), there is such that

( is again the Euclidean norm on ), the smallness assumption on the Lipschitz constant can be expressed by the inequality


Then, if are two solutions of (35), (36), one has

and (39) leads to a contradiction. Nevertheless, the inequality (39) is often not satisfied and then the uniqueness of the nonlinear problem (35), (36) is open.

Now let us investigate the validity of DMPs for problem (35), (36). To this end, one has to relate the properties of the artificial diffusion matrix to the matrix . This can be done in various ways and we shall use the following assumption that generalizes the one used in Kno17 .

Assumption (A2): Consider any and any . If is a strict local extremum of with respect to from (32), (33), i.e.,


Remark 3.

In contrast to linear problems, it is only assumed that off-diagonal entries of the matrix are nonpositive in rows corresponding to indices where strict local extrema of appear. If does not depend on , then Assumption (A2) implies that the first rows of have nonpositive off-diagonal entries, which is a necessary and sufficient condition for the validity of the local DMP under our assumptions on and .

Theorem 2.

Let (18), (19), and (29)–(33) hold and let Assumption (A2) be satisfied. Then any solution of (35) satisfies the local DMP (21) for all . If condition (24) holds, then the stronger local DMP (23) is also valid.


The proof is basically the same as in Kno17 . Since it is short, we repeat it for completeness. Let satisfy (35). Consider any and let . Denoting , it follows from (37) that


If , we want to prove the first implication in (21) for which it suffices to consider since otherwise the implication trivially holds. If , an arbitrary sign of is considered. Let us assume that for all . Then Assumption (A2) implies that each term of the sum in (40) is nonnegative. If , then there is such that since (see (19)). This together with (30) implies that the sum in (40) is positive. If , then . Thus, in both cases, the left-hand side of (40) is positive, which is a contradiction. Therefore, there is such that , which proves the first implication in (23) and hence also in (21). The statements for follow in an analogous way. ∎

Our next aim will be to show that, under the above assumptions, also a global DMP is satisfied. First we prove the following general form of the DMP, which generalizes a result proved in BJK18 .

Theorem 3.

Let (18), (19), and (29)–(33) hold and let Assumptions (A1) and (A2) be satisfied. Consider any nonempty set and denote


Assume that . Then any solution of (35) satisfies the DMP


If, in addition,




The proof is based on the technique used in (Kno10, , Theorems 5.1 and 5.2). Let satisfy (35) and let for all . We denote

Then, according to (31)–(33), (18), (38), (19) and (35), one has



It suffices to consider the case since otherwise the validity of (42) and (45) is obvious. First, let us show that


Let and . For any , define the vector with and for . Then is a strict local maximum of with respect to and hence, in view of Assumption (A2),

Since for , Assumption (A1) implies that

As , it follows that . For , one has , which completes the proof of (50).

Now we want to prove that the relations (47)–(50) imply (42) and (45). If (44) does not hold, it suffices to consider since otherwise (42) trivially holds. Let us assume that (45) does not hold, which implies that . We shall prove that then


Assume that (51) is not satisfied. Then, applying (47) and (50), one derives for any

which gives

Thus, the matrix is singular, which contradicts (48). Therefore, (51) holds and hence, denoting one obtains using (49) and (50)


If (44) holds, then, in view of (47), the right-hand side of (52) equals . Hence, , which is a contradiction to the definition of . If (44) does not hold, then it is assumed that and hence, in view of (50), the inequality (52) implies that . Thus, in view of (47), the right-hand side of (52) is bounded by , which again implies that . Therefore (45) and hence also (42) hold.

The implications (43) and (46) can be proved analogously. ∎

Remark 4.

Note that may contain also indices from the set . The assumption is always satisfied if (44) holds since otherwise, due to (47), the matrix would be singular, which is not possible in view of (48). If satisfies (35) with for all and , then it was shown in the above proof that , which again implies that . The same holds if satisfies (35) with for all and .

Setting in Theorem 3, one obtains the following global DMP.

Corollary 1.

Let (18), (19), and (29)–(33) hold and let Assumptions (A1) and (A2) be satisfied. Then any solution of (35) satisfies the global DMP


If, in addition, the condition (24) holds, then


Finally, let us return to the convection–diffusion–reaction problem (1) and assume that the algebraic problem (16), (17) is defined by (13)–(15) with given by (3) or (10). Recall that a vector can be identified with a function via (8). Then, for index sets defined by (34), Theorem 3 implies that finite element functions corresponding to obeying to (35) satisfy an analog of the continuous maximum principle (4)–(7).

Theorem 4.

Let the assumptions stated in Sect. 1 be satisfied and let the algebraic problem (16), (17) be defined by (13)–(15) with given by (3) or (10). Let the index sets be given by (34). Consider a matrix depending on and satisfying (29)–(31), (33), and Assumptions (A1) and (A2). Consider any nonempty set and define

Let be a solution of (35) and let be the corresponding finite element function given by (8). Then one has the DMP


If, in addition, in , then



where denotes the interior of . Since for any and is piecewise linear, one has


If , then for any and (61) immediately implies the validity of the right-hand sides in the implications (57)–(60). Thus, assume that . Let and be defined by (41). Then, in view of the definition of , one has and . If in , then for any and hence

according to (42). If , then and hence

Consequently, (57) holds due to (61). The implications (58)–(60) follow analogously. Note that if in , then (44) holds since . ∎

Remark 5.

It might be surprising that the local DMP proved in Theorem 2 was not employed for proving the global DMP and instead a much more complicated proof was considered in Theorem 3. However, the global DMP cannot be obtained as a consequence of the local DMPs as the following example shows. Let be values at the vertices of the triangulation depicted in Fig. 1

Figure 1: Local DMP does not imply a global DMP

numbered as in Sect. 2. Let (values at the black interior vertices) and (values at the white boundary vertices). Let the index sets be given by (34). Then the local DMP

is satisfied but the corresponding global DMP (the right-hand sides of the implications (53) and (55) with and ) does not hold.

5 An error estimate

In the previous section, we analyzed the nonlinear algebraic problem (35), (36) on its own, without relating it to some discretization (except for Theorem 4). If the algebraic problem originates from a discretization of the convection–diffusion–reaction problem (1), then a natural question is how well its solution approximates the solution of (1). This question will be briefly addressed in this section.

Let us assume that the algebraic problem (16), (17) corresponds to the variational problem (9) satisfying (11), i.e., it is defined by (13)–(15). Let correspond to the solution of the nonlinear algebraic problem (35), (36) via (8). Our aim is to estimate the error . To this end, it is of advantage to write the nonlinear algebraic problem in a variational form. We denote