C++ in depth (part 3)

August 19th, 2009
My whirlwind tour of the C++ language and the standard C++ library is complete.  Now I just need to practice using it.  Since I wrote the last entry, I have read about...
  • exception handling
  • the C++ I/O system
    • streams
    • advanced file formatting
    • binary I/O
    • I learned that "files" in C++ are thought of the same way as in UNIX, meaning that all hardware devices are considered to be "files". (The details of I/O could easily fill a book, so I've only had an overview)
  • the run-time type ID (RTTI)
  • casting operators
  • namespaces
  • pointers to functions (allows a function to be passed to another function as an argument)
  • the asm keyword (allows assembly code to be incorporated into C++)
  • linkage specification (allows a program to call functions written in another language)
  • an overview of the standard template library (the full details of the STL could fill a book)

C++ in depth (part 2)

August 17th, 2009
Today I learned how C++ implements object oriented programming (OOP).   OOP allows for the implementation of class libraries, which is why I needed to learn more about C++ in the first place (some of the visualization libraries contain class libraries). The topics that I learned about today include
  • operator overloading
  • inheritance
    • public, protected, and private access rights to member data and functions
    • methods of properly dealing with constuctors and destructors of base classes and derived classes
    • virtual base classes
    • virtual functions
  • Polymorphism
  • Templates
    • Generic functions (very useful!)
    • generic classes
I'll soon forget the details of the syntax and some of the important caveats, but I've learned the general ideas and the most basic syntax, which is the important part.  I can always refer to the details as needed.  Over the next two days, I'll learn about exception handling, the Standard Template Library (STL), namespaces, advanced I/O, and assorted "advanced" topics.

C++ in depth

August 16th, 2009
In the course of evaluating the various visualization options, I discovered that I needed to learn a few more things about C and C++ in order to implement a visualization library.  Rather than continuing to learn more features of the language incrementally, I decided that it would be better to just read a description of the entire C++ language!  On Friday I grabbed one of the C++ books from my bookcase and started reading.  The book is "C++ from the ground up" by Herbert Schildt.  I've already learned a lot of things that I previously didn't know existed in the language.  I'm half-way through the book now and I plan to finish reading it Wednesday night.  Once I'm finished with this, I should be able to read C and C++ code much more easily and I'll be able to write code more efficiently because I will know all of the major options available to me. The book will serve as a reference for the detailed syntax until I internalize everything.  Luckily, the language is very intuitive thus far. I'm still hoping to begin implementing a visualization system before I leave for VA. I also need to read a few more items on my reading list to keep a good pace going.  I'll try to read while I'm in VA. A brief summary of the things that I read about (in exhaustive detail) today:
  • the const, volatile, extern, and static type qualifiers
  • the register keyword (which I need to experiment with)
  • enumerations
  • bitwise operators |, &, ^, ~, >>, <<
  • the ? operator
  • the , operator
  • a review of structures
  • bit-fields
  • unions
  • classes and objects
  • constructors  and destructors
  • friend functions
  • the copy constructor
  • the this keyword
I also found a research group in Computer Science and Engineering that may be useful some day: The Riverside Graphics Lab.  They do visualization research and of course they use Unix / Linux.

Bruno's paper and CUDA

August 13th, 2009
I took a break from writing code today and finished reading Bruno's paper on continuous stellar mass loss in N-body simulations. His code used the particle mesh method and included some heuristics for mass loss and star formation. The "sticky particles" method was used in modeling collisions of gas clouds. There seemed to be very little actual physics in the simulation. I'm still evaluating different possibilities for visualizing my N-body simulation output. The Visualization Toolkit (VTK) is the newest possibility I've found. It's the foundation on which many other visualization systems are built. Hopefully I'll be able to make a decision by the end of next week. My interest in CUDA is growing once again.  Once I become more experienced with programming in C / C++, it should be easier to learn CUDA.  I found an interesting little PDF on using CUDA for N-body simulations: here.  I have the necessary hardware and software to develop CUDA programs--I just need to learn a bit more.

Issues resolved!

August 12th, 2009
My difficulties with storing large arrays and parallelization have been resolved! I've also learned a lot more about the GNU Compiler Collection (GCC) and how to pass flags to the compiler using the Make system (this is extremely useful). In the process of reading the GCC documentation, I skimmed through instructions on how to vectorize code using SIMD extensions like 3D-Now!, SSE, and SSE2. That might come in handy some day, so it's good that I know where to look if I want to learn more. I've also switched from KDE's IDE to GNOME's IDE (Anjuta DevStudio) because it turns out to be much more straightforward to use, and thus it makes me more productive.   It appears that creating GUI programs wouldn't be very difficult using Anjuta and Glade along with GTK or wxWidgets.  Evidently programs built using these tools can easily be ported to Mac and Windows. The images below were created on my quad core home workstation.  Figure 1 shows 10 orbits calculated with the newly parallelized code.  The orbits are color-coded by thread number.   Each thread was executed on a different CPU core.  Figure 2 shows the CPU and memory utilization for a run consisting of 10,000 orbits with 25,000 steps per orbit.  The memory usage increases because the position of each particle and corresponding thread number are saved after every step (I was intentionally trying to stress the memory system).  Note that all 4 cores reach 100% CPU utilization and the system runs out of main memory and has to start using the swap partition toward the end of the run.

Figure 1

Figure 2

Figure 1 Frigure 2
Here is an animation of the orbits rotating so that the 3D effect is more obvious:

OpenMP and more C++

August 10th, 2009
The leap-frog integrator has been implemented and it yields exactly the same behavior as discussed in the GADGET-2 paper; It's more reliable than second order Runge-Kutta even though they are both second order methods.  I've modified the code so that it initializes particles according to a simple distribution function. The images below are IFrIT plots of  3, 10, and 100 randomly oriented orbits in a simple \Phi \propto r^{-1} potential, as well as a plot of 25,000 initial positions.
3-orbits10-orbits100-orbits25000-points
The main focus of the last few days has been C++ programming and building the individual pieces that will be needed for the N-body code later this month.  I've been able to use OpenMP to parallelize the main for-loop.  It's surprisingly easy.  There is  noticeable speed increase on my dual processor system ( Note: I have not yet benchmarked the program to measure the actual speed-up factor).  The CPU utilization now reaches 100% rather than the previous 50%.  I'm currently working out the bugs with the shared memory system.  Writing a code for parallel execution requires me to modify my way of thinking a bit.  I am also trying to figure out how to handle large amounts of memory in C++.  Once these issues are resolved, I will try to implement a data visualization library so the program can plot output automatically.  This will be important for the next project because I want to be able to make animations of the N-body simulations. In terms of reading, I just started reading Bruno's  paper dealing with continuous stellar mass-loss in N-body models.

Beginning the projects

August 6th, 2009
I finished reading the description of the GADGET-2 code and I've read a few more annual reviews. I've been spending a lot of time becoming more familiar with C++ and the compiler that I'm using. I've also been looking for different plotting options. I've managed to write a very simple code for tracing the orbit of a particle in a central force field, so technically I'm finished with the first of the three projects, but I want to spend tomorrow and possibly Saturday making improvements to the code. I want to implement a leap-frog integrator and possibly use pointers instead of the messy syntax I currently have. I would also like to have the program automatically plot the data. Currently the output is a CSV file with coordinates. Here's a plot of the output after about 4.5 orbits: Screenshot

Winds, Rings, Shells, GADGET-2

July 31st, 2009
I've finished reading an annual review on galactic winds as well as one on rings and shells / ripples.  Galactic winds alone are a fairly complicated topic, and very important in galactic evolution.  Evidently the publicly released version of the GADGET-2 code does not incorporate the processes necessary to simulate the wind.  The publicly released version is actually missing a lot of the more interesting components which allow feedback.  I'm guessing that the public version is simplified so that only preferred users can easily publish major new results.  The full code may also still be buggy and difficult to run--especially since it uses MPI rather than OpenMP. I've started reading the description of the GADGET-2 code.  The most significant things to note are: (1) The multipole moments used in the TreeSPH gravitational potential calculations are just monopoles rather than quadrupoles or higher. (2) The code now includes the option of using a hybrid Tree-particle mesh (TreePM) technique which uses a mesh and Fourier techniques in the far field while using TreeSPH for local calculations. (3) The SPH formulation is based on entropy, so it explicitly conserves entropy with fully-adaptive smoothing length.  Anomolous entropy changes were evidently an issue with the original GADGET code. (4) Cosmological expansion is included explicitly. (5) As far as I can tell, the gravitational model is essentially Newtonian with infinite propagation speed.  Implementing a consistent gravitational model would probably greatly complicate the code and make it run much more slowly. My goal for the next week is to complete the first program described in my first post. It will be written in C++.  I would also like to be able to create a simple plot of 3D data with IFrIT.

SPH & Hierarchical tree method

July 29th, 2009
I now have a fairly good understanding of the mathematics behind SPH.  It is fundamentally a fluid technique.  The fact that it uses "particles" to represent the fluid does not make it an N-body method. The "particles" are merely tracers--a way to store the data describing the fluid properties at various points in space.  There are two levels of smoothing involved when the method is used for gravitationally interacting fluids.  First there is the ordinary smoothing, characterized by the smoothing length h. This parameter determines the size of the smoothing kernel--essentially the spatial extent of each particle.  The second type of smoothing is called "softening".  It is characterized by the softening length \varepsilon and its purpose it to avoid singularities or near-singularities in the gravitational potential.  Thus SPH-based galaxy simulations codes cannot resolve fine details such as individual stars in galaxies unless the smoothing and softening parameters are set to zero, in which case the technique is really no longer SPH.  The method also has great difficulty handling shocks properly with a finite number of particles. The hierarchical tree method is a way of improving the efficiency of N-body simulations by effectively smoothing the distant particle distribution.  In the TreeSPH code, which was the precursor to GADGET, an oct-tree is created which allows the particles to be sorted spatially.  Multipole moments of blocks of space are calculated in the process of contructing the oct-tree.  for each particle, the gravitational force is calculated by descending down the oct-tree toward that particle and collecting data along the way.  The nearest-neighbor particles are treated as pure particles and gravitational forces are computed directly as in an N-body code.  The field of the next-nearest neighbors is calculated form the multipole moments of small "cells" of space.  As the distance from the particle increases, the cells grow larger.   This approximation greatly improves the efficiency of the force calculation. Other than learning about SPH and TreeSPH, I've been tinkering with C++, just getting some practice and becoming more familiar with the compiler.  I'll read a description of the GADGET-2 code.  I found out that it is written in C and it uses MPI, rather than OpenMP. I've installed IFrIT, which looks like a good tool for visualizing output from the simple programs that I plan to write this summer. sites.google.com/site/ifrithome/Home

Sagittarius A* , Jets, visualization

July 27th, 2009
It turns out that reading Binney and Tremaine has paid off. It's much easier to read research papers now that I know a large portion of the vocabulary. Today I read an annual review about Sagittarius A* at the center of the Milky Way. It was a pretty rich source of information. The most interesting part was the discussion of the processes which might cause the jets associated with AGN (the ideas which had been proposed as of 2001, anyway) I read a great paper by John Dubinski called "Visualizing astrophysical N-body systems". The paper gives an overview of the history of N-body simulations and visualizations, describes visualization fundamentals, describes the state of the art techniques,  and discusses what will be done next.  I'm interested in getting a copy of Dubinski's MYRIAD visualization software library, since it is apparently the best and it's free.  It may be more practical to use one of the simpler codes mentioned in the paper first though.  Unlike the space weather world that I was working in, the galactic simulation world seems to be almost entirely open-source! Dubinski's paper led me to one of the first N-body animations from the 1970s: mindcine.com/toomre/ links to all papers discussed in this research log can be found in the reading list