\documentclass{article}
\usepackage{amsmath}
\usepackage[mathletters]{ucs}
\usepackage[utf8x]{inputenc}
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
\usepackage{graphics}
\usepackage[ps,dvips,xdvi,arrow,frame,matrix,graph,curve,cmtip]{xy}
\usepackage[pdftex]{graphicx}

\usepackage[breaklinks=true]{hyperref}
\usepackage{enumerate}

\setcounter{secnumdepth}{0}
\title{HPP, FHP, and FHCP Implementation}
\author{Abdulmajed Dakkak}
\date{June 28th, 2008}
\begin{document}
\maketitle

\section{Computational Fluid Dynamics with CA}

\subsection{Visualization}

In virtually all the papers read, the lattice grid is viewed as a
vector field. This is done for several reasons:

\begin{itemize}
\item 
  It's easy to visualize on paper whether the algorithm is working.
\item 
  It's the simplest to implement.
\end{itemize}
In NKS, for example, Wolfram uses a $96 \times 96$ subgrid to
approximate a vector direction on a 6 million lattice grid.
Previously in his 1987 paper, Wolfram used \ldots{}. for a 4600(?)
lattice grid.

\section{HPP Model}

We represent each lattice with $4$ vectors, each containing is
stored in its own matrix. This means that \verb!a[i][j]! contains
the \verb!a! vector component of the lattice on the $i^{th}$ row
and $j^{th}$ column. With this representation in place, we perform
the following to compute the next frame.

\begin{itemize}
\item 
  Resolve collisions --- the collisions are performed and stored in a
  temporary collision matrix.
\item 
  Propagation --- the propagation is performed and values are copied
  from the temporary collision matrix into the \verb!a!, \verb!b!,
  \verb!c!, and \verb!d! matrices.
\end{itemize}
\includegraphics[scale=0.1]{graphics/hpp-rep.png}

\subsection{Collision}

Four temporary matrices are created and are called \verb!K_a!,
\verb!K_b!, \verb!K_c!, and \verb!K_d! to signify which direction
they correspond to. The vector \verb!K_a[i][j]! is the updated
based on the HPP collision rules.

One can represent the collision rules as a series of binary
statements. The best example is found in%
\footnote{Wolf-Gladrow, Dieter. Lattice-Gas Cellular Automata and Lattice
Boltzmann Models: an Introduction. Santa Clara: Springer-Verlag
TELOS, 2000. (\href{../papers/files/Wol2000c.pdf}{pdf})}
where the authors represent the collision rules as a binary
statement. Examples of the same binary statement can be found in
other papers, however (\textbf{need to record them})

\$$a = a \otimes [(a \wedge c \wedge \neg (b \vee d)) \vee (b \wedge d \wedge \neg (a \vee c))]$\$

Where the logic operator symbols are

\begin{itemize}
\item 
  $\neg$ is the \emph{not} operator
\item 
  $\otimes$ is the \emph{xor} operator
\item 
  $\vee$ is the \emph{or} operator
\item 
  $\wedge$ is the \emph{and} operator
\end{itemize}
In \verb!C!, the following code implements the above rule:

\begin{verbatim}
for (x = 0; x < LATTICE_WIDTH; x++) {
    for (y = 0; y < LATTICE_HEIGHT; y++) {
        change = (a[x][y] & c[x][y] & ~(b[x][y] | d[x][y])) |
                 (b[x][y] & d[x][y] & ~(a[x][y] | c[x][y]));

        K_a[x][y] = a[x][y] ^ change;
        K_b[x][y] = b[x][y] ^ change;
        K_c[x][y] = c[x][y] ^ change;
        K_d[x][y] = d[x][y] ^ change;
    }
}
\end{verbatim}
The boolean expression allows us to perform the computation needed
without the use of \verb!if-else! statements. It should be also
noted that we will never go outside the bounds of the array, the
same cannot be said about the propagation step. Also note that the
expression

\begin{verbatim}
(a1[x][y] & c1[x][y] & ~(b1[x][y] | d1[x][y])) | 
(b1[x][y] & d1[x][y] & ~(a1[x][y] | c1[x][y]))
\end{verbatim}
Is repeated four times in the binary representation, which is why
we store it in a separate value to reduce three computations.

\subsection{Propagation}

The HPP propagation step is the process of moving velocity vectors
in the same direction they are moving in after the collisions are
resolved. Because of the way we oriented our grid, we made this
step unnecessarily complicated. Also it's not clear if our
representation has any advantages over a regular square grid.
Future implementation (if one is needed) would use a regular square
grid.

\includegraphics[scale=0.2]{graphics/hpp-grid.png}

We use the following code to propagation the vectors in the
lattice

\begin{verbatim}
a[i][j] = K_a[i-1][j-1];
b[i][j] = K_b[i-1][j+1];
c[i][j] = K_c[i+1][j+1];
d[i][j] = K_d[i+1][j-1];
\end{verbatim}
it should be evident that we will run out of bounds if we run the
above piece of code over the width of the lattice array. Since we
do not want to enforce any boundary conditions at this stage of the
process, we place our lattice on a toroidal surface. This can be
easily implement in \verb!C! by noting that the following way to
access and store array elements creates the toroidal surface:

\begin{verbatim}
t[(i+n + LATTICE_HEIGHT) % LATTICE_HEIGHT][(j+m + LATTICE_WIDTH) % LATTICE_WIDTH]
\end{verbatim}
and will also always be in the bounds of the matrix
representation.

\subsection{Optimizations}

Since C has no differentiation between a boolean value and a short,
we can overload a short variable to contain multiple vector values.
The boolean expression does not change, but as you see from the
figure below, the propagation step has to.

\includegraphics[scale=0.2]{graphics/propgrid.png}

\subsection{The Code}

There are a few differences between our implementation and the one
described above:

\begin{enumerate}[1.]
\item 
  The variables are called \verb!c1!, \verb!c2!, \verb!c3!, and
  \verb!c4! rather than \verb!a!, \verb!b!, \verb!c!, and \verb!d!.
  The temporary collision arrays have also been renamed to \verb!k1!,
  \verb!k2!, \verb!k3!, and \verb!k4!.
\item 
  Our computation of blocks of vectors at a time forces us to use
  checkerboard propagation.
\end{enumerate}
The crux of the code for HPP is provided bellow

\begin{verbatim}
int x, y;
short change;

// Resolve Collisions
for (x = 0; x < LATTICE_WIDTH; x++) {
    for (y = 0; y < LATTICE_HEIGHT; y++) {
        change = (a1[x][y] & c1[x][y] & ~(b1[x][y] | d1[x][y])) |
                 (b1[x][y] & d1[x][y] & ~(a1[x][y] | c1[x][y]));

        a2[x][y] = a1[x][y] ^ change;
        b2[x][y] = b1[x][y] ^ change;
        c2[x][y] = c1[x][y] ^ change;
        d2[x][y] = d1[x][y] ^ change;
    }
}
// Propagate
for (x = 1; x < LATTICE_WIDTH - 1; x++) {
    for (y = 1; y < LATTICE_HEIGHT - 1; y += 2) {
        a1[x][y] = (a2[x][y - 1] >> 1) + (a2[x - 1][y - 1] << LAST);
        b1[x][y] = b2[x][y - 1];
        c1[x][y] = c2[x][y + 1];
        d1[x][y] = (d2[x][y + 1] >> 1) + (d2[x - 1][y + 1] << LAST);

        a1[x][y + 1] = a2[x][y];
        b1[x][y + 1] = (b2[x][y] << 1) + (b2[x + 1][y] >> LAST);
        c1[x][y + 1] = (c2[x][y + 2] << 1) + (c2[x + 1][y + 2] >> LAST);
        d1[x][y + 1] = d2[x - 1][y + 2];
    }
}
\end{verbatim}
\subsection{CUDA}

Never implemented.

\section{FHP Model}

\subsection{Collision}

\subsection{Propagation}

\subsection{The Code}

\subsection{Optimizations}

\subsection{CUDA}

\section{FCHC Model}

In the \emph{Face Centered Hyper Cubic} model, you have 24 vectors.
16 correspond to the vertexes of the hyper cube while the other 8
correspond to the center of the squares in the hyper cube.

\subsection{Reduction}

Since you cannot store the entire table, you consider only states
that have a normalized momentum. The coordinates of the momentum
$q$ are defined by

\$$q_\alpha = \sum_{i = 0}^{24} S_i C_{i \alpha}$\$

where $S_i$ is a boolean variable signifying the presence or
absence of a particle with velocity $C_i$. The momentum is
normalized if the coordinates satisfy

\$$q_1 > q_2 > q_3 > q_4 > |q_1 - q_2 - q_3|$\$

Different isometries are then applied if one of the following
conditions are satisfied in the following sequence:

\begin{itemize}
\item 
  if $q_1 < 0$, then $S_1$ is applied to $q$
\item 
  if $q_2 < 0$, then $S_2$ is applied to $q$
\item 
  if $q_3 < 0$, then $S_3$ is applied to $q$
\item 
  if $q_4 < 0$, then $S_4$ is applied to $q$
\item 
  if $q_1 < q_2$, then $P_{12}$ is applied to $q$
\item 
  if $q_3 < q_4$, then $P_{34}$ is applied to $q$
\item 
  if $q_1 < q_3$, then $P_{13}$ is applied to $q$
\item 
  if $q_2 < q_4$, then $P_{24}$ is applied to $q$
\item 
  if $q_2 < q_3$, then $P_{23}$ is applied to $q$
\item 
  if $q_1 + q_4 < q_2 + q_4$, then $\Sigma_1$ is applied to $q$
\item 
  if $q_1 < q_2 + q_3 + q_4$, then $\Sigma_2$ is applied to $q$
\end{itemize}
Where

\begin{itemize}
\item 
  $S_\alpha$ is the symmetry with respect to the plane $x_\alpha = 0$
\item 
  $P_{\alpha \beta}$ is the symmetry with respect to the plane
  $x_\alpha = x_\beta$
\item 
  $\Sigma_1$ is the symmetry with respect to the plane
  $x_1 + x_4 = x_2 + x_3$
\item 
  $\Sigma_2$ is the symmetry with respect to the plane
  $x_1 = x_2 + x_3 + x_4$
\end{itemize}
\end{document}
