15-859E Programming Assignment 1:
N-Body Simulation

I strongly suggest you try to compile and run the starter code by Thursday 24 Sept. to help catch compiler and operating system difficulties early.

Compiling and linking

The starter code is in pubdir/src/nbody, where pubdir = /afs/cs/project/classes-ph/859E/pub . Take a look at it here. You'll need to make your own copy of these files elsewhere before you start modifying them. C++ code and headers use suffixes .cc and .hh, while C uses .c and .h . As of 9/22, I've tested it on SGI's only, but on other UNIX machines, changes to the Makefile will hopefully suffice to get it running. You can find public SGI machines in Wean 5205. On an SGI, run "make -f Makefile.sgi" to compile and link the program. On a Sun, use "make -f Makefile.sun". This should create the executable "nbody".

If you get an error message such as "don't know how to make /usr/include/X11/Xlib.h" then Xwindows (or some other library) is not in the same location on your machine as on the cluster machines for which the Makefiles were configured. First try "make -f Makefile.systype depend" to remake the dependency list at the end of the Makefile, and if that doesn't work, search around for the missing include files and change the Makefile appropriately. On a Sun, if you get "warning: file /usr/local/lib/libstdc++.so: section .stabstr: malformed...", don't worry about it; it's just a warning.

See pubdir/www/software.html for instructions on how to set up environment variables so you can link shared object (.so) files. Do this now.

Running the starter code

• left-mouse picks and moves a particle. You have to click fairly close to a particle.
• shift-left-mouse deletes a particle
• middle-mouse creates a particle
• right-mouse translates the world
• shift-right-mouse (moving left & right) zooms the world
• mass slider sets the mass (and the radius) of the currently selected particle (the one most recently picked)
• log(dt) slider sets the time step (currently unused). the number displayed is the log base 10 of delta_time (this seemed like a reasonable way to control small numbers)
• run button is intended to run the gravitational n-body simulation, but that code hasn't been written yet, so it's currently a no-op. Note that you can edit the world (create, move, & delete particles while the simulation is running.
• if you type a key, the code will tell you what key you pressed. This could be useful to you later.

To facilitate testing and interchange of data between students, I added a parser for a simple "particle2" file format to nbody.cc on 2 Oct. The file format is described in the example file pubdir/src/nbody/orbit.p2. I also added the gravitational constant G to the World structure.

The user interface was implemented using the Xforms library (see pubdir/www/software.html for documentation and downloads for this and other libraries mentioned here). The layout for the user interface was set using the fdesign program described in the documentation. If you run "fdesign ui_form.fd" you can rearrange the user interface. Hitting "save" while in that program rewrites ui_form.c and ui_form.h .

Graphics is done using the OpenGL library. We're only doing 2-D graphics, so we could have done it all with Xlib, I suppose, but OpenGL is much cleaner than Xlib. Double buffering (definition) is used to make the redisplay smooth.

Event handling (mouse and keyboard input) is done using Xlib, the low level library for Xwindows.

The main program, nbody.cc, is written in C++, but stays away from the advanced features of the language. The main features of C++ used that are not present in C are:

• classes - for modularity (encapsulation)
• operator overloading - so we can do arithmetic on 2-D vectors and points very concisely
• more stringent type checking during compilation

• Implement a hierarchical n-body algorithm in two or more dimensions that runs in O(NlogN) or O(N) time. Implementing just the brute force O(N^2) algorithm does not suffice.
• Get your program running on large systems (thousands of particles).
• Also implement the brute force O(N^2) algorithm for reference -- probably do this first.
• Check that your program is working right (that simulations are reasonably accurate). A good check is that the potential energy + kinetic energy is constant.
• Test it on the file circ.p2, which gives the initial conditions for a planet in a (nearly) circular orbit around a sun. If your planet does not make an approximately circular orbit, then it's probably wrong. Small time steps should make the orbit more circular, while large time steps will probably cause the planet to spiral out. Or test it on output of the programs orbit.cc or randsq.cc. The former can create an spinning "galaxy" of n stars. These files are in the pubdir/src/nbody directory. (To compile these, run "make -f Makefile.systype orbit", etc.)
• Finish your code by the date listed at top. To submit your program, copy your source code, executable, and test files (if any) to students/yourlastname/nbody. Include a brief README file explaining what you did, what system your executable runs on, what works / doesn't work, and how you tested. Also report the size of the biggest set of particles you've run a simulation on, the time per iteration for it, and the system you ran this on (e.g. 900 MHz Pentium 4). The date and time on these files matters (it's how I check for lateness) so if you want to add additional notes/corrections later, do it in a subdirectory, and talk to me or email me.
• Demo your program to me and the other students during the class period on Oct. 6, preferably on the machines in Wean 520x.
• When grading, I'll look for correctness first, then speed per iteration, then smart time stepping, then user interface and graphics. When I evaluate speed, constant factors don't matter as much to me as asymptotic complexity (O(N^2) is bad, O(NlogN) is good).

Tips

• I suggest you start with my starter code, but you don't have to.
• The recommended algorithm is Barnes-Hut, as described in Nature, vol. 324, 1986. You'll probably want to put their parameter theta on a slider.
• If you'd rather implement some other hierarchical n-body method, such as the fast multipole method, be my guest, but this will probably be harder. The same for 3-D instead of 2-D.
• I recommend you do the simulations in 2-D, since it's much easier to visualize, and use the following formulas for pairwise potential and force between particles i and j:

  potential energy = phi = G*mi*mj*log(rij)
  force = -gradient(phi) = -G*mi*mj*(Xi-Xj)/rij^2

  where mi and Xi are the mass and 2-D position of particle i, rij=||Xi-Xj|| is the distance between particles i and j, and G is the gravitational constant. This formula was discussed in the Greengard-Rokhlin paper (J. Comput. Phys. 73, 325-348, 1987). Note that this is not Newton's law, but it tends to give more interesting behavior in 2-D, and it obeys Laplace's equation in 2-D. Whatever potential you choose, make sure that your potential and force are consistent (force = -gradient(phi)). The total potential energy is just the sum of all pairwise potentials.
• Regardless of which potential you use, the kinetic energy of a particle is .5*mi*||Vi||^2 where Vi is the velocity of particle i. The total kinetic energy is the sum of particle KE's.
• To cause a "planet" to orbit the "sun" in a circle, make the ratio of masses very large and make the planet's initial velocity be tangential to the circle with speed sqrt(Gm), where m is the mass of the sun. This is a solution to the 2-body problem for log(r) potential. The program orbit.cc and file circ.p2 mentioned earlier follow this prescription. If you use 1/r potential, the tangential speed should be sqrt(Gm/r). You can easily verify these relationships by assuming motionless sun, circular motion, and plugging into the differential equation to solve for circumferential speed.
• I've been using non-physical units and G=1. If you want to work in physical units (meters, grams, etc) then you can use Newton's law and set G to the appropriate value.
• I've used Euler's method for solving the differential equation. You could try Runge-Kutta or other methods, if you like. Let us know how well they work. See a numerical analysis text or the book Numerical Recipes.
• When two particles come very close together, small time steps are necessary to get accurate simulations. Otherwise, near encounters sometimes lead to particles shooting out very fast.
• Implementing adaptive time steps is optional. See the book Numerical Recipes for generic methods.
• If you have time, and you want to improve the graphics, you might consider some options: only redraw every k time steps, draw antialiased points or circles, draw trails behind the particles, draw velocity or acceleration vectors, draw the hierarchy in some way so we can visualize how the algorithm is using it.

15-859, Hierarchical Methods for Simulation
Paul Heckbert