On Adaptive Eulerian–Lagrangian Method for Linear Convection-Diffusion Problems

On Adaptive Eulerian–Lagrangian Method for Linear Convection-Diffusion Problems


In this paper, we consider the adaptive Eulerian–Lagrangian method (ELM) for linear convection-diffusion problems. Unlike the classical a posteriori error estimations, we estimate the temporal error along the characteristics and derive a new a posteriori error bound for ELM semi-discretization. With the help of this proposed error bound, we are able to show the optimal convergence rate of ELM for solutions with minimal regularity. Furthermore, by combining this error bound with a standard residual-type estimator for the spatial error, we obtain a posteriori error estimators for a fully discrete scheme. We present numerical tests to demonstrate the efficiency and robustness of our adaptive algorithm.


This work was supported in part by NSF Grant DMS-0915153, and DOE Grant DE-SC0006903.


Convection-diffusion equations oftentimes arise in mathematical models for fluid dynamics, environmental modeling, petroleum reservoir simulation, and other applications. Among them, the most challenging case for numerical simulation is the convection-dominated problems (when diffusion effect is very small compared with the convection effect) [43]. Dominated convection phenomena could appear in many real world problems; for example, convective heat transport with large Péclet numbers [27], simulation of oil extraction from underground reservoirs [23], reactive solute transport in porous media [45], etc. The solutions of these problems usually have sharp moving fronts and complex structures; their nearly hyperbolic nature presents serious mathematical and numerical difficulties. Classical numerical methods, developed for diffusion-dominated problems, suffer from spurious oscillations for convection-dominated problems. Many innovative ideas, like Upwinding, Method of Characteristics, and Local Discontinuous Galerkin methods, have been introduced to handle these numerical difficulties efficiently; see, for example, [1] and references therein.

For problems with nearly hyperbolic nature, it is nature to explore the idea of the so-called Method of Characteristics [1]; and, this idea has been combined with different spatial discretizations like finite difference (FD), finite element (FE), and finite volume (FV) methods. Along this line of research, the Semi-Lagrangian method (or, in the finite element context, the Eulerian–Lagrangian method) treats the convection and capacity terms together to carry out the temporal discretization in the Lagrangian coordinate. The Eulerian–Lagrangian method (ELM) gives rise to symmetric discrete linear systems, stabilizes the numerical approximation, and the corresponding diffusion problems are solved on a fixed mesh [18]. This method and its variants have been successfully applied not only on the linear convection-diffusion problem [18], but also the incompressible Naiver-Stokes equations and viscoelastic flow problems; see, for example, [42].

Adaptive mesh refinement (AMR) for partial differential equations (PDEs) has been the object of intense study for more than three decades. AMR techniques have been proved to be successful to deal with multiscale phenomena and to reduce the computation work without losing accuracy when solution is not smooth. In general, the adaptive algorithm for static problems generates graded meshes and iterations in the following form:

In the ESTIMATE procedure, we usually employ some computable local indicators to estimate the local error of the approximate solution we obtain from the SOLVE procedure. These indicators only depend on the datum and/or the approximate solution, and show in which part(s) of the domain the local error is relatively too big or too small in the current mesh. We then MARK these subdomains and REFINE or COARSEN them accordingly.

Local error indicators determine whether the whole adaptive procedure is effective or not. Therefore, a reliable and efficient error indicator is the key to a successful AMR method and a posteriori error analysis is often used to obtain such an error indicator in practice [4]. In the context of finite element methods, the theory of a posteriori error analysis and adaptive algorithms for linear elliptic problem is now rather mature. Convergence and optimality of adaptive methods for linear elliptic problems have been proved as the outcome of a sequence of work [17]; see the recent review by Nochetto, Siebert, and Veeser [38] and references therein. On the other hand, for the nonlinear and time-dependent problems, the theory is still far from satisfactory. A posteriori error analysis for nonlinear evolution problems is even more challenging.

Adaptivity time-stepping is very important for time dependent problems because the practical problems sometimes have singularities or are multiscale in time. Uniform time step size cannot capture these phenomena. There are considerable amount of work in the literature devoted to the development of efficient adaptive algorithms for evolution problems. A posteriori error estimators for linear parabolic problems are studied in [29] and are also derived for nonlinear problems; see [35] for example. There have been also some efforts for extending a posteriori error analysis for the time-dependent Stokes as well as the Navier-Stokes equations [5]. In particular, a posteriori error estimates for convection-diffusion problems have been discussed in [25].

It is nature to employ ARM techniques for convection-dominated problems because of the complex structures of the solutions and evolution of these structures in time. We also notice that spatial mesh adaptivity plays an important role in ELM to reduce numerical oscillations and smearing effect when inexact numerical integrations are employed [28]. Early adoption of adaptive characteristic methods has been seen since late 1980’s [16]. A posteriori error estimates for characteristic-Galerkin method for convection-dominated problems have been proposed: Demokowicz and Oden [16] considered the method of characteristics combined with Petrov-Galerkin finite element method for spatial discretization. Houston and Süli [26] give an a posteriori error estimator in -norm for the linear time-dependent convection-diffusion problem using the duality argument. Chen and Ji [11] give sharp a posteriori error estimations of ELM in -norm for linear and nonlinear convection-diffusion problems, respectively. A related a posteriori error bound can be found in Chen, Nochetto, and Schmidt [13] for the continuous casting problem (convection-dominated and nonlinearly degenerated).

In the previous error estimators mentioned above, the time residual, which measures the difference between the solutions of two successive time steps, is employed as a local indicator for temporal error. On the other hand, it is observed that ELM is essentially transforming a convection-diffusion problem to a parabolic problem along the characteristics. In this paper, we consider a posteriori error estimators for ELM, with focus on temporal error estimators. Motivated by Nochetto, Savaré, and Verdi [37], in which adaptive time-stepping scheme for an abstract evolution equation in Hilbert spaces is analyzed, we can obtain an a posteriori error bound for the temporal error along the characteristics. Combined with space adaptivity, adaptive method has been designed and implemented. Our numerical experiments in Section 5 suggest that the new a posteriori error estimator is effective. Moreover, the numerical results also indicate that the proposed error estimators are very efficient and can take advantage of the fact that ELM allows larger time stepsize.

The outline of this paper is as follows. In Section 2, we introduce a model problem and its Eulerian–Lagrangian discretization. In Section 3, we discuss temporal a posteriori and a priori error analysis for the linear convection-diffusion problem. In Section 4, we present a posteriori error estimation for the full-discretization. To simplify the discussion, we only consider a linear convection-diffusion model problem and restrict ourselves to the standard residual-type estimator for the spatial error, although the technique discussed in this paper can be potentially extended to other problems and spatial error estimators. In Section 5, we consider the implementation of our adaptive algorithm and present some numerical experiments.

Throughout this paper, we will use the following notation. The symbol denotes the space of all square integrable functions and its norm is denoted by . Let be the standard Sobolev space of the scalar function whose weak derivatives up to order are square integrable, and, let and denote the standard Sobolev norm and its corresponding seminorm on , respectively. Furthermore, and denote the norm and the semi-norm restricted to the domain , respectively. We also use the notation for the functions that belong to and their trace vanish on . The dual space of is denoted by . We use the notation to denote the existence of a generic constant , which depends only on , such that .

2A model problem and discretization

In this paper, we consider the following linear convection-diffusion model problem

with the initial condition

and the boundary condition

Here () is bounded polygonal domain. We assume that is divergence free and vanishes on for . We assume that and .

2.1The Eulerian–Lagrangian method

In this section, we briefly recall the construction of the Eulerian–Lagrangian method. As usual, we introduce the characteristics (particle trajectory) , where is the Lagrangian coordinate (original labeling) at time of the particle. It is also referred to the material coordinate. And is the Eulerian coordinate at the current time and referred to as reference coordinate. Then for the given velocity field , the characteristics is defind by the following ordinary differential equation:

In the Eulerian coordinates, be chain rule, the material derivative is defined as following

and then we can rewrite as follow

In order to avoid the deformation of the mesh, ELM uses a fixed spatial mesh at each time-step and traces back along the characteristics. The characteristics at each time interval () has the same original labeling . Assume that and

We discretize the material derivative in using the backward Euler method as follows

Here, , , and . It is easy to see that is a function of in the Lagrangian coordinate. We usually refer to as the temporal semi-discretization scheme.

2.2Feet searching

Our a posteriori error estimate base on the assumption that the characteristics are solved exactly, which preserves the determinant of the Jacobian of the flow. In the numerical analysis, it is difficult to do that and it was apparent that a computational realization of preserving the determinant of the Jacobian for the flow map to be one was crucial. Therefore, we discuss how to integrate the following ordinary differential equation for the computation of the characteristic feet. The numerical scheme we discuss here has second order accuracy and preserves the determinant of the Jacobian of the flow.


The equation (Equation 8) is often called the source-free dynamical systems due to (Equation 9). For such a system, the solution is often called the flow map or phase flow and it has the following property.

We shall begin with introducing some popular scheme to compute (Equation 8) and showing that the scheme is indeed volume-preserving scheme for but not for . We will then introduce some volume preserving scheme to solve (Equation 8) for , which is due to Feng and Shang [24].

In literatures, the following second order numerical scheme for solving (Equation 8) is popular and it seems to first appear in [?].

First, we integrate the equation (Equation 8) using the mid-point rule to obtain :

where .

Second, the right hand side is approximated by a second order accurate extrapolation. Namely,

The following approximation shall also be used :

For notational conveniences, let us denote and Hence we have the following implicit approximations :

To see that the mid-point rule (Equation 13) results in the volume-preserving scheme, let us take the derivative with respect to for both sides of (Equation 13). We then obtain the following :



We solve (Equation 14) for and obtain that

Under the assumption that some appropriate finite element space for the velocity field is used so that the divergence free condition of the velocity field is imposed in the discrete sense, we have

From (Equation 15), we conclude that identically.

The main reason why such an algorithm is popular seems that it has the volume preserving property. On the other hand, it is easy to see that the algorithm may not result in the volume-preserving scheme for . Note that for , under the assumption that , for given as follows,

we have that

Our purpose here is that by reviewing the volume-preserving scheme in three dimension developed by Feng and Shang in [24], we wish to make sure such a volume preserving scheme can be devised in three dimension and confirm our numerical scheme can be implemented. As far as the author is concerned, such a special algorithmic detail has not been implemented in the context of the semi-Lagrangian scheme.

The basic idea of constructing the volume preserving scheme for is based upon the following observation. Following the idea of H. Weyl, we have :


The actual expressions for , , and as follows :

From this, it is easy to see that (Equation 16) holds and by construction and given as follows are divergence free :

Let us now denote by the volume preserving scheme for

with , then the following composition is trivially volume preserving :

Moreover, assuming is of second order accurate, it is easy to see that the composition (Equation 18) is of second order. The above idea on the composition is from Feng and Shang, [24].

2.3Full discretization

Discretizing with suitable spatial discretizations, we obtain fully discrete numerical schemes. In this paper, we focus on the finite element method. First, we define the weak forms of and as usual. We define the bilinear form as follow

We denote the potential by

and its Frechet derivative by , i.e.,

It is easy to see that is convex and satisfies the following inequality

In fact, it is well-known that we have the following identity,

where . Furthermore, by taking the test function as in and applying , we obtain that

The weak forms of (Equation 5) and (Equation 7) can be written as



On a shape-regular triangulation of , we introduce the piecewise continuous linear finite element space such that

Let and denote the mesh and the finite element space at time respectively. Then the fully discrete scheme can be written as: Suppose is known, find such that

where is some suitable approximation of .

3Error analysis for temporal semi-discretization

In this section, we focus on the a posteriori error estimation for temporal semi-discretization (Equation 7). For the sake of simplicity, we assume that in this section. Different from the standard time error indicators for ELM, which measure the solution difference along the time direction, our new error indicator measures the difference along the characteristics. We will establish a posteriori error estimation based on the new time error indicator and show its efficiency. We will also show an optimal priori error bound as a byproduct. Here optimality does not only mean the optimal convergence rate which can also be achieved by classic error analysis, but also mean the optimal regularity requirement which have not been proved by standard a priori error analysis.

3.1A posteriori error analysis

We first show that the solution of satisfies an energy identity along characteristics. In fact, by taking inner product with on both side of and applying the following identity

we immediately obtain that

We can see that the convection diffusion equation preserve the total energy along characteristics in continuous level.

On the other hand, by choosing the test function in and employing , we have an discrete energy inequality: For any integer ,

We note that the equality holds in if there is no temporal discretization error. This motivates the following definition:

and we can view the computable quantity as a measure of the deviation of numerical solution from satisfying the energy conservation . We can use as a time error indicator for adaptive time stepping in ELM.

In order to give an a posteriori error bound for our new time error indicator, we construct the following linear interpolation


Since the original labeling of the characteristics does not depend on time, we have

Substituting back to and then applying , we have

By adding and subtracting , for any , we have

where is the residual (which does not depend on the choice of test function):

We notice that is a convexity functional, i.e.,

Hence, can be bounded by

By choosing in , and in , adding these two equations, and using , we end up with the following inequality

Then the following upper bound of the new time error indicator holds:

Integration in time, we can obtain

Since is divergence free, which implies , we have, after changing of variables, that

Then follows directly from summing up above inequality from to .

3.2An optimal priori error estimation requiring minimal regularity

Traditional a priori error analysis for ELM treats the temporal semi-discretization as a finite difference method and derive the error estimation based on the Taylor expansion (see [18]). As a result, we obtain an optimal convergence rate but the regularity requirement is suboptimal. Taking advantage of the posteriori error estimation of new time error indicator , we can derive an optimal order priori error estimation with minimal regularity requirement on the datum and the solution.

From the definition of the new a posteriori error estimator , and , we have

By the definition of and , we have

Substitute back into , we have

Set and use the elementary inequality , we obtain

Now using the assumption that is divergence free and summing up from to , we obtain the following optimal a priori error estimation with minimal regularity requirement:

4A posteriori error estimation for full discretization

In this section, a posteriori error estimation for the fully-discrete scheme is considered. The characteristic feet is defined by and is computed exactly. Besides the new time error indicator, residual type error estimator is used as a spatial error indicator. Let be the interior residual, i.e.,

and be the jump residual also defined, i.e.,

where is the edge/face shared by and , is the unit normal vector of from to .

As the analysis of time semi-discretization, we define the linear interpolation of as following

Similar with , we have

Now, the following a posteriori error estimation holds

For any and , we have

Let in and in , we have

Inserting , we have

Let and in , we have

Then let and where is the so called Clement interpolation operator. Moreover, let in and add to above inequality, we can obtain that

Then by the definitions , and , we have

Here in last two inequalities, we use the interpolation property of the Clement interpolation, the Cauchy Schwartz inequality, and the Young’s inequality. The norm is defined by . Now, we have

For any , by integrating in time from to , we have

For , by integrating in time from to , we have

Using the assumption that is divergence free, sum the above two inequalities from to , we have

Note that is arbitrary, the above inequality implies the upper bound .

5Numerical experiments

In this section, we first introduce the adaptive algorithm we use in numerical experiments and then present some numerical results.

5.1Adaptive Algorithm

Since the main purpose of this paper is analyze the proposed time error estimator, we start with the algorithm for time step size control. Let be the total tolerance allowed for the error introduced by the time discretization, according to Theorem ?, that is

Therefore, a natural way to satisfy is to adjust the time step size such that

Based on this, now we propose the following algorithm to adjust the time step size and control the error under given tolerance .

The above algorithm is for adjusting the time step size and controlling the error introduce by time discretiziation. When we consider the fully discretization, we need to take mesh adaption into account. According to Theorem ?, the space error indicator gives us the elements where the local error is relatively large and should be refined. Let be the tolerance allowed for the error introduced by the spatial discretization, then the usual stopping criterion for mesh refinement is:

As usual, in order to achieve equal distribution of error, can be replaced by

The main difficulty for time-dependent problems is the mesh coarsening. We need a coarsening error indicator to guide the mesh coarsening procedure. Here, we employ the same coarsening error indicator as the one proposed in [10]. Let be a coarsening mesh from mesh , and and be the corresponding finite element spaces. It is nature that . Furthermore, let and be the numerical solution at time satisfy that

Subtracting these two equalities and taking , we have

Using the elementary equality for the two terms on the left hand side of the above equality respectively, we can obtain the following Galerkin orthogonal relation