h14520
s 00435/00000/00000
d D 1.1 95/02/15 19:13:59 mnataraj 1 0
c date and time created 95/02/15 19:13:59 by mnataraj
e
u
U
f e 0
t
T
I 1
\documentstyle [12pt, sane]{article}
\setlength{\textheight}{8.7in}
\setlength{\topmargin}{-1cm}
\setlength{\oddsidemargin}{+7mm}
\setlength{\textwidth}{6.2in}
\pagestyle{fancy}
\input{psfig}

\begin{document}
\title{Case Studies on Parallelization of NASA's Weather Code:\\
       Results on Interpolation Techniques}
\author{{\it Director} : Geoffrey Fox \\
{\it Project Leader} : Miloje Makivic \\
{\it Research Assistant} : Manu Natarajan \\ \\
Northeast Parallel Architectures Center \\
	Syracuse University, 111, College Place \\
	Syracuse, NY, 13244-4100 }


\date{}
\maketitle
\clearpage
\tableofcontents

\newpage
\noindent
\section{Abstract}

The purpose of this project is to investigate suitable domain decompositions
that lead to efficient parallelizations of the interpolation techniques
used in the Semi-Lagrangian advection scheme.  In this report we present
the results of our investigation using the lasagnia decomposition of the
input domain of the bilinear and bicubic interpolations used in NASA's 
2-D Shallow Water Model.  

Our approach to the usage of the lasagnia decomposition reduces the amount
of communication needed for either the bilinear or bicubic interpolations.
The criteria for goodness of the algorithm is based on the ratio of the
useful work done to the amount of communication required.

\section{Introduction}

Weather forecasts are essential for the planning of many activities. They need
to be accurate and available in time.  The basis of weather forecasting are
Numerical Weather Prediction (NWP) models which assimilate data and 
predict future conditions of atmosphere. Given the requirement that forecasts
be available in real-time, the numerical models must be efficient.  There is
also considerable interest in building quality data-sets or model predictions
of weather conditions for research purposes and as case studies.  These involve
processing of voluminous data using the numerical models in an efficient
manner.

One chronic problem with NWP models is that the longest time step that can be
chosen is decided by stability rather than accuracy.   This is alleviated to
a considerable degree in the Semi-Lagrangian Transport Model using a 
semi-implicit time differencing scheme \cite{rob-81,rob-82}. This was 
demonstrated for a limited area shallow water model (one layer, homogeneous
 fluid).  An
excellent review of Semi-Lagrangian integration schemes for atmospheric models 
is treated by \cite{stan-91}.  

Later work addressed the deficiencies associated with polar singularity 
in the latitude-longitude co-ordinate system \cite{mc-89}.  NASA
implemented a modified model based on this paper 

Another notable issue is that for accuracy to be improved, among other
things higher resolutions are required. This calls for efficient
implementations on powerful machines.  Effort has been underway to implement 
these models on parallel computers \footnote{Weather Modelling project 
undertaken at NASA Goddard Laboratory, Maryland}. 


In order to establish the most suitable paradigm, certain modules of NASA's
serial implementation were identified as test-cases. 
Initial studies at Syracuse University 
and NASA have indicated that the interpolation techniques deployed occupied a 
total running time of 25\% of running time in the serial code.  It was also 
observed that the usage of the construct "forall" did not have
a good speed-up in the data parallel prototype on the CM-5.  This can be 
attributed to the irregular, non-local pattern of
data accesses.  In addition the decomposition was carried out in the 
North-South axis as shown in Fig. 1; as one expects winds to be mostly in the 
East-West direction.
This lead to a very poor utilization of the nodes when there was 
a small grid size and  one of the grid axes 
declared as serial or local to a
node.  In such a case the data is distributed among only a few nodes along 
the news 
or parallel axis leaving remaining processors with no work \cite{cm5-ref}
The default news:news (square or ravioli) decomposition 
showed better results.  This is because the run time library is able to 
ensure a better utilzation of nodes with the default decomposition. But 
this decompostion does not take advantage the fact that the winds are of a
East-West nature.


Hence, the message passing model was proposed as an alternate
parallelization technique
as it offers a greater flexibility in choice of domain decomposition.  One 
does not need to depend on the vagaries of run-time libraries or systems
to decide the utilization of nodes.  Scalable decompositions may be 
developed.  Additionally, the decomposition can be adapted to ensure that 
the pole regions are treated differently than regions near the equator.  
This is because, at the poles the winds move around the pole (entire latitude).
This is another disadvantage using data parallel 
techniques which permit only simple decompositions.  Hence, nearest neighbor
communication is not possible at the pole for the default square decomposition
on CM-5 when using a latitude-longitude grid.

A lasagnia decomposition along the N-S axis was adopted as the 
first test case as seen in Fig. 2.  This was postulated to have mostly 
nearest neighbor 
communication in the N-S direction.  The polar region problem is
avoided as the lasagnia strips are in E-W direction. The algorithm 
uses a message-passing version of the data-parallel construct "cshift".
The maximum number of "cshifts" is determined at run-time at every time
-step for each component.  The necessary data from the circulated 
wind "strips" is collected and a local interpolation is performed at
every time step.  Two strategies and their merits are discussed in the
next section.

\subsection{Metrics}

We propose some measurements for assessing the goodness of the algorithm
as well as giving direction to future work :

\begin{itemize}
 \item Speed-up of interpolation w.r.t the serial code
 \item Ratio of calculation to communication.
 \item The maximum displacement in terms of nodal distance in the N-S
       direction.
 \item Accuracy of calculations 
\end{itemize}

The first criterion shows how well the algorithm performs with respect to 
the serial code as the problem size grows.  The numerical value is of not
much interest.  But, a very low figure would indicate that the algorithm
was ill-conceived   \footnote{A more appropriate measurement would be the 
speedup 
obtained by using varying number of processors on a parallel machine. }.
The second criteria
indicates how well the algorithm will scale when a global
 3-D Semi-lagrangian model is used.  The third may give a clue as to what could
be a suitable ratio of the partition dimensions in more complex decompositions.
It would also verify if the assumption that the winds are mostly E-W is true.
The last criterion is important as it indicates the presence of errors other
than roundoff.  This done by finding out if the difference between the 
interpolated values obtained from both the serial and parallel implementations
are homogenous for a given resolution.

\section{Implementation}

The project was implemented on CM-5 using the native message passing library
\cite{cmmd-ref}. For each component, the partitions of wind 
components are circulated among the
nodes in both North and South directions.  The amount of circulation depending
on the maximum displacement obtained by a global reduction of local (to a node) maximum displacement.  The idea is to be able to select the required data for 
 the component which will allow each node to perform the interpolation locally.

\subsection{Pre-Determination of Degree of Circulation}

If the assumption that the winds in the North-South direction mostly result in
nearest neighbor communicationis valid , then the following strategy may be 
used.
For any given global indice (ii, jj), we determine the node it belongs to.
The inter-nodal distance is then simply the difference between the owner node
and user node.  But, this distance has a direction, i.e. the ``owner'' node 
may be
reached in the North or South Direction [Fig. 5].  Hence, this is converted to
 the
inter-nodal distance between the two nodes when going in the North direction.
Moreover, if the inter-nodal distance exceeds half the number of nodes, the
shorter path is in the Southern direction.  For, either interpolation there are
two extreme points in the North-South direction [Fig. 3,4]. The inter-nodal 
distance is computed for both the points separately, for all the indices on a
node.  The maximum inter-nodal distance in a direction (both North and South
are evaluated) is then found by a global reduction. The two maxima in the 
North-South direction govern the circulation.

\subsection{Selection of Data Points}

There are at least two strategies to selecting the data points. The two differ
in how the data is retrieved from the circulated wind strips.

In the first approach, the data is retrieved into one of 'n' (4 for in case of
bilinear or 16 in case of bicubic) buffers. These 'n' buffers correspond to the
number of points used in the interpolation [Fig. 3,4].  In order for this to 
be convenient an overlap region (3 columns in case of bilinear) [Fig. 6] 
has to be constructed in the
northern direction.  This convenience comes at the cost of pre-fetching the
data before selection from the northern neighbor.  This communication increases
when the width of an individual strip is smaller than the overlap region.
 [Fig. 7]  The interpolation is done by indexing the 'n' buffers to obtain the
'n' points.  The bilinear interpolation was implemented by using this approach.

In the second approach, the creation of the overlap region is dispensed with.
The circulated wind strip is appended as shown in Fig. 8.  However, additional
book-keeping needs to be done in order to the relative position of the
partition in the segment and for generation of local indices.  The 
interpolation is then carried out locally. The 'n' points are obtained by the
use of the local indices.  The bicubic interpolation was done by using this
approach.

\section {Performance}
	The timings were obtained on CM-5 on CMOST 7.2 on a 32 node
partition.  The serial machine used for comparision were Sun workstations.
The performance is discussed with respect to the four criteria 
discussed in the Introduction.  They are as follows :

\begin{itemize}

 \item Speed-up of interpolation (with respect to serial code)
 \item Ratio of calculation to communication.
 \item The maximum displacement in terms of nodal distance in the 
       N-S direction.
 \item Accuracy of calculations 

\end{itemize}

The first table shows the non-linear nature of the communication overhead
versus message size.  The overhead decreases as message size increases from 
0 to 8 bytes.  Somewhere between a size of 8 to 40 bytes, the overhead
starts increasing again.

	\paragraph\noindent{
	{\bf Table 1 : Communciation Time vs Message Size} \newline
	\begin {tabular}{||l|l|l|l||} \hline
	Message Size & Communication \\ 
	(bytes)      & Time (ms) \\ \hline
	0 & 0.222 \\ 
        4 & 0.189 \\
        8 & 0.135 \\ 
        40 & 0.146 \\
        240 & 0.170 \\
        400 & 0.165 \\
        2000 & 0.237 \\
        4000 & 0.340 \\
        20000 & 1.378 \\
        40000 & 2.469 \\  \hline
	\end{tabular} \\\\


\subsection {Bicubic Interpolation}


	{\bf Table 2 : Speedup with respect to serial code} \newline
	\begin {tabular}{||l|l|l|l||} \hline
	Resolution &   Total Time    & Total Time & Speedup \\ 
	  (I,J)     &  CM-5 (tp ms)  &  Sun (ts ms) &   (ts / tp) \\ \hline
	(8,64)    & 22.60 &    24.64  & 1.09 \\ 
        (64,64)   & 49.48 &   304.16 & 6.15 \\
	(128,128) & 132.58 &   967.32 & 7.30 \\
        (256,512) & 890.84 &  7716.23 & 8.66 \\ \hline
	\end{tabular}\\\\

	{\bf Table 3 : Performance of Bicubic Interpolation} \newline
	\begin {tabular}{||l|l|l|l||} \hline
	Resolution & Total Commn. & Calculation & Ratio of Calc. Time \\ 
	  (I,J)     & Time (ms)  &  Time (ms) &  to Commn. Time \\ \hline
	(8,64) & 4.02 & 18.57 & 4.61 \\  
        (64,64) & 6.49 & 42.99 & 6.63 \\
	(128,128) & 9.77 & 122.80 & 12.56 \\
        (256,512) & 23.72 & 867.12 & 36.54 \\ \hline
	\end{tabular}\\\\



   The second approach to selecting data points is adopted for the bicubic
   interpolation.
   The speedup with respect to the linear code ranges from a little over 1
   to about 8.66 [Table 2] \footnote { The timings of the parallel code in 
   both the
   bicubic and bilinear interpolations do not include the cost of performing
   the lasagnia decomposition. The cost of decomposition should be amortized 
   in calculating the speedup of the entire parallel code 
   with respect to the serial code and not just the
   speedup of either of the interpolation schemes. }.  From this, we can infer
   that the message passing paradigm is capable of achieving a reasonable
   speedup (at higher resolutions) as far as the bicubic interpolation is
   concerned.  

The ratio of calculation to communication increases from about 4.6 to 36.5
as the resolution increases [Table 3].  This indicates that 
the communication may not be a limiting factor but the 
amount of calculations may have to optimized.  

   The inter-nodal distance seems to be fairly static for both components
   (north, south directions inclusive). On no occasion does this 
   exceed 3 "cshifts" in the Northern direction or 1 "cshift" in the 
   Southern.  Typical values are 1 "cshift" in either direction.  
   This is close to the assumption that the
   winds generate mostly nearest neighbor communication.

  An accuracy of upto (but not including) the fifth decimal place was
   obtained between the parallel and serial implementations.
   
\pagebreak
\subsection {Bilinear Interpolation}


%	{\bf Table 4 : Speedup with respect to serial code - 
%2nd algorithm} \newline
%	\begin {tabular}{||l|l|l|l|l||} \hline
%	Resolution &   Total Time    & Total Time & Speedup \\ 
%	  (I,J)   &  CM-5 (tp ms)  &  Sun (ts ms) &   (ts / tp) \\ \hline
%	(8,64)    & 127.45  & 15.24   & 0.119 \\ 
%        (64,64)   & 149.45  & 89.05   & 0.596 \\
%	(128,128) & 207.19  & 273.10  & 1.318 \\
%        (256,512) & 763.82  & 2187.41 & 2.863 \\ \hline
%	\end{tabular}\\\\
%
%	{\bf Table 5 : Performance of Bilinear Interpolation - 2nd
%algorithm} \newline
%	\begin {tabular}{||l|l|l|l||} \hline
%	Resolution & Total Commn. & Calculation & Ratio of Calc. Time \\ 
%	  (I,J)     & Time (ms)  &  Time (ms) &  to Commn. Time \\ \hline
%        (8,64)  & 4.29   & 123.16 & 28.68 \\
%        (64,64) & 6.28   & 143.16 & 22.77 \\
%	(128,128) & 8.34 &  198.84 & 23.83 \\
%        (256,512) & 23.86 & 739.96 & 31.01 \\ \hline
%	\end{tabular}\\\\



	{\bf Table 4 : Speedup with respect to serial code } \newline
%- 1st algorithm} \newline
	\begin {tabular}{||l|l|l|l|l||} \hline
	Resolution &   Total Time    & Total Time & Speedup \\ 
	  (I,J)   &  CM-5 (tp ms)  &  Sun (ts ms) &   (ts / tp) \\ \hline
	(8,64)    & 119.17  & 15.24   & 0.127 \\ 
        (64,64)   & 140.83  & 89.05   & 0.632 \\
	(128,128) & 180.64  & 273.10  & 1.511 \\
        (256,512) & 675.36  & 2187.41 & 3.239 \\ \hline
	\end{tabular}\\\\

	{\bf Table 5 : Performance of Bilinear Interpolation } \newline
% - 1st algorithm} \newline
	\begin {tabular}{||l|l|l|l||} \hline
	Resolution & Total Commn. & Calculation & Ratio of Calc. Time \\ 
	  (I,J)     & Time (ms)  &  Time (ms) &  to Commn. Time \\ \hline
        (8,64)  & 3.17   & 116.00 & 36.56 \\
        (64,64) & 5.22   & 140.82 & 25.96 \\
	(128,128) & 7.47 &  173.16 & 23.17 \\
        (256,512) & 34.06 & 641.29 & 18.82 \\ \hline
	\end{tabular}\\\\

   The first approach to selecting data points was utilised for the 
   bilinear interpolation.  The speedup with respect to the linear code varies
   from less than 1 to about 3.2 [Table 4]. 
    The reason for the lower speedup of the
   bilinear interpolation is largely due to the fact the division of 
   the computation among the processors was largely offset by the
   cost of parallelization.  As the computations are more in bicubic
   interpolation, the ``division of work'' paid better dividends \footnote{If
   the additions are ignored (which are also more in the bicubic 
   interpolation), there are 32 floating point multiplications in the bicubic 
   interpolation versus 6 for the bilinear for every point of interpolation. }.
%But   as the figures show, the message passing version performs better than
%   the data parallel version at higher resolutions.

  The ratio of calculation with communication decreases as the resolution 
   increases [Table 5].  One would have expected the ratio to have been 
   low at lower resolutions, due to the additional communication overhead 
   in the construction of the northern buffer [Fig. 5].  This needs to be 
   investigated further.

   The inter-nodal distance seems to be fairly static for both components
   (north, south directions inclusive) with neither direction
    requiring more one cshift.  This verifies the assumption that the
    winds generate mostly nearest neighbor communication.

   An accuracy of upto (but not including) the fifth decimal place was
   obtained between the parallel and serial implementations.


\section {Conclusions}

	 A reasonable speedup maybe achieved by the message passing version.
The cost of parallelization (which includes both the communication overhead
 and any computational overhead due to the parallelization) 
should be relatively low and not just the communication overhead for any
message passing algorithm to work.

	It appears that communication overhead may not be as high or as 
menacing as was assumed, for either the bilinear or bicubic interpolations.
The pre-determination of the maximum number of "cshifts" plays an important 
role in limiting the communication.  Additionally, the non-linear communication
overhead with varying buffer size on the CM-5 [Table 1] influences the results
for intermediate grid sizes.

\section {Further Work}

	In the short term further testing needs to be done to verify
the accuracy and stability of the parallel code.  The scaling of the
performance of the parallel code with the number of processors needs
to be observed.

	In the long term more, more scalable decompositions such as
the igloo decomposition have to be studied.  A different approach will
be needed.  The "cshift" operation will no longer be applicable. A key
strategy to make the igloo decomposition feasible would be to find a
near-optimal ratio of the partition dimensions.  This is has to be
ascertained by initial experiments.

\newpage
\begin{thebibliography}{9}

\bibitem{cm5-ref}
Thinking Machines Co., 1993: CM Fortran Reference Guide {\bf Ch. 12}

\bibitem{cmmd-ref}
Thinking Machines Co., 1993: CMMD User's Guide and Reference Manual

\bibitem{mc-89}
McDonald, A.; Bates, J.R., 1989: Semi-lagrangian Integration of a Gridpoint
Shallow Water Model on the Sphere. {\em J. American Meteor. Soc.}, {\bf 117},
130-137

\bibitem{rob-81}
Robert, A.,  1981: A stable numerical integration scheme for the primitive
meteorological equations. {\em Atmos. Ocean}, {\bf 19}, 35-46

\bibitem{rob-82}
Robert, A.,  1982: A semi-Lagrangian and semi-implicit numerical integration
scheme for the primitive meteorological equations. {\em J. Meteor. Soc. Japan},
{\bf 60}, 319-324

\bibitem{stan-91}
Staniforth, A.; Cot'e, J., 1991: Semi-Lagrangian Integration Schemes for
Atmospheric Models - A Review. {\em Monthly Weather Review}, {\bf 119}, 
2206-2203


\end{thebibliography}
\end{document}
E 1
