─юъєьхэЄ тч Є шч ъ¤°р яюшёъютющ ьр°шэ√. └фЁхё юЁшушэры№эюую фюъєьхэЄр : http://geophys.srcc.msu.ru/library/parcfd2003.pdf ─рЄр шчьхэхэш : Tue May 22 13:15:20 2007 ─рЄр шэфхъёшЁютрэш : Mon Oct 1 19:46:24 2012 ╩юфшЁютър: IBM-866

1

Parallel Computations in Problems of Climate Modeling
V. Gloukhov
a a



Research Computing Center, Moscow State University, NIVC MGU, Vorob jovi Gory, 119992 Moscow, Russia The pap er is fo cused on different approaches to parallel implementation of climate mo dels on cluster multipro cessors. We present results of applying either MPI or Op enMP to the mo del comp onents: atmospheric and o ceanic blo cks. A library of communication routines was develop ed for the distributed memory approach. It carries out the basic communication op erations, such as b oundary exchanges and transp ositions of decomp osed data. The library has quite general interface that makes it useful for parallelization of a wide range of scientific applications calculated on structured grids. In particular, we have examined it on INM AGCM as well as on a set of atmospheric b enchmarks which integrate the equations of hydrothermo dynamics of atmosphere under the hydrostatic assumption. An Op enMP implementation of the INM Ocean mo del has b een undertaken and tested on few shared memory systems. INTRODUCTION The generally accepted definition of climate is the average course or condition of weather at a place usually over p erio d of years as exhibited by temp erature, wind velo city, precipitation, etc. Or more rigorously, it is a set of subsequent states of atmosphere, o cean and criosphere, quantified by state variables. Since no direct physical exp eriment on climate system is p ossible, scientists use state-of-the-art numerical mo dels to access the role of anthrop ogenic factors in its formation and changes. Long p erio ds of integration, decades or even centuries, make the problem complex to b e pro cessed by a unipro cessor PC or workstation, necessitating the usage of parallel computers, in particular those of cluster architecture, comprised of SMP no des interconnected with fast network. The concept of parallel pro cessing is quite suitable for climate mo dels b ecause they can b e split up into comp onents, for example, into the atmospheric and o cean parts, which interact each other from time to time, through exchange of energy and substance, and are indep endent in-b etween. Furthermore, the comp onents, also can b e calculated concurrently, providing a substantial degree of parallelism exploitable through shared or distributed memory multitasking. This pap er aims to show few of our exp eriences in trying different parallelization strategies for an Atmospheric and an Ocean mo del develop ed in the Institute of Numerical Mathematics of the Russian Academy of Sciences. We highlight some technical details of these approaches that could b e interesting to application develop ers, not only from the


Supported by the INTAS pro jects 00-189, 01-2132 and RFBR pro jects 01-05-64150a, 03-05-64358a.


2 area of climate studies but from other fields also. While an application of distributed memory paradigm to the full INM atmospheric general circulation mo del was describ ed in [1], here we stress only on distributed memory implementation of its dynamical core resp onsible for integration of the governing equations. And the Op enMP paradigm was tested on the INM Ocean mo del that constitutes the o cean comp onent of a coupled mo del [10]. We kept the following structure of the presentation. In section 1, a library of communication routines (communication kernel) is describ ed that p erforms the basic communication op erations on distributed memory: b oundary exchanges and transp ositions. Written in C the library is intermediate b etween MPI and application. The communicated data may have quite general format which makes it suitable for parallelization of a wide range of b oth C and Fortran scientific applications formulated on structured grids. In particular, the library are invoked by the atmospheric dynamical core to maintain 2D data partitioning over the pro cessors. Extended by the idealized Held-Suarez forcing [2], the dynamics can b e used for b enchmarking purp oses alone. It allows to switch b etween the explicit and semi-implicit time stepping scheme, as well as b etween direct and iterative solvers for the Helmholtz equation arising in the semi-implicit scheme, giving an opp ortunity for direct intercomparison of the metho ds at different grid resolutions. Section 2 contains an example of such b enchmarking carried out on MBC1000M computing system installed at the Joint Sup ercomputing Center, in Moscow. Though MPI is quite common programming environment on clusters, other techniques exist facilitating parallelization. In section 3, we present the results of an Op enMP parallelization of the INM Ocean mo del, which requires less programming efforts than the distributed memory approach, but scalable up to a single SMP no de on most systems. 1. COMMUNICATION KERNEL If a climate mo del runs on a distributed memory machine and domain decomp osition technique is used within its comp onents, each pro cessor computes only particular subdomain of either the atmosphere or o cean. Interactions b etween the mo del comp onents and dep endencies at the b oundaries of sub domains urge the pro cessors to interchange messages. These communications have to b e realized as a set of routines. We are ab out to characterize a library of such routines, called communication kernel, which is based on MPI. This library was designed for the communications p erformed within the mo del comp onents for it has no interp olation features. Moreover, the communicated data are assumed to b e stored in multidimensional arrays which dimensions corresp ond to different spatial dimensions or prognostic variables. Therefore, the library is intended only for the structured grids and rectangular domains which should b e decomp osed along one ore more spatial dimensions. It p erforms two ma jor kinds of communications of distributed data: b oundary exchanges and transp ositions. Boundary exchange is a lo cal communication, in the sense that every pro cessor exchanges values adjacent to the b oundaries of his sub domain only with the neighb ors, while the transp osition is global, or an all-to-all communication, such that each pro cessor within a group interacts with each. For example, b oundary exchanges are inevitable, when derivatives are evaluated by the metho d of finite differences. But if an op eration, for


3 instance integral summation, requires all the data along a decomp osed dimension and is indep endent along another, which is not partitioned, the data can b e redistributed along the second dimension, then the op eration is evaluated lo cally within each pro cessor, and in the end the original distribution is restored. That is exactly what we call transp osition. The b oundary exchange routine has the following interface int P_BExchange( void *arr, int nd_array, int *stride, int *blklen, int nd_procgrid, int n_edges, int *cart2arr, int overlap[][2], MPI_Datatype datatype, MPI_Comm comm ) CALL P_BEXCHANGE (ARR, ND_ARRAY, STRIDE, BLKLEN, ND_PROCGRID, N_EDGES, CART2ARR, OVERLAP, DATATYPE, COMM, IERR) INTEGER ND_ARRAY, ND_PROCGRIG, N_EDGES, DATATYPE, COMM, IERR INTEGER STRIDE(ND_ARRAY), BLKLEN(ND_ARRAY), CART2ARR(ND_PROCGRID) INTEGER OVERLAP(2, ND_PROCGRID) where COMM MPI communicator defining a Cartesian top ology, ARR p ointer to the first element of lo cal sub domain, ND ARRAY numb er of dimensions of array ARR, ND PROCGRID numb er of dimensions of the Cartesian top ology, CART2ARR array mapping dimensions of the Cartesian top ology into dimensions of the array ARR, STRIDE dimensions of array ARR, BLKLEN dimensions of the sub domain, OVERLAP, N EDGES width and adjacency degree of b oundaries to b e communicated, resp ectively, DATATYPE MPI data typ e of array elements, IERR error co de. Figure 1 depicts a lo cal to some pro cessor p ortion of data, decomp osed along directions dir1 and dir2 , with its halo regions (b oundaries). We assume, by definition, that the 00 01 10 0 1 0 1 b oundaries B1 , B1 , B2 and B2 has the first degree of adjacency, while B1,,2 , B1,,2 , B1,,2 11 and B1,,2 the second (the last four will b e exchanged in diagonal directions). The transp osition routine has the following interface int P_Transpose ( ndims, void *arr_source, int dim_source, int *lblks_source, void *arr_dest, int dim_dest, int *lblks_dest, int *stride, int *blklen, int *overlap, MPI_Datatype datatype, MPI_Comm comm, int period ) CALL P_TRANSPOSE (NDIMS, ARR_SOURCE, DIM_SOURCE, LBLKS_SOURCE, ARR_DEST, DIM_DEST, LBLKS_DEST, STRIDE, BLKLEN, OVERLAP, DATATYPE, COMM, PERIOD, IERR) INTEGER NDIMS, DIM_SOURCE, DIM_DEST, DATATYPE, COMM, PERIOD INTEGER IERR, LBLKS_SOURCE(*), LBLKS_DEST(*) INTEGER STRIDE(NDIMS), BLKLEN(NDIMS) where ARR SOURCE, ARR DEST p ointers to lo cal sub domains within the arrays containing data in original and transp osed distribution, resp ectively, NDIMS numb er of dimensions of b oth arrays ARR SOURCE and ARR DEST, DIM SOURCE and DIM DEST distributed dimensions


of the comminicated arrays, resp ectively, LBLKS SOURCE and LBLKS DEST distribution of dimensions DIM SOURCE and DIM DEST, resp ectively, STRIDE dimensions of the communicated arrays, BLKLEN size of lo cal sub domains, COMM MPI communicator determining the pro cessors participating in transp osition, DATATYPE MPI data typ e of the communicated arrays, IERR error co de. Parameters OVERLAP and PERIOD define the interpro cessor overlapping of redistributed data and p erio dicity of dimension DIM DEST. There usage is equiualent to a transp osition without overlapping followed by a b oundary exchange in the direction of dimension DIM DEST. Figure 2 represents data in array ARR SOURCE decomp osed in direction dir s (original distribution), while figure 3 shows the same data in array ARR DEST redistributed in direction dir d (transp osed layout). To transfer data asynchronously the kernel routines call MPI non blo cking send-receive functions MPI Isend and MPI Irecv [4], while to form the blo cks of data to b e transfered they invoke the MPI Type vector function. It is worth to mention, that our kernel functionally is similar to the Nearest Neighb or To ol (NNT) develop ed by Ro driguez, Hart, and Henderson [3].

Figure 1. The lo cal sub domain and its halo regions for a 2D data partitioning over the pro cessors.

стс ст стс стс стст стст стст с с стс стст ст стст

с с с с стс стст ст стст с с с с с стст стст с с с с с стст с ст с с с с стс стст ст стст с с с с стс с стст

4

B

w

1,0 1,2

0 2

blk

B

1 1

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с с с с стст с с с с с с стст стст

2. EULERIAN DYNAMICAL CORE

Robustness of the communication kernel describ ed ab ove has b een tested on the dynamical core, used in INM AGCM [5], which integrates a system of primitive equations

стс ст стс стс стст стст стст с с стс стст ст стст с с с с с с стст с стс ст с с с с с стст стст с с с с с стст стст с с с с стс стст ст стст с с с с стс стст ст стст с с с с с стст с ст

стс стс стс стст стст с с с с с стст стст с с с с с стст стст

стс стстс стс стст ст стст с с с с стс стст ст стст с с с с стс стст ст стст с с с с с стст с ст

стс стстс стс стст стст с с с с стс стст стст с с с с стс стст стст

стс стстс стс стст стст с с с с стс стст стст с с с с стс стст стст

стс стстс стс стст стст с с с с стс стст стст с с с с стс стст стст

стс стстс стс стст стст с с с с стс стст стст с с с с стс стст стст

стс стстс стс стст стст с с с с стс стст стст с с с с стс стст стст

стс стстс стс стст стст с с с с стс стст стст с с с с стс стст стст

сфстсф с сф стст стст сус су стсус ст с с с с с стст стст с с с с с стст стст

тф фф т тт ст тт у у ту т т тт ст тт т тт ст тт тт с т с с

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

стс стс ст стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

с с с с с стст с ст

сф сусфс стсу стсус стсф стс с стс с ст стс ст стс с стс с ст стс ст стс с стст с с стст с с ст с с с с с стст

dir1

B

B

0,0 1,2

0 2

( f 1, f 2)

dir

2

subdomain

B

0 1

dir 2

B B

w

B

0,1 1,2 1,1 1,2 1 2 1 2

w w

blk

0 1 1 1

dir 1


of hydrothermo dynamics of the atmosphere in the Boussinesq approximation under the hydrostatic assumption [6]. Written in vertical -co ordinate, the system is discretized on a staggered Arakawa C-grid [7]. Our realization allows user to cho ose b etween explicit or semi-implicit time integration scheme. Since most calculations p erformed by the AGCM are very dep endent in the vertical direction, we applied checkerb oard partitioning, distributing data along b oth longitude and latitude, over two-dimensional Cartesian top ology of pro cessors. Evaluation of the explicit tendencies at a grid p oint requires the values of prognostic variables from some neighb orho o d of that p oint, necessitating a b oundary exchange. Also a transp osition is required to p erform Fourier filtering of the fields near the p oles. The semi-implicit scheme increases time steps 4э5 times, but leads to a discrete twodimensional Helmholtz-like equation that have to b e solved on every vertical level each time step. The matrix of the arising linear system is structured and can b e represented 2 as I - qk [D M + (D M ) I ], where qk are height dep endent co efficients, I is a unit matrix, D is diagonal, M tridiagonal, and M cyclic tridiagonal. There are few fast direct metho ds of solving the system with structured matrix, like we have, which take advantage of the matrix structure. We exploit a direct metho d involving the fast Fourier transform (FFT) along longitude and tridiagonal Gaussian eliminations along latitude. However, b oth longitudes and latitudes are distributed, and the alternative here is whether to use distributed algorithms of FFT and tridiagonal Gaussian elimination or transp ose the data twice, first, in longitude-latitude plane and then, after completion of FFTs, in latitude-longitude plane to p erform Gaussian sweeps. Our exp erience convinces

с с стс ст с с с стс ст с с с стс ст с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст с с стс ст с с с стс ст с с с стс ст с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст с с стс ст с с с стс ст с с с стс ст с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

Figure 2. Data distribution b efore transp osition.

с с стс ст с с с стст с ст с стст с с с с с стст с ст с стст

с с стс ст с с с стст с с стст с с с с с стст с ст с стст

с с стс ст с с с стс ст с с с стс ст с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стс ст с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стс ст с с с стс ст с с с с стст с

с ст стст с ст стст с ст стст с ст стст с ст с стст с с стс ст с с с стс ст с с с стс ст с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с ст ст с с стстс ст сщссщс стс с стс ст сщс стст с с стстс ст сщс стс с стс ст сщс стст с стс ст сщс с с стс ст сщс стст с стс ст сщс ст с стс с стст сщссщс стс с с стс стс сщс с ст с стс с стст сщссщс

с с с с с стс ст с стс ст с стс ст с стс ст стс с стс с стс с стс с с с с ст с с ст с с ст с с ст стс с с стстс с с стстс с с стстс с с стстс с с с с с стст с стст с стст с стст стстс ст с стстст с с с с с стст с ст с стст с с с с с стст с ст с стст с с с с с стст с ст с стст с с с с стст с ст с стст с стс ст с с с с стс

ст ссщс с ссщ стст сщссщ стс ст ссщс стст сщссщ с сщс стстст с сщссщс стс с сщс с ст сщссщс

ст щ с щ стст щщ стс ст щ стст щщ с щ стстст с щщ стс с щ с ст щщ

ст сшсчсш с сшсч стст счсшсч стс ст сшсчсш стст счсшсч с счсш стстст с счсшсчсш стс с счсш с ст счсшсчсш

ст сшсчсш с сшсч стст счсшсч стс ст сшсчсш стст счсшсч с счсш стстст с счсшсчсш стс с счсш с ст счсшсчсш

ст шчш с шч стст чшч стс ст шчш стст чшч с чш стстст с чшчш стс с чш с ст чшчш

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стст с

с ст стст с ст стст с ст стст с ст стст с ст с стст с с стс ст с с с стс ст с с с стст

т т т т тт

с ст стст с ст стст с ст стст с ст стст с ст с стст

с ст с с стст сусфсусф стс с с сусф стст с с стст сусф стс с с сусф стст с с стст сусфсусф стстс стс с сусф ст с с с ст сусфсусф стс с стс с сусф с ст с с с ст сусфсусф

с с с с стс ст с стс ст с стс ст с с с с с с с стст с стст с стст с с с с с с с с стстс с стстс с стстс с стстс с с с с стст с стст с стст стстс ст с стстст с с с с с стст с ст с стст с с с с с стст с ст с стст с с с с стст с ст с стст с стс ст с с с с стс

ст сфсусф с сфсу стст сусфсу стс ст сфсусф стст сусфсу с сусф стстст с сусфсусф стс с сусф с ст сусфсусф

ст фуф сцсхсц с фу сцсх стст уфу схсцсх стс ст фуф сцсхсц стст уфу схсцсх с уф схсц стстст с уфуф схсцсхсц стс с уф схсц с ст уфуф схсцсхсц

ст сцсхсц с сцсх стст схсцсх стс ст сцсхсц стст схсцсх с схсц стстст с схсцсхсц стс с схсц с ст схсцсхсц

ст цхц с цх стст хцх стс ст цхц стст хцх с хц стстст с хцхц стс с хц с ст хцхц

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

с с стс ст с с с стс ст с с с стст с с с с стс ст с с с с стст

с ст стст с ст стст с ст стст с ст стст с ст с стст

т тт т тт т тт т тт т с тт с с стс ст с с с стс ст с с с стст т т тт

dir

s

f

C

sub

1

d

1

dir

w0 w

d
1

sub

C

2

d

2

w0 w

1

sub

C

d

3

3

sub

2

s

5


us that the second approach is more efficient unless the length of transformed vectors is not big. Thereafter, we incorp orated double transp osition in direct Helmholtz equation solver. To apply iterative solver we made use of the PIM 2.2 package [8], which we found to b e useful. It includes a numb er of iterative metho ds for symmetric and nonsymmetric systems. As our matrix is nonsymmetric, we inclined to the Bi-CGSTAB metho d with the stopping criteria ||rk || < ||b||, where rk is the residual, b is the right hand side, and = 10-2 . The only thing supplied by us to PIM was a matrix-vector multiplication subroutine, that invokes the b oundary exchange communication routine. Table 1 contains an example of b enchmarking p erformed with the Eulerian dynamical core on MBC1000M computing system. One can see, that the semi-implicit scheme with the direct solver more scalable than that with the iterative solver and faster than the explicite scheme, at given resolution.

сфсусусф сфсусф сфсусу сфсусусф сфсусф сусфсу сфсусусф сфсусф сфсусу сфсусусф сфсусф сфсусу сфсусусф сфсусф сфсусу сфсусусф сфсусф сфсу
1

сусф уф ст сусф уф с стст сфсусф фуф стст сусфсу уфу ст сусф уф стс сусф уф с ст сфсусф фуф стстст сусфсу уфу ст сусф уф стст сусф уф стст сфсусф фуф стст сусфсу уфу стст сусф уф стст сусф уф стст сфсусф фуф стс сусфсу уфу стст сусф уф с сусф уф с стст сфсусф фуф стст сусфсу уфу ст сусф уф стс сусф уф с ст сфсусф фуф стст су сусфсу уфу стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стстс с с с стс ст с с с стс ст стстс с с с с стс ст

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стст с ст стст стс стст с с ст

стс ст с с с стст с с стст стс с с стст с с стст стс стстс с ст с с с стст с стс стст с с с стст с стс стст с с с стст с с с ст стстст с с стс с стстст с стстс с с с стстст с стстс с с с стст с с стст стс с с стст с с стст стс стстс с ст с с с стст ст с стс с стст стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стстс с с с стс ст с с с стс ст стстс с с с с стс

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стстс с с с стс ст с с с стс ст стстс с с с с с с с стс

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стст с ст стст стс стстстст

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стст с ст стст стс стст с с ст

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стстс с с с стс ст с с с стс ст стс стс с с стс
d
2

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стстст стс с ст стст с с стст

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стс стс стстс с с с стстс с с с с стс ст с с

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стстст стс с ст стст с стст

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стс стс стстс с с с стстс с с с с стс ст с с стс

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стст с ст стст стс стст с с ст

стс ст с с с стст с с стст стс с с стст с с стст стс стстс с ст с с с стст с стс стст с с с стст с стс стст с с с стст с с с ст стстст с с стс с стстст с стстс с с с стстст с стстс с с с стст с с стст стс с с стст с с стст стс стстс с ст с с с стст ст с стс с стст

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стс стс стстс с с с стстс с с с с стс ст с с

Figure 3. Transp osed data.

6

3. OPENMP IMPLEMENTATION OF THE INM OCEAN MODEL

The INM Ocean mo del [9] is based essentially on the same governing equations as the atmospheric mo del, with the difference that the salinity is incorp orated as a prognostic variable. As a comp onent of the coupled mo del [10] it p erforms ab out 10% of overall computations. The rigid lid condition at the surface allows to intro duce the barotropic stream function. The numerical implementation makes use of splitting into physical processes, as well as along spatial co ordinates, facilitating application of the efficient implicit

схсх хх сх х схсх хх схсх хх сх х схсх хх схсх хх сх х схсх хх схсх хх сх х схсх хх сх хх

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стстст стс с ст стст с стст

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стс стс стстс с с с стстс с с с с стс ст с с стс

стст с ст стст с ст стст стс стст стст стст стст стст стст стст с ст стстст стс с ст стст с ст стст стс стст с с ст

стс ст с с с стст с с стст стс с с стст с с стст стс стстс с ст с с с стст с стс стст с с с стст с стс стст с с с стст с с с ст стстст с с стс с стстст с стстс с с с стстст с стстс с с с стст с с стст стс с с стст с с стст стс стстс с ст с с с стст ст с стс с стст

стс с с с с стс ст с с с стс ст стстс с с с с стс с с с с стс с с с с с с ст с стс с с с стстс с с с стс стс стстс с с с стстс с с с с с стст
0

т т т т т тт т т т т тт т т тт тт тт

тт т тт т тт т тт тт тт тт тт тт тт т ттт т т ттт т т тт тт

сх цц схсх цц схсх ц сх цц схсх цц схсх ц сх цц схсх цц схсх ц сх цц схсх цц схсх ц сх цц схсх цц схсх ц сх цц схсх цц т сх цсх

сцсцсхсх сцсц сцсцсх сцсц сцсхсх сц сцсцсхсх сцсц сцсцсх сцсц сцсхсх сц сцсцсхсх сцсц сцсцсх сцсц сцсхсх сц сцсцсхсх сцсц сцсцсх сцсц сцсхсх сц сцсцсхсх сцсц сцсцсх сцсц сцсхсх сц сцсцсхсх сцсц сцсцсхсх сцсцсх сцсх сц

dir

s

dir
сх х схсх хх схсх хх сх х схсх хх

w

d

sub

g

D

D D

3

1

2

w

sub

sub

sub

2

3

s

s

1

s


7 Table 1 CPU time (ab ove) and sp eed-up (b elow) of one day forecast (in sec) of the Eulerian dynamical core at resolution 256 ╜ 128 ╜ 16 on a given numb er of pro cessors of MBC1000M computing system. EXPL explicit scheme, GAUSPH semi-implicit scheme with the direct solver, ITER semi-implicit scheme with the iterative solver. # of pro cs 1 2 4 8 16 32 64 128 256 EXPL 602. 367. 190. 103. 54. 43. 24. 20. 18. 1.0 1.6 3.2 5.8 11.0 13.9 25.4 30.6 33.5 GAUSPH 233. 156. 100. 69. 53. 34. 17. 11. 10. 1.0 1.5 2.3 3.4 4.4 6.9 13.8 20.7 23.4 ITER 306. 178. 99. 53. 35. 35. 31. 31. 32. 1.0 1.7 3.1 5.8 8.5 8.6 9.7 9.9 9.5

algorithms. Intro duction of the relative depth co ordinate makes the computational domain cylindrical. Though the horizontal grid is uniform in spherical co ordinates and fields are stored in multidimensional arrays, which dimensions corresp onding to spatial co ordinates, only wet p oints are computed and the pro cessed domain has rather complex form in horizontal plane. The shared memory approach is attractive, mainly, b ecause it do es not involve such profound restructuring of co de, as MPI do es. Secondly, it is suitable for the Ocean mo del, b ecause the dynamic scheduling, if enabled, yields go o d load balancing even for an ugly domain. The mo del parallelization was accomplished by supplying the PARALLEL DO directive to ma jor lo ops. Red-black ordering of unknowns was intro duced in the SOR metho d that calculates the barotropic stream function. Ultimately, sp eedup 2.6 was reached on a 4-pro cessor Itanium I I system (Table 2). A distributed memory version of the mo del, obtained by A.S. Rusakov [11] outp erforms the Op enMP version. On 8 CPUs of MBC1000M it revealed sp eed-up of 4.9 which, however, was estimated without writing some output files. Table 2 CPU time of the Op enMP parallel # of CPUs 1 IBM Regatta, p670 16m47s Intel Itanium I I 15m46s Sun UltraSparc I I 1h38m24s

version of 2 11m28s 9m46s 1h5m48s

the INM Ocean mo del one year runs. 4 8 8m39s 7m16s 5m58s 48m48s

ACKNOWLEDGMENTS We would like to thank the Joint Sup ercomputer Center for providing there computing facilities which include the first Russian Teraflop machine MBC1000M. Also supp ort


8 of Russian-Indian Centre for Advanced Computing Research is appreciated. This work would not b e p ossible without assistance of E.M. Volo din, N.A. Diansky, and M.A. Tolstykh, all from INM RAS. REFERENCES 1. V. Gloukhov, Parallel implementation of the INM atmospheric general circulation mo del on distributed memory multipro cessors, Lecture Notes in Computer Science 2329 (2002) 753э762. 2. I.M. Held and M.J. Suarez, A prop osal for the intercomparison of the dynamical cores of atmospheric general circulation mo dels, Bull. Am. Meteorol. So c. 73 (1994) 1825э1830. 3. B. Ro driguez, L. Hart, and T. Henderson, Parallelizing op erational weather forecast mo dels for p ortable and fast execution, Journal of Parallel and Distributed Computing 37 (1966) 159э177. 4. M. Snir, S. Otto, S. Huss-Lederman, D. Walker, J. Dongarra, MPI: The Complete Reference, The MIT Press, Cambridge, MA, 1998 5. V.A. Alexeev, E.M. Volo din, V.Ya. Galin, V.P. Dymnikov, V.N. Lykossov, Simulation of the present-day climate using the INM RAS atmospheric mo del. The description of the mo del A5421 (1997 version) and the results of exp eriments on AMIP I I program. Institute of Numerical Mathematics RAS, Moscow (1998) (reg. VINITI 03.07.98 No. 2086-B98) 6. G.I. Marchuk, V.P. Dymnikov, V.B. Zalesny, V.N. Lykossov, V.Ya. Galin, Mathematical Mo deling of the Atmosphere and Ocean, Gidrometeoizdat , Leningrad, 1984 (in Russian) 7. A. Arakawa, V.R. Lamb, Computational design of the basic dynamical pro cesses of the UCLA general circulation mo del, Metho ds Comput. Phys. 17 (1977) 173э265. 8. R.D. Cunha and T. Hopkins, PIM 2.2 The parallel iterative metho ds packages for systems of linear equations. User's guide. Internal Rep ort 2э94 of the Computing Lab oratory, UKC. 9. N.A. Diansky, A.V. Bagno, and V.B. Zalesny, Sigma mo del of global o cean circulation and its sensitivity to variations in wind stress, Izvestiya, Atmospheric and Oceanic Physics, 38 No. 4 (2002). 10. N.A. Diansky, E.M. Volo din, Repro ducing the present day climate using a coupled atmosphere-o cean general circulation mo del. Izvestia, Atmospheric and o ceanic physics 38 No. 6 (2002). 11. A.S. Rusakov, N.A. Diansky, Parallel o cean general circulation mo del for distributed memory computer systems, Parallel CFD 2003, Bo ok of abstracts.