Archive for the ‘regular update’ Category

OpenMP and more C++

Monday, 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

Thursday, 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

Friday, 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

Wednesday, 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

Monday, 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

Finished with Binney & Tremaine

Friday, July 24th, 2009
I have finished reading Binney and Tremaine's  Galactic Dynamics book. I will read the updated material from the second edition of the text next week. I've compiled a rather long list of reading material here.  The list includes Annual Reviews, papers by Gaby, Bruno, and Phil Hopkins et al. as well as papers on n-body simulation, SPH, and the GADGET code. I am also compiling a list of useful links here. The reading list and various programming projects should occupy my time for the next few months.

Stellar collisions, binaries, core collapse

Wednesday, July 22nd, 2009
I finished reading an introduction to stellar collisions, the formation and effects of binaries, and core collapse.  The topics are interrelated and must be included in a good galaxy simulation.  The problem is that properly modeling the interactions between stars while simultaneously modeling the global galaxy evolution is very computationally expensive.  However it is necessary if the formation and growth of black holes in the nuclei of galaxies and globular clusters is to be consistently modeled.

Kinetic theory & GADGET-2

Tuesday, July 21st, 2009
After a bit of tinkering I was able to compile "Hello World!" programs in Fortran and C++.  Now I need to figure out which one would be better to use for the upcoming programs and I need to learn more about Python. I discovered that Phil Hopkins is a year younger than me and he is now a post-doc at Berkeley. He and many of the other people in this field use the same code called GADGET-2: www.mpa-garching.mpg.de/gadget/. I began reading about the application of kinetic theory to gravitationally bound particles today. I already had an introduction to the subject when I was reading for my proposal assignment, so it was nice to see the theory presented more thoroughly. The most interesting feature, theoretically speaking, is the negative heat capacity caused by gravitation. Removing energy from a system with negative heat capacity makes the system's temperature increase! As with much of the rest of the book, this section contains some silly attempts to derive analytic approximations for things which can only be studied properly through simulation.  Analytic arguments in systems as complex and asymmetric as clusters and galaxies seem to only be useful for building a vocabulary with which to describe phenomena.

Dynamical friction

Sunday, July 19th, 2009
I've begun reading about collisions of galaxies.  The first topic in the chapter was dynamical friction.  There is a very clever derivation in section 7.1 of Binney and Tremaine.  If I taught a class on galactic dynamics, I would definitely include this derivation.

Spiral structure summary and Warps

Saturday, July 18th, 2009
I finished reading the introduction to spiral structure.  This entry is a very brief overview. Observations: The vast majority of spiral arms are observed to be trailing the rotation of the stars (i.e. the relationship between the pattern and rotation direction is the same as in the case of hurricanes  and cyclones on Earth).  Most spiral galaxies containing prominent arms only have two spiral arms, but some have more. The number density of all stellar types is higher in the arms, but the distribution of young blue stars is much narrower than the older stars. Anti-spiral theorem and consequences: If spiral galaxies are in a steady state, the cause of the spiral pattern cannot be solely time-reversible physics. Otherwise, trailing patterns and leading patterns would be equally likely.   If time-reversible processes are the cause of the arms, then spiral galaxies are not in a steady state. Lin-Shu hypothesis: Spiral arms are quasi-static density waves.  The density wave is caused by a combination of orbital resonances in a rotating potential and instabilities which arise in thin (quasi 2D) potentials.  The Swing Amplifier helps to enhance the pattern. Stars and clouds of gas drift into and out of the arms but the overall pattern is relatively fixed.

Swing amplification: A disturbance (density perturbation) which is arranged in a leading, rather than trailing, pattern unwinds, then propagates through the central part of the galaxy and emerges as an amplified disturbance with a trailing pattern.  This was not explained very well in the text.

Some alternatives to Lin-Shu: Spiral arms could be transient and in a state of constant renewal.  There is a "detonation wave" theory that states that star formation induces more star formation and thus the spiral arms are like the front line of a forest fire.  Some spirals could be induced from the recent interaction between galaxies or from the perturbations due to a satellite galaxy.  The spirals could also be caused by the rotating oval or bar-shaped potential at the center of the galaxy.
Warps: Disk galaxies are observed to be "warped" on their outer edges.  If a cylindrical coordinate system is defined at the center of the galaxy with z perpendicular to the disk, the matter in the outer regions is seen to be displaced in the +z direction in some regions and in the -z direction in other regions.  The cause of this warping is not understood.