Optimal error bounds for two-grid schemes applied to the Navier-Stokes equations

# Optimal error bounds for two-grid schemes applied to the Navier-Stokes equations

Javier de Frutos Departamento de Matemática Aplicada, Universidad de Valladolid. Spain. Research supported by Spanish MICINN under grant MTM2010-14919, and by JCyL grant VA001A10-1 (frutos@mac.uva.es)    Bosco García-Archilla Departamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research supported by Spanish MICINN under grant MTM2009-07849 (bosco@esi.us.es)    Julia Novo Departamento de Matemáticas, Universidad Autónoma de Madrid, Instituto de Ciencias Matemáticas CSIC-UAM-UC3M-UCM, Spain. Research supported by Spanish MICINN under grant MTM2010-14919 (julia.novo@uam.es)
###### Abstract

We consider two-grid mixed-finite element schemes for the spatial discretization of the incompressible Navier-Stokes equations. A standard mixed-finite element method is applied over the coarse grid to approximate the nonlinear Navier-Stokes equations while a linear evolutionary problem is solved over the fine grid. The previously computed Galerkin approximation to the velocity is used to linearize the convective term. For the analysis we take into account the lack of regularity of the solutions of the Navier-Stokes equations at the initial time in the absence of nonlocal compatibility conditions of the data. Optimal error bounds are obtained.

## 1 Introduction

In this paper we study two-grid mixed finite-element (MFE) methods for the spatial discretization of the incompressible Navier–Stokes equations

 ut−Δu+(u⋅∇)u+∇p = f, (0) div(u) = 0,

in a bounded domain () with a smooth boundary subject to homogeneous Dirichlet boundary conditions on . In (1), is the velocity field, the pressure, and  a given force field. As in , ,  we assume that the fluid density and viscosity have been normalized by an adequate change of scale in space and time. We approximate the solution and corresponding to a given initial condition

 u(⋅,0)=u0. (0)

Two-grid methods are a well established technique for nonlinear steady problems, see , . The main idea in a two-level method involves the discretization of the equations over two meshes of different size. A nonlinear system over the coarse mesh is solved in the first step of the method. In a second step, a linearized equation based on the approximation over the coarse mesh is solved on the fine mesh. In ,  several two-level methods are considered to approximate the steady Navier-Stokes equations. In these papers, depending on the algorithm, the second step is based on the solution of a discrete Stokes problem, a linear Oseen problem or one step of the Newton method over the fine mesh with the coarse mesh approximation as initial guess.

Several two-level or two-grid schemes have also been considered in the literature to approximate the evolutionary nonlinear Navier-Stokes equations (1)-(1). Again, two approximations to the velocity (and correspondingly two approximations to the pressure), are computed. One of them is defined by a discretization of the nonlinear equations over a coarse mesh, , and another one, , is defined by an appropriate linearization over a fine mesh. Different classes of algorithms can be seen as two level methods. In particular, although they were originally developed from different ideas, the so called nonlinear Galerkin methods, postprocessed and dynamical postprocessed methods, fall into this category.

Postprocessed Galerkin methods were first introduced for Fourier spectral methods in ,  (see also ) and later extended to finite element methods in , , , . In all these works the main idea is the following: one first compute the standard Galerkin approximation to the velocity and pressure over a coarse mesh of size and then compute the postprocessed approximation in a finer mesh at selected times in which one wants to obtain an improved approximation. More precisely, the postprocessed approximation computed at a given time is an approximation in a mesh of size to the following (steady) Stokes problem:

 −Δ~u+∇~p=f−ddtuH(t∗)−(uH(t∗)⋅∇)uH(t∗)\@@LTX@noalign\vskip3.0ptplus3.0ptminus1.0pt\omitdiv(~u)=0⎫⎪ ⎪⎬⎪ ⎪⎭in~{}Ω,\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus3.0ptminus1.0pt\omit~u=0,\,on~{}∂Ω.

Here, , , is the standard MFE approximation computed in the coarse mesh in a time interval and . Note that the computation of , , is completely independent of the computation of in the fine mesh. The postprocessed approximation improves the rate of convergence of the standard Galerkin approximation over the coarse mesh in the following sense. If the rate of convergence of the Galerkin approximation to the velocity in the () or () norm is then the rate of convergence of the postprocessed approximation to the velocity is +. Analogous results are obtained for the pressure. For first order mixed finite element methods the improvement in the rate of convergence of the velocity is only achieved in the norm, . Then, if one wants to achieve the optimal accuracy of the fine level in the norm, one can first compute the Galerkin approximation on a coarse mesh of size and then compute the postprocessed approximation over the fine mesh of size at the desired time levels. For example, one should take and  for linear and quadratic mixed finite elements, respectively. It can be expected that the computational cost of the postprocessed approximation is smaller than that of the Galerkin approximation on the same fine mesh, since in the former method the time evolution is done on the coarse mesh, and only at selected time levels are computations done on the fine mesh. This has been confirmed by the numerical experiments in  (see also  and )

In  a related algorithm, the so-called dynamical postprocessing, is introduced for the Fourier case. In this algorithm, the standard Galerkin approximation, , is computed over a coarse mesh in the first level, as before. For the second level an approximation to a linear evolutionary problem, instead of the steady problem (1), is computed. More precisely, the dynamical postprocessing involves the approximation, at each time step, over a mesh of size of the problem:

 ddt~u−Δ~u+∇~p=f−(uH⋅∇)uH\@@LTX@noalign\vskip3.0ptplus3.0ptminus1.0pt\omitdiv(~u)=0⎫⎪ ⎪⎬⎪ ⎪⎭in~{}Ω,\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus3.0ptminus1.0pt\omit~u=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,on~{}∂Ω.

Note that in the dynamical postprocessing, the computation of and , , is coupled. The rate of convergence of the dynamical postprocessing scheme is proved in  to be the same as the rate of convergence of the standard postprocessing. In the case of highly oscillatory solutions the dynamical algorithm is shown to be more efficient than the standard postprocessing in some one dimensional examples. The dynamical postprocessing method is also considered in , named now as two-level method, in the case of mixed finite elements. In , the author treats the fully discrete case integrating in time with the backward Euler method. A similar two-level scheme is also considered and analyzed in  where the author uses first order mixed finite elements in space, Crank-Nicolson extrapolation for the time integration over the coarse mesh and the backward Euler method for the time integration over the fine mesh.

The so-called nonlinear Galerkin methods are also two-level methods that have been used to compute approximations to (1)-(1). They were first introduced for Fourier spectral methods , , and later extended to mixed finite element methods in . In this work the authors obtain the rate of convergence of the nonlinear Galerkin method in the case of first order elements. The rate of convergence is the same one of the postprocessed method. The main difference between the nonlinear Galerkin methods and the postprocessed or two-grid methods is that in the former the approximation on the coarse mesh takes into account the influence of the fine mesh, whereas in the latter it is just the standard Galerkin method (i.e., computed independently of the fine mesh).

In this paper we analyze two two-grid algorithms in the context of spatial mixed finite element discretizations to approximate the solutions of (1)-(1). The two algorithms we consider are very similar to the dynamical postprocessing method. The difference is the treatment of the nonlinearity in the second level. In the dynamical postprocessing method the nonlinear convective term of the fine level is approximated by the data (see the right-hand-side of (1)). In the two algorithms we consider in the present paper, the approximation to the velocity of the coarse mesh is used to linearize the nonlinear convective term of the fine level. In the first algorithm, the linearized convective term of the fine level is . In the second algorithm is regarded as an initial guess to perform one Newton step in the fine level. For the spatial discretization we consider mixed finite elements of first, second and third order. More precisely, we consider the mini-element and the quadratic and cubic Hood-Taylor elements. The analysis for other mixed finite elements of the same order is similar. As in ,  due to the lack of regularity at of the solution of (1)-(1) no better than error bounds can be expected. For this reason we do not analyze higher than cubic finite element discretizations. For the temporal discretization we use the backward Euler method or the two-step backward differentiation formula. The analysis of the fully discrete methods is similar to the one appeared in  and it is only briefly indicated in this paper.

This is not the first time these two algorithms have been considered. The first algorithm was introduced in , where the authors analyze the semi-discrete in space case for first order finite elements. In  the authors extend this analysis to the fully discrete case and in  the second order Hood-Taylor finite element is used for the spatial discretization and the two-step backward differentiation formula for the time integration. In  the second algorithm is analyzed for the Fourier spectral case while in  the analysis is extended to the case of first order mixed finite elements considering the fully discrete case coupled with the Crank-Nicolson scheme for the time integration. As opposed to the above mentioned works on the same methods, in the present paper we take into account the lack of regularity suffered by the solutions of the Navier-Stokes equations at the initial time. Then, for the analysis in the present paper we do not assume more than second-order spatial derivatives bounded in up to initial time , since demanding further regularity requires the data to satisfy nonlocal compatibility conditions unlikely to be fulfilled in practical situations , . This is the first time these methods are analyzed under realistic regularity assumptions. Also, this is the first time the cubic case is considered and the first time the quadratic case is considered for the second method.

There are some other improvements with respect to previous works. In  the authors get an error bound of order for the fine approximation to the velocity in the norm whenever the following inequality is satisfied and being constants independent of and . With the technique of this paper an error bound of order for the same fully discrete method in the norm can be obtained for and independently chosen. With the new error bound obtained in this paper one can achieve the rate of convergence of the fine mesh in the norm by taking instead of . This fact improves the efficiency of the method compared with the (same order) standard Galerkin method over the fine mesh. Also, the authors of  remark that they have observed the same rate of convergence for the two-grid method with and in the numerical tests they have carried out, which supports the improved rate of convergence we obtain in this paper. We want to remark that in all the numerical experiments of ,  and  the two-grid algorithms improve the efficiency of the standard Galerkin method in the sense that a given error can be achieved with less computational cost with the new algorithms than with the standard Galerkin method. In  a comparison in the Fourier case between the standard postprocessing, the dynamical postprocessing and the second two-grid algorithm is also included. Although the computational cost of the two-grid approximation over the fine mesh is bigger than that of the postprocessed approximations, the two-grid algorithm produces smaller errors in the case of moderate to high Reynolds numbers. Finally, comparing the two algorithms we analyze in this paper we remark that with the second algorithm better error bounds are obtained in terms of . Although this fact could make the choice of the second algorithm preferable for computations, it turns out in practice to be rather inefficient to solve the linear systems accurately. For this reason, some authors suggest solving instead an Oseen problem leading then to the first algorithm, see .

The rest of the paper is as follows. In Section 2 we introduce some preliminaries and notations. In Section 3 we carry out the error analysis of the first two-grid algorithm in the semi-discrete in space case. In Section 4 we consider the analysis of the second two-grid algorithm in the semi-discrete in space case. Finally, in Section 5 we consider the fully discrete case integrating in time with the backward Euler method or the two-step backward differentiation formula.

## 2 Preliminaries and notations

We will assume that is a bounded domain in , not necessarily convex and smooth enough. When dealing with linear elements ( below) may also be a convex polygonal or polyhedral domain. We will consider the Hilbert spaces

 H ={u∈(L2(Ω))d∣div(u)=0,u⋅n|∂Ω=0}, V ={u∈(H10(Ω))d∣div(u)=0},

endowed with the inner product of and , respectively. For integer and , we consider the standard Sobolev spaces of functions with derivatives up to order in , and . We will denote by the norm in , and  will represent the norm of its dual space. We consider also the quotient spaces with norm .

Let us recall the following Sobolev’s imbedding inequalities : For , there exists a constant such that

 ∥v∥Lq′(Ω)d≤C∥v∥Ws,q(Ω)d,1q′≥1q−sd>0,q<∞,v∈Ws,q(Ω)d. (0)

For , (2) holds with .

Let be the projection onto . We denote by the Stokes operator on :

 A:D(A)⊂H⟶H,A=−ΠΔ,D(A)=H2(Ω)d∩V.

Applying Leray’s projector to (1), the equations can be written in the form

 ut+Au+B(u,u)=Πf in Ω,

where for , in .

We shall use the trilinear form defined by

 b(u,v,w)=(F(u,v),w)∀u,v,w∈H10(Ω)d,

where

 F(u,v)=(u⋅∇)v+12(∇⋅u)v∀u,v∈H10(Ω)d.

It is straightforward to verify that enjoys the skew-symmetry property

 b(u,v,w)=−b(u,w,v)∀u,v,w∈H10(Ω)d. (0)

Let us observe that for .

We shall assume that

 ∥u(t)∥1≤M1,∥u(t)∥2≤M2,0≤t≤T,

and, for integer,

 sup0≤t≤T∥∥∂⌊k/2⌋tf∥∥k−1−2⌊k/2⌋+⌊(k−2)/2⌋∑j=0sup0≤t≤T∥∥∂jtf∥∥k−2j−2<+∞,

so that, according to Theorems 2.4 and 2.5 in , there exist positive constants and such that for

 ∥u(t)∥k+∥ut(t)∥k−2+∥p(t)∥Hk−1/R≤Mkτ(t)1−k/2 (0)

and for

 ∫t0σk−3(s)(∥u(s)∥2k+∥us(s)∥2k−2+∥p(s)∥2Hk−1/R+∥ps(s)∥2Hk−3/R)ds≤K2k, (0)

where and  for some . Observe that, for , we can take  and . For simplicity, we will take these values of and . We note that no further than will be needed in the present paper.

Let , , be a family of partitions of suitable domains , where is the maximum diameter of the elements and are the mappings of the reference simplex onto . We restrict ourselves to quasi-uniform and regular meshes .

Let , we consider the finite-element spaces

 Sh,r = {χh∈C(¯¯¯¯Ωh)|χh|τhi∘ϕhi∈Pr−1(τ0)}, S0h,r = {χh∈C(¯¯¯¯Ωh)|χh|τhi∘ϕhi∈Pr−1(τ0),χh(x)=0∀x∈∂Ωh},

where denotes the space of polynomials of degree at most on . Since we are assuming that the meshes are quasi-uniform, the following inverse inequality holds for each (see, e.g., [10, Theorem 3.2.6])

 ∥vh∥Wm,q(τ)d≤Chl−m−d(1q′−1q)∥vh∥Wl,q′(τ)d, (0)

where , , and is an element in the partition .

We shall denote by the so-called Hood–Taylor element [8, 26], when , where

 Xh,r=(S0h,r)d,Qh,r−1=Sh,r−1∩L2(Ωh)/R,r≥3,

and the so-called mini-element  when , where , and . Here, is spanned by the bubble functions , , defined by , if  and 0 elsewhere, where denote the barycentric coordinates of . For these mixed elements a uniform inf-sup condition is satisfied (see ); that is, there exists a constant independent of the mesh grid size such that

 infqh∈Qh,r−1supvh∈Xh,r(qh,∇⋅vh)∥vh∥1∥qh∥L2/R≥β. (0)

The approximate velocity belongs to the discretely divergence-free space

 Vh,r=Xh,r∩{χh∈H10(Ωh)∣(qh,∇⋅χh)=0∀qh∈Qh,r−1}.

We observe that, for the Hood–Taylor element, is not a subspace of . Let be the discrete Leray’s projection defined by

 (Πhu,χh)=(u,χh)∀χh∈Vh,r.

We will use the following well known bounds

 ∥(I−Πh)u∥j≤Chl−j∥u∥l,1≤l≤r,j=0,1. (0)

Assuming that is has a smooth enough boundary, we also have

 ∥∥A−m/2Π(I−Πh)u∥∥0≤Chl+min(m,r−2)∥u∥l,1≤l≤r,m=1,2. (0)

Since , for all , it follows that

 ∥A−1/2hΠhf∥0≤C∥f∥−1. (0)

Moreover it holds for , see :

 ∥A−s/2hΠhf∥0 ≤ Chs∥f∥0+∥A−s/2Πf∥0s=1,2. (0)

Let  denote either or . Notice that both are positive self-adjoint operators with compact resolvent in  and , respectively. Let us consider then for and the operators and , which are defined by means of the spectral properties of  (see, e.g., [11, p. 33], ). An easy calculation shows that

 ∥Aαe−tA∥0≤(αe−1)αt−α,α≥0, t>0, (0)

where, here and in what follows, when applied to an operator denotes the associated operator norm. Also, using the change of variables , it is easy to show that

 ∫t0s−1/2∥∥A1/2he−(t−s)Ah∥∥0ds≤1√2eB(12,12), (0)

where is the Beta function (see, e.g., ).

## 3 Semi-discretization in space. The first two-grid algorithm.

In this section we carry out the error analysis of the first two-grid algorithm for the Hood-Taylor mixed finite element with or . At the end of the section we include the results that can be obtained for the mini-element with a similar but simpler analysis than the one showed along the section.

The first algorithm we consider is the following. Let us choose so that . Then, in the first level we compute the standard mixed finite-element approximation to (1)–(1). That is, given , we compute and , , satisfying, for all and

 (˙uH,ϕH)+(∇uH,∇ϕH)+b(uH,uH,ϕH)+(∇pH,ϕH) =(f,ϕH) (0) (∇⋅uH,ψH) =0. (0)

In the second level we solve a linearized problem on a finer grid and given we compute and , , satisfying, for all and

 (˙~uh,ϕh)+(∇~uh,∇ϕh)+(uH⋅∇~uh,ϕh)+(∇~ph,ϕh) =(f,ϕh) (0) (∇⋅~uh,ψh) =0. (0)

To obtain the error bounds for we will follow the error analysis of  and introduce an auxiliary approximation (see [14, Section 4.1]). For a and solution of (1)–(1) let us consider the approximations and , respectively, solutions of

 (˙vh,ϕh)+(∇vh,∇ϕh)+(∇gh,ϕh) =(f,ϕh)−b(u,u,ϕh) (0) (∇⋅vh,ψh) =0, (0)

for all and , with initial condition . We will also use the following notation:

 zh=Πhu−vh. (0)

Next, we state some lemmas that are needed in the proof of the main theorems. The first one summarizes previous results.

###### Lemma 1

Let be the solution of (1)–(1). There exists a positive constant such that the discrete velocity  defined by (3)-(3) and the Hood–Taylor element approximation to , , satisfy the following bounds for , and :

 ∥vH(t)−uH(t)∥j ≤ Ct(r−2)/2|log(H)|Hr+1−j,3≤r≤4, (0) ∥uH(t)−u(t)∥j ≤ Ct(r−2)/2Hr−j,2≤r≤5, (0) ∫t0∥uH(s)−u(s)∥2jds ≤ CH2(3−j),3≤r≤4. (0)
###### Proof

The bound (1) is proved in Theorems 4.7 and 4.15 in . The case in (1) is proved in [23, Theorem 3.1] and [24, Theorem 3.1]. The case follows from the case by applying (2) and (2), see also Corollaries 4.8 and 4.16 in . Finally, (1) is proved in Lemmas 5.1 and 5.2 in .

For the convenience of the reader, we will reproduce here the following two Lemmas, the first one from  and the second one from .

###### Lemma 2

For any , the following estimate holds for all :

 ∫t0∥∥Ahe−(t−s)AhΠhf(s)∥∥0ds≤C|log(h)|max0≤t≤T∥f(t)∥0.
###### Lemma 3

Let be the solution of (1)–(1). Then, there exists a positive constant such that the error in (3) satisfies the following bound:

 ∥A(−1+j)/2hzh∥0 ≤Ct(r−2)/2hr+1−j, j=0,1,2, r≥3. (0)
###### Lemma 4

For each there exist positive constants and depending on and such that, for , or , and every and satisfying the threshold condition

 ∥∥u(t)−w1h1(t)∥∥j≤αh2−j1,∥∥u(t)−w2h2(t)∥∥j≤αh2−j2,j=0,1, t∈[0,T], (0)

for , , satisfying

 ∥wh(t)∥j≤2αmax(h1,h2)2−j,

the following bounds hold:

 ∥∥\@fontswitchF(wh,w2h2)∥∥0+∥∥\@fontswitchF(w1h1,wh)∥∥0 ≤ K∥wh∥1, (0) ∥∥\@fontswitchF(wh,w2h2)∥∥−1+∥∥\@fontswitchF(w1h1,wh)∥∥−1 ≤ K∥wh∥0, (0) ∥∥\@fontswitchBh(wh,w2h2)∥∥0+∥∥\@fontswitchBh(w1h1,wh)∥∥0 ≤ K∥wh∥1, (0) ∥∥A−1/2h(\@fontswitchBh(wh,w2h2))∥∥0+∥∥A−1/2h(\@fontswitchBh(w1h1,wh))∥∥0 ≤ K∥wh∥0, (0)

where can be either or , and .

###### Proof

The proof of the present lemma can be found in that of Lemma 4.4 in  for and . With obvious changes, the proof is also valid when , as well as when .

Remark 1  We will apply the above inequalities for , and . Let us also remark that the Lemma 4 also holds if either or is replaced by . In what follows we will apply Lemma 4 to and both satisfying the threshold condition (4) for an appropriate value of  (see [14, Remark 4.1]).

###### Lemma 5

For there exists a positive constant such that the following bound holds for :

 aa∥∥A−1Π[\@fontswitchF(v,e)+\@fontswitchF(e,v)]∥∥0≤K∥v−w∥−1, (0)

where can be either or .

###### Proof

The proof of this result when can be found as part of the proof of [7, Lemma 3.4]. With obvious changes, the proof is also valid when .

Let us observe that the approximation over the finer grid and the recently defined satisfy

 ˙~uh+Ah~uh+Πh(uH⋅∇~uh)=Πhf,uh(0)=Πhu0, (0) ˙vh+Ahvh+Πh(u⋅∇u)=Πhf,vh(0)=Πhu0, (0)

respectively. Then satisfies

 ˙eh+Aheh+Πh(uH⋅∇eh)=Πhρh,H,eh(0)=0, (0)

where

 ρh,H=uH⋅∇vh−u⋅∇u.

In the proof of Theorem 1 below we will use the following lemmas.

###### Lemma 6

Let be the solution of (1)–(1). There exists a positive constant such that the following inequality holds for :

 ∥A−1hΠhρh,H∥0≤Ct(r−2)/2|log(H)|Hr+1,t∈(0,T].
###### Proof

Let us write , where

 ρ1h,H=((uH−u)⋅∇vh),ρ2h,H=(u⋅∇(vh−u)). (0)

By applying (2) we have

 ∥A−1hΠhρjh,H∥0≤Ch2∥ρjh,H∥0+∥A−1Πρjh,H∥0,j=1,2.

To bound let us recall Remark 3 and apply (4) to get

 ∥ρ1h,H∥0≤C∥uH−u∥1≤CHr−1t(r−2)/2, (0)

where we have applied (1) from Lemma 1 in the last inequality. Applying (4) we also get

 ∥ρ2h,H∥0≤C∥vh−u∥1≤Chr−1t(r−2)/2, (0)

where in the last inequality we have applied standard bounds for (see (2)) together with the estimates (3) for in Lemma 3. Let us next bound . We will use the decomposition

 ρ1h,H=((u−uH)⋅∇)(vh−u)+((u−uH)⋅∇)u. (0)

Then, we obtain

 ∥A−1Πρ1h,H∥0 = ∥A−1Π(((uH−u)⋅∇)(vh−u))∥0 (0) +∥A−1Π(((uH−u)⋅∇)u)∥0.

To bound the second term in (3) we apply (5) from Lemma 5 to get

 ∥A−1Π(((uH−u)⋅∇)u)∥0≤C∥uH−u∥−1.

Applying then (1) together with (2) and (3) we get

 ∥uH−u∥−1 ≤ ∥uH−vH∥0+∥vH−u∥−1 (0) ≤ Ct(r−2)/2|log(H)|Hr+1+Ct(r−2)/2Hr+1.

To bound the first term in (3) we argue by duality, using (2), we get