%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[12pt]{article}
\setlength\topmargin{-0.0cm}
\setlength\textheight{22.0cm}
\setlength\textwidth{15.0cm}
\addtolength\evensidemargin{-1.cm}
\addtolength\oddsidemargin{-0.5cm}

\usepackage{psfig}
\bibliographystyle{unsrt}

\title{CSE 557 \\ Spring 2000\\
Parallel Implementation of Multigrid Methods for Solving Poisson's Equation}
\author{Anirudh Modi and Anupam Sharma}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\myfigure}[5]
        {
         \begin{figure}
                \centerline{\psfig
                        {figure=#1,angle=0,height=#3,width=#4}}
                \caption{#5}
                \label{fig:#2}
         \end{figure}
        }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\myfigureB}[8]
        {
         \begin{figure}
                \centerline{\psfig
                        {figure=#1,angle=0,height=#5,width=#6}}
                \caption{#7}
                \label{fig:#3}
         \vspace{0.5cm}
                \centerline{\psfig
                        {figure=#2,angle=0,height=#5,width=#6}}
                \caption{#8}
                \label{fig:#4}
         \end{figure}
        }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\maketitle

\abstract
This work attempts to solve the Poisson's Equation using 
accelerator techniques like 
{\em Multigrid} along with several relaxation schemes
such as {\em Jacobi}, {\em Gauss-Siedel} and {\em Successive Over
Relaxation} (SOR). The implementation is parallelized using the
parallel grid communication library, {\em GridComm}. The convergence
rates for the various schemes are studied, and their parallel efficiencies
for various 1D and 2D partitioning and grid sizes are compared.

\section{Introduction}
Analytical solution to physical problems is highly desirable, but is 
almost impossible to obtain for most practical problems. Numerical
techniques provide a convenient way to obtain the solution for very
complex problems by simulating the physics on a computer. Simulation
of any problem usually employs the standard relaxation techniques:
 {\em Jacobi}, {\em Gauss-Siedel}, {\em SOR}, etc. Simulation of most of the
engineering problems are very demanding in terms of memory requirement
and computations. The problems can be so huge and computationally
intensive that the simulations may take unacceptable amount of time.
This presents the need for {\em accelerating} the solution
process. A fast computer can be understood as an {\em accelerator}.
Parallel implementation of a problem is another way to accelerate
the solution process. A number of {\em accelerating} techniques
are available. Among these, the most popular technique is the {\em
Multigrid} technique.

Multigrid employs more than one interacting grids to solve the problem.
The inefficiency of the single-grid solvers is owing to the fact that it
can only efficiently resolve the high-frequency error components, while
the low frequency components take a long time to be resolved. Multigrid,
on the other hand, uses coarse grids to relax the low-frequency error
components. The low-frequency error components can be represented as
high-frequency error components on the coarse grids. These can then be
efficiently relaxed and interpolated back to the fine grid to update
the solution. Implementation of multigrid involves restriction and
interpolation operators to communicate between different grid levels.
The standard relaxation schemes are used to relax the error (solution in
case of fine grid) at every grid level.

This work deals with solving the Poisson's equation in parallel using
Multigrid technique. The Poisson's equation in 2D Cartesian co-ordinate
system is written as:

\begin{equation}
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f
\label{eq:poissons}
\end{equation}

\section{Relaxation Schemes}
The matrix $A$ in the problem $Au=b$ can be re-written as $A=D+L+U$ where $D$
is the diagonal matrix, $L$ is the lower diagonal and $U$ is the upper
diagonal matrix. The standard relaxation schemes used for most numerical
techniques are {\em Jacobi}, {\em Gauss-Siedel} and {\em Successive Over
Relaxation} (SOR). These can be mathematically represented by
Eqns.~\ref{eq:jacobi},~\ref{eq:GS} and~\ref{eq:SOR}.

\begin{equation}
u^{n+1} = -D^{-1}(L+U)u^n + D^{-1}b
\label{eq:jacobi}
\end{equation}

\begin{equation}
u^{n+1} = -(D+L)^{-1}u^{n} + (D+L)^{-1}b 
\label{eq:GS}
\end{equation}

\begin{equation}
u^{n+1} = (D+\omega L)^{-1}\left [(1-\omega)D -\omega U\right ] u^{n} + (D+\omega L)^{-1}\omega b 
\label{eq:SOR}
\end{equation}

Gauss-Siedel relaxation scheme is exactly the same as the Jacobi scheme
except that the Gauss-Siedel utilizes the updated values as soon as
they are available. On the other hand, the Jacobi scheme only uses
the updated values after one full iteration. The Gauss Siedel scheme
therefore converges faster when compared to the Jacobi scheme.

For parallel implementation of Gauss-Siedel, Eq.~\ref{eq:GS} cannot be
used, as all the updated values from the neighboring processors has to
available before $u$ is updated. This problem can be solved by using
{\em Red-Black ordering}. The points for which the sum of the indices ($i+j$)
is {\em even} are called {\em red} cells, and the points for which the sum is 
{\em odd} are called {\em black} cells. The red cells are updated first
using the values of black cells (calculated in the previous iteration),
and then the black cells are updated using the updated values of the red
cells. Message passing between the neighbors is done after updating the
red cells. Red-Black ordering is able to solve the problem because only
the black cells are required to update the red cells and vice-versa.
This may not always be the case, especially when one is using {\em higher
order accurate} methods.

The idea behind SOR is to {\em more than update} the solution in one
iteration. This is done by introducing a relaxation parameter, $\omega$.
The correction to the solution, obtained using Gauss-Siedel approach,
is multiplied with $\omega$ and added to the solution. If $\omega$ is
greater than $1$, the scheme is called over-relaxation; and if it is less
than $1$, it is under-relaxation. The case $\omega=1$, is the Gauss-Siedel
relaxation. The convergence rate of an SOR scheme critically depends on
the choice of $\omega$. For the Poisson's problem, the optimum value of
$\omega$ is equal to $1.5$.

\section{Multigrid as an Accelerator}
For implementation of multigrid one needs toe define the ``restriction''
and the ``interpolation'' operators. The interpolation operator is
represented by I$_{2h}^{h}$ and the restriction operator is represented
by I$_{h}^{2h}$. The restriction operator is used to restrict the
residue on the fine grid to the coarse grid, while the interpolation
operator is used to interpolate the correction to the solution from the
coarse grid to the fine grid. Usually, these operators are simple
weighted-averages.
The algorithm for implementing multigrid can be described as follows.

\begin{enumerate}
    \item Pre-smoothing ($\nu_1$ cycles of relaxation on the fine grid)
    \item Coarse grid correction
    \item Post smoothing ($\nu_2$ cycles on the fine grid)
\end{enumerate}

The coarse grid correction is the step which involves the communication
between coarse grids using the interpolation and restriction operators.
The solution obtained by relaxation on the fine grid is used to find the
residue on the fine grid.

\begin{equation}
r^{h}=b^{h}-A^{h}u^{h}
\end{equation}
The residue is then restricted to the coarse grid and a different problem is solved at the coarse grid.

\begin{equation}
A^{2h}e^{2h}=I_{h}^{2h}r^{h}
\end{equation}

The error obtained at the coarse grid $e^{2h}$ is then interpolated to
the coarse grid and the solution is updated.

\begin{equation}
u^{h} = u^{h} + I_{2h}^{h}e^{2h}
\end{equation}
This completes the coarse grid correction step for a two-grid multigrid.
For multi grid levels, the coarse grid correction step is recursively
used to correct for error at every coarse grid level.

\section{Multigrid in Parallel}
The parallel implementation of the multigrid technique was a trivial
extension of the serial algorithm, thanks to the {\em GridComm} library
written previously. Since the {\em GridComm} library takes care of all
the grid partitioning and communication, the main program simply needs
to initialize and call routines supplied by the library interface.
All the matrices in the original serial algorithm were initialized as
parallel grids of type {\tt GCgrid}, and were updated (via communication
of the ghost cells) each time they were changed, using the library
function {\tt GCupdateGrid()}. Each coarse grid in the multigrid level
is partitioned in the same way as the original fine grid, and is
communicated through separate calls to the library. It should be noted
that although all the local grids at the same multigrid level always
work in parallel, the grids at different levels are still being executed
serially on each processor. The library tries to assign equal number of
nodes to each processor to achieve proper load balancing, and does all
the communication through ghost cells at the boundaries.

The implementation was extremely general with minimal assumptions
being made regarding the partition type or the grid size. As a result,
the program can run with grids of arbitrary size $m\times n$ on
arbitrary distribution of processors $p \times q$. The maximum
size of the problem being solved is limited only by the available system
memory. An $auto-levels$ option is provided which selects the optimal
number of levels for the multigrid solver based on the size of the
grid and the partitioning given. This option ensures that the size
of the coarsest grid is never smaller than a fixed size ($8\times8$
in our case) on any processor, thus avoiding unnecessary communication
(hence latency) for the extremely small grids.
The number of multigrid levels and various other parameters
(e.g., relaxation method to be used, value of $\omega$ for SOR, 
number of pre-smoothing and post-smoothing cycles, maximum error
allowed, maximum number of iterations to run, etc.)
can be entered by the user as command line arguments to the program (see
Appendix~\ref{ch:appendixSyntax}). This made it extremely simple and
hassle-free to run several test cases with varying parameters without
ever recompiling the program. The boundary condition and the initial
condition can be specified using separate user-defined functions {\tt
setBoundary(x,y)} and {\tt setGrid(x,y)}. The source code for the
driving program can be found in Appendix~\ref{ch:appendixmain} 
and~\ref{ch:appendixmultigrid}.

\section{Results and Comparison with Theory}
Several hundred parallel and serial runs were carried 
on COCOA (COst effective COmputing Array,
{\tt http://cocoa.ihpca.psu.edu}) in order to generate
all the plots in this section. COCOA is a 50-processor PII-400 Beowulf
cluster running Linux and connected via 100 Mbit/s fast ethernet switch.
For our test case, the boundary condition was set to $0$ everywhere
and the grid was initialized to $100$. Thus, we were solving the
Laplace equation (with RHS term zero in Eq.~\ref{eq:poissons}). The exact
solution given by the program for this problem (within error $10^{-12}$)
is illustrated in figure~\ref{fig:soln_case1}. All our convergence
and parallel efficiency plots correspond to this case. Some other
test cases were also carried out (figures~\ref{fig:soln_case2}
and~\ref{fig:soln_case3}) for the validation of the program. We used $3$
pre-smoothing and post-smoothing cycles in all our runs, and ran for
multigrid levels $2$ to $5$.

\subsection{Convergence Rates}
For the study of the convergence rates of each method, we did a comparison
of {\em Jacobi}, {\em Gauss-Siedel} and {\em SOR} (with $\omega=1.25$)
 smoothing with various
multigrid levels. All the convergence data  looks linear in the 
semi-log plot suggesting that the convergence is occurring at an exponential
rate, something that matches very well with the theory. According to
the theory, for our test case, the error for the Jacobi case should
decrease by a factor $\rho_J$ during every iteration, which is given by:
$$
\rho_J = \cos \pi h \approx 1 - \frac{(\pi h)^2}{2}
$$
where $h$ is the grid spacing. For Gauss-Siedel, this factor is:
$$
\rho_{GS} = \cos^2 \pi h \approx 1 - (\pi h)^2
$$
and for
SOR, it is:
$$
\rho_{SOR} = \frac{1-\sin \pi h}{1+\sin \pi h} \approx
1 - 2\pi h
$$

For our case with a $200\times 200$ grid, $h = 1/200 = 0.005$ and therefore
$\rho_J = 0.9998766$, $\rho_{GS} = 0.9997533$ and
$\rho_{SOR} = 0.9690712$.
For the single grid case (no multigrid) depicted in
figure~\ref{fig:1level}, we can see that Gauss-Siedel converges almost
exactly twice as fast (i.e., using half the number of iterations) than
the Jacobi relaxation for a fixed error norm. Looking at the above
formulae for $\rho_J$ and $\rho_{GS}$, this is expected as $\rho_{GS}
= \rho^2_J$. SOR took yet $40\%$ fewer iterations, proving to be
the fastest among the three schemes (again, this is expected since
$\rho_{SOR}$ is the smaller value above). However, according to the
theory, a much faster convergence was expected here. The reason for our
slower convergence may be due to the fact that the $\rho_{SOR}$ obtained
above is for the value $\omega_{opt} = 1.5$, whereas $\omega$ for our
SOR run was $1.25$ (the value of $1.25$ was chosen as a precautionary
measure, since some of our test cases diverged for $\omega=1.5$). The
same trend was seen for the other cases (2-level to 5-level) as depicted
in figures~\ref{fig:2level}, \ref{fig:3level}, \ref{fig:4level} and
\ref{fig:5level}. From the comparison of various multigrid levels for
each relaxation scheme as shown in figures~\ref{fig:jac_alllevel},
\ref{fig:gs_alllevel} and ~\ref{fig:sor_alllevel}, it is also clear
that the increasing of the multigrid levels has a great impact on
the convergence rate. Roughly $3$ times less number of iterations is
required for the 2-level multigrid case as opposed to the conventional
1-level (i.e., single grid) case, and the same is seen when going from
the 2-level to the 3-level case.

\subsection{Parallel Efficiency}
In order for any parallel implementation to be usable, it has to have
a reasonably high parallel efficiency. The parallel efficiency of
the multigrid implementation was tested for various 1D and 2D grid
partitions. A fixed 3-level multigrid scheme was used for these
tests. Both fixed problem size (Amdahl's Law)
and scaled problem size (Gustafson's Law) cases were run. 

For the fixed problem size case using 1D partitioning, two grid sizes
of $100\times100$ and $1000\times1000$ were tested for each of the
relaxation schemes. As seen in figure~\ref{fig:1Dfixed}, it is clear
that the parallel efficiency decreases rapidly with the number of
processors for the small problem size of $100\times100$ ($10\%$ for
16-proc case), whereas it only suffers a small decline in the large
problem size case of $1000\times1000$ ($75\%$ for 16-proc case). This is
well in line with theory, since the computation cost is fixed here, and
the communication cost rapidly increases as the number of processors
increase, thus increasing the communication-to-computation ratio
(which is $1.5\%$ for the $1000\times1000$ case, and $15\%$ for
the $100\times100$ case, both for 16 processors), and
hence the parallel efficiency. For the scaled problem size case using
1D partitioning, again two grid sizes of $100\times100 = 10^4$ nodes
and $1000\times1000 = 10^6$ nodes were used, this time for each local
grid (or processor). From figure~\ref{fig:1Dscaled}, it is noted that
the parallel efficiency again reduces greatly with the number of
processors for the $100\times100$ case, but not as rapidly as for the
fixed size case ($50\%$ for 16-proc case). For the $1000\times1000$
case, the parallel efficiency is only slightly reduced with the number
of processors ($80\%$ for 16-proc case). The reduction in parallel
efficiency is much less here, since the communication-to-computation
ratio is much lower and grows at a very small ratein this case 
($0.71\%$ to $3.75\%$ for
$100\times100$ case with 2-16 processors, 
and $0.071\%$ to $0.375\%$ for the $1000\times1000$
case with 2-16 processors).

Similar trend as above is seen for the 2D partitioning case, as
depicted in figures~\ref{fig:2Dfixed} and \ref{fig:2Dscaled}.
For the 2D case with fixed problem size with 16 processors, the
communication-to-computation ratio is $6\%$ for the $100\times100$,
and $0.6\%$ for the $1000\times1000$ case. For scaled problem size,
it is $1.5\%$ for the $100\times100$ case, and $0.15\%$ for the
$1000\times1000$ case. The $100\times100$ grid seemed to perform
slightly worse than the corresponding case for 1D partitioning
with fixed problem size, whereas the opposite is seen for the
$1000\times1000$ grid. For the scaled problem size case, the same
trend was observed. This can be explained easily, as a slightly worse
parallel efficiency is expected here when the problem size is small due
to the communication in two directions for every processor as opposed
to just one in the 1D partitioning case (thus, the network latency will
be double in this case although the transfer time will reduce a bit).
For large grids however, an advantage is to be expected, as the total
length of the shared data decreases in the 2D partitioning (consider,
for example, a $2\times2$ and $4\times 1$ partitioning for a $N\times N$
grid; here the total shared data is $2N$ for the 2D case and $3N$ for
the 1D case). Thus, the loss due to the high network latency is more
than made up with the savings in the transfer time.

The final plot, figure~\ref{fig:parkxk} studies the effect of grid size
($k \times k$) on the parallel efficiency for the fixed partitioning
cases of $9\times1$ processor and $3\times3$ processor. It is clear
from the figure that the parallel efficiency goes up with the increase
in the grid size until a certain critical size is reached ($1000\times
1000=10^6$ in our case, for which the parallel efficiency is $85\%$).
Beyond this size, the parallel efficiency is almost constant.

\section{Conclusions}
A successful parallel implementation of the multigrid technique to solve
the Poisson's equation using various relaxation algorithms has been
carried out. The implementation is both robust and very general, and
efficient. It is seen that SOR is the best relaxation algorithm for all
the cases studied above. It is also very clear from the exercise that
multigrid offers tremendous advantage over the conventional one-grid
method. The more the number of multigrid levels, the faster is the
convergence rate. Thus, SOR with 3-5 levels offers the best advantage
for most cases. A comprehensive study of the parallel efficiency
for various partitioning scenarios and grid sizes has been done on
COCOA. It is found that scaling the problem size with the number of
processors offers maximum parallel efficiency in all cases, and hence
maximum utilization of resources. Fixed problem size cases are best
run on several processors only when the problem size is reasonably
large (around $1000\times1000$ in our case), otherwise the communication
dominates and kills the efficiency. 1D partitioning is recommended for
smaller problem sizes, whereas 2D partitioning is more suitable for
larger problem sizes. While the performance numbers may vary from one
computer/cluster to another, the trend seen here is expected to remain
the same.

% Sample solutions
\myfigure{plots/final_sol_case1.eps}{soln_case1}{8.7cm}{9.9cm}
	{Contour plot of converged solution ($\| Residual \| = 10^{-12}$) for 
a $1000\times1000$ test case (BC: $f_{BC}(x,y)=100$)
using 5-level Multigrid with {\em SOR}, run on a $3\times3$ processor configuration}

\myfigure{plots/final_sol_case2.eps}{soln_case2}{8.7cm}{9.9cm}
	{Contour plot of converged solution ($\| Residual \| = 10^{-12}$) for 
a $1000\times1000$ test case (BC: $f_{BC}(x,y)=\sin 2 \pi x +\sin 2 \pi y$)
using 5-level Multigrid with {\em SOR}, run on a $3\times3$ processor configuration}

\myfigure{plots/final_sol_case3.eps}{soln_case3}{8.7cm}{9.9cm}
	{Contour plot of converged solution ($\| Residual \| = 10^{-12}$) for 
a $1000\times1000$ test case (BC: $f_{BC}(x,y)= x + y$)
using 5-level Multigrid with {\em SOR}, run on a $3\times3$ processor configuration}

% Plots of comparison of convergence rates of Jacobi, GS & SOR 
% for various Multigrid levels 
\myfigure{plots/n200_l1_p1.eps}{1level}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Jacobi}, {\em Gauss-Siedel}
and {\em SOR} without using Multigrid acceleration ($200\times200$ grid on 1 processor)}

\myfigure{plots/n200_l2_p1.eps}{2level}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Jacobi}, {\em Gauss-Siedel}
and {\em SOR} without 2-level Multigrid acceleration ($200\times200$ grid on 1 processor)}

\myfigure{plots/n200_l3_p1.eps}{3level}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Jacobi}, {\em Gauss-Siedel}
and {\em SOR} without 3-level Multigrid acceleration ($200\times200$ grid on 1 processor)}

\myfigure{plots/n200_l4_p1.eps}{4level}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Jacobi}, {\em Gauss-Siedel}
and {\em SOR} without 4-level Multigrid acceleration ($200\times200$ grid on 1 processor)}

\myfigure{plots/n200_l5_p1.eps}{5level}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Jacobi}, {\em Gauss-Siedel}
and {\em SOR} without 5-level Multigrid acceleration ($200\times200$ grid on 1 processor)}

% Plots of comparison of convergence rates of Jacobi, GS & SOR 
% for various Multigrid levels
\myfigure{plots/jacobi_ml_n200_p1.eps}{jac_alllevel}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Jacobi}
with varying levels in Multigrid acceleration ($200\times200$ grid on 1 processor)}

\myfigure{plots/gs_ml_n200_p1.eps}{gs_alllevel}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em Gauss-Siedel}
with varying levels in Multigrid acceleration ($200\times200$ grid on 1 processor)}

\myfigure{plots/sor_ml_n200_p1.eps}{sor_alllevel}{9cm}{10.08cm}
	{Comparison of convergence rates of {\em SOR}
with varying levels in Multigrid acceleration ($200\times200$ grid on 1 processor)}

% Parallel efficiency for 1D partitioning case
\myfigureB{plots/1Dfixed_all.eps}{plots/1Dscale_all.eps}{1Dfixed}{1Dscaled}
	{9cm}{10.08cm}
	{Parallel efficiency for 1D partitioning of processors ($N\times1$)
for 3-level Multigrid scheme with fixed global grid sizes of $100\times100$ and $1000\times1000$}
	{Parallel efficiency for 1D partitioning of processors ($N\times1$)
for 3-level Multigrid scheme with fixed local grid sizes (or scaled global
grid size) of $10^4$ nodes/proc ($100\times100$) and 
$10^6$ nodes/proc ($1000\times1000$)}

% Parallel efficiency for 2D partitioning case
\myfigureB{plots/2Dfixed_all.eps}{plots/2Dscale_all.eps}{2Dfixed}{2Dscaled}
	{9cm}{10.08cm}
	{Parallel efficiency for 2D partitioning of processors ($N\times N$)
for 3-level Multigrid scheme with fixed global grid sizes of $100\times100$ and $1000\times1000$}
	{Parallel efficiency for 2D partitioning of processors ($N\times N$)
for 3-level Multigrid scheme with fixed local grid sizes (or scaled global
grid size) of $10^4$ nodes/proc ($100\times100$) and 
$10^6$ nodes/proc ($1000\times1000$)}


% Parallel efficiency for 2D partitioning case
\myfigure{plots/pareffvsGrid_9all.eps}{parkxk}{9cm}{10.08cm}
	{Parallel efficiency vs. grid sizes $k\times k$
for 1D and 2D partitioning for the $9$ processor case for all the
three relaxation methods}

\appendix
\begin{scriptsize}
\include{prog_syntax}
\end{scriptsize}

\begin{scriptsize}
\include{prog_main}
\end{scriptsize}

\begin{scriptsize}
\include{prog_multigrid}
\end{scriptsize}

\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
