source: trunk/anuga_core/source/anuga_parallel/documentation/parallel.tex @ 8950

\chapter{Running the Code in Parallel}

This chapter looks at how to run the code in parallel. The first section gives an overview of the main steps required to divide the mesh over the processors. The second sections describes some example code and the final section talks about how to run the code on a few specific architectures.


\section {Partitioning the Mesh}\label{sec:part}
There are four main steps required to run the code in parallel. They are;
\begin{enumerate}
\item partition the mesh into a set of non-overlapping submeshes
(\code{pmesh_divide_metis} from {\tt pmesh_divide.py}),
\item build a \lq ghost\rq\ or communication layer of triangles
around each submesh and define the communication pattern (\code{build_submesh} from {\tt build_submesh.py}),
\item distribute the submeshes over the processors (\code{send_submesh}
and \code{rec_submesh} from {\tt build_commun.py}),
\item and update the numbering scheme for each submesh assigned to a
processor (\code{build_local_mesh} from {\tt build_local.py}).
\end{enumerate}
See Figure \ref{fig:subpart}

\begin{figure}[hbtp]
  \centerline{ \includegraphics[scale = 0.75]{figures/domain.eps}}
  \caption{The main steps used to divide the mesh over the processors.}
  \label{fig:subpart}
\end{figure}

\subsection {Subdividing the Global Mesh}\label{sec:part1}

The first step in parallelising the code is to subdivide the mesh
into, roughly, equally sized partitions. On a rectangular mesh this may be
done by a simple co-ordinate based dissection method, but on a complicated
domain such as the Merimbula mesh shown in Figure \ref{fig:mergrid}
a more sophisticated approach must be used.  We use pymetis, a
python wrapper around the Metis
(\url{http://glaros.dtc.umn.edu/gkhome/metis/metis/overview})
partitioning library. The \code{pmesh_divide_metis} function defined
in {\tt pmesh_divide.py} uses Metis to divide the mesh for
parallel computation. Metis was chosen as the partitioner based on
the results in the paper \cite{gk:metis}.

\begin{figure}[hbtp]
  \centerline{ \includegraphics[scale = 0.75]{figures/mermesh.eps}}
  \caption{The Merimbula mesh.}
 \label{fig:mergrid}
\end{figure}

Figure \ref{fig:mergrid4} shows the Merimbula grid partitioned over
four processors. Note that one submesh may comprise several
unconnected mesh partitions. Table \ref{tbl:mer4} gives the node
distribution over the four processors while Table \ref{tbl:mer8}
shows the distribution over eight processors. These results imply
that Pymetis gives a reasonably well balanced partition of the mesh.

\begin{figure}[hbtp]
  \centerline{ \includegraphics[scale = 0.75]{figures/mermesh4c.eps}
  \includegraphics[scale = 0.75]{figures/mermesh4a.eps}}


  \centerline{ \includegraphics[scale = 0.75]{figures/mermesh4d.eps}
  \includegraphics[scale = 0.75]{figures/mermesh4b.eps}}
  \caption{The Merimbula grid partitioned over 4 processors using Metis.}
 \label{fig:mergrid4}
\end{figure}

\begin{table}
\caption{Running an 4-way test of
{\tt run_parallel_sw_merimbula_metis.py}}\label{tbl:mer4}
\begin{center}
\begin{tabular}{c|c c c c c c c c}
CPU & 0 & 1 & 2 & 3\\
Elements & 2757 & 2713 & 2761 & 2554\\
\% & 25.56\% & 25.16\% & 25.60\% & 23.68\%\
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{Running an 8-way test of
{\tt run_parallel_sw_merimbula_metis.py}} \label{tbl:mer8}
\begin{center}
\begin{tabular}{c|c c c c c c c c}
CPU & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\
Elements & 1229 & 1293 & 1352 & 1341 & 1349 & 1401 & 1413 & 1407\\
\% & 11.40\% & 11.99\% & 12.54\% & 12.43\% & 12.51\% & 12.99\% & 13.10\% & 13.05\%\\
\end{tabular}
\end{center}
\end{table}

The number of submeshes found by Pymetis is equal to the number of
processors; Submesh $p$ will be assigned to Processor $p$.

\subsection {Building the Ghost Layer}\label{sec:part2}
The function {\tt build_submesh.py} is the work-horse and is
responsible for setting up the communication pattern as well as
assigning the local numbering scheme for the submeshes.

Consider the example subpartitioning given in Figure
\ref{fig:subdomain}. During the \code{evolve} calculations Triangle
2 in Submesh 0 will need to access its neighbour Triangle 3 stored
in Submesh 1. The standard approach to this problem is to add an
extra layer of triangles, which we call ghost triangles. The ghost
triangles are read-only, they should not be updated during the
calculations, they are only there to hold any extra information that
a processor may need to complete its calculations. The ghost
triangle values are updated through communication calls. Figure
\ref{fig:subdomaing} shows the submeshes with the extra layer of
ghost triangles.

\begin{figure}[hbtp]
  \centerline{ \includegraphics[scale = 0.6]{figures/subdomain.eps}}
  \caption{An example subpartitioning of a mesh.}
 \label{fig:subdomain}
\end{figure}


\begin{figure}[hbtp]
  \centerline{ \includegraphics[scale = 0.6]{figures/subdomainghost.eps}}
  \caption{An example subpartitioning with ghost triangles. The numbers
in brackets shows the local numbering scheme that is calculated and
stored with the mesh, but not implemented until the local mesh is
built. See Section \ref{sec:part4}. }
 \label{fig:subdomaing}
\end{figure}

When partitioning the mesh we introduce new, dummy, boundary edges.
For example, Triangle 2 in Submesh 1 from Figure
\ref{fig:subdomaing} originally shared an edge with Triangle 1, but
after partitioning that edge becomes a boundary edge. These new
boundary edges are are tagged as \code{ghost} and should, in
general, be assigned a type of \code{None}. The following piece of
code taken from {\tt run_parallel_advection.py} shows an example.
{\small \begin{verbatim} T = Transmissive_boundary(domain)
domain.default_order = 2 domain.set_boundary( {'left': T, 'right':
T, 'bottom': T, 'top': T, 'ghost': Non e} )
\end{verbatim}}

Looking at Figure \ref{fig:subdomaing} we see that after each
\code{evolve} step Processor 0  will have to send the updated values
for Triangle 2 and Triangle 4 to Processor 1, and similarly
Processor 1 will have to send the updated values for Triangle 3 and
Triangle 5 (recall that Submesh $p$ will be assigned to Processor
$p$). The \code{build_submesh} function builds a dictionary that
defines the communication pattern.

Finally, the ANUGA code assumes that the triangles (and nodes etc.) are numbered consecutively starting from 0. Consequently, if Submesh 1 in Figure \ref{fig:subdomaing} was passed into the \code{evolve} calculations it would crash. The \code{build_submesh} function determines a local numbering scheme for each submesh, but it does not actually update the numbering, that is left to \code{build_local}.


\subsection {Sending the Submeshes}\label{sec:part3}

All of the functions described so far must be run in serial on
Processor 0. The next step is to start the parallel computation by
spreading the submeshes over the processors. The communication is
carried out by \code{send_submesh} and \code{rec_submesh} defined in
{\tt build_commun.py}. The \code{send_submesh} function should be
called on Processor 0 and sends the Submesh $p$ to Processor $p$,
while \code{rec_submesh} should be called by Processor $p$ to
receive Submesh $p$ from Processor 0.

As an aside, the order of communication is very important. If someone was to modify the \code{send_submesh} routine the corresponding change must be made to the \code{rec_submesh} routine.

While it is possible to get Processor 0 to communicate it's submesh to itself, it is an expensive and unnecessary communication call. The {\tt build_commun.py} file also includes a function called \code{extract_hostmesh} that should be called on Processor 0 to extract Submesh 0.


\subsection {Building the Local Mesh}\label{sec:part4}
After using \code{send_submesh} and \code{rec_submesh}, Processor $p$ should have its own local copy of Submesh $p$, however as stated previously the triangle numbering will be incorrect on all processors except Processor $0$. The \code{build_local_mesh} function from {\tt build_local.py} primarily focuses on renumbering the information stored with the submesh; including the nodes, vertices and quantities. Figure \ref{fig:subdomainf} shows what the mesh in each processor may look like.


\begin{figure}[hbtp]
  \centerline{ \includegraphics[scale = 0.6]{figures/subdomainfinal.eps}}
  \caption{An example subpartitioning after the submeshes have been renumbered.}
 \label{fig:subdomainf}
\end{figure}

Note that ghost triangles are not stored in the domain any differently to the other triangles, the only way to differentiate them is to look at the communication pattern. This means that the {\tt inundation} code is essentially the same whether it is being run in serial or parallel, the only difference is a communication call at the end of each \code{evolve} step and a check to ensure that the ghost triangles are not used in the time stepping calculations.

\section{Some Example Code}

Chapter \ref{chap:code} gives full listings of some example codes.

The first example in Section \ref{subsec:codeRPA} solves the advection equation on a
rectangular mesh. A rectangular mesh is highly structured so a coordinate based decomposition can be used and the partitioning is simply done by calling the
routine \code{parallel_rectangle} as shown below.
\begin{verbatim}
#######################
# Partition the mesh
#######################

# Build a unit mesh, subdivide it over numproces processors with each
# submesh containing M*N nodes

points, vertices, boundary, full_send_dict, ghost_recv_dict =  \
    parallel_rectangle(N, M, len1_g=1.0)
\end{verbatim}

Most simulations will not be done on a rectangular mesh, and the approach to subpartitioning the mesh is different to the one described above, however this example may be of interest to those who want to measure the parallel efficiency of the code on their machine. A rectangular mesh should give a good load balance and is therefore an important first test problem.


A more \lq real life\rq\ mesh is the Merimbula mesh used in the code shown in Section \ref{subsec:codeRPMM}. This example also solves the advection equation. In this case the techniques described in Section \ref{sec:part} must be used to partition the mesh. Figure \ref{fig:code} shows the part of the code that is responsible for spreading the mesh over the processors. We now look at the code in detail.

\begin{figure}[htbp]
\begin{verbatim}
if myid == 0:

    # Read in the test files

    filename = 'merimbula_10785.tsh'

    mesh_full = pmesh_to_domain_instance(filename, Advection_Domain)
    mesh_full.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))

    # Subdivide the mesh

    nodes, triangles, boundary, triangles_per_proc, quantities  =\
            pmesh_divide_metis(mesh_full, numprocs)

    # Build the mesh that should be assigned to each processor.
    # This includes ghost nodes and the communication pattern

    submesh = build_submesh(nodes, triangles, boundary, quantities, \
                            triangles_per_proc)

    # Send the mesh partition to the appropriate processor

    for p in range(1, numprocs):
      send_submesh(submesh, triangles_per_proc, p)

    # Build the local mesh for processor 0

     points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict =\
              extract_hostmesh(submesh, triangles_per_proc)

else:
    # Read in the mesh partition that belongs to this
    # processor (note that the information is in the
    # correct form for the ANUGA data structure

    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
             rec_submesh(0)
\end{verbatim}
  \caption{A section of code taken from {\tt run_parallel_merimbula_metis.py} (Section \protect \ref{subsec:codeRPMM}) showing how to subdivide the mesh.}
 \label{fig:code}
\end{figure}
\newpage
\begin{itemize}
\item
These first few lines of code read in and define the (global) mesh. The \code{Set_Stage} function sets the initial conditions. See the code in \ref{subsec:codeRPMM} for the definition of \code{Set_Stage}.
\begin{verbatim}
    filename = 'merimbula_10785.tsh'
    mesh_full = pmesh_to_domain_instance(filename, Advection_Domain)
    mesh_full.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))
\end{verbatim}

\item \code{pmesh_divide_metis} divides the mesh into a set of non-overlapping subdomains as described in Section \ref{sec:part1}.
\begin{verbatim}
    nodes, triangles, boundary, triangles_per_proc, quantities  =\
            pmesh_divide_metis(mesh_full, numprocs)
\end{verbatim}

\item The next step is to build a boundary layer of ghost triangles and define the communication pattern. This step is implemented by \code{build_submesh} as discussed in Section \ref{sec:part2}. The \code{submesh} variable contains a copy of the submesh for each processor.
\begin{verbatim}
    submesh = build_submesh(nodes, triangles, boundary, quantities, \
                            triangles_per_proc)
\end{verbatim}

\item The actual parallel communication starts when the submesh partitions are sent to the processors by calling \code{send_submesh}.
\begin{verbatim}
    for p in range(1, numprocs):
      send_submesh(submesh, triangles_per_proc, p)
\end{verbatim}


The processors receives a given subpartition by calling \code{rec_submesh}. The \code{rec_submesh} routine also calls \code{build_local_mesh}. The \code{build_local_mesh} routine described in Section \ref{sec:part4} ensures that the information is stored in a way that is compatible with the Domain datastructure. This means, for example, that the triangles and nodes must be numbered consecutively starting from 0.

\begin{verbatim}
    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict=\
             rec_submesh(0)
\end{verbatim}

Note that the submesh is not received by, or sent to, Processor 0. Rather     \code{hostmesh = extract_hostmesh(submesh)} simply extracts the mesh that has been assigned to Processor 0. Recall \code{submesh} contains the list of submeshes to be assigned to each processor. This is described further in Section \ref{sec:part3}.
%The \code{build_local_mesh} renumbers the nodes
\begin{verbatim}
    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict)=\
              extract_hostmesh(submesh, triangles_per_proc)
\end{verbatim}

\end{itemize}

\section{Running the Code}
\subsection{Compiling Pymetis and Metis}
Unlike the rest of ANUGA, Metis and its Python wrapper Pymetis are not built
by the \verb|compile_all.py| script. A makefile is provided to automate the
build process. Change directory to the \verb|ga/inundation/pymetis/| directory
and ensure that the subdirectory \verb|metis-4.0| exists and contains an
unmodified Metis 4.0 source tree. Under most varieties of Linux, build the
module by running \verb|make|. Under x86\_64 versions of Linux, build the
module by running \verb|make COPTIONS="-fPIC"|. Under Windows, build the
module by running \verb|make for_win32|. After the build completes, verify
that the module works by running the supplied PyUnit test case with
\verb|python test_metis.py|.
\subsection{Running the Job}
Communication between nodes running in parallel is performed by pypar, which
requires the following:
\begin{itemize}
\item Python 2.0 or later
\item Numeric Python (including RandomArray) matching the Python installation
\item Native MPI C library
\item Native C compiler
\end{itemize}
Jobs are started by running appropriate commands for the local MPI
installation. Due to variations in MPI environments, specific details
regarding MPI commands are beyond the scope of this document. It is likely
that parallel jobs will need to be scheduled through some kind of queuing
system. Sample job scripts are available for adaptation in section
\ref{sec:codeSJ}. They should be easily adaptable to any queuing system
derived from PBS, such as TORQUE.