Time and space parallelization of the navier-stokes equations

Time and space parallelization of the navier-stokes equations

Isidoro I. Albarreal Núñez^I; M. Carmen Calzada Canalejo^II; José Luis Cruz Soto^II; Enrique Fernández Cara^I; José R. Galo Sánchez^II; Mercedes Marín Beltrán^II

^IDepartamento E.D.A.N., Universidad de Sevilla, Aptdo. 1160, E-41080 Sevilla

^IIDepartamento de Informática y Análisis Numérico, Universidad Córdoba, Campus de Rabanales, Ed. C2-3, E-14071 Córdoba

E-mails: [e-mail]/ [e-mail]/ [e-mail]/ [e-mail]/ [e-mail]/ [e-mail]

ABSTRACT

In this paper, we will be mainly concerned with a parallel algorithm (in time and space) which is used to solve the incompressible Navier-Stokes problem. This relies on two main ideas: (a) a splitting of the main differential operator which permits to consider independently the most important difficulties (nonlinearity and incompressibility) and (b) the approximation of the resulting stationary problems by a family of second-order one-dimensional linear systems. The same strategy can be applied to two-dimensional and three-dimensional problems and involves the same level of difficulty. It can be also useful for the solution of other more complicate systems like Boussinesq or turbulence models. The behavior of the method is illustrated with some numerical experiments.

Mathematical subject classification: 65M06, 35A35, 68W10.

Key words: Navier-Stokes equations, numerical solution, parallel algorithms.

1 Introduction

We will consider here a numerical method for solving the incompressible, time-dependent, Navier-Stokes equations. These equations can be used to model the behavior of a homogeneous, incompressible, viscous newtonian fluid. When we impose Dirichlet conditions on the velocity field, the problem reads

Here, W Ì

^d is a bounded regular domain (d = 2 or 3), u = u(x,t) is the velocity field, p = p(x,t) is the pressure, n> 0 is the kinematic viscosity (a positive constant) and f = f(x,t) is the density function of a field of external forces. For simplicity, we have assumed in (1) that the fluid has unit mass density.

Concerning the solution to (1), it is well known that we can only expect to get numerical approximations. However, it is also well understood nowadays that this is a very difficult task.

A good strategy seems to be the use of parallel computers. Of course, in order to optimize their efficiency, one has to design appropriate algorithms (in general terms, the next generation of processors is expected to multiply the speed of computation by a factor 10; at the same time, new forthcoming parallel algorithms are expected to produce an increase of a factor 100, see [22]). However, up to now, parallelization has been performed almost always at the lowest level, when the task has been reduced to the solution of finite-dimensional linear systems with probably many unknowns.

The goal of this work is to propose a different method which relies on parallelization at the highest possible level and tries to reduce as much as possible the computer time by using a large number of processors. It will be seen that this method leads to difficulties essentially of the same kind in the 2D and 3D settings. The seminal ideas for this approach can be found in [19].

In this paper, we will only consider low or moderate Reynolds numbers (respectively up to 4000 and 1000 in 2D and 3D problems). Recall that the Reynolds number of (1) is given by Re = UL/n, where U and L are characteristic values of the velocity field modulus and the length, respectively. For higher Re, we would need more subtle arguments and methods. The design of appropriate techniques, similar to those in this paper, for the numerical solution of large Reynolds number problems (1) will be the subject of future work.

As usual, the approximation of (1) is performed in two steps. We first discretize in the time variable and, then, we solve numerically the resulting stationary problems by introducing a spatial approximation.

At both steps, we can use a plenty of methods. Among all them, let us mention viscosity splitting methods and in particular q-scheme fractional schemes for the approximation in time and finite element and finite difference methods for the approximation in space. A detailed analysis of the behavior of these and many other methods can be found in [15].

Our interest has focused on the design and analysis of numerical schemes relying on two main ideas: (a) to split or separate in parallel the most important difficulties (nonlinearity and incompressibility) and (b) to approximate the resulting stationary problems by a (large) family of second-order, completely independent, one-dimensional linear systems.

At the 1D level, it will be then easy an adequate to apply finite difference techniques to produce good approximations. In this way, the solution strategy will make possible a very high level or parallelization.

The research described in this paper is a small part of a much larger project concerning parallelization and nonlinear partial differential systems. Up to now, this has led to some publications and PhD Theses. See for instance [1]-[3], [5], [6], [8]-[12]. However, the numerical techniques we present below can be useful for solving many different problems: linear and semilinear elliptic and parabolic systems with nonlinear boundary conditions, Boussinesq systems, one-equation and two-equation turbulence models, fluid-solid interaction models, fully nonlinear equations of the Monge-Ampère kind, etc.

2 The algorithm

Before recalling the formulation of the algorithm, let us introduce some notation:

• J(W) = { j Î (W)^d: Ñ = j in W}; H (resp. V) is the closure of J(W) in the space L²(W)^d (resp. (W)^d). Thus, H (resp. V) is a Hilbert space for the scalar product of L² (W)^d (resp. (W)^d), which will be denoted by (· ,·) (resp. ((· ,·))). The associated norm will be denoted by |·| (resp. ||·||).

• V' is the dual space of V; á· ,·ñ denotes the duality pairing between V' and V.

• We also introduce the trilinear forms

  b(· ,· ,·) and

  (· ,· ,·), with

for any u, v, w Î H¹(W)^d (here, the usual summation convention is used).

The following properties of V and H are well known:

We can now give a rigorous formulation of the unsteady Navier-Stokes problem in W×(0,T):

In (2), u₀Î H and f Î L²(0,T;L²(W)^d). It is well known that (2) possesses at least one solution which is furthermore unique if d = 2. If u is a solution, then u solves, together with some scalar distribution p, the Navier-Stokes equations (1) (for instance, see [18]). One also has

and u|_{t = 0} = u₀ in an appropriate sense.

Notice that

Also, (u,v,v) = 0 for all v Î (W)^d (even when div u ¹ 0). Consequently, the variational evolution equation in (2) can also be written in terms of (· ,· ,·) and this gives the following equivalent formulation:

We are now going to indicate how to approximate in time. Let us divide the interval [0,T] in M subintervals of length k (k = T/M) and let us assume that the parameters s Î (0,1], q,m Î[0,1] and a,b > 0 are given.

We first put

Then, for given m > 0 and u^mÎ (W)^d (an approximation of u at time t_m = mk ), we compute u^m+a, u^m+b and then u^m+1 as follows. We first solve in parallel the elliptic systems

PROBLEM (BP) (Burgers)

and

PROBLEM (SP) (Stokes)

Then, we set

In (5) and (6), f^m+a y f^m+b are appropriate approximations of f. For instance, we can make the following choice:

In (5), several different definitions of u* and u** are possible. Thus, it seems natural to put

u** = au^m+a + (1 - a)u^m

for some a. Actually, the choice of u** is crucial when one tries to establish ''a priori'' estimates of the numerical solutions. On the other hand, the particular u^* we use determines the degree of linearity we conserve in (5).

Using more or less standard arguments, we can deduce existence and uniqueness results for (5) and (6), at least when n is not too small (see for instance [17] and [18]).

In the previous works [5] and [6], we have presented theoretical and numerical results obtained for some parallel schemes of the kind (4)-(7). There, parallelization was performed only at the time approximation level and the stationary problems (BP) and (SP) were solved with finite element techniques. In this work, we are going to extend the parallelization procedure to all the variables.

To this end, we will apply simultaneous directions implicit (SDI) techniques to the previous stationary problems. Thus, let us denote by W_h a finite dimensional Hilbert space determined by a second-order finite difference approximation of (W)^d (h is a parameter that allows to identify the mesh; of course, we pass from the finite-dimensional to the infinite-dimensional problems by letting h ® 0). Let V_h Î W_h be the subspace formed by the functions in W_h with vanishing discrete divergence (see [21] for several possible W_h and V_h and the associated definitions of the discrete divergence).

Then, the spatial approximation of (4)-(7) is the following:

First, is the orthogonal projection of u₀ on W_h for the L² scalar product, i.e.

Then, for any given m > 0 and Î W_h , we compute , and then as follows. We first solve in parallel two independent problems:

Then, we put

Again, we have several possibilities for the choice of and . For instance, we can take

= s + (1 - s)

= or = , etc.

The existence and uniqueness of a solution to (10) is an immediate consequence of Lax-Milgram's lemma. The existence of a solution to (9) is easily implied by Brouwer's fixed point theorem (see for instance [18]).

In [6] and [8], we have deduced convergence and stability results for the completely discretized scheme (in time and space). In some particular situations, we have also deduced error estimates, cf. [2] and [1].

3 A convergence result

We recall in this section a convergence-stability result for the previous numerical method. In the sequel, C denotes a generic positive constant only depending on the data W, T, n, u₀ and f and, possibly, the parameters s, q and m.

There exist optimal quantities d₀ , S(h) and S₁(h) such that

and

More precisely, we have d₀ = 2 where is the smallest size of W in the directions x₁,...,x_d,

(see [21]).

For each time step k and each h > 0, we introduce the functions u_{k h} , v_{k h} , w_{k h} , z_{k h} , _{k h} , _{k h} and _{k h} , given as follows:

Theorem 1. Assume that s Î (,1], a+b = 2, = s + (1-s)and, for instance, = in (9). There exist constants K₀and K₁ , only depending on W, |u₀|, ||f||_L²(0,T;H) , n and s, such that, whenever

we have:

1. There exist subsequences u_k'h' ,... ,

_k'h' that converge strongly in the space L²(0,T;L²(W)^d), weakly in L² (0,T;(W)^d) and also weakly-* in L^¥(0,T;L²(W)^d) to the same function u.

2. The limit of any such subsequence is a solution of (2).

3. Consequently, when d = 2, the whole sequences u_{k h} ,... ,

_{k h} converge (in the above sense) to the unique solution of (2).

4. Finally, if d = 2 and k and h satisfy

we also have strong convergence of the whole sequences in the space L²(0,T;(W) ^d).

The proof of this result is given in [2]. Notice that (15) can be viewed as a stability condition. It means in practice that, for any small h, k cannot be too large.

Remark 1. From the proof of this result we see that, in order to get stability with restrictions as weak as possible, it seems preferable to take q = 1, which is equivalent to leave the whole nonlinear term in (BP). On the other hand, as expectable, we see that the choice of m has no influence on (15). A more detailed analysis shows that the best parameters conserving stability and low computational cost are those satisfying

(see [9]), where h₀ = min_{1 < i < d} h_i . A remarkable fact is that, for small n (the most interesting situation from the realistic viewpoint), the stability requirements (17) become weak.

4 The numerical solution in practice

After time discretization, we must solve independent stationary problems of two kinds:

• Burgers-like problems (BP) that can be linear or not, depending on the definition of u*.

• Generalized (linear) Stokes problems (SP).

We are now going to indicate the way these problems are solved in practice. The main idea is to reduce the task to the solution of a family of Poisson problems. Then, as already said, we will apply SDI techniques. As a result, we will only find (many) independent 1D differential problems.

4.1The numerical solution of Burgers problems

The goal is to solve numerically a system of the kind

where

When u* = v, we find a linear elliptic system. Contrarily, when u* = u, (18) is nonlinear. In both cases, (18) is solved applying an iterative fixed point algorithm leading to standard Poisson equations completed with Dirichlet conditions.

To this end, let us rewrite (18) in the form

with

G(u, v) = F - 2q(u*· Ñ)(su + (1 - s)v) - q(divu*)((1 - s)v + su).

We first take u⁰ = v (in practice, the velocity field we know from the previous time step). Then, for any n > 0, we compute the solution uⁿ⁺¹ to the linear system

and we iterate until the desired precision is reached.

Observe that (20) is a set of d independent scalar Poisson-Dirichlet problems, for which the unknowns are u_i for i = 1,...,d. Consequently, they can be solved in parallel. To each of these problems, SDI techniques will be applied (see subsection 4.3). Thus, we see that at least conceptually (18) reduces to a family of 1D differential problems, many of them independent, all them leading to similar numerical difficulties.

When (18) is nonlinear and the considered Reynolds number is large, the previous fixed point argument does not suffice. In that case, more sophisticated methods are required based for instance on least square reformulations and Newton-like or conjugate gradient algorithms (see [14] for a complete analysis). However, for the low or moderate Reynolds numbers considered in this paper, it is sufficient to argue as before. In fact, the numerical experiments show that good convergence is attained after very few iterates (see [8] for more details).

4.2The solution of the generalized Stokes problems

Now, we deal with the linear problem

where

This generalized Stokes problem has been solved using a conjugate gradient algorithm adapted from the methods in [7]. Its complete description has been given in [9].

This procedure reduces (again) the task to the solution of Poisson-Dirichlet problems of the kind (22) (see below). In practice, in order to improve its behavior, we have to incorporate preconditioners. For more details, see [7] and [8].

Remark 2. At each time step, when we solve (SP), we get in particular a numerical approximation to the pressure p(t_m+1).

4.3The solution to the Poisson-Dirichlet problems with SDI methods

We include here a brief description of the parallel method originally proposed in [19] and later used and analyzed in [8] and [11] to solve the previous Poisson-Dirichlet problems.

Thus, let us consider the system

where f and h are given. Let us write L in the form L = L₁ +...+ L_d , where

We start with an arbitrary U⁰ satisfying

Then, for any m > 0, U^m+1 is found from U^m as follows. First, we compute U^m+1,1,...,U^m+1,d by solving in parallel the ''one-dimensional'' problems

completed with appropriate boundary conditions deduced from (22). Secondly, we set

Here, t and w (that can depend on m) are parameters that must be determined in order to improve convergence properties. For instance, when d = 2, the iterative method (23)-(24) leads to the following problems:

completed with Dirichlet (two-point) boundary conditions.

We see that, for each fixed x₂ , the unknown function

U^m+1,1(×, x₂)

solves in (25) an ordinary differential equation (in (25), x₂ is a parameter). In practice, it seems reasonable to fix a finite set of x₂ values and, then, solve numerically the corresponding problems (25). In this way, we will obtain approximations to the values of the unknown U^m+1,1 in a finite set of grid points (, ).

Of course, similar things can be said for (26). Since both unknowns U^m+1,1 and U^m+1,2 have to be used in (23) for the computation of U^m+1, it is desirable to use the same coordinates and in (25) and (26). Accordingly, we are led to use rectangular grids.

It is important to emphasize that the level of difficulty is not increased in the case d = 3 since, at the end, the task is reduced to the numerical solution of (many) one-dimensional problems like (25) and (26).

Several slight generalizations of (23)-(24) and a convergence analysis have been given in [11].

5 The behavior of the SDI method

We are now going to illustrate the behavior of the SDI method when it is applied to the following test problems:

• The two-dimensional problem (22) in W = (0,1)×(0,1) with a = 0, b = 1 and f º 2. The exact solution is

  u = sinh(px₁) sin(px₂)+x₁(1 - x₁).

• A similar three-dimensional problem in W = (0,1)×(0,1)×(0,1), with a = 0, b = 1 and

  f º 4. Now, the exact solution is

  u = sinh(p x₁) sin(p x₂)+sinh(p x₁) sin(p x₃)+x₁(1-x₁)+x₃(1-x₃) .

Of course, in both cases the imposed boundary conditions are those fitted by u.

The previous algorithm has been implemented in a SGI Origin 2000 computer with 8 processors, using the parallel computing model of OpenMP. In order to measure the performance of the parallel algorithm, let us introduce two parameters: the speed-up S_p, defined by

and the efficiency h = S_p/p. We have obtained results for different meshes and for 2, 4, 6 and 8 processors. In Fig. 1-4, the associated speed-up's and efficiencies are shown.

The observed behavior is very similar in dimensions 2 and 3. For coarse grids, parallelization does not improve the speed-up. In fact, in the case of a 65 ×65 mesh, the results when using 8 processors are worse than those provided by a sequential method. This can be justified because, in this case, the cost of initializing the processors is, probably, greater than the benefit (the computational work of each processor is too small). On the contrary, when the number of nodes is high, the speed-up and the efficiency increase, and we obtain an efficiency of 0.6 for 8 processors, that we think reasonable. The results are analogous for other tests in non rectangular domains. For more details, see [8].

6 Two-dimensional numerical experiments

We are now going to present some numerical results concerning the Navier-Stokes problem (1). We begin with the so called 2D no-flow test. This is related to the motion of a viscous 2D fluid under the action of gravitational forces of the form f = (0,-1). The considered domain W and the pressure isolines are displayed in Fig. 5.

The spatial approximation has been determined by a regular mesh of step h = 0.05. The time step has been k = 0.01. The following parameters were chosen: s = 0.51, q = 0.5 and a = 0.5. We have started the computations from zero velocity field and pressure at time t = 0.

The computed isobars, in accordance with the theoretical prediction, are horizontal straight lines. The exact error for the pressure (in norm L^¥) is less than 3.4 ×10^-7. The L^¥ norm of the computed velocity is less than 3.28 ×10^-11.

We consider now a squared cavity of unit side filled by a fluid. We assume that the upper wall slips with constant velocity and we try to determine which is the effect of this on the fluid (see [13]).

The results we obtain for various different parameters q at time T = 10 have been presented in Fig. 6 and 7. This test has been performed for a Reynolds number Re = 1000, with spatial mesh size h = 0.01 and time approximation step k = 0.01. We have taken a = 0.5 and s = 0.51.

It can be observed that, as q increases, the occurrence of spurious pressures becomes more important. A detailed analysis of the reasons that lead to this phenomenon has been presented in [8] and [10]. There, a discrete filtering operator that eliminates these undesirable fluctuations has been introduced.

The numerical results corresponding to the particular case q = 1 without and with filtering have been displayed in Fig. 7(a) and Fig. 7(b), respectively.

In order to make a comparison with the results furnished by Ghia [13], Botella and Peyret [4] and Griebel [20], we have fixed a regular spatial approximation with 129 points in each direction. We have taken again k = 0.01, Re = 1000 and a final time T = 30. The computed streamlines and isobars are presented in Fig. 8.

Some computed values of the velocity field and the pressure, together with the deviations of these results from those in [13], [4] and [20] are also given in [9]. We can observe that the proposed algorithm leads to numerical solutions whose behavior is correct.

In a similar way, we have displayed in Fig. 9 the streamlines and isobars corresponding computed for this test for Re = 4000 at time T = 50, using a spatial mesh size h = 0.01 and a time approximation step k = 0.005. Of course, in this case we have had to take a smaller k in order to ensure numerical stability of the Burgers problems.

7 Three-dimensional numerical experiments

We have also considered the no-flow and the cavity tests for a 3D fluid.

For the no-flow test, the results obtained at time T = 3 for h = 0.05 and k = 0.01 are presented in Fig. 10. We have used here the following parameter values: q = 0.5, s = 0.51 and a = 0.5. Again, the algorithm has started from a vanishing velocity field and a vanishing pressure. We have displayed the surfaces P = const. noticing that the exact error, in norm L^¥, is less than 10^-5.

For the 3D cavity test, we have compared our results to those furnished by [16]. For instance, the results obtained for Re = 1000 have been presented in Fig. 11-12, respectively. Good agreement is found.

8 Final comments

The proposed parallel algorithm leads to satisfactory numerical results for both 2D and 3D Navier-Stokes fluids with small or moderate Reynolds numbers. The theoretical results have been confirmed by a set of numerical experiments.

Although we have only reported here Reynolds numbers up to 4000 in 2D tests and up to 1000 in 3D tests, other experiments reveal that, with an appropriate choice of the parameters, we can obtain realistic results for larger Re.

The algorithm we have used is very flexible can be easily adapted to many particular flows. Another positive fact is the incorporation of SDI methods for the solution of Poisson problems, which can be combined with multigrid techniques in view of their smoother properties (see [12]).

Let us finally mention that other boundary conditions can be considered. Indeed, SDI methods can be adapted to the numerical solution of Poisson equations with Neumann and Robin boundary conditions (see [8] and [3]).

[20] M. Griebel and F. Koster. NaSt3DGP. http://wwwwissrech.iam.uni-bonn.de/research/projects/koster/NaSt3DGP/index.htm (20 January 2002).

Received: 9/XII/04. Accepted: 5/IV/05.

#622/04.

• [1] I.I. Albarreal. PhD Thesis. University of Sevilla, Spain (2004).
• [2] I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fernández-Cara, J.R. Galo and M. Marín. Convergence analysis and error estimates for a parallel algorithm for solving the Navier-Stokes equations, Numer. Math., 93 (2002), 201-221.
• [3] I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fernández-Cara, J.R. Galo and M. Marín. A new method for the numerical solution of elliptic and parabolic problems with Neumann boundary conditions, to appear.
• [4] O. Botella and R. Peyret. Benchmark spectral results on the lid-driven cavity flow, Computers & Fluids, 27 (1998), 421-433.
• [5] J.L. Cruz, M.C. Calzada, E. Fernández-Cara and M. Marín. A parallel algorithm for solving the incompressible Navier-Stokes problem Comp. Math. with Appl., 25(9): (1993), 51-58.
• [6] J.L. Cruz, M.C. Calzada, E. Fernández-Cara and M. Marín. A convergence result for a parallel algorithm for solving the Navier-Stokes equations, Comp. Math. with Appl., 35(4): (1998), 71-88.
• [7] E. Dean and R. Glowinski. On some finite element methods for the numerical simulation of incompressible viscous flow, in ''Incompressible computational fluid dynamics'', M.D. Gunzburger and Nicolaides RA eds., Cambridge U.P. (USA), (1993), 17-65.
• [8] J.R. Galo. PhD Thesis. University of Sevilla, Spain (2002).
• [9] J.R. Galo, I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fernández-Cara and M. Marín. A parallel algorithm in time and space for the Navier-Stokes equations, to appear.
• [10] J.R. Galo, I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fernández-Cara and M. Marín. Analysis of a finite difference regularizing scheme with the same grid for the velocity and the pressure, to appear.
• [11] J.R. Galo, I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fernández-Cara and M. Marín. The convergence of the parallel SDI method for elliptic problems, to appear.
• [12] J.R. Galo, M.C. Calzada, J.L. Cruz, M. Marín, I.I. Albarreal and E. Fernández-Cara. The smoothing effect of a simultaneous directions parallel method as applied to Poisson problems, to appear in Numer. Methods for Differ. Eqn.
• [13] U. Ghia, K.N. Ghia and C.T. Shin. High-Resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comp. Physics, 48 (1982), 387-411.
• [14] R. Glowinski. Numerical Methods for Nonlinear Variational Problems (2nd edition), Springer-Verlag, New York (1984).
• [15] R. Glowinski. Finite element methods for incompressible viscous flow Handbook of numerical analysis, Vol. IX, 3-1176, North-Holland, Amsterdam (2003).
• [16] K. Kakuda and H. Kalo. Petrov-Galerkin finite element approach using exponential functions for three-dimensional incompressible viscous flow problems, Finite Element in Fluids (1993), 204-213.
• [17] O.A. Ladyzhenskaya. The Mathematical Theory of Viscous Incompressible Flow Gordon and Breach, New York (1969).
• [18] J.L. Lions. Quelques Méthodes de Résolution des Problčmes aux Limites non Linéaires Dunod, Gauthiers-Villars, Paris (1969).
• [19] T. Lu, P. Neittaanmäki and X.C. Tai. A parallel splitting-up method for partial differential equations and its applications to Navier-Stokes equations, RAIRO Modél. Math. Anal. Numér. 26(6): (1992), 673-708.
• [21] R. Témam. Navier-Stokes Equations (2nd edition), North-Holland, Amsterdam (1984).
• [22] S. Turek. Efficient solvers for incompressible flow problems Springer-Verlag, Berlin (1999).

Publication Dates

• Publication in this collection
  20 Apr 2006
• Date of issue
  Dec 2005

History

• Accepted
  05 Apr 2005
• Received
  09 Dec 2004