An adaptive algorithm for the transport equation with time dependent velocity

An a posteriori error estimate is derived for the approximation of the transport equation with a time dependent transport velocity. Continuous, piecewise linear, anisotropic finite elements are used for space discretization, the Crank-Nicolson scheme scheme is proposed for time discretization. This paper is a generalization of Dubuis S, Picasso M (J Sci Comput 75(1):350–375, 2018) where the transport velocity was not depending on time. The a posteriori error estimate (upper bound) is shown to be sharp for anisotropic meshes, the involved constant being independent of the mesh aspect ratio. A quadratic reconstruction of the numerical solution is introduced in order to obtain an estimate that is order two in time. Error indicators corresponding to space and time are proposed, their accuracy is checked with non-adapted meshes and constant time steps. Then, an adaptive algorithm is introduced, allowing to adapt the meshes and time steps. Numerical experiments are presented when the exact solution has strong variations in space and time, illustrating the efficiency of the method. They indicate that the effectivity index is close to one and does not depend on the solution, mesh size, aspect ratio, and time step.


Introduction
Space-time adaptive algorithms are efficient tools to approximate solutions of partial differential equations with accuracy and low computational cost. Whenever possible, the adaptive criteria is based on theoretical error estimates, this is mostly the case for elliptic and parabolic problems, fewer results are available for hyperbolic problems [11,14,18,31,37], nonlinear systems [3,13,30,39,41] or PDEs with variable coefficients [9,24].
The classical theory of a posteriori error analysis for finite element methods was first developed on isotropic meshes [6,16,40], the involved constants were depending on the mesh aspect ratio. However, anisotropic finite elements, that is to say elements with possibly large aspect ratio, have been widely used to approximate phenomena involving boundary or internal layers. The isotropic theory for a posteriori error estimates was therefore updated, see for instance [4,19,22,23,25], and the involved constants were proved to be aspect ratio independent whenever the mesh was aligned with the solution.
The Crank-Nicolson method is a popular second order scheme for time dependent problems. However, most of the a posteriori error estimates are proved for first order methods only, for instance the Backward Euler scheme [8,10,32,38]. Moreover, standard a posteriori proofs lead to suboptimal estimates for second order methods, as reported in [2]. To circumvent this problem, piecewise quadratic reconstructions have been introduced for parabolic problems [1,2,7,28], in [21] for the transport equations and in [27] for the wave equation.
In this paper, the a posteriori error analysis of [21] for the transport equation is extended to the case where the velocity varies in space and time. As in [21], the a posteriori error analysis is valid for anisotropic meshes and a second order time discretization. In [21], the a posteriori error estimate was restricted to the case of a steady velocity while the new estimate is valid for velocity fields that may exhibit strong variations in time. The techniques used should be useful when coupling the transport equation to Navier-Stokes equations, which is the case for instance when considering free surface flows [20].
The outline of this paper is the following. In Sect. 2, the transport equation and the numerical method are presented. An a posteriori error estimate is proposed in Sect. 3, numerical experiments confirm the sharpness of the estimate in Sect. 4. Finally, an adaptive algorithm is presented in Sect. 5 and numerical results are discussed.

Statement of the problem and numerical scheme
Given a polygon Ω , a divergence free velocity field ∈ C 1 Ω × [0, T ] , and an initial data 0 ∈ C 0 Ω , we are looking for ∶ Ω × (0, T ] → ℝ the solution of the transport problem where the inflow boundary is defined by with standing for the unit outer normal of Ω . It is assumed that the velocity field is such that the inflow boundary does not depend on time and that the data , Ω, 0 are smooth enough to justify the forthcoming computations. To approximate in space (1) a finite element stabilization scheme is needed. This scheme corresponds to the stabilized scheme studied in [12], in which the stabilization term is updated to the anisotropic setting, following [11,21,29].
The fully discrete method reads as follows. For every h > 0 , let T h be a conformal triangulation of Ω into triangles K of diameter h K ≤ h . Let V h be the set of continuous, piecewise linear functions on each triangle K ∈ T h , zero valued on Γ − . Let N be a non-negative integer and a partition 0 = t 0 < t 1 < t 2 < ⋯ < t N = T . We denote n+1 = t n+1 − t n the time step, t n+ 1 ∕2 = (t n+1 + t n )∕2 , Here h > 0 stands for an O(h) stabilization parameter that will be specified in the sequel.
A priori error estimates for (2) have been proposed in [12] for isotropic meshes, constant time steps and a transport velocity independent of time. In [21], the analysis is extended to anisotropic meshes and variable time steps, and both a priori and a posteriori error estimates are derived. An a posteriori error analysis involving a time dependent transport velocity is proposed here, so as an adaptive algorithm. The corresponding a priori error estimates can be found in [20], so as additional numerical experiments.
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 [22,23,29], see also [25] for similar results. Let K ∈ T h and T K ∶K ⟶ K be the affine transformation mapping the reference triangle , into K defined by with M K ∈ ℝ 2×2 , K ∈ ℝ 2 . Observe that since M K is invertible it admits a singular value decomposition M k = R T K Λ K P K , 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 correspond to the value of maximum and minimum stretching.
A posteriori error analysis can be obtained using Clément's interpolant [17]. Since anisotropic meshes are considered, the usual regularity assumption is omitted but two additional assumptions are needed. First, we assume that each vertex has a number of neighbours bounded from above, uniformly with respect to h. Second, we assume that for each K, the diameter of T −1 K (ΔK ) (here ΔK is the union of triangles sharing a vertex with K) is uniformly bounded with respect to h. For more details, we refer again to [22,23,29]. In this framework, the following estimation holds where R h is the Clément's interpolant, Ĉ > 0 is a constant depending only on the reference triangle K , and with G K standing for the gradient matrix given by

A posteriori error estimates
To recover an error estimate of order two in time, we follow the idea introduced in [28] for the Crank-Nicolson approximation of the heat equation, and define a piecewise quadratic reconstruction of the numerical solution n h of (2). We introduce the following notations. For any quantity w n and for n = 0, … , N − 1 we note and for n = 1, … , N − 1 Setting = max( 1 , … , N ) , we define the piecewise numerical reconstruction h by for ( , t) ∈ Ω × t n , t n+1 , n ≥ 1, and by for ( , t) ∈ Ω × t 0 , t 1 .
Note that (5) is the unique quadratic polynomial interpolating exactly n+1 h , n h , n−1 h at time t = t n+1 , t n , t n−1 . This quadratic reconstruction was first introduced in [28] to recover an a posteriori error estimate for the Crank-Nicolson method applied to the heat equation which was of order two in time. It is modification of another quadratic reconstruction proposed in [2], where it was indeed observed that linear reconstructions would yield to suboptimal error estimates. Here, for the transport equation (1), the motivation is the following. A priori error estimates for a simplified differential problem indicate that the error is order two, the constant depending on ttt [20], which is formally equal to − tt ( ⋅ ∇ ) . Therefore, some information about the second derivative in time of must be available in the error indicator, which is precisely the role of the quadratic term in (5).
Our error indicator for the time discretization is obtained by inserting the numerical reconstruction (5) (6) into the transport equation.
be the solution of (2). Let h be the numerical reconstruction (5) (6). For all 0 ≤ n ≤ N − 1 , for all v h ∈ V h and for all t ∈ (t n , t n+1 ) , we have where n is given for n ≥ 1 by and for n = 0 by and n is given for n ≥ 1 by for some t in (0, T) and for some t in (0, T). Therefore, n (for n ≥ 1 ) is an approximation of and is thus of optimal order O( 2 ) . We can also approximate n (for n ≥ 1) by thus, since h = O(h) , n is of higher order O(h 2 ). (ii) Note that if and h are independent of time, then n = 0 , n = 0, 1, … , N − 1, and n reduces to the time error indicator given in [21], Lemma 1.
Proof Let n ≥ 1 and t ∈ (t n , t n+1 ) , we have: Observe that I 3 is already part of (7) and order two in time, so it remains to transform I 1 and I 2 into quantity of second order. A straightforward computation for I 1 yields that Using the numerical scheme (2), the first term is in fact zero, remaining the following expression for I 1 : To simplify I 2 , we note that it looks like the discrete derivative of (2). Thus we first compute the difference between (2) taken at two successive steps, that is to say and we divide by ( n+1 + n )∕2 . We obtain that We now use the fact that we can transform the centered finite difference operator ̄ into the backward finite operator through the relation Plugging this last relation into (11), we get By multiplying the last equality by (t − t n+ 1 ∕2 ) it yields for I 2 the following expression Finally, in order to have finite differences of the velocity, we add and substract the terms and We therefore obtain as final expression for I 2 Now, by adding the first term of I 2 with I 3 and the third and the last terms with I 1 yields the result for n ≥ 1.
Finally, for n = 0 , reproducing the steps using to compute I 1 here above yields that ◻ In order to prove an a posteriori error estimate for the numerical method (2) we now need to define the stabilization parameter h in (2). Following [11,29], h is defined for all t ∈ [0, T ] and for K ∈ T h by if (t) is not identically zero on K and by h (t) |K = 0 otherwise.
where is the solution of (1). Let ( n h ) N n=0 be the solution of (2) with h defined by (12) and consider h the numerical reconstruction given by (5) (6). Setting e = − h , there exists Ĉ > 0 depending only on the reference triangle, in particular independent of T , Ω, , , the mesh size, aspect ratio and time step such that where K is defined in (4) , n , n are defined in Proposition 1 and c 0 = 1 , c n = T − 1 , for n ≥ 1.

Remark 2
(i) The a posteriori error estimate (13) is not standard since the exact solution is contained in K (e) . To obtain a computable upper bound, post-processing techniques were advocated in [11], for instance Zienkiewicz−Zhu (ZZ) post-processing [42][43][44].
More precisely, to compute G K ( − h ) , we replace the first order partial derivatives with respect to x i where, for any v h ∈ V h and for any vertex P of the mesh It corresponds to the weighted mean value of the gradient ∇v h over the support of the hat functions centered at vertex P. Superconvergence of the ZZ recovery has been proved for elliptic problems and structured meshes [44] and more recently for unstructured anisotropic meshes [15]. Numerical experiments have shown that the efficiency of ZZ post-processing is better than theoretical predictions, see for instance [11,28,[33][34][35][36].
(ii) Observe that if is independent of time, then due to definitions of n and n , the a posteriori error estimate (13) reduces to the one proven in [21], Theorem 2. (iii) Based on the a priori error estimates in [12,20], it is expected that ‖e(T )‖ 2 L 2 (Ω) = O(h 3 + 4 ) (written in the isotropic setting for simplicity). Numerical experiments will confirm that and The term is a mixed quantity involving both space and time discretizations. It can be shown (see point (iv) below) to be of higher order when 2 = Θ(h 3 ∕2 ) , that is exactly the goal of the adaptive algorithm proposed fur ther. Since the initial error , the leading terms of the estimate (13) reduces to that can be used as an a posteriori estimator for the error ‖e(T )‖ 2 .
(iv) Assuming that the solution is smooth enough [20], it is expected that Since the goal of our adaptive algorithm is to equidistribute the error due to space and time, then 2 ≃ h 3 ∕2 and the above term is O(h 3.5 + 5+1∕3 ) , that is to say of higher order compared the error at final time ‖e(T )‖ 2 L 2 (Ω) = O(h 3 + 4 ).

Proof
In the following, we denote by Ĉ any positive constant, which may depend only on the reference triangle and may change from line to line. In particular, Ĉ is independent of T, Ω, , , the mesh size, aspect ratio and the time step. Let t ∈ t n , t n+1 , n ≥ 1 . First observe that for all v ∈ H 1 (Ω) that is zero on Γ − , we have Therefore, applying (14) with v = e and using the fact that is the solution to (1) we have � Subtracting any v h + h (t) (t) ⋅ ∇v h and using Proposition 1, we obtain Splitting the last inequality into a sum over the triangles, and using the triangle and Cauchy-Schwarz inequalities yields We now choose v h = R h e where we recall that R h stands for the Clément's interpolant. Using the anisotropic interpolation error estimate (3), we can prove that Indeed, using that ∇e = (∇e ⋅ 1,K ) 1,K + (∇e ⋅ 2,K ) 2,K ) and implies that Thus, since we obtain (15) by applying the interpolation error estimate (3). Moreover using (15) and the definition of h|K (t) , we can finally prove that So we obtain that where we have set . The discrete Cauchy-

Schwarz implies then that
Finally, the Young's inequality yields We conclude by using a Gronwall's type inequality. Multiplying by exp(−t∕(T − t 1 )) on both sides and integrating between t 1 and T yields Note that the choice to use exp(−t∕(T − t 1 )) is made in order to eliminate the exponential growth with respect to T in the estimate. Proceeding in the same manner we can obtain an estimate for ‖e(t 1 )‖ 2 Combining both estimates yields finally

Numerical experiments with non-adapted meshes and constant time steps
We denote e(T ) L 2 the L 2 numerical error at final time where we recall that c 0 = 1 , c n = T − 1 , n ≥ 1 , and n is defined as in Proposition 1. As already anticipated in the Remark 2, the space error indicator (17) is made computable by replacing all the derivatives of by their ZZ postprocessing. The sharpness of the error indicator will be investigated by computing the effectivity index ei given by Since we use the ZZ post-processing to make the space error indicator computable, we also check the efficiency of this procedure by computing the true H 1 error and the approximated H 1 error and the ZZ effectivity index is denoted by To test the convergence of the numerical method presented above and check the sharpness of our error indicators, we start with numerical experiments for nonadapted meshes and constant time steps. The numerical experiments should possibly demonstrate the following properties: (i) To validate the equivalence between the error indicator and the true error, the effectivity index ei should remain close to a constant, in particular, the effectivity index ei should not depend on the solution, the time step and the mesh size and aspect ratio. (ii) To ensure that the ZZ post-processing is asymptotically exact and justify its use, the effectivity index ei ZZ should be close to one.
We consider a "1D" problem where the initial condition is given by with C > 0 and we solve the transport equation with Finally, we impose Dirichlet boundary conditions on the left side of Ω that is chosen as (0, 1) 2 and we set T = 0.05 . The exact solution is then given by Observe that the solution is smooth, with small variations, except in a thin layer of width controlled by C. 10t 2 , 0).
The numerical results are reported in Table 1 where we choose = O(h 2 ) . When the error is mainly due to the space discretization, it is observed that ei stays close to a value of 20, as it was already concluded for non-depending on time velocity [21]. ei ZZ is close to one, thus the ZZ postprocessing is asymptotically exact, as predicted by the theory. The numerical errors are respectively ≃ O(h 1.8 ) for the L 2 error at final time and ≃ O(h) for the L 2 (0, T ;H 1 (Ω)) error. When we choose h = O( 2 ) , the numerical error is mainly due to the time discretizazion and the effectivity index is close to 2. This value is already observed in [21] and we note that both numerical errors are O( 2 ) , see Table 2. In all cases, ei is independent of the solution (the same effectivity indices are obtained with C = 60 and C = 240 ), the mesh size and aspect ratio, and the time step as expected.
In order to have effectivity indices that are close to one, we finally perform numerical experiments with the normalized error indicator defined by Numerical experiments will show that indeed, the effectivity indices corresponding to the normalized error indicator will be close to one and independent of the solution, the mesh size and aspect ratio, and the time step. This normalized indicator will be then used as a criterion towards mesh and time step adaptivity. In Table 3, the effectivity indices are reported when using (21)

An adaptive algorithm
We now present an adaptive algorithm. The a priori error analysis contained in [12,20] show that the final error increases as the square root of the final time T. Therefore, our goal is to ensure that the relative error over the time e(T ) L 2 T 1 ∕2 = TOL , where TOL is a preset tolerance chosen by the user. Since the normalized error indicator (21) is shown to be close to the true error, the goal of the adaptive algorithm is to build a sequence of meshes and time steps such that T 1 ∕2 stays close to TOL. Therefore, we would like that at the end of the simulation, verifies A sufficient condition so that (22) holds is to equidistribute the error between the space and time approximations. Therefore, we build a sequence of meshes such that and a sequence of time steps such that Finally, sufficient conditions at every steps n = 0, 1, .., N − 1 to satisfy the above criteria are and (22) Table 1 Convergence results when = O(h 2 ) with C = 60 (rows 1-6) and C = 240 (rows 7-12) The aspect ratio is 50 (rows 1-4 and 7-9) and 500 (rows 5-6 and 10-12). The space error is dominated and ei is computed to be close to 20 If the conditions (26) are not satisfied, we refine or coarsen the time step. If (25) are not satisfied, the mesh is changed and a new anisotropic mesh is generated with the BL2D software [26]. In practice, the new mesh is built by equidistributed A K ,n between the direction 1,K and 2,K , and by aligning each triangle K with the eigenvectors of the gradient matrix G K ( − h ) (that is post-processed with ZZ post-processing). For more details on this procedure, we refer to [20]. Everytime a new mesh has to be built, the old values of the solutions have to be interpolated on the meshes. A detailed discussion and several interpolation operator are presented in [11] and [21]. In both studies, the conservative interpolation [5] demonstrates the best result. Therefore we decide to use it also for the numerical experiments presented below. The main steps of the adaptive algorithm are summarized in Table 4. Table 2 Convergence results when h = O( 2 ) with C = 60 (rows 1-7) and C = 240 (rows [8][9][10][11][12][13][14][15] Thhe aspect ratio is 50 (rows 1-4 and 8-11) and 500 (rows 5-7 and 12-15). The temporal error is dominated and ei is computed to be close to 2  Table 3 Convergence results with the normalized error indicator when h 3 ∕2 = O( 2 ) with C = 60 (rows 1-7) and C = 240 (rows [8][9][10][11][12][13][14] The aspect ratio is 50 (rows 1-4 and 8-11) and 500 (rows 5-7 and 12-14). Here ei = √ A 20 2 + T 2 2 ∕e(T) L 2 and is computed to be around 1 Example 1 (An anisotropic example) We apply the adaptive algorithm to the following example. The domain is still chosen as (0, 1) 2 and we set the final time to T = 0.3 . Dirichlet boundary conditions are still imposed on the left side of Ω and we choose the initial condition as before as with C = 60 or C = 240. We now choose (x 1 , x 2 , t) = (ūt(t), 0) with ūt given by  (25) and (26) are satisfied then accept the current mesh and time step (25) is not satisfied, then adapt the mesh For K ∈ T n h,i , compute G K and align the mesh with the eigenvectors of G K For i = 1, 2 Compute the mesh size i,K If the mesh size is too small in direction x i increase i,K If the mesh size is too big in direction x i decrease i,K Build a new mesh T n h,i+1 using the BL2D software If n > 0 Interpolate the old solutions on the new mesh is too big, decrease the current time step Maximum aspect ratio at final time, the aspect ratio on an element K being 1,K ∕ 2,K ar : Average aspect ratio at final time ūt is a smooth function that is mainly 1 or 10, except in a small boundary layer of width 0.03 around t =t . Since the velocity quickly accelerates after t =t , it is expected that smaller time steps are chosen by the algorithm. For this particular example,we choose t = 0.25. We run the adaptive algorithm for various values of the prescribed tolerance TOL. We investigate the number of vertices, aspect ratio, number of time steps and remeshings. We summarize the notations we used for the analysis in Table 5. We make the following observations: -The number of vertices is multiplied by 1.6 as the tolerance is divided by two, expressing the fact that the method is O(h 3 ∕2 ).
In Tables 6, we present the converge results for several values of TOL. In Fig. 1, we check the evolution of the time step. We observe that the time step follows the evolution of , up to some oscillations occurring when the adaptive algorithm refuses the current time step. In particular, when is 10 times larger, the adaptive algorithm selects a time step that is approximatively 10 times smaller. In Fig. 2, we represent the generated meshes and the solutions for C = 240 and TOL = 0.00125. This example demonstrates the efficiency of anisotropic adaptive finite elements to approximate the behaviour of solutions with boundary layers. When TOL = 0.00125 , the mesh size in the x 1 direction is 0.0013 for C = 60 and 0.0003 for C = 240 . An isotropic adaptive algorithm would need the same mesh size in both x 1 and x 2 directions, thus resulting in million vertices! Morover, numerical experiments show that the effectivity index corresponding to the normalized error indicator (21) remains close to one (say between 0.9 and 1.6).
Example 2 (An isotropic example) We present an isotropic example in order to show that the effectivity index indeed does not depend on the solution. The initial condition is given by where ūt(t) is as in (27) and we choose t = 0.125 . The final time is set to 0.15. Numerical results are reported in Table 7 when runing the adaptive algorithm. In Fig. 3, we represent the generated meshes and the solutions at time t = 0, 0.125, 0.15. The same observations can be made, namely that the effectivity index is close to one (it ranges between 1 and 1.7) and the ZZ post-processing is asymptotically exact.  Since the flow is reversed at t = 2 , we must have (x 1 , x 2 , 4) = 0 (x 1 , x 2 ). We start the adaptive algorithm with an initial grid of mesh size h = 0.1 and an initial time step 1 = 0.001 . Several meshes and numerical solutions are presented in Figs. 4 and 5 when TOL = 0.00125 . In Figs. 6 and 7, convergence of the computed solution at final time is checked for several values of TOL.

Conclusion
We prove an a posteriori error estimate for the space-time approximation of the transport equation in the case where the transport velocity depends on space and time. Anisotropic finite elements are considered and the Crank-Nicolson method is used to advance in time. The corresponding a posteriori error estimate is shown to be of optimal order. Error indicators for space and time are proposed, numerical experiments confirm their sharpness. An adaptive algorithm is then introduced in order to capture with accuracy and low computational cost solutions that have strong variations in space and time. The efficiency of the method is shown for three 2D examples. A few 3D computations are reported in [20].