\chapter{Synchronous Applications I}
\label{synch1}

\section{QCD and the Beginning of C$^3$P}
\label{begin}

The Caltech Concurrent Computation Project started with QCD,\index{QCD} or 
Quantum Chromodynamics, as its first application.  QCD is discussed in more
detail in Section~\ref{apps} and Section~\ref{qcd}, but here we will put
in the historical perspective.  This nostalgic approach is developed in
\cite{Fox:87d}, \cite{Fox:88oo} as well as Chapter~\ref{intro} and
Chapter~\ref{capri} of this book.

We show in Table~\ref{tab1.qcd}, fourteen QCD simulations, labeled by
representative physics publications, performed within C$^3$P using
parallel machines.  This activity started in 1981 with simulations, using
the first four node 8086-8087 based prototypes of the Cosmic Cube, that
were quite competitive in performance with the VAX~11/780 on which we
had started in 1980 our computational physics program within high energy
physics at Caltech.  The 64-node Cosmic Cube was used more or less
continuously from October, 1983 to mid 1984 on what was termed by
Caltech, a ``Mammoth Calculation'' in the press release shown in
Figure~\ref{fig2.begin}.  This is the modest $12^3 \times 16$
four-dimensional lattice calculation reported in line~3 of
Table~\ref{tab1.qcd}.  As trumpeted in Figure~\ref{fig1.begin} and
Figure~\ref{fig2.begin}, this was our first major use of parallel
machines and a critical success on which we built our program.

Our 1983--84 calculations totaled some 2,500 hours on the 64-node
Cosmic Cube and successfully competed with 100 hour CDC Cyber~205
computations that were the state of the art at the time
\cite{Barkai:84b}, \cite{Barkai:84c}, \cite{Bowler:85a},
\cite{DeForcrand:85a}.  We used a four-dimensional lattice size of $12^3
\times 16$ for 27,648 grid points where on each grid point defined eight
gluon field values.  The resultant 221,184 degrees of freedom seems
modest today as QCD practitioners contemplate lattices of order $100^4$
simulated on machines of teraflop performance \cite{Aoki:91a}.  However,
this lattice was comparable to those used on vector supercomputers at
the time.

Hallmarks of this work was the interdisciplinary team building hardware,
software and parallel application.  Further, from the start we stressed
large supercomputer level simulations where parallelism will make the
greatest initial impact.  It was also worth noting that our use of
comparatively high level software paid off---Otto and Stack were able to
code better algorithms \cite{Parisi:83a} than the competing vector
supercomputer teams.  The hypercube could be programmed conveniently
without use of microcode or other unproductive environments needed on
some of the other high performance machines of the time.

Our hypercube calculations used an early C plus message passing
programming approach which later evolved to the {\em Express\/} system
described in the next chapter.  Although not as elegant as data parallel
C and Fortran discussed in Chapter~\ref{data}, our approach was easier
than hand coded assembly that was quite common for alternative high
performance systems of the time.

Figure~\ref{fig1.begin} and Figure~\ref{fig2.begin} shows extracts from
Caltech and newspaper publicity of the time.  We were only a collection
of 64 IBM-PC's.  Was that a good thing (as we thought) or an indication
of our triviality (as a skeptical observer thought).  1985 saw the start
of a new phase as conventional supercomputers and availability increased
in power and NSF and DOE allocated many tens of thousands of CRAY~X-MP
(2, Y-MP) and ETA-10 time to QCD simulations.  Our final QCD hypercube
calculations in 1989 within C$^3$P used a 64-node JPL Mark~IIIfp with
approximately three hundred megaflops performance.  Since this work, we
switched to using the Connection Machine CM-2, which by 1990 was the
commercial standard in the field.  C$^3$P helped the Los Alamos group of
Brickner and Gupta (one of our early graduates!) to develop the first
CM-2 QCD codes, which now (1991) perform at six gigaflops on the full
size $\rm 64\,K$ CM-2.

\begin{figure}
\vspace*{2.0in}
\caption{Caltech Scientists Develop `Parallel' Computer Model}
\label{fig1.begin}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Caltech's Cosmic Cube Performing Mammoth Calculations}
\label{fig2.begin}
\end{figure}

It is not surprising that our first hypercube calculations in C$^3$P did
not need the MIMD structure of the machine.  This was also a
characteristic of Sandia's pioneering use of the 1024 node nCUBE-1
\cite{Gustafson:88a}.  Synchronous applications like QCD are
computationally important and have a simplicity that made them the
natural starting point for our project.

\section{Synchronous Applications}
\label{apps}

Table~\ref{tab1.begin} indicates that 70\% of our first set of applications
were of the class we call Synchronous.  As remarked above, this could be
expected in any such early work as these are the problems with the
cleanest structure that are, in general, the simplest to code and, in
particular, the simplest to parallelize.  As already defined in
Section~\ref{props}, Synchronous applications are characterized by a
basic algorithm that consists of a set of operations which applied
identically in every point in a data set.  The structure of the problem
is typically very clear in such applications, and so the parallel
implementation easier than for the irregular loosely synchronous and
asynchronous cases.  Nevertheless, as we will see, there are many
interesting issues in these problems and they include many very
important applications.  This is especially true for academic
computations that address fundamental science.  These are often
formulated as studies of fundamental microscopic quantities such as the
quark and gluon fundamental particles seen in QCD of Section~\ref{qcd}.
Fundamental microscopic entities naturally obey identical evolution laws
and so lead to synchronous problem architectures.  ``Real world
problems''---perhaps most extremely represented by the battle
simulations of Chapter~\ref{pci} in this book---typically do not involve
arrays of identical objects, but rather the irregular dynamics of
several different entities.  Thus, we will find more loosely synchronous
and asynchronous problems as we turn from fundamental science to
engineering and industrial or government applications.  We will now
discuss the structure of QCD in more detail to illustrate some general
computational features of synchronous problems.

{\small
\begin{table}
\begin{center}
\begin{tabular}{|lllrl|}
\hline
\multicolumn{1}{|c}{Field}&\multicolumn{1}{c}{Project}&
\multicolumn{1}{c}{Reference}&\multicolumn{1}{c}{Section}&
\multicolumn{1}{c|}{Problem}
\\
&\multicolumn{1}{c}{Leader}&&\multicolumn{1}{c}{Number}&
\multicolumn{1}{c|}{Architecture}
\\
\hline
\medStrut
Lattice Gauge Theory&S. Otto&\cite{Otto:83b},\cite{Otto:84d}&\ref{qcd}&
Synchronous \\
&&&& \\
Acoustic Wave Equation&R. Clayton&\cite{Clayton:88a}&\ref{apps}&
Synchronous \\
&&&& \\
Coulomb Gas&F. Fucito&\cite{Fucito:85a}&\ref{apps}&Synchronous \\
Monte Carlo&S. Solomon&\cite{Fucito:85f}&& \\
&&&& \\
Two-Dimensional&D. Meier&\cite{Meier:88a}&\ref{apps}&Synchronous \\
Hydrodynamics &&& \\
&&&& \\
Chemical Reaction&P. Hipes&\cite{Hipes:88a},\cite{Hipes:88b}&\ref{banded}&
Synchronous \\
Dynamics&A. Kuppermann&\cite{Wu:90a}&& \\
&&&& \\
Evolution of Universe&J. Salmon&\cite{Salmon:84a}&\ref{nbody}&Synchronous \\
using FFT&&\multicolumn{2}{l}{(more advanced algorithm)} & \\
&&&& \\
Dynamics of&B. Werner&\cite{Werner:87a}&\ref{micro}&Loosely \\
Sand Grains&P. Haff&\cite{Werner:88a}&&Synchronous \\
&&&& \\
Traveling Salesman&S. Otto&\cite{Felten:85b}&\ref{tsp}&Synchronous \\
Optimization by&E. Felten&&& \\
Simulated Annealing&S. Karlin&&& \\
&&&& \\
Sorting&E. Felten&\cite{Felten:85a}&\ref{sort}&Loosely \\
&S. Karlin&&&Synchronous \\
&S. Otto&&& \\
&&&& \\
Melting of a Two&M. Johnson&\cite{Johnson:86a}&\ref{melting}&Asynchronous \\
Dimensional Solid&&&& \\
&&&& \\
\hline
\end{tabular}
\caption{The Ten Pioneer Hypercube Applications}
\label{tab1.begin}
\end{center}
\end{table}
}

The applications using the Cosmic Cube were well established by 1984 and
Table~\ref{tab1.begin} lists the ten projects which were completed in the
first year after we started our interdisciplinary project in the summer
of 1983.  All but one of these projects are described more or less in
this book while the other will be found in \cite{Fox:88a}.  They covered
a reasonable range of application areas and formed the base on which we
first started to believe that parallel computing works!
Figure~\ref{fig3.begin} illustrates the regular lattice used in QCD and
its decomposition onto 64 nodes.  QCD is a four-dimensional theory and
all four dimensions can be decomposed.  In our initial 64-node Cosmic
Cube calculations, we used the three-dimensional decompositions show in
Figure~\ref{fig3.begin} with the fourth dimension, as shown in
Figure~\ref{fig4.begin}, and internal degrees of freedom stored
sequentially in each node.  Figure~\ref{fig3.begin} also indicates one
subtlety needed in the parallelization\@; namely one needs a so called
red-black strategy with only half the lattice points updated in each of
the two (red and black) phases.  Synchronous applications are
characterized by such a regular spatial domain as shown in
Figure~\ref{fig2.begin} and an identical update algorithm for each
lattice point.  The update makes use of a Monte Carlo procedure
described in Section~\ref{qcd} and in more detail in Chapter~12 of
\cite{Fox:88a}.  This procedure is not totally synchronous as the
``accept-reject'' mechanism used in the Monte Carlo procedure does not
always terminate at the same step.  This s no problem on a MIMD machine
and even makes the problem ``slightly loosely synchronous''.  However,
SIMD machines can also cope with this issue as all systems (DAP, CM-2,
Maspar) have a feature that one can allow processors either to execute
the common instruction or ignore it.

\begin{figure}
\vspace*{2.0in}
\caption{A $12^3$ Problem Lattice decomposed onto a 64-node machine
arranged as a $4^3$ machine lattice.  Points labeled X (``red'') or
$\bullet$ (``black'') can be updated at the same time.}
\label{fig3.begin}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{The 16 time and eight internal gluon degrees of freedom stored
at each point shown in Figure~\protect\ref{fig2.begin}.}
\label{fig4.begin}
\end{figure}

Figure~\ref{fig5.begin} illustrates the nearest neighbor algorithm used
in QCD and very many problems described by local interactions.  We see
that some updates require communication and some don't.  In the message
passing software model used in our hypercube work described in
Chapter~\ref{cros}, the user is responsible for organizing this
communication with an explicit subroutine call.  Our later QCD
calculations and the Spin Simulations of Section~\ref{spins} use the
data parallel software model on SIMD machines where a compiler can
generate their messaging.  Chapter~\ref{data} will mention later
projects which are aimed at producing a uniform data parallel Fortran or
C compiler which will generate the correct message structure for either
SIMD or MIMD machines on such regular problems as QCD.

\begin{figure}
\vspace*{2.0in}
\caption{Template of a local update involving no communication in a) and
the value to be communicated in b).}
\label{fig5.begin}
\end{figure}

The calculations in Section~\ref{qcd} and Section~\ref{spins} used a
wide variety of machines, inhouse and commercial multicomputers, as well
as the SIMD DAP and CM-2.  The spin calculations in Section~\ref{spins}
can have very simple degrees of freedom including that of the binary
``spin'' of the Ising model.  These are naturally suited to the single
bit arithmetic available on the AMT DAP and CM-2.  Some of the latest
Monte Carlo algorithms do not use the local algorithms of 
Figure~\ref{fig4.begin} but exploit the irregular domain structure seen
in materials near a critical point.  These new algorithms are much more
efficient but very much more difficult to parallelize---especially on
SIMD machines.  They are discussed in Section~\ref{tree}.  We also see a 
taste of the ``embarrassingly parallel'' problem structure of 
Chapter~\ref{indepar} in Section~\ref{spins}.  For the Potts simulation, we 
obtained parallelism not from the data domain (lattice of spins) but from 
different starting points for the evolution.  This approach described in more 
detail in Section~\ref{dynamic}, would not be practical for QCD with many 
degrees of freedom as one must have enough memory to store the full lattice 
in each node of the multicomputer.

Table~\ref{tab1.begin} lists the early seismic simulations of the group
led by Clayton.  These solved the elastic wave equations using finite
different methods as discussed in Section~\ref{spatial}.  The equations
are iterated with time steps replacing the Monte Carlo iterations used
above.  This work is described in Chapter~\ref{cros} and
Chapter~\ref{indepar} of \cite{Fox:88a} and represents methods that can
tackle quite practical problems such as predicting the response of
complicated geological structures such as the Los Angeles Basin.  The
two-dimensional hydrodynamics work of Meier \cite{Meier:88a} is
computationally similar using the regular decomposition and local update
of Figure~\ref{fig3.begin} and Figure~\ref{fig5.begin}.  These
techniques are now very familiar and may seem ``old-hat''.  However, it
is worth noting that as described in Chapter~\ref{data} we are only now
(in 1992) developing the compiler technology that will automate these
methods developed ``by-hand'' by our early users.  A much more
sophisticated follow-on to these early seismic wave simulations is the
ISIS interactive seismic imaging system described in Chapter~\ref{complex}.

Chapter~9 of \cite{Fox:88a} explains the well known synchronous
implementation of long range particle dynamics.  This algorithm was not
used directly in any C$^3$P application as we implemented the much
better cluster algorithms described in Sections~\ref{nbody},\ref{vortex}, and 
Section~\ref{tree}.  As explained in \cite{Fox:88a}, one can use the
long range force algorithm in any case where the calculation involves a
set of $N$ points with observables requiring functions of every pair of
which there are $N^2$.  In the language of Chapter~\ref{methodology},
this problem has a system dimension of one, whatever its geometrical
dimension.  This is illustrated in Figure~\ref{fig6.begin} which
represents $f_C$ in the form of Equation~\ref{eqn10.method} and
Equation~\ref{eqn13.method}.  We find $1/\beta\ = d_{\rm num} = 2$ for
the simple two-dimensional decompositions described for the Clayton and
Meier applications for Table~\ref{tab1.begin}.  We increase range of
``interaction''---defined formally by
\begin{eqnarray}
\mbox{Calculation at}\ \underline{x} = & \mbox{function}\  
(\underline{x},\underline{y}) \nonumber \\
& \mbox{for all}\ \left|\underline{x} - \underline{y}\right| < R 
\end{eqnarray}
from small (nearest neighbor) $R$ to $R=\infty$, the long range force.
$f_C$ decreases as $R$ increases with $\beta = d_{\rm num} = 1$
independent from $d_{\rm geom}$ for $R=\infty$.  Noederlinger
\cite{Lorenz:87a} and Theiler \cite{Theiler:87a}, \cite{Theiler:87b}
used such a ``long-range'' algorithm for calculating the correlation
dimension of a chaotic dynamical system.  This measures the essential
number of degrees of freedom for a complex system which in this case was
a time series becoming a plasma.  The correction function involved
studying histograms of $\left|\underline{x} (t_i) - \underline{x}
(t_j)\right|$ over the data points $\underline{x}\,(t_i)$ at time $t_i$.

\begin{figure}
\vspace*{2.0in}
\caption{Communication overhead as a function of increasing radius of
interaction.}
\label{fig6.begin}
\end{figure}

Fucito and Solomon \cite{Fucito:85b}, \cite{Fucito:85f} studied a simple
coulomb gas which naturally had a long range force.  However, this was a
Monte Carlo calculation that was implemented efficiently by an ingenious
algorithm that cannot directly use the analysis of the particle dynamics
(time stepped) case shown in Figure~\ref{fig6.begin}.  Monte Carlo is
typically harder to parallelize than time evolution, as all
``particles'' can be evolved in time together.  However, Monte Carlo
updates can only proceed simultaneously if they are involved disjoint
particle sets.  This implies the red-black ordering of
Figure~\ref{fig5.begin} and requiring a difficult asynchronous algorithm
in the irregular melting problem of Section~\ref{melting}.  Johnson's
application was the technically hardest in our list of pioneers in
Table~\ref{tab1.begin}.

Finally, Section~\ref{grain} uses cellular automata ideas that lead to a
synchronous architecture to grain dynamics, which if implemented
directly as in Section~\ref{grain} would be naturally loosely
synchronous.  This illustrates that the problem architecture depends on
the particular numerical approach.

\section{Quantum Chromodynamics}
\label{qcd}

\subsection{Introduction}

Quantum Chromodynamics\index{Quantum Chromodynamics} is the proposed theory 
of the so called strong interactions that mediate the forces between protons, 
neutrons, and mesons and so control and formation of nuclei.  Quantum 
Chromodynamics is always abbreviated QCD.  The fundamental properties of QCD 
cannot be directly tested, but a wealth of indirect evidence supports this 
theory.

Gauge theories\index{Gauge theories} are ubiquitous in elementary particle 
physics\@: the electromagnetic interaction between electrons and photons is 
described by quantum electrodynamics (QED) based on the gauge group U(1), the 
strong force between quarks and gluons is believed to be explained by quantum
chromodynamics (QCD) based on SU(3), and there is a unified 
description of the weak and electromagnetic interactions in terms of the 
gauge group $\rm SU(2)\times U(1)$.  The strength of these interactions is 
measured by a coupling constant.  This coupling constant is small for QED, so 
very accurate analytical calculations can be performed using perturbation 
theory.  However, for QCD, the coupling constant appears to increase with 
distance so that perturbative calculations are only possible at short 
distances.  In order to solve QCD at longer distances, Wilson 
\cite{Wilson:74a} introduced lattice gauge theory\index{lattice gauge
theory} in which the space-time continuum is discretized to provide a cut-off 
that regulates ultraviolet divergences non-perturbatively.  This 
discretization onto a lattice, which is typically hypercubic, also makes the 
gauge theory amenable to numerical simulation by computer.

{\small
\begin{table}
\begin{center}
\begin{tabular}{|l|l|c|l|}
\hline
\multicolumn{1}{|c|}{Authors}&\multicolumn{1}{c|}{Project}&
\multicolumn{1}{c|}{Hypercube}&\multicolumn{1}{c|}{Reference}
\\
\hline
\medStrut
Brooks, Fox, Otto, & Glueball mass (SU(2)) & Cosmic Cube &
\protect\cite{Brooks:83a} \\
Randeria, Athas, & $4^3 \times 8$ lattice & 4 node & \\
De Benedictis, & & & \\
Newton, Seitz & & & \\
\hline
Otto, Randeria & Glueball mass, modified action & Cosmic Cube &
\protect\cite{Otto:83a} \\
& $4^3 \times 8$ lattice & 4 node & \\
\hline
Otto, Stack & Heavy quark potential & Cosmic Cube & \protect\cite{Otto:84a}
\\
& $12^3 \times 16$ lattice & 64 node & \\
\hline
Otto, Stolorz & Glueball mass, enhanced statistics & Cosmic Cube &
\protect\cite{Otto:85b} \\
& $12^3 \times 16$ lattice & 64 node & \\
\hline
Patel, Otto, Gupta & Monte Carlo renormalization group, & Cosmic Cube &
\protect\cite{Patel:85a} \\
& nonperturbative $\beta$-function & 64 node & \\
\hline
Fucito, Soloman & Chiral symmetry breaking, mass & Mark II &
\protect\cite{Fucito:84a}, \protect\cite{Fucito:85b}, \\
& spectrum, deconfinement transition, & 64 node &
\protect\cite{Fucito:85c}, \protect\cite{Fucito:85d} \\
& pseudo-fermions & & \protect\cite{Fucito:86a} \\
\hline
Flower, Otto & Field distribution & Mark II & \protect\cite{Flower:85a}
\\
& $10^2 \times 12^2$ lattice & 32 node & \\
\hline
Flower, Otto & Heavy quark potential & Mark II &
\protect\cite{Flower:86b} \\
& $20^4$ lattice & 128 node & \\
\hline
Kolawa, & Glueball mass (SU(2)), & Mark II &
\protect\cite{Furmanski:87c} \\
Furmanski & Hamiltonian ``loop'' formalism & 32 node & \\
\hline
Stolorz, Otto & Microcanonical renormalization & Mark I &
\protect\cite{Stolorz:86b} \\
& group (SU(2)) & 64 node & \\
\hline
Flower & Baryons on the lattice $20^4$ & Mark II &
\protect\cite{Flower:87e} \\
& & 128 node & \\
\hline
Chiu & Random block lattice & nCUBE-1 & \protect\cite{Chiu:89b},
\protect\cite{Chiu:88c}, \\
& & 1024 node & \protect\cite{Chiu:88e}, \protect\cite{Chiu:88f} \\
& & & \protect\cite{Chiu:89a} \\
\hline
Ding, Baillie, & Heavy quark potential & Mark IIIfp &
\protect\cite{Ding:90b} \\
Fox & $24^3 \times 10$, $24^3 \times 12$, $12^3 \times 24$, & 32 node & \\
& $16^3 \times 24$ & and 64 node & \\
\hline
Baillie, Brickner, & QCD with dynamical Wilson & TMC CM-2 &
\protect\cite{Baillie:89e} \\
Gupta & fermions $16^4$ lattice & 16K & \\
\hline
\end{tabular}
\caption{Quantum Chromodynamic (QCD) Calculations within C$^3$P}
\label{tab1.qcd}
\end{center}
\end{table}
}

Most of the work on lattice gauge theory has been directed towards solving
lattice QCD and, thus, deriving the hadron mass spectrum from first 
principles.  This would confirm QCD as the theory of the strong force.  Other 
calculations have also been performed, in particular, the properties of QCD 
at finite temperature and/or finite baryon density have been determined.  
Unfortunately, in order to simulate lattice QCD on a computer, one must 
integrate out the quark variables (because they are fermions, 
{\em i.e.,}\@ anticommuting elements of a Grassmann algebra rather than 
numbers) leaving a highly non-local fermion determinant for each flavor of 
quark.  Physically, this determinant arises from closed quark loops.  The 
simplest way to proceed is to ignore these quark loops and work in the 
so-called quenched approximation with no flavors of quark present.  (This may 
be valid for heavy quarks.)  Current state-of-the-art quenched QCD 
calculations are performed on lattices of size $24^3\times 48$.

However, to investigate the physically more realistic fully interacting
QCD, with the inclusion of dynamical quarks, we must go beyond the quenched 
approximation and tackle the problem of the fermion determinant.  We shall 
return to this later.

\subsection{Monte Carlo}

In order to explain the computations for QCD, we use the Feynman path 
integral formalism.  For any field theory described by a 
Lagrangian density ${\cal L}$, the dynamics of the fields $\phi(x)$ are 
determined through the action functional
\begin{equation}
S(\phi) = \int dx {\cal L}(\phi).
\end{equation}
In this language, the measurement of a physical observable represented by an 
operator ${\cal O}$ is given as the expectation value
\begin{equation}
\langle {\cal O}(\phi)\rangle = \frac{1}{Z} \int D\phi {\cal O}(\phi) 
e^{-S(\phi)},
\label{eqn2.qcd}
\end{equation}
where the partition function\index{partition function} $Z$ is
\begin{equation}
Z = \int D\phi e^{-S(\phi)}.
\end{equation}
In these expressions, the integral $D\phi$ indicates a sum over all possible
configurations of the field $\phi$.  A typical observable would be a product 
of fields ${\cal O}=\phi(x)\phi(y)$, which says how the fluctuations in the 
field are correlated, and in turn, tells us something about the particles 
that can propagate from point $x$ to point $y$.  The appropriate correlation 
functions give us, for example, the masses of the various particles in the 
theory.  Thus, to evaluate most any quantity in field theories like QCD,
one must simply evaluate the corresponding path integral.  The catch is that 
the integrals are over an infinite-dimensional space.

To put the field theory onto a computer, we begin by discretizing space and 
time into a lattice of points.  Then the functional
integral is simply defined as the product of the 
integrals over the fields at every site of the lattice $\phi(n)$\@:
\begin{equation}
\int D\phi \rightarrow \int \prod_n \lbrack d\phi(n)\rbrack.
\end{equation}
Restricting space and time to a finite box, we end up with a finite (but 
large) number of ordinary integrals, something we might imagine doing 
directly on a computer.  However, the high multi-dimensionality of these 
integrals renders conventional mesh techniques impractical.  Fortunately, the 
presence of the exponential $e^{-S}$ means that the integrand is sharply 
peaked in one region of configuration space.  Hence, we resort to a 
statistical treatment and use Monte Carlo type algorithms to sample the 
important parts of the integration region.

Monte Carlo algorithms typically begin with some initial 
configuration of fields, and then make pseudo-random changes on the fields 
such that the ultimate probability $P$ of generating a particular field 
configuration $\phi$ is proportional to the Boltzmann factor,
\begin{equation}
P(\phi) = e^{-S(\phi)},
\end{equation}
where $S(\phi)$ is the action associated with the given configuration.
There are several ways to implement such a scheme, but for many theories the 
simple Metropolis algorithm\index{Metropolis algorithm} \cite{Metropolis:53a} 
is effective.  In this algorithm, a new configuration $\phi'$ is generated by 
updating a single variable in the old configuration $\phi$ and calculating 
the change in energy (action)
\begin{equation}
\Delta S = S(\phi') - S(\phi).
\end{equation}
If $\Delta S \le 0$, the change is accepted; if $\Delta S > 0$ the change
is accepted with probability $\exp(- \Delta S)$.  In practice, this is done 
by generating a pseudo-random number $r$ in the interval [0,1] with uniform 
probability distribution and accepting the change if $r < \exp(- \Delta S)$.
This is guaranteed to generate the correct (Boltzmann) distribution of 
configurations, provided ``detailed balance'' is satisfied.  That condition 
means that the probability of proposing the change $\phi\rightarrow\phi'$ is 
the same as that of proposing the reverse process $\phi'\rightarrow\phi$.
In practice, this is true if we never simultaneously update two fields which 
interact directly via the action.  Note that this constraint has important 
ramifications for parallel computers as we shall see below.

Whichever method one chooses to generate field configurations, one updates 
the fields for some equilibration time of $E$ steps, and then calculates the 
expectation value of ${\cal O}$ in equation~\ref{eqn2.qcd} from the next $T$ 
configurations as
\begin{equation}
\langle {\cal O}\rangle = \frac{1}{T} \sum_{i=E+1}^{E+T} {\cal O}(\phi_i).
\end{equation}
The statistical error in Monte Carlo behaves as $1/\sqrt{N}$, where $N$ is 
the number of effectively independent configurations.  $N = T/2\tau$, where 
$\tau$ is the autocorrelation time.  This autocorrelation time can easily be 
large, and most of the computer time spent in generating statistically 
independent configurations.  The operator measurements then become a small 
overhead on the whole calculation.

\subsection{QCD}

Quantum Chromo-dynamics (QCD) is the theory describing the strong interaction 
in high energy physics. This interaction binds quarks together to form 
hadrons---the constituents of nuclear matter---for example, the proton and 
neutron.  Currently, we know of five types (referred to as ``flavors'') of 
quark: up, down, strange, charm and bottom; and expect at least one more 
(top) to show up soon.  In addition to having a ``flavor,'' quarks can carry 
one of three possible charges known as ``color'' (this has nothing to do with 
color in the macroscopic world!); hence, Quantum {\em Chromo\/}-dynamics.
The strong color force is mediated by particles called gluons, just as 
photons mediate light in electromagnetism.  Unlike photons, though, gluons 
themselves carry a color charge and, therefore, interact with one another.
This makes QCD a nonlinear theory which is impossible to solve
analytically.  Therefore, we turn to the computer for numerical
solutions.

\subsection{Lattice QCD}

To put QCD on a computer, we proceed as follows\@:  Four-dimensional 
space-time continuum is replaced by a four-dimensional hypercubic periodic 
lattice, of size $N = N_s\!\times\! N_s\!\times\! N_s\!\times\! N_t$ with the 
quarks living on the sites and the gluons living on the links of the lattice.
$N_s$ is the spatial and $N_t$ is the temporal extent of the lattice.  The 
lattice has a finite spacing $a$.  The gluons are represented by 
$3\!\times\! 3$ complex $SU(3)$ matrices associated with each link in the 
lattice.  The $3$ in $SU(3)$ reflects the fact that there are three colors of 
quarks, and $SU$ means that the matrices are unitary with unit determinant 
({\em i.e.,} ``special unitary'').  This link matrix describes how the color 
of a quark changes as it moves from one site to the next.  For example, as a 
quark is transported along a link of the lattice it can change its color 
from, say, red to green; hence, a red quark at one end of the link can 
exchange colors with a green quark at the other end.  The action functional 
for the purely gluonic part of QCD is
\begin{equation}
S_G = \beta \sum_P\, ( 1 - \frac{1}{3} Re Tr U_P ),
\end{equation}
where
\begin{equation}
U_P = U_{i,\mu}\, U_{i+\mu,\nu}\, U_{i+\nu,\mu}^\dagger\, 
U_{i,\nu}^\dagger
\end{equation}
is the product of link matrices around an elementary square or plaquette on 
the lattice --- see Figure~\ref{fig1.qcd}.  Essentially, all of the time in 
QCD simulations of gluons is spent multiplying these $SU(3)$ matrices 
together.  The main component of this is the $a \times b + c$ kernel (which 
most supercomputers can do very efficiently).  As the action involves 
interactions around plaquettes, in order to satisfy detailed balance, we can 
update only half the links in any one dimension simultaneously, as shown in
Figure~\ref{fig2.qcd} (in two dimensions for simplicity).  The partition 
function for full lattice QCD including quarks is
\begin{equation}
Z = \int D\psi D\bar\psi DU \exp(-S_G - \bar\psi (\not\!\!D +m) \psi),
\end{equation}
where $\not\!\!D +m$ is a large sparse matrix the size of the lattice squared.
Unfortunately, since the quark or fermion variables $\psi$ are anticommuting 
Grassmann numbers, there is no simple representation for them on the computer.
Instead, they must be integrated out, leaving a highly non-local fermion 
determinant:
\begin{equation}
Z = \int DU \det(\not\!\!D +m) \exp(-S_G).
\end{equation}
This is the basic integral one wants to evaluate numerically.

\begin{figure}
\vspace*{1.0in}
\caption{A Lattice Plaquette}
\label{fig1.qcd}
\end{figure}

\begin{figure}
\vspace*{1.0in}
\caption{Updating the Lattice}
\label{fig2.qcd}
\end{figure}

The biggest stumbling block preventing large QCD simulations, including 
quarks, is the presence of the determinant $\det(\not\!\!D +m)$ in the 
partition function.  There have been many proposals for dealing with the 
determinant.  The first algorithms tried to compute the change in the 
determinant when a single link variable was updated \cite{Weingarten:81a}.  
This turned out to be prohibitively expensive.  So instead, the approximate 
method of pseudo-fermions \cite{Fucito:81a} was used.  
Today, however, the preferred approach is the so-called ``Hybrid Monte 
Carlo'' algorithm \cite{Duane:87a} which is exact.  
The basic idea is to invent some dynamics for the variables in the system in 
order to evolve the whole system forward in (simulation) time, and then do a 
Metropolis accept/reject for the entire evolution on the basis of the total 
energy change.  The great advantage is that the whole system is updated in 
one fell swoop.  The disadvantage is that if the dynamics are not correct, 
the acceptance will be very small.  Fortunately, (and this is one of very few 
fortuitous happenings where fermions are concerned) good dynamics can be 
found\@: the Hybrid algorithm \cite{Duane:85a}.  This is a neat combination 
of the deterministic microcanonical method \cite{Callaway:83a}, 
\cite{Polonyi:83a} and the stochastic Langevin method \cite{Parisi:81a}, 
\cite{Batrouni:85a}, which yields a quickly evolving, ergodic algorithm for 
both gauge fields and fermions.  The computational kernel of this algorithm 
is the repeated solution of systems of equations of the form,
\begin{equation}
(\not\!\!D +m) \phi = \eta,
\end{equation}
where $\phi$ and $\eta$ are vectors that live on the sites of the lattice.
To solve these equations, one typically uses a conjugate gradient or one of 
its cousins, since the fermion matrix $(\not\!\!D +m)$ is sparse.  For more 
details, see \cite{Gupta:88a}.  Such iterative matrix algorithms have as 
their basic component the $a \times b + c$ kernel, so again computers which 
do this efficiently will run QCD well.

\subsection{Concurrent QCD Machines}

Where is the computer power needed for QCD going to come from?  Today, the 
biggest resources of computer time for research are the conventional 
supercomputers at the NSF and DOE centers.  The centers are continually 
expanding their support for lattice gauge theory, but it may not be long 
before they are overtaken by several dedicated efforts involving concurrent 
computers.  It is a revealing fact that the development of most 
high-performance parallel computers---the Caltech Cosmic Cube, the Columbia 
Machine, IBM's GF11, APE in Rome, the Fermilab Machine---was actually 
motivated by the desire to simulate lattice QCD.  Geoffrey~Fox and
Chuck~Seitz, at Caltech, built the first hypercube computer, the Cosmic Cube 
or Mark~I, in 1983.  It had 64 nodes, each of which was an 
Intel 8086/87 microprocessor with $\rm 128\,KB$ of memory, giving a total of 
about $\rm 2\,Mflops$ (measured for QCD).  This was quickly upgraded to the
Mark~II hypercube with faster chips, twice the memory per 
node, and twice the number of nodes in 1984.  Now QCD is running on the 
latest Caltech hypercube, the 128-node Mark~IIIfp\index{Mark~IIIfp} (built by 
JPL), at $\rm 600\,Mflops$ sustained \cite{Ding:90b}.  Each node of the 
Mark~IIIfp hypercube contains two Motorola 68020 microprocessors, one for 
communication and the other for calculation, with the latter supplemented by 
one 68881 coprocessor and a 32-bit Weitek floating point processor.  
Norman~Christ and Anthony~Terrano, at Columbia, built their first parallel 
computer for doing lattice QCD calculations in 1984.  It had 16 nodes, each 
of which was an Intel 80286/87 microprocessor, plus a TRW 22-bit floating 
point processor with $\rm 1\,MB$ of memory, giving a total peak performance 
of $\rm 256\,Mflops$.  This was improved in 1987 using Weitek rather than TRW 
chips so that 64 nodes give $\rm 1\,Gflops$ peak.  Recently, at Columbia, 
they have finished building their third machine\@: a 256-node $\rm 
16\,Gflops$ lattice QCD computer.  Don~Weingarton initiated the GF11
project in 1984 at IBM.  The GF11 is a SIMD machine comprising 576 Weitek 
floating point processors, each performing at $\rm 20\,Mflops$ to give the 
total $\rm 11\,Gflops$ peak implied by the name.  First preliminary results 
from this project were report at Lattice `91.  The APE (Array Processor with 
Emulator) computer is basically a collection of 3081/E processors (which were 
developed by CERN and SLAC for use in high energy experimental physics) with 
Weitek floating point processors attached.  However, these floating point 
processors are attached in a special way---each node has four multipliers and 
four adders, in order to optimize the $a \times b + c$ calculations, which 
form the major component of all lattice QCD programs.  This means that each 
node has a peak performance of $\rm 64\,Mflops$.  The first small 
machine---Apetto---was completed in 1986 and had four nodes yielding a peak 
performance of $\rm 256\,Mflops$.  Currently, they have a second generation 
of this machine with $\rm 1\,Gflops$ peak from 16 nodes.  By 1992, the APE 
collaboration hopes to have completed the $\rm 100\,Gflops$ 4096-node 
``Apecento'', based on specialized VLSI chips that are software compatible 
with the original APE\@.  Not to be outdone, Fermilab is also using its 
high-energy experimental physics emulators in constructing a lattice QCD 
machine for 1991 with 256 of them arranged as a $2^5$ hypercube of crates, 
with eight nodes communicating through a crossbar in each crate.  All 
together, they expect to get $\rm 5\,Gflops$ peak performance.  These 
performance figures are summarized in Table~\ref{tab2.qcd}.  (The ``real'' 
performances are the actual performances obtained on QCD codes.)

\begin{table}
\begin{center}
\begin{tabular}{|l|l|r|r|}
\hline
\multicolumn{1}{|c|}{Computer}&\multicolumn{1}{c|}{Year}&
\multicolumn{1}{c|}{Peak}&\multicolumn{1}{c|}{Real}
\\
\hline
\medStrut
Caltech I    & 1983 & 3 & 2" \\
Caltech II   & 1984 & 9 & 6" \\
Caltech III  & 1989 & 2000 & 600" \\
Columbia I   & 1984 & 256 & 20" \\
Columbia II  & 1987 & 1000 & 200" \\
Columbia III & 1990 & 16000 & 6000" \\
IBM GF11     & 1990 & 11000 & 10000* \\
APE I        & 1986 & 256 & 20" \\
APE II       & 1988 & 1000 & 200" \\
APE III      & 1992 & 100000 & 20000* \\
Fermilab     & 1991 & 5000 & 1200* \\
\hline
\end{tabular}
\end{center}
*All real times are measured except these predicted ones.
\caption{Peak and Real Performances in Mflops of ``Homebrew'' QCD
Machines}
\label{tab2.qcd}
\end{table}

Major calculations have also been performed on commercial SIMD machines, 
first on the ICL Distributed Array Processor (DAP) at Edinburgh University 
during the period 1982--1987, and now on the TMC Connection Machine (CM-2); 
and on commercial distributed memory MIMD machines like the nCUBE hypercube 
at Caltech.  In fact, currently, the Connection Machine is the most powerful 
commercial QCD machine available\@:  running full QCD at a sustained rate of 
approximately $\rm 5\,Gflops$ and $\rm >5\,Gflops$ after optimization, on a 
$\rm 64\,K$ CM-2 \cite{Baillie:89e}, \cite{Brickner:91b}.

It is interesting to note that when the various groups began building their
``homebrew'' QCD machines, it was clear they would out-perform all
commercial (traditional) supercomputers; however, now that commercial 
parallel supercomputers have come of age \cite{Fox:89n} the situation is not 
so obvious.  To emphasize this, we describe QCD calculations on both the 
home-grown Caltech hypercube and on the commercially available Connection 
Machine.

\subsection{QCD on the Caltech Hypercubes}

To make good use of MIMD distributed memory machines like hypercubes, one 
should employ domain decomposition.  That is, the domain of the problem 
should be divided into subdomains of equal size, one for each processor in 
the hypercube; and communication routines should be written to take care of 
data transfer across the processor boundaries.  Thus, for a lattice 
calculation, the $N$ sites are distributed among the $P = 2^d$ processors 
using a decomposition, which ensures that processors assigned to adjacent 
subdomains are directly linked by a communication channel in the hypercube 
topology.  Each processor then independently works through its subdomain of 
$n = N/P$ sites updating each one in turn, only communicating with 
neighboring processors when doing boundary sites.  This communication 
enforces ``loose synchronization,'' which stops any one processor from racing 
ahead of the others.  Load balancing is achieved with equal-size domains.  If 
the nodes contain at least two sites of the lattice, all the nodes can 
update in parallel, satisfying detail balance, since loose synchronicity 
guarantees that all nodes will be doing black, then red sites alternately.
The characteristic timescale of the communication, $t_{comm}$, corresponds to 
roughly the time taken to transfer a single $SU(3)$ matrix from one node to 
its neighbor.  Similarly, we can characterize the calculational part of the 
algorithm by a timescale, $t_{calc}$, which is roughly the time taken to 
multiply together two $SU(3)$ matrices.  For all hypercubes built {\em 
without} floating point accelerator chips $t_{comm} << t_{calc}$ and, hence, 
QCD simulations are extremely ``efficient,'' where efficiency is defined by 
the relation
\begin{equation}
\epsilon  =  \frac{T_1}{k\, T_k},
\end{equation}
where $T_k$ is the time taken for $k$ processors to perform the given 
calculation.  Typically, such calculations have efficiencies in the range 
$\epsilon \geq .90$, which means they are ideally suited to this type of 
computation since doubling the number of processors nearly halves the 
total computational time required for solution.  However, as we shall see 
(for the Mark~IIIfp hypercube, for example), the picture changes dramatically 
when fast floating point chips are used; then $t_{comm} \simeq t_{calc}$ and 
one must take some care in coding to obtain maximum performance.

Rather than describe every calculation done on the Caltech hypercubes,
we shall concentrate on one calculation that has been done several times
as the machine evolved---the heavy quark potential\index{quark potential} 
calculation (``heavy'' because the quenched approximation is used).

QCD provides an explanation of why quarks are confined inside hadrons in that
lattice calculations reveal that the inter-quark potential rises linearly as
the separation between the quarks increases.  Thus, the attractive force
between quarks grows linearly with the separation.  This force, called the 
``string tension'', is carried by the gluons, which form a ``string'' between 
the quarks.  On the other hand, at short distances, quarks and gluons are 
``asymptotically free'' and behave like electrons and photons, interacting 
via a Coulomb-like force.  Thus, the quark potential $V$ is written as
\begin{equation}
V(R) = -\frac{\alpha}{R} + \sigma R,
\end{equation}
where $R$ is the separation of the quarks, $\alpha$ is the coefficient of the 
Coulombic potential and $\sigma$ is the string tension.  Experimentally, in 
fitting charmonium data to this Coulomb plus linear potential, Eichten
{\em et al.} \cite{Eichten:80a} estimated that $\alpha=0.52$ and 
$\sigma$=0.18GeV$^2$.  Thus, a goal of the lattice calculations is to 
reproduce these numbers.  Enroute to this goal, it is necessary to show that 
the numbers from the lattice are ``scaling'', that is, if one calculates a 
physical observable on lattices with different spacings then one gets the 
same answer.  This means that the artifacts due to the finiteness of the 
lattice spacing have disappeared and continuum physics can be extracted.

The first heavy quark potential calculation using a Caltech hypercube was 
performed on the 64-node Mark~I in 1984 on a $12^3\times 16$ lattice with 
$\beta$ ranging from $5.8$ to $7.6$ \cite{Otto:84a}.  $\alpha$ was found to 
be $0.25$ and the string tension (converting to the dimensionless ratio) 
$\sqrt{\sigma}/\Lambda = 106$.  The numbers are quite a bit off from the 
charmonium data but the string tension did appear to be scaling, albeit in 
the narrow window $6.0 \le \beta \le 6.4$.

The next time around, in 1986, the 128-node Mark~II hypercube was used on a 
$20^4$ lattice with $\beta = 6.1, 6.3, 6.5, 6.7$ \cite{Flower:86b}.  The 
dimensionless string tension decreased somewhat to 83, but clear violations 
of scaling were observed\@: the lattice was still too coarse to see continuum 
physics.

Therefore, the latest (1989) calculation using the Caltech/JPL 32-node 
Mark~IIIfp hypercube concentrated on one $\beta$ value, $6.0$, and 
investigated different lattice sizes\@:  $24^3\times 10$, $24^3\times 12$, 
$12^3\times 24$, $16^3\times 24$ \cite{Ding:90b}.  Scaling was not 
investigated, however, the values of $\alpha=0.58$ and $\sqrt{\sigma}/\Lambda 
= 77$, {\em i.e.,\/} $\rm \sigma$=0.15\,GeV$^2$ compare favorably with the 
charmonium data.  This work is based on about 1300~CPU hours on the 32-node
Mark~IIIfp hypercube, which has a performance of roughly twice a CRAY~X-MP 
processor.  The whole 128-node machine performs at $\rm 0.6\,Gflops$.  As 
each node runs at $\rm 6\,Mflops$, this corresponds to a speedup of 100, and 
hence, an efficiency of 78\%.  These figures are for the most highly 
optimized code.  The original version of the code written in C ran on the 
Motorola chips at $\rm 0.085\,Mflops$ and on the Weitek chips at $\rm 
1.4\,Mflops$.  The communication time, which is roughly the same for both, is 
less than a 2\% overhead for the former but nearly 30\% for the latter.  When 
the computationally intensive parts of the calculation are written in 
assembly code for the Weitek, this overhead becomes almost 50\%.  This $\rm 
0.9\,msec$ of communication, shown in lines two and three in 
Table~\ref{tab3.qcd}, is dominated by the hardware/software message startup 
overhead (latency), because for the Mark~IIIfp the node to node communication 
time, $t_{comm}$, is given by
\begin{equation}
t_{comm} \simeq ( 150+2*W) \quad \mu\, {\rm sec},
\end{equation}
where $W$ is the number of words transmitted.  To speed up the communication, 
we update all even (or odd) links (eight in our case) in each node, allowing
us to transfer eight matrix products at a time, instead of just sending one 
in each message.  This reduces the $\rm 0.9\,msec$ by a factor  of
\begin{equation}
\frac{8*(150+18*2)}{ 150+8*18*2} = 3.4
\end{equation}
to $\rm 0.26\,msec$.  On all hypercubes with fast floating point chips---and 
on most hypercubes without these chips for less computationally intensive 
codes---such vectorization of communication is often important.  In 
Figure~\ref{fig3.qcd}, the speedups for many different total lattice sizes 
are shown.  For the largest lattice size, the speedup is 100 on the 128-node. 
The speedup is almost linear in number of nodes.  As the total lattice volume 
increases, the speedup increases, because the ratio of 
calculation/communication increases.  For more information on this 
performance analysis, see \cite{Ding:90c}.

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{Programming Level}&\multicolumn{1}{c|}{Calc.}&
\multicolumn{1}{c|}{Comm.}&\multicolumn{1}{c|}{Total}&
\multicolumn{1}{c|}{Mflops} \\
&\multicolumn{1}{c|}{Time}&\multicolumn{1}{c|}{Time}&
\multicolumn{1}{c|}{Time}& \\
\hline
\medStrut
1.\ \ Motorola MC68020/68881 in C & 52\phantom{.}"" & 0.86 & 53\phantom{.}"" &
0.085 \\
2.\ \ Weitek XL all in C             & "2.25 & 0.90 & "3.15 & 1.4"" \\
3.\ \ Weitek XL parts in Assembly    & "0.94 & 0.90 & "1.84 & 2.4"" \\
4.\ \ Weitek XL Assembly, vec. comm. & "0.94 & 0.26 & "1.20 & 3.8"" \\
5.\ \ Weitek XL Assembly, no comm.   & "0.94 & 0.0" & "0.94 & 4.8"" \\
\hline
\end{tabular}
\caption{Link Update Time (msec) on Mark~IIIfp Node for Various Levels
of Programming}
\label{tab3.qcd}
\end{center}
\end{table}

\begin{figure}
\vspace*{2.0in}
\caption{Speedups for QCD on the Mark~IIIfp}
\label{fig3.qcd}
\end{figure}

\subsection{QCD on the Connection Machine}

The Connection Machine Model CM-2 is also very 
well-suited for large-scale simulations of QCD.  The CM-2 is a 
distributed-memory, single-instruction multiple-data (SIMD) massively 
parallel processor comprising up to 65536 ($\rm 64\,K$) processors 
\cite{Hillis:85a} \cite{Hillis:87a}.  Each processor consists of an 
arithmetic-logic unit (ALU), eight or $\rm 32\,Kbytes$ of random-access 
memory (RAM) and a router interface to perform communications among the 
processors.  There are 16 processors and a router per custom VLSI chip, with 
the chips interconnected as a 12-dimensional hypercube.  Communications
among processors within a chip work essentially like a crossbar interconnect.
The router can do general communications, but we require only local ones for 
QCD so we use the fast nearest-neighbor communication software called NEWS\@.
The processors deal with one bit at a time.  Therefore, the ALU can compute 
any two Boolean functions as output from three inputs, and all data paths are
one-bit wide.  In the current version of the Connection Machine (the CM-2),
groups of 32 processors (two chips) share a 32-bit (or 64-bit) Weitek 
floating point chip, and a transposer chip, which changes 32 bits stored 
bit-serially within 32 processors into 32 32-bit words for the Weitek, and 
vice versa.

The high-level languages on the CM, such as *Lisp and CM-Fortran, compile into
an assembly language called Parallel Instruction Set (Paris).  Paris regards 
the $\rm 64\,K$ bit-serial processors as the fundamental units in the 
machine.  However, floating-point computations are not very efficient in the 
Paris model.  This is because in Paris, 32-bit floating point numbers are 
stored ``field-wise''; that is, successive bits of the word are stored at 
successive memory locations of each processor's memory.  However, 32 
processors share one Weitek chip, which deals with words stored 
``slice-wise''---stored across the processors, one bit in each.  Therefore, 
to do a floating-point operation, Paris loads in the field-wise operands, 
transposes them slice-wise for the Weitek (using the transposer chip), does 
the operation, and transposes the slice-wise result back to field-wise for 
memory storage.  Moreover, every operation in Paris is an atomic process; 
that is, two operands are brought from memory and one result is stored back 
to memory, so no use is made of the Weitek registers for intermediate 
results.  Hence, to improve the performance of the Weiteks, a new assembly 
language called CM Instruction Set (CMIS) has been written, which models the 
local architectural features much better.  In fact, CMIS ignores the 
bit-serial processors and thinks of the machine in terms of the Weitek chips. 
Thus, data can be stored slice-wise, eliminating all the transposing back and 
forth.  CMIS allows effective use of the Weitek registers, creating a memory 
hierarchy, which combined with the internal buses of the Weiteks, offers 
increased bandwidth for data motion.

When the arithmetic part of the program is rewritten in CMIS (just like
on the Mark~IIIfp when this was rewritten in assembly code) the
communications become a bottleneck.  Therefore, we need also to speed
up the communication part of the code.  On the CM-2, this is done using
the ``bi-directional multi-wire NEWS'' system.  As explained above, the
CM chips (each containing 16 processors) are interconnected in a
12-dimensional hypercube.  However, since there are two CM chips for
each Weitek floating-point chip, the floating-point hardware is
effectively wired together as an 11-dimensional hypercube, with two
wires in each direction.  This makes it feasible to do simultaneous
communications in both directions of all four space-time directions in
QCD---bi-directional multi-wire NEWS---thereby reducing the
communication time by a factor of eight.  Moreover, the data
rearrangement necessary to make use of this multi-wire NEWS further
speeds up the CMIS part of the code by a factor of two.

Currently, the Connection Machine is the most powerful commercial QCD machine 
available\@: the ``Los Alamos collaboration'' is running full QCD at a 
sustained rate of $\rm 5.3\,Gflops$ on a $\rm 64\,K$ CM-2.  As was the 
case for the Mark~IIIfp hypercube, in order to obtain this performance, one 
must resort to writing assembly code for the Weitek chips {\em and\/} for the 
communication.  Our original code, written entirely in the CM-2 version
of *Lisp, achieved around $\rm 1\,Gflop$ \cite{Baillie:89e}.  As shown in 
Table~\ref{tab4.qcd}, this code spends 34\% doing communication.  When we 
rewrote the most computationally intensive part in the assembly language 
CMIS, this rose to 54\%.  Then when we also made use of ``multi-wire
NEWS'' (to reduce the communication time by a factor of eight) it fell
to 30\%.

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|r|}
\hline
\multicolumn{1}{|c|}{Programming Level}&\multicolumn{1}{c|}{Calc. time}&
\multicolumn{1}{c|}{Comm. time}&\multicolumn{1}{c|}{Total Time}&
\multicolumn{1}{c|}{Mflops}
\\
\hline
\medStrut
All in *Lisp & 8.7 & 4.5 & 13.2 & 900 \\
Inner Loop in CMIS & 3.3 & 3.9 & "7.2 & 1600 \\
Multi-wire CMIS & 1.7 & 0.5 & "2.2 & 5300
\\
\hline
\end{tabular}
\caption{Fermion Update Time (sec) on $\rm 64\,K$ Connection Machine for
Various Levels of Programming}
\label{tab4.qcd}
\end{center}
\end{table}

\subsection{Status and Prospects}

The status of lattice QCD may be summed up as\@: underway.  Already there 
have been some nice results in the context of the quenched approximation, but 
the lattices are still too coarse and too small to give definitive results.
Results for full QCD are going to take orders of magnitude more computer 
time, but we now have an algorithm---Hybrid Monte Carlo---which puts real 
simulations within reach.

When will the computer power be sufficient?  In Figure~\ref{fig4.qcd}, we 
plot the horsepower of various QCD machines as a function of the year they 
started to produce physics results.  The performance plotted in this case is 
the real sustained rate on actual QCD codes.  The surprising fact is that the 
rate of increase is very close to exponential, yielding a factor of ten every 
two years!  On the same plot, we show our estimate of the computer power 
needed to re-do correct quenched calculations on a $128^4$ lattice.  This 
estimate is also a function of time, due to algorithm improvements.

\begin{figure}
\vspace*{2.0in}
\caption{Megaflops for QCD Calculations}
\label{fig4.qcd}
\end{figure}

Extrapolating both lines, we see the outlook for lattice QCD is rather bright.
Reasonable results for the phenomenologically interesting physical 
observables should be available within the quenched approximation in the 
mid-90's.  With the same computer power, we will be able to re-do today's 
quenched calculations using dynamic fermions (but still on today's size of 
lattice).  This will tell us how reliable the quenched approximation is.  
Finally, results for the full theory with dynamical fermions on a $128^4$ 
lattice should follow early in the next century (!), when computers are two 
or three orders of magnitude more powerful.

\section{Spin Models}
\label{spins}

\subsection{Introduction}

Spin models\index{Spin models} are simple statistical models of real systems, 
which exhibit the same behavior and, hence, provide an understanding of the 
physical mechanisms involved.  Despite their apparent simplicity, most of 
these models are not exactly soluble by present theoretical methods.  Hence, 
computer simulation is used.  Usually, one is interested in the behavior of 
the system at a phase transition; the computer simulation reveals where the 
phase boundaries are, what the phases on either side are, and how the 
properties of the system change across the phase transition.  There are two 
varieties of spins\@: discrete or continuously valued.  In both cases, the 
spin variables are put on the sites of the lattice and only interact with 
their nearest neighbors.  The partition function for a spin model is
\begin{equation}
Z = \int Ds \exp(S),
\end{equation}
with the action being of the form
\begin{equation}
S = \beta \sum_{\langle ij\rangle} ( 1 - s_i.s_j ),
\end{equation}
where $\langle ij\rangle$ denotes nearest neighbors.  A great deal of work 
has been done over the years in finding good algorithms for computer 
simulations of spin models; recently some new, much better, algorithms have 
been discovered.  These so-called cluster algorithms are described in
detail in Section~\ref{cluster}.  Here, we shall use them to perform
large-scale Monte Carlo simulations of several spin models---both
discrete and continuous.

\subsection{Ising Model}

The Ising model\index{Ising model} is the simplest model for ferromagnetism 
that predicts phase transitions and critical phenomena.  This model, 
introduced by Lenz in 1920 \cite{Lenz:20a}, was solved in one dimension by 
Ising in 1925 \cite{Ising:25a}, and in two dimensions by Onsager in 1944 
\cite{Onsager:44a}.  However, it has not been solved analytically in three 
dimensions, so Monte Carlo computer simulation methods have been one of the 
methods used to obtain numerical solutions.  One of the best available 
techniques for this is the Monte Carlo Renormalization Group\index{Monte
Carlo Renormalization Group} (MCRG) method \cite{Wilson:80a}, 
\cite{Swendsen:79a}.  The Ising model exhibits a second-order phase 
transition in $d=3$ dimensions at a critical temperature $T_c$.  As $T$ 
approaches $T_c$, the correlation length $\xi$ diverges as a power law with 
critical exponent $\nu$\@:  
\begin{equation}
\xi = \xi_0(T/T_c-1)^{-\nu}
\end{equation}
and the pair correlation function $G(r)$ at $T=T_c$ falls off to zero with
distance $r$ as a power law defining the critical exponent\index{critical 
exponents} $\eta$\@:
\begin{equation}
G(r) = G_0 r^{-(d-2+\eta)}.
\end{equation}
$T_c$, $\nu$ and $\eta$ determine the critical behavior of the 3-D Ising model
and it is their values we wish to determine using MCRG.

In 1984, this was done by Pawley, Swendsen, Wallace and Wilson
\cite{Pawley:84a} in Edinburgh on the ICL DAP computer with high statistics.
They ran on four lattice sizes---$8^3$, $16^3$, $32^3$ and $64^3$---measuring 
seven even and six odd spin operators.  We are essentially repeating their 
calculation on the new AMT DAP computer.  Why should we do this?  First, to 
investigate finite size effects---we have run on the biggest lattice used by 
Edinburgh, $64^3$, and on a bigger one, $128^3$\@.  Second, to investigate 
truncation effects---qualitatively the more operators we measure for MCRG, 
the better, so we have included 53 even and 46 odd operators.  Third, we are 
making use of the new cluster-updating algorithm due to Swendsen and Wang 
\cite{Swendsen:87a}, implemented according to Wolff \cite{Wolff:89b}.  
Fourth, we would like to try to measure another critical exponent more
accurately---the correction-to-scaling exponent $\omega$, which plays an
important role in the analysis.

In order to calculate the quantities of interest using MCRG, one must 
evaluate the spin operators $S_\alpha$.  In \cite{Pawley:84a}, the 
calculation was restricted to seven even spin operators and six odd; we are 
evaluating 53 and 46, respectively.  Specifically, we have decided to 
evaluate the most important operators in a $3\times 3\times 3$ cube 
\cite{Baillie:88h}.  To determine the critical coupling (or inverse 
temperature), $K_c$, one performs independent MC simulations on a large 
lattice $L$ of size $M^3$ and on smaller lattices $S$ of size $(M/2^m)^3$, 
$m = 1,2,\,\ldots$, and compares the operators measured on the large lattice 
blocked $m$ times more than the smaller lattices.  $K = K_c$ when they are 
the same.  Since the effective lattice sizes are the same, unknown finite 
size effects should cancel.  The critical exponents, $y_a$, are obtained 
directly from the eigenvalues, $\lambda_a$, of the stability matrix, 
$T_{\alpha\beta}$, according to $\lambda_a = 2^{y_a}$.  In particular, the 
leading eigenvalue $y_t$ of $T_{\alpha\beta}$ for the even $S_\alpha$ gives 
$\nu$ from $y_5 = 1/\nu$, and similarly $y_h = (d+2-\eta)/2$ from the odd 
eigenvalue $y_h$ of $T_{\alpha\beta}$.

The Distributed Array Processor (DAP) is a SIMD computer consisting of $D 
\times D$ bit-serial processing elements (PEs) configured as a cyclic 
two-dimensional grid with nearest-neighbor connectivity.  The Ising model 
computer simulation is well-suited to such a machine since the spins can be 
represented as single-bit (logical) variables.  In 3-D, the system of spins 
is configured as an $M \times M \times M$ simple cubic lattice, which is 
``crinkle mapped'' onto the $D \times D$ DAP by storing $N \times N$ pieces 
of each of $M$ planes in each PE\@: $M \times M \times M = M \times (N\times 
D) \times (N\times D)$, with $N = M/D$.  Our Monte Carlo simulation uses a 
hybrid algorithm in which each sweep consists of 10 standard Metropolis 
\cite{Metropolis:53a} spin-updates followed by one cluster-update using 
Wolff's single-cluster variant of the Swendsen and Wang algorithm.  On the 
$128^3$ lattice, the autocorrelation time of the magnetization reduces from 
$73 \pm 2$ sets of $100$ sweeps for Metropolis alone to $5.0 \pm 0.2$ sets of 
$10$ Metropolis plus one cluster-update for the hybrid algorithm.  In order 
to measure the spin operators, $S_\alpha$, the DAP code simply histograms the 
spin configurations so that an analysis program can later pick out each 
particular spin operator using a look-up table.  Currently, the code requires 
the same time to do one histogram measurement, one Wolff 
single-cluster-update or $100$ Metropolis updates.  Therefore, our hybrid of 
$10$ Metropolis plus one cluster-update takes about the same time as a 
measurement.  On a DAP~510, this hybrid update takes on average 127
secs (13.5 secs) for the $128^3$ ($64^3$) lattices, so that the total
time used for this calculation is roughly 11,000 hours.  We have performed
simulations on $64^3$ and $128^3$ lattices at two values of the
coupling\@:  $K_{nn} = 0.221654$ (Edinburgh's best estimate) and 
$K_{nn} = 0.221644$.  We accumulated $\rm 140\,K$ measurements for each
of the $128^3$ simulations and $\rm 100\,K$ for the $64^3$.  For error
analysis, this is divided into bins of $\rm 10\,K$ measurements.

In analyzing our results, the first thing we have to decide is the order in 
which to arrange our 53 even and 46 odd spin operators.  Naively, they can be 
arranged in order of increasing total distance between the spins 
\cite{Baillie:88h} (as was done in \cite{Pawley:84a}).  However, the ranking 
of a spin operator is determined physically by how much it contributes to 
the energy of the system.  Thus, we did our analysis initially with the 
operators in the naive order to calculate their energies, then subsequently 
we used the ``physical'' order dictated by these energies.  This physical 
order of the first 20 even operators is shown in Figure~\ref{fig1.spins} with
six of Edinburgh's operators indicated; the seventh Edinburgh operator (E-6) 
is our 21st.  This order is important in assessing the systematic effects of 
truncation, as we are going to analyze our data as a function of the 
number of operators included.  Specifically, we successively diagonalize the 
$1\times 1$, $2\times 2$, $\ldots$, (for even, $46\times 46$ for odd) 
stability matrix $T_{\alpha\beta}$ to obtain its eigenvalues and, thus, the 
critical exponents.

\begin{figure}
\vspace*{2.5in}
\caption{Our Order for Even Spin Operators}
\label{fig1.spins}
\end{figure}

We present our results in terms of the eigenvalues of the even and odd parts 
of $T_{\alpha\beta}$.  The leading even eigenvalue on the first four blocking
levels starting from the $64^3$ lattice is plotted against the number of 
operators included in the analysis in Figure~\ref{fig2.spins}, and on the 
first five blocking levels starting from the $128^3$ lattice in 
Figure~\ref{fig3.spins}.  Similarly, the leading odd eigenvalues for $64^3$ 
and $128^3$ lattices are shown in Figure~\ref{fig4.spins} and
Figure~\ref{fig5.spins} respectively.  First of all, note that there are 
significant truncation effects---the value of the eigenvalues do not settle 
down until at least thirty and perhaps forty operators are included.  We note 
also that our value agrees with Edinburgh's when around seven operators are
included---this is a significant verification that the two calculations are
consistent.  With most or all of the operators included, our values on the 
two different lattice sizes agree, and the agreement improves with increasing 
blocking level.  Thus, we feel that we have overcome the finite size effects 
so that a $64^3$ lattice is just large enough.  However, the advantage in 
going to $128^3$ is obvious in Figure~\ref{fig3.spins} and
Figure~\ref{fig5.spins}\@:  there, we can perform one more blocking, which 
reveals that the results on the fourth and fifth blocking levels are 
consistent.  This means that we have eliminated most of the transient effects 
near the fixed point in the MCRG procedure.  We also see that the main
limitation of our calculation is statistics--the error bars are still
rather large for the highest blocking level.

Now in order to obtain values for $\nu$ and $\eta$, we must extrapolate
our results from a finite number of blocking levels to an infinite
number.  This is done by fitting the corresponding eigenvalues
$\lambda_t$ and $\lambda_h$ according to
\begin{equation}
\lambda(n) = \lambda^{\ast} + a2^{-\omega n}
\end{equation}
where $\lambda^{\ast}$ is the extrapolated value and $\omega$ is the
correction-to-scaling exponent.  Therefore, we first need to calculate
$\omega$, which comes directly from the second leading even
eigenvalue\@: $\omega = -y_{t,2}$.  Our best estimate is in the interval
$\omega = 0.80$--0.85, and we use the value 0.85 for the purpose of
extrapolation, since it gives the best fits.  The final results are
$\nu = 0.624\pm 0.001\pm 0.002$, $\eta = 0.0262\pm 0.0002\pm 0.0030$,
where the first errors are statistical and the second errors are
estimates of the systematic error coming from the uncertainty in
$\omega$.

\begin{figure}
\vspace*{2.0in}
\caption{Leading Even Eigenvalue on $64^3$ Lattice}
\label{fig2.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Leading Even Eigenvalue on $128^3$ Lattice}
\label{fig3.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Leading Odd Eigenvalue on $64^3$ Lattice}
\label{fig4.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Leading Odd Eigenvalue on $128^3$ Lattice}
\label{fig5.spins}
\end{figure}

Finally, perhaps the most important number, because it can be determined the 
most accurately, is $K_{nn}^c$.  By comparing the fifth blocking level on the 
$128^3$ lattice to the fourth on the $64^3$ lattice for both coupling
values and taking a weighted mean, we obtain $K_{nn}^c = 0.221652\pm
0.000003\pm 0.000001$, where again the first error is statistical and
the second is systematic.

Thus, MCRG calculations give us very accurate values for the three critical 
parameters $K_{nn}^c$, $\eta$, and $\nu$, and give a reasonable estimate
for $\omega$.  Each parameter is obtained independently and directly
from the data.  We have shown that truncation and finite size errors at
all but the highest blocking level have been reduced to below the
statistical errors.  Future high statistics simulations on $256^3$
lattices will significantly reduce the remaining errors and allow us to
determine the exponents very accurately.

\subsection{Potts Model}
\label{potts}

The $q$-state Potts model\index{Potts model} \cite{Potts:52a} consists of a 
lattice of spins $\sigma_{i}$, which can take $q$ different values, and whose 
Hamiltonian is
\begin{equation}
H = K \sum_{\langle i,\,j\rangle} \delta_{\sigma_{i}\sigma_{j}} .
\end{equation}
For $q=2$, this is equivalent to the Ising model.  The Potts model is thus a 
simple extension of the Ising model; however, it has a much richer phase 
structure, which makes it an important testing ground for new theories and
algorithms in the study of critical phenomena \cite{Wu:82a}. 

Monte Carlo simulations of Potts models have traditionally used local 
algorithms such as that of Metropolis, {\em et al.} \cite{Metropolis:53a},
however, these algorithms have the major drawback that near a phase 
transition the autocorrelation time (the number of sweeps needed to generate 
a statistically independent configuration) increases approximately as $L^2$, 
where $L$ is the linear size of the lattice.  New algorithms have recently 
been developed that dramatically reduce this ``critical slowing
down''\index{critical slowing down} by 
updating clusters of spins at a time (these algorithms are described in 
Section~\ref{cluster}).  The original cluster algorithm of Swendsen and Wang 
was implemented for the Potts model \cite{Swendsen:87a}, and there is a lot 
of interest in how well cluster algorithms perform for this model.  At 
present, there are very few theoretical results known about cluster 
algorithms, and theoretical advances are most likely to come from first 
studying the simplest possible models.

We have made a high statistics study of the SW algorithm and the single 
cluster Wolff algorithm \cite{Wolff:89b}, as well as a number of variants of 
these algorithms, for the $q=2$ and $q=3$ Potts models in two dimensions 
\cite{Baillie:89m}.  We measured the autocorrelation time $\tau$ in the 
energy (a local operator) and the magnetization (a global one) on lattice 
sizes from $8^2$ to $512^2$.   About 10 million sweeps were required for each 
lattice size in order to measure autocorrelation times to within about 1\%.
From these values, we can extract the dynamic critical exponent $z$, given 
by $\tau(L) \sim L^{z}$, where $\tau(L)$ is measured at the infinite volume 
critical point (which is known exactly for the $2D$ Potts model). 

The simulations were performed on a number of different parallel computers. 
For lattice sizes of $128^2$ or less, it is possible to run independent
simulations on each processor of a parallel machine, enabling us to obtain
100\% efficiency by running 10 or 20 runs for each lattice size in parallel.
These calculations were done using a 32-processor Meiko Computing Surface,
a 20-processor Sequent Symmetry, a 20-processor Encore Multimax, and
a 96-processor BBN GP1000 Butterfly, as well as a network of SUN workstations,
and took approximately 15,000 processor-hours.  For the largest lattice sizes, 
$256^2$ and $512^2$, a parallel cluster algorithm was required, due to the 
large amount of calculation (and memory) required.  We have used the 
self-labeling algorithm described in Section~\ref{cluster}, which gives 
fairly good efficiencies of about 70\% on the machines we have used (an 
nCUBE-1 and a Symult S2010), by doing multiple runs of 32 nodes each for the 
$256^2$ lattice, and 64 nodes for $512^2$.  Since this problem does not
vectorize, using all 512 nodes of the nCUBE gives a performance approximately 
five times that of a single processor CRAY~X-MP, while all 192 nodes of the 
Symult is equivalent to about six CRAYs.  The calculations on these machines 
have so far taken about 1000 hours.

Results for the autocorrelation times of the energy for the Wolff and SW 
algorithms are shown in Figure~\ref{fig6.spins} for $q=2$ and
Figure~\ref{fig7.spins} for $q=3$.  As can be seen, the Wolff algorithm has 
smaller autocorrelation times than SW.  However, the dynamical critical 
exponents for the two algorithms appear to be identical, being approximately 
$0.25(1)$ and $0.57(1)$ for $q=2$ and ${q=3}$ respectively (shown as
straight lines in Figure~\ref{fig6.spins} and Figure~\ref{fig7.spins}), 
compared to values of approximately two for the standard Metropolis 
algorithm. 

\begin{figure}
\vspace*{2.0in}
\caption{Energy Autocorrelation Times, $q=2$}
\label{fig6.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Energy Autocorrelation Times, $q=3$}
\label{fig7.spins}
\end{figure}

Burkitt and Heermann \cite{Heermann:90a} have suggested that the increase in 
the autocorrelation time is a logarithmic one, rather than a power law for 
the $q=2$ case (the Ising model), {\em i.e.,\/} $z = 0$.  Fits to this
are shown as dotted lines in Figure~\ref{fig6.spins}.  These have
smaller $\chi^2$ values than the power law fits, favoring logarithmic
behavior.  However, it is very difficult to distinguish between a logarithm 
and a small power even on lattices as large as $512^2$.  In any case, the 
performance of the cluster algorithms for the Potts model is quite 
extraordinary, with autocorrelation times for the $256^2$ lattice, which are 
hundreds of times smaller than for the Metropolis algorithm.  In the future, 
we hope to use the cluster algorithms to perform greatly improved Monte Carlo 
simulations of various Potts models, to study their critical behavior.

\subsection{XY Model}

The XY (or O(2)) model\index{XY model} consists of a set of continuous valued 
spins regularly arranged on a two-dimensional square lattice.  15 years ago, 
Kosterlitz and Thouless\index{Kosterlitz and Thouless} (KT) predicted that 
this system would undergo a phase transition\index{transitions} as one 
changed from a low temperature spin wave phase to a high 
temperature phase with unbound vortices.  KT  predicted an approximate 
transition temperature, $T_c$, and the following unusual exponential 
singularity in the correlation length and magnetic susceptibility\@:
\begin{equation}
\xi = a_{\xi} e^{b_{\xi}t^{-\nu}},\quad\quad \chi = 
a_\chi e^{b_\chi t^{-\nu}},
\end{equation}
with
\begin{equation}
\nu = \frac{1}{2},\quad\quad \eta = \frac{1}{4},
\end{equation}
where $t = T-T_c$ and the correlation function exponent $\eta$ is defined by 
the relation $\chi = c\, \xi^{2-\eta}$.  

Our simulation \cite{Gupta:88a} was done on the 128-node FPS (Floating Point 
Systems) T-Series hypercube at Los Alamos.  FPS software allowed the use of C 
with a similar software model (communication implemented by subroutine call) 
to that used on the hypercubes at Caltech.  Each FPS node is built around 
Weitek floating point units, and we achieved $\rm 2\,Mflops$ per node in this 
application.  The total machine ran at $\rm 1/4\,Gflops$, or at about twice 
the performance of one processor of a CRAY~X-MP for this application.  We use 
a $1-D$ torus topology for communications, with each node processing a 
fraction of the rows.  Each row is divided into red/black alternating sites of
spins and the vector loop is over a given color.  This gives a natural data 
structure of ($n \times 128$) words for lattices of size $(n \times 256)$.  
The internode communications, in both lattice update and in measurement of 
observables, can be done asynchronously and are a negligible overhead.

Previous numerical work was unable to confirm the KT theory, due to limited 
statistics and small lattices.  Our high-statistics simulations are done on 
$64^2$, $128^2$, $256^2$, and $512^2$ lattices using a combination of 
over-relaxed and Metropolis algorithms which decorrelates as $\tau \approx 
0.15 \xi^{1.2}$.  (For comparison, a Metropolis algorithm decorrelates as 
$\tau \approx 5 \xi^2$.)  Each configuration represents $N_{or}$ over-relaxed 
sweeps through the lattice followed by $N_{met}$ Metropolis sweeps.  
Measurement of observables is made on every configuration.  The over-relaxed 
algorithm consists of reflecting the spin at a given site about $\Sigma$ 
where $\Sigma$ is the sum of the nearest-neighbor spins, {\em i.e.},
\begin{equation}
s_{new} \ \ = \ \ \Sigma \ s_{old}^\dagger \ \Sigma.
\end{equation}
This implementation \cite{Creutz:87a} \cite{Brown:87a} of the over-relaxed 
algorithm is microcanonical, and it reduces critical slowing down even though 
it is a local algorithm.  The ``hit'' elements for the Metropolis algorithm
are generated as $e^{i\theta}$, where $\theta$ is a uniform random number
in the interval $\lbrack - \delta, \delta\rbrack$, and $\delta$ is adjusted to 
give an acceptance rate of 50--60\%.  The Metropolis hits make the algorithm 
ergodic, but their effectiveness is limited to local changes in the energy.
In Figure~\ref{fig8.spins}, we show the autocorrelation time $\tau$ vs. the 
correlation length $\xi$\@; for $N_{or} = 8$,  $N_{met} = 2$ we extract $\tau 
= 0.15\xi^{1.48}$ and for $N_{or} = 15$, $N_{met} = 2$ we get $\tau = 
0.15\xi^{1.20}$.

\begin{figure}
\vspace*{2.0in}
\caption{Autocorrelation Times for the XY Model}
\label{fig8.spins}
\end{figure}

We ran at 14 temperatures near the phase transition and made unconstrained 
fits to all 14 data points (four parameter fits according to equations~(XY)),
for both the correlation length (Figure~\ref{fig9.spins}) and susceptibility 
(Figure~\ref{fig10.spins}).  The key to the interpretation of the data is the 
fits.  We find that fitting programs (e.g., MINUIT, SLAC) move incredibly
slowly towards the true minimum from certain points (which we label spurious
minima), which, unfortunately, are the attractors for most starting points.
We found three such spurious minima (KT1-3) and the true minimum KT4, as 
listed in Table~\ref{tab1.spins}.

\begin{figure}
\vspace*{2.0in}
\caption{Correlation Length for the XY Model}
\label{fig9.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Susceptibility for the XY Model}
\label{fig10.spins}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|l|}
\hline
&\multicolumn{1}{|c|}{KT1}&\multicolumn{1}{c|}{KT2}&\multicolumn{1}{c|}{KT3}&
\multicolumn{1}{c|}{KT4}
\\
\hline
\medStrut
$a_{\chi}$      & "2.155 & "0.270 & 0.784 & 0.2727(3)  \\
$b_{\chi}$      & "1.088 & "2.809 & 1.844 & 2.7734(5)  \\
$T_c$           & "0.701 & "0.911 & 0.867 & 0.8987(1)  \\
$\nu$           & "1.460 & "0.472 & 0.701 & 0.4995(11) \\
$\chi^2/_{dof}$ & 73.0"" & 13.2"" & 1.82" & 1.70       \\
$a_{\xi}$       & "0.708 & "0.228 & 0.441 & 0.2280(3)  \\
$b_{\xi}$       & "0.740 & "1.693 & 1.100 & 1.6725(8)  \\
$T_c$           & "0.689 & "0.910 & 0.867 & 0.8967(2)  \\
$\nu$           & "1.412 & "0.471 & 0.700 & 0.5007(3)  \\
$\chi^2/_{dof}$ & "7.94" & "2.62" & 0.89" & 0.37       \\
$\eta$          & "0.53" & "0.34" & 0.324 & 0.34       \\
\hline
\end{tabular}
\caption{Results of the XY model fits\@: (a) $\chi$ in T, and  (b) $\xi$ 
in T assuming the KT form.  The fits KT1-3 are pseudo-minima while KT4 is 
the true minimum.  All data points are included in the fits and we give the 
$\chi^2/_{dof}$ for each fit and an estimate of the exponent $\eta$.}
\label{tab1.spins}
\end{center}
\end{table}

Thus, our data was found to be in excellent agreement with the KT theory and, 
in fact, this study provides the first direct measurement of $\nu$ from both 
$\chi$ and $\xi$ data that is consistent with the KT predictions.

\subsection{O(3) Model}\index{O(3) Model}

The XY model is the simplest O(N) model, having $N=2$; the O(N) model being a 
set of rotors (N-component continuous valued spins) on an N-sphere.  For $N 
\ge 3$, this model is asymptotically free \cite{Polyakov:75a}, and for $N=3$, 
there exist instanton solutions.  Some of these properties are analogous to 
those of gauge theories in four dimensions; hence, these models are 
interesting.  In particular, the O(3) model\index{O(3) model} in two 
dimensions should shed 
some light on the asymptotic freedom of QCD (SU(3)) in four dimensions.
The predictions of the renormalization group for the susceptibility $\chi$
and inverse correlation length ({\em i.e.,\/} mass gap) $m$ in the O(3) model 
are \cite{Brezin:76a}
\begin{equation}
\chi = C \beta^{-4} \exp\,(4\pi\beta )\, (1 + O(1/\beta))
\end{equation}
and
\begin{equation}
m = \frac{1}{\xi} = C^\prime \beta \exp\,(-2\pi\beta )\,(1 + O(1/\beta)) ,
\end{equation}
respectively.
If $m$ and $\chi$ vary according to these equations, without the correction of
order $1/\beta$, they are said to follow asymptotic scaling.  Previous work 
was able to confirm that this picture is qualitatively correct, but was not 
able to probe deep enough in the area of large correlation lengths to obtain 
good agreement. 
 
The combination of the over-relaxed algorithm and the computational power of 
the FPS T-Series allowed us to simulate lattices of size up to $1024^2$.  We
were thus, able to simulate at coupling constants that correspond to 
correlation lengths up to $300$, on lattices where finite size effects are 
negligible.  We were also able to gather large statistics and, thus, obtain 
small statistical errors.  Our simulation is in good agreement with the 
similar cluster calculations \cite{Wolff:89b} \cite{Wolff:90a}.  Thus, we 
have validated and extended these results in a regime where our algorithm is 
the only known alternative to clustering.

We have made extensive runs at 10 values of the coupling constant.  At the 
lowest $\beta$, several hundred thousand sweeps were collected while for the 
largest values of $\beta$, between $50,000$ and $100,000$ sweeps were made.  
Each sweep consists of between 10 iterations through the lattice at the 
former end and 150 iterations at the latter.  The statistics we have gathered 
are equivalent to about 200 days, using the full 128-node FPS machine.

Our results for the correlation length and susceptibility for each coupling
constant $\beta$ and lattice size, are shown in Table~\ref{tab2.spins}.  The 
autocorrelation times are also shown.  The quantities measured on different 
sized lattices at the same $\beta$ agree, showing  that the infinite volume 
limit has been reached.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{Beta}&\multicolumn{1}{c|}{L,T}&
\multicolumn{1}{c|}{$\tau$}&\multicolumn{1}{c|}{$N_{or}$}&
\multicolumn{1}{c|}{$\chi$}&\multicolumn{1}{c|}{$\xi$}
\\
\hline
\medStrut
1.50 & "256     & 1.3 & "12 & $""176.4 \pm ""0.2$ & "11.1 \\
1.60 & "256     & 2.5 & "12 & $""448.4 \pm ""0.7$ & "19.0 \\
1.70 & "256     & 7.0 & "12 & $"1263\phantom{.}"\pm ""3$\phantom{.}" & "34.4 \\
     & "512     & 1.8 & "35 & $"1264\phantom{.}"\pm ""3$\phantom{.}" & "34.4 \\
1.75 & "768     & 1.8 & "50 & $"2197\phantom{.}"\pm "15$\phantom{.}" & "47.2 \\
1.80 & "512     & 3.5 & "40 & $"3850\phantom{.}"\pm "10$\phantom{.}" & "64.7 \\
     & "768     & 3.0 & "45 & $"3820\phantom{.}"\pm "20$\phantom{.}" & "64.5 \\
1.85 & "768     & 2.3 & "80 & $"6730\phantom{.}"\pm "25$\phantom{.}" & "88.7 \\
1.90 & 1024, 512& 2.0 & 150 & $11600\phantom{.}"\pm "60$\phantom{.}" & 121.5 \\
     & 1024     & 1.5 & 200 & $11870\phantom{.}"\pm "60$\phantom{.}" & 122.7 \\
1.95 & 1024     & 1.7 & 200 & $20640\phantom{.}"\pm 310$\phantom{.}" & 164.8 \\
2.00 & 1024     & 1.7 & 250 & $35100\phantom{.}"\pm 400$\phantom{.}" & 224.3 \\
2.05 & 1024     & 1.7 & 300 & $56220\phantom{.}"\pm 550$\phantom{.}" & 295.6 \\
\hline
\end{tabular}
\caption{Coupling constant, lattice size, autocorrelation time, number of 
over-relaxed sweeps, susceptibility and correlation length for the O(3) 
model.}
\label{tab2.spins}
\end{center}
\end{table}

To compare the behavior of the correlation length and susceptibility
with the asymptotic scaling predictions, we use the ``correlation length
defect'' $\delta_{\xi}$ and ``susceptibility defect'' $\delta_{\chi}$
which are defined as follows\@: $\delta_{\xi} = \beta (e^{-2\pi\beta})\xi$, 
$\delta_{\chi} = \beta^4 (e^{-4\pi\beta})\chi$ so that asymptotic
scaling is seen if $\delta_{\xi}$, $\delta_{\chi}$ go to constants as
$\beta\rightarrow\infty$.  These defects are shown in
Figure~\ref{fig11.spins} and Figure~\ref{fig12.spins} respectively.  It
is clear that asymptotic scaling does not set in for $\beta < 2$, but it
is not possible to drawn a clear conclusions for $\beta > 2$---though
the trends of the last two or three points may be toward constant
behvaior.

\begin{figure}
\vspace*{2.0in}
\caption{Correlation Length Defect vs. the Coupling Constant for the O(3) 
Model}
\label{fig11.spins}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Susceptibility Defect vs. the Coupling Constant for the O(3) Model}
\label{fig12.spins}
\end{figure}

We gauged the speed of the algorithm in producing statistically independent 
configurations by measuring the autocorrelation time $\tau$.  We used this to 
estimate the dynamical critical exponent $z$, which is defined by $\tau = c
\cdot\xi^z$.  For constant $N_{or}=12$, our fits give $z=1.33(1)$.  However, 
we discovered that by increasing $N_{or}$ in rough proportion to $\xi$, we 
can improve the performance of the algorithm significantly.  To compare the
speed of decorrelation between runs with different $N_{or}$, we define a new 
quantity, which we call `effort', $e =  N_{or} \times\tau$.  This measures 
the computational effort expended to obtain a decorrelated configuration.  We 
define a new exponent $z^\prime$  from $e\sim\xi^{z^\prime}$, where $N_{or}$ 
is chosen to keep $\tau$ constant.  We also found that the behavior of
the decorrelation time can be approximated over a good range by
\begin{equation}
\tau  = C^{\prime\prime}\cdot\xi^z \cdot N_{or}^{-z/z^\prime}.
\end{equation}
A fit to the set of points ($N_{or}$, $\xi$, $\tau > 1.0$) gives
$z = 1.301\pm 0.012$, $z^{\prime} = 1.079\pm 0.010$.  Thus, $z^{\prime}$ is
significantly lower than $z$.  Figure~\ref{fig13.spins} shows $\tau$ vs
$N_{or}$ with the fits shown as solid lines.

\begin{figure}
\vspace*{2.0in}
\caption{Decorrelation Time $\tau$ vs. Number of Overrelations Sweeps
$N_{or}$ for Different Values of $\beta$}
\label{fig13.spins}
\end{figure}

\section{An Automata Model of Granular Materials}\index{Granular
materials}
\label{grain}

\subsection{Introduction}

Physical systems comprised of discrete, macroscopic particles or grains,
which are not bonded to one another, are important in civil, chemical, and 
agricultural engineering, as well as in natural geological and planetary 
environments.  Granular systems are observed in rock slides, sand\index{sand} 
dunes, clastic sediments, snow avalanches, and planetary rings. 
In engineering and industry they are found in connection with the processing 
of cereal grains, coal, gravel, oil shale, and powders, and are well-known to 
pose important problems associated with the movement of sediments by streams,
rivers, waves, and the wind.

The standard approach to the theoretical modelling of multiparticle systems 
in physics has been to treat the system as a continuum and to formulate the 
model in terms of differential equations.  As an example, the science of soil 
mechanics has traditionally focused mainly on quasi-static granular systems, 
a prime objective being to define and predict the conditions under which 
failure of the granular soil system will occur.  Soil mechanics is a 
macroscopic continuum model requiring an explicit constitutive law relating, 
for example, stress and strain.  While very successful for the low-strain 
quasi-static applications for which it is intended, it is not clear how it 
can be generalized to deal with the high-strain, explicitly time-dependent 
phenomena which characterize a great many other granular systems of interest.
Attempts at obtaining a generalized theory of granular systems using a
differential equation formalism \cite{Johnson:87a} have met with limited 
success.

An alternate approach to formulating physical theories can be found in
the concept of cellular automata\index{cellular automata}, which was first 
proposed by Von~Neumann in 1948.  In this approach, the space of a physical 
problem would be divided up into many small, identical cells each of which 
would be in one of a finite number of states.  The state of a cell would 
evolve according to a rule that is both local (involves only the cell itself 
and nearby cells) and universal (all cells are updated simultaneously using 
the same rule).

The Lattice Grain Model\index{Lattice Grain Model} \cite{Gutt:89a} (LGrM) we 
discuss here is a microscopic, explicitly time-dependent, cellular automata 
model, and can be applied naturally to high-strain events.  LGrM carries some 
attributes of both particle dynamics\index{particle dynamics} models 
\cite{Cundall:79a}, \cite{Haff:87a}, \cite{Haff:87b}, \cite{Walton:84a}, 
\cite{Werner:87a} (PDM), 
which are based explicitly on Newton's second law, and lattice gas models 
\cite{Frisch:86a}, \cite{Margolis:86a} (LGM), in that its fundamental element 
is a discrete particle, but differs from these substantially in detail.  Here 
we describe the essential features of LGrM, compare the model with both PDM 
and LGM, and finally discuss some applications.

\subsection{Comparison to Particle Dynamics Models}

The purpose of the lattice grain model is to predict the behavior of
large numbers of grains (10,000 to 1,000,000) on scales much larger
than a grain diameter.  In this respect, it goes beyond particle dynamics 
calculations, which are limited to no more than $\sim 10,000$ grains by 
currently available computing resources \cite{Cundall:79a}, \cite{Haff:87a}, 
\cite{Haff:87b}, \cite{Walton:84a}, \cite{Werner:87a}.  The particle dynamics 
models follow the motion of each individual grain exactly, and may be 
formulated in one of two ways depending upon the model adopted for 
particle-particle interactions.

In one formulation, the interparticle contact times are assumed to be of 
finite duration, and each particle may be in simultaneous contact with 
several others \cite{Cundall:79a}, \cite{Haff:87a}, \cite{Walton:84a}, 
\cite{Werner:87a}.  Each particle obeys Newton's law, $F = ma$, and a 
detailed integration of the equations of motion of each particle is performed.
In this form, while useful for applications involving a much smaller number 
of particles than LGrM allows, PDM cannot compete with LGrM for systems 
involving large numbers of grains because of the complexity of PDM 
``automata''.

In the second, simpler formulation, the interparticle contact times are
assumed to be of infinitesimal duration, and particles undergo only binary 
collisions (the hard-sphere collisional models) \cite{Haff:87b}.  Hard-sphere 
models usually rely upon a collision-list ordering of collision events to 
avoid the necessity of checking all pairs of particles for overlaps at each 
time step.  In regions of high particle number density, collisions are very
frequent; and thus in problems where such high density zones appear,
hard-sphere models spend most of their time moving particles through very 
small distances using very small time steps.  In granular flow, zones of 
stagnation where particles are very nearly in contact much of the time are 
common, and the hard-sphere model is therefore unsuitable, at least in its 
simplest form, as a model of these systems.  LGrM avoids these difficulties 
because its time-stepping is controlled, not by a collision list but by a scan 
frequency, which in turn is a function of the speed of the fastest particle 
and is independent of number density.  Furthermore, although fundamentally a 
collisional model, LGrM can also mimic the behavior of consolidated or 
stagnated zones of granular material in a manner which will be described. 

\subsection{Comparison to Lattice Gas Models}

LGrM closely resembles LGM \cite{Frisch:86a}, \cite{Margolis:86a} in some 
respects.  First, for two-dimensional applications, the region of space in 
which the particles are to move is discretized into a triangular 
lattice-work, upon each node of which can reside a particle.  The particles 
are capable of moving to neighboring cells at each tick of the clock, subject 
to certain simple rules.  Finally, two particles arriving at the same cell 
(LGM), or adjacent cells (LGrM), at the same time, may undergo a ``collision'' 
in which their outgoing velocities are determined according to specified 
rules chosen to conserve momentum.

Each of the particles in LGM has the same magnitude of velocity and is
allowed to move in one of six directions along the lattice, so that each
particle travels exactly one lattice spacing in each time step.  The single 
velocity magnitude means that all collisions between particles are perfectly 
elastic and that energy conservation is maintained simply through particle 
number conservation.  It also means that the temperature of the gas is 
uniform throughout time and space, thus limiting the application of LGM to 
problems of low Mach number.  An exclusion principle is maintained in which 
no two particles of the same velocity may occupy one lattice point.  Thus, 
each lattice point may have no more than six particles, and the state of a 
lattice point can be recorded using only six bits.

LGrM differs from LGM in that has many possible velocity states, not just six.
In particular, in LGrM not only the direction but the magnitude of the
velocity can change in each collision.  This is a necessary condition because 
the collision of two macroscopic particles is always inelastic, so that 
mechanical energy is not conserved.  The LGrM particles satisfy a somewhat 
different exclusion principle\@: no more than one particle at a time may 
occupy a single site.  This exclusion principle allows LGrM to capture some 
of the volume-filling properties of granular material, in particular, to be 
able to approximate the behavior of static granular masses.

The determination of the time step is more critical in LGrM than in LGM.
If the time step is long enough that some particles travel several lattice 
spacings in one clock tick, the problem of finding the intersection of 
particle trajectories arises.  This involves much computation and 
defeats the purpose of an automata approach.  A very short time step would 
imply that most particles would not move even a single lattice spacing.
Here we choose a time step such that the fastest particle will move
exactly one lattice spacing.  A ``position offset'' is stored for each of the 
slower particles, which are moved accordingly when the offset exceeds 
one-half lattice spacing.  These extra requirements for LGrM automata imply a 
slower computation speed than expected in LGM simulations\@; but, as a 
dividend, we can compute inelastic grain flows of potential engineering and 
geophysical interest.

\subsection{The Rules for the Lattice Grain Model}

In order to keep the particle-particle interaction rules as simple as
possible, all interparticle contacts, whether enduring contacts or true
collisions, will be modeled as collisions.  Collisions that model 
enduring contacts will transmit, in each time step, an impulse equal to the 
force of the enduring contact times the time step.  The fact that collisions 
take place between particles on adjacent lattice nodes means that some 
particles may undergo up to six collisions in a time step.  For simplicity, 
these collisions will be resolved as a series of binary collisions.  The 
order in which these collisions are calculated at each lattice node, as well 
as the order in which the lattice nodes are scanned, is now an important 
consideration.

The rules of the Lattice Grain Model may be summarized as follows:
\begin{enumerate}
\item The particles reside on the nodes of a two-dimensional triangular 
lattice, obeying the exclusion principle that no node may have more than one 
particle.
\item Each particle has two components of velocity, which may take on any 
value.  At the beginning of each time step, each particle's velocity is 
incremented due to the acceleration of gravity.
\item The size of each time step is set so that the fastest particle will 
travel one lattice spacing in that time step.
\item Two components of a ``position offset'' are maintained for each 
particle.  This offset is incremented after the velocities in each time step 
according to gravitational acceleration and the particle's velocity\@:
\begin{equation}
\Delta q_i = v_i \Delta t + \frac{1}{2} g_i \Delta t^2
\end{equation}
where\@:
\begin{center}
\begin{tabular}{rl}
$i$ & = 1,2; \\
$\Delta q_i$ & = i\@{\em th} component of increment in position offset; \\
$v_i$ & = i\@{\em th} component of particle velocity; \\
$g_i$ & = i\@{\em th} component of gravitational acceleration; \\
$\Delta t$ & = current time step. \\
\end{tabular}
\end{center}
Once the offset exceeds half the distance to the nearest lattice node, and 
that node is empty, the particle is moved to that node, and its offset is 
decremented appropriately.  Also, in a collision, the component of the offset 
along the line connecting the centers of the colliding particles is set to 
zero.
\item The order in which the lattice is scanned is chosen so as not to create 
a coupling between the scan pattern and the particle motions.  Thus, the 
particle position updates are done on every third lattice point of every 
third row, with this pattern being repeated nine times to cover all 
lattice sites.
\item Particle collisions are calculated assuming that they are smooth hard 
disks with a given coefficient of restitution.  Particles on adjacent nodes 
are assumed to collide if their relative velocity is bringing them together.
The following order has been adopted for evaluating possible collisions on 
odd time steps: 3b, 3c, 3f, 2f, 2c, 2b, 4b, 4c, 4f, 1f, 1c, 1b; and for even 
time steps: 1b, 1c, 1f, 4f, 4c, 4b, 2b, 2c, 2f, 3f, 3c, 3b (where the lattice 
numbers and collision directions are defined in Figure~\ref{fig1.grain}).
\item In order to incorporate a container, wall, or other barrier within 
these rules, a second type of particle is introduced---the wall particle.  
This particle is similar to the movable particles, and interacts with them 
through binary collisions (with a separately defined inelasticity), but is 
regarded as having infinite mass.  To allow for the introduction of shearing 
motion from a wall (as in a Couette flow problem), the particles making up the
wall are given a common constant velocity, which is used in the usual fashion 
for calculating the results of collisions.  However, the position of the wall 
particles in the lattice remains fixed throughout the simulation.
\item Though a single particle does not accurately predict the trajectory of 
a single grain, we still regard each particle as representing one grain when 
we are extracting information from the simulation regarding the behavior of 
groups of grains.  Thus, the size of one particle, as well as the spacing 
between lattice points, is taken to be one grain-diameter.
\end{enumerate}

The transmission of ``static'' contact forces within a mass of grains (as in 
grains at rest in a gravitational field) is handled naturally within the 
above framework.  Though a particle in a static mass of grains may be 
nominally at rest, its velocity may be nonzero (due to gravitational or 
pressure forces); and it will transmit the appropriate force (in the form of 
an impulse) to the particles under it by means of collisions.  When these 
impulses are averaged over several time steps, the proper weights and 
pressures will emerge.

\begin{figure}
\vspace*{1.5in}
\caption{Definition of Lattice Numbers and Collision Directions}
\label{fig1.grain}
\end{figure}

\subsection{Implementation on a Parallel Processor Computer}

When implementing this algorithm on a computer, what is stored in the
computer's memory is information concerning each point in the lattice,
regardless of whether or not there is a particle at that lattice point.
This allows for very efficient checking of the space around each
particle for the presence of other particles ({\em i.e.,} information
concerning the six adjacent points in a triangular lattice will be found at 
certain known locations in memory).  The need to keep information on empty 
lattice points in memory does not entail as great a penalty as might be 
thought; many lattice grain model problems involve a high density of 
particles, typically one for every one to four lattice points, and the memory 
cost per lattice point is not large.  The memory requirements for the 
implementation of LGrM as described here are five variables per lattice
site\@: two components of position, two components of velocity, and one 
status variable, which denotes an empty site, an occupied site, or a bounding 
``wall'' particle.  If each variable is stored using four bytes of memory, 
then each lattice point requires 20 bytes.

The standard configuration for a simulation consists of a lattice with a
specified number of rows and columns bounded at the top and bottom by two 
rows of wall particles (thus forming the top and bottom walls of the problem 
space), and with left and right edges connected together to form periodic 
boundary conditions.  Thus, the boundaries of the lattice are handled 
naturally within the normal position updating and collision rules, with very 
little additional programming.  (Note\@: since the gravitational acceleration 
can point in an arbitrary direction, the top and bottom walls can become side 
walls for chute flow.  Also, the periodic boundary conditions can be broken 
by the placement of an additional wall, if so desired.)

Because of the nearest-neighbor type interactions involved in the model, the 
computational scheme was well-suited to an nCUBE parallel processor.  For the 
purpose of dividing up the problem, the hypercube architecture is unfolded 
into a two-dimensional array, and each processor is given a roughly 
equal-area section of the lattice.  The only interaction between sections 
will be along their common boundaries; thus, each processor will only need to 
exchange information with its eight immediate neighbors.  The program itself 
was written in C under the Cubix/CrOSIII operating system.

\subsection{Simulations}

The LGrM simulations performed so far have involved from $\sim 10^4$ to
$10^6$ automata.  Trial application runs included 2D, vertical, 
time-dependent flows in several geometries, of which two examples are given 
here\@: Couette flow\index{Couette flow} and flow out of an hourglass-shaped 
hopper.

The standard Couette flow configuration consists of a fluid confined between 
two flat parallel plates of infinite extent, without any gravitational 
accelerations.  The plates move in opposite directions with velocities that 
are equal and that are parallel to their surfaces, which results in the
establishment of a velocity gradient and a shear stress in the fluid.  For 
fluids that obey the Navier-Stokes equation, an analytical solution is 
possible in which the velocity gradient and shear stress are constant across 
the channel.  If, however, we replace the fluid by a system of inelastic 
grains, the velocity gradient will no longer necessarily be constant across 
the channel.  Typically, stagnation zones or plugs form in the center of the 
channel with thin shear-bands near the walls.  Shear-band formation in 
flowing granular materials was analyzed earlier by Haff and others
\cite{Haff:83a}, \cite{Hui:84a} based on kinetic theory models.

The simulation was carried out with 5760 grains, located in a channel 60 
lattice points wide by 192 long.  Due to the periodic boundary conditions at 
the left and right ends, the problem is effectively infinite in length.  The 
first simulation is intended to reproduce the standard Couette flow for a 
fluid; consequently, the particle-particle collisions were given a coefficient 
of restitution of 1.0 ({\em i.e.,} perfectly elastic collisions) and the 
particle-wall collisions were given a .75 coefficient of restitution.  The 
inelasticity of the particle-wall collisions is needed to simulate the 
conduction of heat (which is being generated within the fluid) from the fluid 
to the walls.  The simulation was run until an equilibrium was established in 
the channel (Figure~\ref{fig2.grain}(a)).  The average $x$- and 
$y$-components of velocity and the second moment of velocity, as functions of 
distance across the channel, are plotted in Figure~\ref{fig2.grain}(b).

\begin{figure}
\vspace*{4.0in}
\caption{(a) Elastic Particle Couette Flow\ \ (b) X-component (1), 
y-component (2), and Second Moment (3) of Velocity}
\label{fig2.grain}
\end{figure}

The second simulation used a coefficient of restitution of .75 for both the 
particle-particle and particle-wall collisions.  The equilibrium results are 
shown in Figure~\ref{fig3.grain}(a) and Figure~\ref{fig3.grain}(b).  As can 
be seen from the plots, the flow consists of a central region of particles 
compacted into a plug, with each particle having almost no velocity.  Near 
each of the moving walls, a region of much lower density has formed in which 
most of the shearing motion occurs.  Note the increase in value of the second 
moment of velocity (the granular ``thermal velocity'') near the walls, 
indicating that grains in this area are being ``heated'' by the high rate of 
shear.  It is interesting to note that these flows are turbulent in the sense
that shear stress is a quadratic, not a linear, function of shear rate.

\begin{figure}
\vspace*{4.0in}
\caption{(a) Inelastic Particle Couette Flow\ \ (b) X-component (1), 
y-component (2), and Second Moment (3) of Velocity.}
\label{fig3.grain}
\end{figure}

In the second problem, the flow of grains through a hopper or an hourglass 
with an opening only a few grain diameters wide, was studied; the driving 
force was gravity.  This is an example of a granular system which contains a 
wide range of densities, from groups of grains in static contact with one 
another to groups of highly agitated grains undergoing true binary collisions.
Here, the number of particles used was 8310; and the lattice was 240 points 
long by 122 wide.  Additional walls were added to form the sloped sides of 
the bin and to close off the bottom of the lattice to prevent the 
periodic boundary conditions from reintroducing the falling particles back 
into the bin (Figure~\ref{fig4.grain}(a)).  This is a typical feature of 
automata modeling.  It is often easier to configure the simulation to 
resemble a real experiment---in this case by explicitly ``catching'' spent 
grains---than by re-programming the basic code to erase such particles.

The hourglass flow, Figure~\ref{fig4.grain}(b), showed internal shear zones, 
regions of stagnation, free-surface evolution toward an angle of repose, and 
an exit flow rate approximately independent of pressure head, as observed
experimentally \cite{Tuzun:82a}.  It is hard to imagine that one could solve 
a partial differential equation describing such a complex, multiple-domain, 
time-dependent problem, even if the right equation were known (which is not 
the case).

\begin{figure}
\vspace*{5.0in}
\caption{(a) Initial Condition of Hourglass\ \ (b) Hourglass Flow after 2048 
Time Steps}
\label{fig4.grain}
\end{figure}

\subsection{Conclusion}

These exploratory numerical experiments show that an automata approach to 
granular dynamics problems can be implemented on parallel computing machines.
Further work remains to be done to assess more quantitatively how well such 
calculations reflect the real world, but the prospects are intriguing.

\subsection{Acknowledgements}

This work was supported in part by USDOE [DE-FG22-86-PC90959] and by
USARO [DAAL03-86-K-0123] and [DAAL03-89-K-0089].  We thank Geoffrey Fox for 
making available to us the resources of the Caltech Concurrent Computation 
Project and Tom Tombrello for his interest and support.

