•

An Upwind Finite Element Scheme for the Navier-Stokes Equations
and its Domain Decomposition Algorithm( -j- r'f l' . I . A ~ - 7 A }J t.¥~O) ~,Q ~JJt!:ff~H~* A ~ - A c!:: f- O)~~5t~u ~ 7 }l.t -::1'1) ;1:..' A )Thesis for the Degree of Doctor of Science, Hiroshima University
Copyright by Shoichi FUJIMA!v1arch 1997
Department of Mathematics, Faculty of Science, Hiroshima University
1-3-1 Kagamiyama, Higashi-Hiroshima 739, JAPANApproved by Supervisory Committee:
Chairperson: Professor Masahisa T ABA T A *Member: Professor Osamu MIY AMURA **Member: Professor Shinnosuke OHARU** Department of Mathematics, Faculty of Science , Hiroshima University** Department of Physics, Faculty of Science, Hiroshima University

An Upwind Finite Element Scheme for
the Navier-Stokes Equations and its
Domain Decomposition Algorithm

Shoichi FUJIMA

Thesis, Hiroshima University
1~997

Acknowledgements

Cfirst of all, 9 wish to express my sincere gratitude to CProfessor Wasahisa %hata of
9/iroshima CZlnioersify. 9te has encourag~'d me constantly and made oaluahle discussions
throughout the whole period of the preparation for the thesis. CBesides he takes the lahor of thechainnan of the refreeing commiffee ofthe thesis.

9 also wish to express my thanks to memhers of the refreeing commiffee of the thesis,
CProfessor ()samu Wlyamura and CProfessor (5hinnosuke ()haru, who read the thesis while theirtightjoh schedule.
9lpp reciation goes to CProfessor ~~ruo CZlshijima of the CZlnioersify of elecfroCommunications, who had supported me while 9 had been at the CZlnioersify.

9 would like to thank CProfessor Cfatsu!:,(uki 9oakaki, who haoe conlrfhuted to improoe thethesis through helpful discussions.

9 would also like to thank Cj)octor 9tiroshi CXanayama of FU-.fIFACOM Corporation
for application of the upwind scheme to engine Qring problems and gioing me useful infonnations.

9 wish to thank Wr. CYasuji Cfukasawa, who had assisted me in three-dimensional
computation.

Contents

General Introduction12

1

Mixed Finite Element Method for the N avier-Stokes Equations

5

1.1

Introduction......................

5

1.2

Variational principle and the finite element method

61.3Standard mixed finite element method for the Navier-Stokes equations with iso-P2
PI/PI elements . . . . . . . . . . . . . . . .

10

1.4

Consistent discretized pressure Poisson equation

15

1.5

Decomposition into sub-tetrahedra

17

1.6

Conclusion

20

Upwind Finite Element Scheme Based on the Choice of U p- and Down-wind Points 21

2.1

Introduction..........................................

21

2.2

Upwind finite element approximation based on the choice of up- and down-wind points

23

2.3

Search algorithm for the upwind points

27

2.4

2.3.1

Heuristic- and blind-search algorithms for U .

27

2.3.2

Heuristic search algorithm for W . . . . . . .

32

2.3.3

A topology oriented method in heuristic algorithms.

33

Finite element scheme based on the choice of up- and downwind points for the NavierStokes equations

2.5

2.6

. . . . . . . . . . . . . .

34

Stability regions with respect to ex and T .

35

2.5.1

Numerical experiments. . .

35

2.5.2

A theoretical consideration

39~umerical results. . . . . . . .

40

2.6.1

Flows in a cubic cavity.

40

2.6.2

Flows past a circular cylinder

42

11

2.7

3

Conclusion..........................................

Domain Decomposition Algorithm

57
58

3.1

Introduction . . . . . . . . . . . . .

3.2

Two data types for data on the nodal points.

60

3.2.1

Conventional programs and data-parallelization

60

3.2.2

Two data types for data on the nodal points

62

3.2.3

Implementation of the conjugate gradient method

65

3.3

3.4

3.5

58

Mortar element method

66

3.3.1

Review

3.3.2

Mortar element scheme for the Stokes problem

66

71

Domain decomposition algorithm for the upwind finite element scheme

73

3.4.1

Domain decomposition algorithm

73

3.4.2

Numerical results.

76

Conclusion

84

General Conclusion

85

Bibliography

88

11l

General Introd u4ction
In this thesis , numerical methods for flow problems described by the incompressible Navier-Stokes
Equations are considered. The aim is to present an upwind finit e element scheme based on the choiceof up- and downwind points for the equations and to apply a domain decomposition algorithm to thisscheme. );Iain parts have been and will be published in the following papers :• M. Tabata and S. Fujima, Finite-element-analysis of high Reynolds number flows past a circularcylinder, J. Compo Appl. Math., 38 , 4:11-424, 1991.
• S. Fujima, :vI. Tabata and Y. Fukasawa, Extension to three-dimensional problems of the upwind
finite element scheme based on the choice of up- and downwind points, Compo Meth. Appl.

Mech. Engrg. , 112 , 109-131 , 1994.
• S. Fujima, Implementation of mortar element method for flow problems in the primitive variables,
to appear in Int. J. Compo Fluid Dyn .

Upwind finite element schemes for the N avier-Stokes equations

vVith the development of computer performance in speed and storage size, a wide class of flow problems
have become solvable by numerical methods . The finite element method, which has advantages such as
geometrical flexibility and easy treatment of natural boundary conditions, is one of the most powerful
numerical methods for flow problems. Even if the shape of a domain is complicated, one can divide
the domain into a union of finite elements and approximate equations as well as boundary conditions
with good accuracy.
As is well known , however , the standard Galerkin finite element method leads to a numericalinstability for flow problems at high Reynolds numbers. Several improved finite element schemes havebeen proposed[65, 52]. Here we refer to [9 26] for the stream upwind/Petrov-Galerkin scheme [13, 14]
for Galerkin/least squares methods or stabilized methods, [30] for a first-order least squares method ,[37] for a third-order upwind scheme based on the Petrov-Galerkin method , and [63] for a scheme

1using a rotatio nal form.
The idea of choosi ng upwind point(s ) was origina lly introdu ced
approx imatio n by Tabata [52].

into an upwind finite elemen t

Althou gh that is based on the choice of an upwind elemen t,

essenti ally same to choose an upwind point.

it is

The convec tion term is approx imated in first-o rder

accura cy. Bristea u et al. [8] propos ed an approx imatio n based on

the choice of two upwind points, thatis in second -order accura cy. The upwind finite elemen t approx imatio
downw ind points was propos ed in [55] for two-di mensio nal proble
and [17]. The approx imatio n presen ted in this thesis is based onn based on the choice of up- andms, and have been develo ped in [56]the choice of two upwind points andtwo downw ind points.
Combi ning the upwind approx imatio n with the iso-P2 PI/PI

finite elemen ts for the velocit y andthe pressu re[6]' and the proced ure of a consis tent discret ized pressu
finite elemen t scheme for the flow proble rns at high Reyno lds

re Poisso n equatio n[22], an upwind

numbe rs is obtain ed, which has thefollowing advant ages.
(i) It is a kind of finite elemen t scheme and no curvili near coordi
therefo re subdiv ide any part of the domai n, includ ing bound ary

nate systems[66] are used. We canlayers, as freely as we like by virtueof the geome trical flexibility of the finite elemen t metho d.
(ii) It has the potent ial to approx imate the convec tion term inthird-o rder accura cy. Under a suitabl e

subdiv ision of bound ary layers, the convec tion term can beapprox imated in third-o rder accura cy

essentially. Conseq uently, the behavi or of high Reyno lds number flows inbound ary layers can be

captur ed well.
(iii) The deriva tion of the scheme is natura l from the mathe matica
bound ary condit ions such as the ;\Ieum ann condit ion in the pressu
therefo re applica ble to any flow proble ms, inner flows or outerl point of view and no artifici alre Poisso n equati on are used. It is

flows.

(iv) Iso-P2 PI/PI finite elemen ts satisfy the inf-sup condition[20]. Conseq uently, no spurio us pressu remode appear in solving the consis tent discret ized pressu re Poisso
n equation.

Doma in decom positi on algori thm for parall elism

For practic al compu tations high perfor mance implem entatio nin terms of faster compu tation

speed
and smalle r storage size is require d. In recent years, paralle
l compu ters have change d techniq ues
to solve proble ms in various kinds of fields. Domai n-deco mposit
ion is an efficient techni que to solve
partial differe ntial equati ons on paralle l compu ters of distrib uted
a compu tationa l domai n into a union of subdom ains, and eachmemor y type. The techniq ue divides

proces sor of a paralle l compu ter take

charge of a differe nt subdom ain.
In the domain -decom positio n situati on, we can consid er two types
conven tional progra ms of the finite elemen t metho d for non-paof metho ds. The first type is thatrallel compu ters are tuned to paralle lcompu ters by distrib uting data on a subdom ain to a proces sor.

This is called data-p arallel ization(based on the domain -decom positio n). The numer ical schemes execut ed by the tuned progra ms areequiva lent to the origina l scheme s [5, 3, 36]. The other type impose
on interfa ces of subdom ains, and solves proble ms iterativ ely

s an auxilia ry bound ary condit ion

in each sub domai n until the auxilia ry

bound ary condit ion converges to attain a slnoot h numer ical solution [40,47 ,21,51 ]. This is called thedomai n decom positio n metho d or the domai n decom positio n algorit

hm. It is possib le to use differe nt

mathe matica l models or differe nt numer ical metho ds in each subdom
metho d has the potent ial to decrea se commu nicatio n betwee n procescompu tations in each subdom ain is higher than the data paralle

ain. The domai n decom positio nsors becaus e the indepe ndency oflizatio n.

Bernar di et al. [7] propos ed a mortar elemen t metho d, that is akind of the domai n decom positio nmetho d. It is a Lagran ge multip lier metho d with a constr aintof contin uity
When the system is solved with the conjug ate gradie nt(CG )betwee n subdom ains.

metho d, the proced ure can be highly

paralle lized . A flow solver in the stream functio n - vortici ty
form by a mortar
propos ed[l, 2], but the way of the stream functio n - vortici ty

elemen t metho d isform made it difficult to treat natura lbound ary condit ions.
In this thesis we apply a domai n decom positio n algorit hm tothe upwin d finite elemen t scheme

based on the choice of up- and downw ind points. At first, we
introdu ce the idea of two data types,'pre-as sembly type' and 'restric tion type' for data on nodal points.

They simplif y the implem entatio n

and reduce the amoun t of commu nicatio n betwee n proces sors. Second
metho d to the Stokes equatio ns in the primiti ve variabl e form,

ary, we apply the mortar elemen twhich mainta ins the natura l treatm entof stress bound ary conditi ons. Its implem entatio n can be describ

ed concise ly by the idea of the twodata types. Finally , a domai n decom po ition algorit hm is applied

to the finite elemen t scheme basedon the choice of up- and downw ind points .

Orga nizati o n of t h e thesis
In Chapt r 1 the theory of the variati onal princip le and the formul
elemen t metho d for the incomp res ible 1 - avier-S tokes equatio

ation of the standa rd mixed finite

ns are review ed accord ing to [54, 59].

The former a sure the olvabil ity of the schem e discuss ed in
Section 2.4 and also plays

a key role in

the analys is in the domai n decom positio n m.etho d. The proced
ure of con istent discret ized pre sure3

Poisson equation by Gresho et al. [22] is introduced and we show a remark in order to apply iso-P2
PI/PI elements in boundary layers.
In Chapter 2, the upwind finite element approximation based on the choice of two upwind and two
downwind points are presented. We present efficient search algorithms of heuristic type and of blind
type for upwind-point search problems. The relation between one of the upwind-point search problems
and a point search problem in a subdivision of a lower-dimensional space is discussed. Afterwards weformulate the finite element scheme for the Navier-Stokes equations using the upwind approximation.
The choice of two parameters in the scherne are investigated numerically and theoretically. Somenumerical results are shown to demonstrate that the solver can produce numerical results possessing
good Reynolds-number dependency and capturing three-dimensionality of flows well.

In Chapter 3 we discuss the parallelization of the finite element scheme based on the choice of upand downwind points. At first we present two data types for data on nodal points. Next, we review
the mortar element method and apply it to the Stokes problem of primitive variable form. At last, weapply a domain decomposition algorithm to the finite element scheme based on the choice of up- anddownwind points.

Chapter 1

Mixed Finite EI(~ment Method for the
Navier-Stokes Equations
1.1

Introduction

In this chapter, we review the standard mixed finite element method for the incompressible NavierStokes equations and the procedure of consistent discretized pressure Poisson equation.
Let D be a bounded domain in Rd(d = 2,3) and T be a positive constant. We denote the spacecoordinates by x = (Xl"'" Xd) and the time by t. We consider motion of incompressible viscous fluidin D that is described by the incompressible Navier-Stokes equations:

au/ at + (u . grad)u + grad p(1/ Re)b..u + f

°

div u

in D x (0, T),

in D x (O,T),(1.1 )(1.2)

where the velocity, u : D x (0, T) -+ R d , and the pressure, p : D x (0, T) -+ R, are unknown functions,the body force, f : D x (0, T) -+ R d , is a known function, Re is the Reynolds number and(a=1' ... ' a=J T
div(a=1 ' .. , a=J 'ab.. = div . grad = L a

grad

d

2

i=l

2'
xi

(1.3)
(1.4)

(1.5)

Let f D( =I: 0) and f IV be two parts of the boundary aD,

aD = fD Uf,(1.6)

We consider the velocity boundary condition,

U=gD

onfDx(O,T),5(1. 7)

and the stress boundary condition,(Tn

on rN x (0, T),

= gN

(1.8)

where (T is the stress tensor

2-pO'ij + Re Dij( u),~ (~aui + aU j ) ,
2·

aXj

aXi

i, j = 1, ... ,d,i, j = 1, ... ,d,(1.9)(1.10)

n is the unit outward normal to r N, gD : r D x (0, T) - t Rd is a given velocity and gN : r N x (0, T) - t Rd

is a given surface force.
The initial condition isin n

u = UIC

at t = 0,(1.11)

where U/C : n - t Rd is a given initial velocity.

R e mark 1.1 The Navier-Stokes equations have two factors that make numerical computations moredifficult . They are the convection term (the second term in the left hand side of (1.1)) and theincompressibility condition (1.2). The equations in which the convection term is omitted:

au/ at + grad p
div u(1/ Re)~u + f

in n x (0, T),° n (O,T),
inx(1.12)(1.13)are called the Stokes equations.

1.2

Variational principle and the finite element method

The fin ite element method is based on the variational principle. In this section we briefly review it
according to [54] . P roofs of the theorems in this section are shown in [10, 20].

Variational principle of minimization type
Let .X and Y be Hilbert spaces . A mapping a : ~~ X Y - t R is called to be a continuous bilinear formon ~~ x Y if it satisfiesi) VUi E ~Y,Vv E Y,VCi E R(i = 12),a(clul +C2U2 v) = cla(ul,v) +C2 a (U2,V),

ii ) Vu E .X, VVi E Y VCi E R (i = 1,2),a(u, CIVI + C2V2) = cla(u, vd + C2a(U, V2),iii ) T here exist a constant Co such thatVu E X, Vv E Yla(u v)1 ::; collullxllvlly.
6(1.14)

If Y = X, it is called to be a continuous bilinear form on X. Let V be a closed subspace of X. A
continuous bilinear form a on X is called to be coercive on V if there exists a constant a > 0 suchthata(v ,v) 2:: allvll~.

"Iv E V,

(1.15)A continuous bilinear form a on X is called to be symmetric if

VU,VEX ,a(u,v)=a(v,u).

(1.16)A mapping 1 : X --+ R is called to be a continuous linear functional on X if it satisfiesii) There exists a constant Co such that

I(1, v) I ::; Co II v II x .

vv EX,( 1.17)The set of all continuous linear functionals on X is denoted by X'.
Now, let X be a Hilbert space, V a closed subspace a a continuous bilinear form on X which is
coercive on V. For 1 E X', we consider the following problem:
Problem 1.1 Find u E V such thata(u, v) = (1, v).

Vv E lV,(1.18)The existence and uniqueness of Problem 1.1 is shown as follows.
Theorem 1.1 Problem 1.1 has a unique solution. The mapping from 1 to u is an isomorphism fromV' onto V.
This result is known as Lax-:vIilgram s theorem when V = X. When a is symmetric, we also consider
Problem 1.2 (Variational problem of nninimization type) Find u E V such that

J [u] = min J[vJ,
vEV(1.19)

where1

J[v] == -a (v, v) - (f, v) .
2

(1.20)

is a functional on X. We already know that Problem 1.1 and 1.2 make no difference:
Theorem 1.2 (Variational principle of lminimization type) Problem 1.2 has a unique solution.

This solution coincides with the solution of Problem 1.1.
7

Ritz-Galerkin method
Let V be a subspace of a Hilbert space X. Choose a set of func t ions, {¢d i= 1,... ,N 1 eX , such that{¢di is linearly independent. Let X h be a linear span of {¢d i=1 ,.. .,N 1 • Let
N

Vh = {v h ;L Ci ¢i, Ci E R}

Vh =(1.21)

i=lbe a N(~ N"d dimensional subspace of V.
Let us consider a discrete version of Problem 1.1:

Problem 1.3 Find Uh E Vh such that(1.22)
Since Vh is a subspace of V, a is coercive on Vh if so is on V. Thus , these exists a unique solution of
Problem 1.3. ivloreover, if a is symmetric , the solution minimizes the functional J on Vh .
The solution
N

LCj¢j

Uh =

(1.23 )

j=l

can be computed by solving the system of linear equations ,N

L aijCj=b i ,i=l, . . . , N ,(1. 24)

j=l

where

(1.25)(1.26)
This method is called Ritz-Galerkin metho d .

Variat ional princ iple of saddle-point typ e
In the case where the solution subspace V C X is described by a constraint it may be difficult to apply
finite element approximations to the variational principle of minimization type. The incompressible
condition , that appears in the (l\"avier- )Stokes problems is an example of such a constraint. Therefore,variational principle of saddle-point type is used.
Let ~y and ilvI be Hilbert spaces, let a : .X X X--+ Randb : X x i.'vI --+ R be continuous bilinear

forms and let f : .X --+ Rand 9 : ..'vI --+ R be continuous linear functionals. Consider the problem:

8

Problem 1.4 Find a pair (u, p) E X x NI such that

Vv EX,a(u, v) + b(v,p) = (j, v),(1.27)Vq E NI,b( u, q) = (g, q).

(1.28)

By setting

V(g)

v

{v E X; Vq E iYI, b(v,q) = (g,q)},(1.29)

V(O),(1.30)

we have
Theorem 1.3 If a is coercive on V and if b satisfies the inf-sup condition, that is, there exists a

positive constant /3 > 0 such thatb( v, q)_> /3,

inf supIIvllxllqllM

qEM vEXthen Problem 1·4 has an unique solution.

(1.31)The mapping from data (f, g) E X' X lvI' to solution(u, p) E X x M is an isomorphism.
In the case where a is symmetric, by Theorem 1.4 stated below, Problem 1.4 corresponds toProblem 1.5 (Minimization problem under a constraint) Find u E V(g) such that

J[u] =min J[vJ,vEV(g)(1.32)

where J is the functional on X defined by (1.20).
Theorem 1.4 The solution u of Problem 1.4 is also a solution of Problem 1.5.

In order to avoid the constraint which appears in the definition of V(g), we introduce anotherfunctional L : .X x i'v1 --* R ,L(v, q) = J [v] + b(v, q) - (g, q),(1.33)and consider
Problem 1.6 (Variational problem of saddle-point type) Find a pair (u, p) E X x lvI such that

V(v, q) E X x i'v1,

L(u, q) ::; L(u,p) ::; L(v p).

(1.34)For the solution of Problem 1.6, it follows
Theorem 1.5 (Variational principle of saddle-point type) If a is semi positive definite on X,

that is,VvEX,

a(v,v) 20,(1.35)

then Problem 1.6 has a unique solution. This solution coincides with the solution of Problem 1.4.
9Mixed finite element approximation
Let X h C X and NIh C M be finite dimensional subspaces. Consider a discrete version of Problem
1.4,

Problem 1. 7 Find a pair (Uh' Ph) E X h X Mh such that
(1.36)(1.37)

Proposition 1.1 There exists a unique solution of Problem 1.7 if(i) The bilinear form a is coercive on(1.38)(ii) The bilinear form b satisfies the inf-sup condition with respect to X h x M h , that is, there exists a
constant f3h > 0 such that
(1.39)In the case where a is symmetric and senli-positive definite on X h , Problem 1.7 can be transformed
to a minimization problem under a constraint or a variational problem of saddle-point type. It is tobe noted that Vh is not a subspace of V in general.
In order to estimate the difference between the solutions of Problem 1.4 and the solutions ofProblem 1.7 the next theorem can be used .

T he o rem 1. 6 If the constant O'.h with respect to the coerciveness on Vh of a and the constant f3h withrespect to the inf-sup condition of b are independent of the discretization parameter h,(1.40)

holds, where c > 0 is independent of h.

1.3Standard mixed finite element method for the Navier-Stokes
equations with iso-P2 PI/PI elements

In this section we review the tandard mixed finite element method for the 1 - avier-Stokes equations
according to [59].

Variational formulation
We transform equations (1.1),(1.2),(1.7), (1.8) and (1.11) into a variational form. We introduce function spaces,

V(gD)
Q

{v E (H1(D))d;gD on f D },V =(1.41)

L2(D) ,

(1.42)

for the velocity and the pressure, respectively, and let V = V(O).

Remark 1.2 In the case of f N = 0, if a pair (u,p) is a solution of (1.1), so is (u,p + c), where c is aconstant. Thus we employ(1.43)
rather than L 2 (D).

Let v be a function in V. We integrate the product of (1.1) and v over D, apply Gauss-Green's
theorem and substitute the boundary condition (1.8). We also integrate the product of (1.2) andq E Q over D. Then the following variational form is derived:
Problem 1.8 Find (u,p) : (0, T) -+ V(gD) x Q such thatu(O)\Iv E V,(~~ (t), v) + al (u(t), u(t), v) + ao(u(t), v) + b(v,p(t))
Vq E Q,b(u(t),q)UIC,(1.44)(j(t), v),(1.45 )

O.

(1.46)

Notations that appear in the above problern are as follows.

(u,v)

In

U· vdx.

(1.4 7)i(1.48)dal(w,u,V)

I:

i,j=l

ao( u, v)

D

2 d
If
I:

aUi
wjavidx,
Xje i ,j=lb( v, q)

(j,v)

-In

i

i

Dij(u) Dij(v)dx,

q divv dx,

j· vdx+D(1.49)

DifN(1.50)

gN·vds.

(1.51 )The bilinear form (1.49) is employed for the viscosity term so that we can naturally treat the stressboundary conditions (1.8) .
11

Galerkin approximation
In the following of this section we assume gD =°for the simplicity so that V(gD) = V.

We approximate V and Q by finite dirnensional subspaces, Vh and Qh. Practical choice of thesesubspaces will be discussed later. Let(1.52)(1.53)

cPi : n -+ R d , i = 1, ... ,Nu , be base functions of Vh, and 'l/Ji : n -+ R , i = 1, ... , Np be base functionsof Qh. We set(1.54)
j=l

Np

Ph(t)

=

LPj(t)'l/Jj.

(1.55)

j=l

We employ the Galerkin method, that is, the same base functions are used for the discretization of v
and q in (1.45),(1.46). We approximate Problem 1.8 by a discrete problem:

Problem 1.9 Find (Uh' Ph) : (0, T) -+ Vh X Qh such thatddt (Uh(t), cPi) + a1(uh(t), Uh(t), cPi)+ao( Uh (t), cPi) + b( cPi, Ph (t))(1.56)b(V'h(t), 'l/Ji)

i = 1, ... , N p ,0,(Uh (0) , cPi)

(1.57)(1.58)

This is a system of ordinary differential-algebraic equations to find Uj, Pj : [0, T] -+ R , that is called asemi-discrete approximation.

Discretization in time
Let T be a time increment. We discretize Problem 1.9 in time. Several methods have been proposed
for the discretization, e.g. an implicit scherne[12]' a semi-explicit scheme[ll] and a three-step explicit
cheme[29]. vVe here employ a semi-explicit method[22].
Approximating the time derivative in (1.56) by a forward difference,d-Uhn('T ) f'.J

dt

U n+l - un

h

h

T

(1.59)

and evaluating other terms in (1.56) at time nT, (1.57) at time (n +- l)T, we obtain the following
system,Un+ 1 _ un)
( h T h, ¢i +- al(uh, Uh' ¢i)+- ao (Uh' ¢i) +- b(¢i , Ph)(1.60)
0,i = 1, ... , N p ,(1.61)(1.62)

where uh' Ph and jn (n = 0, ... , [TIT] == NT) denotes uh(nT), Ph(nT) and j(nT), respectively.
Setting vectors un, pn, pn as,
n
n ]T
['U1
,··· ,uNu'(1.63)]T ,(1.64)

r, n

~Pl'"

n

. ,PNp(1.65)
and a vector UIC as,(1.66)
the system is rewritten by

u n+ 1 _ un
M
T+- Al(un )U n +- AoUn +- BT pnpn,0,n=0, ... ,NT -1,
n = 0, ... , NT - 1,

(1.67)(1.68)(1.69)

where lvI Ai (un) and Ao are Nu x Nu-matrix who e components are(1. 70)[.4. 1 (Un)]ij(1. 71)[Ao]ij(1. 72)

They are called the mass matrix, the convection matrix and the viscosity matrix, respectively. The
components of the lVp x Nu-matrix Bare
(1.73)

It is called the divergence matrix.

13

Iso-P2 PI/PI finite elements
We divide 0 into a union of d-simplexes (triangles when d = 2 and tetrahedra when d = 3):(1.74)
An ej is called an ' element '. We suppose that ej is a open set. If i -I- j , ei n ej has to be one of anempty set , a common vertex, a whole edge or a whole face (a whole face is only in the case d = 3).
We call (1.74) as a subdivision of the domain and set(1.75)
where h represents the maximum of the diameters of the elements. We assume the triangulation isregular [10], i.e. there exists a constant Bo > 0 such that the angle of two arbitrary faces (two edges
in the two-dimensional case) of an arbitrary element is larger than or equal to Bo.
The space Qh is made by the PI (continuous piecewise linear) element:(1.76)
where pI is a space of linear functions.
Furthermore , we decompose each d-sim.plex into 2d sub-d-simplexes.

There are three ways todecompose a tetrahedron into 8 sub-tetrahedra. The choice of the way is discussed in Section 1.5.
The space Vh is made by the PI element on the sub-d-simplexes (see Figure 1.1):(1. 77)The bilinear form b : Vh x Qh -+ R defined by (1.50) satisfies the inf-sup condition[6]'(1. 78)where (3 > 0 is a constant independent of h.

Remark 1.3 Let us consider the standard mixed finite element method for the steady Stokes problem : Find (Uh Ph) E Vh x Qh such that ,(1. 79)o.

(1.80)

By virtue of Proposition 1.1 and Theorem 1.6, there exist a unique solution and the error estimation(1.40) holds .
14

p

u

Figure 1.1: Iso-P2 PI/PI elements, d = 3
If u E (H 2 (D))d and p E Hl(D), interpolation errors are estimated byinf

Ilu-Vhllv(1.81)

inflip - qhllQ <(1.82)

vhEVh
qhEQh

Substituting them into (1.40), we have the following error estimations.
(1.83)

1.4

Consistent discretized pressure Poisson equation

In order to solve the system (1.67),(1.68), a procedure of the consistent discretized pressure Poisson
equation[22] is efficient to reduce the computation cost.
In the scheme, the lumping technique is applied to the time difference term in (1.67), that is, themass matrix !vI is replaced by the lumped [nass matrix, !VI,
(1.84)
where(u, V)h =L u(P)v(P)measDp(1.85)

pis corresponding to (1.47). In (1.85) P rune; over all the nodal points of Vh , and D p is the barycentric domain with respect to P. Figure 1.2 illustrates the barycentric domain of a nodal point in atriangulation.

15

Figure 1.2: Barycentric domain
Introducing an intermediate variable, [;, we can separate (1.67) into
pn ,(1.86)o.

(1.87)

Operating BJVI- 1 to (1.87), and substituting (1.68), we derive(BM'-lBT)pn = ~B[; .
T(1.88)

Thus we can solve pn separately by (1.88) and aftenvards compute U n+ 1 by (1.87):(1.89)
We can consider [; as the velocity which is predicted from un by effects of convection and viscosity,n1and U + as the velocity which is corrected as to atisfy the incompressibility. The size of the linear
equation system to be solved is equal to iVp ., which i a relatively small number compared with Nu in
the iso-P2 PI/PI elements.
Proposition 1. 2 Since b satisfies the inf-s ~up condition. equation (l.88) is uniquely solvable(S6j.
Remark 1.4

NIas lumping is required to maintain par eness of the matrix (B~\lI-l B T ).

Remark 1.5 Equation (1.88) is corresponding to a pressure Poisson equation. Let us assume thatfunctions in (1.1)-(1.2) are sufficiently smooth. If we operate (-div) to (1.1) and substitute (1.2) , we
get a pressure Poisson equation,- ~p =: div( u . gradu) + div f.

(1. 90)This operation corresponds to the above derivation of (1.88). It is possible to apply the finite element
method directly to (1.90) but some artificial boundary conditions for the pressure are required in suchmethods.

1.5

Decomposition into sub-tetrahedra

In the three-dimensional iso-P2 PI/PI elements, each tetrahedron is decomposed into 8 sub-tetrahedra.
Then there exist three ways of the decomposition (see Figures 1.3 and 1.4). Obtained element matrices
are different depending on the three ways. The way that produces the best subdivision in the following
sense is recommended. Four sub-tetrahedra containing the vertices, e~ , ... , e~, are identical for thethree ways (see Figures 1.3(b) and 1.4(b)), while the other four sub-tetrahedra, e~, ... ,e~, are different
in general (see Figures 1.3( d-f) and 1.4( d-f)). We compute

q=

max

~la?C (rJek grad/\ · gradAjdX)(1.91)

k=5,6,7,8 t"tJ

for each way, where Ai( i = 1, ... ,4) are the barycentric coordinate on each sub-tetrahedron. A positiveq implies the existence of some obtuse subtetrahedron, which is considered to be a bad subdivision.
Thus we choose the way where the minimurn q is attained.
Figure 1.3 shows three ways of the decomposition of a tetrahedron whose sides PQ, QR and RS
of the original tetrahedron are orthogonal to each other. The way shown in Figure 1.3( d) is chosenbecause the way yields q = 0 and the other two ways make q > O. Figure 1.4 shows another example,
where the lengths of the sides PQ , PR and PS are 2a, 2b and 2c (a < b < c), re pectively, andthose sides are orthogonal to each other. Although the original tetrahedron is of weakly acute type,every group of sub-tetrahedra includes some obtuse sub-tetrahedra and q = ab/(6c), ac/(6b), bc/(6a),respectively. The decomposition of Figure 1.4(d) is cho en, which yields q = ab/(6c). When the
aspect ratio of the tetrahedral element i large, that is a « b or b «c such a situation often occurs
in boundary layers, the value q becomes pretty large in the ca e of Figures 1.4( e) and 1.4( f).

17

s

( c )

( b )

( a )

:I

~
.

!

( d)

.

.

.......... .

( e )

( f )

Figure 1.3: Example of a decomposition into 8 sub-tetrahedra

18b

R ~-----+-<
b

( a )

( b )

( c )

( d)

( e )

( f )

Figure 1.4: Example of a decomposition into 8 sub-tetrahedra

19

1.6

Conclusion

In this chapter, we have reviewed the standard mixed finite element methods for the incompressible
Navier-Stokes equations and the procedure of the consistent discretized pressure Poisson equation.
The fact that the iso-P2 PI/PI elements satisfies the inf-sup condition assures the unique solvability of the consistent discretized pressure Poisson equation and gives an error estimation of orderh for a class of flow problems. Thus one of the difficulty to solve the Navier-Stokes equations, theincompressibility condition, is overcomed by the iso-P2 PI/PI elements.
A remark for the suitable way to decoDlpose a tetrahedron into 8 sub-tetrahedra has been given
in Section 1.5. That is useful to use the three-dimensional iso-P2 PI/PI elements in the boundary
layers.

20

Chapter 2

Upwind Finite E~lement Scheme Basedon the Choice of' Up- and Down-wind
Points2.1

Introduction

At first we review necessity of upwind approximations. For detail and more examples, we refer to [59].
The Reynolds number is defined by

Re = UL,v(2.1)

where U is the representive velocity, L is the representive length and D is the kinematic viscosity. If
the Reyno lds number is large, the convection term, (u . grad)u, is dominant rather than the viscosity
term. Then it is difficult to obtain numerical solution with the Galerkin method because it is anapproximation of central type . In order to extract the essential aspect, let us consider an initial value
problem,a¢,at + (U· grad)¢(2.2)¢(x,O)(2.3)

where u E Rd is a constant vector. The solution of this problem is¢(x, t) = ¢IC(X - ut).

(2.4)

By this solution we can see that the convection is the phenomena of the temparature (or a quantity etc.)
on the flow. T hus upwind information is important to decide the phenomena. The Galerkin method
is the approximation that use information of every neighbor with respect to a point. YIethods of thistype is known to be unstable for (2.2) (2.3)[48] . In order to obtain stable numerical solutions to highReynolds-number problems we need to employ approximations that use more upwind information.

21

In the finite difference field , upwind approximat ions are easily attained by one-sided differences.
In the finite element field , since the Galerkin method is of central type , much effort has been paid:· Heinrich-Huyakorn-Zienkiewicz-Mitchell( [25],1977)· Tabata(Approximation based on the choice of an upwind element)([52]'1977)· Kanayama([32]'1980)· Baba-Tabata([4]'1981)· Brooks and Hughes(Stream upwindjPetrov-Galerkin scheme) ([9 , 26]'1982)· Ikeda([27],1983)·Ohmori-Ushijima([44]'1984)· Kondo-Tosaka-)Jishimura([37]'1990)· Tanahashi-Okanaga-Saito([63]'1990)· Kakuta-Tosaka([31]'1992)· Jiang-Lin-Povinelli([30],1992)· Franca-Frey-Hughes(The stabilized methods)([13, 14]1992)
In the finite element approximation based on the choice of an upwind element[52]' the simplicial
piecewise linear element is used. Let £ (0) be the set of elements that have 0 as a vertex,£( 0) = {e E 7h; 0 E e}.

(2.5)

The value of the convection term (u . grad¢h) at a nodal point 0 is approximated by(u · grad¢h)uo(O) = u(O) . grad¢h(eu),(2.6)

where eu, the upwind element of 0, is defined as an element which locates the upwind-side of 0 withrespect to u(O) (see Figure 2.1):

eu E £(0),{ 3E>0,0-EU(0) Ee u,(2.7)

Since grad¢h( eu) is constant, (2.6) is equivalent to(u· grad¢h)uo(O) = lu(O)1¢h(O) - ¢h(U)length OU '(2.8)

where U is the point that attains the maxinlal E in (2.7). Bilinear form al is replaced bya?(u, ¢h, h) = ~(u. grad¢h)uo(P)?J;h(P)measDp,

(2.9)

p

where P runs over all nodal points . The truncation error of this scheme is evaluated by(u· grad¢)uo(O) = (u· grad¢)(O) + chlu(0)la2¢ja~2(0) + 0(h 2),22(2.10)u(O)uFigure 2.1: Finite element approximation based on the choice of upwind element
for every sufficiently smooth function ¢, where c is a constant independent of h and ~ is a coordinate
in the flow-direction defined in the next section.
In order to capture Reynolds-number dependency of flows well, higher order approximations are
required . Bristeau et al. extended (2.8) to an approximation by choosing two upwind points, WandU (the choice of W is different from one defined later). The truncation error becomes(2.11)In this chapter we present a finite element approximation based on choosing not only upwind points
but also choosing downwind points.

2.2

Upwind finite element approximation based on the choice ofup- and down-wind points

We explain t he upwind approximation to the convection term u . grad¢, where u : n -+ Rd is a given
velocity and ¢ : n -+ R is a scalar function .
Let 0 be a nodal point. When u(O) = 0, we set (u· grad¢h(O) = o. W hen u(O) ¥ 0 we consider
a li ne l passing through 0 in the u( O)-direction and introduce a coordinate ~ on l with origin 0 andunit vector u( 0) / lu( 0) I. Then the convection term can be written as(u · grad¢)(O) = lu(O)la¢/a~(O).

(2.12)In orde r to app roximate a¢/a~(O) we choose four other points U, W, D and B on I such that U
and ltV are upwind point of 0 with respect to the flow u( 0) and D and B are downwind points (seeFigure 2.2 and 2.3). T he method of practical choice of these points is presented in the follow ing
section . Setting
h = length OU

23(2 .13)

Figure 2.2: Velocity u(o) and points W, U, 0, D, B d = 2

Figure 2.3: Velocity u( 0) and points W, U, 0, D, B, d = 3
we denote the ~-coordinate of P E {W, U, 0, D, B} by ~P = ~~h (~o = 0, ~~ = -1). We approximatea¢/a~(O) by the values ¢ at the five points W, U, 0, D and B as follows:I:(u· grad¢)h(O) = lu(O)1 {1P¢(P) } / h,

(2.14)

PE{W,U,O,D,B}

where

1p

II (~~ - ~Q)'
{ II (~o - ~Q) + a}/Q¢P- L 1Q,P =f: 0,(2 .15 )Q¢~O

10(2.16)Q¢O

and a is a non-negative parameter.
vVe approximate the trilinear form,(2.17)
24

by(2.18)
pwhere P runs over all nodal points.
The following property had been shown in [56].

Lemma 2.1 Suppose that there exist positive constants Cl and C2 independent of h such that

(2.19)
Then, for every sufficiently smooth function ¢, it follows that

Proof: We set 'YP = 'Y~ + 'Y~, where{ II °(~o - ~Q) } / II (~~ - ~Q)' P
Q=/:-P,

1'Yo=1= 0,(2.21)

Q=/:-P- L 'Yb,(2.22)

Q=/:-Oa/ II (~~ -- ~Q).

(2.23)

Q=/:-P

Then we have'Y~¢(P)/h(2.24)

PE{W,U,O,D,B}

L'Y~¢(P)/h

PE{W,U,O,D ,B}

where P4(~; ¢(P), P = W, U, 0, D , B) is the fourth order Lagrange interpolation and¢[W, U, 0, D, B] is the fourth-order divided difference. Hence we obtain (2.20).

R e mark 2.1

The cond it ion (2.19) is required in order to avoid a cancellation of significant digits.

In fact if the distance between any two of the five points tends to 0, the absolute values of several
coefficients 'Y tends to +00.

Remark 2.2(i) When a = 0, (2.14) is the fourth order central approximation.

(ii) Equat ion (2. 14) becomes a third-order upwind approximation using W,U,O and D whena = --II (~o - ~Q).
Q=/:-BO

25(2 .26)

We have reported this case briefly in [55].
(iii) Equation (2.14) becomes a third-order upwind approximation using W ,U ,O and B whenIT (~o - ~Q)'a = -(2.27)

Qf:-D ,O(iv) Suppose the case where W, U, 0, D and B are located uniformly, that is, the lengths of BD , DO ,OU and UvV are equal to h. Then (2.14) becomes the third order Kawamura et al.-type approximation[35] when a = 6. ~ote that a = 2 in (2.26) and a = 4 in (2 .27) when W, U, 0, D and B are locateduniformly.
The parameter a stabilizes the scheme and we call it the stabilizing parameter. We choose abetween 0 and 6 from the above remark . V\Te will investigate the choice of a and T in Section 2.5.
Since we use the piecewise linear element , the values of Uh at upwind or downwind points are
obtained by the linear interpolation of nodal values. Therefore (2.14) is not of third-order accuracy ingeneral, but it is when these points are located suitably, e.g. when they coincide with nodal points.
In the neighborhood of the boundary where the points W or B cannot be chosen , we employ the
first-order upwind approximation (2.8) or the third-order scheme (2.26) instead of (2.14) .

Choice of up- and down-wind points
vVe choose the upwind and downwind points in the following manner.
We denote by F the collection of all the (closed) faces of the elements ( only in the three-dimensionalcase ), all the (closed) edges of the elements and all the vertices. We denote by C(F) all elements in

F which have no intersection with F E F:

C(F) = {A E F;An F = 0} .

(2.28)The upwind point U(~u) on l is defined by~u = max{~;~ < 0,P(O E F E C(O)}.

(2.29)

Let dim(F) be the dimension of F . vVe de:fine Fu as the element F in C(O) which includes U andattainsdim(Fu) = min{dim(F); U E F E C(O)}.

(2.30)

This means that F u = U when U coincides with a nodal point F u i an edge when U i just on theedge otherwise Fu is a face in the three-dimensional case. In the two-dimensional case, Fu = U when26

U coincides with a nodal point, otherwise Fu is an edge. The second upwind point W(~w) on l isdefined by~w = max{~;~ < ~u, P(O E F E C(Fu)}·(2.31)The downwind points D and B are defined as the upwind points U and W corresponding to -u( 0).
By virtue of the assumption that the finite element subdivision is regular, we have:
Proposition 2.1 The upwind and downwind points chosen as above satisfy the relation {2.19}.

2.3

Search algorithm for the upwind points

In order to evaluate (2. 14), we need to determine two upwind points and two downwind points definedby (2 .29) ,(2.31) with respect to all the nodal. points of the velocity. In the computation of unstationary
flow problems, the velocity field changes every time step so that we need to determine such points atevery time step. This determination needs to solve point search problems:
Problem 2.1 Search the upwind point U defined by {2.29} with respect to a nodal point 0 and thevelocity u( 0) on a finite element subdivision Th of a domain.

Problem 2.2 Search the upwind point W defined by {2.31} with respect to a nodal point 0 and thevelocity u( 0) on a finite element subdivision Th of a domain.

Fast search algorithms are required in order to apply the upwind finite element approximation withlow computer cost. In this section we discuss algorithms for the search problems.

2.3.1

Heuristic- and blind-search algorithms for U

Blind search algorithm

At first let us consider Problem 2.1. It is essentially same to choose an upwind element.
------t------t

Let e be an element in [( 0) and Pi , ... ) P d be vertices of e other than o. The vectors, 0 Pi) . .. 0 P d ,are linearly independent so that we can represent the vector u( 0) as the linear combination of them:du(O) =

LCiM.

(2 .32)i=i

If all the coefficients are not positive
Ci ::; 0i = 1, . . . , d(2.33)

the element e is located in the upwind side of O. Hence a search algorithm of blind type is written as
follows.
27

Algorithm 2.1 (Blind search algorithm)
For j = 1, ... ,Ne:Let Po, ... , P d be the vertices of ej.
For i = 0, ... , d:If ej is located in the upwind side of Pi with respect to U(Pi):
Record ej as the upwind element of Pi.

Heuristic search algorithm
A search algorithm of heuristic type is written as follows (see also Figure 2.4).

Algorithm 2.2 (Heuristic search algorithm)
1. Choose an element e E £(0).

2. Let PI, ... , Pd and 0 be the vertices of e.
----+

----+

3. Factorize u(O) by OPI , ... , OPd as (2.32).

4. For i = 1, ... , d:
If Ci > 0:
Replace e by the element in £(0) that shares all the vertices except Pi,
Go back to Step 2.
o. e is the upwind element of O.
In unstational flow computations, change of the velocity at a nodal point can be expected small.
Thus if we choose the upwind element at the last time step as e in Step 1, we can expect to reduce
iteration times of the loop Step 2 - Step 4.
Let (0, AI, . .. , Ad) be the barycentric coordinates of U in the upwind element:
----+

d
~'

---+

OU = L, AiOPi.

(2.34)

i=l----+

Substituting (2.32) into the definition of the ~-coordinate of U, OU = ~uu( 0) Ilu( 0) I, the relations,~u

Ai = Ci lu(O)I'i = 1. ... d,(2 .35)t(2.36)

holds . Summation of this equation yields

lu(O)1 =:
Ci.
~u
i=1

Values ~u and Ai, i = 1, ... d, are used to compute the quantities in (2.14),(2.15),(2.16). Thus all
the result of computations in the point-searching are also u ed to evaluate the convection term if noiteration occur in the heuristic search algorithm.
28

u(O)

u
Figure 2.4: Heuristic search of point U
Table 2.1: Comparison of the heuristic- and the blind-search algorithm on a scalar- and a vectorcomputer
Time
Step
1
2
3
4
5

Scalar
Heuristic
109.
146.
66.6
31.4
30.2

Blind
138.
138.
138.
138.
138.

VectorVelocity Field
Heuristic Blind
64.5
1.89
u 1 = (X2 - 0.51, -Xl + 0.501, -1.0)
2
92.0
1.90
u = (-X2 + 0.51, Xl - 0.501,1.0)
3
32.0
1.90
u = (-X2 + 0.05 , Xl - 0.041,1.0)
4
4.88
1.94
u = (-X2 + 0.01, Xl - 0.001,1.0)
uS = (-X2 + 0.01, Xl - 0.001,1.0)
3.68
1.93
(10- 3 sec.)

Comparison of two algorithms

The heuristic search algorithm is expected to attain less computation step while the blind search canbe written by a program of simple structure. Performance of both algorithms on a scalar computerand a vector computer was checked numerically in [16].
The benchmark problem is the upwind point search at 5 time step in a 10 x 10 x 10 x 6 tetrahedral
subdivision of a cube domain.

-irtual velocity fields , u I , ... , u 5 , are given at the time steps as shownin Table 2.1. Change of the velocity fields at the second. the third , the fourth and the fifth time step
are stated as u 2 = -u l , lu 3 -u2 1is large lu4 - u3 1is small and u 5 = u 4 , respectively. Scalar machineis HITAC YI-880 and the vector one is HITAC S-:3800. Their computation times are also listed inTable 2.1.

29

We can see that the computation time of the heuristic search algorithm much depends on thetime step. At the first time step, initial upwind elements in the algorithm are given in random sothat it needs relatively long computation time. At the later time steps, if change of the velocity fieldis smaller, computation time becomes smaller. In the practical flow computation, the change of thevelocity field is considered to be small so that performance at the 4th or the 5th time step may be reallyattained. Consequently, we can obtain good performance of scalar computers by the heuristic searchalgorithms, while the blind search algorithm constantly shows good performance of vector computers
by virtue of a simple program structure.

Relation with a point search problem in a field
Problem 2.1 is related to the following point search problem in the (d - 1 )-dimensional field.

Problem 2.3 Let Th be a (d - I)-simplex subdivision of R d - 1 , and P be a point in R d - 1 . Find a(d - 1) -simplex e E Th such that PEe.
This problem appears in some numerical methods for flow problems, e.g., computations using marker
particles, visualization by streaklines. For Problem 2.3 a search algorithm of heuristic type had been
proposed by Lawson[39]:Algorithm 2.3 (Lawson)
1. Choose a e(O) E Th . Set n = O.
) b e t h e vertIces.
f (n) .
2 . L et p.0(n) , ... , p(nd0 e3. Represent P by the barycentric coordinates of e(n) :d~

P = L )~~n) p/ n ) .

(2.37)

i=O

4. For i = 0, ... , d:

If A~n) < 0:
Set e(n+l) the d-simplex that shares all the vertices except Pi(n) with e(n),n:= n + 1,

Go back to Step 2.

5. e(n) is the element containing P.

Proposition 2.2 vVhen d = 2, execution of the Lawson's algorithm always terminates, if the numberof segments is finite.

Figure2.5: Example of an infinite loop in execution of the Lawson's algorithm, d = 3
Remark 2.3 When d = 3, there are situations where execution of the Lawson 's algorithm falls into

an infinite loop and never terminates (see Figure 2.5 ) .
The relation between Problem 2.1 and Problenn 2.3 is explained as follows (see also Figure 2.6).
The upwind face (edge when d = 2) is on a m-polygon that contains point°

as an inner point.

Consider a unit sphere(circle) whose center is 0 , and move all the vertices of the m-polygon onto thesphere( circle) wi thou t changing the angle from the point 0. Let the points of intersection between
the line l and the sphere(circle) be P and Q (~p > 0, ~Q < 0). Project the m points onto a plane
which is tangential to the sphere(circle) at point Q by a point source of light at P. These operation
gives a subdivision of a plane(line). Thus, the signs of Ci in (2.32) and the signs of A~n ) in (2.37) for apair of corresponding elements are same.
Consequently, for the upwind-point search problem, the following holds.
Proposition 2.3 When d = 2, ifTrt is regular and 0 is an inner nodal point, execution of the heuristic
search algorithm always terminates

From Remark 2.3 we can see that there also exist exceptional cases for the heuristic search algorithm
when d = 3. We have not yet experienced such a case in practical computations. The algorithm may
always terminate under some kind of regularity condition of a subdivision.

Figure 2.6: Relation between the search problem of U and a point search problem, d = 2

2.3.2

Heuristic search algorithm for W

We assume that the other upwind point W is contained in a face (edge in 2D) F E C(Fu). We callFu and F the upwind face and the second upwind face, respectively. The second upwind face can besearched by the following algorithm:

Algorithm 2.4 (Heuristic search algorithm for VV)
1. If the upwind face was not changed from the previous time step:

If the second upwind face at the previous tinae step, PI ... P d, is still located at the
upwind direction of 0 go to Step 9.
2. Denote Fu by PI . . . Pd· Set the logical variables h, ... , fd true, j = O.
3. Let PI ... PdQ be the neighbor element which is located at the opposite side of point O.
4. Do j := j + 1; until fj is false.
5. Decompose u( 0) such thatu(O) = CjoQ +L ciM.

(2.38)

ii=j

If one of Cl, ... ,Cd is positive j := j + 1, otherwise go to Step 7.
6. In the three-dimensional case, decompose u( 0) such thatu(O) = CjoQ -+L CiM.
if:-j

If one of CI, ... ,Cd is positive j := j + 1.
7. Replace Pj by Q . Set fj to false.
32(2.39)

8. If one of il,' .. , id is true, go to Step 3.

9. Pl··· Pd is the second upwind face.
When the algorithm terminates the values cl, ... , Cd are used to evaluate the upwind approximationin the similar manner as in(2.35),(2.36).

Remark 2.4 Algorithm 2.4 is constructed for almost all cases. It is possible theoretically to construct
an algorithm that covers all degenerated cases, for example, the case when Fu is just a nodal point. But
it is difficult to implement the ideal algorithm because of finite-precision arithmetic of real computers.

2.3.3

A topology oriented method in heuristic algorithms

In the heuristic algorithms in the previous sections, we need to express u( 0) by the linear combinationof vectors,du( 0) =

2= CiM.

(2.40 )

i==l

It is not good to solve this equation by Gauss's elimination method, since the method may produce
some wrong signs in a solution (Cl, ... , Cd) because of the finite-precision arithmetic. Thus we have arisk that an execution of a heuristic algorithm may fall into an infinite loop. The risk is overcomedby the following method, which essentially use the Cramel's formula.
The location of a point P is expressed with the barycentric coordinates of P, (Ao, .. . , Ad) (~t=o Ai =1), byd

DP= L:AiM.

(2.41)

i=l

We therefore recognize that ((1 - ~t=l Ci), cl, . . . , Cd) is the barycentric coordinates of the point 0 +u(O) .
Relative barycentric coordinates are efficient tools to compute that.

We explain in the three-dimensional case. Let (~, 7J, () be the Cartesian coordinates . Let the vertices of a tetrahedron be
Pi(~(i),7J(i),((i)),i = 0,1,2,3 where Po = O. The barycentric coordinates of P(~,7J,() are

Ai

a(i)

+ b~i) ~ + b~i)71 + b~i\,a(i)( _l)iVb( i)((_l)i 1
V

~(j)

~(k)

~(l)

7J(j)

7J

(k)

7J(l)

((j)

C(k)

((l)

1

1

7J(j)

7J(k)

1
7J( I)

((j)

((k)

((I)

33

i=0,1,2,3,(2.42)(2.43)(2.44)

( -1 )i+ l

b(i )
1]

V
(_I)i+ l

b(i)
(

V

1

1

1

((j )

( (k )( (l )

~(j )

~(k)

~(l)

1

1

1

~(j )

~(k)

~(l )

7](j )

7] (k)

7] ( l)(2.45)(2.46)

where (i,j , k, l) is a cyclic permutation of (0, 1, 2, 3) and V is the volume of the tetrahedron. Since thebarycentric coordinates of Po = 0 is (1 , 0,0 , 0), we have
(2.4 7)
Subtracting this equation from (2.42) , we obtain'·(P) = (b(i)b(i) b(i)) . O----tp ,~ , 1] , (. = 1, 2, 3 .

(2.48)(b(i)b(i) b(i)) . u (O) , 1 ,. = 1, 2, 3 .
~ ' 1] , ((2.49)/\t1,

Hence we can compute coefficients Cl, C2 , C3 by
Ct· =i
In the definitions of bi ), b~i) , b~i), i = 0, ... , d, t he coordinate of Pi does not appear so that thesigns of Ci is not affected by Pi. We can avoid infinite loops by this fact. We explain the reason in an
example. The correct signs are Cl > 0, C2 <°at elernent el and Cl , C2 <°at element e2 in the upwind

point search in Figure 2.4. If only the sign of C2 in e2 were computed wrong, the execution of thealgorithm would fall into an infinite loop . But because of above fact , Cl in el and C2 in e2 always have
opposite signs regardless of round-off error so that either el or e2 is decided as an upwind element andthe execution of the algorithm terminates. This method is a kind of topology-oriented method[50J.

R e mark 2 .5 The first de rivatives of the barycentr ic coordinates , bi ) , b~i), b~i), i = 0, ... , d , are worthistoring in the memory since they can be also used to compute the element stiffness matrices [16J.

2.4

Finite element scheme based on the choice of up- and downwind points for the Navier-Stokes equations

We consider a finite element scheme for the ~av ier-Stokes equations by combining• Standard mixed finite element formulation for the ~avier-Stokes equations with the iso-P2 PI/PI
elements (except the followings) ,• Lumping technique for the time-difference terrD ,34• Procedure of the consistent discretized pressure Poisson equation,• Upwind finite element approximation based on the choice of up- and downwind points for theconvection term, that is, replacement of al byda?(wh, Uh, Vh) = L

L(Wh . graduhih(P)vhi(P)measDp.

(2.50)P i=l

Therefore the finite element scheme solves
Problem 2.4 Find (uh,ph)n=O, ... ,NT E Vh x Qh such that

(2.51)

0,(2 .52)(2.53)We use an iterative method (ICCG, see [41]) as the matrix solver when d = 3, while we use a direct
method (the Cholesky decomposition for band matrices) when d = 2. The employment of the direct
method in the three-dimensional code is not realistic because it requires huge memory (or storage).

2.5

Stability regions with respect to a and T

2.5.1

Numerical experiments

There are two parameters a and T in the scheme. Since (2 .20) and in order to reduce a computation
time, it is better to choose a small a and a large T in the range of stable computation. In this section,
we find stability regions numerically in the a - T plane where we can compute problems stably. We
will show that for two problems. One is a three-dimensional test problem and the other is a flow
problem past a circular cylinder in the two-dimensional space.
Let n be a unit cube. We decompose n into the union of 6 x 6 x 6 cubes and each cube into 6
tetrahedra. Let[Xl(Xi-3X~)/3 ]UP =X2(X~ - 3x~)/3
X3(X~ - 3xi)/32
xl + X22 + X32 - 1 ,(2 .54)(2.55 )

be an exact solution and obtain the corresponding external force of the ~avier-Stokes equations. We
impose velocity boundary conditions gD = U on fD(= an) and set an initial condition u/c = O. For
35

0.175

0.150

0.125

0.100

0.075
0.050

~
Re=l000.025

o~--~~----------------------------a
1
2
3
4
5
6
Figure 2.7: Stability regions for Re = 100 , 400 and 1000 in a test problem
the increments ~a = 0.1 and ~T = 0.005 we compute the problem and check if the solutions converge
or diverge while time marching.
The regions below the lines in Figure 2.7 show the stability regions where the solution converged.
When Re = 100, there exists a stability region even ex = 0, and it i recognized that we can compute theproblem in this case without an upwind approximation. When Re = 400, the computation diverged
for every T if a < 0.9. We have to choose a ~ 0.9. 'When Re = 1000, we have to choose a ~ 1.4.
'rVe next treat a more practical problem, that is two-dimensional flows past a circular cylinder. Let

Lin Lout and ltV be constants and n c R2 be ,n = { x-Ixl > ~ x E (-Lin, Lout) X ( - W, W) } .
The problem is to find a pair (u.p) which satisfies (1.1) and (1.2) in(2.56)n x (0, T) under no external force.

'rVe set the initial condition u = O. The boundary conditions and the finite element subdivision of thedomain are shown in Figure 2. : where Lin = ltV = 7.5 and Lout = 22.5, and (Tl,"" Td) is the stressvector on the boundary:
d

Ti =

L
j=l

36

(5ijnj.

(2.57)

0

0

II

II

C'1

f1....

~

....

r-1

0

II

II

r-i

r-i

~

l-

Figure 2.8: Two-dimensional flO\v problem around a circular cylinder, boundary conditions( u = 0 on
the cylinder surface) and a divi ion of the domain
The total element number (of p) i --1320 and the minimum element size (of p) is hp

0.006. We

compute the problem at Re = 3200 and Re = 10000.
Figures 2.9 and 2.10 shows the result. where symbols 0 mean the computations were stable while
symbols x mean the computations diverged. In all stable computations, we obtained Karman vortex
sheddings. We can observe that \H' have to choose 0: ;::: 2.25 in the case of Re = 3200 and 0: ;::: 2.5 inthe case of R e = 10000.
The stability regions of the two problem have imilar shapes to each other. The figures of theregions show the follO\ving. (i) In orner to choose a larger T \ve need to choose a smaller 0:. (ii) But
the choice of too small 0: makes the computation u.n table regardless of T. (iii) The lower bound ofa increases as Reynolds numb r increa ~ . Thu we recognize nunlerically that a suitable choice of 0:
stabilizes the scheme .
Although we examined only in t\\·o sample problems. uch a tendency of stability is considered to
be valid for a wide claof flo\,- problem . In every computation in this thesis we choose a pair of 0:

and T which is con idered to be near t he top of the tability region (T is almost equal to the maximum
stable time step).

37

To. 010
x
x
x
x
x
x
x
x
x
x

O. 008x
O. 006

Xx x x x
x x x
x x x
x x
x x x
x x
x x x

x

x

x
x

0

0

O. 004

L
0

I
1

2

3

4

I

I

5

6

a

>

Figure 2.9: Stability regions, Re = 3200, in the flow problem around a circular cylinder

T

x

x

o. 010

x
x
x

x

x

x x

x
O. 008x x
x0

O. 006

x

x

0

1

2

3

0

O. 004

o

4

5

6

a

>

Figure 2.10: Stability regions, Re = 10000, in the flow problem around a circular cylinder

38

2.5.2A theoretical consideration

The above result (i) is suggested from the stability condition(2.58)of the forward-time and centered-space (FTCS) type finite difference scheme for the one-dimensional
linear scalar equation,(2.59)
where K, is a positive constant. From (2 .20) we see K, corresponds to (aluI/4!)h 3 , which implies(2.60)
This means the stability region is below a hyperbolic curve in the a - T plane, which is observed forlarger a in Figures 2.7, 2.9 and 2.10.
We refer to [43] for linear stability analysis of higher-order upwind finite difference approximationsof one-dimensional convection-diffusion equations. Stability analysis of such a type of the upwind
scheme for the Navier-Stokes equations is an open problem.

Derivation of condition (2.58) by the von Neumann stability analysis[49]Let consider (2.59) in R x (0, T) with a suitable initial condition at t = O. With a uniform gridsize h > 0 and a time increment T, FTCS type finite difference scheme is(2 .61)
where ¢i denotes the value ¢( ih, nT). Setting d = K,T I h4, we can rewrite the scheme as(2.62)A Fourier component of the numerical solution whose wave number is k (wavelength is 27r I k) is writtenby(2.63)
where (j;n is the amplitude function at the nth time step (also depends on k) and I = vCI. Define
the phase angle e = kh. Replacing ¢i in (2 .62) by(2.64)

39

we obtain(2.65)
Canceling the common term eIiO gives(pn {I - d(e- 2IO ._ 4e- 10 + 6_ 4e IO + e 2IO }(2.66)(pn[l - d{6 - 4(2 cos 8) + 2 cos 28}]

(2.67)¢;n{l - 4d(1 - cos 8)2}.

(2.68)

The amplification factor Go, defined by(2.69)
is the form of

Go = 1 - 4d(1 - cos 8)2.

(2.70)

If solutions are to remain bounded, we have

IGol:S 1(2.71)- 1 ~ 1 - 4d(1 - cos 8)2 ~ 1(2.72)

d ~ 1/8.

(2.73)

for all 8. Thus we havethat yields

2.6

Numerical results

2.6.1

Flows in a cubic cavity

We compute flows in a lid-driven cubic cavity as the first model problem.
Let D be a unit cube and consider the problem to find a pair (u,p) which satisfies the Navier-Stokes
equations subject to the boundary data on the lid,gD = (1,0 0) on r Dl == {x E 0S1'x3 = I}

(2.74)and the no-slip boundary condition on the other walls,
9D =

0 on

r DO := aD \ r Dl.

The external force is f = o.
40

(2.75)

This is a typical bench-mark problem for solvers of the ~avier-Stokes equations and several computation have been done. In [38, 19, 18, 34]' it has been reported that the flow becomes stationary at

Re = 10 3 and in [28] so is at Re = 2 x 10 3 . As for pairs of Taylor-G6rtler-like vortices, they are observed at Re = 2000 in [28] and at higher Reynolds numbers in [19, 18, 34, 45]. Recently Jiang[30] has
shown on the uniform 50 x 50 x 50 subdivision that the flow becomes unstationary even at Re = 400
and obtained pairs of Taylor-G6rtler-like vortices even at Re = 1000.
Here we compute the cases of Re = 100,400, 1000 and 2000. We analyze the problem in a halfdomain imposing the symmetric conditions
U2 =(2.76)

0,on the plane X2 = 0.5. The finite element subdivision of the domain is made by dividing the domain
into a union of 20 x 10 x 20(16 x 8 x 16 when Re = 100) non uniform bricks and subdividing each bricks
into six tetrahedra (see Figure 2.11). The minimum element size for the pressure is hp :::: 0.0172(0.0379
when Re = 100). We impose u = 0 as the initial condition in the case of Re S; 1000, while the steadystate solution of Re = 1000 is employed in the case of Re = 2000. We solve the Navier-Stokes equations
evolutionally. The parameters of the scheme are listed in Table 2.2.
In the cases of Re = 100,400 and 1000 the nurnerical solutions converged to steady solutions inthe sense satisfying the inequality

max{llun+l_unll,lIpn+l_pnll} < 10- 3 ,T max{lIunll, IIpnll}(2.77)

where II . II is the maximum norm. The times of convergence are also listed in Table 2.2.

Re = 2000, the values of the left-hand side of inequality (2.77) remained between 0.8

/"'.J

\Vhen

1.0 and did

not converge. This fact may suggest that the flow is unsteady at Re = 2000.
Figures 2.12, 2.13 and 2.14 show the velocity vectors of the steady flows obtained in the case of

Re = 100, 400 and 1000, respectively. Figure 2.15 shows the velocity vectors of the flow at t = 50 in
the case of Re = 2000. Those vectors are projected orthogonally to the three planes, X2 = 0.5, Xl = 0.5
and X3 = 0.5 and the length of the vectors are doubled in the two latter planes to observe flows clearly.
vVe can observe on the symmetric plane, X2 = 0.5, in Figures 2.12-2.15 that the center of the primaryvortex moves towards the center of the cube as the Reynolds number increases and secondary vorticeappear in the two lower corners. This behavior is similar to that in the two-dimensional lid-driven
square cavity. But there does not appear a econdary vortex near the upper left corner. At Re = 400,
we can observe a pair of secondary vortices near the bottom wall on the plane Xl = 0.5 and near the41

Table 2.2: Parameters of the scheme and convergence times

Re
100
400
1000
2000

a
1.5
1.5
2.0
2.0

T

0.003
0.004
D.OOS
0.010

t
5.61
15.94
27.2S

upstream wall on the plane X3 = 0.5. When Re =: 1000 another pair of secondary vortices appears
near the top wall on the plane Xl = 0.5. When Re = 2000 a pair of Taylor-Gortler-like vortices canbe observed around (O.OS, 0.5 ± 0.1,0.07) (see Figure 2.15).
Figure 2.16 depicts paths of marker particles starting from (0.5, 0.75, 0.95) in a steady flow at

Re = 100 ,400 and 1000. The path seems to remain on the surface of a torus when Re = 100. But
at higher Reynolds numbers the paths become distorted and much complicated. This behavior agrees
well with the results of [2S].
Figures 2.17 and 2.1S show the graphs of UI (0.5,0.5 .) and U3(-, 0.5, 0.5) in the case of Re = 400
and 1000 with other results obtained by(a) Fuseki[19]' finite difference method, 62 x 62 x 62 (R e = 400,1000),(b) Fukuyama[lS], finite volume method, 20 x 20 x 20 (Re = 400),(c) Kato et al.[34], finite element method, 20 x 20 x 20 (Re = 400),30 x 30 x 30(Re = 1000),(d) Ku et al.[3S], quasi spectral method, 24 x 12 x 24(Re = 400) 30 x 15 x 30(Re = 1000). Our
results are almost the same as or a little better than the results of similar subdivisions but they looka little more diffusive than the results of finer subdivisions.
From these observations, we can recognize that our scheme captured three-dimensional flow behavior well although we used relatively rough meshes.

2.6.2

Flows past a circular cylinder

We next compute flows past a circular cylinder. In the computation of the two-dimensional problem
explained in Section 2.5 we obtained Karman vortex sheddings. But physical experiments[23] showthat there appear three dimensional movements behind the cylinder over Re = 200 and this movements
grow as Reynolds number increases. It means that two-dimensional flows are unstable in the threedimensional space for high Reynolds numbers. vVe refer to [24 42, 61, 62] for numerical investigationsof such instability and comparisons between two- and three-dimensional results. We here compute theproblem in the case of Re = 1000 and observe the effect to three-dimensional movement by the non-42

Figure 2.11: Fini te element ubdi\'i ion of a half cube
uniformity of the subdivision and the transition from a two-dimensional flo\v to a three-dimensionalflow.
Let H > 0 be a constant. The comput ation domain 0 C R :3 i defined by

O=02 x (0.H ).

(2.78)

where O 2 is defined by the right hand sid e of (:2 ..56 ). Th e external force is J = O. The boundary
conditions corresponding to those of two-dimensional ca e a re shov;n in Figure 2.19. We assume thatthe phenomena is periodic of period H in the T:3- di r ection.

(2.79)

p(-,·. X:3+ H : ·)(2.80)As the initial datum we used a Karman vortrx -head ing rpsult of th corre. pondi ng two-dimensional
computation at a specific time. A small d i. t urbancp () {I - co.-( iTt)} / 2 \Va - added to the X3- com ponentof the velocity on the cylinder surface for 0 :s; t :s; :? to ' e .he 'fi nal ' ' tate of th solution more quickly.
vVe observed the solution ufficiently afterwards \\'hp n rh effect of he eli t urbance was considered todisappear. (vVe have checked that the re ult is <llmo:'t ind pendent of a (0 '" 0.5) later enough.)-1:3

(a) X2= 0.5

X3- ---_ ._ --.-,.

"'~'"'''''''''''''''--'--------~/////III'''''. . . ...............
...... '" " ' ' ' " , .... , ............... _---------',...". "'"" .. . .
•

•

•

..

. .
... ..
....

..

'

.

•

•

.

. .

. . ..
4'

..

..

..

..

..

..

. .

...

•

..

.

. . . . . . .. ..... . ... .. . .

. . .

" i 1 ~ ~~~~ ii:::::::::~~~~~~~ l\'::;~~~ZZZZZZ 7.

:d~~1.%~~~~~ 1. .

:Un~~~~~~~ ~ ~~ ~~~~~~~n~::~~~~~\~~~~~::~~~~~~~~~nl::\ \\\\\\\\\\\\1,",/1/////////11
01111/1// / / / / I I I \ \ \ \\\\\\\\\\1,111111 1 1/ //1 1111 1 \ \ \ \ \\\\\\\11,",11111111111• \ \ ' \ \\ \\\1111,. . • \ , \ \ \ \ 1/11,

I II I I (

I

I

I

,

" 11 I I I I

I

,

,

.

·

•

,

I

I

"It 1 1 \ I

,

•••

•

••

,

, I 1 I I I ,.

01

.Il l

"

.1 1 '\. ,, "., , ,

\
\

\

\

,

\

.....

\

\

"

....

,

" . .... . .

•

,

•

, " .. .. . .. ..

1 1 I II I ,I' " 1,

•

...

,

I

I

..

....

,

I

I

I

f

10

..

-

....

.

,

",.

to

I'

.... .. . . . . .. .
. . . . . . . .

.... . , ' ----_._-----...., .,. .........
...
.. , ' " ~.----~--------_-----------,----_._------.. " ,_._--_._------.......

......

••(c) XJ= 0.5

X2

•
••

-------.. " ;".~.~.----...-.•. , I, ,,,.,.~.__- - -..- - - - - - -------------.. I ' .".",;.,- _ _- - . - - - - - - -...
...."" , .......
, , - ---_._----... . . .
- - - - - ---------- -_ ...
~Figure 2.12: Velocity vectors, Re = 100. on the planes (a)x2 = 0.5, (b)Xl = 0.5 and (C)X3 = 0.5

44(a) x 2= 0.5

X3:.

.:...; .:.

; .:.

.

.

.

. ...

.

:":\\\z\\\\\Ill ;;::' j' ~IU
llll!!r:,::::;~/~z1.~~1,1. z~ ~~ l\~ \'''~~:::
"

~

;

Z f.

::::~~IIIrf;II./'/,/,//~~ I~!.,.dll////;,11111/////%/ /•,llfZ\\\\\1 •.
\\\\10:\, \\\1 .

/ / / / / / / / I r \ \ \\."''''''.. ,\\11'~\ \11 ,,IIt!!.1 ///////'/' / I r I \ \ \ " ' ' - ,,III/ / / / / , / / / / 1 1 I 1 " , ............' ' - "Idllili/////---,III/ / /' ..... -- __ '" ,
I

~

"II

X3

n~I~:\~~\~ ;~~ii::

\\ \
\\\\\
\ \ \ \"""// // / ! I,III I !

(b) XI= 0.5

)

"II

,1111
,1111

.1

I

I

\

\

"

.... _ _ _, "\

I

,

,

,

,

_ _ - _ ..... ,{// / ___ . . . ,_-\\___ ... __ . . . '\\ \,'• , ;' / '" ~ _ , \ \

\

\

\

I

I

I

I

I'

\

\

\

I

I

I

I

I

•,11 \ \ \ \ \ \ ,_

I

\

•

\

,

... _~:;::.oI'I\\\"

,

I

I

\

I

\11 ,

1\11,

\ 11 11,tllll'I,I I
11111 •
_ , I I 1111 II": : ; ; ~:_"'I !lIII",":::!!~~~~~::: ~ ~ :

'

\

\\1.

1 II!,\\\ 'r I ;' _ , '" " "II/ I •
\ I I,'\ \•.

~ ~ ~ ~ : :::~ ~~~~~::::::::~~~~~~:= =;: ~ ~~ ~ :::>~nH ~~ ~ ~ ~ ~==:~~~~~~:: : :~ ~ ~ ¥~~ ~;;;:; :: :::::::: :::::: : :,..... . .. . .. .. .... .,.,,,,.
,.,-..-"""--- ...
,.., ,,,. ""......-.--. .
:::::::: ;: : : ..
.. " , '" ,. ,... //",~~---~....---, , "',. ///-------"'-----------.;~---.,..".

~~• ••• " ' 1

I

\

~

...

_,.,.

.",

;..-os.----/////"."".,..".__.•... ,1' \ \ ..... ___ ;' /////",.",,.,....,,..---=.........
, _--.••• 1' \\ \" .. ___ .-//////./...,..,----._._ - - •••• , 1 \ ' "...... _ _ _ _ _ "••••• 1 \ "... _ _ _ _ _ _ " " " . " . . , . / / / " ' " . . . . _ - - - -. / . / ".,.",./ . " " , - " " - - . . _ - - - ----_._ ----..... \\' ...... _---- . . . . .//"",..""..,..,.. *""'...,..,-..--._ - - - -.... .. ...... -------- --~................

,."-'" -'"(c) X3= 0.5• • ••••"... _ _ _ _ _ _ '"'"-0.

...............

~~.......

_ _ _ _ _ _ __--------------__. _--_----_._---.... ,"""------~, ............................... "-~-----X2.... ,"'''''----_-..., ............... """~---'"-- - - _....................................
.................. , ..............
" , ,- - _ .......... , ...... ,:'- ........

..........

-.....:.--::: : ~ : :. :. : : : ~ ::: ::: :::::::::::::::;~--:::::====, ........................................................---~...

. . - . .

- .

: ~ ::: :: ~:;-:~:::~-~:--~-~~-~--~- -------- -..

.

-

.. .

.

..

Figure 2.13: Velocity vectors R e = 400, on the planes (a)x2 = 0.5. (b) Xl = 0.5 and (C)X3 = 0.5

45(a) x 2= 0.5

x)

L

Xl

~1\~~~~\'7} \~ \\~ \\ i I; II; /' II i {: j: I: l~l~r~rr~::~I:f..... ' , I \I I:: :: :: ~ ~ ~ / I I I I

I\.. ... ,////11/1 I I I I 'L\ •. ..,,,(1 1\ \ \ \ \ \ \ \ \ \ " •. :::( 1 I \
\ \ \ \ \ \ \ " " ....

,

r

I

\ \ \ " " " " , .. ..

r

1 I \ \ \ " " " , ' \ \ .. ..

I

I
I

I

I

I

I

.. , 1 1 1 1 / / / / / " " 1

I• ,1111////// .. "I"

\ \ \ \ \ \ \ " " ...

"I

I

/

1 1\

,r ( 1 1\ \ \ \ , \ \ \ " " , ...

... , , / / / / / / / 1 / 1 1(b) Xl= 0.5

\

I'

I

....• ",1// / / / / I I...• , , / 1 1 / / / / / I I... • d l i l i / / / / I I... , / / 1 1 / / / / I\~ ~ ~ I , ( 1 1 : : ~ \ \ \ \ \ \ ~ : : : ...

r
r

\

\

\

\

\

... " " " " ' \ \ " ..

, ... " " , ,'\\\ ...

X2"'I~III:"_::~~~j\)};: ~;:::~;:::~~::::,{((, ~ .. \ '}I I ; ' / / / / - - - s::.. ,'::" \, ~; ; ; I / ~ ~ ~: ~ ~ -;; -:;-:;-:;-~-----\/ / / // / //", ..-'''', :: , - .
.. ,,' \, " - ,.,,1 I\ , _ ,

\

/

",

/
/

/

/

/

I

/

/

/

/

I'

""'\\\"--"//111"_---... , ' \ \ ' , .... _ - " //.... ,\', ...

~

/

-

JI'

. . . . . ............ , - - -II"/ / ", ", ........
------.. - - - - - - - -

" , ,
"

I

;'

(c) Xl= 0.5.. .... ",,--- - - ....

.....

,"'~--_ ....

X2

L

Xl

Figure 2.14: Velocity vectors R e1000. on the planes (a)x2

460.5. (b )Xl

0.5 and (C)X3

0.5

:, \\~~ \ r. Zll.l.~ ~;..::z/"."....,-_--________~ -;::;:;::::::::=---"""'1~... ,,'
Zl.%~/// ..~-------------~\::::: \: zz'I.'/. '/ ~ ~ ~ -;~ ~ ~ ~ -:: -:::: =----' . ~~\,..of III, I.z!'z Z!. /,'/.1, /1 / / /' .,..
~ ~ ~'.• ",11/./.",,1.
/ / / ' " .... - - - .•• ," 1/1 Z'f
1.( / I / /' " .- .... __ : : = : : : : : 'I,;0.. • I-

Z/. I.

.",IIIf,"."
1111 ,.

-

-

-

-

Ii

:: :

1

)

,

##_-"':: ~:I'llI,11111'"I I I I / /::::~~~!! IIf{~I.,If II I.,111111\I(a) X2= 0.5-

I,'[I., / I I " '" '" .. - - tI / I.

"
tt
1
lYJ

\\::!jIU

;;'
I,

I

,

'

\

I

,

•

'

I

I

IJ

J /

I
I I I 1J j
\ \ , \ , \ , . . . . . . : : ; : 1/ J ,J ::,tllllff \'\\\\ \ \ \ \, . . . . . , "III~,;7!l. ~·1I",til" 1/ .\ \ \ \ , \ , ,• •II/I UJ/;f, I".,1111 II \ \ \ \ \ \ \ , , , ,. • , / ; : / /,,ij!~ ~ •.

,tllIlIl

I

Xl•,,1111/ \\\•• ,,"//11\

L... ,"111\\~ ~~///~/'I»IJt;j'/I"
1)/. /1" "\ \ " " " " ... , "'-"",.----.-/'/" " - - - - - ___~//:"--.'------~/::---..._____. ..

,,",'\\........
,'1\'-/1"/....... ,' \ "~----.;:::::/ ,.,.
,. ..

/~~ :

.. ..

:::::~~~~~~~~~,(~P •• ·•Xl... . _--_ .. .

.\\. \,

L-,~~~"~~~~~~:~\

I.... ,. ,11111

:

I

I

I

.... ' , ' 1 / I.. .. ' , " I I I I
I I

,
,

I,
I
I

., ••• ~ , I I I I III. .• •• - , ~ I! II {' rn:-,,:~,: :: \if!;:'!!::

,

"

I

I
I

I
,

,
,

I

I

I

\
,

I

I

I

I

I

I

I

I

,
I

(b) Xl= 0.5,., III"." ____

,

I

I

I

I

I

I

I

I

,

, \ \ ' \1 \" ... ..
,\\ 11\ \"... ..
~

I• • ••• - - " 1 1

I

I

\

\

\

\

\

\ \ , , , .... .
" ' , ...... .

I

..

Xl

L

X2
::;:.:. :. :.::..=-=...=- =- =.. :.. :. :.

~ ~ ~~ -=~= :::~~~:~~~~---::~~ ~ ~. . "\ \';:::::~~~:=-, ,•• 11" I , •. • " " ; J ; I•• 11", • . .
II/I//~,.,,, I, , • .
"" I I I • • •""""l ..

••• 'f\\,.....

... . ........ ,......

--'tl'" ., , , .... -(c) Xl= 0.5

I

I

I

/

I

I

I

/

,

I
I

I
I

I
I

I
I

I

I

,

- " .....

,

I

- ,

.• ,1" _. . , , __, •• ,

f.

~ ,

I ,

"1, ""

L

1......

,, -''''''" ...
·t"'l
l

l

'

<00

..

, \ ,
\

\

\

•

, ,

\

\

, ,

, , ,
\

\

\

, ,
\
\

Xl

Figure 2.15: Velocity vectors, R e(C)X3 = 0.5

2000 , t

50, on the planes (a)x2

0.5 , (b )Xl

0.5 andFigure 2.16 : Paths of marker particle. (a)Re = 100, (b)Re = 400, (c)Re = 1000

1.0

UJ

0.0 ~-----------,liJtr-I-:::..:l'b~:-----------;I

-1.0
-1.0

0.0

0

0
D.

UI

1.0

present
FusekiFukuyama
Kato et al.
Ku et al.

Figure 2.17: Graphs of Ul (0. 5, 0.5 , .) and U3(- , 0. 5, 0. 5) , R e = 400

49

1.0

U3

0.0 J--------------,~9-..:J:t:c:~oo:;;::::-------__1-1.0 L -_ _ _ _ _ _ _ _ _ _~.-------------'
-1.0
0.0
1.0
UI

o
6

present
Fuseki
Kato et al.
Ku et al.

Figure 2.18: Graphs of ul(0.5, 0.5 , ·) and U3(·, 0.5 , 0.5) , R e = 1000

50

Figure 2.19: Three-dimensional flow problem around a circular cylinder, boundary conditions (u = 0
on the cylinder surface) and a divi ion of the domain
The sizes of the computational domain D are Lin = Lout = vV = 7.5 and H = 1. The finite
element subdivision of the domain i constructed by dividing the domain into a union of triangularprisms and further subdividing each triangular pri m into three tetrahedra. In this process, a triangularsubdivision of [22 (the total element number i 2.J60 and hp = 0.00 ) and the division of (0, H) into mequi-intervals (m = 16) are combined ( ee Figure 2.19 ). The subdivi ion obtained is symmetric withrespect to the planes X3 = iH / m (i = O..... m .- 1). w h re ymmetric solutions with respect to thoseplanes may appear. But the subdivi ion i not uniform in the x3-direction because of the characterof the tetrahedral element (see the decompo ition of thide wall in Figure 2.26). We examine theeffect of this fact later.
The two parameter of the scheme arc .' t as T = 0.00 and 0' = 2. Figures 2.20,2.21, 2.22 and 2.23
shows the velocity at a specific time on horizontal plane. X3 = 0.0.0.25,0.5 and 0.75, respectively.
Figure 2.24 depict the velocity on a vrrtical plan ..r2 = O. at the arne time. We can observe that
there appear vertical velocity components on many pointb hind the cylinder and that horizontal

flow patterns are slightly different depending on the planes.
I n order to observe three-dimen ionality qllantitati\' 1y. we computed the deviations from mean

.51

~-

values in the x3-direction,max(max lUi .- Uil)/ max IUjl,
02

02J

X3

i = 1,2,3,

max(max Ip _. 151)/ max 1151,
02

02

X3

(2.81)(2.82)

where 'Uj and 15 are the mean values in the x3-direction at each point in O 2 , e.g.,

Io.If

V = -1H
VdX3 .

0

We show the time history of e in the upper graph of Figure 2.25, where the vertical axis is measured
logarithmically. The deviations e remain in alm.ost the same levels until t = 30, grow in t = 40 rv 45
and are almost stationary after t = 45. As mentioned above, the finite element subdivision of thedomain is not uniform in the x3-direction. To examine the effect of the nonuniformity, we also computevalues e' defined by the nodal points lying on the even number layers, that is, equations (2.81),(2 .82)
are evaluated only on the layers,

X3= i/m,

i = 0,2,4 ... , m.

They are depicted in the lower graph of Figure 2.25. The deviations e' are pretty small (rv 0(10- 5 ))until t = 10 . Then they grow up exponentially and become almost stationary (rv 0(10- 1 )) after t = 45.
The values e and e' are almost the same after t = 35. We recognize that the three-dimensionality
caused by the non-uniformity of the subdivision is dominant until about t = 30. Later the difference isovercomed by the genuine three-dimensionality. Thus, the instability of the two-dimensional solution
at Re = 1000 is verified.

R e mark 2.6 We also examined the difference of stability between flows at Re = 100 and Re = 1000
qualitatively using two relative rough subdivisions . The sizes of the computational domain 0 are

Lin = 2.25, Lout = 6.75, ltV = 2.165 and H = 1. The finite element subdivisions are made in themethod mentioned above (see Figure 2.26). The total element number is 678 and hp = 0.25 in thesubdivision of the two-dimensional domain O 2 . The subdivisions in the x3-direction are m = 4 and6. The two parameters are set as T

= 0.04 and a = 2.25 . Table 2.3 shows values e at sufficiently

elapsed times . \Ve can see that in the case of Re = 100, these values e decreases as m increa es, whichsuggest that the finite element solutions converge to a solution uniform in the x3-direction. On theother hand, in the case of Re = 1000, these values e increases a m increases, which suggests thatthree-dimensionality appears.

52

Table 2.3: Deviatiions from x3-mean values

Re

m

t

100

4
6
4
6

40
40
40
90

1000

e1
4.04 x 10- 2
3.57 X 10- 2
4.09 X 10- 2
3.66 X 10- 1

e2
6..44 x 10- 2
6 . 33 X 10- 2
6.95 X 10- 2
3.35 X 10- 1

e3
2.32 x 10- 2
1.86 X 10- 2
2.23 X 10- 2
1.49 X 10- 1

e4
5.11 x 10- 2
4.05 X 10- 2
5.22 X 10- 2
8.46 X 10- 2

Figure 2.20: Velocity on the plane X3 = 0.00. Vectors of (Ul' U2) and contours of U3 = ±0.02, ±0.04, ...

53

,....-----,

Figure 2.21: Velocity on the plane X3±0.02, ±0.04, ...

0.25. Vectors of (Ul, U2) and contours of U3

---

---

---

...........

...........

""-...
'-...

"

Figure 2.22 : Velo ity on th plane X3

0.50. Vector of (Ul . U2) and contours of U3"-±0.02, ±0 .04 , ...

Figure 2.23: Velocity on the plane X3

Figure 2.24 : Velo ity on the .r2

0.75. Vectors of (Ul' U2) and contours of U3

o planvectors of (Ul . U3) and contours of U2±0.02, ±0.04 , ...

0, ±O.l, ±0.2, ...

QQ~_ _ _ _ _i~Q~--~2~Q----

_3~Q~

Q

-3

o
6

o

-2

o
'V

o

e'1
e''1
e'3
e'<\

-3

-5

F ig u re 2.2
5: G ro w thof d e v ia ti o

F ig u re 2.2
6: F in it een s from X m e a n

le m e n t s u

3-

b d iv i ion ,valuesm =' 6

2.7

Conclusion

In this chapter, we have presented the upwind finite element approximation based on the choice of upand downwind points. It has a potential to approximate the convection term in third-order accuracy.
Combining this approximation with iso-P2 Pl/Pl elements and the procedure of a consistent discretized pressure Poisson equation, an upwind finite element solver for the two- and three-dimensional
incompressible Navier-Stokes equations at high Reynolds numbers has been presented. The scheme
has two scalar parameters , the stabilizing parameter and the time increment. We have investigatedstability regions in the a -Tplane by two sarnple problems and have studied the shapes of the sta-bility regions. Some numerical results have been shown to demonstrate that the solver can produce
numerical results possessing good Reynolds-number dependency and capturing three-dimensionalityof phenomena well.
We have also presented a blind search algorithm and heuristic search algorithms for upwind-pointsearch problems. In the search problem of U, the heuristic search algorithm shows good performanceon scalar (conventional) computers while the blind one is good on vector computers. We discussed therelation between the point-search problem of U and a point-search problem in the (d - I)-dimensional
space. In order to avoid a risk of an infinite loop caused by finite-precision arithmetic, a topologyoriented method is presented, where the Cramel's formula is used rather than the Gauss's elimination
method to solve a small problem in the heuristic search algorithms.
In order to use finer subdivisions of a domain or to compute more complicated flow phenomena,
large- cale computations are required.

Thus development of high performance implementation isimportant in engineering point of view . One effort has been paid for high performance implementationof the upwind-point search in this chapter. In Chapter 3 we consider parallel execution of the scheme.

57Chapter 3

Domain Decomposition Algorithm
3.1

Introduction

In this chapter, we discuss parallelization of the finite element method in massively parallel computersof distributed memory type (see Figure 3.1). In the computers data can be shared by communication
procedures called message-passing, whose speed is slower than that of computations in a processor.
From a practical point of view , it is important to reduce amount of message-passings.

Policy to use domain-decomposition
We decompose a domain into K non-overlapping subdomains , and assign data and computationsin the kth su bdomain to the kth processor. As is discussed in General Introduction there are two
ways in the domain-decomposition situation. They are data-parallelization[5 , 3, 36] and the domain
decomposition method[47, 21 , 51].
The later, the domain decomposition method, has been developed from two aims[7]. One is totreat problems in complicated-shape domains.

In such problems , one can divide the domain into

subdomains , apply the spectral element methods etc. in several rectangular subdomains, and use thefinite element method in other subdomains. The other aim is to achieve parallel computations.
vVe here use the domain decomposition method as a tool of parallel computations.

Notations on domain-decomposition
In this chapter v. e restrict ourselves to two-dimensional problems.
Let D be a bounded domain in R2. We divide the domain into a union of finite elements, Th ={edi=l ,...,N e - This subdivision is supposed to be regular. vVe denote by w the set of all nodal points
in D, and set Np = IwlLet {1k,h} k=l, ...,K be a partition of Th into a union of K non-overlapping subsets. Associated with58

Figure 3.1: Massively parallel computer of distributed memory typethe partition, il is decomposed into K non-overlapping subdomains,(3.1)
We assume that the partition is made uch that each subdomains is simply connected. We denote by
nk the unit outward normal vector on fJr2 k . and by Wk

Nek

C W the set of all nodal points in ilk, We set= 11k,hl and N pk = IWkl.
If ilk n ill (k # I) includes an edge of an element, \ve ay an interface of the subdomains appears.

We denote all interfaces by r m, m = L ... , .11. Let us define integers fi,_(m) and fi,+(m) by(3.2)
The sets Wk, k = 1, .. . , K are overlapping subsets of w. \\-e denote intersection sets corresponding tor m by "'1m = W,c(m)nwK+(m), and other intersection sets by ~/m = WK_(m)nWK+(m), m = NI +1 ... , !'vir.
The latter corresponds to the case where r2k n ill (k # l) do not include an edge but a vertex (orvertices) .

R e mark 3.1 The numbers of element a signed to the processors hould be the same, if it is possible

in order to balance loads of the proce ors. It is also prefered reducing the number of nodal points on
interfaces of the subdomains on each processors because of the message-passing. We refer to [46] forthe recursive spectral bisection method. that respond to such requirements.

59

~-

Organization of this chapter
The presence of overlaps between nodal points of subdomains (rm defined above) makes parallelizationof programs of the finite element method more complicated. In order to overcome that we first present
the idea of two data types for data on nodal points in the finite element method. That is useful forboth of data-parallelization and the domain decomposition method.
Secondary, we introduce a mortar element Inethod [7], that is a kind of the domain decomposition
method. We review a scheme of mortar element method[l, 2] and apply the scheme to the Stokes
equations of primitive variables.
Finally, we develop the domain decomposition algorithm of the upwind finite element scheme basedon the choice of up- and downwind points described in Section 2.4.

3.2

Two data types for data on the nodal points

3.2.1

Conventional programs and data-parallelization

Conventional programs
For simplicity, let us consider the Helmholtz equation with a Dirichlet boundary condition,fu - j.L.6.u
u

in D,

(3.3)

o on an,(3.4)

where u is an unknown function, f is a given function and j.L is a positive constant. Defining V =HJ (0,)

anda(u,v)(f , v)kkjl'Vu· 'Vvdx +

k

uvdx,

f vdx,(3.5)(3 .6)

this problem is reduced to Problem 1.1.
In order to apply the Ritz-Galerkin method described in Section 1.2, we employ the finite element
method with the triangular PI element. vVe then have a system of linear equations (1.24).
A program to solve the system consists of two stages: (i) to compute components of the matrix
and the right-hand-side vector in (1.24) , (ii) to solve the system. In the first stage, the element matrix,a;{3 =

1eJ(j.LV¢l(o: ,j) . V ¢l({3,j) + rPl(o:,j)'Pl({3,j))dx,0', ,3 = 1"',,3, j = 1"", lYe,

(3 .7)and the element vector,

bd = /. . f¢l(o:,j )dX0:= 1, ... ,m, j = L ... ,iVe ,J

60(3 . )

are computed, where {¢di is a set of the base functions of PI element, and l (a, j) is a global nodal
point number corresponding to the local nodal point number a in the element ej.

If they wereassembled over all the elements, the global stiffness matrix and the global vector would be obtained.
The program is shown in Program 3.1. In the second stage, the system of linear equations is solved bya conjugate gradient method. It has three computational components, that is, the linear combinationof two vectors, the inner product of two vectors, and the product of a matrix and a vector. Programsof the latter two are shown in Program 3.2 and 3.3, respectively.
Program 3.1 Computation of the element matrices and the element vectors

REAL A(1:3,l:3,l:Ne),B(1:3,l:Ne)
DO j=l,Ne
Compute Iej(j.L'V¢l(fr ,j)· 'V¢l({3 ,j) + ¢l(fr,j)¢l({3,j))dx(a,f3 = 1,··· ,3)and store them into A{1:3,l:3,j)Compute Iej f¢l(fr,j)dx(a = 1,···,3) and store them into B(1:3,j)

END DO
Program 3.2 Computation of an inner product, a = (b, c)

REAL A,B{1:NP),C(l:NP)

A = 0.0
DO i=l,NP
A = A + B(i)*C(i)END DOSum up A over all the subdomain in domain-decomposition cases
Program 3.3 Computation of a matrix-vector product, c = Ab

REA L A (1 : 3,1: 3,1: N e) ! the element matrices
REAL B(l:NP), C(l:NP) ! vectors
Clear C(l:NP)
DO j=l,Ne
DO a=l,3
Add L~=lA(a, (3,j)*B(l(f3,j)) to C(l(a,j))
END DOEND DO61

Data parallelization

We consider here data parallelization of the above program in a massively parallel computer of distributed memory type (see Figure 3.1).
If all the Ne (resp. NP) are replaced by Nek (resp. NPk) in Program 3.1, the subdomain stiffness

matrix (a N pk x Npk-matrix) and the subdomlain vector (a Npk-vector) are computed and stored inparallel in every subdomain.
With similar replacements, Program 3.2 and 3.3 can be executed in parallel. In order to makethe results valid, we need to take care to handle the data, because nodal points on interfaces of the
subdomains are shared with the neighbor subdomains. If the vector b has the same values at theshared nodal points in all the subdomains, Program 3.3 can compute the matrix-vector multiplication
correctly, but the product c has different values at those nodal points in the subdomains. In fact ,Program 3.2 can compute the inner product of such a product vector and another vector which has
the same values at the shared nodal points correctly. We can see that there are two types for the data.
It is useful to distinguish them to simplify the implementation of the finite element method based on
the domain-decomposition. We formulate the two data types exactly in the next subsection.

3.2.2

Two data types for data on the nodal points

Assembly mapping and restriction mapping

We consider the lVp-dimensional real vector space defined on w,
V={v;v:w---+R},(3.9)

and the N Pk-dimensional real vector spaces defined on Wk,
Vk = {v;v: wk ---+ R},(3.10)

for k = 1" . . ,K. We define by V* the product space,(3.11)
which is a (Lf=l lVPk)-dimensional real vector space. For these vector spaces, the usual inner products
are introduced. For example , for u* = (Ul,'" ,UK),V* = (VI'" ,VK) E V*,K(U*,v*)v*= L(Uk,Vk)Vk,
k=l62(3.12)

where(3.13)For any 1 :S k :S K, we define by rv: Vk ~ V the zero-extension operator, that is, for Vk E Vk ,(3.14)We define two mappings between the spaces V and V*. One is an 'assembly mapping' A : V* ~ V,that is, for v* = (VI, ... , VK) E V*, the zero-extensions of Vk E Vk are summed,
K(3.15 )

A v* = LVk'
k=1

The other is a 'restriction mapping' R : V ~ V~ defined by(3.16)The range of the mapping is the subspace of V*,(3.17)
Proposition 3.1 (i) A is a continuous linear mapping from V* into V. (ii) R is a continuous linear

mapping from V into V*. (iii) R is an isomorphism from V onto V*.
For computation of an inner product, we have,
Proposition 3.2(3.18)

Remark 3.2 In general,(u*,V*)v.

(u*, v* E V*)

do not holds.
In order to consider matrix operations, we also introd lice matrix spaces in a similar manner as
above

v

{H;

H:wxw~R}(3.19)( k = 1 , ... , K) ,(3.20)(3.21)

63

We can also define the assembly mapping A and the restriction mapping R for these spaces. The
range of the restriction mapping is the subspace of V *,H K_ (m )1 ITn X ,Tn -- H K+ (m )1 ITn X ,Tn ' (1 :S m :S M 1 )} .

(3.22)(3.23)
It follows that
Proposition 3.3
A(H* . v*)(AH*)('R,-lV*)

A(H* . v*)(R- H*)(Av*)

1(H* E V* ,v* E V*),

(3.24)(H* E V*, v* E V*).

(3.25 )

Data of the pre-assembly type and the restriction type

In an implementation based on domain-decomposition, no processor can handle data on the entire
domain, that is, we cannot do any operation in V or V. But it is possible to treat subdomain-wise
vectors and matrices, that is, vectors in V* and rnatrices in V *. There are two types for data. A vector
v* E V* is considered of 'pre-assembly type' and is identified with Av* E V. A vector v* E V* Cv:

can be considered of 'restriction state' and be identified with R -lv* E V. In the following of thissection, we denote by subscript '*' the quantities in the pre-assembly type and by superscript '*' thequantities in the restriction type to distinguish them clearly.
Proposition 3.2 and 3.3 suggest operations in V* and V* equivalent to operations in V or V. The
prod uct of H* E V * and v* E V* gives the corresponding product in V and V by Proposition 3.3.
These operations do not need message-passings. The inner product of two vectors in V is equal to theinner product of corresponding vectors in the space V* by Proposition 3.2, provided that one vector
is prepared in the pre-assembly type and the other is prepared in the restriction type. After products
in all the subdomains are computed, a summation is needed over all the subdomains. It needs amessage- passing over all the processors (we call it global message-passing).
If we need to compute the inner product of two vectors in the pre-assembly type, one vector v* E V*

is transformed to RAv* E V*. For this transformation, every processor collects data at nodal points on
the interfaces from neighbor subdomains and su:ms up them. This operation needs message-passings
among processors sharing interfaces (we call thern local message-passings).

643.2.3

Implementation of the conjugate gradient method

Let us consider to solve a system of linear equations, Ax = b, where A E V is a symmetric positive
definite matrix, b E V is a known vector and :X E V is an unknown vector. A procedure of a preconditioned conjugate gradient(PCG) method is described by Algorithm 3.1. It is required that A
and b are given in the pre-assembly type and that a pre-conditioning matrix P ~ A is available in the
restriction type. The solution is produced in the restriction type. This procedure needs one local andtwo global message-passings in an iteration.

Algorithm 3.1[PCG method]

x* := an initial vector;

while p ~ E :x* := x* + Q7f*;

We next consider a system of linear equations(~(3.26)

We note that a problem with constraintsBuj,(3.27)

0,(3.28)

is a special case of (3.26) .
If A is regular we obtain by elimination of x in (3.26),(3.29)
After y is solved by this equation, x is computed by(3.30)
65

Assuming the matrix C - BA -1 BT is symmetric positive definite, we can solve (3.29) by Algorithm
3.1. Imbedding the computation of g- BA -1 f in (3.29) and the computation of (3.30) in the procedure,
we obtain the following.
Algorithm 3.2 :y* : = an initial vector;x* .=. A *-l(j* - BTy*).
*,r* := B*x* - (g* - C*y*);7r*:= RA(P*-l r *);p:= (r*,7r*);

while p ~ E do :y* := y* + a7r*;
x* := x* - aw*;

It is required that A, B, C and (j, g) are given in the pre-assembly type and a pre-conditioning matrix

P is available in the restriction type. Solutions (x, y) are generated in the restriction type.

3.3

Mortar element method

In data-parallelization of the finite element method , there may remaIn two bottlenecks to obtainhigh performance of parallel computers: (i) Message-passings between the processors occur at every
iteration step in the CG solver.

(ii) Even if two subdomains share only one nodal point, whichcorresponds to 1m, m = !vI + 1, ... ,M1 , a local rnessage-passing appears. The domain decomposition
method has the potential to overcome the bottlenecks.

3.3.1

Review

Here, we review a mortar element method[l, 2], that is a kind of the domain decomposition method.
Let us consider to apply a mortar element method to the two-dimensional Helmholtz problem(3.3),(3.4). In the mortar element method, interfaces of subdomains should be straight segments, but
the total numbers and the positions of the nodal points on a interface may differ between two subdomains sharing the interface. (In fact , the basic idea is to divide the entire domain into quadrilateralsubdomains first , and divide each subdomain into a union of elements independently.)66

...---We introduce a Hilbert space,with a norm //vl/y. ={Lf=ll/vkl/~l(nk)} 1/2. y* is equivalent to HJ(D). We therefore obtain the

weak form:

Problem 3.1 Find U E y* such that
(3.32)
This is equivalent to the weak form discussed in Section 3.2.l.
The Galerkin method is adopted for the space discretization. Each subdomain Dk is divide into aunion of quadrilateral elements: ~ , h. Defining finite element function spaces subdomainwise by(3.33)
we construct the finite element space on entire domain byK

II Xk ,h,X h ==(3.34)

k=l

where Q1 is the space of isoparametric bilinear functions.
If T",±(m),h has more vertices than T"'~(m),h OIn r m, we define(3.35)
where Trk : H1(Dk) - 7 H1 / 2(8Dk) is a trace operator. If we omit the degrees of freedom at bothedge points of r m, Wm,h reduces to Wm ,h (for the base functions of this space see Figure 3.2). vVe set
HAm =dim VVm ,h. We adopt their product space,M

vV; = II Wm ,h,(3.36)

m=l

as a finite element function space for the mortar function.
Let a* : ";'(h x X h - 7 R be a symmetric bilinear form,K

2:= ak( k,h, Vk,h),(3.37)

1L

k=l

J11n V~Uk,h· VVk ,hdx + Jnr Uk ,hVk ,hdx ,
k

k

67

(3.38)

Subdomain A

Figure 3.2: The bases of the mortar function space
j* : X;' x W; ---+ R be a bilinear form,M

L {jm ,+(VK:+(m),h, J-Lm ,h) + jm ,-(VK:_(m),h, J-Lm,h)}- 1fr (±VK:±(m),h)J-Lm,hdS ,(3.39)

m=l(3.40)

mand f* : X;' ---+ R be a linear functional,(3.41 )
k=l(fk, Vk,h)==10r fkVk ,hdx.

(3.42)kRemark 3.3 Substituting (3 .40) into (3.39), we obtain

L Irfm (VK:+(m),h - VK:_(m),h)J-Lm,h ds .
M

j*(Vh, J-Lh) = -(3.43 )

m=l

We can define the discretized space corresponding to y* by(3.44)
A weak form of a discrete version is
Problem 3.2 Find Uh E Y h* such that
(3.45 )
Referencing Theorem 1.4, it is actually treated as a saddle type problem as follows. Let A* : .Xh ---+eX;')' and J* : ..Y;' ---+ ("vV;) , be linear operators corresponding to a* and j* , respectively,(3.46)(3.47)
Then

Yt = {Vh E ~Yh; J*Vh = O} holds and the weak form is rewritten by6

Problem 3.3 Find (Uh' Ah) E Xi; x ltV; such thatin (X::)',(3.48)o in (W;)'.

(3.49)

fh

We represent elements in Xk,h and vVm,h by
N pk

L Uk,i¢k,i (E Xk,h),(3.50)

i=l
Nt\tn

Ak ,h

=

L .Am,iwm,i (E vVm,h),(3.51 )

i=l

respectively, where ¢k,i, i = 1, ... , N pk and wm,i, i = 1, ... , N m are the base functions of Xk,h and

W m,h, respectively. The system of linear equations obtained from Problem 3.3 is(3 .52)where notations areu(3 .53)(3 .54)A(3.55)(3.56)

F(3.57)(3 .58)

A is a block- d iagonal matrix,(3 .59)and each block is defined by(3.60)The matrix J, called the Jump matrix, is a sparse block matrix,

J =

[ 0

... 0

compo nents of each block are

Jm,K_(m)0(3 .61)Equation (3.52) can be solved via eliminating U,(3.63)
with the CG method. In its procedure, the mult iplication A-I G appears , where G corresponds to anelement in (X:J'. Since each block of A corresponds to a Helmholtz problem in a subdomain with a~eumann boundary condition on the interfaces, the multiplication just solves a sub-problem in each

subdomain. It can be executed in parallel.
The quantity Am,h E VVm ,h corresponds to the flux by diffusion ,(3.64)

Comparison with data-parallelization

Here we compare the data parallelization and the mortar element method in a one-dimensional
Helmholtz problem,

u-u"

1,

x E (0,1),(3 .65 )u(O)0,u'(l) =0.

(3.66)We divide the domain into a union of 4 equi-Iength subdomains, and further divide each subdomain
uniformly into 25 elements . In the both methods , the systems of linear equations are solved by the
CG method.
Since there are 101 nodal points, the CG method converges at most in 100 iterations in the dataparallelization. In each iteration step, one local and two global message-passings occur.
The mortar element method has two nested CG loops. The outer loop is used to solve (3.63) andthe inner is used to solve A.-I . The outer loop converges in 3 iterations. Iteration times of the inner
loop in each subdomain are shown in Table 3.1. The shape of the solution at each step in the outer
loop is shown in Figure 3.3. The total number of iterations are almost the same as the former method,
but the number of iteration with the message-passing is only 3.
Thus, the mortar element method has advantages: (i) Independency of problems in subdomains isrelatively high and we can reduce the amounts of message-passings , (ii) Two subdomains which hare
only one nodal point do not need to communicate each other by the local message-passings. Therefore:
the mortar element method has the potential to perform more efficiently in the parallel computers.

lll~

0.3

·············,················t········X-·············t ............. -::0.2····I·u . ; ~._ . . u.

~

0.1

predicted~n~~ ~~bl?~c5'p

converged

o I!S-"..............I.-..I.--L..~~-L.....I'--'---L..L.....J.-L...L..J--L..J......J

o

0.4 x

0.2

0.6

0.8

Figure 3.3: Shapes of the solution at each step in the outer CG loop (u = 1.0 in the interval where
they are not plotted)
Table 3.1: Iteration times of the inner CG loop
CG step
prediction
1
2
3

3.3.2

11
25
25
25
25

subdomain
12 13 14
1
1
1
26
0
0
26 26
0
26 26 26

Mortar element scheme for the Stokes problem

As a application of the mortar element method to flow problems, a stream function-vorticity formulation is proposed[l, 2]. They reduced the problem to two Poisson problems and two Helmholtz
problems. In this formulation, however, it is not easy to treat stress boundary conditions. We here
apply the mortar element method to the Stokes problem in primitive variable form. We impose thecontinuity constraint only on the velocity.
For the velocity and the pressure, we adopt the triangular iso-P2 PI/PI elements. The discretizedfunction spaces on each subdomains, .Xk ,h and Qk ,h, are constructed on the subdivision.

By thedefinition of Xk,h, the function space for the m.ortar function ~Vm,h is made in a similar way as in
Section 3.3.1. The function spaces on the ntire domain are constructed by
K.X"h = II .X"k,h,
k=l

K

Q

h = II Qk ,h,
k=l

M

~Vh =

II VVm,h,(3.67)

m=l

for the velocity, the pressure and the mortar function. respectively. For time discretization we employ
an implicit scheme, that is, we evaluate the time derivative term by a backward difference,(3.6 )

71

and evaluate other terms at time nT.
Consequently we have the following problem. in matrix form:(

Ao +J!vI IT

JT BT) ( un)
0
0
An
BOO
pn

(pn
)
0,

(3.69)

0where M and Ao are the mass matrix and the viscosity matrix concerning Xh, respectively, B is thedivergence with respect to X h and Q"h, and J is the jump matrix with respect to X h and Who
Elimination of un gives(3.70)
where A = Ao + NI IT. An algorithm to solve this system by the conjugate gradient method is shownin Algorithm 3.3. In the algorithm, multiplications of the matrix in the left hand side of (3.70) and avector corresponding to W h x Q"h appear. It can be computed in parallel without message-passings.
Each block of A is symmetric positive definite even if 80. k n 80. = 0.
Algorithm 3.3 The CG method of a mortar scheme for the Stokes problem), *, p; {initial vectors}

u := A -l(f - J;)'* - BT p),.

r)..*:= J*u ,·rp:= Du,·

while p 2: f :W := A-1(J;),* := ),*+ 07r;,·/3 := J-LI p,.

Remark 3 .4

7r; + BT 7r p),.

q)..*:= J*w,·

p:= p + 07rp ,·u := u -qp:= Dw ,·OW,·

P := J-L,.

The idea of two data types discus ed in Section 3.2 is applied to the mortar functions

and the jump matrix. Scripts '* in Algorithm 3.3 show a type of data a in Section 3.2. For thevelocity and the pressure, we do not need to introduce the idea of two data types.
72

Remark 3.5

This scheme is also applicable to the ='Javier-Stokes equations if the convection term

is treated explicitly in the time discretization. The convection term contributes to the vector pn in(3.69) in such a case.

3.4

Domain decomposition algorithm for the upwind finite element
scheme3.4.1

Domain decomposition algorithrn

In this section we present a domain decomposition algorithm of the upwind finite element scheme
described in Section 2.4.
Defining iso-P2 PI/PI finite element function spaces for the velocity and the pressure subdomainwise by
Vie E (pl(e))2, e E Tx:,h/2, V = 0 on aD k n f D },{v E (C(Dk))2;{q E C(D k );qle E pl(e), e E 1k,h},(3.71 )
(3.72)

respectively, we construct the finite element spaces by
K

1/;*

h

II Vk,h,=(3.73)k=lK

II Qk ,h .

Qh =(3.74)

k=l

Concerning weak continuity of the velocity between subdomains, we employ the Lagrange multiplieron the interfaces.

For the discretization of the spaces of the Lagrange multiplier defined on f m(1 :S m :S lv'I), we investigate three cases:
Case 1. The conventional iso-P2 PI element, that is defined by(3 .75 )where(3.76)
Case 2. A modified iso-P2 PI element having no freedoms at both edges of interfaces.
Case 3. The conventional PI element, that is defined by(3.77)
The finite element space VV;: is defined byMw; = II VVm,h.
,m=l

73

(3 .7 )

Introducing the Lagrange multiplier, (2.51),(2.52) are transformed toProblem 3.5 Find (u~+l,Ph' Ah ) E (Vh*' Qh' W;) such that

(3.79)\/qh E Qh,b*(u~+l, qh)

0,(3.80)\/P,h E liV;,j*(u~+l,P,h)0,(3.81 )

where

1
L1
K

(u, v)*L

k=l

Uk' Vk dx ,(3 .82 )(Wk' graduk)vkdx,(3 .83)D.kK

ai(w,u,v)ao( u, V)

k=l

D.k

2

K

(

2

-R e k=l
L) j=l
L Dij(Uk )Dij(Vk)dx,(3 .84)D.kKb*( v, q)-{; kkMj*( v, p,)(3.85 )

qkdivvkdx,

1- L r (Vx:+(m) - vx:_(m))P,m ds ,m=l(3.86)m

L(j f· Vk dx + lK(j*,v)

k=l

D.k

8D.k

nr N

9N . Vk ds ),(3.87)(,)h denotes the mass-lumping corresponding to (, )*, and a?* is a trilinear form which corresponds toapplication of the upwind finite element scheme based on the choice of up- and downwind points to

The quantity Am,h corresponds to () . nx:+(m) ""'1m 'Remark 3.6 (Parallelization for the up- and downwind point search)

In order to evaluatea?*(u h u h,Vh), we need to find two upwind points and two downwind points for each nodal point.
vVhen the domain decomposition method is used , some of these up- and downwind points for nodal
points near interfaces may be included in the neighbor subdomains. In order to treat it, each processor
corresponding to a su bdomain has geometry information of all elements which share a point at least
with the subdomain. The proce sors exchange each other the value of u h before the evaluation of theconvection term. Hence the evaluation is completely parallelized.

We rewrite this problem by a matrix form as follows.

Tn
M BT J ) ( U +
BOO
pn
(
J
0
0
An

1) (
=

Fn )

0,
0

(3.88)

where !VI is the lumped-mass matrix which corresponds to the form (, )'h' B is the divergence matrix
which corresponds to the form b* , J is the jump matrix which corresponds to the form j* , and Fn isa known vector. Eliminating U n + 1 from the equation, we get a system written by(AJ f)(pn)
I An
/vt

(~),

(3.89 )

whereA

BM - lB T ,J

JM-lB T ,M

JM-lJ T ,:F

BJVI- l F n ,()JM-1F n .

The matrix A is a K - block diagonal matrix. Each block corresponds to the matrix of the consistent
discretized pressure Poisson equation in a subdomain. The matrix J can be considered as a kind ofjump matrix with respect to the pressure:(3.90)The system (3.89) is reduced into

(3.91)
Applying CG method to this equation, the procedure can be executed by the domain composition
method.
R e mark 3.7 When r N = 0, the pressures should be sought in(3.92)
which corresponds to L6(D). Let us define(3.93)
75

If we eliminate Ah from Problem 3.5, it reduces to a problem of Problem-2.4-type. The problem hasa solution (Uh,Ph) E Yh* x Q-;;'o' (If b* satisfies the inf-sup condition with respect to

yt x Q-;;'o, thesolution is unique.) Substituting this solution into Problem 3.5, we can obtain Ah. The space,(Uh,Ph, Ah) is a solution of Problem 3.5}(3.94)

is a proper subspace of W; and the codimension of WhO in W; corresponds to that of Q-;;'o in Q-;;'.
For the simplicity of the procedure, we seek a Lagrange multiplier inof Problem 3.5 is not regular, a solution (u~+l, Ph, ~h) E

W;. Although the matrix

yt x Q-;;' x W; can be found by virtue of

employing the conjugate gradient method. The obtained pressure is Ph t/: Q-;;'o in general. In order toobtain the pressure in Q-;;'o, we therefore need to normalize the pressure and the Lagrange multiplier
by

c

Ph
A~,h

K

.

k=l

nk

-(~= jPhdx)/measD,

Ph + C,(3.95)(3.96)

Am,h - cnK.+(m)1 :S m :S lvI.

(3.97)

It is easy to verify that if (u~+l ,Ph, ~h) is a solution of Problem 3.5, so is (u~+l 'Ph' Ah ).

3.4.2

Numerical results

At first we compute a test problem. Let D = (0, 1) X (0, 1) and r D = aD (r N = 0). The boundary
data is(3.98)

the external force is(3 .99 )the ini tial function is(3.100)and the Reynolds number is set to 400.
The exact solution, which is stational is(3 .101 )xi +x~ - 2/3.

(3 .102 )

/
//

:/ /

/~
/~

/,/

:'
~

~:~/
//
~
~/:'

:~;i
/)///

<: m
~: :: :
, , , ,
..
..
~

,

,

"

,

"

'\

,

,

"-=-"O""';"~"'-"-J

Figure 3.4: Domain-decomposition( 4 x 4) and the triangulation( n = 32)
We have divided n into a union of uniform 1'1 x n x 2 triangular elements, where n = 4, 8, 16 or
32. vVe have computed in two domain-decomposed ways, where the number of sub domains in eachdirection is 2 or 4. Figure 3.4 shows the domain-decomposition and the triangulation in the casen = 32 and 4 x 4 subdomains.
Since the velocity in the exact solution does not belong to the finite element space of the iso-P2
PI element, the initial condition in the discretized problem is given by a projection (1.62) and thenumerical solution is expected to converge to the stational solution in time-marching. If
max I ukn i-Ukn-11/i
T < 10-.'5
k ,l

'

,

(3.103)is satisfied, we judge that the numerical solution has converged and stop the computation. We denoteby tconv the time when (3.103) is satisfied . The two parameters of the upwind scheme are set as
a = 1.0 and T = 0.24/n.

Table 3.2 and Figures 3.5, 3.6 show errors bet.ween the obtained numerical solutions and the exact
solutions . They are measured by(3.104)(3.105)Results of the non-domain-decompo ition case are also listed lil the table.

We can observe thatthe errors of the velocity and the pressure realize the optimal onvergence rate of the iso-P2 PI/PI

77

elements regardless of choice of"VVm,h. In the case where W m,h is generated by the conventional iso-P2
PI element, the maximum error of the Lagrange multiplier does not converge to 0 when h tends too. It may indicate the appearance of some spurious Lagrange multiplier modes, since the degree of
freedom of the Lagrange multiplier is larger than that of jump of the velocity in the choice. In thecase of the modified iso-P2 PI element for Wm,h, we can see that the error of the Lagrange multiplieris O(h) but the errors of the velocity and the pressure become larger as the number of subdomainincreases. In the case of the conventional PI element for W m,h, the error of the Lagrange multiplierconverges to 0 but the rate of convergence depends on the number of subdomains.
Since the conventional PI element has the slmallest degree of freedom of the Lagrange multiplier,it can decrease the amount of computation steps in a iteration time in the conjugate gradient solver.
Hence we adopt PI element for Wm,h in the following.
We next computed the two-dimensional lid-driven cavity flow problem. The Reynolds number is400. The domain n = (0,1) x (0,1) is divided into a uniform 24 x 24 x 2 triangular subdivision.
We chose T = 0.01 and ex = l. vVe computed in the cases of 4 x 4, 4 x 2, 2 x 2, 2 x 1 and 1 x 1
domain-decompositions. The total number of elements, the degree of freedoms, the convergence time
when the stationarity criterion was satisfied and the computation time are listed in Table 3.3. \Ve seethat the computation time becomes shorter as the number of subdomains (i.e. processors) increases,
except for the case of 2 x 1 subdomains.
The velocity vectors and the pressure contours of the computed stationary flow in 4 x 4 subdomains
are shown in Figures 3.7 and 3.8.

We can observe that the flow is captured well in th

domaindecomposition algorithm .
The computations have been done on Intel Paragon XP /S in Institute for ~umerical Simulations
and Applied Mathematics, Faculty of Science, Hiroshima u niver ity. The system has 56 processors
and each processor has 16 :vIega byte memory.

Table 3.2: Errors in a test problem

Wm,hIso-P2 PI
element

Decomposition

n

4x4

4
8
16
32
Cony. rate
4
8
16
32
Cony. rate
816
32
Cony. rate
4
8
16
32
Cony. rate
4
8
16
32
Cony. rate
4
8
16
32
Cony. rate
4
8
16
32
Cony. rate

2x 2

Modified
iso-P2 PI
element

4x 4

2x 2

PI element

4x 4

2x 2

Non-DDM

1x 1

IUh - ulv*
lulv*
0.082160.043160.018890.009355
1.014
0.10250.04368
0.01892
0.009355
1.016
0.41180.11450.03457
1.728
0.43590.17580.050540.01506
1.747
0.08697
0.04846
0.01956
0.009411
1.055
0.13200.046830.019070.009362
1.026
0.1066
0.043430.018910.009355
1.015

Ilph - pIIQ*

max P\h -IlpIIQ"
0.027760.007125
0.002175
0.0008618
1.336
0.027200.0071300.002175
0.0008617
1.336
0.012400.003676
0.0009287
1.985
0.032620.007728
0.002254
0.0008706
1.372
0.02736
0.007088
0.002168
0.0008611
1.332
0.027210.007096
0.002170
0.0008612
1.333
0.027290.007126
0.0021720.0008611
1.335

maxipi0.56890.53480.52430.5196
0.39670.38320.37500.37060.02338
0.01133
0.005755
0.977
0.034800.02172
0.01087
0.005470
0.991
0.06966
0.002982
0.001184
0.0007348
0.688
0.009901
0.0031210.001144
0.0004529
1.337

Table 3.3: Domain-decomposition vs. computation times
Decomposition
4x 4
4 x 2
2x 2
2x 1
1x 1

Nek

lVuk

72
144
288
576
1152

169
325
625
1225
2401

~Dk
49
91
169
325
625

79

tconv

32.28
32.25
32.17
32.16
32.14

Computation times(sec. )
415.4
510.9
720.1
1231.
1090.

)'1

4

8

16

n

32

4

8

16

n

32

0.1
- -0- - Iso-P2 PI / 4x4- -D- - Iso-P2 PI / 2x2
0

Error

- - -cr - - Modified Iso-P2 PI /4x4
- - -0- - - Modified Iso-P2 PI /2x2--O-PI / 4x4

"

Error

0
\

'
\ ,

b
"

0.1

9,
0.01

--O-PI / 2x2--+-- Non-DDM

,'

b, '

q
"

0.01

0.001

(a)

(b)

Figure 3.5: Errors in a test problem. (a) Velocity, (b) Pressure

80

4

8

16

32 n

0.1

Error

0.01

I

....

"

""

""

- - -0- - - Modified Iso-P2

PI / 4x4- - -0- - - Modified Iso-P2
PI / 2x2---O-PI /4x4

0.001--O-Pl / 2x2

0.0001Figure 3.6: Errors of the Lagrange multiplier in a test problem.

81

Figure 3.7: Velocity vec t ors of th lid-drin'll c8xity flo \\" p roblem . R e = -100 , on a uniform 24 x 24 x 2
triangular subdivision and a -± x -± domaill-d comp o. ition82

......

-

Figure 3.8: Pressure contour lines for the lid-driven cavity flow problem , R e = 400 , on a uniform
24 x 24 x 2 triangular subdivision and a 4 x 4 dOlmain-decomposi t ion

83

3.5

Conclusion

In this chapter we have discussed parallelization of the finite element method in massively parallel
computers of distributed-memory type. The parallelization is based on decomposition of a domain.
We have started from presenting the two data types , 'pre-assembly type' and 'restriction type, fordata on nodal points. They simplify the implementation and reduce the amount of message-passings
between sub domains in both of the data parallelization and the domain decomposition method. We candescribe concisely the procedures of conjugate gradient method for some systems of linear equations
by the idea.
In order to reduce message-passings further a mortar element method was studied. In the methoda Lagrange multiplier in terms of weak continuity between subdomains is introduced. We have applied
the mortar element method to the Stokes equations of primitive variable form. Its procedure can be
described concisely by the idea of two data types.
Finally, we have considered a domain decomposition algorithm of the finite element scheme forflow problems described in Section 2.4. In the scheme sub domain-wise finite element spaces by iso-P2
PI/PI elements are constructed and weak continuity of the velocity between subdomains are treatedby a Lagrange multiplier method. This domain decomposition algorithm has advantages such as:
1. :v1ost part of the computation is completely parallelized.

2. Each sub domain-wise problem is a consistent discretized pressure Poisson equation so that it isregular.
3. The size of a system of linear equations to be solved by the conjugate gradient method is smaller
than that of the original consistent discretized pressure Poisson equation.
For the discretization of the Lagrange multiplier, we compared three cases: the conventional i o-P2
PI element, a modified iso-P2 PI element having no freedoms at both edges of interfaces, and theconventional PI element. In every case, we checked numerically in a sample problem that the scheme
could produce solutions which converged to the exact solution at the optimal rates for the velocity andthe pressure. In the latter two cases we have also observed the convergence of the Lagrange multiplier.
Employing the conventional PI element, we have computed the lid-driven cavity flow problem. The
computation time becomes shorter when the number of processor increases. ~Ioreover, we have madea remark on simultaneous normalization for the pressure and the Lagrange multiplier.

84

-

General Conclusiol1.
We have presented the upwind finit e element approximation based on the choice of up- and downwind
points. It has a potential to approximate the convection term in third-order accuracy. Combining thisapproximation with the iso-P2 PI/PI elements and the procedure of a consistent discretized pressure
Poisson equation, the upwind finite element scherne for the two- and three-dimensional incompressible
Navier-Stokes equations at high Reynolds numbers was presented. About this scheme we have got thefollowing results throughout this thesis.
• A remark for the suitable way to decompose a tetrahedron into 8 sub-tetrahedra has been given.
That is useful to the calculation by the three-·dimensional iso-P2 PI/PI elements in the boundary
layers.
• vVe presented a blind search algorithm and heuristic search algorithms for upwind-point search
problems. In the search problem of U , the heuristic search algorithm shows good performance on
scalar (conventional) computers while the blind one is good on vector computers. The relation
between the point-search problem of U and a point-search problem in the (d - I )-dimensional
space is discussed. In order to avoid a risk of an infinite loop caused by finite-precision arithmetic ,a topology-oriented method is presented , where the Cramel's formula is used rather than theGauss 's elimination method to solve a small problem in the heuristic search algorithms.
• Shape of stability regions with respect to the two parameter of the scheme , a stabilizing parameter a and the time increment T , have been investigated numerically. We recognize that asuitable choice of a stabilizes the cherne: (i) In order to choose a larger T, we need to choose asmaller a (ii) A choice of too small a makes the computation unstable regardless of T, (iii ) The
lower bound of a increases as Reynolds number increases.
• Some numerical results was shown to demonstrate that the solver can produce numerical results
possessing good Reynolds-number dependency and capturing three-dimensionality of phenomena
well. For other results of this scheme, we refer to [57] for computation of two dimensional flows85past a circular cylinder and for observation of the drag coefficient of the cylinder at high Reynolds
numbers, [56] for computations of a lid-driven cavity flow problem at high Reynolds numbers
and for a discussion of contribution of rounding errors to flow simulation, [60] for simulation of
streaklines emerging from the surface of a circular cylinder, and [33] for an air flow computation
around a vehicle in a clean room.
• For parallel executions of the finite element Dlethod in massively parallel computers of distributedmemory type, two data types, 'pre-assembly type' and 'restriction type', for data on nodal points
have been presented. They simplify the irnplementation and reduce the amount of messagepassings between subdomains in both of the data parallelization and the domain decomposition
method. We can describe concisely the procedures of conjugate gradient method for some systems of linear equations by the idea .
• A method to apply the mortar element method to the Stokes equations of primitive variables
form has been presented. Its implementation can be described concisely by virtue of the idea ofthe two data type .
• We have considered a domain decomposition algorithm of the finite element scheme ba ed on
the choice of up- and downwind points for flow problems. In the scheme, subdomainwise finite
element spaces by the iso-P2 PI/PI elements are used for the velocity and the pressure, respectively. A Lagrange multiplier concerning the weak continuity of the velocity between subdomains
is introduced. This domain decomposition algorithm has advantages such as:

1. NIost part of the computation is completely parallelized.
2. The size of the system of linear equations to be solved by the conjugate gradient method
is smaller than that of the original consistent discretized pressure Poisson equation.
3. Each subdomain-wise problem is a consistent discretized pressure Poisson equation so thatthey are regular.
For the discretization of the Lagrange multiplier, we compare three cases: (i)the conventional
iso-P2 PI element, (ii)a modified iso-P2 PI element having no freedoms at both edges of interfaces and (iii)the conventional PI element. In every case, we checked numerically in a sample
problem that the scheme could produce solutions which converged to the exact solution at the
optimal rates for the velocity and the pre sure. In the latter two cases we have also ob erved
the convergence of the Lagrange multiplier. Employing the conventional PI element we have86computed the lid-driven cavity flow problenl. The comp'utation time becomes shorter when thenumber of processor increases. Moreover, we have made a remark on simultaneous normalization
for the pressure and the Lagrange multiplier.

87

Bibliography[1] Y. Achdou and Y. A. Kuznetsov, Algorithm for a non conforming domain decomposition method, Tech.
Rep. 296, Ecole Poly technique, 1994.

[2] Y. Achdou and O. Pironneau, A fast solver for l -avier-Stokes equations in the laminar regime using mortar
finite element and boundary element methods, SIAM. J. Numer. Anal., 32, 985-1016, 1995.

[3] S. K. Aliabadi and T. E. Tezduyar, Parallel fluid dynamics computations in aerospace applications, Int. J.
Num. Meth. Fluids, 21, 783-805, 1995.
[4] K. Baba and M. Tabata, On a conservative upwind finite element scheme for convective diffusion equations,
R. A. 1. R. 0., Analyse Numerique/Numerical Analysis, 15, 3-35, 1981.

[5] M. Behr and A. Johnson, Computation of incompressible flows with implicit finite element implementations
on the Connection Machine, Compo Meth. Appl. Mech. Engrg., 108,99-118, 1993.

[6] M. Bercovier and O. Pironneau, Error estimates for finite element method solution of the Stokes problem
in the primitive variable, Numer. Math., 33, 211-224, 1979.

[7] C. Bernardi, Y. Maday and A. Patera, A new nonconforming approach to domain decomposition: themortar element method, Nonlinear Partial Differential Equations and their Applications XI, H. Brezis and
J. L. Lions (eds.), Longman Scientific & Technical, Essex, UK, 13-51, 1994.

[8] :vI. O. Bristeau, R. Glowinski, B. Mantel, J. Peri.aux and O. Pironneau, A finite element approximationof I avier-Stokes equations for incompressible viscous fluids, Approximation Methods for Navier-Stokes
Problems, Lecture Notes in Mathematics 771, R. Rautmann (ed.), Springer, Berlin, 78-128, 1980.

[9] A. N. Brooks and T. J. R. Hughes, Stream upwind,lPetrov-Galerkin formulations for convection dominatedflows with particular emphasis on the incompressible ~ avier-Stokes equations, Compo Meth. Appl. Mech.
Engrg., 32, 199-259, 1982.
[10] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, -orth-Holland, Amsterdam, 1978.
[11] J. Donea, S. Giuliani and H. Laval, Finite element solution of the unsteady Navier-Stokes equations by afractional step method, Compo Meth. Appl. Mech. Engrg., 32, 53-73, 1982.
[12] M. Fortin and F. Thomasset, Mixed finite-element methods for incompressible flow problems, J. Compo
Phys., 31, 113-145 1979.
[13] 1. P. Franca, S. L. Frey and T. J. R. Hughes , Stabilized finite element methods: I. Application to theadvective-diffusive model, Compo Meth. Appl. Meek Engrg., 95, 253-276, 1992.
[14] L. P. Franca and S. L. Frey, Stabilized finite element methods: II. The incompressible Navier-Stokes
equations, Compo Meth. Appl. Mech. Engrg., 99, 2'09-233, 1992.
[15] S. Fujima, Implementation of mortar element method for flow problems in the primitive variables, to appear
in Int. J. Compo Fluid Dyn.
[16] S. Fujima and ?v!. Tabata, High performance implementation of the upwind finite element scheme basedon the choice of up- and downwind points, Computational Fluid Dynamics '94, Proc. Second European
Computational Fluid Dynamics Conference, S. v\ agner, E.H. Hirschel, J. Periaux and R. Piva (eds.),
Wiley, Chichester, UK, 129-134, 1994.
[17] S. Fujima, :V!. Tabata and Y. Fukasawa, Extension to three-dimensional problems of the upwind finite
el ment scheme based on the choice of up- and downwind points, Compo Meth. Appl. Mech. Engrg., 112,109-131 1994.

88

[18] Y. Fukuyama, Flows in a three-dimensional cubic cavity, Proc. 3rd Workshop on Computational Fluid
Dynamics, Japan Soc. Mechanical Engineers, 91-98 (in Japanese) , 1990.
[19] T. Fuseki, Flows in a three-dimensional cubic cavi ty, Proc. 3rd Workshop on Computational Fluid Dynamics, Japan Soc. Mechanical Engineers, 85-90 (in Japanese), 1990.
[20] V. Girault and P. A. Raviart , Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms,
Springer, Berlin, 1986.
[21] R. Glowinski, T.-W. Pan and J. Periaux, A one shot domain decomposition/fictitious domain method forthe ~avier-Stokes equations, Domain Decomposition Methods in Scientific and Engineering Computing,Proc. 7th Int. Conf. on Domain Decomposition, D. E. Keyes and J. Xu (eds.), Contemporary Mathematics
180, A. :vI. S., Providence, Rhode Island, 211-220, 1994.
[22] P. M. Gresho, S. T. Chan, R. L. Lee and C. D. Upson, A modified finite element method for solving thetime-dependent, incompressible Navier-Stokes equations, Part 1: Theory, Int. J. Num. Meth. Fluids, 4,557-598,1984.
[23] F. R. Hama, Three-dimensional vortex pattern behind a circular cylinder, J. Aeronaut. Sci., 24, 156-159,
1957.
[24] M. O. L. Hansen, J. N. Sorenson and V. A. Barker, A numerical investigation of 3-D flow past an infinite cylinder, Computational Fluid Dynamics '9~~, Proc. First European Computational Fluid Dynamics
Conference, C. Hirsch, J. Periaux and W. Kordulla (eds.), Elsevier, Amsterdam, 375-381, 1992.
[25] J. C. Heinrich, P. S. Huyakorn, O. C. Zienkiewicz and A. R. litchell, An 'upwind' finite element scheme
for two-dimensional convective transport equation, Int. J. Num. Meth. Engrg., 11, 131-143, 1977.
[26] T. J. R. Hughes and A. Brooks, A theoretical framework for Petrov-Galerkin methods with discontinuous weighting functions: Application to the streamline-upwind procedure, Finite Elements in Fluids 4
R. H. Gallagher et al. (eds.), Jown Wiley & Sons, ew-York, 47-65, 1982.
[27] T. Ikeda, Maximum Principle in Finite Element: Models for Convection-Diffusion Phenomena, Lecture
Notes in Numerical and Applied Analysis 4, North-Holland, Tokyo, 1983.
[28] R. Iwatsu, K. Ishii, T. Kawamura, K. Kuwahara and J. M. Hyun, Simulation of transition to turbulence ina cubic cavity, AIAA-Paper 89-0400, Proc. AIAA 27th Aerospace Sciences Meeting, Reno, ~ evada, 1989.
[29] C. B. Jiang, M. Kawahara, K. Hatanaka and K. Kashiyama A three-step finite element method forconvection dominated incompressible flow, Compo Fluid. Dyn. J., 1,443-462, 1983.
[30] B. Jiang, T. Lin and L. A. Povinelli, Least squares finite element solution of 3D incompressible NavierStokes problems, Computational Fluid Dynamics 'g2, Proc. First European Computational Fluid Dynamics
Conference, C. Hirsch, J. Periaux and "V. Kordulla (eds.), Elsevier, Amsterdam, 391-395, 1992.
[31] K. Kakuta and ~ . Tosaka, Exponential function type Petrov-Galerkin finite element method for unsteadyincompressible viscous flow problems, Research reports in structure field in Japan Soc. Construction 439,
189-198(in Japanese) 1992.
[32] H. Kanayama, Discrete models for salinity distribution in a bay: Conservation laws and maximum principle,
Theoretical and Applied Mechanics 28, "C niv. Tok) 0 Press, 559-579, 1980.
[33] H. Kanayama, K. Toshigami Y. Tashiro :V!. Tabata and S. Fujima, Air flow computation around an
automatic guided vehicle, Int. J. Compo Fluid. Dyn., 7, 155-162, 1996.
[34] Y. Kato and T. Tanahashi, Flows in a three-dimensional cavity, Proc. 3rd Workshop on Computational
Fluid Dynamics, Japan Soc. :vIechanical Engineers, 99-106 (in Japanese), 1990.
[35] T. Kawamura, H. Takami and K. Kuwahara, ~ew higher-order upwind scheme for incompressible. -avierStokes equations, Numerical Methods in Fluid Dynamics, Soubbaramayer and J. P. Boujot (eds.) Lecture
Notes in Physics 218, Springer, Berlin, 291-295, 1985.
[36] J. G. Kennedy, :V!. Behr, V. Kalro and T. E. Tezduyar, Implementation of implicit finite element methods
for incompressible flows on the C. 1-5, Compo Meth. Appl. Mech. Engrg., 119,95-111 1994.
[37]~. Kondo ~. Tosaka and T. ~ishimura, ~umericaJ simulation of viscous flows by the third-order upwind
finite element method, Theoretical and Applied Mechanics 39 , "Cniv. Tokyo Press, 237-250, 1990.

9

[38] H. C . Ku , R. S. Hirsh and T. D. Taylor , A pseudospectral method for solution of the three-dimensional
incompressible Navier-Stokes equations , J. Comp o Phys ., 70 , 439-462 , 1987.
[39] C. L . Lawson , Software for Cl surface interpolation , Mathematical Software III , J. Rice (ed .), Academic
Press, New-York, 161-194, 1977.
[40] P. L. Lions, On the Schwarz alternating method. I , First Int ernational Symposium on Domain D ecomposition Methods for Partial Differential Equations , R. Glowinski, G . H. Golub , G. A. Meurant and J. Periaux(eds.) , S. 1. A. M. , Philadelphia, pages 1-42 , 1988.
[41] J. A. :\IIeijerink and H. A. van der Vorst , An iterative solution method for linear systems of which thecoefficient matrix is a symmetric M-matrix , Math. Comp. , 31 , 148-162 , 1977.
[42] K. Mori and Y. Doi , Three-dimensional flow simulation around a two-dimensional circular cylinder Proc.
4th Symp. on Computational Fluid Dynamics, Chuo-University, Tokyo, 503-506(in Japanese) , 1990.
[43] S. ~agoya, Stability Analysis and Error Estimate for Difference Problems Derived from a Spatially On e
Dimensional Convective Diffusion Equation, Thesis , The University of Electro-Communications , Tokyo ,1996.
[44] K. Ohmori and T. Ushijima, A technique of upstream type applied to a linear nonconforming finite element
approximation of convective diffusion equations , R.A.I.R. O. Analyse Numerique/Numerical Analysis, 18 ,309-332, 1984.
[45] C.- Y. Perng and R. L. Street , Three-dimensional unsteady flow simulations: alternative strategies for a
volume-averaged calculation , Int. J. Num. Meth. Fluids , 9 , 341-362 , 1989.

[46] A. Pothen, H. D. Simon and K.-P. Liou , Partitioning sparse matrices with eigenvectors of graphs , SIAM
J. Matrix Anal. Appl., 11 , 430-452 , 1990.
[47] A. Quarteroni, Domain decomposition and parallel processing for the numerical solution of partial differential equations Surv. Math. Ind. , 1 , 75-118 , 1991.
[48] R. D. Richtmyer and K. W. :Nlorton , Difference Methods for Initial Value Problems, Interscience, ~ew- York,1967.
[49] P. J. Roache , Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, 1976.
[50] K. Sugihara, Design of geometric algorithms in finite-precision arithmetic , Bulletin of the Japan S. 1. A.
M., 4, 280-299 , 1991.
[51] A. Suzuki , Implementation of domain decomposition methods on parallel computer ADEN ART, Parallel
Computational Fluid Dynamics: New Algorithms and Applications, 1 -. Satofuka, J. Periaux and A. Ecer(eds.) , Elsevier , 1995.
[52]

1. Tabata, A finite element approximation corresponding to the upwind finite differencing, Mem. Numer.
Math., 4 47-63, 1977.

[53] TvI. Tabata, A theoretical and computational study of upwind-type finite element methods , Patterns andWav es, T. l -ishida et al. (eds .). Kinokuniya/Korth-Holland , Tokyo , 319-356 , 1986.
[54] :\11. Tabata, Finite element method , Encyclopaedia of Math ematical Science , H. Hironaka et al. (eds. ), (as asubsection 'Partial Differential Equation (2) ' in Section ' ~umerical Computation') , Osaka-Shoseki , 0 aka,
930-937 (in Japanese). 1991.
[55] :\II. Tabata and S. Fujima, A finite element scheme for the ~avier-Stokes equations using a third-order
upwind approximation , Finite Element Analysis in Fluids , Proc. 7th Int. Conf. on Finite Element Methods
in Flow Problems T. J. Chung and C. R. Karr (eds.) , "Cniv. Alabama in Huntsville Press , Huntsville1034-1039 , 1989.
[56] YI. Tabata and S. Fujima, An upwind finite element scheme for high-Reynolds-number flows Int . J. Num .
Meth. Fluids, 12 , 305-322. 1991.
[57]

I. Tabata and S. Fujima Finit -element-analysi of high Reynolds number flows past a circular cylinder,
J . Comp o Appl. Math. , 38 -111-424 , 1991.

[58] :\II. Tabata and S. Fujima, Some applications of the third order upwind finite element approximation ,Computational Fluid Dynamics '92, Proc . First European Computational Fluid Dynamics Conference,
C. Hirsch J. Periaux and vV. Kordulla (eds.) , Elsevier Amsterdam , 139-145 1992.

[59] M. Tabata and S. Fujima, Numerical analysis of flow problems by the finite element method, Mathematics,
48, 22-36 (in Japanese), 1996.
[60] M. Tabata, S. Fujima and K. Itakura, Numerical simulation of streaklines emerging from the surface of acircular cylinder, PIOC. First China-Japan Seminar on Numerical Mathematics, Z.-C. Shi and T. Ushijima(eds.), World Scientific, Singapore, 146-154, 1993.
[61] T. Tamura and K. Kuwahara, Direct finite difference computation of turbulent flow past a circular cylinder,
PIOC. Int. Symp. on Computational Fluid Dynami'cs (at Nagoya), 701-706, 1989.
[62] T. Tamura, 1. Ohta and K. Kuwahara, On the reliability of two-dimensional simulation for unsteady flows
around a cylinder-type structure, J. Wind Engrg Ind, Aerodynamics 35, 275-298, 1990.
[63] T. Tanahashi, H. Okanaga and T. Saito, GSMAC finite element method for unsteady incompressible
Navier-Stokes equations at high Reynolds numbers, Int. 1. Num. Meth. Fluids, 11,479-499,1990.
[64] T. E. Tezduyar, S. Mittal, S. E. Ray and R. Shih, Incompressible flow computations with stabilized bilinear
and linear equal-order-interpolation velocity-pressure elements, Compo Meth. Appl. Mech. Engrg., 95, 221242, 1992.
[65] F. Thomasset, Implementation of Finite Element Methods for Navier-Stokes Equations, Springer, ~ew
York, 1981.
[66] J. F. Thompson, Z. U. A. Warsi and C. W. Mastin, Numerical Grid Generation, Foundations and Applications, North-Holland, Amsterdam, 1985.

91

List of Tables
2.1
2.2
2.3

3.1
3.2

3.3

Comparison of the heuristic- and the blind-search algorithm on a scalar- and a vector-computer.
Parameters of the scheme and convergence times
Deviations from x3-mean values . . .

29
42

Iteration times of the inner CG loop
Errors in a test problem . . . . . . .
Domain-decomposition vs . computation times.

71
79
79

9i.,'""

53

List of Figures
1.1
1.2
1.3
1.4

Iso-P2 P1/P1 elements, d = 3 . . . . . . . . . . . .
Barycentric domain . . . . . . . . . . . . . . . . . .
Example of a decomposition into 8 sub-tetrahedra
Example of a decomposition into 8 sub-tetrahedra

2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8

Finite element approximation based on the choice of upwind element
Velocity u(O) and points W, U, 0, D, B , d = 2 .
Velocity u(O) and points W, U, 0, D, B, d = 3 . . . . . . . . . . . .
Heuristic search of point U . . . . . . . . . . . . . . . . . . . . . . .
Example of an infinite loop in execution of the Lawson's algorithm, d = 3
Relation between the search problem of U and a point search problem, d = 2
Stability regions for Re = 100, 400 and 1000 in a test problem . . . . . . . . . '.
Two-dimensional flow problem around a circular cylinder, boundary conditions(u = 0 on thecylinder surface) and a division of the domain. . . . . . . . . . . . . . . . . .
2.9 Stability regions, Re = 3200, in the flow problem around a circular cylinder .
2.10 Stability regions, Re = 10000, in the flow problem around a circular cylinder
2.11 Finite element subdivision of a half cube. . . . . . . . . . . . . . . . . . . . .
2.12 Velocity vectors, Re = 100, on the planes (a)x2 = 0.5, (b)Xl = 0.5 and (C)X3 = 0.5
2.13 Velocity vectors, Re = 400, on the planes (a)x2 = 0.5, (b)Xl = 0.5 and (C)X3 = 0.5
2.14 Velocity vectors, Re = 1000, on the planes (a)J~2 = 0.5, (b)Xl = 0.5 and (C)X3 = 0.5 .
2.15 Velocity vectors, Re = 2000, t = 50, on the planes (a)x2 = 0.5, (b)Xl = 0.5 and (C)X3 = 0.5
2.16 Paths of marker particles, (a)Re = 100, (b)Re = 400, (c)Re = 1000
2.17 Graphs of Ul (0.5, 0.5,,) and U3(" 0.5, 0.5), Re == 400 . . . . . . . .
2.18 Graphs of Ul (0 .5 ,0.5 , .) and U3 (·,0 .5 ,0.5), Re == 1000 . . . . . . . .
2.19 Three-dimensional flow problem around a circular cylinder, boundary conditions (u = 0 on thecylinder surface) and a division of the domain. . . . . . . . . . . . . . . . . . . . . . . . .
2.20 Velocity on the plane X3 = 0.00. Vectors of (Ul' U2) and contours of U3 = ±0.02, ±0.04,. . .
2.21 Velocity on the plane X3 = 0.25. Vectors of (Ul, U2) and contours of U3 = ±0.02, ±0.04, . . .
2.22 Velocity on the plane X3 = 0.50. Vectors of (Ul U2) and contours of U3 = ±0.02, ±0.04,...
2.23 Velocity on the plane X3 = 0.75. Vectors of (Ul' U2) and contours of U3 = ±0.02, ±0.04, . . .
2.24 Velocity on the X2 = 0 plane: vectors of (Ul' U3) and contours of U2 = 0, ±O.l, ±0.2, . . .
2.25 Growth of deviations from X3-mean values
2.26 Finite element subdivision, m = 6. . . . . . . . . . . . .
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8

:YIassively parallel computer of distributed memory typeThe bases of the mortar function space. . . . . . . . . .
Shapes of the solution at each step in the outer CG loop (u = 1.0 in the interval where they are
not plotted) . . . . . . . . . . . . . . . . . . . . . . . . . . .
Domain-decomposition(4 x 4) and the triangulation(n = 32)
Errors in a test problem. (a) Velocity, (b) Pressure. . . . .
Errors of the Lagrange multiplier in a test problem.
Velocity vectors of the lid-driven cavity flow problem, Re = 400, on a uniform 24 x 24 x 2
triangular subdivision and a 4 x 4 domain-decomposition . . . . . . . . . . . . . . . . . . . . . .
Pressure contour lines for the lid-driven cavity How problem Re = 400, on a uniform 24 x 24 x 2
triangular subdivision and a 4 x 4 domain-decomposition . . . . . . . . . . . . . . . . . . . . . .

15
16
18
19
23
24
24
29
31
32
36
37
38
38
43
44
45
46
47
48
49
50
51
53
54
54
55
55
56
56
59
68
71
77
80
81
82
83

--

(1)

Masahisa TABATA and Shoichi FUnMA (83Yfffi1£j\, Airl3~ ~-),Finite-element analysis of high-Reynolds nUlmber flows past a circular cylinder(pjif.f ~ j@1 ~t,~ ~ v -1 / )v ;Z:'i5£VTEh0)1r~R~~M1JT),

Journal of Computational and Applied Mathlematics, Vol. 38, pp. 411-424, 1991.

(2)

Shoichi FUnMA, Masahisa TABATA and Yasuji FUKASAWA(Ri rl3~ ~ - , 83 Yfffi 1£)\,?* i~ ~ PJ ) ,Extension to three-dimensional problems of the upwind finite element scheme based on
the choice of up- and downwind points(J~LtJOO. T/8~1R~:£·':5 < JOO.J:.~1f~R~~ A ~ - A 0) 3 *jGFp~M~O)tJt~),

Computer Methods in Applied Mechanics and Engineering, Vol.l12, pp.109-131,
1994.

(3)

Shoichi FUnMA (Ri rl3~ ~ -),Implementation of a mortar element method for flow problems ill the primitive
variables(mr ~B~~~Jf~0)1JTEhFp~~~0)-t)v ~ )v~~1:t0)~~ ),To appear in the International Journal of Conlputational Fluid Dynamics.

(1)Masahisa T ABA TA and Shoichi FUnMA ( EB YffiI lE:A., ~ F£3, J§ ~ ) ,An upwind finite element scheme for high-Reynolds-number flows(~ v1 / )v ;(it1JTEhO)t-=~0)~~~ ~I:.~:ff~N~~A~-A),

International Journal for Numerical Method.s in Fluids, Vol. 12, pp. 305-322, 1991.

(2)

Masahisa TABATA, Shoichi FUnMA and Kazuhiro ITAKURA( E81iffiLiE:A., ni F£3, J§ ~ , t& ~*D 1E) ,Numerical simulation of streaklines emerging from the surface of a circular cylinder inhigh-Reynolds-number flows(~ v 1 /)v ;('itifrEh91: FI3 tf£ZITff 7J)~) 1NtB T ~ 1JTElPlFff'7'F.O) it1i[ Y ~ ::L V - Y :3/ ) ,Proceedings of the First China-Japan Seminar on Numerical Mathematics, pp.146-154,
1993.

(3)

Shoichi FUnMA and Masahisa TABATA (~ F£3' J§~, EBYffiIlE:A.) ,High performance implementation of the upwind finite element scheme based on the
choice of up- and downwind points( ~ J:. ~ T ,8~tR~ = £ 0<~ J:.~:ff ~N ~~~ A ~ - A O)~ tE~~ ~~) ,Computational Fluid Dynamics '94, Proceedings of the 2nd European Computational
Fluid Dynamics Conference, pp.129-134, 1994.

(4)

Hiroshi KANAYAMA, Kouji TOSillGAMI, Yukari TASHIRO, Masahisa TABATA
and Shoichi FUnMA (i: LlJ Jl, 5tf J:.1E P]l'

EB 1-t ~ 7J) ~ , EB YffiI lE:A., ni F£3, J§ ~ ) ,Air flow computation around an automatic guided vehicle(§ JJJ~ *~mJ illO) ~VfEBt $),International Journal of Computational Fluid Dynamics, Vol. 7, pp.155-162, 1996.