An Adaptive Algorithm for the Time Dependent Transport Equation with Anisotropic Finite Elements and the Crank–Nicolson Scheme

Anisotropic finite elements and the Crank–Nicolson scheme are considered to solve the time dependent transport equation. Anisotropic a priori and a posteriori error estimates are derived. The sharpness of the error indicator is studied on non-adapted meshes and time steps. An adaptive algorithm in space and time is then designed to control the error at final time. Numerical results show the accuracy of the method.


Introduction
Anisotropic adaptive meshes with large aspect ratio have proved to be extremely efficient for partial differential equations with free boundaries or boundary layers, see for instance [1][2][3] for applications in computational fluid dynamics. In most cases, the adaptive criteria are based on heuristics or interpolation error estimates rather than rigorous a posteriori error estimates. This is particularly the case when the time dependent transport equation is involved, since few a posteriori error estimates are available [4].
In [5], anisotropic a posteriori error estimates were derived for the time dependent transport equation discretized in space only. An adaptive algorithm with meshes having large aspect ratio was proposed. Since the error due to time discretization was not considered, small time steps were used. In this paper, the error due to time discretization is taken into account. More precisely, the order two Crank-Nicolson scheme is used and an appropriate piecewise quadratic time reconstruction is advocated, as in [6]. The quality of our error estimator is first validated on non-adapted meshes and constant time steps. An adaptive algorithm is then B Marco Picasso marco.picasso@epfl.ch Samuel Dubuis samuel.dubuis@epfl.ch 1 Institute of Mathematics, Station 8, EPFL, 1015 Lausanne, Switzerland proposed, with goal to build a sequence of anisotropic meshes and time steps, so that the final error is close to a preset tolerance. Numerical results on adapted, anisotropic meshes and time steps show the efficiency of the method.
It is well known that the classical Galerkin formulation is unsuitable for the transport equation and that some stabilization techniques are necessary. Assume that Ω is a polygon and Γ − is the union of edges lying on ∂Ω. For any h > 0, let T h be a conformal triangulation ofΩ into triangles K of diameter h K less than h. Let V h be the set of continuous piecewise linear functions on each triangle of T h , with zero value on Γ − . A possible finite element discretization in space is to search for u h : Ω × (0, T ) −→ R such that u h (·, 0) = r h u 0 (r h is the Lagrange interpolant) and, for all 0 ≤ t ≤ T where δ h > 0 is a stabilization parameter that will be specified later on. A numerical study of (2) with anisotropic finite elements has already been proposed in [5], our goal is to take into account an order two in time discretization, namely the Crank-Nicolson scheme. Let N be a non-negative integer and consider a partition 0 = t 0 < t 1 < · · · < t N = T. We denote by τ n+1 = t n+1 − t n the time step, n = 1, 2, . . . , N − 1. Starting from u 0 h = r h u 0 , for n = 0, 1, . . . , N − 1, we are looking for u n+1 h ∈ V h such that Ω ∂u n+1 h + β · ∇u where ∂u n+1 We will only consider one mesh T h for the theoretical analysis of the scheme (3). Comments are added in the case of dynamic meshes in Sect. 4.2.

Anisotropic Finite Elements
In this paper, anisotropic finite elements will be used, that is to say meshes with possibly large aspect ratio. We will use the notations and results of [8][9][10], see also [11] for similar results. Let K ∈ T h and T K :K −→ K be the affine transformation mapping the reference triangleK into K defined by where R K and P K are orthogonal matrices and We note where r 1,K , r 2,K are the unit vectors corresponding to directions of maximum and minimum stretching, respectively, so that λ 1,K , λ 2,K are the value of maximum and minimum stretching. With these notations, the following interpolation results holds for the Lagrange interpolant r h [8,9]: where C > 0 is a constant depending only on the reference triangleK , and where H (v) is the Hessian matrix defined by The considerations that follow require also the use of Clément's interpolant. Since anisotropic meshes are considered, we assume that each vertex has a number of neighbours bounded from above, uniformly with respect to h. Moreover, we suppose that for each K , the diameter of where K is the union of triangles sharing a vertex with K , is uniformly bounded, independently of the mesh geometry. For more details, we refer again to [8][9][10]. In this framework, the following estimation holds where R h is the Clément's interpolant, C > 0 is a constant depending only on the reference triangleK , and with

A Priori Error Estimates
We now prove that the solution of the numerical method (3) converges to that of problem (1) for anisotropic meshes. This has already been proved for isotropic meshes in [12]. The stabilization parameter is kept constant in time and space. The key ingredient of the proof consists in taking test functions of the form v + δ h ∂v ∂t [12] and to observe that, since β is Since P 1 finite elements and the Crank-Nicolson method are used, it is expected that the error at final time in the isotropic settings. Theorem 1 Assume that β is not identically zero on Ω. Let u be the solution of (1) and let u N h be the solution of (3) with a constant δ h defined by Let Assume that the data T , Ω, f , β, u 0 are such that u ∈ H 1 (0, T ; H 2 (Ω)) and ∂ 3 u h ∂t 3 ∈ L 2 (0, T ; L 2 (Ω)). Let e(t n ) = u(t n ) − u n h , n = 0, . . . , N . Then, there exists C > 0 independent of the data T , Ω, f , β, u 0 , the mesh size, aspect ratio and the time step such that Remark 1 In the case of isotropic meshes, λ 1,K λ 2,K h K and L 2 K (u) ≤ Ch 4 K |u| 2 H 2 (K ) , where C is independent of the mesh size but can depend on the mesh aspect ratio. Thus, in these settings, (11) reduces to where h.o.t. stands for higher order terms.
Remark 2 As already explained in [5], estimate (11) is optimal with respect to the space discretization parameter for anisotropic meshes. Indeed, assume that the solution u depends only on one variable and that the mesh is aligned with the solution, then the estimate (11) reduces to and max K ∈T h λ 2,K → 0 is sufficient to ensure the convergence of the numerical method.

Remark 3
We have not been able to prove that dt is bounded independently of h and τ . The proof is not obvious, even for parabolic problems [13], and out of the scope of the present paper. It should be noticed that an a priori error estimate can also be proved introducing the anisotropic equivalent of the hyperbolic projector used in [12]. In this case, only derivatives of the exact solutions , appear in the error bound instead of dt. We have completed the proof, which is not presented here since it is significantly longer than the one below.
Remark 4 As in [12], a similar analysis can be performed if div β = 0 under restrictions on h and τ , with the price to pay that all the constants involved depends exponentially on the final time and the divergence of β.
where u h (t) is the solution of (2). Then we apply Theorem 3.1 of [5] on I 1 and we obtain where C > 0 depends on the reference triangleK only. In particular, C is independent of Ω, f, β, u 0 , u, T, N , the mesh size and aspect ratio and the time step.
We now have to estimate I 2 . By using several times the Fundamental Theorem of Calculus, one can derive that where In particular, we observe that In the sequel, we will note e n h = u h (t n )−u n h . By using (2), (3) and (13), the following relation holds for the numerical error Choosing and using (9), we therefore obtain Using Cauchy-Schwarz and Young's inequality yields Multiplication by 2τ n+1 and use of Cauchy-Schwarz, triangle and Young's inequalities yield after summing from 0 to N − 1 Here we use the fact that e 0 h = 0 and we have set τ 0 = τ N +1 = 0. Finally, we use the discrete Gronwall's Lemma (see [14], Lemma 5.1) and we get μ n 1−μ n ≤ 1 and using (14), we obtain Estimates (16) and (12) together yield the result.

A Posteriori Error Estimate
We now prove an a posteriori error estimate involving time and space discretization for problem (3). As in [5], the following choice for the stabilization parameter δ h is advocated.
else δ h|K is set to zero. As proposed in [6], we introduce a piecewise quadratic reconstruction of the computed solution in order to recover an O(τ 2 ) error estimator. We shall use the following notation x ∈Ω, n ≥ 1.
and for n = 0, Observe that (18) is a Newton polynomial; for every n ≥ 1, u hτ is the unique quadratic polynomial in time that equals u n−1 h , u n h , u n+1 h , at time t n−1 , t n , t n+1 , respectively. We first prove the following lemma : and Proof We start with n = 0. Then The result is then obtained noticing that We then take the difference of (3) with superscript n and (3) with superscript n − 1 to obtain Inserting into (23) yields the result.
We are now ready to prove our a posteriori error estimate.
Proof Let t ∈ t n , t n+1 , n ≥ 1. Using (9), (1), (3) and Lemma 1, we have The triangle and Cauchy-Schwarz inequalities imply Choosing v h = R h e, using estimation (7) and definition of δ h|K , we have see [5] for details. Therefore we have where we have set and where C denotes a positive constant, independent of T , Ω, f , β, u 0 , the mesh size, aspect ratio and the time step, which may change from line to line. Using the discrete Cauchy-Schwarz and Young's inequalities we therefore obtain where ε is any positive number. Multiplying by 2e −εt and integrating between time t 1 and T , we get Finally, we choose ε = 1 T so that the exponential growth in time is eliminated: In order to estimate e(t 1 ) 2 L 2 (Ω , we proceed in the same manner to obtain The desired estimate is obtained plugging (28) into (27). (24) is not a standard a posteriori estimate since the exact solution u is contained in ω K (e). However, post-processing techniques can be applied in order to approximate G K (e), for instance Zienkiewicz−Zhu (ZZ) post-processing. More precisely, we will replace the first order partial derivatives with respect to x i

Remark 5 Estimate
where, for any v h ∈ V h , for any vertex P of the mesh Numerical results already presented in [5,6,[15][16][17] showed the efficiency of ZZ post-processing for anisotropic meshes for elliptic, parabolic, and hyperbolic problems.

Remark 6
The following a posteriori error estimate can also be proved. Starting form (25), we have Cauchy-Schwarz inequality implies that which yields Estimate (29) was already pointed out in [4] and is valid for non-smooth solutions. Numerical experiments (not reported here) have shown that (29) is suboptimal for smooth solutions, thus estimate (24) should be preferred for smooth solutions.

Remark 7
We have not been able to prove a lower bound corresponding to estimate (24), this being also the case for parabolic problems with anisotropic finite elements [6]. However, for elliptic problems [16], we have been able to prove a lower bound provided that that is to say that the error is equidistributed in both directions r 1,K , r 2,K . (24) can be generalized in the case div β = 0 under the assumptions of the Theorem 2. In this case, the constant involved in (24) depends exponentially on the final time T and div β L ∞ (Ω) . Indeed, (25) becomes

Remark 8 Estimate
We conclude the proof using the same techniques as in Theorem 2, using the Gronwall's Lemma to control 1 2 Ω (div β)e 2 dx. Therefore, (24) becomes where C and c n are as in Theorem (2).

A Posteriori Error Indicators
We now define our error indicator Here the anisotropic error indicator in space η A is defined by

The error indicator in time is defined by
with θ given in Lemma 1 and c n as in Theorem 2. The reader should note that the other terms in (24) have not been considered since they are of higher order. In order to check the sharpness of these error indicators, we will compare them to the true errors. To this end, we introduce the effectivity indices ei and ei Z Z defined by Here ei measures the sharpness of our space-time error indicator, whereas ei Z Z measures the quality of our Zienkiewicz−Zhu post-processing.

Numerical Experiments on Non-adapted Meshes with Constant Time Steps
We now investigate the sharpness of our indicators by performing numerical experiments on nonadapted meshes with constant time steps. Problem (1) is considered in the unit square Ω = (0, 1) 2 , with T = 0.5, f = 0, β = (1, 0) T , the initial condition is given by the smooth function Γ − is the left boundary of Ω, thus the exact solution u(x 1 , x 2 , t) is given by The solution is smooth with small variations, except in a thin layer of width controlled by C, the larger C, the smaller the layer, the larger the error for a given mesh size. Several experiments have been performed on anisotropic meshes with aspect ratio varying from 50 to 500, where we keep the time step constant. In what follows, h 1 -h 2 denotes the mesh size in the directions x 1 , x 2 and τ is the time step. We first investigate the sharpness of the anisotropic error indicator in space η A , choosing τ = O(h 2 ) so that the error due to time discretization is negligible, see Tables 1 and 2. It is observed that the L 2 (Ω) error at final time is O(h 1.8 ) while the L 2 (0, T, H 1 (Ω)) error is O(h). The post-processed ZZ gradient is asymptotically exact, while the effectivity index ei converges to a value close to 20. These results agree with those of [5].  Table 3 Convergence results when h = O(τ 2 ) with C = 60 and aspect ratio 50 (rows 1-4) and 500 (rows 5-7) Then, we check that the quadratic reconstruction in (18) and (19) yields an error indicator of optimal second order in time. We choose h O(τ 2 ) so that the error due to the space discretization is negligible. The numerical results presented in the Tables 3 and 4 show that both the L 2 (Ω) error at final time and the time indicator η T are O(τ 2 ). The effectivity index tends to a value close to 2. Note that in this case, ei Z Z is away from 1, which implies that the post-processing included in our error indicator in space η A is not accurate; but this is unimportant since η A is much smaller than the error indicator in time η T .
In order to check that the effectivity index does not depend on Ω and T , we reproduce the same experiment on a domain Ω = (0, 10) × (0, 1) for several values of the final time Table 4 Convergence results when h = O(τ 2 ) with C = 240 and aspect ratio 50 (rows 1-5) and 500 (rows 6-8) T . The corresponding results are presented in Tables 5 and 6 for C = 60 and meshes with aspect ratio 50. The effectivity index remains close to the values obtained previously. In order to obtain an effectivity index close to one, we divide the space indicator η A by 20 and the time indicator η T by 2. We report the result obtain in Tables 7 and 8 where we consider the normalized error indicator The corresponding effectivity index is shown to be near a value of 1 when h 3 /2 = O(τ 2 ).

An Adaptive Algorithm in Space and Time
Although the analysis in Sect. 3 is restricted to a single mesh T h , we now present an adaptive space-time algorithm which involves several meshes. Then the question of interpolation between meshes is discussed. The goal of the adaptive space-time algorithm is to control e(T ) L 2 (Ω) . Given a prescribed tolerance TOL, we want to ensure that A sufficient condition is to ensure that, for n = 0, 1, 2, . . . , N − 1 and 0.75 2 TOL 2 τ n+1 2 ≤ The main steps of the adaptive algorithm are summarized in Fig. 1. At each time step, a new mesh is built, whenever needed. Then, the previous finite element approximation, u n h , has to be interpolated in order to compute the current one, u n+1 h . More precisely, if we denote Fig. 1 Adaptive algorithm. The index i denotes the number of remeshing required to build an acceptable mesh at current time, starting from the mesh accepted at previous time by T n h,i and T n h,i+1 two successive meshes generated at time t n+1 , and by V n h,i , V n h,i+1 the associated finite elements spaces, we consider the interpolation operator If a new mesh has to be built, then we interpolate the values of u n h from V n h,i to V n h,i+1 and compute u n+1 for all v h ∈ V n h,i+1 Five interpolation operators have been considered. -the linear Lagrange interpolation, -the exact L 2 projection [5,18], -the conservative algorithm of [19], -the Ritz hyperbolic projection [12], -the modified hyperbolic projection defined below.
We give more details on the last choice. For g ∈ H 1 (Ω), we define π n h,i+1 : The projection π n h,i+1 clearly satisfies the following property in (34). Using (9) yields We conclude by multiplying on each side by 2τ n+1 and using (36).
This stability result has a little interest in practice since δ h is not constant for adapted meshes.
In the numerical experiments, the best results have been obtained using the conservative algorithm of [19], the other four interpolation operators are shown to be less accurate. The BL2D software [20] is used in order to build anisotropic meshes, the indicator (η A K ) 2 being equidistributed in the directions of maximum and minimum stretching r 1,K , r 2,K . Each triangle K is aligned with the eigenvectors of the error gradient matrix G K (e), where ZZ post-processing is used in order to approximate ∂e/∂ x i . We shortly describe this remeshing procedure. Since the BL2D software uses informations about vertices, we need to translate the error indicator η A from triangles to vertices. We define, for all K ∈ T h , the anisotropic error indicator in direction r i,K , by and for all vertex P Then, a sufficient condition to ensure (32) is the following. For all vertex P of the mesh, η A,i P,n should satisfy Hereabove, the factor 3 is due to the fact that summing over all vertices is equivalent to summing 3 times over all triangles; the factor 1/2 to the fact that the error is equidistributed in both directions r 1,K , r 2,K . Also, N v denotes the number of vertices of the current mesh. The remeshing procedure is then the following : for every P, we set If (37) is not satisfied, we modify λ i,P by a factor β, else we keep it as is. Based on these stretching values, a new mesh will be generated by the BL2D software. The results for example (31) hereafter have been obtained setting β = 2 3 , while β = 1 2 was set for example (38).

Numerical Results with Adapted Meshes and Adapted Time Steps
We now analyse the efficiency of the adaptive algorithm of Fig. 1. We first consider example (31) with C = 60. The initial mesh is an isotropic mesh with mesh size h = 0.01, while the initial time step is taken as τ 1 = 0.002. The mesh and solution at final time are shown in Figs. 2 and 3 when using conservative interpolation [19], C = 240 and TOL = 0.001.   We investigate the number of vertices, aspect ratio, number of time steps and remeshings, for various values of the prescribed tolerance TOL. The notations are summarized in Table 9, the results in Table 10. The observations are the following when using conservative interpolation.
-The error at final time is approximatively divided by 2 when TOL is divided by 2.
-Both effectivity indices ei and ei Z Z are close to one, -The number of remeshing depends on the exact solution u (the larger C, the larger N m ).
-Since the solution depends only on the x 1 variable, the total number of vertices at final time is only doubled as the tolerance is divided by two (it should be multiplied by four with isotropic meshes). -The total number of time steps is multiplied by √ 2 as the tolerance is divided by 2, which confirms the second order convergence of the error indicator in time η T .
Linear interpolation, the L 2 projection, the Ritz hyperbolic projection and the modified hyperbolic projection (35) yield worse results, for instance the ZZ effectivity index is away of one. This has already been been observed in [5,17] for hyperbolic problems, whereas interpolation between meshes seems not to be an issue for parabolic problems [6].
The exact solution is not known, however, since the flow is reversed at t = 2, we must have u(x 1 , x 2 , 4) = u 0 (x 1 , x 2 ). This example is not covered by our our theory, since the velocity field β depends on time. Although the anisotropic error indicator in space η A remains valid even for a time dependent velocity field β, the error indicator in time η T should be modified. However, this is beyond the scope of the present study. Nevertheless, the use of the time indicator η T defined by (30) yields good results. Several meshes and numerical solutions are presented in Fig. 4 when    Stretching of a circle in a vortex flow. Convergence results for the adaptive algorithm with C = 60 (rows 1-4), C = 240, (rows 4-6). Conservative interpolation between meshes is used Stretching of a circle in a vortex flow. Convergence results with non-adapted meshes and constant time steps (τ 2 = O(h 3/2 )) with C = 60 (rows 1-3), C = 240, (rows 4-6) TOL = 0.025 and conservative interpolation is used. In Figs. 5, 6, 7, and 8 and Table 11 we have checked convergence of the computed solution at final time for several values of TOL. For comparison, we present in Table 12 results with non-adapted uniform meshes and constant time steps. In Fig. 9, we compare the solution computed on a non-adapted meshes with the one obtained with the largest value TOL = 0.1 of the adaptive algorithm. Clearly, the coarsest adapted solution is more accurate than the finest non-adapted one. Note that the number of vertices of the non-adapted mesh is 200 larger than that of adapted meshes.