Spectral Element Methods: Algorithms and Architectures 

/fl/l J' 3/ 

+fi ~I''fl1 

/flJ /k/-C& 

Paul Fischer, Einar M. Rdnquist 

Daniel Dewey, Anthony T. Patera 

Department of Mechanical Engineering 

Massachusetts Institute of Technology 

Cambridge, MA 02139 

/ 

U.S.A. 

Running title: Spectral Element Algorithms and Architectures 

(liA5a-CR- 182701) SPBClIiAL ELEUEII UETIIODS: M 8-2ca78 

ALGOR IT H ES A NL C 6 L ElP E C I U E E 5 4 Rass ach use t ts 

Itst. cf lech,) 33 F CSCL 09B 

Unclas 

~3/61 0134769 
Abrtract. Spectral element methods are high-order weighted 

residual techniques for partial differential equations that 

combine the geometric flexibility of finite element methods 

with the rapid convergence of spectral techniques. In this 

paper we describe spectral element methods for the 

simulation of incompressible fluid flows, with special 

emphasis on implementation of spectral element techniques on 

medium-grained parallel processors. Two parallel 

architectures are considered: the first, a connnercially 

available message-passing hypercube system; the second, a 

developmental reconfigurable architecture based on GeometryDefining Processors. High parallel efficiency is obtained 

in hypercube spectral element computations, indicating that 

load balancing and conmumication issues can be successfully 

addressed by a high-order technique/medium-grained processor 

algorithm-architecture coupling. 

1. Introduction. Spectral element methods are high-order 

weighted-residual techniques for partial differential 

equations th’at exploit both the common foundations and 

competitive advantages of h-type finite element methods 

(Strang and Fix, 1973) and p-type spectral techniques 

(Gottlieb and Orszag, 1977). In the spectral element 

discretization (Patera, 1984; Maday and Patera, 19871, the 

computational domain is broken up into macro-spectral 

elements, and the dependent and independent variables are 

approximated by Nth order t e n s o r - p r o d u c t p o l y n o m i a l 

expansions within the individual subdomains. Variational 

projection operators and Gauss numerical quadrature are used 

to generate the discrete equations, which are then solved by 

*Department of M e c h a n i c a l E n g i n e e r i n g , M a s s a c h u s e t t s 

Institute of T e c h n o l o g y , C a m b r i d g e , MA 02139 USA. 

1 
w 2 

direct or iterative procedures based on tensor-product sumfactorization techniques. Convergence to the exact solution 

is achieved by increasing the degree, N, of the polynomial 

approximation, while keeping fixed the number and identity 

of the underlying spectral elements. 

domain decomposftion in that an individual spectral element 

is the ”basic unit” as regards strong coupling between 

degrees-of-freedom. In the spectral element discretization 

the elements serve as the basic building blocks for mesh 

generation and physics, as well as the basic units of 

locally structured mesh. The latter is critical as regards 

the computational efficiency of the method, as it allows for 

the sum-factorized tensor-product matrix-vector products 

that render h i g h - o r d e r m e t h o d s competitive with low-order 

techniques. 

Although spectral element methods have proven 

competitive with h-type finite element methods for 

incompressible flow simulations on serial and vector 

computers (Ghaddar, Magen, Mikic and Patera, 19861, it is 

clear that in the future techniques must be judged on the 

basis of their performance on parallel machines. 

F o r t u n a t e l y , w h e r e a s global spectral methods are not 

obviously parallelizable due to the strong coupling of all 

points in the computational domain, spectral element methods 

have an intrinsic ” d o m a i n - d e c o m p o s i t i o n ” g r a n u l a r i t y that in 

turn leads to natural and efficient implementation on 

medium-grained parallel processors. 

of spectral element algorithms. In Section 2 we describe 

the basic spectral element discretization as applied to 

elliptic problems. In Section 3 we discuss the parallel 

implementation of spectral element elliptic algorithms on a 

message-passing hypercube system, and present results for 

parallel efficiency on a commercially available Intel 

machine. The viability of the h i g h - o r d e r m e t h o d / m e d i u m - 

grained approach to parallel computing is discussed, and the 

issues of communication latency and load-balancing are 

addressed. Lastly, in Section 4 we present a developmental 

special-purpose reconfigurable architecture, GeometryDefining Processors, and discuss the advantages of parallel 

implementation of spectral element and substructured finite 

element algorithms on this novel configuration. 

The spectral element discretization is an example of 

In this paper we focus on the parallel implementation 

2. Spectral E l e m e n t M e t h o d a . 

Ellintic Equatiopg. To illustrate the basic spectral 

element concepts we first consider the following onedimensional elliptic Helmholtz problem: Find u(x) defined 

over A=]-1,1[ such that 

xEA , 2 (2.1a) -(pu’)’ + x u = f 

(2.lb) u(-l) = u(l) = 0 
Y 3 

Here prime denotes differentiation, and we assume that p(x) 

’ Po solving (211) is based on the variational formulation: 

Find u E Ho(d) such that 

(2.2) a(u,v> = (f,v) V~EH&O , 

where 

> 0, and E R. The spectral element discretization for 

Here Hi is the space of all functions which satisfy the 

homogeneous boundary conditions (2.lb). which are square 

integrable, and whose derivatives are also square 

integrable. 

The spectral element discretization proceeds by 

breaking up the interval A into K subintervalgh(spectral 

elements) dl,..,%, where the length of the k interval is 

Lk. The space of approximation for the solution y is then 

takef to be a high-order polynomial subspace of H (A), 

Xh-Ho(A)nP .(A), where PN K(A)=UkP,(Ak), and P (9,) is the 

space of arl polynomials of degree d4 on the inferval nk. 

To accurately evaluate the integrals associated with the 

bilinear forms in (2.3) we use Gauss-Lobatto Legendre 

numerical quadrature with quadrature points and weights 

given by (ti,pi), i-0, ...., N. 

The abstract formulation of the resulting discrete 

problem is: Find uh E Xh such that 

where subscript h,GL refers to Gauss-Lobatto numerical 

quadrature of the bilinear forms in (2.3). Convergence of 

the spectral element solution uh to the exact solution u is 

then achieved by increasing the degree of the approximation, 

N, for fired K. It is shown in Maday and Patera (1987) for 

this spectral (p-type) convergence strategy that the error 

goes to zero exponentially fast with N for sufficiently 

smooth solutions u. This should be contrasted to finite 

element h - t y p e m e t h o d s (N fixed, K=+oo) for which algebraic 

convergence results. 

with a nodal basis, that is, any function w in X is 

written in terms of elemental Lagrangian inperpolants hi(r) 

through the Gauss-Lobatto Legendre points, 

To proceed further, we represent our polynomial space 

h 

(2.5a) wh(r) k = Nk wihi(r) 

i -0 

xEAk =+ rEA , 

2 ( 2.5 b ) h iEPN(d) , h ( ) =b ij Vi,j€{O, ...., N} , 
4 

where wk is the value of wh at theOlocal .point ti in the 

interval Ak. 

the boundary conditions we further require that 

In order to ensure C -continuity and satisfy 

(2.6a) k k+l VkE(1,. . . . ,K-1} , "N = wO 

and 

(2.6b) w 1 'WN30 K . 0 

Expressing u and v in terms of the nodal basis (2.61, and 

choosing each test !unction to be nonzero at only one global 

collocation point, we arrive at the following set of 

discrete equations: 

k 2Lk K,2N N 

E [---E E P pkD .D .u + -p.ui] 1 - 

k-1 Lkj=O q-0 q q 41 qJ j 2 

(2.7) 

K, Lk k E - Pifi P 

k-1 2 11 

where pk and fk are the values of p(x) and f(x) at the local 

point rgc in , Di -dh.(Ei)/dr, and denotes the direct 

stiffness's-a$ion drhicd incorporates (2.6). 

with p(x1-e , f(x)=e (cosx-sinx), x=O on x=]O,T[, for which 

the solution is u(x)=-sinx. The particular spectral element 

discretization corresponds to division of the domain into 

K-2 similar spectral elements, each of order N. Fig. 1 shows 

a plot of the maximum nodal error llu-uhlIGkfLm as a 

function of the total number of degrees-of- eedom N =2N+1. 

Exponential convergence is achieved due to the fact that the 

solution is analytic. 

The spectral element discretization of multidimensional elliptic equations corresponds to a tensorproduct extension of the one-dimensional method. We 

consider he,re the two-dimensional Poisson equation on a 

domain nwith homogeneous boundary conditions on the domain 

boundary a?, 

(2.8a) -Au = f in f? 

(2.8b) u = o on LU?. 

As a numerical example we consider the problem (2.1) X X 

In this case the domain is broken up into K quadrilateral 

elements f? l,...,nK, and the corresponding discrete 

formulation is given by: Find uh E Xh such that 

1 where Xh-Ho(n)nPN .<n>. Here PN .(n)=U P (n 1, and P (n 

the space of all fiiecewise polybomialskoy dkgree IN Citk is 
respect to each spatial variable x and y. For an error 

analysis of the discretization (2.9) we refer to Funaro 

(1986) where it is shown that.spectra1 convergence obtains 

as N=+W for fixed K. 

A basis for the space Xh is constructed by expressing 

any function w 

element nk ( (k,y> E nk=+ (r,s) E AXA 1, 

in Xh in tensor-product form within each 

NN 

i-0 j=O 

wk(r,s) k = wk hi(r)h.(s) 

ij J (2.10) 9 

where h. are the one-dimensional Gauss-Lobatto Lagrangian 

interpofants defined in (2.5b). The solution u , the 

testfunctions vh, and the (isoparametric) geomeiry (x 

are all expressed in terms of the nodal basis (2.10>,hink 

the integrals in (2.9) are evaluated using Gauss-Lobatto 

Legendro X Gauss-Lobatto Legendre quadrature in (r,s). 

Choosing appropriate "delta" testfunctions v we arrive at 

the set of linear equations, 

y ) 

h 

(2.11) Au-Bf , 

where A is the discrete Laplace operator, B is the 

(diagonal) m a s s m a t r i x , and underscore denotes nodal 

unknowns. 

(2.8a) on the domain ~=~xE]O,l[,yE]0,1+~sin~x[~ with f-0 and 

u=sinx.~-~. Fig. 2 shows the isoparametric spectral element 

discretization (K-41, while in Figure 3 we plot the maximum 

nodal error Ilu-u II 00 as a function of the number of 

degrees-of-freedom in dne spatial direction. Note that 

although the domain is relatively deformed compared to the 

mapped rectilinear problem, exponential convergence to the 

analytic solution is obtained (Ronquist and Patera, 1987). 

As a two-dimensional example w e c o ~ s i d e r the problem 

h GL L 

tic Solvers. In this section we address the problem of 

solution of the linear positive-definite symetric systems 

resulting from spectral element discretization of selfadjoint elliptic equations. For large three-dimensional 

problems the most attractive approach is iterative solvers, 

both in terms of operation count and storage requirements. 

In order for a n y s o l u t i o n m e t h o d to be efficient in this 

context, fast matrix-vector evaluation schemes are needed; 

for high-order methods this is effected by tensor-product 

sum factorization, which we now discuss. 

Consider the solution of the two-dimensional Poisson 

equation (2.8) in a rectilinear domain n, resulting in the 

spectral element discretization (2.11). The matrix-vector 

product A& can be written in terms of an elemental sum 

2 Vi,j€{O, ...., N} . k K, N 

k-1 m,n 

(2.12) A&= *ijmn kU mn 
6 

4 Naive evaluation of (2.12) reqtires O(N 

element,kand also requires O(N storage as the elemental 

matrix A is, in general, full. Fortunately, due to the 

tensor-pfJWct representation (2.101, a typical elemental 

term in (2.12) can be f5ctoZed as (consider here only the 

term corresponding to a /ax 

operations per 

k 2 N N 

(2.13) ED .[ P P.[ CDqmumj 111 Vi,j€{O,.. .,N} , q-0 q1 J m-0 

5 yielding an operation c2unt of O(N 

requirement of only O(N 1. 

first forming the elemental matrices Ak knd then evaluating 

the elemental matrix-vector produc s Ai results in an 

operation count and storage of O(N per element, while the 

tensor-product sum factorization yields an operation count 

of O(N and storage requirement of O(N 1. It should be 

noted that for a ge3fial d-dimensional problem, the 

operation count O(N persists even for isoparametric 

discretization of non-separable problems. 

Given the efficient tensor-product evaluation, the 

system (2.11) is solved by preconditioned conjugate gradient 

iteration (Golub and Van Loan, 1983). First, we construct 

A, an h-type finite element approximation defined on the 

spectral element quadrature points (Orszag, 1980; Deville 

and Mund, 1985). The matrix A is then further approximated 

by its incomplete Cholesky decomposition, A , which in turn 

is used as a preconditioner for A (Meijerink and Van der 

Vorst, 1977; Kershaw, 1978). As a test problem to 

demonstrate the preconditioned conjugate gradient algorithm, 

we consider solut'on of th two-dimensional Poisson equation 

(2.8a), with f-e 

(xE~O,l[,yE]O,1[). The domain is discretized by an array of 

K-K spectral elements, each of order N. In Figure 4 we 

plot the convergence history as a function of number of 

iterations for K -4 and N-5,7,9, and 11. The results 

demonstrate the significant savings due to effective 

preconditioning. 

and a storage 

For three-dimensional problems thg naive approach of 

6 

4 3 

, ~-~e "', and n defined by 

1 

Stokes Discretizatiom. To illustrate how the spectral 

element discretization of the Poisson equation extends to 

the full Stokes equations we consider here the following 

problem: Find a velocity u=(u,v) and a pressure p in the 

domain n=]-l,l[x]0,27r[ such that 

(2.14a) -VAU + Vp = f , 

(2.14b) -divu = 0 , 

subject to the semi-periodic boundary conditions 

(2.15a) VyE]0,27r[, u(-l,y) = U(1,Y) = 0 9 

(2.15b) VxE]-1,1[, u(x,o) = u(x,27r) 
7 

Here V is the kinematic viscosity of the fluid, and f is the 

prescribed force. 

Due to periodicity the dependent variables can be 

written in terms of Fourier series in the y-direction, 

where i = fl. The Fourier representation (2.16) results in 

a set of decoupled equations for each Fourier mode n, the 

variational statement of which is given by (dropping the 

Fourier superscfipt ijnd carat notaiion): For n E N* find 

u-(u;v) in X=[Ho(n>l and p in M-L (A) such that 

(2.17b) (q,uX+inv)* = 0 ,Vq€M , 

where 4=]-1,1[, L2 is the space of square-integrable 

functions, * refers to complex conjugate, and ($,$'I = In 

d$'*dx . For simplicity, we do not consider the case of n-0. 

In the same fashion as for the elliptic problem (2.11, 

the spectral element discretization proceeds by breaking up 

the interval A into K subintervals and searching for a 

solution in polynomial subspaces of X and M. However, 

choosing the polynomial degree N to be the same for the 

velocity and the pressure leads to spurious modes in the 

pressure. It can be shown that an optimal strategy is to 

use a staggered mesh for the pressure for which the 

associated polynomial degree is N-2 (Bernardi and Maday, 

1987; Maday, Patera and Rdnquist, 1987a; Maday, Patera and 

Rdnquist, 

(2.17) can then be stated as: Find uN in XN-[Ho(n)nPN,K(A)l 

and PN in %='N-2,K (A) such that 

1987b). The discrete formulation corfesponding to2 

where the subscripts N,G and N,GL denote numerical 

quadrature based on the Gauss and Gauss-Lobatto points, 

respectively. Theoretical error estimates for the 

approximation error due to the discretization (2.18) can be 

found in (Maday and Patera, 1987; Maday, Patera and 
8 

Rdnquist, 1987a), in which spectral convergence is 

demonstrated as N=+m for fixed K. 

element nk the velocity uN is expanded in terms of N 

Lagrangian interpolants h through the Gauss-Lobatto i Legendre points i$ i-0, ...., N , 

We now define the bases for the space XN%. Int&ach 

order 

i’ 

(2.19) uN(r) k = Nk uihi(r) 

i-0 

xEAk =+ rEA , 

while the pressure is exkanded in terms of (N-2lth order 

Lagrangian interpolants h through the Gauss Legendre points 

Si, i=l, ..., N-1, i 

Note the Gauss points are naturally suited for the pressure, 

which need not be continuous across elemental boundaries. 

The expansions for the the velocity and the pressure are 

inserted into (2.181, appropriate ”delta” testfunctions 

wN E XN and q E % are chosen, and we arrive at the 

discrete saddye problem, 

(2.21a) AX - Rp-Bf * 

(2.21b) -Qx*Q 

T 

Here the matrix A is the discrete Laplacian, B is the mass 

matrix, R is the discrete gradient o p e r a t o r , a n d underscore 

refers to nodal unknowns. 

problemwith V=1.0, n=l, 

As a numerical example we consider the following test 

on a domain A which is divided into K-2 similar spectral 

elements. We show in Figure 5 the error in the velocity and 

the pressure as the order N of the polynomial expansions is 

increased. As expected, exponential convergence is achieved 

due to the fact that the solution is analytic. In the case 

when the solution is less regular we obtain an algebraic 

convergence rate reflecting all the regularity of the 

solution (Maday, Patera and Rdnquist, 1987a). Note that we 

obtain spectral convergence with only limited continuity 

between elements due to proper choice of the weak form. 

Wes Solvers. In this section we consider solution of the 

algebraic system of equations (2.21) resulting from the 

spectral element discretization of the Stokes problem 

(2.14). It should be noted that the algorithms presented 

here are equally appropriate for other types of variational 

discretizations, and are in fact extensions of the classical 
9 

Uzawa algorithm used in finite element analysis (Girault and 

Raviart, 1986; Bristeau, Glowinski and Periaux, 1987). 

Our approach (Maday, Meiron, Rdnquist and Patera, 1987) 

to solving (2.21) is a global iterative procedure in which 

the original coupled saddle problem is decoupled into two 

positive-definite symnetric forms, 

(2.23a) -AX+ QTp- Bf , 

(2.23b) Sp = QA-lBf , 

where 

(2.24) S-QA-Q . 1T 

In the solution process (2.23b) is first solved for p and 

then (2.23a) is solved for uwith p known. 

The equation for the velocities -IL correspond to 

standard Laplacian solves, and the inversion of the matrix A 

is done by conjugate gradient iteration. Although the 

matrix S in thf equation for the pressure is full due to the 

presence of A- , the matrix is extremely well-conditioned. 

In particular, it can be shown that the matrix S is 

spectrally close to the variational equivalent of the 

identity matrix L, namely the m a s s m a t r i x B defined on the 

pressure mesh. This suggests an inner/outer conjugate 

gradient iteration procedure in which the outer iteration 

corresponds to inverting the matrix S preconditioned with 

the diagonal mass matrix B, and the inner iteration 

corresponds to inverting the discrfte Laplacian, A. As the 

condition number for the matrix B- S is of order unity, the 

above algorithm requires only order unity elliptic solves, 

and is therefore an ideal decoupling of the Stokes problem. 

conditioning of the matrix S for the semi-periodic model 

pioblem (2L14) and (2.15). 

Xk(n> of B S for the spectral element discretization with 

K-4, N-7, and wave number n-1. The clustering of the 

spectrum around unity is apparent, and a comparison with the 

continuous operator spectrum is seen to be virtually exact 

(Maday, Mei-ron, Rdnquist and Patera, 1987). For multidimensional Stgkes prpblem it can be shown that the 

conditioning n of B- S is still of order unity, allowing 

for rapid convergence of the outer ite ation. It should be 

noted that the good conditioning of B S is directly related 

to good approximation of the pressure (Maday and Patera, 

1987). 

The discretization and solution of the steady Stokes 

equations (2.14) can be readily extended to solve the 

unsteady Stokes equations. The implicit Stokes algorithm 

can be further extended to solve the unsteady Navier-Stokes 

equations by explicit treatment of the nonlinear term using 

a t h i r d - o r d e r A d a m s - B a s h f o r t h scheme. The resulting NavierStokes algorithm can be viewed as an implicit Stokes scheme 

with an "augmented" force which includes the explicit 

We now present numerical results demonstrating the good 

Figure 6 shows the spectrum 

-1 
10 

convective contributions. Numerous spectral element 

calculations of unsteady flows have been performed to date 

(Ghaddar, Korczak, Mikic and Patera, 1986; G h a d d a r , M a g e n , 

Mikic and Patera. 1986). 

3. Implementation on a Hypercube. 

Introduction. It is clear that the new direction of high 

performance architectures will entail large scale 

parallelism, and that successful algorithms must be able to 

exploit these new technological advances. The spectral 

element discretization and solution algorithms described in 

the previous section have been designed with parallel 

architecture in mind; in this section, we outline some 

parallel algorithm concepts, and describe the implementation 

of the spectral element method on the Intel iPSC Hypercube 

parallel processor. 

The Intel iPSC Hypercube is typical of the general 

class of parallel architectures for which our algorithms are 

appropriate. The Intel machine is a true multiple 

instruction-multiple data parallel processor, with M-2 

processors arranged in a classical hypercube fashion (Saad 

and Schultz,1985). Each processor has local source code 

operating on local data, with execution proceeding 

asynchronously until data is needed from another processor 

in the system. Information is exchanged between processors 

via a message passing scheme, in which data-send and 

data-receive statements are issued by ?he participating 

processors. Efficient use of all the processors can be 

achieved only if the ratio of communication time to 

computation time is small, and if all the processors have 

the same amount of work, that is, are balanced as regards 

their computational load. 

the Navier-Stokes problemwhich allow for successful 

parallelization are the following: 

The key aspects of the spectral element formulation of 

a) use of a domain decomposition approach; 

b) limited coupling between domains; 

c) high order discretization within domains; 

d) iterative solution techniques. 

Point (a) addresses the central issue of how one divides the 

problem among the available processors. Following the 

domain d e c o m p o s i t i o n a p p r o a c h , one or more spectral elements 

are assigned to each processor, with all the coefficients, 

geometry, and data associated with each element being 

locally resident. For our problem, the most natural 

decomposition is geometric, that is, adjacent elements 

reside on the same processor so that inter-processor 

comunication can be minimized. An alternative approach is 

to use a random element placement such that incremental work 

associated with local mesh refinement is distributed among 

several processors (Fox and Otto, 1985). 
11 

Points (b) and (c) address the requirement that there 

be minimal comnunication/computation if high parallel 

efficiencies are to be attained. 

condition at elemental interface~~impliej that the interelement communication will be O(N ) in R , while the highorder discr8tization implies that the computational work 

will be O(N (see Section 2). In essence, high-order 

methods imply a large amount of computation per point, thus 

resulting in a low ratio of communication to computation. 

Point (d) essentially addresses the load balance issue; 

iterative procedures such as conjugate gradient iteration 

are more r e a d i l y a m e n a b l e to node-simultaneous calculation 

than direct methods such as Gaussian elimination. By the 

same argument standard conjugate gradient preconditioning 

algorithms, such as incomplete Cholesky factorization, are 

not readily parallelizable; work is underway to find 

suitable alternatives. 

The Co continuity 

AlPorithmDescriptiu. At the heart of the parallel NavierStokes algorithm are the elliptic solvers which employ 

(unpreconditioned) conjugate gradient iteration (Golub and 

Van Loan, 1983). To illustrate the parallel implementation 

we describe the central step in this procedure, which is 

formation and evaluation of matrix vector products of the 

type: 

(3.1) 3-42 9 

where2 is the desired global result vector, p is an 

intermediate s e a r c h d i r e c t i o n , and A is the global Laplacian 

operator defined in Section 2. By definition of the 

variational statement (2.2-2.31, (3.1) can be written as: 

(3.2) 

where Ak is the local Laplace operatorkcoupling the nodal 

unknowns within the kth element, and p represents the 

values of p in element k. 

To compute r. in parallel, we first calculate 

(3.3) kk 

r. -k=Ap 

and then compute 

(3.4) 

k=1,2,. . . ,K , 

Here direct stiffness corresponds to summation of all 

e l e m e n t a l n o d a l v a l u e s corresponding to the same global 

degree-of-freedom. Note that in practice (3.4) is not a 

global operation, but corresponds to local exchange of face 

data between neighboring spectral elements. The sequence of 

steps (3.3-3.4) is depicted graphically in Figure 7. 

The only other computational step in the conjugate 

gradient algorithm r e q u i r i n g i n t e r - p r o c e s s o r . c o m m u n i c a t i o n 
12 

is the evaluation of inner products of the intermediate 

result vectors. This is readily handled using a 

substructuring technique, whereby local contributions to the 

inner product are computed on individual processors and then 

sumned via a binary-tree-like gather and sum. The final 

result is then redistributed to all the processors via the 

same binary tree, resulting in a total communication count 

of 2D, where D is the dimension of the hypercube network 

(Saad and Schultz, 1985). The communication associated with 

the inner products is therefore small compared to that 

associated with direct stiffness sumnation. 

Since the direct stiffness summation (3.4) invo ves 

only the faces of the spectrgl elements it is an O(N 

operation (per elementl4in R , whereas the matrix-vector 

product (3.3) is an O(N ) operation. The ratio of 

calculation to comunication is therefore O(N 1, which will 

typi.cally be quite large and will thus result in noncomnunication bound simulations. Note also that the direct 

stiffness summessage length of N2 will be relatively large, 

and thus startup costs associated with each data-send will 

be amortized over a large number of words, particularly on 

pipelined configurations. Thus, in theory, the spectral 

element/conjugate gradient algorithm suffers little 

comnunication overhead when implemented on a distributed 

memory parallel processor. 

practice, we plot in Figure 8 parallel gfficiency 77 vs. D 

for solution to a Poisson equation in R . Parallel 

efficiency is defined here as: 

1 

As an example of the comunication overhead incurred in 

(3.5) T1 77- . 

where M=2D is the number of processors, and T is the time 

required to solve a given problem on m processors. The 

efficiencies range from 0.96 to 0.99 for D ranging from 0 to 

5. The K-32 spectral element discretization used in this 

example is o f m o d e r a t e o r d e r , N=6; higher order solutions 

will show a further increase in parallel efficiency. 

strong function of the particular hardware used. The 

results of Figure 8 were obtained on a D-5 iPSC Hypercube 

without high performance vector processors. Increases in 

floating point performance will increase the absolute 

calculation speed, but will at the same time make the 

parallel efficiency worse. The question remains as to 

whether comunications capability will continue to keep pace 

with the rapid advances in floating point hardware. 

rn 

The actual efficiency obtained will of course be a 

Virtual Parallel Processor. Central to the development of 

the parallel spectral element algorithm is the treatment of 

each element as a separate entity, that i s , a s a virtual 

processor. All data structures are set up on an element by 

element basis, and all operations are executed element by 
13 

element on this element-local data. In this way, 

parallelism is inherent in the a l g o r i t h m , a n d not overly 

dependent upon the particular machine on which the code is 

run. 

The practical implementation of this virtual parallel 

processor concept is relatively simple. The critical 

feature is that all outer loops in the code run over the 

number of elements; on a serial machine, loops extend over 

the entire domain, while on a parallel machine, loops extend 

over the subset of elements associated with a given 

processor. The source code on e a c h o f the processors is 

therefore identical. The spectral element-virtual processor 

association further insures that the few operations 

requiring comunicat ion between processors (direct stiffness 

sunmation, inner product evaluation) can be easily effected 

by device drivers which translate s e n d a n d receive 

subroutines into hardware-compatible datcsend and data 

receive comnands. 

Alternative AlgPritWArchitecture Coup-. We now 

consider our particular choice of architecture and algorithm 

in light of other available alternatives. Our algorithm is 

well suited for what we define as medium-grained parallel 

processing. Defining Nt to be the total number of degreesof-freedom, M to be the number of processors, and p = Nt/M, 

we classify systems as 

coarse-grained: Nt =+m, M fixed, P + 0 , 

medium-grained: N, =+m, M=+oo, p<<l, (p =& 0) ; 

C o a r s e - g r a i n e d m a c h i n e s do not exploit the full parallel 

potential of our algorithm and are consequently not of 

interest. It is also fairly clear, u priori. that highorder methods will not readily parallelize on fine-grained 

s y s t e m s d u e to comunication considerations. We therefore 

restrict our comparison to high-order/medium-grained 

couplings versus low-orderlfine-grained couplings. 

We first note that a low-order approach is necessarily 

going to require many more points to achieve the same level 

of accuracy as the spectral element method. Hence, the 

fine-grained system will need many more processors than the 

number of degrees-of-freedom in the spectral element 

solution. Second, for e a c h o f the fine-grained processors 

the ratio of comunication time to computation time will be 

relatively large, as the number o f w o r d s transmitted during 

direct stiffness summation will be of the same order as the 

number of calculations required to evaluate the local 

residual. Third, the messages will of necessity be short on 

a fine-grained processor, and the advantages of pipelined 

communication may not be realized. 

Next, we consider the issue of mesh refinement. 

Adaptive resolution is readily implemented in the spectral 

elementlmedium-grained algorithm. This is effected by 

fine-grained: Nt+m, Maw, p =+ 1 
14 

increasing the order of the polynomial ap.proximation in the 

vicinity of the desired high resolution area (Mavripilis); 

due to spectral accuracy, an order o f m a g n i t u d e improvement 

can be obtained at the cost of only a few additional points 

in each s p a t i a l d i r e c t i o n . B e c a u s e of this exponential gain 

in accuracy, and because the work per processor is already 

large, the incremental work incurred due to refinement will 

be a f r a c t i o n o f the total work on the processor. 

Therefore, the processor load imbalance resulting from this 

localized refinement will not be significant. In contrast, 

mesh refinement on a fine-grainedllow-order implementation 

will probably require at least a doubling of the number of 

points on the relevant processors to improve local 

resolution. The incremental work increase will then be of 

the same order as the initial work, and the parallel 

efficiency for an initially well balanced problemwill be 

reduced by a significant factor unless a large-scale load 

redistribution is carried out. 

associated algorithmmust be coupled to economic 

considerations which are difficult to predict. The above 

discussion should only be considered as plausibility 

arguments for a h i g h - o r d e r m e t h o d / m e d i u m - g r a i n e d coupling. 

Ultimately, the choice of parallel architecture and 

4. Geometry Defining Processors (GDPs). 

tivation. The use of the hypercube for solution of finite 

or spectral element equations as described in Section 3 

results in efficient parallel solution of' partial 

differential equations (PDE) in three dimensional 

geometries. However, when applied to this class of 

problems, the general-purpose hypercube comunications 

network is too flexible, and, at the same time, non-optimal. 

In particular, as more processors are applied to the 

solution of a problem, the hypercube's log(M) connections 

will grow, and eventually exceed, the number required for 

nearest neighbor comunication; these unnecessary 

connections represent a significant increase in per 

processor cost and complexity. Furthermore, even with the 

available lpg(M) connections, computational problems in 

complex geometries will require some inter-processor 

conmunication to be routed through intermediate nodes, 

potentially resulting in increased communication latency. 

parallel processing issue that has not yet been explicitly 

considered in this paper. This is the mapping problem 

(Bokhari, 1979)., which addresses how spectral (or finite) 

elements should be distributed onto a particular processor 

configuration. This mapping is typically performed by a 

(serial) front-end, which is also responsible for the 

related activity of mesh generation. A number of computeraided design systems have been designed to facilitate the 

transfer of a problem geometry from the physical domain to a 

computer representation, however the input and d i s p l a y o f 

three dimensional objects and d o m a i n s v i a essentially two 

The issue of inter-processor communication raises a 
15 

dimensional tools is often cumbersome, and contributes 

little to solution of the associated parallel processing 

mapping problem. 

and convenient mesh generation appear superficially 

unrelated, it can be seen that they are, in fact, intimately 

connected. In essense, mesh generation is related to the 

problem of specifying domain and connection topology, which 

is, in turn, a central issue in the subsequent parallel 

solution of the governing PDEs. In this section we briefly 

describe a special purpose parallel architecture that 

simulateously addresses the mesh generation and parallel 

solution problems by exploiting their underlying 

similarities (Dewey and Patera, 1987). 

Although the problems of efficient parallel solution 

Processors. Geometry-Defining Processors 

(GDPs) are microprocessors housed in physical geometrfc 

packages which can be easily manually assembled or 

reconfigured to define any particular system geometry. An 

individual GDP is aware of the parameters of its physical 

package, is able to communicate with neighboring GDPs as 

well as with a host processor, and is capable of performing 

independent numerical computations. To solve a PDE, that 

is, to simulate a physical system, GDPs are assembled in a 

"scale" model of the actual geometry, as shown in Figure 9. 

The GDP assembly provides a means for geometry input, and 

subsequently functions as an optimally-connected dedicated 

parallel processor for the particular problem at hand. 

GDP is shown schematically in Figure 10 for the simple case 

of a rectangular two-dimensional GDP. A communication port 

on each face of the GDP allows for two types of 

communication within the GDP assembly. First, each GDP is 

able to communicate with neighboring GDPs through local 

busses. The local bus is reconfigurable, in that the lines 

emanating from different faces can be internally connected 

to route messages through the GDP; this allows high-speed, 

parallel, non-nearest-neighbor communication (e.g., for 

vector reduction operations (Lin and Sips, 1986)). The 

second type of comnunication supported by the GDP faces is a 

global bus which all GDPs and the host computer can access. 

This bus is primarily used for communication between the 

host computer and one or more of the GDPs. 

Each G D P - r e s i d e n t m i c r o p r o c e s s o r , in addition to 

controlling communications, performs local calculations 

related to mesh generation or equation solution using 

programs and parameters downloaded from the host computer 

via the global bus. The GDPs in an assembly are 

synchronized by the global bus, however they otherwise 

operate independently, performing autonomous calculations on 

locally-available data. 

following sections present details of the o p e r a t i o n o f GDP 

assemblies when applied to mesh generation and parallel 

solution. 

The basic processor/communication structure of a single 

With the above description of an individual GDP, the 
16 

Mesh Generation with GDPs. GDP mesh generation proceeds 

t h r o u g h s e v e r a l s t e p s of comunication and computation. 

After (manual) assembly of the desired geometry, the host 

interrogates the assembly of GDPs to determine the number, 

names, and geometry types of the constituent GDPs. Next, 

each GDP is instructed to determine and report any 

connections between itself and its neighbors. By virtue of 

their physical p r o x i m i t y , a n d hence communication proximity 

through the local busses, the information needed to 

determine these connections is locally available to each 

GDP. F i n a l l y , b y "integrating" the geometry-type and 

connection information obtained from the GDP assembly, the 

host-resident software can determine the absolute geometry 

of the assembly. 

With a philosophy similar to that of the graphics 

tablet, the GDPs allow for data input directly in physical 

space, with all symbolic, computational-space processing 

performed transparently to the user. Both the minimal 

learning curve associated with initial use of GDPs, and the 

minimal time required for subsequent particular realizations 

of GDP assemblies, should make the devices a useful addition 

to conventional software geometry-input techniques. 

P a r a l l e l S o l u t i o n with GDPs. Geometry-Defining Processor 

solution of partial differential equations proceeds as 

follows. The physical domain is first (manually) 

constructed as shown in Figure 9, and the system geometry is 

then determined by the host as described in the previous 

subsection. It should be noted that the geometry as input 

by the GDPs need only represent a "rough-cut" of the actual 

system geometry, as conventional computer-aided design tools 

can be used to subsequently tailor the GDP shapes in 

software. 

a substructure of finite or spectral elements. Note that 

the mapping problem has been completely eliminated: with no 

computational overhead the processor/elements are now 

optimally connected in that neighboring elements in problem 

space are represented by connnunicating processors in 

computational space. In addition to this nearest-neighbor 

coxnnunicatiOn many solution algorithms, such as conjugate 

gradient iteration, require global summation and broadcast; 

these vector reduction operations can be efficiently 

provided through the use of the reconfigurable local bus. 

Thus, the GDP assembly next determines a (time-dependent) 

re-routing of the local busses that will allow these 

operations to be performed in log(M) time (Dewey and Patera, 

1987). Finally, the discrete equations, physical 

p r o p e r t i e s , a n d boundary conditions for each GDP are 

downloaded from the host processor. 

preliminaries performed, the actual parallel s+olution is 

then initiated and subsequently synchronized by the host 

through the global bus. The operations carried out by the 

GDPs in this solution stage are essentially identical to 

Next, each GDP is assigned a single, or more generally, 

With these mesh generation and equation definition 
those now implemented on the hypercube processor as 

described in Section 3, with the advantage of an optimal 

connection topology that results in reduced communication 

latency and greatly reduced hardware and software 

complexity. 

Development of Prototvpe GDPs. We have developed prototype 

hardware GDPs and associated host software to perform simple 

three-dimensional geometry input. The GDP shapes chosen for 

these prototypes are the cubical and wedge elements shown in 

the example of Figure 9. These prototype GDPs have faces 

outfitted with optical comunication ports for the local and 

global busses; a typical face is diagramed in Figure 11. 

Optical emitters and detectors (rather than 

physical/electrical connections) are used in order to have 

rotationally independent, easily reconfigurable interfaces. 

A central port on each face is used for global 

co&unication, while the four ports along the periphery of 

each face allow for local comunication and the 

determination of the rotational orientation of communicating 

faces. A 16 bit microprocessor with nominal amounts of RAM 

and EPROM in each GDP controls the communications ports and 

runs the host-downloaded mesh generation program. 

In operation, these prototypes allow simple threedimensional geometries to be conveniently input and 

reconfigured. With these prototypes as a starting point, 

continued research aims to develop sufficient hardware and 

software to allow for a meaningful evaluation of the GDP 

concept as a mesh-generation and parallel solution 

architecture. 

Conclusion. In this paper, we have shown how high-order 

spectral element algorithms coupled to medium-grained 

message-passing or geometry-defining architectures can 

result in significant improvements in performance over 

conventional serial s o l u t i o n m e t h o d s . The results presented 

demonstrate that the use of parallel processors for the 

solution of complex problems is no longer only a promise, 

but in fact a reality. Future gains in computational power 

will be the result of continued emphasis on parallelism and 

algorithm-architecture coupling. 

Acknowledgements. We would like to acknowledge our 

extensive collaborations with Yvon Maday of the Universite 

Pierre et Marie Curie, and Daniel I. Meiron of the 

California Institute of Technology. This work was supported 

by the ONR and DARPA under Contract N00014-85-K-0208, by the 

NSF under Grant DMC 87-04357, by NASA under G r a n t N A G 1 - 5 7 4 , 

by the Massachusetts Institute of Technology, by Intel 

Scientific Computers, and by Nektonics, Incorporated. 

REFERENCES 

(1) S.H. BOKHARI, Qn the maDDinp problem , in Proc. 1979 

Int. Conf. on Parallel Processing, 1979, p. 239-248. 
18 

(2) C. BERNARD1 AND Y. WAY, A collocation method over 

in Fluids, to appear, 1987. 

ered prids for the Stdces problem, Int. J. Num. Meth. 

.. (5) D. DEWEY AND A.T. PATERA, Geometrv-deflnlng processors 

for partial differential eauatiu, in Architecture and 

performance of specialized computer systems, ed. B.J. Alder, 

Academic Press. 1987. 

(6) G.C. FOX AND S.W. OTTO, Concurrent cqmputation and the 

Sheorv of cqmplex svst-, in Hypercube Multiprocessors, ed. 

by M.T. Heath, SIAM, 1986. 

(7) D. FUNARO, Altidornain spectral aIrgroXimation of 

 tic eau-, Nu. Meth. for Partial Differential 

Equations, 2 (19861, p. 187. 

(8) N.K. GHADDAR, K.Z. KORCZAK, B.B. MIKIC AND A.T. PATERA, 

els. Part 1: Stabilitv and self-sustained osclllatu, 

rical investi~ation - of incpmpressible flow in grooved .. 

J. Fluid Mech., 163 (19861, pp. 99. 

(9) N.K. GHADDAR, M. MAGEN, B.B. MIKIC AND A.T. PATERA, 

rical investigation of incqm~ressible flow in grooved 

els. Part 2: Resonance and oscillatorv heat transfer, 

J. Fluid Mech., 168 (19861, p. 541. 

(10) V. GIRAULT AND P.A. RAVIART, Finite el1986. 

tion of the Navier-Stokes equationa, Springer, 

(11) G.H. GOLD AND C.F. VAN LOAN, Matrix hnputatiom, 

Johns Hopkins University Press, 1983. 

(12) D. GOTTLIEB AND S.A. ORSZAG, Numerical anaLvsis of 

-methods, SIAM, 1977. 

(14) H.X. LIN AND H.J. SIPS, A patallel vector reduction 

architecture, in Proc. 1986 Int. Conf. on Parallel 

Processing, 1986, p. 503-510. 
19 

(15) Y. MADAY AND A.T. PATEIU, Spectral element methods for 

oe iwressible Navler-Stokes equations, in State-of-theArt Surveys in Computational Mechanics (ed: A.K. Noor), 

ASME, 1987. 

(17) Y. MADAY, A.T. PATERA AND E.M. RBNQUIST, Qglimal 

perlodlc problem, in preparation, 1987a. 

(18) Y. MADAY, A.T. PATERA AND E.M. RBNQUIST, ODtimal 

Le p end1 e s pe c t t1 t 1 -d mlonal 

Stok.keem, in preparation, 1987b. 

thods for the Stokes ScIXlL 

(19) C. MAVRIPLIS, Ph.D. Thesis, Dept. of Aero. & Astro., 

Massachusetts Institute of Technology, in progress. 

.. (20) J.A. MEIJERINK AND H.A. VAN DER VORST, An iterative 

solution method for 11near sv-s of which the coeffrclelu 

tric M-matrk, Math. Comp., 31 (19771, p. 

148. 

(21) S.A. ORSZAG, Spectral methods for problemsinex 

-, J. Comput. Phys., 37 (19801, p. 70. 

(22) A.T. PATERA, Apectral element method for fluid 

Phys., 54 (19841, p. 468. 

cs: unar flow in chamel expanslou, J. Comput. 

(23) E.M. R0NQUIST AND A.T. PATERA, A LeEendre spectral 

Engr., to appear, 1987. 

(24) Y. SAAD AND M.H. SCHULTZ, Topological properties of 

h-, Research Report YALEU/DCS/RR-389, Yale 

University, New Haven, June 1985. 

nt method for the Stefan problem, Int. J. Num. Meth. in 

.. (25) G. STKANG AND G.J. FIX, Anvsis of the FULL 

Element, Prentice-Hall, 1973. 
20 

FIGURE CAPTIONS 

Figure 1. A plot of the maximum nodal error in 

the spectral element solution to the differntial 

equation (2.1) as a f u n c t i o n o f the total number 

of degrees-of-freedom, Nt. In this problem A-0, 

p(x)-e , f(x>-e (cosx-sinx), for which the 

solution is u--sinx on x E ]O,T[. The domain is 

divided into K-2 spectral elements of equal 

length. Exponential convergence is achieved as 

the polynomial degree of the fixed elements is 

increased. 

X X 

Figure 2. The domain and spatial discretization 

for s o l u t i o n o f the Poisson equation (2.8a) with 

f-0. Dirichlet boundary conditions are imposed 

such that the solution is ~==sinx.e-~. The domain 

is divided into K-4 spectral elements. 

Figure 3. A plot of the maximum nodal error in 

the Legendre spectral element solution of the 

Poisson equation (2.8a) as a function of the 

total number of degrees-of-freedom in one 

spatial direction, N,. Exponential convergence 

is achieved as the degree of the the elements is 

increased. 

Figure 4. Convergence history for preconditioned 

(closed symbols) and unpreconditioned (open 

symbols) conjugate gradient solution of a 

spectral element Poisson discretization. In all 

cases K1=4, with N-5(A), 7(0), 9(0), and ll(v). 

Figure 5. The error in the velocity uN=(uN,vN) 

and the pressure p as a f u n c t i o n o f the total 

numbe r of degree -op- f reedom (Gaus s -Loba t to 

Legendre points) in the x-direction, N,, when 

solving the test problem (2.22). The total 

interval d-1-l,l[ is divided into K-2 spectral 

elements A,-]-l,O[ and A2=]0,1[. Exponential 

convergence is obtained. 

Figure 6. A plot of thz-fpectrum A (n) of the 

preconditioned matrix B 

pressure matrix given in (2.24) and E is the 

mass matrix defined on the Gauss pressure mesh. 

The spectrum (A) corresponds to a spectral 

element discretization K-4, N-7 for a wavenumber 

n-1; the agreement with the continuous operator 

spectrum (0) is very good. 

S 

S, where 5 is the 

Figure 7. Caption on figure. 
Figure 8. Parallel efficiency vs. hypercube 

dimension for solution to three-dimensional 

Poisson equation using sixteen P spectral 

elements. 6 

Figure 9. Geometry-Defining Processor assembly 

corresponding to a physical domain D, which in 

this case is a block cooled by a heat-sinking 

fin. The GDP subsystem consists of the assembly 

of GDPs, a GDP interface block, and a controller 

through which the host communicates with the 

GDPs . 

Figure 10. Schematic of a (rectangular, twodimensional) Geometry-Defining Processor. A GDP 

consists of a package, comnunication ports on 

each face, a microprocessor, a floating point 

unit (FPU), RAM and EPROM, a communications 

controller, a global bus (GB), and a 

reconfigurable local bus (LB). 

Figure 11. Layout of a typical face of the 

prototype GDPs, showing the arrangement of 

e m i t t e r s a n d d e t e c t o r s u s e d for local and global 

comnunication. The use of four local detectors 

allows rotational orientation to be determined. 
22 

PRIGINAL PAGE IS 

03 POOR QUALITY 

1 oo 

lo-' 

1 o-2 

1 o-~ 

1 o-~ 

1 o-~ 

1 o-6 

1 o-' 

1 o-8 

1 o-~ 

I 1 I I I I I I I I I I I I I I I I I 

10-'O 

i= k lo-" 

lo-'* 

t 

0 

0 

0 

0 

0 

0 

0 

0 

0 

I I I I I I I I I I I I I I 1 I I I I I I I 

2 4 6 8 10 12 14 16 18 20 22 24 10-1~ I 

N, 

Figure 1. 
23 

x 

x 
24 

1 o-2 

1 o-~ 

1 o-+ 

1 o-s 

1 o-6 

1 o-' 

1 o-8 

1 o-~ 

1 o-'O 

lo-" 

1 0-l2 

0 

0 

0 

0 

0 

0 

0 

0 

t 

Figure 3. 
10' 

1 oo 

lo-' 

1 o-2 

1 o-~ 

1 o-~ 

1 o-s 

1 o-6 

1 o-' 

1 o-8 

CONVERGENCE HISTORY - N=5,7,9,1 1 K1 14 

I I I I I I I I I I I I I I I I I I I I 

L 

b 

V 

Ll v 0 

0 0 

0 A 

V 

V 

0 

V 

IO 

V 

0 n 

0 

0 

0 

U 

0 

0 

0 

A 

0 

0 

n 

V 

V 

V 

V 

0 

0 

0 

V 

V 

-1 

in 

Figure 4. 
26 

1 oo 

lo-' 

1 o-2 

1 o-~ 

1 o-~ 

1 o-~ 

1 o-6 

1 o-' 

1 0-a 

1 o-~ 

1 o-'O 

lo-" 

0 

0 Ago 

0 

BA 

00 

a 

A 

0 

0 

n 

0 

A 

00 

0 

A 

0 

A 

IIV - v N 

I IP - Px 

00 

a 

AA 

0 

0 

0 

0 

c 

4 6 8 10 12 14 16 18 20 22 24 26 25 30 

Figure 5. 
27 

1.3 I 1 I I I I I I I I I I I I I I 1 I I I I I 

I I I I I I I I I I I I I I I I I I I I I 

0 2 4 6 8 10 12 14 16 18 20 22 0' 

k 

Figure 6. 
28 

p.' A- 

-1 * ir:i 

global local 

(a) 

-k o form r 

Compu t a t 1 

\, 

I I 

I I 

I I 

o d'rect stiffnes sum 

r in parallel 

Communicatfon O(1) 

d 1 

(C) 

Figure 7. One dimensional example of parallel matrix-vector 

multiply, r=&, for three spectral elements. (a) Global and 

local representations of intermediate s e a r c h d i r e c t i o n p, 

where * indicates nodab values at element interfaces which 

must be equal by the C continuity requirement. (b) Global 

form of the discrete Laplace operator, A, where + indicates 

sunmation of overlapping contributions from local matrices. 

(E) Expansion of global product-into three local products, 

r , evaluated in parallel. The r are th n direct stiffness 

sumned to form the final result vecior, L . Final ratio of 

computation to conxnunication is O(N in all space 

dimensions. 

E 
29 

D 

Figure 8. 
30 

c 

\ \ \ 

\ 

z 

0 

A 

W 

m 

ICE 

\ 

\ 

Z - 

a 

z 

0 

0 

J 

0 

cn 

>- 

I 

a. 

a 
U 

P 

31 

cn 

7 

0 - 

c 

W 

0 

LL 

a 

\ 

m 

c3 

a 

m 

0 

-J 

m 

J 

- 

r 
S 

32 

)E VIEW FRONT VIEW 

...... q& 

P 

I 

\ / 

i\ I // 

I \ /’ 

I \ I&’\, I 

I 

I 

I 

rI 

\ 

/\\ 

\ / 

\ ‘a; 

‘f / 

/\ 

1 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

1 

I 

I 

I 

I 

I 

/ / \] 

\ \ 

GLOBAL EMITTER “LOCAL EMITTER 

AND 

BEAM SPLITTER /DETECTOR /J 

LOCAL DETECTOR (4) 

Figure 11.