ERROR ANALYSIS OF THE SUPG FINITE ... - Semantic Scholar

Comment

Report 2 Downloads 38 Views

ERROR ANALYSIS OF THE SUPG FINITE ELEMENT DISCRETIZATION OF EVOLUTIONARY CONVECTION-DIFFUSION-REACTION EQUATIONS VOLKER JOHN∗ AND JULIA NOVO† Abstract. Conditions on the stabilization parameters are explored for different approaches in deriving error estimates for the SUPG finite element stabilization of time-dependent convectiondiffusion-reaction equations. Exemplarily, it is shown for the SUPG method combined with the backward Euler scheme that standard energy arguments lead to estimates for stabilization parameters that depend on the length of the time step. The stabilization vanishes in the time-continuous limit. However, based on numerical experience, this seems not to be the correct behavior. For this reason, the main focus of the paper consists in deriving estimates in which the stabilization parameters do not depend on the length of the time step. It is shown that such estimates can be obtained in the case of time-independent convection and reaction. An error estimate for the time-continuous case with the standard order of convergence is derived for stabilization parameters of the same form as they are optimal for the steady-state problem. Analogous estimates are obtained for the fully discrete case using the backward Euler method and the Crank–Nicolson scheme. Numerical studies support the analytical results. Key words. Evolutionary convection-diffusion-reaction equation, Streamline-Upwind Petrov– Galerkin (SUPG) finite element method, backward Euler scheme, Crank–Nicolson scheme, timecontinuous problem, error analysis AMS subject classifications. 65M12, 65M60

1. Introduction. Evolutionary convection-diffusion-reaction equations model the transport and reaction of species. In applications, typically the size of the diffusion is much smaller than the size of the convective term and solutions develop sharp layers. In this case, it is well known that standard finite element methods perform poorly and exhibit non-physical oscillations. Stabilization techniques are required in order to get physically sound numerical approximations. This paper studies one of the currently most popular finite element stabilizations, the Streamline-Upwind Petrov–Galerkin (SUPG) method that was introduced for steady-state equations in [9, 2]. In [10], the semi-discrete SUPG formulation was extended to the fully discrete space-time formulation using a discontinuous Galerkin method in time to treat linear convection-diffusion systems. Although time-space elements might be a natural setting to develop stabilized methods for evolutionary equations, effective algorithms for treating such problems are usually designed by separating temporal and spatial discretization. In particular, the increased cost with respect to the number of degrees of freedom for coupled time-space formulations is a significant drawback. Here, the approach of separating spatial and temporal discretization will be considered. Concerning the numerical analysis of this approach, the transient convection equation without diffusive and reactive term is considered in [3]. It is shown that a finite element discretization in space coupled with the backward Euler, the Crank–Nicolson or the second order backward differentiation formula in time leads to the classical error bound for the SUPG method in the L2 norm (suboptimal by an order of one half) ∗ Weierstrass

Institute for Applied Analysis and Stochastics (WIAS), Mohrenstr. 39, 10117 Berlin, Germany and Free University of Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany. † Departamento de Matem´ aticas, Universidad Aut´ onoma de Madrid, Instituto de Ciencias Matem´ aticas CSIC-UAM-UC3M-UCM, Spain. Research supported by Spanish MEC under grants MTM2007-60528 and MTM2010-14919. 1

and to an optimal error bound in the norm of the material derivative. The results are obtained under certain regularity conditions on the data and with stabilization parameters that depend only on the mesh size in the space variable, i.e. δ = O(h), where δ is the SUPG stabilization parameter. However, an optimal error bound for the error in the streamline derivative is not obtained. If the data are not sufficiently smooth or if the velocity field is non-solenoidal, then the bound for the backward Euler method is valid under the condition δ 2 = O(k) and the bound for the Crank–Nicolson scheme holds for δ = O(k), where k is the length of the time step. On the other hand, as before, optimal error bounds in space require to take δ = O(h). An analogous hypothesis for δ, i.e., δ = O(k), was also applied in [15] where a Galerkin least-squares method in space coupled with a θ-scheme in time is analyzed. The analysis of [15] excludes the case θ = 1/2. Finally, the stability of the SUPG finite element method for transient convection-diffusion equations is studied in [1]. However, as it is shown already in [3], the coercivity result of [1] leads to suboptimal global estimates in time. Numerical studies of the SUPG method, together with a discussion on relations to other stabilized finite element methods, can be found in [6]. In [12, 13], the SUPG method was compared comprehensively with other stabilized finite element methods. The approach in these studies was as follows: 1) discretize the equation in time, 2) consider the equation in each discrete time as a steady-state convection-diffusionreaction equation, 3) discretize this equation in space with a stabilized method and apply a parameter choice that is appropriate for this type of steady-state equation. This methodology leads to a parameter that is (in the notation of formula (3.1) below) proportional to the length of the time step, see formulae (8) and (11) in [12]. The numerical results with this approach show large spurious oscillations compared with other methods. Such oscillations can be observed also if the SUPG method with this parameter choice is used in coupled systems coming from applications, as in [11]. Altogether, the numerical results obtained so far are not at all satisfactory. We think that a reason for this is the choice of the stabilization parameters depending on the length of the time step. In the case of a uniform time-space grid, i.e. k ∼ h, this choice gives the size of the stabilization parameters as it is known from the steadystate case. However, approaching the time-continuous limit, i.e. using an anisotropic grid in time-space, the SUPG stabilization vanishes and one reverts to the Galerkin method that is unstable for the convection-dominated case, see also [8] for the same opinion. Very small time steps might be necessary if, e.g., problems with very fast reactions are simulated. In addition, it is desirable to have balanced temporal and spatial errors. This balance is usually given by the orders of the methods. However, additional conditions coming from the coupling of the length of the time step and the mesh width via the stabilization parameters might contradict the balances given by these orders. In [8], another approach for deriving the fully discrete equation is considered: 1) discretize the equation in space with a stabilized method, 2) choose standard stabilization parameters for this equation, 3) discretize the equation in time. Because the temporal discretization is performed after the choice of the stabilization parameters, these parameters cannot depend on the time step. Numerical studies in [8] show that this approach leads to much more stable results for small time steps compared with the approach from [12, 13]. In addition, another parameter choice is proposed in [8] that, e.g., does not depend on the length of the time step if a steady-state solution is approached, the so-called element-vector-based parameter. The goal of the present paper consists in exploring the conditions on the stabilization parameters for different approaches in the numerical analysis for deriving 2

error estimates. In particular, error estimates that do not lead to a dependency of the stabilization parameters on the length of the time step are of interest. To our best knowledge, error estimates of this kind for the SUPG method applied to evolutionary convection-diffusion-reaction equations are not yet available. The main difficulty in the analysis of the method comes from the fact that the time derivative has to enter the stabilization term in order to ensure consistency. This adds a non-symmetric term that cannot be easily bounded using standard energy arguments with stabilization parameters that do not depend on the length of the time step. In the first part of the paper, exemplarily the backward Euler scheme is considered as temporal discretization. This is just for illustrating the difficulties in the standard analysis in the most simple setting. In Sections 3 and 4, stability bounds and error estimates are derived based on energy arguments. Two different ways to argue lead to error estimates under the conditions δ = O(k) and δ = O(k 1/2 h), respectively. These conditions arise in the stability bounds from the stabilization term with the discretization of the time derivative. In both choices, the stabilization parameters tend to zero on a fixed spatial grid as the length of the time step approaches zero. As discussed above, this seems not to be the correct choice. This is also seen in numerical studies, e.g., in Example 6.2 below. Altogether, the limit of the time-continuous case could not be treated so far satisfactorily by standard energy arguments. To obtain some insight in the time-continuous case, Section 5 studies a special problem, where the convection field and the reaction do not depend on time and the SUPG method is applied on a uniform grid. The stabilization parameters are chosen to be the same on all mesh cells, depending only on the coefficients of the equation and on h: δ = O(h). Under certain regularity assumptions on the solution and extending the analysis of [3], an error estimate for the L2 norm and the norm of the material derivative is derived with the standard order of convergence. In the next step, based on this result, an estimate for the error in the norm of the streamline derivative is proven with the same order of convergence. To our best knowledge, this is the first result that proves standard order of convergence for the SUPG method applied to evolutionary convection-diffusion-reaction equations with a parameter choice that is essentially the same as in the steady-state case. A similar result was obtained in [7] for subgrid viscosity stabilizations of Galerkin approximations. Then, the analysis is extended to the fully discrete case using the backward Euler method and the CrankNicolson scheme as time integrators, respectively. The last part of the paper, Section 6, presents some numerical studies. In a first example, the derived error estimates are supported. Then, a rotating body problem is studied for the P1 finite element, on a given spatial grid, and for a very small length of the time step. The results show clearly that in this situation a choice of the stabilization parameter independently of the length of the time step has to be preferred. The paper concludes in Section 7 with a summary of the results and an outlook to open questions. 2. The SUPG Method and Preliminaries of the Analysis. Throughout this paper, standard notations are used for Lebesgue and Sobolev spaces. Generic constant that do not depend on the mesh width or the length of the time step are denoted by C. A linear time-dependent convection-diffusion-reaction equation is given by ut − ε∆u + b · ∇u + cu = f in (0, T ] × Ω, u = 0 on [0, T ] × ∂Ω, u(0, x) = u0 (x) in Ω, 3

(2.1)

where Ω is a bounded open domain in Rd , d ∈ {1, 2, 3}, with boundary ∂Ω, b(t, x) and c(t, x) are given functions, ε > 0 is a constant diffusion coefficient, u0 (x) are given initial data and T is a given final time. For simplicity, the case that Ω is a convex polygonal or polyhedral domain is considered. In the following, it is assumed that there is a constant µ0 > 0 such that 1 0 < µ0 ≤ µ(t, x) = c − ∇ · b (t, x), ∀ (t, x) ∈ [0, T ] × Ω. (2.2) 2 This is a standard assumption in the analysis of equations of type (2.1), [16]. Let V = H01 (Ω). A variational form of (2.1) reads as follows: Find u : (0, T ] → V such that (ut , v) + (ε∇u, ∇v) + (b · ∇u + cu, v) = (f, v) ∀ v ∈ V,

(2.3)

and u(0, x) = u0 (x). Here, (·, ·) denotes the inner product in L2 (Ω)d , d ∈ {1, 2, 3}. In numerical simulations, V is replaced by a finite dimensional (finite element) space Vh,r , where h indicates the fineness of the underlying triangulation Th and r ∈ N the degree of the local polynomials. This paper considers the case of a conforming finite element method, i.e. Vh,r ⊂ V . The time-continuous finite element problem aims to find a function uh ∈ Vh,r that fulfills a problem of form (2.3) for all test functions from Vh,r with an appropriate approximation of u0 (x) at the initial time. Using some temporal discretization, a finite element Galerkin method for solving (2.3) is obtained. It is well known that in the case of small diffusion, in particular compared with the convection, the Galerkin method is instable and leads to solutions that are globally polluted with huge spurious oscillations. A stabilization of the Galerkin method becomes necessary. The probably most popular stabilized finite element method is the SUPG method. This residual-based method adds artificial diffusion along the streamlines of the solution. It has the form (time-continuous case): Find uh : (0, T ] → Vh,r such that X (uh,t , vh ) + aSUPG (uh , vh ) + δK (uh,t , b · ∇vh )K K∈Th

= (f, vh ) +

X

δK (f, b · ∇vh )K

∀ vh ∈ Vh,r ,

K∈Th

with uh (0, x) being an appropriate approximation of u0 (x) and aSUPG (uh , vh ) = ε(∇uh , ∇vh ) + (b · ∇uh , vh ) + (cuh , vh ) X + δK (−ε∆uh + b · ∇uh + cuh , b · ∇vh )K .

(2.4)

K∈Th

Here, K ∈ Th denotes the mesh cells of the triangulation, (·, ·)K the inner product in L2 (K), and {δK } are local parameters which has to be chosen appropriately. Next, preliminaries for the analysis are introduced. The elliptic projection πh : V → Vh,r is defined by (∇(u−πh u), ∇vh ) = 0 for all vh ∈ Vh,r . Note that the functions of Vh,r do not depend on time. Hence, if follows (πh u)t = πh (ut ) = πh ut .

(2.5)

Assuming that the meshes are quasi-uniform, the following inverse inequality holds for each vh ∈ Vh,r , see, e.g., [4, Theorem 3.2.6], l−m−d q10 − q1

kvh kW m,q (K) ≤ cinv hK

4

kvh kW l,q0 (K) ,

(2.6)

where 0 ≤ l ≤ m ≤ 1, 1 ≤ q 0 ≤ q ≤ ∞, hK is the size (diameter) of the mesh cell K ∈ Th and k · kW m,q (K) is the norm in W m,q (K). The following interpolation error estimate for u ∈ V ∩ H r+1 (Ω) is well known, [5, 18] ku − πh uk0 + hku − πh uk1 ≤ Chr+1 kukr+1 ,

(2.7)

where k · kr denotes the norm in H r (Ω) with H 0 (Ω) = L2 (Ω). In particular, stability estimates for u ∈ H01 (Ω) of the following form can be derived kπh uk0 ≤ ku − πh uk0 + kuk0 ≤ Chkuk1 + kuk0 ≤ Ckuk1 .

(2.8)

It is assumed that the space Vh,r satisfies the following local approximation property: for each u ∈ V ∩ H r+1 (Ω) there exists a u ˆh ∈ Vh,r such that ku − u ˆh k0,K + hK k∇(u − u ˆh )k0,K + h2K k∆(u − u ˆh )k0,K ≤ Chr+1 K kukr+1,K

(2.9)

for all K ∈ Th . For example, this property is given for Lagrange finite elements on mesh cells which allow an affine transform to a reference mesh cell. Lemma 2.1. With (2.9) follows for all u ∈ V ∩ H r+1 (Ω) X k∆(u − πh u)k20,K ≤ Ch2r−2 kuk2r+1 . (2.10) K∈Th

Proof. The result is easily obtained using the triangle inequality, the local approximation property (2.9), the inverse inequality (2.6), the quasi-uniformity of the mesh, and the interpolation error estimate (2.7). The coercivity of the bilinear form aSUPG (·, ·) under the condition that the parameters {δK } are appropriately bounded from above is a well-known result, see, e.g., [17, Lemma 10.3]. Lemma 2.2. Coercivity of aSUPG (·, ·). Let (2.2) be satisfied. If the SUPG parameters are chosen such that δK ≤

µ0 , 2kck2K,∞

δK ≤

h2K , 2εc2inv

(2.11)

then the bilinear form aSUPG (·, ·) associated with the SUPG method satisfies aSUPG (uh , uh ) ≥

1 kuh k2SUPG 2

(2.12)

with !1/2 kuh kSUPG :=

εk∇uh k20

+

X

δK kb ·

∇uh k20,K

1/2

+ kµ

uh k20

.

K∈Th

For linear finite elements, the condition δK ≤ h2K /(2εc2inv ) can be omitted. 3. Stability for stabilization parameters depending on the length of the time step. This section studies a fully discrete method for solving (2.3). Besides the finite element SUPG discretization (2.4), the temporal derivative is approximated with the backward Euler scheme. For this simple setting, standard energy arguments are applied to derive stability bounds. It turns out that this analysis proposes parameter choices in the SUPG method that depend on the length of the time step. 5

Consider the case of a fixed time step k = ∆t. The fully discrete solution at time tn = nk will be denoted by Uhn . The backward Euler/SUPG method reads as follows: For n = 1, 2, . . . find Uhn ∈ Vh,r such that n Uh − Uhn−1 , ϕ + ε(∇Uhn , ∇ϕ) + (b · ∇Uhn , ϕ) + (cUhn , ϕ) = (f n , ϕ) k n X Uh − Uhn−1 + δK f n − + ε∆Uhn − b · ∇Uhn − cUhn , b · ∇ϕ (3.1) k K K∈Th

for all ϕ ∈ Vh,r and Uh0 (x) = uh (0, x). Method (3.1) can be written equivalently in the form X (Uhn − Uhn−1 , ϕ) + kaSUPG (Uhn , ϕ) = k(f n , ϕ) + k δK (f n , b · ∇ϕ)K K∈Th

−

X

δK (Uhn − Uhn−1 , b · ∇ϕ)K .

(3.2)

K∈Th

Theorem 3.1. Stability, stabilization parameters proportional to the length of the time step. Let (2.2) and (2.11) be fulfilled. With the additional condition δK ≤

k 4

∀ K ∈ Th ,

(3.3)

the solution of (3.1) satisfies at tn = nk kUhn k20 +

X n n kX j 2 2 kUh kSUPG ≤ kUh0 k20 + k kf j k20 . +k 2 j=1 µ0 j=1

Proof. The proof starts in the usual way by setting ϕ = Uhn . This gives with (3.2) X (Uhn − Uhn−1 , Uhn ) + kaSUPG (Uhn , Uhn ) = k(f n , Uhn ) + k δK (f n , b · ∇Uhn )K K∈Th

−

X

δK (Uhn

− Uhn−1 , b · ∇Uhn )K .

K∈Th

From (Uhn − Uhn−1 , Uhn ) =

1 2

kUhn k20 − kUhn−1 k20 + kUhn − Uhn−1 k20 follows with (2.12)

k 1 kUhn k20 − kUhn−1 k20 + kUhn − Uhn−1 k20 + kUhn k2SUPG (3.4) 2 2 X X ≤ |k(f n , Uhn )| + k δK (f n , b · ∇Uhn )K + δK (Uhn − Uhn−1 , b · ∇Uhn )K . K∈Th

K∈Th

The first two terms on the right hand side are estimated using the Cauchy–Schwarz inequality and Young’s inequality

n 2

f k 1/2 n 2 k n 2 k 1/2 n 2 n n

|k(f , Uh )| ≤ k

µ1/2 + 4 kµ Uh k0 ≤ µ0 kf k0 + 4 kµ Uh k0 , 0 X X k X δK (f n , b · ∇Uhn )K ≤ 2k δK kf n k20,K + δK kb · ∇Uhn k20,K . k 8 K∈Th

K∈Th

K∈Th

6

The estimate of the last term on the right hand side of (3.4) uses condition (3.3) on the stabilization parameters X n−1 n n δK (Uh − Uh , b · ∇Uh )K K∈Th

≤

2 X k X δK kUhn − Uhn−1 k20,K + δK kb · ∇Uhn k20,K k 8 K∈Th

K∈Th

1 k X δK kb · ∇Uhn k20,K . ≤ kUhn − Uhn−1 k20 + 2 8 K∈Th

Inserting all estimates leads to kUhn k20 +

X k n 2 2k n 2 kUh kSUPG ≤ kUhn−1 k20 + kf k0 + 4k δK kf n k20,K . 2 µ0

(3.5)

K∈Th

Summation of the time steps j = 1, . . . , n, and using once more condition (3.3) gives the statement Pnof the theorem. Note k j=1 kUhj k2SUPG is an approximation of kUh k2L2 (0,T ;SUPG) by a Riemann sum using as node in the quadrature rule always the right end of the time intervals. Theorem 3.1 covers the case that the stabilization parameter is proportional to the length of the time step. On a fixed spatial grid, the stabilization becomes small for small time steps and it vanishes in the time-continuous limit. This behavior does not seem to be correct, see the discussion in the introduction. The desired situation in the convection-dominated regime, δK ∼ hK , is obtained if spatial and temporal mesh width are proportional h ∼ k. Note that for the mesh width and the time step being of the same order, the parameter choice of [12, 13] leads also to δ ∼ k ∼ h. Theorem 3.2. Stability, stabilization parameters proportional to some function of the length of the time step. Let (2.2) and (2.11) be fulfilled. With the choice δK =

σ(k)hK kbk∞,K cinv

0 < σ(k) ≤

with

1 4

∀ K ∈ Th ,

(3.6)

where σ(k) is a function to be specified later, the solution of (3.1) satisfies at tn = nk kUhn k20 +

n kX j 2 kU k 2 j=1 h SUPG 

≤ (1 + 2σ 2 (k))n kUh0 k20 + 2k

n X j=1

! X 1 j 2 kf k0 + δK kf j k20,K  . µ0

(3.7)

K∈Th

Proof. The proof starts exactly as the proof of Theorem 3.1 until estimate (3.4) is reached. The first two terms on the right hand side of (3.4) are estimated also in the same way as in the proof of Theorem 3.1 k n 2 k 1/2 n 2 |k(f n , Uhn )| ≤ kf k0 + kµ Uh k0 , µ0 4 X X k X δK (f n , b · ∇Uhn )K ≤ k δK kf n k20,K + δK kb · ∇Uhn k20,K . k 4 K∈Th

K∈Th

K∈Th

7

The last term on the right hand side of (3.4) will now not be absorbed into k2 kUhn kSUPG . It is estimated by using the inverse inequality (2.6) and Young’s inequality X n−1 n n δK (Uh − Uh , b · ∇Uh )K K∈Th

≤

X

δK

K∈Th

+

X

kbk∞,K cinv n kUh − Uhn−1 k20,K hK

δK kbk∞,K kUhn − Uhn−1 k0,K k∇Uhn−1 k0,K

K∈Th

≤

X

δK

K∈Th

K∈Th

+

X 1 kbk∞,K cinv n kUh − Uhn−1 k20,K + kU n − Uhn−1 k20,K hK 4 h

X

2 δK kbk2∞,K k∇Uhn−1 k20,K

K∈Th

X X kbk∞,K cinv kbk2∞,K c2inv n−1 2 1 2 kUhn − Uhn−1 k20,K + kUh k0,K . + δK ≤ δK hK 4 h2K K∈Th

K∈Th

The first term can be absorbed into the left hand side of (3.4) if δK

kbk∞,K cinv 1 hK 1 =⇒ δK ≤ + ≤ . hK 4 2 4kbk∞,K cinv

Set the stabilization parameter as in (3.6), then it follows X 1 n−1 n n δK (Uh − Uh , b · ∇Uh )K ≤ kUhn − Uhn−1 k20 + σ 2 (k)kUhn−1 k20 . 2 K∈Th

Collecting all estimates leads to the recursion X k 2k kUhn k20 + kUhn k2SUPG ≤ (1 + 2σ 2 (k))kUhn−1 k20 + kf n k20 + 2k δK kf n k20,K . (3.8) 2 µ0 K∈Th

Now, one obtains by induction kUhn k20 +

k n 2 kU k 2 h SUPG 2

≤ (1 + 2σ (k))

n

kUh0 k20

+ 2k

n X

n−j 1 + 2σ (k) 2

j=1

 ≤ (1 + 2σ 2 (k))n kUh0 k20 + 2k

n X j=1

X kf j k20 + δK kf j k20,K µ0 K∈Th !

X kf j k20 + δK kf j k20,K  . µ0 K∈Th

Summation of (3.8) gives kUhn k20 +

n n−1 X j kX j 2 kUh kSUPG ≤ 2σ 2 (k) kUh k20 + (1 + 2σ 2 (k))kUh0 k20 2 j=1 j=1 ! n X X kf j k20 j 2 +2k + δK kf k0,K . µ0 j=1 K∈Th

8

!

(3.9)

Inserting (3.9) and applying some estimates for the sake of simplifying the representation lead to n kX j 2 kUhn k20 + kU k 2 j=1 h SUPG (1 + 2σ 2 (k))n − (1 + 2σ 2 (k)) 2 + 1 + 2σ (k) kUh0 k20 ≤ 2σ 2 (k) 1 + 2σ 2 (k) − 1   ! n−1 n j 2 X X X kf k 0 +2k 1 + 2σ 2 (k) (1 + 2σ 2 (k))j  + δK kf j k20,K  µ 0 j=1 j=1 K∈Th  ! n j 2 X X kf k 0 ≤ (1 + 2σ 2 (k))n kUh0 k20 + 2k + δK kf j k20,K  . µ 0 j=1 K∈Th

Consider a finite time interval [0, T ] and a fixed length of the time step. Then, Theorem 3.2 gives stability with the desired stability parameter (in the convectiondominated regime) δK = O(hK ) without a coupling of the mesh width to the time step by choosing σ(k) = const ≤ 1/4. However, the stability bound blows up for σ(k) = const in the time-continuous limit k → 0. Given a length of the time step k, the number of time steps to solve the equation in [0, T ] is n = T /k. The stability 2 1/k estimate will not blow is bounded uniformly. A possible √ up for k → 0 if (1 + σ (k)) choice is σ(k) = δ0 k leading to the stabilization parameter √ khK δK = δ0 , (3.10) kbk∞,K cinv √ where δ0 has to be chosen such that δ0 k ≤ 1/4. For fixed h and sufficiently small k, the parameter from (3.10) is larger than the parameter from (3.3). 4. Error estimates for stabilization parameters depending on the length of the time step. For the following error analysis, it is assumed that all functions are sufficiently regular. Summaries of these assumptions are given below in Theorem 4.1. The error analysis for (3.1) starts by decomposing the error into an interpolation error and the difference of the interpolation and the solution Uhn − u(tn ) = (Uhn − πh u(tn )) + (πh u(tn ) − u(tn )) . The interpolation error can be estimated with (2.7). For brevity, denote πhn u := πh u(tn ) and enh = Uhn − πh u(tn ). Straightforward calculations yield the following error equation n n (enh − en−1 , ϕ) + kaSUPG (enh , ϕ) = k(T˜zero , ϕ) + k(Tconv , ϕ) Xh X n−1 n n ˜ +k δK (Tstab,K , b · ∇ϕ)K − δK (eh − eh , b · ∇ϕ)K , K∈Th

K∈Th

with π n u − πhn−1 u n T˜zero = (ut (tn ) − πhn ut ) + c (u(tn ) − πhn u) + πhn ut − h , k n Tconv = b · ∇(u(tn ) − πhn u), n n n T˜stab,K = T˜zero + Tconv + ε∆(πhn u − u(tn )) . K

9

Using integration by parts and assuming δK > 0, the convective term can be distributed to the term with the zeroth order derivatives (with respect to space) and the stabilization term n X πh u − u(tn ) n n (Tconv , ϕ) = −((∇ · b)(πh u − u(tn )), ϕ) − δK , b · ∇ϕ . δK K K∈Th

Redefining the zeroth order and the stabilization term n n Tzero = T˜zero − (∇ · b)(πhn u − u(tn )),

π n u − u(tn ) n n Tstab,K = T˜stab,K − h , δK

leads to the error equation n (enh − en−1 , ϕ) + kaSUPG (enh , ϕ) = k(Tzero , ϕ) + k h

X

n δK (Tstab,K , b · ∇ϕ)K

K∈Th

−

X

δK (enh

− ehn−1 , b · ∇ϕ)K .

(4.1)

K∈Th

This error equation is similar to equation (3.1), only the arguments on the first two terms on the right hand side are not the same. Deriving error estimates from (4.1) starts essentially in the same way as the derivation of the stability bounds. After this, the arising terms have to be bounded by norms of the solution of the continuous equation (2.3). Since the stability bounds derived in Theorems 3.1 and 3.2 are similar, the detailed analysis for the error estimates is presented here only for the case that was considered in Theorem 3.1. For applying the techniques of the stability estimate to (4.1), only the last two terms cannot be combined in the summation of the analog to (3.5). One gets kenh k20 +

n n n X X 2k X j kX j 2 j keh kSUPG ≤ ke0h k20 + kTzero k20 + 4k δK kTstab,K k20,K . 2 j=1 µ0 j=1 j=1 K∈Th

(4.2) Using the triangle inequality and (2.7), one obtains j kTzero k20 ≤ Ch2r+2 kut (tj )k2r+1 + kck2L∞ (0,T ;L∞ ) ku(tj )k2r+1

2

πhj u − πhj−1 u

2 2 +k∇ · bkL∞ (0,T ;L∞ ) ku(tj )kr+1 + C πh ut (tj ) −

.

k 0

The last term is in essence the approximation error of ut (tj ) by a backward finite difference, hence an estimate of O(k) can be expected. The derivation of this estimate uses Taylor’s formula with remainder in integral form, the application of (2.5), the Cauchy–Schwarz inequality, and the stability estimate (2.8)

2

Z

2

t πhj u − πhj−1 u 1

j

(t − tj−1 )πh utt dt

πh ut (tj ) −

= 2

k k tj−1 0 0  !1/2 Z !1/2 2 Z tj tj 1   (t − tj−1 )2 dt kπh utt k20 dt ≤ 2 k tj−1 tj−1 Z tj ≤ Ck kutt k21 dt = Ckkutt k2L2 (tj−1 ,tj ;H 1 ) . (4.3) tj−1

10

Summation over the time steps, taking into account that the number of time steps n is inverse proportional to the length of the time step, and assuming that all norms are uniformly (in time) bounded give k

n X

j kTzero k20 ≤ Cknh2r+2 + Ck 2 kutt k2L2 (0,tn ;H 1 ) ≤ C h2r+2 + k 2 .

j=1

The estimate of the first term can be applied, in combination with (2.10), to obtain an estimate for the second term on the right hand side of (4.2) X j 2 δK kTstab,K k0,K ≤ C max δK h2r+2 kut (tj )k2r+1 + ku(tj )k2r+1 K∈Th

K∈Th

+kkutt k2L2 (tj−1 ,tj ;H 1 ) + kbk2L∞ (0,T ;L∞ ) h2r ku(tj )k2r+1 −1 2 2r−2 2 +ε h ku(tj )kr+1 + C min δK h2r+2 ku(tj )k2r+1 . K∈Th

Hence, k

n X X

j δK kTstab,K k20,K

≤C

h2r+2 + k 2 + h2r + ε2 h2r−2

max δK

K∈Th

j=1 K∈Th

+

!

−1

min δK

2r+2

h

K∈Th

.

Inserting all estimates into (4.2) and applying the triangle inequality leads to the following error estimates. Theorem 4.1. Error estimates for the stabilization parameter obeying (3.3). Suppose b ∈ L∞ (0, T ; (L∞ )d ), ∇ · b, c ∈ L∞ (0, T ; L∞ ) for the coefficients in (2.3) and u, ut ∈ L∞ (0, T ; H r+1 ), utt ∈ L2 (0, T ; H 1 ) for the solution of (2.3). Let the stabilization parameters {δK } fulfill (2.11), (3.3) and δK > 0 for all K ∈ Th . Denote δ = max∈Th δK . Then, the error Uhn − u(tn ) satisfies " n kUh − u(tn )k0 ≤ C hr+1 + k + hr−1 δ 1/2 h2 + h + ε #

hr+1

+

1/2

(minK∈Th δK ) and  k

n X

1/2 kUhj

− u(tj )k2SUPG 

+ kπh u0 −

Uh0 k0

,

(4.4)

"

≤ C hr (ε1/2 + δ 1/2 + h) + k + hr−1 δ 1/2 h2 + h + ε

j=1

+

#

hr+1 1/2

(minK∈Th δK )

+ kπh u0 −

Uh0 k0

,

(4.5)

where the constants C depend on u, ut , utt , b, ∇ · b and c. Applying the analysis of Theorem 3.2 to estimate (4.1) and using (3.7) leads essentially to (4.2), only with an additional factor of (1 + 2σ 2 (k))n on the right hand 11

side. The same analysis as in the proof of Theorem 4.1 gives error estimates of form (4.4) and (4.5) with this additional factor. n An alternative approach, without integration by parts of Tconv , leads to a subopr timal estimate with respect to space of order h . However, δK does not appear in the denominator of some term in this estimate. Hence, this suboptimal estimate shows that the error is bounded even as δK → 0. 5. Error analysis with stabilization parameters depending not on the length of the time step. The numerical analysis presented so far is only valid if, for a constant mesh and a small time step, the stabilization parameters are sufficiently small. In the time-continuous limit, the SUPG stabilization even vanishes. As discussed in the introduction and as demonstrated in the numerical studies, Example 6.2, we think that this is not the correct asymptotic of the stabilization parameters. In this section, optimal error estimates are derived with stabilization parameters proportional to the mesh width. In the first part, the time-continuous problem is considered while the second and third part extend the analysis to fully discrete cases using the backward Euler and the Crank-Nicolson scheme, respectively. 5.1. The time-continuous case. In the first step, an error estimate for the material derivative is derived, Lemma 5.1. The analysis of this step uses some ideas from [3], like the application of a special test function to obtain (5.9). Extensions of the analysis from [3] were necessary to include diffusion and reaction and the case of a non-solenoidal vector field b. Based on the estimate for the material derivative, an error estimate for the streamline derivative is proven in a second step. Consider problem (2.1) with the following assumptions: • bt (t, x) = 0, ct (t, x) = 0, i.e., b = b(x), c = c(x) and µ = µ(x), • the mesh is uniform with mesh width h, • the stabilization parameters are the same for all mesh cells, i.e. δK = δ. Furthermore, it is assumed that all functions are sufficiently smooth such that all norms appearing below are well defined. In addition, only the convection-dominated regime is considered, i.e. it is assumed that ε ≤ h. Then, the stabilization parameter is set to be ( ) ) ( 1/2 1/2 1/2 µ0 µ0 µ0 1 µ0 h , , min , , , kbk1/2 , 1 . (5.1) δ = min ∞ 4cinv kbk∞ 2 4 4kck∞ 2kck∞ kck1/2 ∞ Hence, the stabilization parameter is proportional to the mesh width and it is bounded from above by a constant. In addition, it is assumed that ε ≤ δ. Consider a finite time interval [0, T ] and let t ∈ [0, T ]. In the analysis of this section, a formally steady-state problem derived from (2.1) is used. Let Πh u(t) ∈ Vh,r be the solution of aSUPG (Πh u(t), vh ) = (f (t) − ut (t), vh ) + δ(f (t) − ut (t), b · ∇vh ) ∀ vh ∈ Vh,r . (5.2) The corresponding continuous equation is solved by u(t). Hence, firstly the Galerkin orthogonality of the SUPG method gives aSUPG (Πh u(t), vh ) = aSUPG (u(t), vh ) ∀ vh ∈ Vh,r .

(5.3)

Secondly, error estimates of the form ku(t) − Πh u(t)kSUPG ≤ Chr+1/2 ku(t)kr+1 12

t ∈ [0, T ],

(5.4)

can be proven, see [16]. A direct calculation, using the linearity of the equation and the time-independency of convection, reaction and the test functions, shows (Πh u(t))t = Πh (ut (t)) = Πh ut .

(5.5)

For brevity, the dependency on time will be omitted from now in the notations. Let uh : (0, T ] → Vh,r be the finite element solution of the continuous-in-time SUPG method (uh,t , vh ) + aSUPG (uh , vh ) = (f, vh ) + δ(f − uh,t , b · ∇vh ) ∀ vh ∈ Vh,r

(5.6)

with uh (0) given. For the error analysis, the following norms in Vh,r are introduced kvh kb := kvh k20 + δ 2 kb · ∇vh k20

1/2

,

kvh kmat := δ 1/2 kvh,t + b · ∇vh k0 .

The expression in the second norm is the material derivative. Note, k · kb is equivalent to the L2 norm, since by using the inverse inequality and the definition (5.1) of the stabilization parameter, one obtains √ 17 2 2 2 2 −2 2 1/2 kvh k0 . kvh k0 ≤ kvh kb ≤ kvh k0 + δ kbk∞ cinv h kvh k0 ≤ 4 Denote the error between the continuous-in-time finite element solution and the solution of the steady-state problem by eh = uh − Πh u and let Ttr = ut − Πh ut be the truncation error. An error equation is obtain by subtracting (5.2) from (5.6) (eh,t , vh )+aSUPG (eh , vh ) = (Ttr , vh )+δ(Ttr , b·∇vh )−δ(eh,t , b·∇vh ) ∀ vh ∈ Vh,r . (5.7) Setting vh = eh in (5.7) gives 1 d keh k20 + εk∇eh k20 + δkb · ∇eh k20 + kµ1/2 eh k20 + δ(eh,t , b · ∇eh ) 2 dt X = (Ttr , eh + δb · ∇eh ) − δ(ceh , b · ∇eh ) + δε(∆eh , b · ∇eh )K .

(5.8)

K∈Th

Analogously, one obtains for vh = eh,t in (5.7) ε d δ d 1 d 1/2 k∇eh k20 + (b · ∇eh , eh,t ) + kb · ∇eh k20 + kc eh k20 2 dt 2 dt X 2 dt = (Ttr , (eh + δb · ∇eh )t ) − δ(ceh , b · ∇eh,t ) + δε(∆eh , b · ∇eh,t )K (5.9)

keh,t k20 +

K∈Th

−δ(eh,t , b · ∇eh,t ). The addition of δ times (5.9) to (5.8) leads to 1 d εδ d δ d 1/2 keh k2b + εk∇eh k20 + kµ1/2 eh k20 + keh k2mat + k∇eh k20 + kc eh k20 2 dt 2 dt 2 dt = (Ttr , eh + δb · ∇eh ) + δ(Ttr , (eh + δb · ∇eh )t ) X −δ(ceh , b · ∇eh ) − δ 2 (ceh , b · ∇eh,t ) + δε(∆eh , b · ∇eh )K (5.10) K∈Th

+

X

2

2

δ ε(∆eh , b · ∇eh,t )K − δ (eh,t , b · ∇eh,t ),

K∈Th

13

where the definition of k · kb and δ keh,t k20 + 2(eh,t , b · ∇eh ) + kb · ∇eh k20 = keh k2mat have been used. Now, the goal consists in estimating all terms on the right hand side of (5.10) and absorbing the terms with eh and eh,t into the left hand side. Note, all terms with eh,t have to be absorbed into the norm of the material derivative since this is the only term on the left hand side which includes the temporal derivative. The estimate of all terms on the right hand side of (5.10) starts with the Cauchy– Schwarz inequality. Observe that keh + δb · ∇eh k0 ≤ keh k0 + δkbk∞ cinv h−1 keh k0 −1/2

≤ µ0

−1/2

kµ1/2 eh k0 + µ0

δkbk∞ cinv h−1 kµ1/2 eh k0 .

Assuming (5.1) gives for the first term on the right hand side of (5.10), together with the application of Young’s inequality, 1 −1/2 (Ttr , eh + δb · ∇eh ) ≤ kTtr k0 µ0 kµ1/2 eh k0 + kTtr k0 kµ1/2 eh k0 8 1 1 1/2 −1 2 ≤ 4µ0 + kTtr k0 + kµ eh k20 . 16 8

(5.11)

For estimating the second term on the right hand side of (5.10), a norm of the temporal derivative of the error has to be bounded. With (5.1), one obtains δkeh,t + δb · ∇eh,t k0 ≤ δ(keh,t k0 + δkbk∞ cinv h−1 keh,t k0 ) ≤ 2δkeh,t k0 ≤ 2δ(keh,t + b · ∇eh k0 + kb · ∇eh k0 ) −1/2

≤ 2δ 1/2 keh kmat + 2δkbk∞ cinv h−1 µ0 1 ≤ 2δ 1/2 keh kmat + kµ1/2 eh k0 . 8

kµ1/2 eh k0

It follows

1 1/2 δ(Ttr , (eh + δb · ∇eh )t ) ≤ kTtr k0 2δ keh kmat + kµ eh k0 (5.12) 8 1 1 1 ≤ 16δ + kTtr k20 + keh k2mat + kµ1/2 eh k20 . 16 16 16 1/2

The third term on the right hand side of (5.10) is estimated simply by 1/2 δ(ceh , b · ∇eh ) ≤ δkck∞ kbk∞ cinv h−1 µ−1 eh k20 ≤ 0 kµ

1 1/2 kµ eh k20 , 16

(5.13)

where (5.1) has been used. Reasoning in the same way for the fourth term gives −1/2

δ 2 (ceh , b · ∇eh,t ) ≤ δ 2 kck∞ kbk∞ cinv h−1 µ0

kµ1/2 eh k0 keh,t k0

−1/2

≤ δ 3/2 kck∞ kbk∞ cinv h−1 µ0

kµ1/2 eh k0 keh kmat

1/2 +δ 2 kck∞ kbk2∞ c2inv h−2 µ−1 eh k20 . 0 kµ

Assuming (5.1) yields δ 1/2 1/2 1 kµ eh k0 keh kmat + kµ1/2 eh k20 8 16 δ 1 1 ≤ + kµ1/2 eh k20 + keh k2mat 16 16 16 1 1/2 1 ≤ kµ eh k20 + keh k2mat . 8 16

δ 2 (ceh , b · ∇eh,t ) ≤

14

(5.14)

For the fifth term on the right hand side of (5.10), one obtains X X ε δε(∆eh , b · ∇eh )K ≤ δεcinv h−1 kbk∞,K k∇eh k20,K ≤ k∇eh k20 , 2 K∈Th

(5.15)

K∈Th

if (5.1) is fulfilled. The estimate of the sixth term starts as follows X X δ 2 ε(∆eh , b · ∇eh,t )K ≤ δ 2 εc2inv h−2 kbk∞,K k∇eh k0,K keh,t k0,K K∈Th

K∈Th

≤

X

δ 3/2 εc2inv h−2 kbk∞,K k∇eh k0,K keh kmat,K

K∈Th

+

X

δ 2 εc2inv h−2 kbk∞,K k∇eh k0,K kb · ∇eh k0,K

K∈Th

≤

X

δ 3/2 εc2inv h−2 kbk∞,K k∇eh k0,K keh kmat,K

K∈Th

+

X

δ 2 εc2inv h−2 kbk2∞,K k∇eh k20,K .

K∈Th

Now, assuming (5.1) and ε ≤ δ leads to X δ 2 ε(∆eh , b · ∇eh,t )K K∈Th

≤

X ε X εδ −1/2 k∇eh k0,K keh kmat,K + k∇eh k20,K 16 16 K∈Th

K∈Th 2

≤

ε ε 3 1 ε k∇eh k20 + keh k2mat + k∇eh k20 ≤ εk∇eh k20 + keh k2mat . (5.16) 32δ 32 16 32 32

Finally, the last term on the right hand side of (5.10) is estimated as follows δ 2 (eh,t , b · ∇eh,t ) ≤ δ 2 kbk∞ h−1 cinv keh,t k20 ≤ 2δ 2 kbk∞ h−1 cinv (keh,t + b · ∇eh k20 + kb · ∇eh k20 ) 1/2 ≤ 2δkbk∞ h−1 cinv keh k2mat + 2δ 2 kbk2∞ h−2 c2inv µ−1 eh k20 . 0 kµ

Assuming (5.1) gives δ 2 (eh,t , b · ∇eh,t ) ≤

1 1 keh k2mat + kµ1/2 eh k20 . 4 8

(5.17)

Inserting now (5.11) – (5.17) into (5.10) leads to 1 d 13 1 19 εδ d δ d 1/2 keh k2b + εk∇eh k20 + kµ1/2 eh k20 + keh k2mat + k∇eh k20 + kc eh k20 2 dt 32 2 32 2 dt 2 dt ≤ CkTtr k20 , where C = 4µ−1 0 + 16δ + 1/8 and the conditions on the stabilization parameter are given in (5.1). Redefining C, this inequality can be written in the form d keh k2b + εk∇eh k20 + kµ1/2 eh k20 + keh k2mat dt d d +εδ k∇eh k20 + δ kc1/2 eh k20 ≤ CkTtr k20 . dt dt 15

(5.18)

Integration in (0, t) yields keh (t)k2b + εk∇eh k2L2 (0,t;L2 ) + kµ1/2 eh k2L2 (0,t;L2 ) + keh k2L2 (0,t;mat) +εδk∇eh (t)k20 + δk(c1/2 eh )(t)k20 ≤ keh (0)k2b + εδk∇eh (0)k20 + δk(c1/2 eh )(0)k20 + C

t

Z

kTtr k20 dτ.

(5.19)

0

Now, the terms on the right hand side of (5.19) have to be bounded. It is keh (0)k2b + εδk∇eh (0)k20 + δk(c1/2 eh )(0)k20 17 εδc2inv 2 2 + δkck + ≤ ∞ keh (0)k0 ≤ Ckeh (0)k0 16 h2

(5.20)

since ε ≤ h is assumed. The next step of the error analysis uses that convection and reaction do not depend on time. Hence (5.2) can be differentiated with respect to time. Using (5.5), one obtains steady-state SUPG problems for Πh ut (t) with corresponding error estimates of type (5.4) kTtr (t)kSUPG ≤ Chr+1/2 kut (t)kr+1 =⇒ kTtr (t)k0 ≤ C

hr+1/2 1/2

µ0

kut (t)kr+1 , t ∈ [0, T ].

(5.21) Summarizing all constant into a generic constant, the following lemma is proven. Lemma 5.1. Error estimate for the time-continuous case involving the material derivative. Let t ≤ T < ∞ and let ut ∈ L2 (0; T ; H r+1 (Ω)). Then, the error eh = uh − Πh u satisfies 1/2 keh (t)kb + εk∇eh k2L2 (0,t;L2 ) + keh k2L2 (0,t;mat) + kµ1/2 eh k2L2 (0,t;L2 ) +δ 1/2 ε1/2 k∇eh (t)k0 + k(c1/2 eh )(t)k0 ≤ C keh (0)k0 + hr+1/2 kut kL2 (0,t;H r+1 ) , (5.22) where C depends on the coefficients of the problem and on cinv . An estimate for u − uh is now obtained by applying the triangle inequality and using (5.4) for estimating the terms with u − Πh u. In the second step, an estimate with the stronger SUPG norm δkb · ∇eh kL2 (0,t;L2 ) instead of keh kL2 (0,t;mat) is derived. To this end, insert once more vh = eh into the error equation (5.7) and apply a standard analysis by using the coercivity (2.12) 1 d 1 kTtr k20 kµ1/2 eh k20 keh k20 + keh k2SUPG ≤ + + 2δkTtr k20 2 dt 2 µ0 4 kb · ∇eh k20 +2δkeh,t k20 + δ . 4 The second and the last term can be absorbed into the left hand side. The first and the third are estimated by (5.21). Integration in (0, t) gives keh (t)k20 +keh k2L2 (0,t;SUPG) ≤ keh (0)k20 +C h2r+1 kut k2L2 (0,t;H r+1 ) + δkeh,t k2L2 (0,t;L2 ) . (5.23) 16

The estimate of the last term uses once more that convection and reaction are functions independent of time. Hence, (5.2) and (5.6) can be differentiated with respect to time, giving the same type of equations. Now, the error analysis for eh leading to (5.22) can be carried out in the same way for eh,t . One obtains δ δkeh,t k2L2 (0,t;L2 ) ≤ 0 kµ1/2 eh,t k2L2 (0,t;L2 ) ≤ Cδ keh,t (0)k20 + h2r+1 kutt k2L2 (0,t;H r+1 ) . µ (5.24) Next, δkeh,t (0)k20 has to be bounded in terms of eh (0) and Ttr (0) since it is not clear how to control eh,t (0) by an appropriate choice of uh (0). To this end, eh,t (t) is inserted into the error equation (5.7) leading to keh,t k20 = −aSUPG (eh , eh,t ) + (Ttr , eh,t + δb · ∇eh,t ) − δ(eh,t , b · ∇eh,t ).

(5.25)

Applying the Cauchy–Schwarz inequality and the inverse inequality, using the assumption ε ≤ h, and (5.1) yield aSUPG (uh , vh ) εδc2inv kbk∞ εcinv 1/2 k∇uh k0 + kb · ∇uh k0 + kck1/2 uh k0 + k∇uh k0 ≤ ∞ kc h h2 ! 1/2 δcinv kbk∞ kck∞ δcinv kbk∞ kb · ∇uh k0 + kc1/2 uh k0 kvh k0 + h h ≤ C k∇uh k0 + kb · ∇uh k0 + kc1/2 uh k0 kvh k0 , where C depends on cinv and the coefficients of the problem. Applying this estimate to (5.25) and using (5.1) give 1 keh,t k20 ≤ C k∇eh k0 + kb · ∇eh k0 + kc1/2 eh k0 + kTtr k0 keh,t k0 + keh,t k20 . 8 From this follows, with (5.21), δkeh,t k20 ≤ Cδ k∇eh k20 + kb · ∇eh k20 + kc1/2 eh k20 + h2r+1 kut (t)k2r+1 . Using this estimate for t = 0 in (5.24), inserting then (5.24) into (5.23), and applying the triangle inequality leads to the following error estimate. Theorem 5.2. Error estimate for the time-continuous case involving the SUPG norm. Let t ≤ T < ∞, let ut (t) ∈ H r+1 (Ω) for all t ∈ [0, T ], and u, ut , utt ∈ L2 (0; T ; H r+1 (Ω)). Then, the error estimate k(u − uh )(t)k0 + ku − uh kL2 (0,t;SUPG) h i ≤ C keh (0)k0 + δ 1/2 k∇eh (0)k0 + k(b · ∇eh )(0)k0 + k(c1/2 eh )(0)k0 h +Chr+1/2 ku(t)kr+1 + δ 1/2 kut (0)kr+1 + kukL2 (0,t;H r+1 ) + kut kL2 (0,t;H r+1 ) i +δ 1/2 kutt kL2 (0,t;H r+1 ) (5.26) holds. The constants depend on the coefficients of the problem and on cinv . Choosing the initial finite element solution uh (0) in the way that uh (0) solves aSUPG (uh (0), vh ) = (f (0) − ut (0), vh ) + δ(f (0) − ut (0), b · ∇vh ) = (−ε∆u0 + b · ∇u0 + cu0 , vh + δb · ∇vh ) ∀ vh ∈ Vh,r leads to eh (0) = 0 such that all terms with eh (0) vanish in (5.26). 17

5.2. Temporal discretization with the backward Euler scheme. The analysis of the backward Euler scheme (3.1) is based on the same ideas that were used in the time-continuous case. For this reason, only its most important steps are presented. Let enh = unh − Πh u(tn ), enh,τ = (enh − en−1 )/k for n ≥ 1, and h Πh u(tn ) − Πh u(tn−1 ) n Ttr,eul . = (ut (tn ) − Πh ut (tn )) + Πh ut (tn ) − k A straightforward calculation, using (5.3), gives the error equation n n (enh,τ , vh ) + aSUPG (enh , vh ) = (Ttr,eul , vh ) + δ(Ttr,eul , b · ∇vh ) − δ(enh,τ , b · ∇vh ) (5.27)

for all vh ∈ Vh,r . Taking first the inner product with enh and then with δenh,τ , one obtains (enh,τ , enh ) + δ 2 (b · ∇enh,τ , b · ∇enh ) + εk∇enh k20 + kµ1/2 enh k20 + kenh k2mat,eul +εδ(∇enh , ∇enh,τ ) + δ(cenh , enh,τ ) n n = (Ttr,eul , enh + δb · ∇enh ) + δ(Ttr,eul , (enh + δb · ∇enh ),τ )

(5.28)

−δ(cenh , b · ∇enh ) − δ 2 (cenh , b · ∇enh,τ ), X X + δε(∆enh , b · ∇enh )K + δ 2 ε(∆enh , b · ∇enh,τ )K − δ 2 (enh,τ , b · ∇enh,τ ), K∈Th

K∈Th

where the material derivative for the Euler scheme is defined by kvh kmat,eul := δ 1/2 kvh,τ + b · ∇vh k0 . The right hand side of (5.28) has principally the same form as the right hand side of (5.10). The second term of (5.28) can be estimated in an analogous way as (5.12), with the material derivative replaced by its discrete counterpart. Applying the same techniques as for the time-continuous case, one arrives at the discrete version of (5.18) (enh,τ , enh ) + δ 2 (b · ∇enh,τ , b · ∇enh ) + εk∇enh k20 + kµ1/2 enh k20 + kenh k2mat,eul n +εδ(∇enh , ∇enh,τ ) + δ(cenh , enh,τ ) ≤ CkTtr,eul k20 .

1 kenh k20 − ken−1 k20 + kenh − ehn−1 k20 and similar expressions Note that (enh,τ , enh ) = 2k h hold for all other inner products on the left hand side. Summation of the discrete times j = 1, . . . , n, and neglecting terms of the form kejh − ehj−1 k20 on the left hand side give kenh k2b +

n X k εk∇ejh k20 + kµ1/2 ejh k20 + kejh k2mat,eul + εδk∇enh k20 + δkc1/2 enh k20 j=1

≤ ke0h k2b + εδk∇e0h k20 + δkc1/2 e0h k20 + C

n X

j kkTtr,eul k20 .

(5.29)

j=1

As in the time-continuous case, the errors at the initial time can be bounded by ke0h k2b + εδk∇e0h k20 + δkc1/2 e0h k20 ≤ Cke0h k20 . 18

Using (5.4) and reasoning as in (4.3) yield n X

j kkTtr,eul k20

j=1

2 !

Πh u(tn ) − Πh u(tn−1 )

≤2 k kut (tn ) − + Πh ut (tn ) −

k 0 j=1 Z tn kΠh utt k20 dt ≤ C knh2r+1 max kut (ts )k2r+1 + k 2 (5.30) 0≤s≤n t0 ≤ C knh2r+1 kut k2L∞ (0,tn ;H r+1 ) + k 2 kΠh utt k2L2 (0,tn ;L2 ) . n X

Πh ut (tn )k20

Inserting the last two estimates into (5.29) gives an estimate for the material derivative of the error. The convergence analysis for the error in the SUPG norm starts by taking vh = enh in (5.27) and applying similar estimates as in the time-continuous case. Absorbing terms leads to an estimate of the form (5.23)   n n n X X X j j j kenh k20 + k keh k2SUPG ≤ C ke0h k20 + kkTtr,eul k20 + kδkeh,τ k20  . j=1

j=1

j=1

The second term on the right hand side is estimated in (5.30). Set zhn = enh,τ and assume for simplicity e0h = 0. Subtracting (5.27) for n − 1 from (5.27) for n, as n ≥ 2, and dividing by k give for n ≥ 2 n n n n (zh,τ , vh ) + aSUPG (zhn , vh ) = (Ttr,eul,τ , vh ) + δ(Ttr,eul,τ , b · ∇vh ) − δ(zh,τ , b · ∇vh ) n−1 n n for all vh ∈ Vh,r , where Ttr,eul,τ = (Ttr,eul − Ttr,eul )/k. Then, the same error analysis with respect to the material derivative as for enh can be applied to zhn , leading to the analog of (5.29)

µ0

n X

kδkzhj k20 ≤ µ0 kδkzh1 k20 +

j=1

n X

kδkµ1/2 zhj k20 ≤ Cδkzh1 k20 + C

j=2

n X

j kδkTtr,eul,τ k20 ,

j=2

where a contribution to the first term on the right hand side comes also from the error analysis for the material derivative. It is ut (tj ) − ut (tj−1 ) j − utt (tj ) + (utt (tj ) − Πh utt (tj )) Ttr,eul,τ = k Πh u(tj ) − 2Πh u(tj−1 ) + Πh u(tj−2 ) + Πh utt (tj ) − . k2 The first and second term can be bounded using the technique applied in (4.3) and the estimate (5.4). The third term is an O(k) approximation to Πh utt (tj ) such that an analog estimate to (4.3) can be performed. One obtains n X

h j kδkTtr,eul,τ k20 ≤ Cδ knh2r+1 max kutt (ts )k2r+1 0≤s≤n

j=2

+k 2 kuttt k2L2 (0,tn ;L2 ) + kΠh uttt k2L2 (0,tn ;L2 ) 19

i .

Now, it remains to bound δkzh1 k20 . To this end, the same approach can be used as in the time-continuous framework to bound keh,t (0)k0 . One gets 2 ε 1 2 1 2 1/2 1 2 1 2 δkzh1 k20 ≤ Cδ k∇e k + kb · ∇e k + kc e k + kT k h 0 h 0 h 0 tr,eul 0 h2 δkck∞ 1/2 1 2 δε 1 2 1 2 1 2 ≤C εk∇e k + δkb · ∇e k + kµ e k + δkT k h 0 h 0 h 0 tr,eul 0 . h2 µ0 In addition to the time-continuous case, the first three terms on the right hand side have to be estimated. Taking n = 1 and vh = e1h in (5.27), and recalling e0h = 0, standard estimates lead to ke1h k2SUPG ≤

1 1 2 1 ke k + ke1h k2SUPG ≤ CkTtr,eul k20 . k h 0

1 A direct estimate gives kTtr,eul k20 ≤ C(h2r+1 + k 2 ), where the constant depends on kΠh utt (t1 )k0 and kut (t1 )kr+1 . Now, the assumption ε ≤ h and the definition (5.1) of the stabilization parameter lead to δkzh1 k20 ≤ C h2r+1 + k 2 . Theorem 5.3. Error estimate for the backward Euler scheme involving the SUPG norm. Let tn = T < ∞, let u, ut ∈ L∞ (0, T ; H r+1 (Ω)) and Πh utt , uttt , Πh uttt ∈ L2 (0; T ; L2 (Ω)). Then, the error estimate

k(u − uh )(tn )k20 + k

n X

k(u − uh )(tj )k2SUPG ≤ C h2r+1 + k 2

(5.31)

j=1

holds with the stabilization parameter defined in (5.1). The constant in (5.31) depends on norms of u and Πh u in the spaces given above, on the coefficients of the problem (2.1) and on cinv . Applying the triangle inequality and (5.4), the constant can be expressed in norms of time derivatives of u instead of Πh u. 5.3. Temporal discretization with the Crank–Nicolson scheme. The Crank– Nicolson/SUPG scheme has the form: For n = 1, 2, . . . find Uhn ∈ Vh,r such that n n Uh + Uhn−1 f + f n−1 ,ϕ = k ,ϕ (Uhn − Uhn−1 , ϕ) + kaSUPG 2 2 n X X f + f n−1 + δK , b · ∇ϕ − δK (Uhn − Uhn−1 , b · ∇ϕ)K . 2 K K∈Th

K∈Th

The error analysis for this method is to the most part similar to that of the backward Euler scheme. Using the Galerkin orthogonality (5.3), one can derive the error equation n n n (enh,τ , vh ) + aSUPG (en,∗ h , vh ) = (Ttr,CN , vh ) + δ(Ttr,CN , b · ∇vh ) − δ(eh,τ , b · ∇vh ) (5.32) n−1 n for all vh ∈ Vh,r , where en,∗ )/2, and the truncation error is h = (eh + eh n Ttr,CN =

1 ((ut (tn ) − Πh ut (tn )) + (ut (tn−1 ) − Πh ut (tn−1 ))) 2 1 Πh u(tn ) − Πh u(tn−1 ) + (Πh ut (tn ) + Πh ut (tn−1 )) − . 2 k 20

Taking first vh = en,∗ in (5.32), then vh = δenh,τ , and adding both equations yield h n,∗ n,∗ 2 2 n 1/2 n,∗ 2 2 (enh,τ , en,∗ eh k0 + ken,∗ h ) + δ (b · ∇eh,τ , b · ∇eh ) + εk∇eh k0 + kµ h kmat,CN n,∗ n n +εδ(∇en,∗ h , ∇eh,τ ) + δ(ceh , eh,τ ) n,∗ n n n n = (Ttr,CN , en,∗ h + δb · ∇eh ) + δ(Ttr,CN , (eh + δb · ∇eh )τ ) n,∗ n,∗ 2 n −δ(cen,∗ h , b · ∇eh ) − δ (ceh , b · ∇eh,τ ), X X n,∗ n 2 n n + δε(∆en,∗ δ 2 ε(∆en,∗ h , b · ∇eh )K + h , b · ∇eh,τ )K − δ (eh,τ , b · ∇eh,τ ), K∈Th

K∈Th

1/2 n where ken,∗ keh,τ + b · ∇en,∗ h kmat,CN := δ h k0 . This equation has the same form as n (5.28), with en,∗ instead of e . Thus, using analogous estimates gives h h n,∗ n,∗ 2 2 n 1/2 n,∗ 2 2 (enh,τ , en,∗ eh k0 + ken,∗ h ) + δ (b · ∇eh,τ , b · ∇eh ) + εk∇eh k0 + kµ h kmat,CN n,∗ n n n 2 +εδ(∇en,∗ h , ∇eh,τ ) + δ(ceh , eh,τ ) ≤ CkTtr,CN k0 .

Note that (∇enh,τ , ∇en,∗ h )=

1 k∇enh k20 − k∇en−1 k20 h 2k

(5.33)

and similarly for the other inner products. Summation of the time steps j = 1, . . . , n gives kenh k2b +

n X j,∗ 2 2 1/2 j,∗ 2 n 2 1/2 n 2 k εk∇ej,∗ k + kµ e k + ke k eh k0 0 0 mat,CN + εδk∇eh k0 + δkc h h h j=1

≤ ke0h k2b + εδk∇e0h k20 + δkc1/2 e0h k20 + C

n X

j kkTtr,CN k20 .

(5.34)

j=1

The initial error is bounded as in (5.20). For the truncation error, the first part is estimated with (5.4) 1 2 k(ut (tn ) − Πh ut (tn )) + (ut (tn−1 ) − Πh ut (tn−1 ))k0 2 ≤ Ch2r+1 kut (tn )k2r+1 + kut (tn−1 )k2r+1 . For the second term, one obtains with a direct calculation, using integration by parts, and with the Cauchy–Schwarz inequality

Πh ut (tn ) + Πh ut (tn−1 ) Πh u(tn ) − Πh u(tn−1 ) 2

−

2 k 0

Z

2

tn 1 1

(tn − t)(t − tn−1 )Πh uttt dt = 2

k tn−1 2 0  ! !1/2 2 1/2 Z tn Z tn 1  1  ≤ 2 (tn − t)2 (t − tn−1 )2 dt kΠh uttt k20 dt k tn−1 4 tn−1 Z tn 3 kΠh uttt k20 dt. ≤ Ck tn−1

21

Altogether, one gets n X

j kkTtr,CN k20 ≤ Cknh2r+1 max kut (ts )k2r+1 + Ck 4 kΠh uttt k2L2 (0,tn ;L2 ) . 0≤s≤n

j=1

(5.35)

Inserting this estimate into (5.34) leads to an error estimate involving the material derivative. Setting vh = en,∗ in (5.32) gives with the coercivity (2.12), a relation of form h (5.33), and the absorption of terms similar to the time-continuous case,   n n n X X X j 2 ke0h k20 + kenh k20 + k kej,∗ kkTtr,CN k20 + kδkejh,τ k20  . h kSUPG ≤ C j=1

j=1

j=1

The second term of the right hand side is estimated in (5.35) and it remains the estimate of the third term. Setting zhn = enh,τ and assuming for simplicity e0h = 0 leads analogously to the analysis of the backward Euler scheme to n n n n (zh,τ , vh ) + aSUPG (zhn,∗ , vh ) = (Ttr,CN,τ , vh ) + δ(Ttr,CN,τ , b · ∇vh ) − δ(zh,τ , b · ∇vh ) n−1 n n − Ttr,CN )/k and n ≥ 2. In the same way as = (Ttr,CN for all vh ∈ Vh,r , where Ttr,CN,τ n n for eh , an error estimate for zh is derived. In contrast to the backward Euler scheme, the terms kµ1/2 zhj,∗ k20 are on the left hand side of this estimate and not kµ1/2 zhj k20 , see (5.34). Hence, different terms are needed to continue the estimate. Using the equivalence of k · kb and the L2 norm, and T = kn give   n n n X X X j 2 j j kδkzh k0 ≤ kδkzh1 k20 + C kδkzh k2b ≤ CT δkzh1 k20 + kδkTtr,CN,τ k20  . j=1

j=2

j=2

The estimate of the terms in the last line proceeds in the same manner as for the backward Euler scheme. It is ut (tn ) − ut (tn−2 ) n − utt (tn−1 ) + (utt (tn−1 ) − Πh utt (tn−1 )) Ttr,CN,τ = 2k Πh u(tn ) − 2Πh u(tn−1 ) + Πh u(tn−2 ) + Πh utt (tn−1 ) − . k2 Similarly to (4.3), it can be shown that the first and third term are of order k 2 and that the factor depends on kutttt kL2 (tn−2 ,tn ;L2 ) and kΠh utttt kL2 (tn−2 ,tn ;L2 ) , respectively. The second term is estimated with (5.4) to be of order hr+1/2 and the constant depends on kutt kL2 (tn−2 ,tn ;H r+1 ) . Finally, δkzh1 k20 ≤ C(h2r+1 + k 4 ) is proven analogously as for the backward Euler scheme, where C depends on kΠh uttt (t1 )k0 and kut (t1 )kr+1 . Theorem 5.4. Error estimate for Crank–Nicolson scheme involving the SUPG norm. Let tn ≤ T < ∞, let u, ut ∈ L∞ (0, T ; H r+1 (Ω)), utt ∈ L2 (0, T ; H r+1 (Ω)), and Πh uttt , utttt , Πh utttt ∈ L2 (0; T ; L2 (Ω)). Then, the error estimate

n X

u(tj ) + u(tj−1 ) uh (tj ) + uh (tj−1 ) 2

k(u − uh )(tn )k20 + k − ≤ C h2r+1 + k 4

2 2 SUPG

j=1

(5.36) holds with the stabilization parameter defined in (5.1). The constant in (5.36) depends on the coefficients of the problem (2.1), on norms of u and Πh u in the spaces given above, on cinv , and linearly on T . 22

6. Numerical studies. Two examples will be presented in the numerical studies. The first one, possessing a given smooth solution, serves as support for the orders of convergence that were proven in the previous sections. The second example is the well-known rotating body problem from [14]. It demonstrates the superiority of the parameter choice from Section 5 compared with the choices from Sections 3 and 4 for small time steps on a fixed, rather coarse, spatial mesh. Example 6.1. Smooth solution. Consider (2.3) with Ω = (0, 1)2 , T = 1, different values of ε, b = (1, −1), c = 1, and the right-hand side was chosen such that u(t, x, y) = esin(2πt) sin(2πx) sin(2πy) is the solution of (2.3). The simulations were performed with ε = 10−8 in the convection-dominated regime and with ε = 1 in the diffusion-dominated regime. Uniform triangular grids were used with the coarsest grid (level 0) obtained by dividing the unit square with a diagonal from (0, 0) to √ (1, 1). The mesh width h was defined by dividing the diameters of the mesh cells by 2. To prevent superconvergence, the convection field was chosen such that it is not parallel to any grid line. Consider the error estimates (4.4) and (4.5). First, optimal scalings of the mesh width h and the length of the time step k were derived from these estimates. Then, the error estimates lead to only one asymptotic order of convergence that serves as criterion. The stabilization parameters for the estimates under the assumptions of Theorem 4.1 were set to be δK = δ = k/4, according to condition (3.3). In the convectiondominated regime, ε h, the terms O(k) and O(hr+1 δ −1/2 ) = O(hr+1 k −1/2 ) have to be balanced to obtain an optimal L2 -error estimate (4.4). This leads to the scaling k = O(h2(r+1)/3 ). The same reasoning applies for the SUPG error (4.5). If the final time T = 1 was not obtained exactly with the chosen time steps, the simulations were stopped at the first discrete time larger than T . In the diffusion-dominated regime, h ≤ ε, the terms O(k),O(k 1/2 hr−1 ε), and O(hr+1 k −1/2 ) need to be balanced. This leads to k = O(h2(r+1)/3 ) or k = (h2 /ε). If h ε, the second scaling gives a better order of convergence for r = 1 (piecewise linear elements). Note, in this case, δ = k = O(h2 /ε) is a standard choice of the stabilization parameter in the diffusiondominated regime for steady-state problems. For r = 2, both scalings are essentially the same and for r ≥ 3, the first scaling leads to a higher order of convergence. For the SUPG estimate, the same terms have to be balanced. In addition, the order of convergence is bounded by the term O(ε1/2 hr ), such that for h ε only first order convergence can be expected for r = 1. Figure 6.1 presents the orders of convergence for the P1 , P2 , and P3 finite element. It can be seen that all orders match the predictions from the analysis. √ For the estimates which can be obtained with δK = khK /(4kbk2 ), one can also observe the expected orders of convergence. For brevity, this is not presented in detail. Next, estimate (5.26) for the time-continuous case is considered, from which one can expect convergence for the L2 norm and the L2 (0, T ; SUPG) norm of order r + 1/2 for Pr finite elements if the initial condition is sufficiently well approximated. The simulations were performed with a very small time step, k = 10−6 . As initial condition, the Lagrange interpolant of u(0, x, y) was used. The results are presented in Figure 6.2. The observed order of convergence in the L2 norm is even higher than the prediction from the analysis. As example for the order of convergence obtained for a fully discrete scheme, estimate (5.36) of the Crank–Nicolson scheme is considered. The equilibration of the terms on the right hand side of (5.36) gives the length of the time step k = O(hr/2+1/4 ) 23

Fig. 6.1. Example 6.1, orders of convergence for the estimates (4.4) and (4.5); left: convectiondominated regime; right: diffusion-dominated regime.

and the expected order of convergence is r + 1/2. This order can be well observed in the numerical results in the right picture of Figure 6.2.

Fig. 6.2. Example 6.1, left: orders of convergence for the estimate (5.26), right: orders of convergence for the estimate (5.36) for k = hr/2+1/4 , both convection-dominated regime.

Example 6.2. Rotating body problem. This problem was studied numerically for finite element discretizations of convection-diffusion equations already in [12]. Here, exactly the same setting is used. The aim of this example is to illustrate that the choice of the stabilization parameter δK = O(hK ) from Section 5 is much better than the choices δK = O(k), δK = O(k 1/2 hK ) from Sections 3 and 4 in the presence of very small time steps. Let Ω = (0, 1)2 , ε = 10−20 , b = (0.5 − y, x − 0.5)T , and c = f = 0. The initial condition, consisting of three disjoint bodies. is presented in Figure 6.3, see [12] for a precise definition of the bodies. The rotation of the bodies occurs counter-clockwise. A full revolution takes t = 6.28 ≈ 2π. With the extremely small diffusion, the solution after one revolution is essentially the same as the initial condition. Homogeneous Dirichlet boundary conditions were imposed. In the simulations, a uniform grid consisting of 128 × 128 triangles was used. This leads to 16 641 degrees of freedom for the P1 finite element method, including Dirichlet nodes. For brevity, only results obtained with the backward Euler scheme 24

Fig. 6.3. Example 6.2, initial condition and ideal solution after one rotation.

will be presented in detail. Computational studies were performed for the Galerkin finite element method (δK = 0 for all mesh cells), the choice of the stabilization parameter from [12, formulae (8)√and (11)], that results in δK = k for the backward Euler scheme, the choice δK = khK /4 and δK = hK /4. Analogously to [12], a measure for the under- and overshoots is given by var(t) := max uh (t; x, y) − min uh (t; x, y), (x,y)∈Ω

(x,y)∈Ω

with the optimal value var(t) = 1 for all t. The backward Euler scheme with k = 10−6 is considered. It can be observed that by far the best result was obtained with δK = hK /4, see Figures 6.4 and 6.5. However, the computed solution with these parameters possesses still non-negligible spurious oscillations. Using the stabilization parameters from the analysis of Sections 3 and 4 leads for very small time steps on a fixed spatial grid to similar results as for the Galerkin finite element method. Only a slight damping of the spurious oscillations can be observed, see the ranges of the finite element solutions in Figure 6.5.

Fig. 6.4. Example 6.2, under- and overshoots measured with var(t), backward Euler scheme, k = 10−6 .

7. Summary and Outlook. This paper studied different ways to obtain error estimates for the SUPG finite element method applied to evolutionary convectiondiffusion-reaction equations. 25

Fig. 6.5. Example 6.2, computed solutions after one revolution with the backward Euler√ scheme and k = 10−6 : Galerkin finite element method, SUPG with δK = O(k), SUPG with δK = O( khK ), SUPG with δK = O(hK ); left to right, top to bottom.

Exemplarily for the backward Euler scheme, it was shown that standard energy arguments for the fully discrete problem yield error estimates under conditions that couple the choice of the stabilization parameters to the length of the time step. In particular, the SUPG stabilization vanishes in the time-continuous limit. Numerical evidence shows that this is not the correct behavior. For this reason, the special case of convection and reaction being independent of time was considered. On uniform grids, optimal error estimates could be proven with the choice of the stabilization parameter δ = O(h) in the convection-dominated regime. The proofs were given for the time-continuous case and the fully discrete method with the backward Euler scheme as well as with the Crank–Nicolson scheme. The analysis of the general time-continuous problem, with time-dependent coefficients and using non-uniform meshes, is open. An extension of the analysis from Section 5 seems to be hard, since this analysis uses several times that the original equation can be differentiated with respect to time yielding essentially the same equation. To answer the question, if the general time-continuous problem allows at all estimates with stabilization parameters independent of the length of the time step, requires further research. With respect to the usage of the SUPG finite element method in time-dependent convection-diffusion-reaction equations, the results of Section 5, Example 6.2 and other numerical studies from the literature strongly suggest to define the stabilization parameters in the convection-dominated regime in the classical way by δK = O(hK ), 26

despite the open analysis for the general time-continuous problem. REFERENCES [1] P.B. Bochev, M.D. Gunzburger, and J.N. Shadid, Stability of the supg finite element method for transient advection-diffusion problems, Comput. Methods Appl. Mech. Engrg., 193 (2004), pp. 2301 – 2323. [2] A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations., Comput. Methods Appl. Mech. Engrg., 32 (1982), pp. 199 – 259. [3] E. Burman, Consistent SUPG-method for transient transport problems: Stability and convergence, Comput. Methods Appl. Mech. Engrg., 199 (2010), pp. 1114 – 1123. [4] P.G. Ciarlet, The finite element method for elliptic problems, North-Holland Publishing Company, Amsterdam – New York – Oxford, 1978. , Basic error estimates for elliptic problems, in Handbook of Numerical Analysis II, P.G. [5] Ciarlet and J.L. Lions, eds., North-Holland Amsterdam, New York, Oxford, Tokyo, 1991, pp. 19 – 351. [6] R. Codina, Comparison of some finite element methods for solving the diffusion-convectionreaction equation, Comput. Methods Appl. Mech. Engrg., 156 (1998), pp. 185 – 210. [7] J.L. Guermond, Subgrid stabilization of Galerkin approximations of linear contraction semigroups of class C 0 in Hilbert spaces, Numer. Methods Partial Differential Equations, 17 (2001), pp. 1 – 25. [8] M. C. Hsu, Y. Bazilevs, V. M. Calo, T. E. Tezduyar, and T. J. R. Hughes, Improving stability of stabilized and multiscale formulations in flow simulations at small time steps, Comput. Methods Appl. Mech. Engrg., 199 (2010), pp. 828 – 840. [9] T.J.R. Hughes and A.N. Brooks, A multidimensional upwind scheme with no crosswind diffusion, in Finite Element Methods for Convection Dominated Flows, AMD vol.34, T.J.R. Hughes, ed., ASME, New York, 1979, pp. 19 – 35. [10] T.J.R. Hughes, L.P. Franca, and M. Mallet, A new finite element formulation for computational fluid dynamics: VI. Convergence analysis of th3 generalized SUPG formulation for linear time-dependent multidimensional advective-diffusive systems, Comput. Methods Appl. Mech. Engrg., 63 (1987), pp. 97 – 112. [11] V. John, T. Mitkova, M. Roland, K. Sundmacher, L. Tobiska, and A. Voigt, Simulations of population balance systems with one internal coordinate using finite element methods, Chem. Engrg. Sci., 64 (2009), pp. 733 – 741. [12] V. John and E. Schmeyer, Stabilized finite element methods for time–dependent convection– diffusion–reaction equations, Comput. Methods Appl. Mech. Engrg., 198 (2008), pp. 475 – 494. [13] , On finite element methods for 3d time–dependent convection–diffusion–reaction equations with small diffusion, in BAIL 2008 – Boundary and Interior Layers, vol. 69 of Lecture Notes in Computational Science and Engineering, Springer, 2009, pp. 173 – 182. [14] R.J. LeVeque, High–resolution conservative algorithms for advection in incompressible flow, SIAM J. Numer. Anal., 33 (1996), pp. 627 – 665. [15] G. Lube and D. Weiss, Stabilized finite element methods for singularly perturbed parabolic problems, Appl. Numer. Math., 17 (1995), pp. 431 – 459. [16] H.-G. Roos, M. Stynes, and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, vol. 24 of Springer Series in Computational Mathematics, Springer, 2nd ed., 2008. [17] M. Stynes, Steady-state convection-diffusion problems, in Acta Numerica, A. Iserles, ed., Cambridge University Press, 2005, pp. 445 – 508. [18] L.B. Wahlbin, Superconvergence in Galerkin Finite Element Methods, no. 1605 in Lecture Notes in Math., Springer-Verlag, Berlin, 1975.

27