© Hindawi Publishing Corp. THE hp VERSION OF EULERIAN-LAGRANGIAN MIXED DISCONTINUOUS FINITE ELEMENT METHODS FOR ADVECTION-DIFFUSION PROBLEMS

We study the hp version of three families of Eulerian-Lagrangian mixed discontinuous finite element (MDFE) methods for the numerical solution of advectiondiffusion problems. These methods are based on a space-time mixed formulation of the advection-diffusion problems. In space, they use discontinuous finite elements, and in time they approximately follow the Lagrangian flow paths (i.e., the hyperbolic part of the problems). Boundary conditions are incorporated in a natural and mass conservative manner. In fact, these methods are locally conservative. The analysis of this paper focuses on advection-diffusion problems in one space dimension. Error estimates are explicitly obtained in the grid size h, the polynomial degree p, and the solution regularity; arbitrary space grids and polynomial degree are allowed. These estimates are asymptotically optimal in both h and p for some of these methods. Numerical results to show convergence rates in h and p of the Eulerian-Lagrangian MDFE methods are presented. They are in a good agreement with the theory.


Introduction.
The development of discontinuous Galerkin finite element methods for the numerical solution of partial differential equations dates back to the early 1970s.They were first introduced in [34] for the neutron transport equation (a first-order linear differential equation).Their extensions, stability, and convergence analyses were later carried out by many researchers [19,26,27,28,29,35].
Discontinuous Galerkin methods have also been developed for second-order elliptic problems.A discontinuous method for these problems, the so-called global element method, was introduced in [21].The matrix arising from the space discretization of diffusion terms by this method is indefinite.Also, its stability and convergence properties are not available.The interior penalty method in [2,22,37] used the bilinear form of the global element method augmented with an interior penalty term.The drawback of this penalty method is that its stability and convergence depend on penalty parameters.The penalty idea was introduced in [31] by replacing boundary multipliers by normal fluxes for the finite element solution of elliptic problems.Similar approaches were utilized in [4,5] for fourth-order differential problems.
A slightly different discontinuous Galerkin method has recently been introduced in [32].Their method is different from the global element method in that the sign in front of the so-called symmetric term was altered.The disadvantage of this method is that the bilinear form is nonsymmetric even for a symmetric differential problem.Also, only suboptimal error estimates in the L 2 -norm for numerical solutions can be obtained.This method and its penalized versions have been studied in [36].For the relationship between all the existing discontinuous Galerkin methods, we refer to [13].
There have been a lot of applications of the above two families of discontinuous Galerkin methods and their stabilized versions to practical problems.Among them are the applications to fluid flows in porous media [11,17], to the Euler equations [7], to nonlinear acoustic waves [30], and to semiconductor device modeling [15,16], for example.
Discontinuous finite element methods in mixed form were utilized in [15,16] and also developed in [6].While the lowest-order Raviart-Thomas space [33] was used in [15,16], Lagrange multipliers on interior boundaries were exploited to relax the continuity of normal derivatives of functions in the vector space, so all finite elements are discontinuous.These methods in mixed form have been recently studied in [20] for the convection-diffusion equations and in [8] for the Laplace equation.
The objective of this paper is to study the hp version of the three families of Eulerian-Lagrangian mixed discontinuous finite element (MDFE) methods and their equivalent Galerkin versions for the numerical solution of advectiondiffusion problems.These methods are based on a space-time mixed formulation of the advection-diffusion problems.In space, they use discontinuous finite elements, and in time, they approximately follow the Lagrangian flow paths (i.e., the hyperbolic part of the problems).Boundary conditions are incorporated in a natural and mass conservative manner.In fact, these methods are locally conservative.A characteristic treatment to handle the advection term in time is similar to those exploited in [1,10]; it is known that this technique is suitable for advection-dominated problems.Also, the stationary version of the MDFE methods (i.e., as the MDFE methods applied to elliptic problems) has been studied in [12].The first family (see (3.1)) was developed in [13,14] and extends the mixed discontinuous methods originally proposed in [6,15,16,20], the second family (see (4.2)) is a generalization of the local discontinuous Galerkin method in [9] (which is a particular case of the local discontinuous method introduced in [20]), and the third family (see (5.1)) is a modification of the second family.
The analysis of this paper focuses on advection-diffusion problems in one space dimension.Error estimates are explicitly obtained in the grid size h, the polynomial degree p, and the solution regularity; arbitrary space grids and polynomial degree are allowed.For the first family, a suboptimal rate of convergence in the L 2 -norm for both the solution and its derivative in space is obtained for odd p and an optimal rate is derived for even p, see Theorems 3.2 and 3.3.For the second and third families, the optimal rate of convergence can be achieved for certain boundary conditions or by choosing appropriately a parameter at the boundary of the domain, see Theorems 4.2 and 5.2.Numerical results to show convergence rates of the Eulerian-Lagrangian MDFE methods in p and h are presented.They are in a good agreement with the obtained theoretical results.
As observed in [13,14], when discontinuous finite element methods are defined in a mixed form, they not only preserve good features of these methods, but also have some advantages over the classical Galerkin form such as they are more stable and more accurate (in space).While an auxiliary variable is introduced in the mixed formulation, the MDFE methods can be reduced to the Galerkin version in one of the two variables because of their local property and can be implemented (if desired) in nonmixed form [13,14].
The rest of the paper is outlined as follows.In Section 2, we describe the continuous problem and its mixed weak formulation.The first, second, and third families of the Eulerian-Lagrangian MDFE methods are defined and discussed in Sections 3, 4, and 5, respectively.The proof of the convergence results is presented in Section 6.Some generalizations of these methods are mentioned in Section 7. Numerical experiments for showing convergence rates in both p and h are presented in Section 8.

The continuous problem. We consider the advection-diffusion equation for u on a bounded interval I = (c, d) with boundary ∂I
where , and u 0 (x) ∈ L 2 (I) are given functions (the standard Sobolev spaces H k (I) = W k,2 (I) with the usual norms are used in this paper), and ν I is the outer unit normal to ∂I (while we call it a normal, it is, in fact, a scalar in one dimension; i.e., it is 1 or −1).To define the mixed weak formulation, we rewrite this equation as follows: (2.2) 1.An illustration of the unit normal ν.
Namely, an auxiliary variable σ is introduced, which has a physical meaning in applications such as the electric field in semiconductor modeling [15,16] or the velocity field in petroleum simulation [23,25], as mentioned in the introduction.In the next four sections, we study the case where the following assumption holds: with a * and a * being constants.The situation without this assumption will be addressed in Section 7. Also, to demonstrate the idea of the proposed approximation scheme, we consider the case where φ and b are constants; the nonconstant case will be addressed in Section 7 as well.
For h > 0, let I h be a partition of I into subintervals h denote the set of all interior nodes of I h , Ᏽ b h the set of the end points on ∂I, and We tacitly assume that Ᏽ o h = ∅.For l ≥ 0, define (2.4) With each x i ∈ Ᏽ h , we associate a unit normal ν.For x i ∈ Ᏽ b h , ν = ν x i is just the outer unit normal to ∂I; that is, for the left end, ν = −1 and for the right end, ν = 1.For x i ∈ Ᏽ o h , it is chosen pointing to the element with lower index (see Figure 2.1); that is, ν = −1 at all interior points.This is just for notational convenience; other choices are possible.
For v ∈ H l (I h ) with l > 1/2, we define its average and jump at x i ∈ Ᏽ o h as follows (see Figure 2.2): (2.5) , the following convention is used (from inside I): For each positive integer ᏺ, let 0  (2.7) The origin of our approximation scheme can be seen by considering (2.2) in a space-time framework.For any x ∈ I and two times 0 ≤ t n−1 < t n ≤ T , the hyperbolic part of problem (2.1), φ∂u/∂t + b∂u/∂x, defines the characteristic xn (x, t) along the interstitial velocity ϕ = b/φ (it emanates backward from x at t n , see Figure 2.3): Similarly, we can modify the characteristic at the left boundary For some element I i ∈ I h , let Ǐi (t) indicate the traceback of I i to time t, t ∈ J n : Ǐi (t) = x ∈ I : x = xn (y, t) for some y ∈ I i .
(2.10) Also, let I n i be the space-time region that follows the characteristics: For an element on either the boundary , an analogous definition can be given (see Figure 2.3).Let (•, •) S denote the L 2 (S) inner product (we omit S if S = I).Multiplying the first equation of (2.2) by v ∈ H 1 (I n i ), integrating the resulting equation on I n i , and using Green's formula, we see that (2.12) where we use the fact that (b, φ) • ν = 0 on the space-time sides (2.13) Note that v(x, t n−1,+ ) = vn−1,+ (x) follows the characteristics forward from t n−1 to t n to become v n (x).If we use this type of test functions in (2.12), the last term in the right-hand side of (2.12) becomes zero.With this, summing the resulting equation over I i ∈ I h , applying the boundary condition in (2.2), and assuming the continuity of σ in I, we see that where with each characteristic a unit normal ν is associated, the characteristic xn (x j ,t) emanates back from the boundary 3), and v is extended by zero outside Ī.Note that we can formally write ) and using the Dirichlet boundary condition in (2.2), (2.16) becomes (2.17) The jump terms in the left-hand side of (2.14) are called the consistent terms (which come from the Green formula), while the corresponding terms in (2.17) are termed the symmetric terms (which are added).Finally, notice that the strips I n i are oriented along the Lagrangian flow paths (i.e., the characteristic paths), and a fixed partition is utilized at the time level t n .These features suggest the name Eulerian-Lagrangian method [1,10].

The first mixed discontinuous method
3.1.The definition.Let V h ×W h be the finite element spaces for the approximation of σ and u, respectively.They are finite dimensional and defined locally on each element

Neither continuity constraint nor boundary values are imposed on
Based on (2.14) and (2.17), the first Eulerian-Lagrangian mixed discontinuous method for (2.1) is: find where v is extended by zero outside Ī.The initial approximation is given by where ᏼ h denotes the We define the bilinear forms (3.4) Then system (3.1) is of the form It can be seen [14] that the stiffness matrix arising from this system is positive definite for any V h and W h .Also, this system is symmetric after changing a sign in either of the two equations.But this would alter the positive definiteness.
As seen in the subsequent analysis, if (3.1) is written in a nonmixed form (or the standard Galerkin version), the positive definiteness and symmetry can be preserved simultaneously.Thus, in terms of implementation, it is desirable to write (3.1) in a nonmixed form.However, as pointed out in the introduction, we emphasize that the mixed formulation naturally stabilizes the discontinuous finite element method, see Section 3.4.That is why we start with the mixed discontinuous method.For the advantages of the present mixed discontinuous methods over the classical mixed finite element methods, see [13,14].Let I i ∈ I h be any element.Take v = 1 on I i and zero elsewhere in the first equation of (3.1) to see that This equation expresses local conservation of mass along the characteristics if the flux across ∂I i is defined as {σ h • ν}.This is one of the advantages of the Eulerian-Lagrangian approach over the classical modified method of characteristics [14,24], which has an inherent difficulty in conserving mass locally.

Existence and uniqueness.
Stability in terms of data for (3.1) has been considered in [14].For completeness, we state existence and uniqueness of the solution to (3.1).Below, C (with or without a subscript) indicates a generic constant independent of h and p (the polynomial degree), which may take on different values in different occurrences.

Proposition 3.1. System (3.1) has a unique solution under (2.3).
Proof.Since (3.1) is a finite system, it suffices to show the uniqueness of the solution.Setting f = g D = g N = u 0 = 0, it follows from (3.1) that We complete the proof by an induction argument.Note that u 0 h = 0 by assumption and by (3.2).Let u n−1 h = 0.The addition of the two equations in (3.7) with v = u n h and τ = σ n h yields which implies, by (2.3), that σ n h = u n h = 0. Thus, the desired result follows.

Convergence. From now on, let
where P p i (I i ) is the set of polynomials of degree at most p i on I i .Note that p i can be different on different elements.
For v ∈ H l (I i ), l ≥ 0, there exists v h ∈ P p i (I i ) and a positive constant C dependent on l but independent of v, h, and p i such that [3] (3.10) We now state error estimates which will be shown in Section 6. Set Note that the estimate in Theorem 3.2 gives a suboptimal order in h of convergence, but it is sharp for (3.1) in the general case, as indicated by numerical experiments.In the case where p i for all I i ∈ I h is even, we can prove an optimal order in h of convergence for (3.1).Theorem 3.3.Under (2.3), if (σ , u) and (σ h ,u h ) are the respective solutions to (2.2) and (3.1) and if p i is even for all I i ∈ I h , then, for ∆t sufficiently small, The proof of this theorem will be carried out in Section 6 as well.Notice that the error estimate in Theorem 3.3 is optimal in h.

Implementation.
While method (3.1) is in mixed form, it can be implemented (if desired) in a nonmixed form [14].We introduce the coefficientdependent projections P n h : for w ∈ L 2 (I), and R n h : for v ∈ H 1 (I h ).Using (3.15) and (3.16), (3.1) can be rewritten as follows [14]: with σ h given by To see the relationship of (3.17) with traditional discontinuous Galerkin finite element methods, we consider the case where a is piecewise constant.In this case, (3.17) becomes: find Note that without the terms involving R h , (3.19) (as applied to elliptic problems) is just the global element method introduced in [21].With a positive sign in front of the fourth term in the left-hand side of (3.19) (and without the R h terms), it is the method recently introduced in [32].The R h term in the lefthand sides of (3.17) and (3.19) naturally comes from the mixed formulation, which stabilizes the classical method and leads to higher convergence rates (in space).For more information on the relationship between the present MDFE method and other earlier methods, we refer to [13,14].
Although R h appears, (3.19) can be evaluated virtually in almost the same amount of work as in the evaluation of the global element method.This is due to the definition of R h in (3.16),where the matrix associated with the left-hand side can be diagonal if the basis functions of V h are appropriately chosen.Also, (3.16) is totally local.The stiffness matrices arising from (3.17) and (3.19) are symmetric and positive definite.Numerical experiments will be given in Section 8.

The second mixed discontinuous method
4.1.The definition.For v ∈ H 1 (I h ), define the one-sided limits at the nodes At the endpoint c, for τ ∈ V h , we use the convention τ(c + ) = τ(c − ); a similar meaning can be given at x = d.The second Eulerian-Lagrangian mixed discontinuous method for (2.1) is: find The second method differs from the first one in that the averaged quantities in (3.1) are replaced by the right-hand side limits.This method is an extension of the local discontinuous Galerkin method proposed in [9] for advectiondiffusion problems.

Convergence.
The finite element spaces V h and W h are defined as in (3.9).The convergence result in Theorem 3.2 holds for (4.2) as well.As shown in [12] for elliptic problems, the optimality of convergence in h and p depends on the type of boundary conditions.We consider the case where c ∈ Γ D and d ∈ Γ N .In this case, we have the following optimal convergence result. 3) The proof of this theorem will be carried out in Section 6.The error estimate in Theorem 4.2 is optimal in both h and p.For other cases of boundary conditions, the error estimate in Theorem 3.2 is sharp for (4.2).However, with an addition of appropriate boundary terms at x = c, d to (4.2), the resulting scheme can generate optimal error estimates.We will not pursue this; for more information in this direction, we refer to [12] for the treatment of elliptic problems.

Implementation. System (4.2) can be implemented (if desired) in the same fashion as in (3.1). The operator P n
for v ∈ H 1 (I h ).Then can be reduced to: find u h : {t 1 ,...,t ᏺ } → W h such that with σ h given as in (3.18).Again, in the case in which a is piecewise constant, (4.5) becomes: find The stiffness matrices arising from (4.5) and (4.6) are symmetric and positive definite.Numerical results for (4.2) will be given in Section 8 as well.

The third mixed discontinuous method
5.1.The definition.The third method is analogous to the second one.Recall that the averaged quantities in (3.1) are replaced by the right-hand side limits in (4.2).In the third method, they are replaced by the left-hand side limits.That is, the third Eulerian-Lagrangian mixed discontinuous method for (2.1) is: find (5.1)

Existence and uniqueness.
The analysis for (5.1) can be done in the same manner as for (4.2).We just state the corresponding results.Proposition 5.1.Under (2.3), ( 5.1) has a unique solution.

Convergence.
The finite element spaces V h and W h are defined as in (3.9).Again, the convergence result in Theorem 3.2 holds for (5.1).Similar to the second method, the optimality of convergence in h and p depends on the type of boundary conditions.The case where c ∈ Γ N and d ∈ Γ D has the following optimal convergence result.Theorem 5.2.Under (2.3), if (σ , u) and (σ h ,u h ) are the respective solutions to (2.2) and (5.1), c ∈ Γ N , and d ∈ Γ D , then, for ∆t sufficiently small, ( The estimate is optimal in both h and p.For other cases of boundary conditions, the error estimate in Theorem 3.2 is sharp for (5.1).As for the second method, with an addition of appropriate boundary terms at x = c, d to (5.1), the resulting scheme can produce optimal error estimates [12].System (5.1) can be similarly implemented (if desired) in nonmixed form as in the second method; the right-hand side limits in (4.5) and (4.6) are replaced by the left-hand side limits.

Proof of Theorem
To prove Theorem 3.2, we need three lemmas whose proofs can be found in [1,18].Lemma 6.1.With definition (2.13), for each n, Moreover, for ∆t n small enough, there are positive constants C 1 and C 2 such that We define the average Note that (2.14) can be written as follows: We are now in a position to prove Theorem 3.2.Recall that Below, is a positive constant independent of h and p, as small as we please.
The error equations follow from the subtraction of (3.1) from (6.8) and (2.17): Set 1 and τ = η n 1 in the first and second equations of (6.9), respectively, and add the resulting two equations to see that with the obvious definition of E n i , i = 1,...,9.Each of the terms in (6.11) is estimated as follows.

Proof of Theorem 4.2.
In the proofs of Theorems 3.2 and 3.3, we have utilized the standard L 2 -projections ᏼ h : L 2 (I) → W h and Π h : L 2 (I) → V h .To prove Theorem 4.2, we use a different pair of operators.For χ ∈ H 1 (I h ), on each interval and Π h χ ∈ V h is given by With this pair of operators, we have the next lemma.Lemma 6.5.Let χ ∈ H l+1 (I i ) with I i = (x i−1 ,x i ) and l ≥ 0. There exits a constant C(l) dependent only on l such that
For simplicity, let the initial approximation u 0 h be given as in (3.2) with ᏼ h being defined in (6.24).Proof theorem 4.2.As for (6.11), it follows from (2.17), (6.8), and (4.2) that Now, by the definition of Π h and ᏼ h , we see that All other terms can be bounded as in the proof of Theorem 3.2.

Proof of Theorem 5.2.
For the proof of Theorem 5.2, the operators ᏼ h : H 1 (I h ) → W h and Π h : H 1 (I h ) → V h are accordingly modified as follows.For χ ∈ H 1 (I h ), on each interval and Π h χ ∈ V h is given by (6.30) Lemma 6.5 remains valid for these new operators, and Theorem 5.2 can be shown using the same technique as for Theorem 4.2.
The corresponding weak formulation is defined as follows: With this, the three families of Eulerian-Lagrangian mixed discontinuous methods in the earlier sections can be accordingly defined.Moreover, all the stability and convergence results hold by combining the present techniques and those in [13,14] for treating the case of a ≥ 0. Note that (7.3) applies to the case where a ≡ 0. In this case, only the inflow boundary condition is needed.

The case of variables φ and b.
For the purpose of demonstration, we have considered the case of constants φ and b.We now generalize to the variable case.Again, for any x ∈ I and two times 0 ≤ t n−1 < t n ≤ T , the hyperbolic part of problem (2.1), φ∂u/∂t + b∂u/∂x, defines the characteristic xn (x, t) along the interstitial velocity ϕ = b/φ: xn x, t n = x. (7.4) In general, we cannot exactly follow the characteristic in (7.4), it can be followed only approximately.There are many ways to solve the first-order ordinary differential equation (7.4).We consider only the Euler method [1].For   (7.9) The second and third methods can be modified in the same manner.The existence and uniqueness result previously obtained holds for (7.9) as well, and the convergence result holds provided the following assumption is satisfied: φ, ϕ ∈ L ∞ J; W 1,∞ (I) .(7.10) 8. Numerical results.We have observed numerical results for the three presented methods similar to those obtained for the corresponding stationary problem in [12] when ∆t is sufficiently small.As an example, we just report numerical results for the second mixed discontinuous method with c ∈ Γ D and d ∈ Γ N .More numerical results can be found in [12].The coefficients in (2.1) are constant (φ = 1), the domain is the unit interval (0, 1), and the exact solution is chosen as follows: u(x, t) = exp(−at) 1 − cos 2π(x − bt) . (8.1) The numerical results for u h at T = 0.1 are shown in Table 8.1 in terms of the space mesh size h and the polynomial degree p.In these results, we have randomly chosen the values of a and b (they have the same magnitude) and have considered uniform grids; analogous results have been observed on nonunigrids.These results are in a good agreement with Theorem 4.2.That is, the estimates are optimal in both h and p.Also, similar numerical results are obtained for σ h .

1 .
The case of a 0. So far, we have assumed that (2.3) holds.When a is nonnegative only, let κ ≥ 0 satisfy a = κ 2 .(7.1)Now, (2.2) is written in the form t n−1,m = t n−1 + m∆t n /M n , with m = M n ,...,1, set xn−1,m−1 n (x) = xn−1,m (x) − ϕ xn−1,m (x), t n−1,m ∆t n M n , xn−1,M n n (x) = x.(7.5)Note that the number of steps M n can depend on n.Denote by xn (x, t) the piecewise linear interpolant in time of xn−1,m (x).If ∆t n /M n is sufficiently small (depending upon the smoothness of ϕ), the approximate characteristics do not cross each other, which is assumed here.Again, indicate the inverse of xn (•,t) by xn (•,t).