\documentclass{article}

\title{Implementation of a Lattice Gas Simulation}
\author{Tony H. Kim}
\date{Spring 2007}

\usepackage{amsmath}
\usepackage{graphicx}


\setlength\topmargin{0in}
\setlength\headheight{0in}
\setlength\headsep{0in}
\setlength\textheight{9in}
\setlength\textwidth{6.5in}
\setlength\oddsidemargin{0in}
\setlength\evensidemargin{0in}
\setlength\parindent{0in}
\setlength\parskip{0.1in}

\begin{document}
\maketitle

\section{Introduction}

In this paper I discuss the implementation of a lattice-gas simulation on a hexagonal grid. The application is written in the C++ language using the Windows and DirectX APIs. Due to the nature of the project, following the initial introduction to the lattice gas concept, the discussion has somewhat the feel of a software manual. After documenting my application in Section 3, I spend some time discussing the technical difficulties associated with attempting to simulate a large number of (albeit simple) automata; and the strategies employed in order to reduce the cost of computations. In the last section are screens and comments from the various demos I have produced using the engine.

\section{The Lattice Gas}

\subsection{A Simple Description}

\begin{figure}
  \centering
    \includegraphics[width=0.5\textwidth]{intro.png}
  \caption{The lattice gas system as visualized in my application}
  \label{fig:mylattice}
\end{figure}

The lattice gas particle is an entity that hops from point to point on the discrete sites of a lattice with each time step (Figure \ref{fig:mylattice}). These spatial and temporal discretizations imposed on the lattice gas particle make the model a highly artificial system to begin with, but there are additional simplifications that we further assume. For instance:

\begin{itemize}
	\item Owing to the hexagonal grid, at each lattice site (or ``node'') the particle can move in only one of six directions.
	\item There cannot exist more than one particle at each node heading in the same direction. This implies that each node can support up to a maximum of six particles at a particular instant.\footnote{Physics joke: the lattice gas particle must have half-integer spin.}
	\item The magnitude of the particle velocity is fixed.
	\item A particle cannot remain stationary. (Note: This is relaxed in more involved implementations which may explicitly allow for stationary particles. However, the restrictions presented here correspond to the version I have produced.)
\end{itemize}

\begin{figure}
  \centering
    \includegraphics[width=0.9\textwidth]{von_karman.png}
  \caption{At large scales, the lattice gas produces the qualitative features of real flows}
  \label{fig:vonkarman}
\end{figure}

Given these severe, unphysical assumptions of the lattice gas model, it is questionable whether we are actually ``simulating'' anything -- let alone the complex interactions and behaviors of real gases and fluids. However, despite the alarmingly simple framework, it is remarkable that at large enough scales the lattice gas reproduces some of the features associated with real flows. A well-known example of this is shown above in the lattice gas reproduction of the von Karman vortex street phenomenon (Figure \ref{fig:vonkarman}). This and other surprising successes have led some to characterize the lattice gas concept as one that ``has nearly nothing in common with real fluids except for one special property -- at a macroscopic scale it flows just like them!''\cite{rothman1}

\subsection{Advantages of the lattice gas}

If our goal is to attempt the simulation of physical particles, one may reasonably ask: given the performance of modern computers and special-purpose machines, why not simply create a more realistic application that models the interactions of particles in continuous space, and continuous time?\footnote{As ``continuous'' as they can be, implemented in a computer program.} While straightforward, the computational order of growth associated with such an approach prevents simulations at a size large enough to be interesting. 

Intuitively, the problem scales as $O(N^2)$, where N denotes the number of particles being simulated. At every frame, for each of the N particles, we must query $N-1$ others for a possible collision. This is an underestimate, however, for we must take into account three-particle collisions, reflections off of walls, and other auxiliary calculations. Of course, we have not yet even considered the momentum computation itself, which, for just two spherical particles in the plane, will involve many parameters such as the angle of incidence or the size of the particles in addition to the initial velocities.

Hence, to feasibly simulate a large number of particles, we must be willing to make simplifications to the physical model. The lattice gas represents one extreme set of such simplifications. In particular, it is of interest because of its unusual $O(n)$ order of growth in the computational effort. As we will find later, this owes to the fact that the relevant parameter for the ``size'' of the lattice gas simulation is \emph{not} the number of particles present, but rather the number of nodes; and that these nodes, in turn, have (more-or-less) \emph{fixed} neighbors. We will return to this point in Section 2.7, after discussing the means to describe analytically the lattice gas system.

\subsection{The analytic description of the lattice network}

In order to construct the lattice gas, we begin with a coordinate plane. Consider a node at the vector position $\vec{x}_0$ as in Figure \ref{fig:nodegeneration}(a). For simplicity, we require that there exist only one node per location on the plane. We may then use the vector $\vec{x}_0$ to refer to the \emph{node} at the location $\vec{x}_0$. This is convenient notation, and I will be utilizing it throughout the discussion.

Once we have a single node at $\vec{x}_0$, we may construct the rest of the lattice \emph{with respect to} $\vec{x}_0$ with the help of ``displacement'' vectors $\vec{c}_i$ whose directions in the plane are given by:

\begin{center}
	$\vec{c}_i = (cos(\pi i/3),sin(\pi i/3))$, where $i = 0, 1, 2, 3, 4, 5;$
\end{center}

In addition, the magnitude of $\vec{c}_i$ denotes the desired spacing between adjacent nodes of the lattice. Using this construct, we may place a node at $\vec{x}_0+\vec{c}_0$ as in Figure \ref{fig:nodegeneration}(b).

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{node_gen.png}
  \caption{Creating nodes in the coordinate plane}
  \label{fig:nodegeneration}
\end{figure}

Having repeating the above process for each direction $i$, we expect the node $\vec{x}_0$ to be connected to the six nodes at locations $\vec{x}_0+\vec{c}_i$. (see Figure \ref{fig:nodegeneration}(c).) This is so because the ``connections'' of the lattice are made based on the relative positions of the nodes on the coordinate plane.

While the discussion above may be trivial, it is still of some importance for my application, because the process of node generation described here is actualy quite versatile; and it is \emph{exactly} the same procedure that the program follows in creating the node networks. For more information on how to procedurally generated nodes, please consult the files NodeGeometry.h/.cpp.

\subsection{Denoting particles on the lattice}

As mentioned previously, each node can support up to six particles, heading in each direction. More specifically, this is the case when the node actually has six neighbors, as $\vec{x}_0$ in the previous diagrams. In general, if a node has $0 \leq n \leq 6$ neighbors, then it can support upto $n$ particles.

The restriction of one particle per one direction of a site gives rise to a convenient notation for denoting the occupancy of a node. Let $\vec{n}(\vec{x})$ be a six-component vector function of a node $\vec{x}$ whose components $n_i(\vec{x})$ denote the presence ($n_i = 1$) or absence ($n_i = 0$) of a particle heading in the i-th direction.

So, for the example of $\vec{x_0}$ in Figure \ref{fig:nodeoccupation}, the occupation vector may be written:

\begin{figure}
  \centering
    \includegraphics[width=0.4 \textwidth]{particles.png}
  \caption{A possible particle occupation for the node $\vec{x}_0$}
  \label{fig:nodeoccupation}
\end{figure}

\begin{equation}
	\vec{n}(\vec{x}_0) = (n_0,n_1,n_2,n_3,n_4,n_5) = (1, 1, 0, 0, 1, 0)
\end{equation}

Note that each of the entries $n_i$ are ``Boolean'' variables in the sense that they possess only two possible values. The fact that such little information is sufficient to capture the complete state of the system owes to the remarkable simplicity of the lattice gas model. Later, we exploit this simplicity to achieve a very efficient method for collision calculations.

\subsection{The ``equation of motion'' for the lattice gas}

The state of the system is specified for some particular time if $\vec{n}$ is given for every node $\vec{x}_j$ present in the lattice. Given this information, how may we produce the next time frame? The answer to this question is captured in the following equation:

\begin{equation}
	n_i(\vec{x}+\vec{c}_i,t+1) = n_i(\vec{x},t) + \Delta_i[\vec{n}(\vec{x},t)]
	\label{eqn:evolutioneqn}
\end{equation}

As defined previously, the first two terms denote the presence or absence of a particle heading in a particular direction $i$; they are evaluated at adjacent nodes $\vec{x}$ and $\vec{x}+\vec{c}_i$ and at successive times. The second term on the right $\Delta_i$ represents the ``collision operator'' that acts on the complete state of a node occupation at some time. The momentum operator can output a range of values {-1,0,1} and hence can modify the ``trajectory'' of a particle.

To illustrate this point, consider the simple case below:

Suppose that at time $t$, the node at $\vec{x}_0$ posseses a sole particle heading in the $i=1$ direction. Then we may ask whether there will be a particle at the node $\vec{x}_0+\vec{c}_1$ at time $t+1$ heading in the $i=1$ direction. Physical intuition (Another reminder that we are in fact simulating a physical system!) would indicate that the particle would keep moving along its path; i.e. that $n_1(\vec{x}_0+\vec{c}_1) = 1$. We can see how this result follows from the dynamical equation. For $i=1$ we have from Eq. \ref{eqn:evolutioneqn}:

\begin{equation}
	n_1(\vec{x}+\vec{c}_1,t+1) = n_1(\vec{x},t) + \Delta_1[\vec{n}(\vec{x},t)]
\end{equation}

Presumably, for a single particle, there will be no ``collision'' and the collision operator will yield $\Delta_1[\vec{n}(\vec{x},t)] = 0$. This in turn gives:
\begin{eqnarray*}
	n_1(\vec{x}+\vec{c}_1,t+1) &=& n_1(\vec{x},t)\\
														 &=& 1
\end{eqnarray*}

In other words, the particle does in fact continue on its path.

In contrast, suppose now we have the following three-body heads-on collision scenario (Figure \ref{fig:threebody}), i.e.:

\begin{figure}
  \centering
    \includegraphics[width=0.4 \textwidth]{threebodycollision.png}
  \caption{A possible particle occupation for the node $\vec{x}_0$}
  \label{fig:threebody}
\end{figure}

\begin{equation}
	\vec{n}(\vec{x}) = (1,0,1,0,1,0)
\end{equation}

In such a case, the collision operator will yield $\Delta_1[\vec{n}(\vec{x},t)] = -1$, indicating that the particle does \emph{not} continue on its original path in the case of a three-body heads-on collision as shown. (Since now: $n_1(\vec{x}_0+\vec{c}_1) = 0$) Simply put, the collision operator $\Delta_i$ is a mapping from the set of occupation states at a node (all the possible values of $\vec{n}$) to the integers {-1,0,1} that modify the trajectories of particles.

\subsection{More on the collision operator}

In the lattice gas approach, the collision operator $\Delta_i$ alone incorporates the possible ways in which particles interact with one another. The particular form of the collision operator chosen differs among implementations, and determines the complexity of each model. In my application, I employ the simplest collision operator, that handles two-body and three-body collisions.\footnote{There were slight modifications made to allow for two distinguishable types of particles.} This form of the collision operator is known as the FHP-I model (Frisch, Hasslacher, Pomeau).\cite{rothman1}

\subsection{Remark on the computational order of growth}

Consider for a moment a more physical simulation of N-particle dynamics based on Newton's laws of motion and gravitation. For that approach, the relevant equation of motion in two-dimensions is:

\begin{equation}
	\ddot{\vec{x_j}} = \sum_{k=1}^N(G m_k/(r_k-r_j)^2)
\end{equation}

for the j-th particle.

An important feature to note is that, here, the relevant index $j$ is over the \emph{particles}. Hence, the number of particles being simulated is the relevant measure of the ``size'' of the simulation.

In contrast, for the lattice gas, the analogous equation reads:

\begin{equation}
	n_i(\vec{x_j}+\vec{c}_i,t+1) = n_i(\vec{x_j},t) + \Delta_i[\vec{n}(\vec{x_j},t)]
\end{equation}

which is to be computed over all \emph{nodes} $\vec{x}_j$. Hence, the size parameter in this case is the number of nodes in the lattice. It is amusing to note that virtually the same number of computations will be needed to simulate a lattice of \emph{no particles} as it would to simulate a completely saturated lattice!

Finally, note that in the case of the Newtonian simulation, the number of terms on the right-hand side increases proportionally to the size of the problem. However, in the case of the lattice gas, the computations on the RHS are fixed in number. Without detailed analysis, this is the main observation that leads us to conclude that the order of growth in the lattice gas simulation is $O(n)$, rather than $O(n^2)$.

\section{Program documentation}

\subsection{Feature list}

\begin{figure}
  \centering
    \includegraphics[width=0.9\textwidth]{features.png}
  \label{fig:features}
\end{figure}

In developing the application, I have kept the following objectives and features in mind:

\begin{itemize}
	\item Ability to generate and modify the lattice network at runtime. This allows for great control in creating customized ``terrains'' for the particles. (Figure \ref{fig:features}(a))
	\item Ability to simulate at least $N > 1000$ particles at a reasonable speed.\footnote{My benchmark is my own laptop - 1.6GHz, 256MB memory and a practically nonexistent graphics card!} (Figure \ref{fig:features}(b))
	\item Simulate two distinctly colored particles. (Figure \ref{fig:features}(b))
	\item Implementation of various behaviors at the boundaries.
\end{itemize}

\subsection{Program control}

Interaction with the software requires both the keyboard and the mouse. The commands are as follows.

\textbf{Mouse:}
\begin{itemize}
	\item Left-click: Select (or unselect) a node;
	\item Right-click: Create a new node in the clicked location;
	\item Mouse-scroll: Zoom in/out;
\end{itemize}

\textbf{Keyboard:}
\begin{itemize}
	\item W,S,A,D: Move the camera up/down/left/right respectively if no node is selected; If a node is selected, moves the location of the selected node.
	\item R,F: Zoom the camera in/out;
	\item Q: Delete the selected node;
	\item C: Clear all particles from the nodes;
	\item V: Randomly fill the node network with particles;
	\item 1,2,3,4,5,6,0: If a node is selected, then toggles a blue particle in the particular direction. ($i=6 \Leftrightarrow i=0$). If left-shift is held while the number is pressed, then will toggle a \emph{red} particle;
	\item Space-bar: Unselect a node
\end{itemize}

\subsection{Code organization}

The following briefly describes the source files  that comprise my application and their dependencies:

\begin{figure}
  \centering
    \includegraphics[width=0.9\textwidth]{code_layout.png}
\end{figure}


\subsection{Further details on the development environment}

The program was written using the C++ language of the Microsoft Visual Studio package, using the Windows and DirectX APIs. (So I'm quite confident the code will not run on a Mac or on Athena.) Furthermore, I have used a very recent release of the DirectX software development kit (SDK February 2007) so that the runtime software may need to be updated in order to run the demonstration executables.

In addition, on some machines I have encountered a ``d3dx9\_32.dll missing'' error. This issue is discussed at and can be resolved by consulting:

\begin{verbatim}
  http://www.toymaker.info/Games/html/d3dx_dlls.html
\end{verbatim}

But I understand if you would not want to install additional software just to get my lattice simulations to run. (That's why I've included Section 5!)

\section{Implementation details}

\subsection{Runtime generation of the node network}

Having decided on runtime lattice generation, the pertinent programming challenge is: How does a node, acting independently, find its neighbors and make connections?

To resolve this issue, we must turn to the underlying coordinate system since that is the only framework that we may rely on initially. My solution relies on a data structure dubbed the ``NodeManager'' whose main purpose is to rapidly retrieve, for a given input coordinate $(x_0,y_0)$, the node that is sufficiently closest to that target.

\begin{figure}
  \centering
    \includegraphics[height=300px]{node_manager.png}
  \caption{Representation of the ``NodeManager'' data structure}
  \label{fig:node_manager}
\end{figure}

As seen in Figure \ref{fig:node_manager}, the NodeManager is comprised of two vector and list containers. The actual node data -- which is represented as a C struct (in node.h) -- is housed in the list. Throughout the execution, the list of node data lies dormant, in the sense that its location in memory is relatively constant. The primary operations that we perform on this list of ``node-data'' are addition and deletion, both of which are well suited for the list.\cite{dray}

The two vectors, on the other hand, contain pointers to the elements of node-data; and are sorted by the properties of the pointed data. One vector is sorted by the x-coordinate of the pointees; and the other by the y-coordinate. The two vectors are kept ordered at all times through careful insertions and modifications.

The purpose of this construct is the following. Suppose we have a collection of N nodes at arbirary points of the coordinate plane, and were interested in the node that is located at the point $(x_0,y_0)$ within some tolerance. Without proper design, to answer this question, we would have to test all N nodes for their distances from the target point. The scaling properties of this solution isn't bad; however my solution is far superior.
 
The approach is the following:
\begin{enumerate}
	\item For some tolerance $t$,
	\item Determine the subset of the x-pointers whose pointees have a x-coordinate within the range $[x_0-t,x_0+t]$. This occurs in $O(log(n))$, since we take advantage of the binary search capabilites of the vector.
	\item Determine the subset of the y-pointers whose pointees have a y-coordinate within the range $[y_0-t,y_0+t]$. This occurs in $O(log(n))$.
	\item Find the intersection of the two resulting sets. This gives the pointers of all nodes whose x- and y-coordinates lie within the prescribed tolerance of the input target. While the intersection operation proceeds at $O(n)$, the previous subset operations have drastically reduced the size of the candidate nodes, so that we may consider this step as taking place in nearly constant time.\footnote{I simply mean to say, that as long as the particle distribution is not becoming more \emph{dense} with respect to the tolerance $t$, the results of the previous two operations should give subsets whose sizes are relatively independent of the actual number of nodes present in the lattice.}
	
	\item Cycle through the elements of the intersection and find the node that is closest to the target. Again, since the size of this set is so small in comparison to the entire collection of nodes, we can consider this to scale at $O(1)$.
\end{enumerate}

With this structure we can determine whether there are nodes at a particular coordinate location in effectively $O(log(n))$ time. This is ideal performance and allows my application to rely on the NodeManager extensively even for very large lattices. Hence, the runtime node network generation revolves around each node querying the NodeManager for nodes in their vicinity. (More specifically, at the locations $\vec{x}+\vec{c}_i$ where $\vec{x}$ denotes its own position.)

\subsection{Rapid collision calculations}

While an efficient algorithm for establishing node neighbors is helpful, it is in the overall operation of the program only a ``minor'' issue since it is invoked infrequently; e.g. it is called only when the network is first generated or modified. As implemented, each node, after establishing its neighbors, records the relevant identifiers so that it does not require the NodeManager apparatus in the future for later computations.

On the other hand, for \emph{every} frame we must perform the collision computations as described in Section 2.5. Hence, \emph{this} is the process that needs to be optimized for maximum performance.

The lattice gas simulation is particularly well-suited to implementation in the C language because C allows manipulations of individual bits in memory. Because the state of the lattice gas system can be represented by a sequence of \emph{booleans} i.e. 0 or 1 (recall that knowledge of $\vec{n}(\vec{x}_j)$ for all nodes j gives the complete state of the system), a natural method to represent each node occupancy information is to use the six-bits of a byte. Unfortunately, because the natural ``chunk'' of memory in most systems is a byte (8-bits) I will simply leave the two highest-order bits unused. This is rather unfortunate since I cannot then rely on built-in bit-manipulation operations which assume the full 8-bit width. (On the other hand, many of the built-in bit operations are machine-dependent and are not good to rely on anyway.)

Consider again the following node geometry and occupation. For the node at $\vec{x}_0$ we may write the occupancy as follows:

\section{Gallery}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{511.png}
  \caption{Experiment 1: The particles are initially concentrated in the left container}
  \label{fig:511}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{512.png}
  \caption{Experiment 1: The particles have diffused evenly throughout both containers}
  \label{fig:512}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{521.png}
  \caption{Experiment 2: Two distinguishable sets of particles begin from opposite containers}
  \label{fig:521}
\end{figure}


\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{522.png}
  \caption{Experiment 2: It takes a considerable amount of time for the two sets to mix to any extent. (Compared to the diffusion time of experiment 1)}
  \label{fig:522}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{531.png}
  \caption{Experiment 3: Attempted to see if interesting boundaries form along the edges of a pin-like obstacle in the way of laminar flow}
  \label{fig:531}
\end{figure}


\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{541.png}
  \caption{Experiment 4: A group of blue particles with net momenta to the right collides with a collection of red particles whose velocities are distributed randomly.}
  \label{fig:541}
\end{figure}


\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{542.png}
  \caption{Experiment 4: The aftermath of such a collision}
  \label{fig:542}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{551.png}
  \caption{Experiment 5: A well-ordered alternating-color rows of particles all move to the right towards obstacles}
  \label{fig:551}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{552.png}
  \caption{Experiment 5: The resulting mess after the main body makes past the obstacle.}
  \label{fig:552}
\end{figure}

\begin{thebibliography}{9}

				\bibitem{dray}
				In developing the NodeManager algorithm used in the application, I benefited greatly from discussions with high school friend and current CS/math student at the University of British Columbia, Karl Dray. (karldray@interchange.ubc.ca)

				\bibitem{Josuttis}
				Josuttis, Nicolai M.,
				\emph{The C++ Standard Library}.
				Addison-Wesley, Indianapolis,
				1999

				\bibitem{luna}
				Luna, Frank,
				\emph{Introduction to 3D Game Programming with DirectX 9.0c}.
				Wordware Publishing, Inc. Plano, Texas,
				2006

        \bibitem{rothman1}
          Rothman and Zaleski,
          \emph{Lattice-Gas Cellular Automata}.
          CUP, Cambridge,
          1997.

			
				\bibitem{prata}
				Prata, Stephen,
				\emph{C Primer Plus}.
				Sams Publishing, Indianapolis,
				2005
				
\end{thebibliography}

\end{document}