% Week 3 Progress Report % Abdulmajed Dakkak % June 27th, 2008 Saturday/Sunday =============== Read a few papers about the LBM, HPP, and FHP CA. Lattice Boltzmann Methods ------------------------- A list differences between Lattice Boltzmann Methods (LBM) and regular numerical methods are found in Chen and Dollen's paper. The differences include: 1. The convection operator of the LBM in phase space is linear. If advection is the same as convection, then in CFD (primarily Stam's Stable Fluids method), one needs to use either a semi-Lagrangian scheme (which is unstable), or, as Stam suggests, a more numerically stable technique such as the _method of characteristics_. Advection is the term $-(\vec{u} \bullet \nabla) u$ in the Navier-Stokes equations. 2. The incompressible Navier-Stokes equations can be derived from the nearly incompressible limit of the LBM. 3. The LBM method utilizes a minimal set of velocities in phase space. HPP --- One of the first attempts to use CA to approximate NSE is called the _HPP model_ after "Hardy-Pomeau-Pazzis", the authors of a 1973 paper describing the model. The idea of the model is quite simple: * Discretize the domain to a square lattice * Only one particle found in the same vertex on the grid, * Move particles in the same previous direction * If a collision occurs, then the particles move in a direction perpendicular to their entry point. The image bellow summarizes the gist of the model \includegraphics[scale=0.1]{graphics/hpp.png} While such a model is simple and ideal for 1973 (when parallel computing was the fad), if we consider its limit, it does not approximate the NSE. This is partly because the inadequate rotational symmetry of the square lattice (having only $\pi/2$ symmetry). The HPP model is still used to simulate sand mixed with water and other rigid liquids (**need to find paper that mentioned this**). FHP --- In 1986 it was realized that the only two dimensional lattice that remotely approximates the NSE was the hexagonal lattice. This model was realized by _Frish-Hasslacher-Pomeau_ and a few months later independently by Wolfram (who references the 1986 article in CA fluids 1). The model was extended to 3 space by FHP in 1987. The FHP model is similar to the HPP model. Again, only one particle can be at any position on the grid. The update rules are a bit more complex, and described by the following figure (appears in "Lattice-gas cellular automata and lattice Boltzmann by Dieter (springer) and page 378 in NKS). \includegraphics[scale=0.1]{graphics/fhp.png} An output is selected randomly if multiple outputs can be realized. Papers Read ----------- * Somers, 1993. J.A. Somers, Direct simulation of fluid flow with cellular automata and the lattice-Boltzmann equation. Appl. Sci. Res. 51 (1993), pp. 127–133. * CHEN, S., and G.D. DOOLEN, 1998. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech, 30, 329–364. * Something on Rule 184 and its relation to traffic flow. Lectures Watched ---------------- * Wolfram's NKS [lecture][mit_lecture] at MIT in 2003. * Irving, G., Guendelman, E., Losasso, F., and Fedkiw, R. 2006. Efficient simulation of large bodies of water by coupling two and three dimensional techniques. In ACM SIGGRAPH 2006 Papers (Boston, Massachusetts, July 30 - August 03, 2006). SIGGRAPH '06. ACM, New York, NY, 805-811. DOI= http://doi.acm.org/10.1145/1179352.1141959 * Treuille, A., Lewis, A., and Popović, Z. 2006. Model reduction for real-time fluids. In ACM SIGGRAPH 2006 Papers (Boston, Massachusetts, July 30 - August 03, 2006). SIGGRAPH '06. ACM, New York, NY, 826-834. DOI= http://doi.acm.org/10.1145/1179352.1141962 [mit_lecture]: http://mitworld.mit.edu/video/149/ Monday ====== My first attempt to implement the FHP CA did not go anywhere because I had no idea how to generate the hexagonal lattice. Most of the time was thus devoted figuring out the optimal way to generate and represent the lattice. Representing the Hexagonal Lattice ---------------------------------- While the CA is simple, I wanted to look at how to optimize my grid code as early as possible, since it would be a bottle neck if the lattice is scarce. The first question was how to represent that lattice as an array. A Naive was is to think of the two dimensional array in terms of $3 \times 3$ blocks. Each of these blocks represents a single lattice; shown bellow. \includegraphics[scale=0.05]{graphics/repr.png} It should be easy, however, to see that such representation would waste 3 positions that the model does not care about. An alternative approach is to remove the middle row from the representation, and concentrate on the first and last row. We thus get the following representation 6 1 5 2 4 3 based on the following numbering scheme \includegraphics[scale=0.07]{graphics/g5.png} A $3 \times 2$ block. If we then rewrite our representation in of the FHP model using this representation we get \includegraphics[scale=0.2]{graphics/fhp2.png} It should be noted that if we want $n^2$ lattices ($n$ rows with $n$ columns) then we do not need a $n*3 \times n*2$ array. This can is evident if we stack more lattices in this representation. 6 1 51 26 462 351 531 462 426 5131 35 42 4 3 Which represents the orange region found in the following pictorial representation \includegraphics[scale=0.1]{graphics/hex.png} This $7 \times 2$ represents on three lattice rows with one lattice per row. We can generalize the size of the matrix to be 2*num_of_lattice_rows + 1 x num_of_lattice_columns + 1 > Note: Originally I was representing this as a $2 \times 3$ matrix > > 6 1 2 > 5 4 3 > > because of some complications writing this section, however, I chose to > switch to a $3 \times 2$ coding scheme. Papers Read ----------- * Ilachinski, A., (2001). Cellular Automata. City: World Scientific Publishing Company. (Chapter on Fluids) Tuesday ======= Revised the representation of the lattice. Started coding the FHP model. Code ---- * [fhp.c][tuesday-fhp](does not work) [tuesday-fhp]: code/tuesday-fhp.c Wednesday ========= After talking to Dr. Francis who explained some the history behind I was doing, I decided to go back and implement the simpler HPP model in Python, in the end, however, it did not work. One question that I and Dr. Francis had was what exactly did Wolfram try on the Connection Machine. After a few hours of research, I found some interesting information. The following is a snippet of an email exchange between me and Dr. Francis on the topic. This < http://www.stephenwolfram.com/publications/articles/ca/88-cellular/1/text.html> is about the only time wolfram mentions his work on the connection machine in the 1980s. As to what type of method he used, one has to look elsewhere from < http://www.longnow.org/views/essays/articles/ArtFeynman.php > " Cellular automata started getting attention at Thinking Machines when Stephen Wolfram, who was also spending time at the company, suggested that we should use such automata not as a model of physics, but as a practical method of simulating physical systems. Specifically, we could use one processor to simulate each cell and rules that were chosen to model something useful, like fluid dynamics. For two-dimensional problems there was a neat solution to the anisotropy problem since [Frisch, Hasslacher, Pomeau] had shown that a hexagonal lattice with a simple set of rules produced isotropic behavior at the macro scale. Wolfram used this method on the Connection Machine to produce a beautiful movie of a turbulent fluid flow in two dimensions. Watching the movie got all of us, especially Feynman, excited about physical simulation. We all started planning additions to the hardware, such as support of floating point arithmetic that would make it possible for us to perform and display a variety of simulations in real time." The above quote comes partly from a book called "Feynman and Computation," and that is the only reference to CA and fluids (or wolfram) in that book. A more (auto)bibliographical account can be found in NKS < http://www.wolframscience.com/nksonline/page-880b-text > "I had always thought that cellular automata could be a way to get at foundational questions in thermodynamics and hydrodynamics. And in mid-1985, partly in an attempt to find uses for the Connection Machine, I devised a practical scheme for doing fluid mechanics with cellular automata (see page 378 ). Then over the course of that winter and the following spring I analyzed the scheme and worked out its correspondence to the traditional continuum approach. By 1986, however, I felt that I had answered at least the first round of obvious questions about cellular automata, and it increasingly seemed that it would not be easier to go further with the computational tools available. In June 1986 I organized one last conference on cellular automata--then in August 1986 essentially left the field to begin the development of * Mathematica*." He thus used the FHP model that I am trying to get to work right now. Also, the 4D fluid simulation is called FHCP < http://arxiv.org/abs/chao-dyn/9508001v1 > It should also be mentioned that Wolfram patented not only the hexagonal method for solving fluid dynamics, but the entire lattice gas method as a whole. See patent number 4809202 < http://www.google.com/patents?id=FF4WAAAAEBAJ&dq=U.S.+Patent+Number+4,809,202> Code ---- * [hpp.py][hpp.py](does not work) [hpp.py]: code/hpp.py Thursday ======== I realized that my interpretation of the model was incorrect the entire time. After I gave a presentation, received some helpful tips, and learned more, I decided to start from scratch. After a few hours, I was able to correctly implement the HPP model. The implementation, along with the algorithm, is detailed in [HPP document][hpp.doc] posted on the website. Code ---- * [hpp.c][hpp.c] [hpp.c]: code/hpp.c [hpp.doc]: ../lattice-gas/lattice-gas.html Friday ====== Looked at how Jared Schaber, and wrote [notes][dbz]. Talked to Jonathan Manton who described how to implement the FHP model (I wish I had that talk before, since I figured out much of what he said the day before). Jonathan also explained the concepts behind CUDA, compiler optimization, and a slew of other things. He also gave me tips on how to optimize my program, while stressing that clarity is the most important objective at the beginning. [dbz]: ../dbz/dbz.html Questions ========= In _Cellular Automata Super Computers for Fluid Dynamics Modeling_ (by Margolus) the author writes the following about implementing CA in conventional architecture (not CAM): > two dimensional planes are processed serially (with a substantial amount > of pipelining); a third dimension is achieved by staking planes and > operating on them in parallel. The question is: how can this be achieved without running into race conditions.