new

[book_gpu.git] / BookGPU / Chapters / chapter14 / ch14.tex

diff --git a/BookGPU/Chapters/chapter14/ch14.tex b/BookGPU/Chapters/chapter14/ch14.tex
index 916717df33fc168f80bffd43948fef61a3cf5b08..2ce9331300ab9c7fbaa4e4182c9edf49f10b96c0 100755 (executable)
--- a/BookGPU/Chapters/chapter14/ch14.tex
+++ b/BookGPU/Chapters/chapter14/ch14.tex
@@ -1,25 +1,26 @@
-\chapterauthor{Alan Gray and Kevin Stratford}{EPCC, The University of Edinburgh}
+\chapterauthor{Alan Gray and Kevin Stratford}{EPCC, The University of Edinburgh, United Kingdom}

 
 

-\chapter{Ludwig: multiple GPUs for a complex fluid lattice Boltzmann
+\chapter[Ludwig: multiple GPUs for a fluid lattice Boltzmann
+application]{Ludwig: multiple GPUs for a complex fluid lattice Boltzmann

 application}
 
 %\putbib[biblio]
 
 
 \section{Introduction}
 application}
 
 %\putbib[biblio]
 
 
 \section{Introduction}

-The lattice Boltzmann (LB) method (for an overview see, e.g.,
+The lattice Boltzmann (LB) method \index{Lattice Boltzmann method} (for an overview see, e.g.,

 \cite{succi-book}) has become a popular approach to a variety of fluid
 \cite{succi-book}) has become a popular approach to a variety of fluid

-dynamics problems.  It provides a way to solve the incompressible,
+dynamics problems.  It provides a way to solve the incompressible

 isothermal Navier-Stokes equations and has the attractive features of
 being both explicit in time and local in space. This makes the LB
 isothermal Navier-Stokes equations and has the attractive features of
 being both explicit in time and local in space. This makes the LB

-method well-suited to parallel computation. Many efficient parallel
+method well suited to parallel computation. Many efficient parallel

 implementations of the LB method have been undertaken, typically using
 a combination of distributed domain decomposition and the Message
 Passing Interface (MPI). However, the potential
 implementations of the LB method have been undertaken, typically using
 a combination of distributed domain decomposition and the Message
 Passing Interface (MPI). However, the potential

-performance benefits offered by GPUs has motivated a new `mixed-mode'
+performance benefits offered by GPUs has motivated a new ``mixed-mode''

 approach to address very large problems. Here, fine-grained
 parallelism is implemented on the GPU, while MPI is reserved for
 approach to address very large problems. Here, fine-grained
 parallelism is implemented on the GPU, while MPI is reserved for

-larger-scale parallelism.  This mixed mode is of increasing interest
+larger scale parallelism.  This mixed mode is of increasing interest

 to application programmers at a time when many supercomputing services
 are moving
 toward clusters of GPU accelerated nodes. The design questions which
 to application programmers at a time when many supercomputing services
 are moving
 toward clusters of GPU accelerated nodes. The design questions which


@@ -55,12 +56,12 @@ and particle suspensions, and typically require additional physics
 beyond the bare Navier-Stokes equations to provide a full
 description~\cite{aidun2010}. The representation of this extra
 physics raises additional design questions for the application
 beyond the bare Navier-Stokes equations to provide a full
 description~\cite{aidun2010}. The representation of this extra
 physics raises additional design questions for the application

-programmer. Here, we consider the \textit{Ludwig} code \cite{desplat},
+programmer. Here, we consider the \textit{Ludwig} code \cite{desplat}\index{Ludwig code},

 an LB application developed specifically for complex fluids
 (\textit{Ludwig} was named for Boltzmann, 1844--1906).
 We will present the steps
 required to allow \textit{Ludwig} to exploit efficiently both a single
 an LB application developed specifically for complex fluids
 (\textit{Ludwig} was named for Boltzmann, 1844--1906).
 We will present the steps
 required to allow \textit{Ludwig} to exploit efficiently both a single

-GPU, and also many GPUs in parallel.  We show that
+GPU and many GPUs in parallel.  We show that

 \textit{Ludwig} scales excellently to at least the one thousand GPU level
 (the largest resource available at the time of writing) with
 indications that the code will scale to much larger systems as
 \textit{Ludwig} scales excellently to at least the one thousand GPU level
 (the largest resource available at the time of writing) with
 indications that the code will scale to much larger systems as


@@ -68,7 +69,7 @@ they become available. In addition, we consider the steps required
 to represent fully-resolved moving solid particles (usually referred
 to as colloids in the context of complex fluids). Such particles
 need to have their surface resolved on the lattice scale in order
 to represent fully-resolved moving solid particles (usually referred
 to as colloids in the context of complex fluids). Such particles
 need to have their surface resolved on the lattice scale in order

-to include relevant surface physics, and must be able to move: e.g.,
+to include relevant surface physics and must be able to move, e.g.,

 to execute Brownian motion in response to random forces from the
 fluid. Standard methods are available to represent such particles
 (e.g., \cite{ladd1994,nguyen2002}) which are amenable to effective
 to execute Brownian motion in response to random forces from the
 fluid. Standard methods are available to represent such particles
 (e.g., \cite{ladd1994,nguyen2002}) which are amenable to effective


@@ -81,7 +82,7 @@ In the following section we provide a brief overview of the lattice
 Boltzmann method, and mention some of the general issues which can
 influence performance in the CPU context.
 In Section \ref{ch14:sec:singlegpu}, we then describe the alterations
 Boltzmann method, and mention some of the general issues which can
 influence performance in the CPU context.
 In Section \ref{ch14:sec:singlegpu}, we then describe the alterations

-which are required to exploit the GPU architecture effectively, and
+which are required to exploit the GPU architecture effectively and

 highlight how the resulting code differs from the CPU version. In
 Section \ref{ch14:sec:parallelgpu}, we extend this description to
 include the steps required to allow exploitation of many GPUs in
 highlight how the resulting code differs from the CPU version. In
 Section \ref{ch14:sec:parallelgpu}, we extend this description to
 include the steps required to allow exploitation of many GPUs in


@@ -127,8 +128,8 @@ are chosen to capture the correct symmetries of the Navier-Stokes
 equations. A typical choice, used here, is the so-called D3Q19
 basis in three dimensions where there is one velocity such that
 $\mathbf{c} \Delta t$ is zero, along with six extending to the nearest
 equations. A typical choice, used here, is the so-called D3Q19
 basis in three dimensions where there is one velocity such that
 $\mathbf{c} \Delta t$ is zero, along with six extending to the nearest

-neighbour
-lattice sites, and twelve extending to the next-nearest neighbour sites
+neighbor
+lattice sites, and twelve extending to the next-nearest neighbor sites

 ($\Delta t$ being the discrete time step). The fundamental object
 in LB is then the distribution function $f_i (\mathbf{r};t)$ whose
 moments are related to the local hydrodynamic quantities: the fluid
 ($\Delta t$ being the discrete time step). The fundamental object
 in LB is then the distribution function $f_i (\mathbf{r};t)$ whose
 moments are related to the local hydrodynamic quantities: the fluid


@@ -138,16 +139,16 @@ function is described by a discrete Boltzmann equation
 f_i(\mathbf{r} + \mathbf{c}_i \Delta t; t) - f_i(\mathbf{r}; t) 
 = - {\cal L}_{ij} f_j(\mathbf{r};t).
 \end{equation}
 f_i(\mathbf{r} + \mathbf{c}_i \Delta t; t) - f_i(\mathbf{r}; t) 
 = - {\cal L}_{ij} f_j(\mathbf{r};t).
 \end{equation}

-It is convenient to think of this in two stages. First, the right hand
+It is convenient to think of this in two stages. First, the right-hand

 side represents the action of a collision operator ${\cal L}_{ij}$,
 which is local to each lattice site and relaxes the distribution toward
 a local equilibrium at a rate ultimately related to the fluid viscosity.
 side represents the action of a collision operator ${\cal L}_{ij}$,
 which is local to each lattice site and relaxes the distribution toward
 a local equilibrium at a rate ultimately related to the fluid viscosity.

-Second, the left hand side represents a propagation step (sometimes referred
+Second, the left-hand side represents a propagation step (sometimes referred

 to as streaming step), in which each element $i$ of the distribution is
 displaced $\mathbf{c}_i \Delta t$, i.e., one lattice spacing in the
 appropriate direction per discrete time step. 
 
 to as streaming step), in which each element $i$ of the distribution is
 displaced $\mathbf{c}_i \Delta t$, i.e., one lattice spacing in the
 appropriate direction per discrete time step. 
 

-More specifically, we store a vector of 19 double precision floating
+More specifically, we store a vector of 19 double-precision floating

 point values at each lattice site for the distribution function
 $f_i(\mathbf{r};t)$.
 The collision operation ${\cal L}_{ij}$, which is local at each lattice
 point values at each lattice site for the distribution function
 $f_i(\mathbf{r};t)$.
 The collision operation ${\cal L}_{ij}$, which is local at each lattice


@@ -156,18 +157,34 @@ ${\cal M}_{ij}f_j$ is used to transform the distributions into the
 hydrodynamic quantities, where ${\cal M}_{ij}$ is a constant 19x19
 matrix related to the choice of
 $\mathbf{c}_i$. The non-conserved hydrodynamic quantities are then
 hydrodynamic quantities, where ${\cal M}_{ij}$ is a constant 19x19
 matrix related to the choice of
 $\mathbf{c}_i$. The non-conserved hydrodynamic quantities are then

-relaxed toward their (known) equilibrium values, and are transformed
+relaxed toward their (known) equilibrium values and are transformed

 back to new post-collision distributions via the inverse transformation
 back to new post-collision distributions via the inverse transformation

-${\cal M}^{-1}_{ij}$. This gives rise to the need for a minimum of 2x19$^2$
+${\cal M}^{-1}_{ij}$. This gives rise to the need for a minimum of $2\times 19^2$

 floating point multiplications per lattice site. (Additional operations are
 required to implement, for example, the force $\mathbf{f}(\mathbf{r})$.)
 In contrast, the
 propagation stage consists solely of moving the distribution values
 floating point multiplications per lattice site. (Additional operations are
 required to implement, for example, the force $\mathbf{f}(\mathbf{r})$.)
 In contrast, the
 propagation stage consists solely of moving the distribution values

-between lattice sites, and involves no floating point operations.
+between lattice sites and involves no floating point operations.
+
+
+\begin{figure}[tbp]
+\centering
+\includegraphics[width=12cm]{Chapters/chapter14/figures/decomphalo}
+\caption[The lattice is decomposed between MPI tasks.]{Left: the lattice is decomposed between MPI tasks. For
+  clarity we show a 2D decomposition of a 3D lattice, but in practice
+  we decompose in all 3 dimensions. Halo cells are added to each
+  subdomain (as shown on the upper right for a single slice) which store
+  data retrieved from remote neighbors in the halo exchange. Lower
+  right: the D3Q19 velocity set resident on a lattice site;
+  highlighted are the 5 ``outgoing'' elements to be transferred in a
+  specific direction.}
+\label{ch14:fig:decomphalo}
+\end{figure}
+

 
 In the CPU version, the \textit{Ludwig} implementation stores one time
 level of distribution values $f_i(\mathbf{r}; t)$. This distribution
 
 In the CPU version, the \textit{Ludwig} implementation stores one time
 level of distribution values $f_i(\mathbf{r}; t)$. This distribution

-is stored as a flat array of C data type \texttt{double}, and laid out
+is stored as a flat array of C data type \texttt{double} and laid out

 so that all elements of the velocity are contiguous at a given site
 (often referred to as ``array-of-structures'' order). This is
 appropriate for sums over the distribution required locally at the
 so that all elements of the velocity are contiguous at a given site
 (often referred to as ``array-of-structures'' order). This is
 appropriate for sums over the distribution required locally at the


@@ -176,24 +193,27 @@ Listing~\ref{ch14:listing1}: the fact that consecutive loads are from
 consecutive memory addresses allows the prefetcher to engage fully.
 (The temporary scalar \texttt{a\_tmp} allows caching of the
 intermediate accumulated value in the innermost loop.)  A variety of
 consecutive memory addresses allows the prefetcher to engage fully.
 (The temporary scalar \texttt{a\_tmp} allows caching of the
 intermediate accumulated value in the innermost loop.)  A variety of

-standard sequential optimisations are relevant for the collision stage
+standard sequential optimizations are relevant for the collision stage

 (loop unrolling, inclusion of SIMD operations, and so on
 \cite{wellein2006}).  For example, the collision stage in
 \textit{Ludwig} includes explicit SIMD code, which is useful if the
 (loop unrolling, inclusion of SIMD operations, and so on
 \cite{wellein2006}).  For example, the collision stage in
 \textit{Ludwig} includes explicit SIMD code, which is useful if the

-compiler on a given platform cannot identify it.  The propagation
-stage is separate, and is organised as a ``pull'' (again, various
-optimisation have been considered, e.g.,
-\cite{pohl2003,mattila2007,wittmann2012}).  No further optimisation is
+compiler on a given platform cannot identify the appropriate vectorization.
+The propagation
+stage is separate and is organized as a ``pull'' (again, various
+optimizations have been considered, e.g.,
+\cite{pohl2003,mattila2007,wittmann2012}).  No further optimization is

 done here beyond ensuring that the ordering of the discrete velocities
 allows memory access to be as efficient as possible. While these
 done here beyond ensuring that the ordering of the discrete velocities
 allows memory access to be as efficient as possible. While these

-optimisations are important, it should be remembered that for some
+optimizations are important, it should be remembered that for some

 complex fluid problems, the hydrodynamics embodied in the LB
 calculation is a relatively small part of the total computational cost
 complex fluid problems, the hydrodynamics embodied in the LB
 calculation is a relatively small part of the total computational cost

-(at the 10\%-level in some cases).  This means optimisation effort may
+(at the 10\% level in some cases).  This means optimization effort may

 be better concentrated elsewhere.
 
 be better concentrated elsewhere.
 

+
+

 \begin{lstlisting}[float, label=ch14:listing1,
 \begin{lstlisting}[float, label=ch14:listing1,

-caption = Collision schematic for CPU.]
+caption = collision schematic for CPU.]

 /* loop over lattice sites */
 for (is = 0; is < nsites; is++) {
   ...
 /* loop over lattice sites */
 for (is = 0; is < nsites; is++) {
   ...


@@ -214,44 +234,32 @@ for (is = 0; is < nsites; is++) {
 \end{lstlisting}
 
 
 \end{lstlisting}
 
 

-\begin{figure}[!t]
-\centering
-\includegraphics[width=12cm]{Chapters/chapter14/figures/decomphalo}
-\caption{Left: the lattice is decomposed between MPI tasks. For
-  clarity we show a 2D decomposition of a 3D lattice, but in practice
-  we decompose in all 3 dimensions. Halo cells are added to each
-  sub-domain (as shown on the upper right for a single slice) which store
-  data retrieved from remote neighbours in the halo exchange. Lower
-  right: the D3Q19 velocity set resident on a lattice site;
-  highlighted are the 5 ``outgoing'' elements to be transferred in a
-  specific direction.}
-\label{ch14:fig:decomphalo}
-\end{figure}
+

 
 
 The regular 3D decomposition is illustrated in Fig.~\ref{ch14:fig:decomphalo}.
 
 
 The regular 3D decomposition is illustrated in Fig.~\ref{ch14:fig:decomphalo}.

-Each local sub-domain is surrounded by a halo, or ghost, region of width one
+Each local subdomain is surrounded by a halo, or ghost, region of width one

 lattice site. While the collision is local, elements of the distribution
 must be exchanged at the edges of the domains to facilitate the propagation.
 To achieve the full 3D halo exchange, the standard approach of shifting the
 lattice site. While the collision is local, elements of the distribution
 must be exchanged at the edges of the domains to facilitate the propagation.
 To achieve the full 3D halo exchange, the standard approach of shifting the

-relevant data in each co-ordinate direction in turn is adopted. This
-requires appropriate synchronisation, i.e., a receive in the the first
-co-ordinate direction must be complete before a send in the second direction
+relevant data in each coordinate direction in turn is adopted. This
+requires appropriate synchronization, i.e., a receive in the the first
+coordinate direction must be complete before a send in the second direction

 involving relevant data can take place, and so on. We note that only
 ``outgoing'' elements of the distribution need to be sent at each edge.
 For the D3Q19 model, this reduces the volume of data traffic from 19 to
 5 of the $f_i(\mathbf{r};t)$ per lattice site at each edge. In the CPU
 version, the necessary transfers are implemented in place using
 involving relevant data can take place, and so on. We note that only
 ``outgoing'' elements of the distribution need to be sent at each edge.
 For the D3Q19 model, this reduces the volume of data traffic from 19 to
 5 of the $f_i(\mathbf{r};t)$ per lattice site at each edge. In the CPU
 version, the necessary transfers are implemented in place using

-a vector of appropriately strided MPI datatypes for each direction.
+a vector of MPI datatypes with appropriate stride for each direction.

 
 
 
 

-\section{Single GPU Implementation}\label{ch14:sec:singlegpu}
+\section{Single GPU implementation}\label{ch14:sec:singlegpu}

 
 
 In this section we describe the steps taken to enable \textit{Ludwig}
 for the GPU. There are a number of crucial issues: first, the
 
 
 In this section we describe the steps taken to enable \textit{Ludwig}
 for the GPU. There are a number of crucial issues: first, the

-minimisation of data traffic between host and device; second, the
-optimal mapping of available parallelism onto the architecture and
+minimization of data traffic between host and device; second, the
+optimal mapping of available parallelism onto the architecture; and

 third, the issue of memory coalescing. We discuss each of these in
 turn.
 
 third, the issue of memory coalescing. We discuss each of these in
 turn.
 


@@ -273,12 +281,12 @@ parallelism inherent in the GPU architecture, particularly for
 those matrix-vector operations within the collision kernel.
 The GPU architecture features a hierarchy of parallelism. At the
 lowest level, groups of 32 threads (warps) operate in
 those matrix-vector operations within the collision kernel.
 The GPU architecture features a hierarchy of parallelism. At the
 lowest level, groups of 32 threads (warps) operate in

-lock-step on different data elements: this is SIMD style vector-level
+lock-step on different data elements: this is SIMD-style vector-level

 parallelism. Multiple warps are combined into a thread block (in which
 parallelism. Multiple warps are combined into a thread block (in which

-communication and synchronisation are possible), and multiple blocks can
+communication and synchronization are possible), and multiple blocks can

 run concurrently across the streaming multiprocessors in the GPU
 (with no communication or
 run concurrently across the streaming multiprocessors in the GPU
 (with no communication or

-synchronisation possible across blocks).  To decompose on the GPU, we
+synchronization possible across blocks).  To decompose on the GPU, we

 must choose which part of the collision to assign to which level of
 parallelism.
 
 must choose which part of the collision to assign to which level of
 parallelism.
 


@@ -289,14 +297,14 @@ to all levels of parallelism, i.e., we use a separate CUDA thread for
 every lattice site, and each thread performs the full matrix-vector
 operation. We find that a block size of 256 performs well on current
 devices: we therefore decompose
 every lattice site, and each thread performs the full matrix-vector
 operation. We find that a block size of 256 performs well on current
 devices: we therefore decompose

-the lattice into groups of 256 sites, and assign each group to a block
+the lattice into groups of 256 sites and assign each group to a block

 of CUDA threads. As the matrix ${\cal M}_{ij}$ is constant, it is
 assigned to the fast {\it constant} on-chip device memory.
 
 For the propagation stage, the GPU implementation adds a second time
 of CUDA threads. As the matrix ${\cal M}_{ij}$ is constant, it is
 assigned to the fast {\it constant} on-chip device memory.
 
 For the propagation stage, the GPU implementation adds a second time

-level of distribution values. The data-dependencies inherent in the
+level of distribution values. The data dependencies inherent in the

 propagation mean that the in-place propagation of the CPU version
 propagation mean that the in-place propagation of the CPU version

-cannot be parallelised effectively without the additional time level.
+cannot be parallelized effectively without the additional time level.

 As both time levels may remain resident on the GPU, this is not a
 significant overhead.
 
 As both time levels may remain resident on the GPU, this is not a
 significant overhead.
 


@@ -307,7 +315,7 @@ is only achieved when data are structured such that threads within a
 segment in a single transaction: this is memory coalescing. It can be
 seen that the array-of-structures ordering used for the distribution
 in the CPU code would not be suitable for coalescing; in fact, it would
 segment in a single transaction: this is memory coalescing. It can be
 seen that the array-of-structures ordering used for the distribution
 in the CPU code would not be suitable for coalescing; in fact, it would

-result in serialised memory accesses and relative poor performance.
+result in serialized memory accesses and relative poor performance.

 To meet the coalescing criteria and allow consecutive threads to read
 consecutive memory addresses on the GPU, we transpose the layout of the
 distribution so that, for each velocity component, consecutive sites
 To meet the coalescing criteria and allow consecutive threads to read
 consecutive memory addresses on the GPU, we transpose the layout of the
 distribution so that, for each velocity component, consecutive sites


@@ -315,7 +323,7 @@ are contiguous in memory (``structure-of-arrays'' order). A schematic of
 the GPU collision code is given in Listing~\ref{ch14:listing2}.
 
 \begin{lstlisting}[float, label=ch14:listing2,
 the GPU collision code is given in Listing~\ref{ch14:listing2}.
 
 \begin{lstlisting}[float, label=ch14:listing2,

-caption = Collision schematic for GPU.]
+caption = collision schematic for GPU.]

 /* compute current site index 'is' from CUDA thread and block */
 /* indices and load distribution from "structure-of-arrays" */
 
 /* compute current site index 'is' from CUDA thread and block */
 /* indices and load distribution from "structure-of-arrays" */
 


@@ -333,18 +341,18 @@ defer some further comments on software engineering aspects of the code to
 the summary.
 
 
 the summary.
 
 

-\section{Multiple GPU Implementation}\label{ch14:sec:parallelgpu}
+\section{Multiple GPU implementation}\label{ch14:sec:parallelgpu}

 
 To execute on multiple GPUs, we use the same domain decomposition
 and message passing framework as the CPU version. Within each
 
 To execute on multiple GPUs, we use the same domain decomposition
 and message passing framework as the CPU version. Within each

-sub-domain (allocated to one MPI task) the GPU implementation
+subdomain (allocated to one MPI task) the GPU implementation

 proceeds as described in the previous section. The only additional
 complication is that halo transfers between GPUs must be staged
 proceeds as described in the previous section. The only additional
 complication is that halo transfers between GPUs must be staged

-through the host (in future, direct GPU to GPU data transfers via
+through the host (in the future, direct GPU-to-GPU data transfers via

 MPI may be possible, obviating the need for these steps). This
 MPI may be possible, obviating the need for these steps). This

-means host MPI sends must be preceded by appropriate device to host
-transfers and host MPI receives must be followed by corresponding host
-to device transfers.
+means host MPI sends must be preceded by appropriate device-to-host
+transfers and host MPI receives must be followed by corresponding
+host-to-device transfers.

 
 %To execute on multiple GPUs in parallel we decompose the problem using
 %the existing high-level parallelisaton framework: each MPI task acts
 
 %To execute on multiple GPUs in parallel we decompose the problem using
 %the existing high-level parallelisaton framework: each MPI task acts


@@ -379,18 +387,18 @@ to device transfers.
 
 In practice, this data movement requires additional GPU kernels to
 pack and unpack the relevant data before and after corresponding MPI
 
 In practice, this data movement requires additional GPU kernels to
 pack and unpack the relevant data before and after corresponding MPI

-calls. However, the standard shift algorithm, in which each co-ordinate
+calls. However, the standard shift algorithm, in which each coordinate

 direction is treated in turn, does provide some scope for the
 overlapping of different operations.
 direction is treated in turn, does provide some scope for the
 overlapping of different operations.

-For example, after the data for the first co-ordinate direction have been
+For example, after the data for the first coordinate direction have been

 retrieved by the host, these can be exchanged using MPI between
 hosts at the same time as kernels for packing and retrieving of
 retrieved by the host, these can be exchanged using MPI between
 hosts at the same time as kernels for packing and retrieving of

-data for the second co-ordinate direction are
+data for the second coordinate direction are

 executed. This
 executed. This

-overlapping must respect the synchronisation required to ensure
-that data values at the corners of the sub-domain are transferred
+overlapping must respect the synchronization required to ensure
+that data values at the corners of the subdomain are transferred

 correctly.
 correctly.

-We use a separate CUDA stream for each co-ordinate direction:
+We use a separate CUDA stream for each coordinate direction:

 this allows some of the host-device communication time to be
 effectively ``hidden'' behind the host-host MPI communication,
 resulting in an overall speedup.
 this allows some of the host-device communication time to be
 effectively ``hidden'' behind the host-host MPI communication,
 resulting in an overall speedup.


@@ -414,7 +422,7 @@ CPUs per node, and with nodes interconnected using Cray Gemini technology.
 For GPU results, a Cray XK6 system is used: this is very similar to the
 XE6, but has one CPU per node replaced with an NVIDIA X2090 GPU. Each
 node in the GPU system therefore features a single Opteron CPU acting
 For GPU results, a Cray XK6 system is used: this is very similar to the
 XE6, but has one CPU per node replaced with an NVIDIA X2090 GPU. Each
 node in the GPU system therefore features a single Opteron CPU acting

-as a host to a single GPU. The inter-node interconnect architecture is
+as a host to a single GPU. The internode interconnect architecture is

 the same as for the Cray XE6. The GPU performance tests use a prototype
 Cray XK6 system with 936 nodes (the largest available at the time of writing).
 To provide a fair comparison, we compare scaling on a \textit{per node}
 the same as for the Cray XE6. The GPU performance tests use a prototype
 Cray XK6 system with 936 nodes (the largest available at the time of writing).
 To provide a fair comparison, we compare scaling on a \textit{per node}


@@ -437,7 +445,7 @@ a simulation in terms of accounting budgets, or electricity.
 \begin{figure}[!t]
 \centering
 \includegraphics[width=10cm]{Chapters/chapter14/figures/two_graphs_ud}
 \begin{figure}[!t]
 \centering
 \includegraphics[width=10cm]{Chapters/chapter14/figures/two_graphs_ud}

-\caption{The weak (top) and strong (bottom) scaling of \textit{Ludwig}.  Closed
+\caption[The weak (top) and strong (bottom) scaling of \textit{Ludwig}.]{The weak (top) and strong (bottom) scaling of \textit{Ludwig}.  Closed

   shapes denote results using the CPU version run on the Cray
   XE6 (using two 16-core AMD Interlagos CPUs per node), while open
   shapes denote results using the GPU version on the Cray XK6 (using a
   shapes denote results using the CPU version run on the Cray
   XE6 (using two 16-core AMD Interlagos CPUs per node), while open
   shapes denote results using the GPU version on the Cray XK6 (using a


@@ -464,7 +472,7 @@ a simulation in terms of accounting budgets, or electricity.
 
 Figure \ref{ch14:fig:scaling} shows the results of both weak and
 strong scaling tests (top and bottom panels, respectively). For
 
 Figure \ref{ch14:fig:scaling} shows the results of both weak and
 strong scaling tests (top and bottom panels, respectively). For

-weak scaling, where the local sub-domain size is kept fixed
+weak scaling, where the local subdomain size is kept fixed

 (here a relatively large 196$^3$ lattice), the time taken by
 an ideal code would remain constant. It can be seen that while
 scaling of the CPU version shows a pronounced upward slope, the
 (here a relatively large 196$^3$ lattice), the time taken by
 an ideal code would remain constant. It can be seen that while
 scaling of the CPU version shows a pronounced upward slope, the


@@ -474,9 +482,9 @@ has a factor of 32 fewer MPI tasks per node, so communications
 require a significantly smaller number of larger data transfers.
 The performance
 advantage of the GPU version ranges from a factor of around 1.5 to
 require a significantly smaller number of larger data transfers.
 The performance
 advantage of the GPU version ranges from a factor of around 1.5 to

-around 1.8. Careful study by other authors \cite{williams2011} have
-found the absolute performance (in terms of floating point
-operations per second, or floating point operations per Watt) to
+around 1.8. Careful studies by other authors \cite{williams2011} have
+found the absolute performance (in terms of floating-point
+operations per second, or floating-point operations per Watt) to

 be remarkably similar between architectures.
 
 Perhaps of more interest is the strong scaling picture (lower
 be remarkably similar between architectures.
 
 Perhaps of more interest is the strong scaling picture (lower


@@ -492,18 +500,18 @@ In contrast, the GPU version shows
 reasonably robust scaling even for the smaller system sizes and good
 scaling for the larger systems.
 
 reasonably robust scaling even for the smaller system sizes and good
 scaling for the larger systems.
 

-\section{Moving Solid Particles}\label{ch14:sec:particles}
+\section{Moving solid particles}\label{ch14:sec:particles}

 
 \begin{figure}[t]
 \centering
 %\includegraphics[width=10cm]{Chapters/chapter14/figures/bbl}
 
 \begin{figure}[t]
 \centering
 %\includegraphics[width=10cm]{Chapters/chapter14/figures/bbl}

-\includegraphics[width=10cm]{Chapters/chapter14/figures/colloid_new}
-\caption{
+\includegraphics[width=10cm]{Chapters/chapter14/figures/colloid_new_gray}
+\caption[A two-dimensional schematic picture of spherical particles on the lattice.]{

 A two-dimensional schematic picture of spherical particles on the lattice.
 Left: a particle is allowed
 to move continuously across the lattice, and the position of the
 A two-dimensional schematic picture of spherical particles on the lattice.
 Left: a particle is allowed
 to move continuously across the lattice, and the position of the

-surface defines fluid lattice sites (light blue) and solid lattice
-sites (dark red). The discrete surface is defined by links
+surface defines fluid lattice sites (gray) and solid lattice
+sites (black). The discrete surface is defined by links

 where propagation would intersect the surface (arrows). Note the
 discrete shape of the two particles is different. Right: post-collision
 distributions are reversed at the surface by the process of bounce-back
 where propagation would intersect the surface (arrows). Note the
 discrete shape of the two particles is different. Right: post-collision
 distributions are reversed at the surface by the process of bounce-back


@@ -534,15 +542,15 @@ The introduction of moving solid particles poses an additional hurdle
 to efficient GPU implementation of an LB code such as \textit{Ludwig}.
 In this section, we give a brief overview of the essential features
 of the CPU implementation, and how the considerations raised in the
 to efficient GPU implementation of an LB code such as \textit{Ludwig}.
 In this section, we give a brief overview of the essential features
 of the CPU implementation, and how the considerations raised in the

-previous sections --- maximising parallelism and minimising both
-host to device and GPU to GPU data movement --- shape the design decisions
+previous sections---maximizing parallelism and minimising both
+host-to-device and GPU-to-GPU data movement---shape the design decisions

 for a GPU implementation. We restrict our discussion to a simple fluid
 in this section; additional solid-fluid boundary conditions (e.g.,
 wetting at a fluid-fluid-solid contact line) usually arise elsewhere
 for a GPU implementation. We restrict our discussion to a simple fluid
 in this section; additional solid-fluid boundary conditions (e.g.,
 wetting at a fluid-fluid-solid contact line) usually arise elsewhere

-in the calculation and broadly independent of the hydrodynamic boundary
+in the calculation and are broadly independent of the hydrodynamic boundary

 conditions which are described in what follows.
 
 conditions which are described in what follows.
 

-Moving solid particles (here, spheres) are defined by a centre position
+Moving solid particles (here, spheres) are defined by a center position

 which is allowed to move continuously across the space of the lattice,
 and a fixed radius which is typically a few lattice spacings $\Delta r$.
 The surface of the particle is defined by a series of \textit{links}
 which is allowed to move continuously across the space of the lattice,
 and a fixed radius which is typically a few lattice spacings $\Delta r$.
 The surface of the particle is defined by a series of \textit{links}


@@ -559,11 +567,11 @@ and torque on the sphere. The particle motion can then be updated in a
 molecular dynamics-like step.
 
 For the CPU implementation, a number of additional MPI communications
 molecular dynamics-like step.
 
 For the CPU implementation, a number of additional MPI communications

-are required: (1) to exchange the centre position, radius, etc of each
-particle as it moves and (2), to allow the accumulation of the force
+are required: (1) to exchange the center position, radius, etc.\ of each
+particle as it moves and (2) to allow the accumulation of the force

 and torque for each particle (the links of which may be distributed
 between up to 8 MPI tasks in three dimensions). Appropriate marshaling
 and torque for each particle (the links of which may be distributed
 between up to 8 MPI tasks in three dimensions). Appropriate marshaling

-of these data can provide an effective parallelisation for relatively
+of these data can provide an effective parallelization for relatively

 dense particle suspensions of up to  40\% solid volume fraction
 \cite{stratford2008}. As a final consideration, fluid distributions
 must be removed and replaced consistently as a given particle moves
 dense particle suspensions of up to  40\% solid volume fraction
 \cite{stratford2008}. As a final consideration, fluid distributions
 must be removed and replaced consistently as a given particle moves


@@ -573,12 +581,12 @@ With these features in mind, we can identify a number of competing
 concerns which are relevant to a GPU implementation:
 \begin{enumerate}
 \item
 concerns which are relevant to a GPU implementation:
 \begin{enumerate}
 \item

-Minimisation of host-device data transfer would argue for moving the
+Minimization of host-device data transfer would argue for moving the

 entire particle code to the GPU. However, the code in question involves
 largely conditional logic (e.g., identifying cut surface links) and
 irregular memory accesses (e.g., access to distribution elements around
 a spherical particle). These operations seem poorly suited to effective
 entire particle code to the GPU. However, the code in question involves
 largely conditional logic (e.g., identifying cut surface links) and
 irregular memory accesses (e.g., access to distribution elements around
 a spherical particle). These operations seem poorly suited to effective

-parallelisation on the GPU. As an additional complication, the sums
+parallelization on the GPU. As an additional complication, the sums

 required over the particle surface would involve potentially tricky
 and inefficient reductions in GPU memory.
 \item
 required over the particle surface would involve potentially tricky
 and inefficient reductions in GPU memory.
 \item


@@ -596,11 +604,12 @@ particle information.
 \end{enumerate}
 
 
 \end{enumerate}
 
 

-We have implemented the second option as follows. For each sub-domain,
-a list of boundary-cutting links is assembled on the CPU which includes
-the identity of the relevant element of the distribution. This list,
+We have implemented the second option as follows. For each subdomain,
+a list of boundary-cutting links is assembled on the CPU; the list includes
+the identity of the relevant element of the distribution in each case.
+This list,

 together with the particle information required to compute the correct
 together with the particle information required to compute the correct

-bounce-back term, are transferred to the GPU. The updates to the relevant
+bounce-back term, is transferred to the GPU. The updates to the relevant

 elements of the distribution can then take place on the GPU. The
 corresponding information to compute the update of the particle dynamics
 is returned
 elements of the distribution can then take place on the GPU. The
 corresponding information to compute the update of the particle dynamics
 is returned


@@ -617,14 +626,14 @@ at solid lattice points by consulting a look-up table of solid/fluid
 status. On the GPU, it is perhaps preferable to perform the collision
 stage everywhere instead of moving the look-up table to the GPU and
 introducing the associated logic.
 status. On the GPU, it is perhaps preferable to perform the collision
 stage everywhere instead of moving the look-up table to the GPU and
 introducing the associated logic.

-Ultimately, the GPU might favour other boundary methods which treat solid and
+Ultimately, the GPU might favor other boundary methods which treat solid and

 fluid on a somewhat more equal basis, for example, the immersed boundary
 fluid on a somewhat more equal basis, for example, the immersed boundary

-method \cite{ch14:immersed,ch14:immersed-lb} or smoothed profile method
+method \cite{ch14:immersed1,ch14:immersed2,ch14:immersed-lb} or smoothed profile method

 \cite{ch14:spm}.
 However, the approach adopted here  allows us to exploit
 \cite{ch14:spm}.
 However, the approach adopted here  allows us to exploit

-the GPU for the intensive fluid simulation whilst maintaining the complex
+the GPU for the intensive fluid simulation while maintaining the complex

 code required for particles on the CPU. Overheads of CPU-GPU transfer are
 code required for particles on the CPU. Overheads of CPU-GPU transfer are

-minimised by transferring only those data relevant to the hydrodynamic
+minimized by transferring only those data relevant to the hydrodynamic

 interaction implemented via bounce-back on links.
 
 % END NEW SOLID TEXT
 interaction implemented via bounce-back on links.
 
 % END NEW SOLID TEXT


@@ -768,7 +777,7 @@ keep the entire problem resident on the device.
 \section{Summary}
 \label{ch14:sec:summary}
 
 \section{Summary}
 \label{ch14:sec:summary}
 

-We have described the steps take to implement, on both single and
+We have described the steps taken to implement, on both single and

 multiple GPUs, the \textit{Ludwig} code which was originally
 designed for complex fluid problems on a conventional CPU
 architecture. We have added the necessary functionality using
 multiple GPUs, the \textit{Ludwig} code which was originally
 designed for complex fluid problems on a conventional CPU
 architecture. We have added the necessary functionality using


@@ -776,14 +785,14 @@ NVIDIA CUDA and discussed the important changes to the main
 data structures. By retaining domain decomposition and message
 passing via MPI, we have demonstrated it is possible to scale
 complex fluid problems to large numbers of GPUs in parallel.
 data structures. By retaining domain decomposition and message
 passing via MPI, we have demonstrated it is possible to scale
 complex fluid problems to large numbers of GPUs in parallel.

-By following the key design criteria of maximising parallelism
-and minimising host-device data transfers, we have confirmed
+By following the key design criteria of maximizing parallelism
+and minimizing host-device data transfers, we have confirmed that

 the mixed MPI-GPU approach is viable for scientific applications.
 For the more intricate problem of moving solid particles, we find
 it is possible to retain the more serial elements related to
 particle link operations on the CPU, while offloading only the
 parallel lattice-based operations involving the LB distribution
 the mixed MPI-GPU approach is viable for scientific applications.
 For the more intricate problem of moving solid particles, we find
 it is possible to retain the more serial elements related to
 particle link operations on the CPU, while offloading only the
 parallel lattice-based operations involving the LB distribution

-to the GPU. Again, this minimises host-device movement of data.
+to the GPU. Again, this minimizes host-device movement of data.

 
 From the software engineering viewpoint, some duplication of code
 to allow efficient implementation on both host and device is
 
 From the software engineering viewpoint, some duplication of code
 to allow efficient implementation on both host and device is


@@ -796,7 +805,7 @@ codes which also exhibit portable performance (perhaps in conjunction
 with automatic tuning approaches).
 
 So, while the challenges in designing portable and efficient scientific
 with automatic tuning approaches).
 
 So, while the challenges in designing portable and efficient scientific

-applications remain very real, this work provides some hope the large
+applications remain very real, this work provides some hope that large

 clusters of GPU machines can be used effectively for a wide range of
 complex fluid problems.
 
 clusters of GPU machines can be used effectively for a wide range of
 complex fluid problems.
 


@@ -928,7 +937,7 @@ complex fluid problems.
 
 
 % use section* for acknowledgement
 
 
 % use section* for acknowledgement

-\section*{Acknowledgments}
+\section{Acknowledgments}

 AG was supported by the CRESTA project, which has received
 funding from the European Community's Seventh Framework Programme
 (ICT-2011.9.13) under Grant Agreement no. 287703. KS was supported
 AG was supported by the CRESTA project, which has received
 funding from the European Community's Seventh Framework Programme
 (ICT-2011.9.13) under Grant Agreement no. 287703. KS was supported


@@ -958,8 +967,8 @@ by UK EPSRC under grant EV/J007404/1.
 % set second argument of \begin to the number of references
 % (used to reserve space for the reference number labels box)
 
 % set second argument of \begin to the number of references
 % (used to reserve space for the reference number labels box)
 

-
-\begin{thebibliography}{1}
+\putbib[Chapters/chapter14/biblio14]
+%\begin{thebibliography}{1}

 
 %\bibitem{IEEEhowto:kopka}
 %H.~Kopka and P.~W. Daly, \emph{A Guide to \LaTeX}, 3rd~ed.\hskip 1em plus
 
 %\bibitem{IEEEhowto:kopka}
 %H.~Kopka and P.~W. Daly, \emph{A Guide to \LaTeX}, 3rd~ed.\hskip 1em plus


@@ -967,22 +976,22 @@ by UK EPSRC under grant EV/J007404/1.
 
 
 
 
 
 

-\bibitem{succi-book}
-S. Succi, \textit{The lattice Boltzmann equation and beyond},
-Oxford University Press, Oxford, 2001.
+%% \bibitem{succi-book}
+%% S. Succi, \textit{The lattice Boltzmann equation and beyond},
+%% Oxford University Press, Oxford, 2001.

 
 %\bibitem{mpi-standard}
 %Message Passing Interface Forum, http://www.mpi-forum.org 
 
 
 %\bibitem{mpi-standard}
 %Message Passing Interface Forum, http://www.mpi-forum.org 
 

-\bibitem{desplat}
-Desplat, J.-C., I. Pagonabarraga, and P. Bladon,
-\textit{LUDWIG: A parallel lattice-Boltzmann code for complex fluids}.
-Comput. Phys. Comms., \textbf{134}, 273, 2001.
+%% \bibitem{desplat}
+%% Desplat, J.-C., I. Pagonabarraga, and P. Bladon,
+%% \textit{LUDWIG: A parallel lattice-Boltzmann code for complex fluids}.
+%% Comput. Phys. Comms., \textbf{134}, 273, 2001.

 
 

-\bibitem{aidun2010}
-C.K. Aidun and J.R. Clausen,
-\textit{Lattice Boltzmann method for complex flows},
-Ann. Rev. Fluid Mech., \textbf{42} 439--472 (2010).
+%% \bibitem{aidun2010}
+%% C.K. Aidun and J.R. Clausen,
+%% \textit{Lattice Boltzmann method for complex flows},
+%% Ann. Rev. Fluid Mech., \textbf{42} 439--472 (2010).

 
 %\bibitem{bray1994}
 %A.J. Bray,
 
 %\bibitem{bray1994}
 %A.J. Bray,


@@ -1001,10 +1010,10 @@ Ann. Rev. Fluid Mech., \textbf{42} 439--472 (2010).
 %fluid systems},
 %Phys. Rev. E, \textbf{54}, 5041--5052 (1996).
 
 %fluid systems},
 %Phys. Rev. E, \textbf{54}, 5041--5052 (1996).
 

-\bibitem{stratford2008}
-K. Stratford and I. Pagonabarraga,
-Parallel domain decomposition for lattice Boltzmann with moving particles,
-\textit{Comput. Math. with Applications} \textbf{55}, 1585 (2008).
+%% \bibitem{stratford2008}
+%% K. Stratford and I. Pagonabarraga,
+%% Parallel domain decomposition for lattice Boltzmann with moving particles,
+%% \textit{Comput. Math. with Applications} \textbf{55}, 1585 (2008).

 
 %\bibitem{xe6}
 %Cray XE6 Product Brochure, available from
 
 %\bibitem{xe6}
 %Cray XE6 Product Brochure, available from


@@ -1017,146 +1026,146 @@ Parallel domain decomposition for lattice Boltzmann with moving particles,
 %http://www.cray.com/Products/XK6/XK6.aspx (2011)
 
 
 %http://www.cray.com/Products/XK6/XK6.aspx (2011)
 
 

-\bibitem{wei2004}
-X. Wei, W. Li, K. M\"uller, and A.E. Kaufman,
-\textit{The lattice Boltzmann method for simulating gaseous phenomena},
-IEEE Transactions on Visualization and Computer Graphics,
-\textbf{10}, 164--176 (2004).
+%% \bibitem{wei2004}
+%% X. Wei, W. Li, K. M\"uller, and A.E. Kaufman,
+%% \textit{The lattice Boltzmann method for simulating gaseous phenomena},
+%% IEEE Transactions on Visualization and Computer Graphics,
+%% \textbf{10}, 164--176 (2004).

 %% Apparently first LB via GPU; a serious contribution using single fluid
 %%d3q19 single relaxation time.
 
 %% Apparently first LB via GPU; a serious contribution using single fluid
 %%d3q19 single relaxation time.
 

-\bibitem{zhu2006}
-H. Zhu, X. Liu, Y. Liu, and E. Wu,
-\textit{Simulation of miscible binary mixtures based on lattice Boltzmann method},
-Comp. Anim. Virtual Worlds, \textbf{17}, 403--410 (2006).
+%% \bibitem{zhu2006}
+%% H. Zhu, X. Liu, Y. Liu, and E. Wu,
+%% \textit{Simulation of miscible binary mixtures based on lattice Boltzmann method},
+%% Comp. Anim. Virtual Worlds, \textbf{17}, 403--410 (2006).

 %% Single relaxation (MRT mentioned) time d3q19 apparently sound although
 %%not a lot of detail. Pre-cuda so graphics code.
 
 %% Single relaxation (MRT mentioned) time d3q19 apparently sound although
 %%not a lot of detail. Pre-cuda so graphics code.
 

-\bibitem{zhao2007}
-Y. Zhao,
-\textit{Lattice Boltzmann based PDE solver on the GPU},
-Visual Comput., doi 10.1007/s00371-0070191-y (2007).
-
-\bibitem{toelke2010}
-J. T\"olke,
-Implementation of a lattice Boltzmann kernel using the compute unified
-device architecture developed by nVIDIA,
-Comput. Visual Sci. 13 29--39 (2010).
-
-\bibitem{fan2004}
-Z. Fan, F. Qiu, A. Kaufman, and S. Yoakum-Stover,
-\textit{GPU cluster for high performance computing},
-Proceedings of ACM/IEEE Supercomputing Conference, pp. 47--59,
-IEEE Computer Society Press, Pittsburgh, PA (2004).
-
-\bibitem{myre2011}
-J. Myre, S.D.C. Walsh, D. Lilja, and M.O. Saar,
-\textit{Performance analysis of single-phase, multiphase, and multicomponent
-lattice Boltzmann fluid flow simulations on GPU clusters},
-Concurrency Computat.: Pract. Exper., \textbf{23}, 332--350 (2011).
-
-\bibitem{obrecht2011}
-C. Obrecht, F. Kuznik, B. Tourancheau, and J.-J. Roux,
-\textit{Multi-GPU implementation of the lattice Boltzmann method},
-Comput. Math. with Applications,
-doi:10.1016/j.camwa.2011.02.020 (2011).
-
-\bibitem{bernaschi2010}
-M. Bernaschi, M. Fatica, S. Melchionna, S. Succi, and E. Kaxiras,
-\textit{A flexible high-performance lattice Boltzmann GPU code for the
-simulations of fluid flow in complex geometries},
-Concurrency Computat.: Pract. Exper., \textbf{22}, 1--14 (2010).
-
-\bibitem{xian2011}
-W. Xian and A. Takayuki,
-\textit{Multi-GPU performance of incompressible flow computation by
-lattice Boltzmann method on GPU cluster},
-Parallel Comput., doi:10.1016/j.parco.2011.02.007 (2011). 
-
-\bibitem{feichtinger2011}
-C. Feichtinger, J. Habich, H. K\"ostler, G. Hager, U. R\"ude, and
-G. Wellein,
-A flexible patch-based lattice Boltzmann parallelization approach
-for heterogeneous GPU-CPU clusters,
-\textit{Parallel Computing} \textbf{37} 536--549 (2011).
-
-\bibitem{wellein2006}
-G. Wellein, T. Zeiser, G Hager, and S. Donath,
-On the single processor performance of simple lattice Boltzmann kernels,
-\textit{Computers and Fluids}, \textbf{35}, 910--919 (2006).
-
-\bibitem{pohl2003}
-T. Pohl, M. Kowarschik, J. Wilke, K. Igelberger, and U. R\"ude,
-Optimization and profiling of the cache performance of parallel
-lattice Boltzmann code,
-\textit{Parallel Process Lett.} \textit{13} 549--560 (2003).
-
-\bibitem{mattila2007}
-K. Mattila, J. Hyv\"aluoma, T. Rossi, M. Aspn\"as and J. Westerholm,
-An efficient swap algorithm for the lattice Boltzmann method,
-\textit{Comput. Phys. Comms.} \textit{176} 200-210 (2007).
-
-\bibitem{wittmann2012}
-M. Wittmann, T. Zeiser, G. Hager, and G. Wellein,
-Comparison of different propagation steps for lattice Boltzmann methods,
-\textit{Comput. Math with Appl.} doi:10.1016/j.camwa.2012.05.002 (2012).
-
-\bibitem{walshsaar2012}
-S.D.C. Walsh and M.O. Saar,
-Developing extensible lattice Boltzmann simulators for general-purpose
-graphics-processing units,
-\textit{Comm. Comput. Phys.}, \textbf{13} 867--879 (2013).
-
-
-\bibitem{williams2011}
-S. Williams, L. Oliker, J. Carter, and J Shalf,
-Extracting ultra-scale lattice Boltzmann performance via
-hierarchical and distributed auto-tuning,
-\textit{Proc. SC2011}.
-
-
-\bibitem{ch14:stratford-jsp2005}
-K. Stratford, R. Adhikari, I. Pagonabarraga, and J.-C. Desplat,
-\textit{Lattice Boltzmann for Binary Fluids with Suspended Colloids},
-J. Stat. Phys. \textbf{121}, 163 (2005).
-
-\bibitem{ladd1994}
-A.J.C. Ladd,
-Numerical simulations of particle suspensions via a discretized
-Boltzmann equation. Part 1. Theoretical foundation,
-\textit{J. Fluid Mech.} \textbf{271} 285--309 (1994);
-Part II. Numerical results,
-\textit{ibid.} \textbf{271} 311--339 (1994).
+%% \bibitem{zhao2007}
+%% Y. Zhao,
+%% \textit{Lattice Boltzmann based PDE solver on the GPU},
+%% Visual Comput., doi 10.1007/s00371-0070191-y (2007).
+
+%% \bibitem{toelke2010}
+%% J. T\"olke,
+%% Implementation of a lattice Boltzmann kernel using the compute unified
+%% device architecture developed by nVIDIA,
+%% Comput. Visual Sci. 13 29--39 (2010).
+
+%% \bibitem{fan2004}
+%% Z. Fan, F. Qiu, A. Kaufman, and S. Yoakum-Stover,
+%% \textit{GPU cluster for high performance computing},
+%% Proceedings of ACM/IEEE Supercomputing Conference, pp. 47--59,
+%% IEEE Computer Society Press, Pittsburgh, PA (2004).
+
+%% \bibitem{myre2011}
+%% J. Myre, S.D.C. Walsh, D. Lilja, and M.O. Saar,
+%% \textit{Performance analysis of single-phase, multiphase, and multicomponent
+%% lattice Boltzmann fluid flow simulations on GPU clusters},
+%% Concurrency Computat.: Pract. Exper., \textbf{23}, 332--350 (2011).
+
+%% \bibitem{obrecht2011}
+%% C. Obrecht, F. Kuznik, B. Tourancheau, and J.-J. Roux,
+%% \textit{Multi-GPU implementation of the lattice Boltzmann method},
+%% Comput. Math. with Applications,
+%% doi:10.1016/j.camwa.2011.02.020 (2011).
+
+%% \bibitem{bernaschi2010}
+%% M. Bernaschi, M. Fatica, S. Melchionna, S. Succi, and E. Kaxiras,
+%% \textit{A flexible high-performance lattice Boltzmann GPU code for the
+%% simulations of fluid flow in complex geometries},
+%% Concurrency Computat.: Pract. Exper., \textbf{22}, 1--14 (2010).
+
+%% \bibitem{xian2011}
+%% W. Xian and A. Takayuki,
+%% \textit{Multi-GPU performance of incompressible flow computation by
+%% lattice Boltzmann method on GPU cluster},
+%% Parallel Comput., doi:10.1016/j.parco.2011.02.007 (2011). 
+
+%% \bibitem{feichtinger2011}
+%% C. Feichtinger, J. Habich, H. K\"ostler, G. Hager, U. R\"ude, and
+%% G. Wellein,
+%% A flexible patch-based lattice Boltzmann parallelization approach
+%% for heterogeneous GPU-CPU clusters,
+%% \textit{Parallel Computing} \textbf{37} 536--549 (2011).
+
+%% \bibitem{wellein2006}
+%% G. Wellein, T. Zeiser, G Hager, and S. Donath,
+%% On the single processor performance of simple lattice Boltzmann kernels,
+%% \textit{Computers and Fluids}, \textbf{35}, 910--919 (2006).
+
+%% \bibitem{pohl2003}
+%% T. Pohl, M. Kowarschik, J. Wilke, K. Igelberger, and U. R\"ude,
+%% Optimization and profiling of the cache performance of parallel
+%% lattice Boltzmann code,
+%% \textit{Parallel Process Lett.} \textit{13} 549--560 (2003).
+
+%% \bibitem{mattila2007}
+%% K. Mattila, J. Hyv\"aluoma, T. Rossi, M. Aspn\"as and J. Westerholm,
+%% An efficient swap algorithm for the lattice Boltzmann method,
+%% \textit{Comput. Phys. Comms.} \textit{176} 200-210 (2007).
+
+%% \bibitem{wittmann2012}
+%% M. Wittmann, T. Zeiser, G. Hager, and G. Wellein,
+%% Comparison of different propagation steps for lattice Boltzmann methods,
+%% \textit{Comput. Math with Appl.} doi:10.1016/j.camwa.2012.05.002 (2012).
+
+%% \bibitem{walshsaar2012}
+%% S.D.C. Walsh and M.O. Saar,
+%% Developing extensible lattice Boltzmann simulators for general-purpose
+%% graphics-processing units,
+%% \textit{Comm. Comput. Phys.}, \textbf{13} 867--879 (2013).
+
+
+%% \bibitem{williams2011}
+%% S. Williams, L. Oliker, J. Carter, and J Shalf,
+%% Extracting ultra-scale lattice Boltzmann performance via
+%% hierarchical and distributed auto-tuning,
+%% \textit{Proc. SC2011}.
+
+
+%% \bibitem{ch14:stratford-jsp2005}
+%% K. Stratford, R. Adhikari, I. Pagonabarraga, and J.-C. Desplat,
+%% \textit{Lattice Boltzmann for Binary Fluids with Suspended Colloids},
+%% J. Stat. Phys. \textbf{121}, 163 (2005).
+
+%% \bibitem{ladd1994}
+%% A.J.C. Ladd,
+%% Numerical simulations of particle suspensions via a discretized
+%% Boltzmann equation. Part 1. Theoretical foundation,
+%% \textit{J. Fluid Mech.} \textbf{271} 285--309 (1994);
+%% Part II. Numerical results,
+%% \textit{ibid.} \textbf{271} 311--339 (1994).

  
 
  
 

-\bibitem{nguyen2002}
-N.-Q. Nguyen and A.J.C. Ladd,
-Lubrication corrections for lattice Boltzmann simulations of particle
-suspensions,
-\textit{Phys. Rev. E} \textbf{66} 046708 (2002).
-
-\bibitem{ch14:immersed}
-C.S. Peskin,
-Flow patterns around heart valves; a numerical method,
-\textit{J. Comp. Phys.}, \textbf{10}, 252--271 (1972);
-C.S. Peskin,
-The immersed boundary method,
-\textit{Acta Nummerica} \textbf{11} 479--517 (2002).
-
-\bibitem{ch14:immersed-lb}
-Z.-G. Feng and E.E. Michaelides,
-The immersed boundary-lattice Boltzmann method for solving
-fluid-particles interaction problem,
-\textit{J. Comp. Phys.}, \textbf{195} 602--628 (2004).
-
-\bibitem{ch14:spm}
-Y. Nakayama and R. Yammamoto,
-Simulation method to resolve hydrodynamic interactions in colloidal
-dispersions,
-\textit{Phys. Rev. E}, \textbf{71} 036707 (2005).
-
-
-\end{thebibliography}
+%% \bibitem{nguyen2002}
+%% N.-Q. Nguyen and A.J.C. Ladd,
+%% Lubrication corrections for lattice Boltzmann simulations of particle
+%% suspensions,
+%% \textit{Phys. Rev. E} \textbf{66} 046708 (2002).
+
+%% \bibitem{ch14:immersed}
+%% C.S. Peskin,
+%% Flow patterns around heart valves; a numerical method,
+%% \textit{J. Comp. Phys.}, \textbf{10}, 252--271 (1972);
+%% C.S. Peskin,
+%% The immersed boundary method,
+%% \textit{Acta Nummerica} \textbf{11} 479--517 (2002).
+
+%% \bibitem{ch14:immersed-lb}
+%% Z.-G. Feng and E.E. Michaelides,
+%% The immersed boundary-lattice Boltzmann method for solving
+%% fluid-particles interaction problem,
+%% \textit{J. Comp. Phys.}, \textbf{195} 602--628 (2004).
+
+%% \bibitem{ch14:spm}
+%% Y. Nakayama and R. Yammamoto,
+%% Simulation method to resolve hydrodynamic interactions in colloidal
+%% dispersions,
+%% \textit{Phys. Rev. E}, \textbf{71} 036707 (2005).
+
+
+%\end{thebibliography}