src/gpu/pennant/doc/pennantdoc.tex

%%%%%%%%%% BEGIN PREAMBLE %%%%%%%%%%
\documentclass[11pt,letterpaper]{article} % 11pt font = 12pt Word font
\usepackage{graphicx}
\usepackage{fancyhdr,lastpage}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage[english]{babel}
\usepackage{helvet}
\usepackage{hyperref}
\usepackage{pdfpages}
\usepackage{times}
\usepackage{comment}
\usepackage{geometry}

\geometry{pdftex, portrait,
          headsep=8pt, headheight=8pt, footskip=16pt,
          top=88pt, bottom=88pt,
          left=72pt, right=72pt
}

\begin{document}

\title{The PENNANT Mini-App}
\author{Charles R. Ferenbaugh \\
        Los Alamos National Laboratory \\
        {\tt cferenba@lanl.gov}}
\date{Version 0.5 -- May 2013}
\maketitle

PENNANT is an unstructured mesh physics mini-app designed for advanced
architecture research.
It contains mesh data structures and a few physics algorithms adapted
from the LANL rad-hydro code FLAG, and gives a
sample of the typical memory access patterns of FLAG.
%It is capable of running several
%standard test problems in shock physics.

\section{Building and running the code}

\subsection{Building}

A simple {\tt Makefile} is provided in the top-level directory for building
the code.  Before using it, you may wish to edit the definitions of CXX
and CXXFLAGS to specify your desired C++ compiler and flags, and to
choose between optimized/debug and serial/OpenMP builds.  Then
a simple ``{\tt make}'' command will create a {\tt build} subdirectory and
build the {\tt pennant} binary in that directory.

PENNANT has been tested under GCC 4.6.1, PGI 11.10, and Intel 12.1.2.
Building under other compilers should require only minor changes.

\subsection{Running tests}

Several test problems are provided in subdirectories under the {\tt test}
directory.  The command line
\begin{quote}
{\tt pennant} {\em testname}{\tt.pnt}
\end{quote}
is used to run a test.

The available test problems are listed in Table~\ref{tbl:tests}.
The smaller problems run quickly and are useful for debugging and
regression tests; gold standard files are provided (see next section).
The larger tests take longer to run and are suitable for timing tests.

\begin{table}
\centering
\caption{Test problems provided with PENNANT.}
\label{tbl:tests}
\begin{tabular}{lll}
    \hline
    name & \# zones & zone type \\
    \hline
    Leblanc problems~\cite{leblanc}: \\
    leblanc    &   900 & square \\
    leblancbig & 57600 & square \\
    \hline
    Noh problems~\cite{noh}: \\
    nohsmall   &    40 & triangle/quad \\
    nohsquare  & 32400 & square \\
    nohpoly    & 22801 & mostly hexagons \\
    \hline
    Sedov problems~\cite{sedov}: \\
    sedovsmall &    81 & square \\
    sedov      &  2025 & square \\
    sedovbig   & 72900 & square \\
    \hline
\end{tabular}
\end{table}

\subsection{Test inputs and outputs}

Each test problem directory contains two input files.  The {\tt.pnt}
file is a small text file containing input parameters for the test.
The {\tt.gmv} file contains the input mesh in LANL GMV format; its name
is specified as one of the entries in the {\tt.pnt} file.

PENNANT generates output files of two kinds.  The {\tt.xy} output file
is a text file containing the per-zone values of zone density, energy,
and pressure.  (It is modeled after a similar file generated by FLAG.)
For the smaller tests, a gold standard file {\em testname}{\tt.xy.std}
is provided for reference.

There are also several graphics output files in Ensight Gold format:
the main file has suffix {\tt.case}, and it refers to auxiliary files
with suffixes {\tt.geo}, {\tt.ze}, {\tt.zp}, and {\tt.zr}.  These
can be viewed by the proprietary
Ensight\footnote{\url{http://www.ensight.com}}
viewer, or by open-source viewers such as
ParaView\footnote{\url{http://www.paraview.org}} and
VisIt\footnote{\url{https://wci.llnl.gov/codes/visit/home.html}}.
Sample outputs are shown in Figure~\ref{fig:output}.
\begin{figure}
    \centering
    \begin{tabular}{ccc}
    \includegraphics[width=0.40\textwidth]{noh-result.png} &
    \hspace{0.06\textwidth} &
    \includegraphics[width=0.40\textwidth]{sedov-result.png} \\
    (a) && (b) \\
    \end{tabular}
    \caption{Final state of (a) {\tt nohsquare} and (b) {\tt sedovbig}
    problems, colored by zone density.}
    \label{fig:output}
\end{figure}

For users who would like to generate additional test cases of varying
sizes (e.g., for scaling studies), the utility script {\tt gmvrect.py}
is provided in the {\tt tools} directory.  This generates rectangular
meshes containing square zones, with user-specified dimensions.
Usage is:
\begin{verbatim}
Usage:  gmvrect.py NZX [NZY [LENX [LENY]]]
where nzx, nzy   = number of zones in x, y directions
                   (no default for nzx; default nzy = nzx)
      lenx, leny = total length in x, y directions
                   (default for both = 1.0)
\end{verbatim}

Note that changing the number
of zones in a problem typically requires a corresponding change to the
{\tt dtinit} parameter in the input file.  As a rule of thumb, if the
resolution of the problem is increased by a factor of $r$, {\tt dtinit}
must decrease by a factor of $r$.

\section{Data structure details}

\subsection{Mesh data structures}

PENNANT is designed to use standard finite-volume meshes similar to
those used by many common physics solvers.  In particular, PENNANT
supports 2-D unstructured meshes composed of arbitrary polygons.

The PENNANT mesh data structures are a subset of those used by FLAG.
These are implemented in the {\tt Mesh} class.
FLAG supports 1-, 2-, and 3-D meshes with various geometries; for
simplicity, PENNANT is restricted to the 2-D, cylindrical geometry
case.

The PENNANT terminology for entities within a mesh is shown in
Figure~\ref{fig:mesh}.  The basic mesh elements in 0, 1, and 2 dimensions
are called {\em points, edges,} and {\em zones} respectively.
PENNANT also uses two types of sub-zone entities.  Within any given
zone, a {\em side} is a triangle whose vertices are two consecutive
boundary points of the zone together with the zone center.  A {\em corner}
is a quadrilateral whose vertices are one boundary point of a zone,
the midpoints of the two adjoining edges, and the zone center.

\begin{figure}
    \centering
    \includegraphics[width=0.50\textwidth]{mesh-entities.png}
    \caption{PENNANT terminology for mesh entities.}
    \label{fig:mesh}
\end{figure}

For each entity type, the first letter of its name is used as an
identifier for variables associated with it.  For example, {\tt px}
is an array of point coordinates, and {\tt zvol} is an array
of zone volumes.  There are scalar variables of the form
{\tt numX} which give the number of {\tt X}'s in the problem:
{\tt nump} is the number of points, and so on.

PENNANT also stores various mapping arrays between different entity
types.  Most of these are side-based, as shown in Figure~\ref{fig:side}.
Given a side $s$, the following mapping arrays are available:
\begin{itemize}
\item {\tt mapsz} gives the zone $z$ of which $s$ is a subregion.
\item {\tt mapse} gives the edge $e$ on the boundary of $s$.
\item {\tt mapsp1} and {\tt mapsp2} give the two mesh points $p_1$ and
$p_2$ on the boundary of $s$.  It is assumed that the mesh is oriented
according to a right-hand rule, so that the edge from $p_1$ to $p_2$
is always in a counter-clockwise direction relative to the zone.
\item {\tt mapsc1} and {\tt mapsc2} give the two corners $c_1$ and
$c_2$ which overlap with $s$, where $c_1$ is before $c_2$ in
a counter-clockwise traversal of the zone.
\item {\tt mapss3} and {\tt mapss4} give the two sides $s_3$ and
$s_4$ on either side of $s$, where $s_3$ is before $s$ and $s_4$ is
after it in a counter-clockwise traversal of the zone.
\footnote{If you're wondering why there isn't a {\tt mapss2}:
FLAG uses $s_2$ to denote the side that is across the edge $e$ from
side $s$.  This variable isn't currently needed by PENNANT, so I didn't
implement it or include it in the figure.  It may be implemented in
a future PENNANT release.}
\end{itemize}
There are also two edge-based arrays {\tt mapep1} and {\tt mapep2},
which give the two endpoints of an edge~$e$; and two corner-based
arrays {\tt mapcz} and {\tt mapcp}, which give the associated zone
and point, respectively, of a corner~$c$.

\begin{figure}
    \centering
    \includegraphics[width=0.40\textwidth]{side-maps.png}
    \caption{Various side map arrays supported by PENNANT.}
    \label{fig:side}
\end{figure}

The {\tt Mesh} class also has methods for computing other
geometry-related variables, such as edge and zone centers, lengths,
volumes, and surface vectors.  The surface vector for a side $s$, shown
by the vector {\bf S} in Figure~\ref{fig:side}, is used by force
computations in the hydro algorithms described later.

\subsection{Chunk processing}

The PENNANT {\tt Mesh} class has been set up to support computation
on chunks of the mesh in parallel.  This is used to implement an OpenMP
version of the code and should lend itself to other task-based approaches
as well.

The maximum chunk size is controlled by the input file parameter
{\tt chunksize}.  If {\tt chunksize} is larger than the total number of
sides in the mesh, the entire mesh is treated as a single chunk
(this is the default).

In the current chunking approach, the lists of points and zones are
simply divided into chunks of size {\tt chunksize} (except for the final,
leftover chunk).  The list of sides is handled similarly, except that
the size of each individual chunk is rounded down slightly if necessary
so that each zone has all of its sides in the same chunk.  For each chunk,
the {\em first} and {\em last} indices of the chunk are stored.
(Note that {\em last} is actually one index beyond the end of the chunk,
in a similar manner to STL iterators, so that the sides in a side chunk
are those with $s_{first} \le s < s_{last}$.)  The total numbers of each
type of chunk ({\tt numXch}, where {\tt X} may be {\tt p}, {\tt s}, or
{\tt z}) are also stored.

Then, nearly all of the routines in the main hydro cycle have been
modified to take as input first and last indices of the appropriate
mesh entity.  This allows the hydro processing to be divided into
five phases, two on point chunks, two on side chunks, and one on zone
chunks.  Within each phase, all chunks are independent and can be
processed in parallel.
See {\tt Hydro::doCycle()} for the complete code flow.

A few of the helper routines, particularly in the {\tt QCS} class,
use scratch arrays the size of the chunk currently being processed.
The prefix {\tt s0} is used for an array with one entry per side
in the current chunk, with prefixes {\tt c0} and {\tt z0} used similarly
for corners and zones respectively.

\section{Physics details}

\subsection{Basic hydro algorithms}

PENNANT provides a subset of the compatible Lagrangian staggered grid
hydrodynamics (SGH) algorithms implemented in FLAG and described
in~\cite{hydro}.  These are implemented in the {\tt Hydro} class.
An outline of the main steps is given in Table~\ref{tbl:dataflow}.

\def\RR{\raggedright}

\begin{table}
\caption{Basic data flow for FLAG/PENNANT hydrodynamics.}
\label{tbl:dataflow}
\begin{tabular}{lp{108pt}p{180pt}}
    \hline
    step & main inputs & main outputs \\
    \hline
    Predictor step: \\
    1. Update mesh & point velocity, position & point position (half-advanced) \\
    1a. Update mesh geometry & point position & {\bf {\RR side and zone volume; zone density; \\ side surface vector } } \\
    2. Compute point masses & {\RR zone density, volume; \\ side mass fraction} & {\bf point mass} \\
    3. Update thermodynamic state & {\RR zone density, specific\\ energy, work rate} & zone pressure, sound speed \\
    4. Compute forces & {\RR side surface vector; \\ zone pressure} & side and {\bf point force} \\
    4a. Apply boundary conditions & point force, velocity & point force, velocity (constrained) \\
    5. Compute accleration & point force & point acceleration \\
    \hline
    Corrector step: \\
    6. Update mesh & {\RR point acceleration, \\ velocity, position} & point velocity, position (fully advanced) \\
    6a. Update mesh geometry & point position & {\bf side and zone volume} \\
    7. Compute work & {\RR point position, velocity; \\ side force} & {\bf zone work, work rate, total energy} \\
    8. Update zone state variables & {\RR zone volume, mass,\\ total energy} & zone density, specific energy \\
    \hline
\end{tabular}
\end{table}

The PENNANT hydro algorithm is a {\em Lagrangian} method, meaning that
the computational mesh moves with the material as the problem state
advances.  This implies that the mass and material type within each zone
are constant throughout the problem, but the zone's position and shape
will change over time.

It is also a {\em staggered-grid} method, meaning that mesh positions and
related variables (velocity, acceleration, etc.) are stored on points,
while most state variables (density, energy, pressure, etc.) are stored
on zones.  Therefore, the calculation must frequently use values of
zone-based variables to compute point-based results, or vice-versa.
In Table~\ref{tbl:dataflow}, such results are shown in {\bf bold}.
(Note that this is true for five of the 11 steps shown.)

To facilitate these calculations, many of the calculation loops are done
over sides, and some intermediate variables are stored on sides.
This works since each side can be easily correlated to its
corresponding zone and points using the mapping arrays in the {\tt Mesh}
object.

PENNANT hydro uses a {\em predictor-corrector} time integration
method.  Each cycle can be broken into two steps, shown in the table.
The cycle begins with all problem state defined for the beginning of the
timestep.  In the predictor step, some variables are advanced to the
middle of the timestep, in order to compute half-advanced point
acceleration values.  In the corrector step, the new accelerations are
then used to advance all variables to the end of the timestep.

To implement the predictor-corrector scheme, it is necessary to store
multiple values of some of the problem variables.  This is done
using the following notation convention:
\begin{itemize}
\item suffix $0$ = the beginning of the timestep (``cycle n'')
\item suffix $p$ = half-way through the timestep (``cycle n + 1/2'')
\item no suffix = completion of the timestep    (``cycle n + 1'')
\end{itemize}
Some examples are shown in Table~\ref{tbl:timestep}.
Note that some entries in the table are blank, since not all quantities
are needed at all times.

\begin{table}
\centering
\caption{Examples of PENNANT variables and their dependence on timestep.}
\label{tbl:timestep}
\begin{tabular}{lccc}
    \hline
    & \multicolumn{3}{c}{Part of time step} \\
    quantity & begin & middle & end \\
    \hline
    point coordinate       &   px0     &   pxp     &   px  \\
    point velocity         &   pu0     &           &   pu    \\
    point acceleration     &           &   pap  \\
    point force            &           &   pf  \\
    point mass             &           &   pmaswt  \\
    \\
    zone center coordinate &           &   zxp      &   zx  \\
    zone mass              &           &           &   zm      \\
    zone volume            &   zvol0   &   zvolp   &   zvol       \\
    zone density           &           &   zrp      &   zr      \\
    zone specific energy   &           &           &   ze  \\
    \\
    side volume            &           &   svolp    &   svol    \\
    side surface vector    &           &   ssurfp \\
    \hline
\end{tabular}
\end{table}


\subsection{Material model}

PENNANT provides finite-volume, arbitrary-polygon cells with a
gas material model, implemented in the {\tt PolyGas} class.
This class includes code to compute a simple gamma-law gas
equation of state, and to compute the resulting pressure-based
forces.


\subsection{Subzonal pressures}

PENNANT provides the Temporary Triangular Subzoning (TTS) algorithm
described in~\cite{szp,tts}.  This is implemented in the {\tt TTS}
class.  This prevents certain kinds of distortions of zones, such as
``hourglassing,'' by estimating a pressure for each side, and adding
a force to each side based on the difference
between the zone and side pressures.

Note to FLAG users:  The FLAG implementation of TTS contains, in addition
to the subzonal pressure treatment, an artificial viscosity algorithm
based on the subzonal pressures; the artificial viscosity is not part of
the standard TTS description in the references.  Only the subzonal pressure
part of TTS is implemented in PENNANT.


\subsection{Artificial viscosity}

PENNANT provides the tensor artificial viscosity algorithm of Campbell and
Shashkov, described in~\cite{qcs}.  This is implemented in the {\tt QCS}
class.  (The symbol $q$ is traditionally used to denote artificial
viscosity, hence the {\tt Q} prefix on the class name.)
Artificial viscosity is a fictitious term commonly introduced into
fluid flow equations to correctly handle
shock regions with large discontinuities in
the problem state variables.


\section*{Acknowledgements}

Thanks to Mikhail Shashkov and the ASCR ``Mimetic Methods for PDEs''
project, and the ASC Hydrodynamics project, for providing support for
this work.
Thanks also to the Lagrangian Applications Project members who have
contributed to the FLAG code; parts of the PENNANT code and documentation
are adapted from their work.
And finally, thanks to the participants in the Intel EPOCH workshop
in the summer of 2012; several optimization ideas from that workshop
have been incorporated into subsequent PENNANT versions.


\begin{thebibliography}{9}
\bibliographystyle{alpha}

\bibitem{leblanc}
D.J.~Benson.
\newblock Momentum advection on a staggered mesh.
\newblock {\em J. Comput. Phys}, 100:143--162, 1992.

\bibitem{qcs}
J.~Campbell and M.~Shashkov.
\newblock A tensor artificial viscosity using a mimetic finite difference
  algorithm.
\newblock {\em J. Comput. Phys}, 172:739--765, 2001.

\bibitem{hydro}
E.J. Caramana, D.E. Burton, M.~Shashkov, and P.P. Whalen.
\newblock The construction of compatible hydrodynamics algorithms utilizing
  conservation of total energy.
\newblock {\em J. Comput. Phys.}, 146:227--262, 1998.

\bibitem{szp}
E.J. Caramana and M.J. Shashkov.
\newblock Elimination of artificial grid distortion and hourglass-type motions
  by means of {L}agarangian subzonal masses and pressures.
\newblock {\em J. Comput. Phys.}, 142:521--561, 1998.

\bibitem{noh}
W. F. Noh,
\newblock Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux.
\newblock {\em J. Comput. Phys.}, 72:78, 1987.

\bibitem{sedov}
L. I. Sedov,
\newblock Similarity and Dimensional Methods in Mechanics.
\newblock Academic Press, New York, 1959.

\bibitem{tts}
K.B. Wallick.
\newblock Temporary triangular subzoning ({TTS}), in {REZONE}: A method for
  automatic rezoning in two-dimensional lagrangian hydrodynamics problems.
\newblock {Technical {R}eport LA-10829-MS}, Los Alamos National Laboratory, Los
  Alamos, NM 1987.

\end{thebibliography}

\appendix
\section{Version Log}

\begin{description}
\item[0.5] May 2013.  Further optimizations.
\item[0.4] January 2013.  First open-source release.
Fixed a bug in QCS and added some optimizations.
Added Sedov and Leblanc test problems, and some new input keywords to
support them.
\item[0.3] July 2012.  Added OpenMP pragmas and point chunk processing.
Modified physics state arrays to be flat arrays instead of STL vectors.
\item[0.2] June 2012.  Added side chunk processing.  Miscellaneous minor cleanup.
\item[0.1] March 2012.  Initial release, internal LANL only.
\end{description}

\section{Copyright and License Information}

Copyright \copyright 2012, Los Alamos National Security, LLC.
All rights reserved.

Copyright 2012. Los Alamos National Security, LLC.
This software was produced under U.S. Government contract
DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is
operated by Los Alamos National Security, LLC for the U.S. Department
of Energy. The U.S. Government has rights to use, reproduce, and
distribute this software. NEITHER THE GOVERNMENT NOR LOS ALAMOS
NATIONAL SECURITY, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is
modified to produce derivative works, such modified software should be
clearly marked, so as not to confuse it with the version available from
LANL.

Additionally, redistribution and use in source and binary forms, with
or without modification, are permitted provided that the following
conditions are met:

\begin{enumerate}
\item Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
\item Redistributions in binary form must reproduce the above
   copyright notice, this list of conditions and the following
   disclaimer in the documentation and/or other materials provided
   with the distribution.
\item Neither the name of Los Alamos National Security, LLC, Los Alamos
   National Laboratory, LANL, the U.S. Government, nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
\end{enumerate}

THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND
CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS
NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

\end{document}