Stabilized bi-grid projection methods in Finite Elements for the 2D incompressible Navier-Stokes

We introduce a family of bi-grid schemes in finite elements for solving 2D incompressible Navier-Stokes equations in velocity and pressure $(u,p)$. The new schemes are based on projection methods and use two pairs of FEM spaces, a sparse and a fine one. The main computational effort is done on the coarsest velocity space with an implicit and unconditionally time scheme while its correction on the finer velocity space is realized with a simple stabilized semi-implicit scheme whose the lack of stability is compensated by a high mode stabilization procedure; the pressure is updated using the free divergence property. The new schemes are tested on the lid driven cavity up to $Re=7500$. An enhanced stability is observed as respect to classical semi-implicit methods and an important gain of CPU time is obtained as compared to implicit projection schemes.


Introduction
Multigrid methods have been widely developed since more than 40 years and were proposed as fast solvers to the numerical solution of elliptic problems, see e.g. [26] for steady linear and nonlinear Dirichlet Problems in finite differences, but also for Steady Navier-Stokes Equations [8,19].
The method is first defined for two levels of discretization (bi-grid case). The two key ingredients are both the separation of the high and the low mode components of the error provided by the use of coarse and fine grids (or spaces) V H and V h respectively, and also the concentration of the main computational effort on the coarset (lower dimensional) subspace V H ⊂ V h ; high modes can be represented in V h while only low modes can be captured in V H . This leads to a drastic save in CPU computation time. The correction to the fine space components belonging to V h is usually realized by a simple and fast numerical scheme and is associated to a high mode smoothing. Then the scheme is recursively applied on a set of nested grids (or spaces).
When considering nonlinear dissipative equations, we can distinguish two approaches for bi-grid methods: In the first one, the low and the high mode components are explicitly handled: thanks to the parabolic regularization, it is expected that the low mode components which carry the main part of the energy of the (regular) solution have a different dynamics from the high mode components, that can be seen as a correction, see [15] and the references therein. Also, a way to speed up the numerical integration is then to apply different schemes to these two sets of components, concentrating the effort on the computation of the low modes, that belong in V H , see [13]. When dealing with spectral methods, the separation in frequency is natural but it is not the case when Finite Differences or FEM methods are used for the spatial discretization. To separate the modes, a hierarchical basis approach is used [43]: using a proper interpolation (or projection) operator between V H and V h , one builds a transfer operator which defines a pre-conditioner for the stiffness matrices but it allows also to express the solution in terms of main part, associated to the low mode components, and of a fluctuant part, of lower magnitude, and associated to high modes components; we refer the reader to [40,41,42] and the references therein for Finite Elements discretizations, [12,34] for Finite Differences and [6,17] in Finite Volumes. These schemes showed to be efficient, however they necessitate to build and manipulate hierarchical bases.
In the second approach, the mode separation is not used explicitly and the methods consist, at each time step, in first computing the coarse approximation u H to the fine solution u h by an unconditionally implicit stable scheme and then to update the fine space approximation by using a linearized scheme at an extrapolated valueũ h of u H in V h . These schemes allow to reduce the computational time with an optimal error as respected to the classical scheme when choosing accurately the mesh size of V h as respect to the one of V H . It is to be underlined thatũ h represents the mean part of the solution while z h = u h −ũ h ∈ V h represents the fluctuant part (carrying the high mode components of the solution u h ∈ V h ) but is not simulated in the schemes. Such methods were developed and applied to the solution of time dependent incompressible Navier-Stokes equations [2], [3], [20] and [27]. Since a linearization is used on V h the matrix to solve at each time iteration changes and it can be costly in transient regime.
The numerical methods we propose here are somehow in between and are inspired from the approach developed in [1] for Allen-Cahn Equation: as previously described, the use of two levels of discretization allows to concentrate the effort on the coarset, yet lower dimensional, space and at the same time to decompose the solution into its mean and fluctuant part u h =ũ h + z h . The z h are not explicitly simulated, however, they are used implicitly for a high mode stabilization as follows: consider the time integrations of the reaction-diffusion equation, in its variational form: where V is a proper Hilbert space.
Given two finite dimensional subspaces of V , W H and V h , with dim(W H ) << dim(V h ), we define the bigrid-scheme as Algorithm 1 Bi-grid Stabilized scheme for Reaction-Diffusion 3: for k = 0, 1, · · · do 4: Step 1 (Coarse Space Implcit Scheme) 5: Step 2 (Fine Space semi-implicit Scheme) 7: It is not necessary to have W H ⊂ V h , an inf-sup like compatibility condition has to be satisfied for defining uniquely the prolongationũ, [1]. The competition between the terms τ (u k+1 h − u k h , φ h ) and τ (u k+1 H − u k H , φ h ) are interpreted as a high mode filter and τ > 0 is the stabilization parameter [1] and allow to compense the explicit treatment of the nonlinear term. Furthermore, this stabilization has few influence on the global dynamics, [1].
The aim of the present work is to adapt this two-grid scheme to Navier-Stokes equations. To this end, we consider a projection method that splits the time resolution into two steps: firstly a parabolic equation on the velocity and secondly an elliptic equation on the pressure. The bi-grid stabilization will be applied to the first step.
The paper is organized as follows: in Section 2 we first present briefly the principle of the bi-grid framework (reference scheme, separation of the modes, high mode stabilization), then, after recalling the definition and some properties of the projection scheme (reference scheme), we propose different bi-grid projection schemes. Section 3 is dedicated to the numerical results. We consider the classical benchmark lid-driven cavity problem. When computing the steady states, the results we obtain agree with those of the literature, also an important gain of CPU time is obtained as respect to classical methods, for a comparable precision. The article ends in Section 4 with concluding remarks. All the computation have been realized using FreeFem++ [18].

Derivation of the Bi-grid Projection Schemes
We here build the bi-grid high mode stabilized projection schemes. For that purpose, we first recall the different approaches of the bi-grid schemes in finite elements and then describe the stabilization procedure, for reaction diffusion problems. Then, we present briefly the projection schemes in finite elements we will start from to introduce the new Bi-grid Projection Schemes.

Principle of the bi-grid approach
As stated in the introduction, when considering nonlinear dissipative equations, a well known way to obtain a gain in CPU time is to consider two levels of discretization and to concentrate the main computational effort on the lower dimensional approximation space W H , while the higher accurate approximation to the solution on the fine space V h is updated by using a simple semi-implicit (yet linear) scheme. The general pattern of a bi-grid scheme for the reaction diffusion equation (1)  3: for k = 0, 1, · · · do 4: Step 1 (Coarse Space Implicit Scheme) 5:

6:
Step 2 (Fine Space semi-implicit Scheme) 7: This approach was proposed by [20,27,29]. However, f (u k+1 H ), the linearized of f at u k+1 H must be computed at each time step, changing the matrix to solve at each iteration. A way to avoid this drawback is to consider a classical semi-implicit scheme and to compense the lack of stability by adding, e.g., a first order damping term τ (u k+1 h − u k h ), obtaining as second step: 1: Step 2 (Stabilized Fine Space semi-implicit Scheme) This stabilization procedure allows to obtain unconditionally stable time scheme for large values of τ > 0, that can be tuned. The resulting scheme is fast, however, it can slow down drastically the dynamics, particularly convergence to steady states occurs in longer times. This is due to the fact that the damping acts on all the mode components, including the low ones which are associated to the mean part of the function and carry the main part of the L 2 energy, when considering a Fourier-like interpretation. A way to overcome this drawback is to apply the stabilization to the only high mode components z h of u h which correspond to a fluctuent part of u h ; the stability of a scheme is indeed related to its capability to contain the propagation of the high frequencies. Hence the scheme becomes 1: Step 2 (High Mode Stabilized Fine Space semi-implicit Scheme) At this point we can distinguish two main ways to decompose u h as a sum of its mean part u h and its fluctuent part z h in finite elements: is the coarse FEM space that can capture only low modes, and W h is the complementary space, generated by the basis functions of V h that do not belong in V H , and that capture high mode components. When using P 1 finite elements, the components of W h of a function u h ∈ V h are built as proper local interpolation error from u H ∈ V H . The functions of approximation V h can be uniquely decomposed as It has been showed that the linear change of variable S : u h → (u H , z h ) provides an efficient preconditioner for the stiffness matrices but also proceeds to a scale separation [9,10,30,31,43].
The spatial discretization reads then to a differential system satisfied by u H and z h , namely New time marching methods are obtained applying two different schemes for u h and for z h ; particularly, the separation of the scales allows to use a simple and fast scheme for the z h components and to save computational time, with a good accuracy. Approximations and simplifications can be applied to each equation, particularly using an approached law linking high and low mode components of type Φ(u H ) = z h , [30,31].
The approach followed in [32] is technically different but is based on the same principle: Denoting by P H the L 2 orthogonal projection from V h on V H , the decomposition u h = u H + w h with u H = P H u h and w h = (Id − P H )u h is used, making u H and w h orthogonal. They proposed and analyzed schemes of Nonlinear Galerkin type, says in which an asymptotic approach law Φ(u H ) = w h ) is implemented together with a time marching scheme on u H . A drawback is that one must build a basis for (Id − P H )V h . Finally let us mention similar developments in finite differences [11,12] and spectral methods [13,15,16,28] • A L 2 -like filtering, used in [1] and to which we will concentrate: we define first the prolongation operator P : Then, settingũ h = P(u H ), we write We now apply the high mode stabilization to the z h components and get, after usual simplifications 1: Step 2 (High Mode Stabilized Fine Space semi-implicit Scheme) we can write the high mode stabilized bi-grid scheme as 3: for k = 0, 1, · · · do 4: Step 1 (Coarse Space Implicit Scheme) 5: Step 2 (Fine Space semi-implicit Scheme) 7: We now consider incompressible Navier-Stokes equations. FEM Bi-grid methods for incompressible NSE as proposed, e.g., in [2,3,20,27,29], are built using the pattern of scheme 2, and written as following. Given two pairs of FEM spaces (X H , Y H ) and (X h , Y h ) satisfying the discrete inf-sup condition, with X H ⊂ X h and Y H ⊂ Y h , one defines the bi-grid iterations as Algorithm 4 Bi-grid Scheme for Navier-Stokes Equation 3: for k = 0, 1, · · · do 4: Step 1 (Coarse Space Implicit Scheme) 5: Step 2 (Fine Space semi-implicit Scheme) 8: We have used here the simple linearization of the nonlinear term (u k+1 H .∇u k+1 h , φ h ) proposed in [20], but other linearisations can be considered. At each step the method needs to solve first a mixed nonlinear FEM problem on the coarse pair FEM spaces (X H , Y H ), then a linear mixed FEM problem on the fine FEM spaces (X h , Y h ); in this last step, the underlying matrix changes at each iteration. To avoid this drawback and to apply the stabilized bi-grid method, as described in Algorithm 3, we will use a projection method. In such a case the computation of the velocity is decoupled from the one of the pressure, the intermediary velocity satisfies a nonlinear convection-diffusion equation. The main idea is then to apply a stabilized bi-grid scheme to this equation to speed up the resolution. We present below several options of this approach.

Bi-grid Projection Schemes in Finite Elements
First of all, let us recall the framework of the Projection Schemes in Finite elements and then derive the new bi-grid stabilized methods.
The mixed variational formulation of the motion of a viscous and incompressible fluid in a domain Ω is described by the unsteady Navier-Stokes equation where Ω is a domain of R d , d = 2, 3 of lipschizian bound ∂Ω, − → n the normal outside unit and a time interval [0, T ], T > 0.
Here is a choice of the approximation spaces X h and Y h for the velocity and the pressure The degrees of freedom for the velocity are the vertices of the triangulation and the midpoints of the edges of the triangles of the triangulation T h . The degrees of freedom for the pressure are the vertices of T h assumed to be uniformly regular. It is well known that because of the constraint of incompressibility, the choice of X h and Y h is not arbitrary. They must satisfy a suitable compatibility condition, the "inf-sup" condition of Babuska-Brezzi (cf. [4,7]): where β * is independent of h. All the results presented in this paper are done with the Taylor-Hood finite element P 2 /P 1 .
Let us consider the semi-discretization in time and focus on marching schemes. Let u k u(x, k∆t) be a sequence of functions; ∆t is the time step. We compare our method to the projection method applied on the fine grid. It is a fractional step-by-step method of decoupling the computation of the velocity from that of the pressure, by first solving a convection-diffusion problem such that the resulting velocity is not necessarily zero divergence; then in a second step, the latter is projected onto a space of functions with zero divergence in order to satisfy the incompressibility condition (cf. [14,23,24,25,35,37,38]). We start by presenting the implicit reference scheme [36]: This will be our reference scheme, used on the coarse FEM Space. It enjoys of unconditional stability properties, see [36], chapter III, section 7.3.

Remark 2.1
We could use also an incremental projection method introduced by Goda [21], see also [23]. Goda has proven that adding a previous value of the gradient of the pressure (∇p k ) in the first step of the projection method then rectify the value of the velocity in the second step will improve the accuracy, in other words, he proposed the following algorithm: The results obtained are similar to the one produced by the non incremental scheme to which we focus for a sake of simplicity.
We can now derive the bi-grid schemes. Our first approach consists in separating in scales the intermediate velocity u * h which we introduce according to the projection method of Chorin-Temam [14,35,37,38]. First we compute u * H on the coarse space X H and then stabilize the high frequencies of u * h on the fine space X h . Based on the free divergence condition, we find the pressure p h whose mean is equal to zero. We end up finding u h and restricting it to the coarse grid to get u k H for the next iteration.
Algorithm 7 Two grids Algo1 1: for k = 0, 1, · · · do 2: We propose in the second algorithm to replace (u k h .∇)u k h by (u * H · ∇)u * h . We obtain: Algorithm 8 Variant of the Two grids Algo1 : Two grids Algo2 1: for k = 0, 1, · · · do 2: The second approach that we consider in the following is to find directly u k+1 H by the projection method applied on the coarse grid. By this technique, we do not need to restrict u k+1 h to X H to define u k+1 H .

The Lid Driven Cavity
We simulate the Navier-Stokes equations in 2D by a stabilization technique. We use the unit square mesh [0, 1] × [0, 1] and vary the dimensions of the coarse and fine FEM spaces. We choose as spaces of approximation the well-known Taylor-Hood element P 2 /P 1 for the velocity and the pressure respectively. Here are the results of the driven cavity flow obtained by adopting the variational formulation provided by the projection technique (Chorin-Temam) proposed by the two preceeding steps. In this case the velocity is imposed only in the upper boundary with u = (1, 0) and zero Dirichlet conditions are imposed on the rest of the boundary (see Fig. 1), below  These flow patterns in Figures 2, 3, 4, 5, 6, 7 fit with earlier results of Bruneau et al. [8], Ghia et al. [19], Goyon [22], Pascal [33] and Vanka [39]. Also, the numerical values and the

CPU Time reduction
In Figures 8 -10 we observe that the bi-grid schemes are faster in computation time than reference the implicit scheme (Algorithm 5) applied on the fine space V h . However, the gain in CPU time is mainly obtained in the transient phase. Indeed, since a stationary solution is computed, the reference scheme will need only one nonlinear iteration at each the time step in the neighborhood of the steady state: it means that only a linear system has then to been solved at each iteration, exactly as for the semi-implicit scheme. Therefore, in a neighborhood of the steady state, an iteration of the bi-grid method needs an additional implicit iteration on the coarse grid and becomes then more expensive in CPU time.

CPU time
One grid Two grids Algo1 Two grids Algo2 Two grids Algo3 Two grids Algo4

CPU time
One grid Two grids Algo1 We also compute the number of nonlinear iterations as a function of time and compute the L 2 norm of ∂u ∂t ; we use the approximation ∂u ∂t | t=t k u k+1 − u k ∆t .We identify the time from which the reference scheme (Algorithm 5) reduces to only one linear iteration at each time step and we define θ s as the associated threshold value of ∂u ∂t L 2 (Ω . We then can obtain an heuristic criteria to define a simple strategy to save more CPU time and we modify the bi-grid scheme as follows: as long as ∂u ∂t > θ s we apply the bi-grid scheme and as soon as ∂u ∂t L 2 (Ω ≤ θ s we

CPU time
One grid Two grids Algo1 Figure 10: From left to right: CPU time of the scheme on a fine grid (Algorithm 5) and that of the scales separation method (Algorithm 7) for Re = 5000, 7500 respectively.
apply the reference scheme on the fine space V h . We represent in Tables 1, 2 Table 4: Non incremental Navier-Stokes for Re = 2000, 3200, 5000, 7500 and P 2 /P 1

Enhanced Stability : comparison with a semi-implicit scheme
As presented in the previous sections, the bi-grid methods are faster than the unconditionally stable reference scheme (Algorithm 5). To illustrate the stabilization properties of the bi-grid schemes, we make now a comparison with the semi-implicit scheme (applied on whole the fine space V h ) and which consists in treating the nonlinear term (u k · ∇)u k explicitly. Particularly we compare the maximum time steps ∆t that can be taken for each of the schemes and compare the respective final CPU time on the time interval [0, T f ]. We give hereafter in Tables 5, 6 and 7 the results comparing the stability of the bi-grid Algorithms 7, 8, 9, 10 and the following semi-implicit scheme: Algorithm 11 Non-incremental semi-implicit scheme 1: for k = 0, 1, · · · do 2:    Table 6: We observe in Figure 11 that the dynamics of the convergence to the steady state of the bi-grid method is similar to the one of the implicit and semi-implicit schemes. We note also that the stability of the bi-grid scheme is guaranteed for large values of τ . For example, for Re = 1000 and ∆t = 0.007, Algorithm 7 was unstable for τ = 0.5 but taking a τ = 30 we can have a stable scheme without deteriorating the history of the convergence to the steady state. In order to locate our gain at time step ∆t, we present in Table 8 a comparison between the maximum time step allowing the stability of the semi-implicit scheme and our Algorithm 7 for the minimum value of τ necessary for stabilization.   For small Reynolds number (Re = 100), we do not obtain an enhanced stability (larger ∆t) as compared to the semi-implicit scheme. The latter is stable even for large values of the time step ∆t. When considering larger values of Re, says Re = 400 andRe = 1000, the bi-grid scheme appears as 10 times to 50 times more stable than the semi-implicit one. As a consequence, the gain in CPU time to compute the steady state increases: the stability of the semi-implicit scheme is limited by ∆t = 0.005 while the bi-grid algorithm remains stable for τ = 100 and ∆t = 0.1 and 0.5.

Precision test: Comparison with an exact solution
In order to show the accuracy of the bi-grid schemes, we simulate the Navier-Stokes equations for a Bercovier-Engelman type solution [5], say and with the corresponding term f . We can observe in Figures 12 and 13 that the solutions coincide with the reference one at final time T = 56 and that the L 2 norm of the error bi-grid schemes is of the order of ∆t on all interval time.

Concluding Remarks and perspectives
The new bi-grid projection schemes introduced in this work allow to reduce the CPU time when compared to fully implicit projection scheme, for a comparable precision. Their stability is enhanced as compared to semi-implicit projection scheme (larger time steps can be used). The stabilization we used on the high mode components is simple and does not deteriorate the consistency, the numerical results we obtained on the benchmark driven cavity agree with the ones of the literature, particularly the dynamics of the convergence to the steady state is close to the one of classical methods. This is a first and encouraging step before considering a strategy of multi-grid or multilevel adaptations of our methods, using more than 2 nested FEM spaces. We have used here FEM methods for the spatial discretization, but the approach is applicable to others discretizations techniques such as finite differences or spectral methods.