The Semi-implicit DLN Algorithm for the Navier Stokes Equations

Dahlquist, Liniger, and Nevanlinna design a family of one-leg, two-step methods (the DLN method) that is second order, A- and G-stable for arbitrary, non-uniform time steps. Recently, the implementation of the DLN method can be simplified by the refactorization process (adding time filters on backward Euler scheme). Due to these fine properties, the DLN method has strong potential for the numerical simulation of time-dependent fluid models. In the report, we propose a semi-implicit DLN algorithm for the Navier Stokes equations (avoiding non-linear solver at each time step) and prove the unconditional, long-term stability and second-order convergence with the moderate time step restriction. Moreover, the adaptive DLN algorithms by the required error or numerical dissipation criterion are presented to balance the accuracy and computational cost. Numerical tests will be given to support the main conclusions.


Introduction.
In the simulation of time-dependent fluid models, various time-stepping schemes have been constructed based on stability and consistency. The backward Euler method, unconditionally stable and easily implemented, can only have first-order accuracy [22,37,51,52]. The trapezoidal rule or two-step backward difference method (BDF2) are both second-order accurate and widely used in computational fluid dynamics [3, 11, 18, 19, 26-29, 43, 47, 54]. However, the trapezoid rule with some unfavorable combinations of time steps leads to instability of the numerical solutions [16,53]. The variable-step BDF2 method only has conditional stability if the time step ratio is small enough [8,9,20,21,30,31].
Dahlquist, Liniger, and Nevanlinna propose a one-parameter family of one-leg, two-step method (thus the DLN method herein) which is G-stable (non-linear stable) [13][14][15][16] and second-order accurate under arbitrary time grids. To our knowledge, the DLN method is the only time-stepping algorithm possessing these two properties under arbitrary time step sequence. Hence its essential properties of stability and consistency have been carefully studied and explored in [45,46]. Recently the variable step DLN method has been applied to the unsteady Stoke/Darcy model and Navier Stokes equations (NSE) and performs well in specific test problems [44,49,50].
Given the initial value problem: Here {0 = t 0 < t 1 < · · · < t N −1 < t N = T } N n=0 are the time grids on interval [0, T ] and y n is the DLN solution to y(t n ). The coefficients in (DLN) are The step variability ε n = (k n − k n−1 )/(k n + k n−1 ) is the function of two step size. k n = α 2 k n − α 0 k n−1 is the average time step. Given sequence {z n } ∞ n=0 , we denote z n,β := β α ℓ y n−1+ℓ = k n f (t n,β , y n,β ).
Let u h n and p h n are the numerical solutions to velocity u(x, t n ) and pressure p(x, t n ) respectively on certain finite element space with diameter h, the fully-implicit DLN algorithm for NSE approximate the non-linear term u ·∇u at t n+1 by u h n,β ·∇u h n,β [44], which results in two main disadvantages of the algorithm: i. the rigorous time step restriction like ∆t ≤ O(ν 3 ) for convergence, ii. the non-linear solver in each time step computation. The above time step restriction in error analysis arises from the use of the discrete Grönwall inequality [35, p.369] and would be very strict even under moderate viscosity value (like ν = 1.e − 2) 1 . Fixed point iteration and Newton's iteration are two common choices for non-linear solvers. Fixed point iteration is easily implemented while Newton's iteration possesses fast convergence. However, they usually cost more than solving a linear system and have the risk of divergence if the initial value for the iteration is poorly guessed.
To address the two issues, we extend Baker's idea [1] and propose the semiimplicit DLN scheme for NSE. The essence of the idea is to extrapolate the first u h n,β in the non-linear term u h n,β · ∇u h n,β by its second-order extrapolation in time (the inner product. The space H r with norm · r and semi-norm | · | r denotes the Sobolev space with p = 2. The velocity space X and the pressue space Q for the NSE in (1.2) are The divergence-free space for the velocity is For any function v ∈ X ∩ H r d , the norm · r and semi-norm | · | r are define . X ′ is the dual space of X with the dual norm We need the Bochner space on the time interval [0, T ] where the corresponding discrete norms are The discrete norm | · | p,r,β in (2.2) is the form of Riemann sum in which the function We apply divergence theorem and integration to ( Thus if u ∈ V , b(u, v, w) = (u · ∇v, w). We have the following lemma about the bounds of the operator b.
For spatial discretization, X h ⊂ X and Q h ⊂ Q are certain finite element spaces for velocity and pressure respectively based on the edge-to-edge triangulation of domain Ω (with the maximum diameter of the triangles h > 0). X h and Q h satisfy the discrete inf-sup condition, i.e.
where r and s are the polynomial degrees of X h and Q h respectively. The inverse inequality for X h is The discrete divergence-free space is For any pair (u, p) ∈ V × Q, the Stokes projection (P (u) The Stokes projection has the following approximations (see [25,38] for proof) 3. The algorithm. Let u h n ∈ X h and p h n ∈ Q h be the numerical solutions of velocity u(x, t n ) and pressure p(x, t n ) respectively. Then the semi-implicit DLN algorithm for the NSE in (1.2) is: given two previous solutions where the second-order, linear extrapolation u h n for u h n,β is The above semi-implicit DLN algorithm in (3.1) can be simplified by the following refactorizaion process (See [45] for the proof of equivalence) Step 1. Pre-possess: Step where the coefficents in the refactorizaion process are

Numerical Analysis.
For numerical analysis, we need the following two lemma about the stability and consistency of the DLN method.
Lemma 4.1. Let Y be the inner product space over R with the inner product (·, ·) Y and the induced norm · Y . For any sequence and the coefficents γ By the above identity in (4.1), the whole family of variable step, one-leg (DLN) methods are G-stable (see [13, p.2] for the definition).
Proof. The proof of identity in (4.1) (implicit in [16]) is an algebraic calculation.
If the mapping u(t) is smooth enough about variable t, then for any θ ∈ [0, 1]

5)
Moreover, if there exists constants C L , C U > 0 such that the ratio of time steps satisfies Proof. Using Taylor theorem and expanding u n+1 , u n and u n−1 at time t n,β .

Stability Analysis. Theorem 4.3.
If the body force f in NSE satisfies f ∈ L 2 (0, T ; X ′ )∩ℓ 2,β (0, N ; X ′ ), the semi-implicit DLN algorithm for NSE in (3.1) satisfies the following unconditional, long-time energy bounds: for any integer N > 1, where {γ Proof. We set v h = u h n,β in (3.1) and use the skew-symmetry property of the operator b We apply (2.1) and Young's inequality to (f n,β , u h n,β ) Then we use the G-stability identity in (4.1) and Lemma 4.2 By the definition of G(θ)-norm in (4.2) and the notations in (2.2), we sum (4.8) over n from 1 to N − 1 to yield (4.7).

Error Analysis.
Let r and s be the polynomial degree of X h and Q h respectively and Let u n and p n be the exact velocity and pressure of the NSE in (1.2) at time t n . We need the following upper and lower bound for the ratio of time steps: there are positive constants C L and C U such that Theorem 4.4. Suppose the velocity u ∈ X, the pressure p ∈ Q and the body force f (x, t) of the NSE in (1.2) satisfy Under the time step bounds in (4.9), the numerical solutions of the semi-implicit DLN scheme in (3.1) satisfy Proof. The exact solutions of NSE at time t n,β satisfy: where u n is second-order, linear extrapolation of u n,β in (4.4) and the truncation error τ n is and decompose the error of velocity e u n to be We restrict v h ∈ V h in (3.1) and subtract (4.11) from the first equation of (3.1) We set v h = φ h n,β in (4.13) and use the G-stability identity in (4.1), (4.14) By Cauchy-Schwarz inequality, Poincaré inequality, Young's equality, (2.8) and (2.10) By Holder's inequality We combine (4.15) and (4.16), By the definition of the Stokes projection, ν ∇η n,β , ∇φ h n,β = 0. We set The non-linear terms in (4.14) become By (2.4), (2.5), (2.9) and (2.10), Poincaré inequality and step requirement in (4.9) We apply Young's inequality to all non-linear terms in (4.18) By (2.8), triangle inequality and (4.5) in Lemma 4.2 We set q h to be L 2 projection of p(t n,β ) onto Q h in (4.14) and use (2.8) Now we treat τ n (φ h n,β ): by (4.5) and (4.6) in Lemma 4.

Since
where and F 1 is in (4.29). We apply (4.50) and discrete Gronwall inequality to (4.49) By the time-diameter condition in ( By (4.51) and (4.53) and triangle inequality which implies (4.35). Theorem 4.6. Suppose the velocity u ∈ X and pressure p ∈ Q of the NSE in (1.2) satisfy and the body force f ∈ L 2 (0, T ; X ′ ∩ L 2 ), then under the time step bounds in (4.9) and the time-diameter condition in (4.33), the pressure component by the algorithm in (3.1) satisfy  α ℓ e u n−1+ℓ , e u n = u h n − u n .
We let v h ∈ X h in (4.11) and subtract (4.11) from the first equation of (3.1)

Implementation of Adaptive DLN Algorithm.
We present two ways of time adaptivity for the whole family of DLN methods. The first way is to use the local truncation error (LTE) criterion: we apply the revised AB2 method (herein AB2-like method) to estimate the error of the DLN scheme for NSE and adjust the time step according to the ratio of the required tolerance and the estimator. The second way is to adapt the time step to control the numerical dissipation.

Local Truncation Error Criterion. Given four previous solutions
and the estimators are where We refer to [46] for the derivation of the AB2-like method in (5.1) and the estimator in (Absolute estimator) and (Relative estimator). We use the step controller proposed by Hairer and Wanner in [33] to adjust the next time step k n+1 = k n · min 1.5, max 0.2, κ Tol/ T n+1 where Tol is the required tolerance and the safety factor κ ∈ (0, 1] is selected to minimize the number of step rejections. If T n+1 > Tol, then the DLN solution is rejected and the current step k n is adjusted by (5.3) for recomputing. We summarize the above adaptive DLN algorithm in Algorithm 1  1) ; compute the AB2-like solution u h,AB2 n+1 by (5.1) ; use k n , k n−1 , k n−2 , k n−3 to update ε n , ε n−1 , ε n−2 ; compute G (n) , R (n) by (5.2) ; ; // adjust step by (5.3) else // adjust current step to recompute k n ⇐ k n · min 1.5, max 0.2, κ Tol

Numerical Dissipation Criterion.
The algorithm proposed by Capuano, Sanderse, De Angelis and Coppola [10] calibrates the step size to ensure the ratio of numerical dissipation and viscosity under the required value and its effect on the fullyimplicit DLN scheme has been tested in [44]. Given the tolerance Tol, the maximum time step k max , the minimum time step k min , we compute the DLN solution and the ratio of numerical dissipation and viscosity If χ n+1 ≤ Tol, we accept the current solutions and double the time step. Otherwise, we halve the time step for recomputing. We summarize the algorithm in Algorithm 2.

Algorithm 2: Adaptivity by Numerical Dissipation Criterion
Input: tolerance Tol, two previous solutions u h n , u h n−1 and p h n , p h n−1 , current time step k n , previous time step k n−1 , ; compute the DLN solution u h,DLN n+1 and p h,DLN n+1 by (3.1) ; compute numerical dissipation E ND n+1 and viscosity E VD n+1 ; // double the step else k n ⇐ max 0.5k n , k min ; // halve current step to recompute 6. Numerical Tests. We apply the semi-implicit DLN algorithm in (3.1) with θ = 2/3, 2/ √ 5, 1 for all numerical tests. θ = 2/3 is suggested in [16] to balance the stability and local truncation error. θ = 2/ √ 5 is recommended in [41] for stability at infinity (a property similar to L-stability). θ = 1 reduces to the midpoint rule. We use software FreeFem++ and Taylor-Hood (P 2 − P 1) finite element space for programming.
We set the parameters ω = 1, τ = 1/ν = 100. The initial value, boundary condition and source function f are determined by the exact solutions in (6.1). We require that the constant time step k and mesh diameter h are the same to satisfy the timediameter condition in (4.33). We simulate the problem over the time interval [0, 1]. The convergence rate R is calculated by The results of the semi-implicit DLN algorithm (with constant time step) are given in Tables 6.1 to 6.6. We see that the semi-implicit DLN has third-order convergence in velocity and second-order convergence in pressure for all three θ values. Hence the semi-implicit DLN scheme has much better performance in the Taylor-Green problem than the theories in Subsection 4.2 suggest. Then we apply the fully-implicit DLN scheme to the same problem and use fixed point iteration to solve the non-linear system at each time step. The error and convergence rate are given in Tables 6.7 to 6.12. From the above tables, we observe that the semi-implicit DLN algorithm outperforms the fully-implicit algorithm for all three θ values in this test problem: The two schemes have almost the same error magnitude but the fully-implicit scheme takes twice the time to finish the simulation since the fully-implicit scheme takes two iterations on average at each time step.
We set ω = 1 and τ = 1/ν = 2500. The exact solutions in (6.2) make the problem more difficult since the Reynolds number is much larger and the energy has an increasing pattern. We use both Algorithm 1 and Algorithm 2 to solve the problem over the time interval [0, 60]. For Algorithm 1, we use the relative estimator in (Relative estimator) and set tolerance Tol = 1.e − 7 and the safety factor κ = 0.95. For Algorithm 2, we set Tol = 1.e − 14 for χ. The value of Tol is chosen to balance accuracy and efficiency. For both algorithms, we set the minimum time step k min = 0.0005, the maximum time step k max = 0.05, the initial time step k 0 = 0.0005 and the mesh diameter h = 1/180. The initial value, boundary value and body force are decided by the exact solutions. We measure the performance of two algorithms by evaluating energy, error of energy, numerical dissipation E ND n+1 and viscosity E VD n+1 . Since E ND n+1 vanishes for the DLN method with θ = 1, we test the two adaptive algorithms with θ = 2/3 and θ = 2/ √ 5. Fig. 6.1 shows the performance of two adaptive DLN algorithms and Table 6.13 tells us the number of steps.
Algorithm 1 surpasses Algorithm 2 in terms of accuracy and efficiency: Fig. 6.1(b) shows that the error magnitude of energy is much smaller for Algorithm 1 while Algorithm 2 takes more number of time steps for both θ values. Two algorithms have similar patterns for numerical dissipation and viscosity in Figs. 6.1(c) and 6.1(d). For both algorithms, T n+1 and χ n+1 are kept below the required tolerance after the first few steps in Fig. 6.1(e) and time steps oscillate between k max and k min in Fig. 6.1(f).

2D Offset Circles Problem.
We use the 2D offset circles problem proposed by Jiang and Layton [37] to verify the stability of the DLN scheme under any arbitrary sequence of time steps and the efficiency of the adaptive algorithms in Section 5. The domain Ω ⊂ R 2 is Ω = {(x, y) : x 2 + y 2 ≤ 1 and (x − 0.5) 2 + y 2 ≥ 0.01}. For both algorithms, T n+1 and χ n+1 are kept below the required tolerance after the first few steps and time steps oscillate between k max and k min . The flow in the domain Ω is driven by the rotational body force with the no-slip boundary condition on both circles. We set the Reynolds number Re = 1/ν = 200 and simulate the problem over time interval [0, 60]. We use the relative estimator of LTE in (Relative estimator) for Algorithm 1 and set the tolerance Tol = 0.001 and safety factor κ = 0.95. For Algorithm 2, we set Tol = 0.01 for ratio χ n+1 . For both adaptive algorithms, the initial time step k 0 = 0.005, the maximum time step k max = 0.05 and the minimum time step k min = 0.0005. The domain triangulation is generated by 80 nodes on the boundary of the inner circle and 320 nodes on the boundary of the outer circle. Since the exact solutions are unknown, we use the constant step DLN algorithm (ε n = 0) with a small time step (k = k min ) and refined mesh (100 nodes on the boundary of inner circle and 400 nodes on the boundary of outer circle) for reference. Figs. 6.2(a) and 6.2(b) show two domain triangulations. The number of time steps is presented in Table 6.14. Fig. 6.3(a) shows that the energy of all algorithms is increasing at the start and then come to the steady level 23 at time t = 8. We deduce that Algorithm 2(θ = 2/ √ 5) has worse performance for this problem because the energy level of this algorithm is low compared to that of other adaptive algorithms. In addition, the number of time steps is least for Algorithm 2 (θ = 2/ √ 5) while the ratio χ goes above the required tolerance value (1.e − 2) for many times in the simulation. Then we compare Figs. 6.3(b) to 6.3(d) and observe that the energy of Algorithm 1(θ = 2/ √ 5) is closer to the energy of reference algorithms with less number of time steps. From Fig. 6.4(a), we can see that the numerical dissipation of Algorithm 1 is at a level as low as that of reference algorithms while that of Algorithm 2 is much larger. All the algorithms have similar viscosity patterns in Fig. 6.4(b). From Figs. 6.4(c) and 6.4(d), T is always below the required tolerance (1.e − 3) thus the time steps of Algorithm 1 never reach k min . However the ratio χ goes above the required tolerance 1.e − 2 frequently and k n = k min occurs very often. The primitive time step controller (doubling and halving time steps) in Algorithm 1 reduces the number of time steps and may cause inaccuracy.      The numerical dissipation of Algorithm 1 is at a level as low as that of reference algorithms while that of Algorithm 2 is much larger. All the algorithms have similar viscosity patterns. T is always below the required tolerance (1.e − 3) thus the time steps of Algorithm 2 never reach k min . The ratio χ goes above the required tolerance 1.e − 2 frequently and k n = k min occurs very often.

Conclusions.
We propose the semi-implicit DLN scheme for the NSE and avoid non-linear solvers at each time step. G-stability of the DLN method results in the long-term, unconditional stability of the numerical solutions. In the error analysis, we prove that both the velocity and pressure of the variable time-stepping, semi-implicit scheme converge in second order under very moderate time conditions. Two adaptive algorithms based on local truncation error and numerical dissipation criteria are presented to improve time efficiency in practice. The advantage of the semi-implicit DLN scheme is observed in numerical tests in Subsection 6.1: the semiimplicit scheme obtains the same accuracy as the fully-implicit scheme and reduces the simulation time by half. Subsection 6.2 shows that two adaptive DLN algorithms obtain enough accuracy in energy and negligible numerical dissipation even the problem with a large Reynolds number has an increasing energy pattern. We verify in the 2D offset problem that the semi-implicit DLN scheme is unconditional, long-time stable in energy with any arbitrary sequence of time steps, and the adaptive DLN algorithm is much more efficient than constant time-stepping DLN scheme (taking less number of time steps and attaining similar magnitude in energy, numerical dissipation and viscosity).

Acknowledgement.
The author thanks Professor Catalin Trenchea (Department of Mathematics, University of Pittsburgh) for very helpful suggestions and discussions.