Error analysis and iterative solvers for Navier–Stokes projection methods with standard and sparse grad-div stabilization

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19 www.elsevier.com/locate/cma

Error analysis and iterative solvers for Navier–Stokes projection methods with standard and sparse grad-div stabilization Abigail L. Bowers a,1, Sabine Le Borne b, Leo G. Rebholz a,⇑,1 b

a O-110 Martin Hall, Department of Mathematical Sciences, Clemson, SC 29634 USA Hamburg University of Technology, Department of Mathematics, Schwarzenbergstrasse 95E, 21073 Hamburg

Received 9 December 2013; received in revised form 6 February 2014; accepted 28 February 2014 Available online 12 March 2014

Abstract This paper shows that use of a recently introduced sparse grad-div stabilization can increase the accuracy of projection methods for solving the Navier–Stokes equations. Sparse grad-div stabilization has recently been introduced as an alternative to standard grad-div stabilization which has a sparser matrix representation. For both sparse and standard grad-div stabilized projection methods, we prove error estimates and provide numerical experiments which reveal that both stabilizations can cause a significant decrease in the error. We then compare iterative solvers for the linear systems of equations arising from the use of both of the stabilizations. A theoretical analysis of a simplified model problem as well as numerical tests show that iterative solvers perform better for systems arising from sparse grad-div compared to standard grad-div stabilized systems. Ó 2014 Elsevier B.V. All rights reserved.

Keywords: Navier–Stokes equations; Projection methods; Sparse grad-div stabilization; Mass conservation; Iterative solvers

1. Introduction We consider an efficient stabilization of grad-div type within the framework of projection method timestepping for approximating solutions of the incompressible Navier–Stokes equations (NSE). On a domain X  Rd (d = 2 or 3), the non-dimensional form of the NSE reads as ut þ u  ru þ rp  Re1 Du ¼ f ;

ð1:1Þ

r  u ¼ 0;

ð1:2Þ

⇑ Corresponding author. Tel.: +1 8646561840. 1

E-mail addresses: [email protected] (A.L. Bowers), [email protected] (S. Le Borne), [email protected] (L.G. Rebholz). Partially supported by NSF grant DMS1112593.

http://dx.doi.org/10.1016/j.cma.2014.02.021 0045-7825/Ó 2014 Elsevier B.V. All rights reserved.

2

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

where u represents the velocity, p the pressure, f an external forcing, and Re the Reynolds number. The Reynolds number is defined to be Re :¼

UL ; m

where L and U are the characteristic length and velocity for a particular flow problem, and m is the kinematic viscosity. It is well known that small Reynolds numbers produce smooth, laminar flows, while large Reynolds number flows are associated with turbulence. For flow simulation, the NSE must be equipped with appropriate initial and boundary conditions. For simplicity of analysis, we will consider wall-bounded flows, so that the velocity satisfies u ¼ 0 on the boundary @X; however, our numerical experiments will use non-homogeneous Dirichlet boundary conditions. The difficulties associated with approximating solutions to the NSE are widely known. The presence of small physical scales leads to the need for fine meshes and thus to discrete systems with (at least) millions of degrees of freedom for most modern problems of interest. How to solve the resulting linear systems at each timestep is a very difficult problem, and even getting an answer to a fluid flow problem can be a major hurdle. Perhaps the most successful temporal discretization methods of the NSE, in terms of their ability to yield more easily solvable linear systems, are projection methods. These methods were originally developed by Chorin and Temam [6,35], and lead to more attractive linear systems because they decouple the pressure and solenoidal constraint from the momentum conservation equation. For simplicity, we will restrict our study to first order methods, although our results and proof techniques have direct analogs for BDF2 methods [17]. The linearized first order projection method (of ‘pressure-correction’ type) for the NSE is given by: Given an initial guess u0 ¼ e u 0 , a timestep 0 < Dt 6 T , for n=0; 1; . . . ; M  1 ð¼ T =Dt  1Þ, Step 1: Find unþ1 such that unþ1  e un þ un  runþ1  Re1 Dunþ1 ¼ f ðtnþ1 Þ; Dt unþ1 j@X ¼ 0:

ð1:3Þ ð1:4Þ

Step 2: Find ðe u nþ1 ; pnþ1 Þ such that e u nþ1  unþ1 þ rpnþ1 ¼ 0; Dt re u nþ1 ¼ 0; e u

nþ1

 nj@X ¼ 0:

ð1:5Þ ð1:6Þ ð1:7Þ

This method has been well studied [18,16,32–34,17], and it can easily be seen that it yields two linear systems that are easier to solve than a coupled system, although the problem to be solved in Step 1 can still be challenging for large Re. It is discussed in a number of works, e.g., [17], how the problem to be solved in Step 2 can be reformulated as a pressure Poisson problem, and thus Step 2 is usually not considered problematic. There are drawbacks to projections methods, however. Compared to coupled methods, projection methods tend to be less accurate due to the error that arises from the splitting of the NSE system. Also, there are now two velocity solutions, and it is not clear what either of them mean from a physical point of view. The velocity solution unþ1 of Step 1 satisfies the boundary condition, but mass conservation is not imposed on it. The velocity solution e u nþ1 of Step 2 is incompressible, but only satisfies the boundary condition in the normal component. Furthermore, from (1.5), we observe that rp  nj@X ¼ 0, which is not a physical condition and can lead to numerical boundary layers in the pressure [17]. Recently, it has been shown that stabilizations of grad-div type increase accuracy in coupled (i.e., without operator splitting) finite element simulations of Stokes and NSE, by reducing the divergence error of the velocity, and furthermore to reduce the scaling of the velocity error with the pressure discretization error [20,30,31,10,22,8,5]. Although its intent was somewhat different, It has appeared at least as far back as the development of the streamline upwind Petrov Galerkin formulation of [3], with the interpretation of a least

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

3

squares penalty term. While the authors are aware of the use of grad-div stabilization by CFD practitioners, much less can be found in the literature about grad-div stabilization applied to projection methods. One interesting recent result is from [25], where it is proven that in some particular finite element settings, if the term crðr  unþ1 Þ is added to Step 1, then the Step 1 velocity solution converges to the solution of a typical coupled scheme as c ! 1, which is known to have optimal convergence rate in both space and time [21]. In [33], the term brðr  ut Þ was added to Step 1, along with some other modifications to the scheme, that increased accuracy and enforced mass conservation in the Step 1 solution. A potential drawback to both of these methods is that in a typical finite element discretization, the term ðr  uh ; r  vh Þ2 will appear in the Step 1 momentum equation with a large coefficient. Since the matrix that arises from this term is singular, large parameters typically slow down the convergence of iterative solvers in the solution of Step 1. However, in both of these works, accuracy improvements were still evident even for smaller coefficients, such as O(1), which are small enough so that the extra difficulty that arises in solving the Step 1 linear system will not become prohibitive. Interestingly, despite its name, grad-div stabilization is typically not thought of as a stabilization operator as much as it is considered a consistent term added to finite element discretizations that decreases divergence error and reduces the effect of the pressure error on the velocity error. Untying stability from accuracy can be difficult, but we note there are several examples in the literature of grad-div stabilization acting to remove significant oscillations that appear in time dependent Navier–Stokes channel flow simulations over a step [22] and around a cylinder [7], and in steady Rayleigh-Bernard convection [11]. The purpose of this paper is to study the projection method together with a new sparse grad-div stabilization. This new sparse grad-div stabilization was proposed by the authors in [26] as a more efficient alternative to standard grad-div stabilization; details of this new stabilization are given in section two. For completeness and for comparison, we will also consider projection methods together with standard grad-div stabilization. We will prove that the error in projection methods can be decreased with the use of either one of these stabilizations, and provide some benchmark numerical tests that demonstrate the improvement. We will also discuss and compare iterative solvers for the sparse and standard grad-div stabilization applied to Step 1, which demonstrate the efficiency advantage offered by the sparse stabilization. This paper is arranged as follows. In Section 2, we introduce notation and both the standard and sparse grad-div stabilization. We summarize some properties of the sparse grad-div stabilization that will be important in the subsequent analysis. In Section 3, we present an a priori error analysis for projection methods including no stabilization, the standard and the sparse grad-div stabilization. The theoretical results are supported by numerical tests for the errors obtained with the three methods. Section 4 addresses the iterative solution of the linear system of equations to be solved in Step 1. We provide a convergence analysis for a simplified model problem and show numerical tests using various common preconditioners for the systems of the two stabilized approaches. 2. Notation and preliminaries Throughout this work, we will consider the domain X  Rd , d = 2 or 3, to be connected and either have smooth boundary or be a convex polygon. The L2 ðXÞ norm and inner product will be denoted by k  k and ð; Þ, respectively. All other norms will be clearly labeled. The natural function spaces for the methods we consider are defined by X :¼ H 10 ðXÞ ¼ fv 2 H 1 ðXÞ; vj@X ¼ 0g; Y :¼ fv 2 L2 ðXÞ; v  nj@X ¼ 0; r  v ¼ 0g;   Z q¼0 : Q :¼ L20 ðXÞ ¼ q 2 L2 ðXÞ; X

The spaces X and Y are velocity spaces that will contain solutions for Steps 1 and 2, respectively, and the space Q will contain pressure solutions.

2

We denote the L2 ðXÞ norm and inner product by k  k and ð; Þ.

4

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

We will study timestepping methods with an endtime T, timestep Dt ¼ T =M where M denotes the number of timesteps. Approximations to a function / at time tn :¼ nDt are denoted by /n . For continuous in time variables, we use the notation k/k1;k :¼ max k/ðtn ÞkH k ðXÞ 16n6M

and for discrete in time variables, k/k1;k :¼ max k/n kH k ðXÞ : 16n6M

2.1. Grad-div stabilization Grad-div stabilization has existed as a divergence stabilization for at least two decades [10,8]. It arises from the addition of 0 ¼ crðr  uÞ to the momentum equation of the incompressible Stokes or NSE (or related equations), and then appears as a nonzero penalty term in discretizations that do not strongly enforce mass conservation. Our intent herein is to consider this and a related stabilization applied to Step 1 of the projection method. 2.2. Sparse grad-div stabilization We introduce now the sparse grad-div stabilization operator recently proposed in [26], and review some of its important properties. Definition 2.1. The sparse grad-div stabilization (i.e., divergence penalization) operator g is defined by ! ðu2 Þxy ; ð2:1Þ g2d ðuÞ ¼ rðr  uÞ þ ðu1 Þxy 0 1 ðu2 Þxy  ðu3 Þxz B C ðu1 Þxy g3d ðuÞ ¼ rðr  uÞ þ @ ð2:2Þ A: ðu1 Þxz Lemma 2.1. The operator g satisfies the following, for u; v 2 H 10 ðXÞd : 1. The inner product of gðuÞ with v can be written as ðg2d ðuÞ; vÞ ¼ ðr  u; r  vÞ  ððu1x ; v2y Þ  ðu2y ; v1x ÞÞ;

ð2:3Þ

ðg3d ðuÞ; vÞ ¼ ðr  u; r  vÞ  ððu1x ; v2y Þ  ðu2y ; v1x Þ þ ðu1x ; v3z Þ  ðu3z ; v1x ÞÞ;

ð2:4Þ

2. g satisfies in 2d and 3d 2

ðgðuÞ; uÞ ¼ kr  uk :

ð2:5Þ

3. If r  u ¼ 0, then in 2d and 3d ðgðuÞ; vÞ ¼ ðu1x ; r  vÞ:

ð2:6Þ

Proof. These properties are proven in [26]. h Remark 2.1. The sparsity advantage in a finite element discretization of the g operator compared to standard grad-div stabilization becomes clear from (2.3) and (2.4). The variational form of the standard grad-div stabilization has interactions between each ui and vj , and thus produces a block-full matrix. However, stabilization using the g operation in 2d has no interaction between u1 and v2 , and in 3d has no interaction between u1 and v2 and u1 and v3 . Thus in 2d, the g operator produces a block upper triangular matrix, and in 3d it produces a matrix with empty (2,1) and (3,1) blocks as illustrated in Fig. 1 (right).

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

5

Fig. 1. Shown above are typical sparsity plots for the matrices arising in Step 1 of the projection method for standard grad-div (left) and sparse grad-div (right) stabilization.

The results of Lemma 2.1 reveal another interpretation for the g operator: it results from the variational formulation of the sum of grad-div stabilization and ru1x , since for v 2 H 10 ðXÞ and divergence-free u, it holds that rðr  uÞ; vÞ þ ðru1x ; vÞ ¼ ðr  u; r  vÞ  ðu1x ; r  vÞ ¼ ðgðuÞ; vÞ: In NSE-type systems, gradient functions are often combined with the true pressure to create a modified pressure (P ¼ p þ cu1x in this case). When such a change of variables is made, the velocity solution remains unchanged. The most common example of this is the use of Bernoulli pressure in formulations of the rotation form of the NSE [15]. Thus, g has the interpretation of grad-div stabilization combined with a pressure alteration that creates a more efficient linear algebraic system compared to the standard grad-div stabilization. Furthermore, we note that in 2d, u2y , and in 3d, either u2y or u3z , could be used instead of u1x . If it is expected that one of these will be “smaller” for a particular application problem, then this would cause a reduced modification of the pressure. The results presented herein would be identical, and the proofs would be analogous. The resulting matrices would of course be different, with the zero blocks occurring in different locations: matrices would be lower triangular in 2d, and two different off-diagonal blocks would be zero in 3d. Thus, the efficiency gain would be the same. 3. Error reduction in projection methods with stabilizations of grad-div type In this section, we will prove basic error estimates for the velocity in standard and sparse grad-div stabilized projection methods. These estimates show that the use of either one of these stabilizations can reduce the total H 1 error, since the divergence part of the total error will have a larger coefficient on the left hand side if c > m. We then give numerical results for two experiments which show the improvement in accuracy that these two stabilizations provide. u 0 2 Y \ H 2 ðXÞ, and that the forcing We will assume that the initial condition satisfies u0 ¼ e 1 2 f 2 L ð0; T ; L ðXÞÞ. Solutions to the systems below are understood to be in the distributional (weak) sense. The three projection methods we consider differ in the computation of Step 1 in terms of their grad-div stabilization. They are given as follows: Unstabilized Step 1: Find unþ1 2 X satisfying unþ1  e un þ un  runþ1  Re1 Dunþ1 ¼ f ðtnþ1 Þ: Dt

ð3:1Þ

6

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

Standard grad-div stabilized Step 1: Find unþ1 2 X satisfying unþ1  e un þ un  runþ1  Re1 Dunþ1  crðr  unþ1 Þ ¼ f ðtnþ1 Þ: Dt

ð3:2Þ

Sparse grad-div stabilized Step 1: Find unþ1 2 X satisfying unþ1  e un þ un  runþ1  Re1 Dunþ1 þ cgðunþ1 Þ ¼ f ðtnþ1 Þ: Dt

ð3:3Þ

Step 2: Find ðe u nþ1 ; pnþ1 Þ 2 ðY ; QÞ satisfying e u nþ1  unþ1 þ rpnþ1 ¼ 0; Dt re u nþ1 ¼ 0; e u

nþ1

ð3:4Þ ð3:5Þ

 nj@X ¼ 0:

ð3:6Þ

We note that the linearized algorithms are used only for simplicity of analysis. The analogous nonlinear schemes would produce very similar results, with the main difference being a timestep restriction on the convergence results of Theorem 3.1 which would result from the use of a different Gronwall inequality. 3.1. A priori error analysis of grad-div stabilized projection methods We now analyze the temporal error in the projection method when it is used with sparse grad-div stabilization. For comparison, we will also present results for cases of no stabilization and standard grad-div stabilization. Theorem 3.1. Let ðv; qÞ solve the NSE on ð0; T   X for a given initial condition v0 2 X \ V and forcing f 2 L1 ð0; T ; L2 ðXÞÞ, with smoothness v 2 L1 ð0; T ; X \ L1 ðXÞÞ;

vt 2 L1 ð0; T ; L2 ðXÞÞ;

vtt 2 L1 ð0; T ; L2 ðXÞÞ;

q 2 L1 ð0; T ; H 1 ðXÞÞ:

Then the error in the projection methods satisfy, 8Dt > 0, Unstabilized Step 1 + Step 2: M   X kvðT Þ  uM k2 þ Re1 Dt krðvðtn Þ  un Þk2 6 C  ðv; Re; T Þ ReDt2 þ Dtkqk21;1 :

ð3:7Þ

n¼1

Standard Grad-div Stabilized Step 1 + Step 2: M M X X 2 2 2 kvðT Þ  uM k þ cDt kr  ðvðtn Þ  un Þk þ Re1 Dt krðvðtn Þ  un Þk n¼1

n¼1

  2 6 C ðv; Re; T Þ ReDt2 þ Dtkqk1;1 : 

ð3:8Þ

Sparse Grad-div Stabilized Step 1 + Step 2: with the additional assumption that v 2 L1 ð0; T ; H 2 ðXÞÞ, 2

kvðT Þ  uM k þ cDt

M X

2

kr  ðvðtn Þ  un Þk þ Re1 Dt

n¼1

  2 6 C ðv; Re; T Þ ReDt2 þ Dtkq  cv1x k1;1 : 

M X

krðvðtn Þ  un Þk

2

n¼1

ð3:9Þ

Remark 3.1. From Lemma 3.8 of [12] it is known that if X is a convex polygon or convex with a C 2 boundary, then for any / 2 H 10 ðXÞ, kr/k2 ¼ kr  /k2 þ kr  /k2 : Because Step 1 does not enforce mass conservation, the divergence part of the error of the Step 1 solution can potentially be a significant part of the total error. Thus the improvement in the error estimates of the stabilized

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

7

P n n 2 methods comes from the term cDt M n¼1 kr  ðvðt Þ  u Þk appearing in the left hand side of the error equation: if c ¼ Oð1Þ and m is small, which is the typical case, then this estimate shows that the divergence part of the energy error can be significantly reduced, which in turn can significantly reduce the total error. This phenomena is clearly observed in a numerical experiment later in this section, where we observe that in the unstabilized case, the divergence error makes up more than half of the total H 1 error of the Step 1 solution, and that if stabilization is added (either sparse or standard), the total H 1 error is reduced by almost exactly the size of the divergence error. Remark 3.2. We observe from Theorem 3.1 that sparse grad-div stabilization can have a negative effect on the error if q  cv1x  q, which is not surprising since this stabilization term’s derivation includes adding the term cðrv1x Þ to the NSE. Hence the size cv1x relative to q is important for the choice of c in sparse grad-div stabilization. Proof. For the unstabilized Step 1 + Step 2 estimate, subtract (3.1) from the NSE at time tnþ1 , and use Taylor’s theorem to get 1 nþ1 ðe  e e n Þ þ en  rvðtnþ1 Þ þ un  renþ1  Re1 Denþ1 ¼ rqðtnþ1 Þ þ Dtðvtt ðtn Þ þ vt ðtn Þ  rvðtnþ1 ÞÞ; Dt ð3:10Þ where en :¼ vðtn Þ  un ; e e n :¼ vðtn Þ  e u n ; tn 6 tn 6 tnþ1 , and tn 6 tn 6 tnþ1 . Multiplying by enþ1 and integrating over the domain gives 1 2 2 2 2 ðkenþ1 k  ke e n k þ kenþ1  e e n k Þ þ Re1 krenþ1 k 2Dt ¼ ðen  rvðtnþ1 Þ; enþ1 Þ þ Dtððvtt ðtn Þ þ vt ðtn ÞÞ  rvðtnþ1 Þ; enþ1 Þ  ðrqðtnþ1 Þ; enþ1 Þ: Since r  e e n ¼ 0, we can use Cauchy–Schwarz and Young’s inequalities to write ðrqðtnþ1 Þ; enþ1 Þ ¼ ðrqðtnþ1 Þ; enþ1  e enÞ 6

1 2 2 kenþ1  ee n k þ Dtkrqðtnþ1 Þk : 4Dt

For the first nonlinear term, Holder’s inequality, Sobolev embeddings, Poincare’s inequality, and Young’s inequality provide the bound 2

2

ðen  rvðtnþ1 Þ; enþ1 Þ 6 Cken kkrvðtnþ1 Þk1 kenþ1 k 6 CRekrvðtnþ1 Þk1 ken k þ

Re1 2 krenþ1 k : 4

The second nonlinear term is majorized using Holder’s and Young’s inequalities, yielding Dtððvtt ðtn Þ þ vt ðtn ÞÞ  rvðtnþ1 Þ; enþ1 Þ ¼ Dtððvtt ðtn Þ þ vt ðtn ÞÞ  renþ1 ; vðtnþ1 ÞÞ 6 Dtkvtt ðtn Þ þ vt ðtn Þkkrenþ1 kkvðtnþ1 ÞkL1 ðXÞ   Re1 2 2 2 6 CDt2 Re kvtt ðtn Þk þ kvt ðtn Þk kvðtnþ1 ÞkL1 ðXÞ þ krenþ1 k 4   2 2 2 6 CDt2 Re kvtt kL1 ðtn ;tnþ1 ;L2 ðXÞÞ þ kvt kL1 ðtn ;tnþ1 ;L2 ðXÞÞ kvðtnþ1 ÞkL1 ðXÞ þ

Re1 krenþ1 k: 4

Combining estimates and multiplying both sides by 2Dt gives   1 2 2 2 2 kenþ1 k  ke e n k Þ þ kenþ1  e e n k þ Re1 Dtkrenþ1 k 2   6 2Dt2 krqðtnþ1 Þk2 þ CDt3 Re kvtt k2L1 ðtn ;tnþ1 ;L2 ðXÞÞ þ kvt k2L1 ðtn ;tnþ1 ;L2 ðXÞÞ kvðtnþ1 Þk2L1 ðXÞ þ CReDtkrvðtnþ1 Þk21 ken k2 :

ð3:11Þ

8

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

Next, from 3.4, 3.5 and 3.6 it follows that ðe e nþ1  enþ1 Þ ¼ Dtrpðtnþ1 Þ;

r  ee nþ1 ¼ 0

and so multiplying by e e nþ1 and integrating gives kee nþ1 k 6 kenþ1 k. Combining this with (3.11) and dropping a positive term on the left hand side provides the estimate   1 2 2 2 2 kenþ1 k  ken k Þ þ kenþ1  e e n k þ Re1 Dtkrenþ1 k 2   2 2 2 2 6 2Dt2 krqðtnþ1 Þk þ CDt3 Re kvtt kL1 ðtn ;tnþ1 ;L2 ðXÞÞ þ kvt kL1 ðtn ;tnþ1 ;L2 ðXÞÞ kvðtnþ1 ÞkL1 ðXÞ þ CReDtkrvðtnþ1 Þk21 ken k2 : Summing over time steps, and noting that e0 ¼ 0, produces keM k2 þ Re1 Dt

M X n¼1

krenþ1 k2 6 CReDt

M X

krvðtn Þk21 ken1 k2

n¼1

  2 2 2 þ CTReDt2 kvtt kL1 ð0;T ;L2 ðXÞÞ þ kvt kL1 ð0;T ;L2 ðXÞÞ kvðtnþ1 ÞkL1 ð0;T ;L1 ðXÞÞ þ 2Dt2

M X 2 krqðtn Þk : n¼1

Finally, applying the discrete Gronwall lemma from [19] with the assumption that v 2 L1 ð0; T ; L1 ðXÞÞ, we 2 get the stated result for all Dt > 0 (the result is independent of Dt because there is no keM k term on the right hand side after summing, a case which is discussed in [19]), ! M M   X X 2 2 2 2 2 keM k þ Re1 Dt krenþ1 k 6 C TReDt2 kvtt kL1 ð0;T ;L2 ðXÞÞ þ kvt kL1 ð0;T ;L2 ðXÞÞ þ Dt2 krqðtn Þk : n¼1

n¼1 1

2

1

1

Finally, the assumptions that vt and vtt are in L ð0; T ; L ðXÞÞ, and q 2 L ð0; T ; H ðXÞÞ, produces the stated result. The proofs for the grad-div stabilized methods work in a similar manner, and thus we present them next in a shortened version. To simplify them further, we apply the smoothness assumptions of the true solution as terms appear, instead of waiting until the last step as we did in the proof above. For the standard grad-div stabilized Step 1 + Step 2 result, subtract (3.2) from the NSE at time tnþ1 to get 1 nþ1 e n Þ þ en  rvðtnþ1 Þ þ un  renþ1  Re1 Denþ1  crðr  enþ1 Þ ðe  e Dt ¼ rqðtnþ1 Þ þ Dtðvtt ðtn Þ þ vt ðtn Þ  rvðtnþ1 ÞÞ: Multiplying by enþ1 and integrating over the domain, and using analysis similar to that performed above yields  Re1 1  nþ1 2 2 2 2 2 ke k  kee n k þ kenþ1  e krenþ1 k þ ckr  enþ1 k enk þ 2Dt 2 1 6 CReken k2 þ CDt2 Re þ Dtkrqðtnþ1 Þk2 þ kenþ1  ee n k2 : ð3:12Þ 4Dt From here, a similar analysis as in the unstabilized case finishes the proof. Lastly, for the case of sparse grad-div stabilized Step 1 + Step 2, add cgðvðtnþ1 ÞÞ to both sides of the NSE, then subtract the discrete scheme to get 1 nþ1 ðe  e e n Þ þ en  rvðtnþ1 Þ þ un  renþ1  Re1 Denþ1 þ cgðenþ1 Þ Dt ¼ rqðtnþ1 Þ þ cgðvðtnþ1 Þ þ Dtðvtt ðtn Þ þ vt ðtn Þ  rvðtnþ1 ÞÞ:

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

9

Since r  vðtnþ1 Þ ¼ 0, the first two terms on the right hand side combine to form rqðtnþ1 Þ þ cgðvðtnþ1 Þ ¼ rðqðtnþ1 Þ  v1x ðtnþ1 ÞÞ: Making this substitution, then multiplying by enþ1 and integrating over the domain, using properties of the g operator from Lemma 2.1 and using an analysis similar to the unstabilized Step 1 case provides  Re1 1  nþ1 2 2 2 2 2 ke k  ke krenþ1 k þ ckr  enþ1 k e n k þ kenþ1  ee n k þ 2Dt 2 1 6 CReken k2 þ CDt2 Re þ Dtkrðq  v1x Þðtnþ1 Þk2 þ kenþ1  e ð3:13Þ e n k2 : 4Dt From here, analysis similar to the unstabilized case finishes the proof, with the main difference being that we must additionally assume here that v 2 L1 ð0; T ; H 2 ðXÞÞ. h We now provide two numerical experiments which show the positive impact of the stabilization. 3.2. Numerical experiment 1: error comparison for projection methods with standard, sparse, and no grad-div stabilization Our first numerical experiment compares the error for the unstabilized, standard and sparse grad-div stabilized projection methods for a test problem with a known analytical solution. The problem we consider is the 3 Ethier–Steinman exact Navier–Stokes solution from [9] on the domain X ¼ ð0; 1Þ . For chosen parameters a; d and viscosity m, the exact solution of the NSE is 2

u1 ¼ aðeax sinðay þ dzÞ þ eaz cosðax þ dyÞÞemd t ; 2

u2 ¼ aðeay sinðaz þ dxÞ þ eax cosðay þ dzÞÞemd t ; az

ay

u3 ¼ aðe sinðax þ dyÞ þ e cosðaz þ dxÞÞe

md 2 t

;

ð3:14Þ ð3:15Þ ð3:16Þ

2

a 2ax ðe þ e2ay þ e2az þ 2 sinðax þ dyÞ cosðaz þ dxÞeaðyþzÞ þ 2 sinðay þ dzÞ cosðax þ dyÞeaðzþxÞ 2 2 þ 2 sinðaz þ dxÞ cosðay þ dzÞeaðxþyÞ Þemd t :

p¼

ð3:17Þ

We compute approximations to (3.14)–(3.16) and (3.17) with a = d=1, Re :¼ m1 ¼ 100, using no grad-div stabilization (c ¼ 0), standard grad-div stabilization with parameter c ¼ 1, and sparse grad-div stabilization with parameter c ¼ 1. We use a standard direct finite element spatial discretization of Step 1 + Step 2 (see, e.g., [21]) using a uniform tetrahedral mesh with ððP 2 Þ3 ; P 1 Þ Taylor–Hood velocity–pressure elements that provides 107,811 velocity degrees of freedom and 4913 pressure degrees of freedom. We choose a time step of Dt ¼ 0:001, and endtime T ¼ 1. In Fig. 2, plots of the L2 and H 1 Step 1 velocity errors, L2 Step 2 velocity errors, and Step 1 divergence errors of the respective solutions of the projection methods are shown with respect to time tn . We observe that the solutions of the two stabilized methods are significantly more accurate than the solution of the unstabilized method, and that the errors of the two stabilized methods are approximately the same. The divergence error in the stabilized methods’ Step 1 solutions is about 66% compared to the unstabilized Step 1 solution, and the error in velocity (in all norms shown) is cut by over 50%. We also observe from the plots that in the unstabilized case, the divergence error accounted for more than half of the total H 1 error, and thus it is the reduction of this part of the error that reduced the total H 1 error. This last statement is in complete agreement with the a priori estimates in Theorem 3.1. 3.2.1. Efficiency of direct solvers in numerical experiment 1 An interesting and important feature of sparse grad-div stabilization compared to standard grad-div stabilization is that the linear system matrix can be explicitly decoupled, since they take the forms shown in Fig. 1. Hence, it is expected that if direct solvers are applied to the linear systems, then the decoupled system solves

10

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

Fig. 2. Shown above are error versus time for the Ethier–Steinman problem with no stabilization, standard and sparse grad-div stabilization, for (TL) L2 error for Step 1 solution, (TR) L2 error for Step 2 solution, (BL) divergence error for Step 1 solution, and (BR) H 1 error for Step 1 solution.

Table 1 Timings of Matlab’s direct solver ‘backslash’ on the sparse and standard grad-div stabilized systems in numerical experiment 1, on several meshes. For the sparse stabilization, the system decouples into two smaller solves, and we give timings for each, and add them up for a total solve time. h

dof

Solve time for sparse grad-div

Solve time for standard grad-div

1/4 1/6 1/8 1/10 1/12

14,739 46,875 107,811 206,763 352,947

.54 + .13 = .67 3.84 + .86 = 4.70 18.59 + 3.32 = 21.91 72.13 + 12.29 = 84.42 193.88 + 32.49 = 226.37

1.04 10.47 49.82 186.52 546.86

arising from sparse grad-div stabilization should total significantly less than standard grad-div stabilization. In the test below, we verify this is the case. Using the same test problem as numerical experiment 1, we time a direct solver (Matlab’s backslash) on the fifth timestep, using several meshes, and using both standard and sparse grad-div stabilization, and display results in Table 1. Here we observe that the direct solver applied to the subproblems from sparse grad-div stabilization is over twice as fast compared to standard grad-div stabilized system. 3.3. Numerical experiment 2: comparison of solutions for 2d channel flow over a step We now test the proposed methods on the benchmark 2d problem of channel flow over a step. The domain is a 40  10 rectangle, with a 1  1 step that is five units into the channel at the bottom. We prescribe no-slip T conditions on the walls and step, and enforce a parabolic inflow and outflow given by ðyð10  yÞ=25; 0Þ . The correct behavior with f ¼ 0 and m ¼ 1=600 is a smooth velocity profile, with eddies forming and detaching behind the step [23].

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

No stabilization

10

1

5

0 0 10

0.5

10 20 30 Sparse grad-div stabilization, γ = 0.1

40

10

0.5

10 20 30 Standard grad-div stabilization, γ = 0.1

40

0 1

5

0 0

0 1

5

0 0

11

0.5

10

20

30

40

0

Fig. 3. Shown above are plots of velocity streamlines over speed contours for the T = 40 Step 1 solutions of the three methods.

We discretize spatially using the standard direct finite element implementation of the three methods. We use ððP 2 Þ2 ; P 1 Þ Taylor–Hood elements for velocity and pressure, and on a mesh that provides 11,774 velocity degrees of freedom and 1501 pressure degrees of freedom. We choose the stabilization parameter c ¼ 0:1 for both stabilized methods and a timestep Dt ¼ 0:1. We compute each method to an endtime of T = 40. Plots of the velocity solutions for each method at T = 40 are shown in Fig. 3, and we observe that the unstabilized method fails to resolve the new eddy after the first one detaches, while both stabilized methods do correctly predict the formation of the second eddy. Divergence errors were also calculated for the Step 1 solutions for each method at T = 40, and were found to be: kr  unostab k ¼ 0:763 h kr  ustandardgd k ¼ 0:111 h kr  usparsegd k ¼ 0:187 h Thus we see that both stabilized methods compute approximations with significantly better mass conservation which therefore have more physical relevance than the solutions of the unstabilized scheme. For a further comparison, we rerun the problem with the same algorithms, but now with a coarser discretization. We again use ððP 2 Þ2 ; P 1 Þ Taylor–Hood elements, but now we use a timestep of Dt ¼ 0:2, and a mesh that provides 2984 velocity degrees of freedom and 398 pressure degrees of freedom. Results for no grad-div stabilization, standard grad-div stabilization with c ¼ 0:1, and sparse grad-div stabilization with c ¼ 0:1 are shown in Fig. 4. Oscillations are clearly present in the unstabilized solution, however in both stabilized methods the oscillations are removed and the eddies do shed behind the step. We again compute the divergence error at T = 40, and get

12

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

No stabilization (coarse mesh) 10

1

8

0.8

6

0.6

4

0.4

2

0.2

0 0

5

10

15

20

25

30

35

40

0

Sparse grad-div stabilization, γ = 0.1 (coarse mesh) 10

1

8

0.8

6

0.6

4

0.4

2

0.2

0 0

5

10

15

20

25

30

35

40

0

Standard grad-div stabilization, γ = 0.1 (coarse mesh) 10

1

8

0.8

6

0.6

4

0.4

2

0.2

0 0

5

10

15

20

25

30

35

40

0

Fig. 4. Shown above are plots of velocity streamlines over speed contours for the T = 40 Step 1 solutions of the three methods on the coarser mesh.

kr  uhnostab k ¼ 0:849; kr  ustandardgd k ¼ 0:161; h kr  usparsegd k ¼ 0:187: h Hence, again we see that both grad-div stabilization operators reduce the divergence error and provide much better solutions. 4. Efficiency of iterative linear solvers for Step 1 with stabilization In the previous section, we demonstrated that a significant increase in accuracy and physical fidelity arises from using either the standard or the sparse grad-div stabilization in Step 1 of the projection method. However, there is a cost to these stabilizations, which comes in the form of making the linear solves more difficult. As we discussed in Section 2, the matrix that arises from Step 1 is block diagonal, and if standard grad-div stabilization is used, then the matrix becomes block full. However, in the case of sparse grad-div stabilization, the matrix becomes block upper triangular in 2d and will have two blocks that are zero in 3d (see Fig. 1). In this section, we analyze and compare preconditioned iterative solvers for the linear systems to be solved in Step 1 for both the standard and sparse grad-div stabilizations. In a discretized setting, Step 1 requires the solution of linear systems of equations of the type ðA þ cMÞu ¼ b:

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

13

The matrix A 2 R3n3n is common to all three versions and is of block diagonal form, i.e., 0 1 A11 B C A¼@ A22 A A33 with Aii 2 Rnn . Typically, it holds that A11 ¼ A22 ¼ A33 is the discrete convection–diffusion-reaction operator and as such might be non-symmetric and ill-conditioned, but will be non-singular. The parameter c is the same one that occurs in the standard or sparse grad-div stabilizations (3.2) and (3.3), respectively. (In the unstabilized case, one simply sets c ¼ 0.) The matrix M 2 R3n3n is a singular matrix and represents the grad-div stabilization. In order to describe this matrix, we write the discrete divergence operator in its three components BT ¼ ðBT1 ; BT2 ; BT3 Þ 2 Rm3n corresponding to the derivatives in the three spatial directions. In this notation, and given an additional non-singular (and typically diagonal) scaling matrix W 2 Rmm – one may think of the (diagonally lumped) pressure mass matrix – the matrix M has the following form for the stabilized versions: M gd :¼ M grad-div :¼ BW 1 BT ; 0 B1 W 1 BT1 B M sgd :¼ M sparse-gd :¼ @ 0

2B1 W 1 BT2 B2 W 1 BT2

0

B3 W 1 BT2

1 2B1 W 1 BT3 C B2 W 1 BT3 A: B3 W 1 BT3

Since the size of these systems prohibits a direct solver, an investigation of the properties of iterative solvers for the systems ðA þ cMÞu ¼ b with M 2 fM gd ; M sgd g is warranted. In particular, we are interested in whether it is easier (i.e., more efficient with respect to time and storage) to iteratively solve the sparse compared to the standard grad-div stabilized system. 4.1. A simplified model to compare sparse and standard grad-div stabilization We begin with a simple but illustrative example. Let A ¼ I 2 R3n3n (identity), W ¼ I 2 Rmm and Bt ¼ ð1; 1; . . . ; 1Þ 2 R13n . In this simplified setting, the two stabilized systems we need to solve have the respective forms 0 1 I þ cN cN cN B C I þ cN cN Au ¼ b; Agd u ¼ @ cN 0

cN

I þ cN B Asgd u ¼ @ 0 0

cN 2cN I þ cN cN

I þ cN

1 2cN C cN Au ¼ b I þ cN

where N ¼ ðnij Þ; nij ¼ 1 8i; j 2 f1; . . . ; ng denotes the matrix with all entries equal to one. The following Fig. 5 shows the convergence histories for solving these systems by a block (backward) Gauß–Seidel method for various values of stabilization parameter c and fixed problem size 3n ¼ 300. In these experiments, we choose the right hand side b so that the exact solution satisfies u 2 NullðMÞ (In MATLAB notation: the exact solution u 2 R‘ is chosen as uð1 : ‘Þ ¼ linspaceð1; 1; ‘Þ, uð‘ þ 1 : 3  ‘Þ ¼ linspaceð1; 1; 2  ‘Þ). We started the iterations with the zero vector u0 ¼ 0. The convergence histories are similar in terms of qualitative behavior: As c ! 0, the solvers converge rapidly since the system matrix approaches the identity matrix. As c ! 1, the convergence rate approaches 1 as the system matrices become almost singular. Quantitatively, the convergence is faster for the (sparser) Asgd system compared to the (standard) Agd -system. In the following, we will explain the observed convergence behavior through an analysis of the spectra of the error propagation matrices. We begin by collecting a few straightforward results in the following Lemma.

14

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19 12

8 gamma=0.01 gamma=0.1 gamma=1.0 gamma=10.0 gamma=100.0

10

gamma=0.01 gamma=0.1 gamma=1.0 gamma=10.0 gamma=100.0

7

6

8 5

6

4

3 4 2 2 1

0

0

2

4

6

8

10

12

14

16

18

20

0

0

2

4

6

8

10

12

14

16

18

20

Fig. 5. Convergence histories: the x-axis represents the number of iteration steps, the y-axis shows the error energy norm kuk  u kA obtained with the block Gauß–Seidel preconditioner for solving Agd u ¼ b (left) and Asgd u ¼ b (right) for various stabilization parameters c (Agd ; Asgd 2 R300300 ).

Lemma 4.1. Let Bt ¼ ð1; 1; . . . ; 1Þ 2 R13n . It holds that Bt  B ¼ 3n; ðI þ cN Þ1 ¼ I 

1 cN ; 1 þ cn

N  N ¼ nN cN ðI þ cN Þ

1



1

¼ ðI þ cN Þ cN ¼

 1 cN : 1 þ cn

Proof. The first equation, Bt  B ¼ 3n, follows directly from the particular choice of B. Since N is the rank-1 matrix N ¼ 1  1T with 1T ¼ ð1; 1; . . . ; 1Þ 2 R1n , the third equation N  N ¼ 1  |{z} 1T 1 1T ¼ nN is proven. The ¼n

second equation follows from  ðI þ cN Þ I 

1 cN 1 þ cn



  1 1 1 cn 2 cN  c N N ¼I þ 1  ¼ I þ cN  cN ¼ I: 1 þ cn 1 þ cn 1 þ cn 1 þ cn

The remaining equalities follow directly from this second equation.

h

Since (block backward) Gauß–Seidel preconditioners are among the most frequently used preconditioners in the context of grad-div stabilized systems [1,2], we define such block preconditioners and analyze the error propagation matrices in the following lemma. Lemma 4.2. Using block Gauß-Seidel preconditioners 0 B P GS gd :¼ @

I þ cN

cN

cN

1

0

I þ cN

cN

C A;

B P GS sgd :¼ @

I þ cN the error propagation matrices have the form

I þ cN

2cN

2cN

1

I þ cN

cN

C A;

I þ cN

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

0

ð2þcnÞcn ð1þcnÞ3

cn ð1þcnÞ3 cn ð1þcnÞ2

0

15

1

 1 B C 1 0C I  P GS Agd ¼ B gd @  ð1þcnÞ2 A ðcN Þ 1 1  1þcn  1þcn 0 0 1 2cn 0 ð1þcnÞ 0 3  1 B C cn C I  P GS Asgd ¼ B sgd @ 0 ð1þcnÞ2 0 A ðcN Þ; 1 0 0  1þcn with standard matrix Kronecker product . Their spectra are given by ) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ( 2  1 3cn þ 2ðcnÞ cnð3cn þ 4Þ GS r I  P gd Agd ¼ 0;  ðcnÞ ; 3 2ð1 þ cnÞ ( ) 2   GS 1 ðcnÞ Agd ¼ 0; r I  P rs ð1 þ cnÞ2 with algebraic (and geometric) multiplicity equal to one for the non-zero eigenvalues. Proof. The preconditioners can be written as 0 1 0 1 1 1 1 B C B GS 1 1 P GS ¼ I þ ðcN Þ; P ¼ I þ @ A @ gd rs 1

2 1

2

1

C 1 A ðcN Þ: 1

Straightforward but tedious computations show 00 ð2þcnÞcn 1 1 0 0 1 1 1 1 00 ð2þcnÞcn cn cn 0 0 1 1 1 ð1þcnÞ3 ð1þcnÞ3 ð1þcnÞ3 ð1þcnÞ3 BB C C @ @ BB C C cn cn 1 1 0 A ðcN ÞA 1 1 A ðcN ÞA  @@  ð1þcnÞ P GS Iþ 2 gd @@  ð1þcnÞ2 ð1þcnÞ2 0 A ðcN ÞA ¼ ð1þcnÞ2 1 1 1 1 1  1þcn 0  1þcn 0  1þcn  1þcn 0 1    0 0 0 1 GS GS ¼ @ 1 0 0 A ðcN Þ ¼ P GS  A ¼ P I  P A gd gd gd gd gd 1 1 0 and analogously one checks that 00 1 1 2cn 0 ð1þcnÞ 0 3    1 BB C C cn GS BB GS GS C C P sgd  @@ 0 ð1þcnÞ2 0 A ðcN ÞA ¼ P sgd I  P sgd Asgd ; 1 0  1þcn 0 proving the identities for the error propagation matrices. The spectra of these matrices consist of the products of eigenvalues of the two respective factors in the Kronecker product. h In either case, it holds that the spectral radius approaches 1 as c ! 1, leading to very slow convergence. This is illustrated in Fig. 6. An important observation from Fig. 6 is that the spectral radius of the error propagation matrix for sparse grad-div stabilization is always smaller than that of standard grad-div stabilization. This suggests that when solving systems arising in projection methods for the Navier–Stokes equations with a block Gauß–Seidel (preconditioned) method, we can expect to have similar or slightly better convergence properties for sparse grad-div stabilized systems compared to standard grad-div stabilized systems. Since sparse grad-div stabilization leads to systems with significantly fewer nonzero entries, matrix vector multiplications and hence individual iteration steps are less costly compared to the systems arising from the standard grad-div stabilization, leading to overall faster methods. We will illustrate this expectation with the following numerical example.

16

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19 1 sparse grad−div standard grad−div

0.9 0.8

spectral radius

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

2

4

6

8

10 12 gamma*n

14

16

18

20

Fig. 6. Spectral radii of the error propagation matrices of the block Gauß–Seidel method depending on cn (c is the stabilization parameter and 3n is the problem size, i.e., n unknowns per spatial direction).

4.2. Numerical experiment 3: comparison of iterative solvers on a test problem We now test the efficiency of some common iterative solvers on linear systems arising for the stabilized Step 1 ((3.2) or (3.3), respectively). The test problem we choose is identical to our numerical experiment 1: 3 the Ethier–Steinman problem on the unit cube ð0; 1Þ with Re ¼ 100, discretized temporally with the projec3 tion method and spatially with ðP 2 Þ velocity and P 1 pressure finite elements. We use two refinement levels of uniform tetrahedral meshes, the first (coarser) one leading to 107,811 velocity degrees of freedom, and the second (finer) mesh yielding 823,875 velocity degrees of freedom. The linear systems we consider for timings are those that occur at the fifth timestep. Since the true solution is observed to decay in time, the solves become easier as time progresses. Also we did not wish to choose the first few timesteps, since the smoothness in the initial condition might help the solver. On the coarse mesh, the matrices in Step 1 have 4.64e+6 non-zeros when using sparse grad-div and 5.78e+6 non-zeros when using standard grad-div stabilization. On the finer mesh, the respective numbers of non-zeros are 40.6e+6 for sparse and 50.6e+6 for standard grad-div stabilization. Hence, on either mesh, the matrix resulting from the sparse stabilization has about 20% fewer non-zeros than the matrix from the standard stabilization. Plots of their sparsity structures are shown in Fig. 1. We solve the linear systems using a preconditioned Bicgstab method. As a stopping criterion, we use a reduction of the residual by a factor of 1e6. We compare performance of three different types of preconditioners: (1) Jacobi, (2) Gauß–Seidel, and (3) approximate block Gauß–Seidel where we solve the diagonal blocks with varying accuracies. Preconditioners (1) and (2) do not require any set-up time. For preconditioner (3), we computed approximate LU factorizations of the three diagonal blocks Aii þ Bi W 1 BTi using the technique of (domain decomposition based) hierarchical (H-) matrices [14,13,2]. H-matrices provide a means to compute an approximate LU-factorization with almost optimal computational and storage requirements. The efficiency results from using blockwise low-rank approximations in certain (off-diagonal) matrix blocks within a hierarchically structured partition of the matrix into subblocks. The accuracy of such approximations can be controlled (blockwise) through a parameter dH 2 ð0; 1Þ. The smaller dH , the closer are the H-LU factors to the exact LU factors (but the higher are the costs to compute and store the H-LU factors). In our experiments, we used relative accuracies dH 2 f0:1; 0:01g. We also performed tests for dH ¼ 0:001 which did not lead to better results compared to dH ¼ 0:01, indicating that a highly accurate H-LU factorization is already obtained for dH ¼ 0:01. We note that the (approximate) solution of the problems involving the diagonal blocks in the block Gauß–Seidel method through H-LU factorization could be replaced with alternative, possibly more efficient methods. Here, our goal is to illustrate the advantage of investing into such a solver compared to the simpler Jacobi or (pointwise) Gauß–Seidel preconditioners.

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

17

Table 2 Iteration counts and timings for iterative solvers applied to the sparse and standard grad-div stabilized systems in numerical experiment 3. Preconditioner

Standard grad-div Iterations

Coarse mesh results (107,811 dof) HLU, dh = 1e1 HLU, dh = 1e2 Gauß–Seidel Jacobi

8 6 15 26

Fine mesh results (823,875 dof) HLU, dh = 1e1 HLU, dh = 1e2 Gauß–Seidel Jacobi

14 7 23 45

Sparse grad-div Time 0.91 0.92 0.51 0.45 15.1 10.7 6.8 6.7

Iterations 8 5 13 27 11 6 21 44

Time 0.88 0.76 0.37 0.41 11.5 9.5 5.3 5.6

In Table 2, we show iteration counts and times for each of these solvers applied to the linear systems arising from Step 1 of the standard and sparse stabilized projection methods. In all tests, we used the stabilization parameter c ¼ 1. We observe that the iteration counts are the same or slightly better for the systems arising from sparse grad-div stabilization, and the iteration times are reduced on average by about 20% for sparse versus standard grad-div stabilization due to the sparser matrix structure. As expected, the block Gauß–Seidel preconditioner requires fewer (butmore expensive) iteration steps than the pointwise Gauß–Seidel or Jacobi preconditioners. However, comparing the iteration counts for the two problem sizes, we observe that only the block Gauß–Seidel method using a highly accurate H-LU factorization (with dH =1e2) to solve for the diagonal blocks leads to an (almost) optimal method where the number of required iteration steps remains (almost) constant as we increase the problem size. These results are in agreement with the analysis obtained for the simplified model problem in the previous section 4.1 which suggested that iterative solvers for linear systems with sparse grad-div stabilization will converge with slightly fewer iterations than those for standard grad-div systems. 5. Conclusions and future work We analyzed both standard and sparse grad-div stabilized projection methods, and proved that they offer an improvement in the divergence part of the H 1 ðXÞ spatial error. Moreover, we provided two numerical experiments which demonstrated the improvement in solution accuracy that can be obtained by using either one of these stabilizations. The convergence behavior of iterative linear solvers applied to the discrete linear systems of equations resulting from both standard and sparse grad-div stabilized Step 1 was also studied herein. We found both analytically and numerically that the convergence properties of preconditioned bicgstab using some common preconditioners are slightly better in terms of iteration counts for sparse grad-div stabilized systems, and significantly better in terms of overall iteration time due to the sparsity advantage of sparse grad-div stabilization. An interesting followup study would be to study optimal stabilization parameters in this context. Some work on optimal parameter for standard grad-div stabilization has been done recently for Stokes [20], nothing has yet been done for sparse grad-div stabilization. Furthermore, although it appears extensions of the optimal parameter work to coupled schemes for Navier–Stokes is relatively straight-forward, there will certainly be more technical details to work out if projection method timestepping methods are used. An important future direction is to study grad-div stabilization, standard and sparse, for turbulent and higher Reynolds number flows, as there appears to be little in the literature on this topic both for projection and coupled timestepping schemes. For these important flows, inaccuracies will still arise from poor mass conservation interacting with large/complex pressures, but at the same time there will be instabilities arising from inertial effects. Some important studies related to this are done in [27,24]. In [24], it is found that pointwise

18

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

divergence free elements (Scott–Vogelius) are more accurate than Taylor–Hood for a two dimensional high Reynolds number flow in a star-shaped domain. Since it is proven in [5] that grad-div stabilized Taylor–Hood solutions of Navier–Stokes equations converge to that of Scott–Vogelius as c ! 1, perhaps a finite c is enough for better accuracy, and it is likely! that the optimal c would be finite. In [27], Lube and Ro¨he studied turbulent flows using grad-div stabilization, and found it did not have much effect when the Smagorinsky-type models were being used. Since Smagorinsky damps the entire gradient, it is possible that its implicit damping of the divergence prevents any addition grad-div stabilization from having an effect. Other types of turbulence models that use indicators to detect laminar/smooth flows (e.g., that employ Vreman’s laminar flow detector [36]) would likely still benefit from grad-div stabilization if they contain high speed shear flows, since they would not apply the stabilization in these areas. Hence, a detailed study in this direction is certainly warranted. A related topic that deserves future study is the connection between grad-div stabilization and filtering based models for large eddy simulation. For variational multiscale methods, it is known from [29] that there is a strong connection to grad-div stabilization, as grad-div can be considered as a variational multiscale method in a sense. Such a connection is seemingly not known for filter based models, however. It is shown in [28] for an NS-a type model that grad-div stabilization can play an important role, however no turbulent computations were done in that paper. Yet another interesting future work could be the study of continuous interior penalty methods, such as that studied by Burman and Hansbo in [4], in the context of projection methods. This stabilization is related to grad-div in that it gives control over the divergence in the L2 sense, and hence it is possible that it could also have a positive effect in the projection method setting. References [1] M. Benzi, M.A. Olshanskii, Z. Wang, Modified augmented Lagrangian preconditioners for the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 66 (4) (2011) 486–508. [2] S. Bo¨rm, S. Le Borne, H-LU factorization in preconditioners for augmented Lagrangian and grad-div stabilized saddle point systems, Int. J. Numer. Methods Fluids 68 (2012) 83–98. [3] A. Brooks, T. Hughes, Steamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 32 (1982) 199–259. [4] E. Burman, P. Hansbo, Edge stabilization for the generalized Stokes problem: a continuous interior penalty method, Comput. Methods Appl. Mech. Eng. 195 (2006) 2393–2410. [5] M. Case, V. Ervin, A. Linke, L. Rebholz, A connection between Scott–Vogelius elements and grad-div stabilization, SIAM J. Numer. Anal. 49 (4) (2011) 1461–1481. [6] A.J. Chorin, Numerical solution for the Navier–Stokes equations, Math. Comput. 22 (1968) 745–762. [7] E. D’Agnillo, L. Rebholz, On the enforcement of discrete mass conservation in incompressible flow simulations with continuous velocity approximation. Recent advances in scientific computing and applications, in: Jichun Li, Eric Macharro, Hongtao Yang (Eds.), Proceedings of the 8th International Conference on Scientific Computing and Applications, AMS Contemporary Mathematics, vol. 586, 2013, pp. 143–152. [8] O. Dorok, W. Grambow, L. Tobiska, Aspects of finite element discretizations for solving the Boussinesq approximation of the Navier–Stokes equations, in: F.-K. Hebeker, R. Rannacher, G. Wittum (Eds.), Proceedings of the International Workshop held at Heidelberg, October 1993, Notes on Numerical Fluid Mechanics: Numerical Methods for the Navier-Stokes Equations, vol. 47, 1994, pp. 50–61. [9] C. Ethier, D. Steinman, Exact fully 3d Navier–Stokes solutions for benchmarking, Int. J. Numer. Methods Fluids 19 (5) (1994) 369– 375. [10] L. Franca, Incompressible flows based upon stabilized methods, in: Proceedings of the International Workshop on Numerical Methods for the Navier–Stokes equations, Heidelberg, October, 1993, Notes on Numerical Fluid Mechanics, 1993. [11] K. Galvin, A. Linke, L. Rebholz, N. Wilson, Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection, Comput. Methods Appl. Mech. Eng. 237 (2012) 166–176. [12] V. Girault, P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, 1986. [13] L. Grasedyck, R. Kriemann, S. Le Borne, Domain decomposition based H-LU preconditioning, Numer. Math. 112 (2009) 565–600. [14] L. Grasedyck, S. Le Borne, H-matrix preconditioners in convection-dominated problems, SIAM J. Math. Anal. 27 (2006) 1172–1183. [15] P. Gresho, R. Sani, Incompressible Flow and the Finite Element Method, vol. 2, Wiley, 1998. [16] J. Guermond, Un resultat de convergence d’ordre deux en temps pour l’approximation des equations de Navier–Stokes par une technique de projection incrementale, Math. Modell. Numer. Anal. 33 (1) (1999) 169–189. [17] J. Guermond, P. Minev, J. Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Eng. 195 (2006) 6011–6045.

A.L. Bowers et al. / Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19

19

[18] J.-L. Guermond, Some practical implementations of projection methods for Navier–Stokes equations, M2AN 30 (1996) 637–667. [19] J. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier–Stokes problem. Part IV: Error analysis for the second order time discretization, SIAM J. Numer. Anal. 2 (1990) 353–384. [20] E. Jenkins, V. John, A. Linke, L. Rebholz, On the parameter choice in grad-div stabilization for incompressible flow problems, Adv. Comput. Math. (2013) (to appear). [21] W. Layton, An Introduction to the Numerical Analysis of Viscous Incompressible Flows, SIAM, Philadelphia, 2008. [22] W. Layton, C. Manica, M. Neda, M.A. Olshanskii, L. Rebholz, On the accuracy of the rotation form in simulations of the Navier– Stokes equations, J. Comput. Phys. 228 (9) (2009) 3433–3447. [23] W. Layton, C. Manica, M. Neda, L. Rebholz, Numerical analysis and computational testing of a high accuracy Leray-deconvolution model of turbulence, Numer. Methods Partial Differ. Equ. 24 (2) (2008) 555–582. [24] A. Linke, Collision in a cross-shaped domain – a steady 2d Navier–Stokes example demonstrating the importance of mass conservation in CFD, Comput. Methods Appl. Mech. Eng. 198 (41–44) (2009) 3278–3286. [25] A. Linke, M. Neilan, L. Rebholz, N. Wilson, Improving efficiency of coupled schemes for Navier–Stokes equations by a connection to grad-div stabilized projection methods, 2013, Submitted for publication. [26] A. Linke, L. Rebholz, On a reduced sparsity stabilization of grad-div type for incompressible flow problems, Comput. Methods Appl. Mech. Eng. 261 (2013) 142–153. [27] G. Lube, L. Ro¨he, Analysis of a variational multiscale method for large-eddy simulation and its application to homogeneous isotropic turbulence, Comput. Methods Appl. Mech. Eng. 199 (37–40) (2010) 2331–2342. [28] C. Manica, M. Neda, M.A. Olshanskii, L. Rebholz, Enabling accuracy of Navier–Stokes-alpha through deconvolution and enhanced stability, M2AN Math. Modell. Numer. Anal. 45 (2011) 277–308. [29] M. Olshanskii, G. Lube, T. Heister, J. Loewe, Grad-div stabilization and subgrid pressure models for the incompressible Navier– Stokes equations, Comput. Methods Appl. Mech. Eng. 198 (2009) 3975–3988. [30] M.A. Olshanskii, A low order Galerkin finite element method for the Navier–Stokes equations of steady incompressible flow: a stabilization issue and iterative methods, Comput. Methods Appl. Mech. Eng. 191 (2002) 5515–5536. [31] M.A. Olshanskii, A. Reusken, Grad-div stabilization for the Stokes equations, Math. Comput. 73 (2004) 1699–1718. [32] A. Prohl, Projection and Quasi-Compressibility Methods for Solving the Incompressible Navier–Stokes Equations, Teubner-Verlag, Stuttgart, 1997. [33] A. Prohl, On pressure approximation via projection methods for nonstationary incompressible Navier–Stokes equations, SIAM J. Numer. Anal. 47 (1) (2009) 158–180. [34] R. Rannacher, On Chorin’s projection method for the incompressible Navier–Stokes equations, The Navier–Stokes Equations II – Theory and Numerical Methods, vol. 1530, Springer, 1992. [35] R. Temam, Sur l’approximation de la solution des equations de Navier–Stokes par la methode des pas fractionnaires (II), Arch. Ration. Mech. Anal. 33 (1969) 377–385. [36] A.W. Vreman, An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications, Phys. Fluids 16 (2004) 3670–3681.