Archive for the ‘Programming’ Category

I realise I’ve been neglecting the blog a bit the last few months, with only a few posts that link to other things, but that’s because I’ve been solving a nasty problem for my research and I think it’s now time to update what’s going on with a ‘real’ post.

For the past year or so, we’ve been working on modelling a galaxy using a set of models developed by my supervisor, Larry Widrow, and his collaborator, John Dubinski, who is at Toronto. The models assume a standard exponential disc for spiral galaxies, but the curious thing is that, because discs appear to contain a luminous bulge at the centre (which is generally considered to be a dynamically distinct entity) there is degeneracy in how surface brightness profiles can be broken down: typically, we break them down into a part due to the bulge (which is often assumed to be a de Vaucouleurs-type r^1/4 profile) and a part due to the disc, which is assumed to be exponential. In the centre, the bulge and the disc contribute to the total light seen in a galaxy. There is some theory behind assuming that discs form exponential density distribution, but in principle other functional forms are allowed for disks. One of these is called a Kormendy disc (named after John Kormendy who proposed this in 1977).

A Kormendy disc is basically an exponential disc with a hole in the centre. So the centre of the disc emits no light, because there are no stars in the centre. This sounds weird, but there is nothing preventing us from attributing all the central light to the bulge, and letting the disc take over as the bulge falls off – we can find fits using the Kormendy functional form as well as the standard exponential. In fact, a siginificant minority of galaxies appear to be better fit by a bulge + Kormendy disc than by a bulge + exponential disc.

The presence of a hole is highly counterintuitive, and you may wonder if ‘seeing’ a hole is merely a measurement artifact or a real measurement due to some other effect. Neither of these appears to be the case; measurement artifacts have been ruled out, and effects such as central dust (which could obscure disc light) are unlikely to be the cause. But, there are good reasons to think it reflects a real physical effect; many of the galaxies that display holes are known to be barred (i.e., they have a central bar of stars out of which typically emanate spiral arms), and bars could be responsible for the effect, but even some that don’t have bars are well fit by holes.

So what has this got to do with my research? Well, the galaxy we’re looking at might in fact be fit better by a central hole than by the exponential disc in our models. Originally we accounted for this by pretending that the real mass distribution was exponential while the light distribution had a Kormendy hole. This is reasonable, because the galaxy appears to have some dust at the centre, although, as I mentioned above, holes are not likely due to dust. So we decided to fundamentally change the structure of the models to generate a Kormendy disc instead, and that’s what I’ve been doing for the past two months… and it’s a pain, because the code is written in Fortran (the F-word again) 77, and Fortran 77 is just plain annoying.

Fortunately, the code is structured in such a way that the radial profile can be changed separately without altering anything else, so it was in principle just a matter of changing the density wherever it was called in the code. Unfortunately, the code itself has been changed several times over the past 15 years and the syntax and variable names used across different subroutines are inconsistent. Moreover, just changing the density profile at first did not work, partly because there is a mathematical (though not a physical) singularity at the centre; to rectify this, I just redefined the central density (and the first two derivatives) to be zero. Debugging was annoying because outputting important quantities in subroutines ran into fortran’s absurd I/O recursion error, which really seems to be an annoying bug that has just never been fixed (as evidenced by the fact that it would pop up sometimes and not at other times when I had made no changes to the code). Finally, the code has a bothersome habit of giving different results for the same input when compiled on different machines… I’ve quietly decided to ignore that since it seems to work on my department machine (but on my home desktop or my laptop). Regardless, after two months of fighting through oceans of NANs and INFs, outputting dozens of different quantities at all points of the program, and banging my head against the wall (not literally, but it wasn’t going to stay not literal for long), I tried an arbitrary but well-behaved radial profile and found that it worked perfectly, which indicated that the main problem was likely a coding bug or something having to do with the centre of the Kormendy disc. I finally traced the problem to the calculation of certain spline quantities at the centre… where the density was zero. Instead of fighting my way through tedious, possibly buggy numerical subroutines, I changed the density (and the first two derivatives) to a very small nonzero quantity… and presto! We have self-consistent Kormendy discs.

So now we have real Kormendy discs, and we can model galaxies with either a Kormendy profile or an exponential profile, and the sun is shining, and life is grand (well, sort of, but that’s another story :P).


Read Full Post »

I haven’t posted in along time, and the reason is that I’ve been trying to get a production run going for the research on HPCVL (big computing cluster) and it’s always trying to port your code to another machine with a different compiler. So here’s what happened.

Three weeks ago we decided to get a big MCMC run going on HPCVL, which required me to compile code on a Solaris machine using the Sun compiler. My programs use automake, which means I had lots of fun figuring out how to configure and install them into my home directory on HPCVL (aha! ‘prefix’ flag required!), and had my memory refreshed on changing environment variables/linking libraries many times. So after running configure, ‘make’ choked in several different places, not all of which seemed logical:

  • Problem: CC can’t figure out what ‘sqrt’, ‘sin’, etc mean. Fine, I need a #include <cmath> statement in a header file. Makes sense, but it’s still not clear to me why it compiles on my computer at Stirling. Possibly something to do with gcc 3.2, because it won’t compile without the aforementioned directive on my home computer with gcc 4. In any case, that’s an easy fix.
  • Problem: There’s some kind of scope problem involving trying to define a member of a class definition and putting an extraneous scope operator. I think. I still have no idea what the error message meant, but removing the extraneous scope operator made it go away. Good enough.
  • Problem: I can’t compile with the cosmology header file as it is because the compiler won’t accept initializing const static int (the Hubble constant) members in the ‘protected:’ section of a class definition. Since this pushes the limits of my knowledge of C++, I try to initialize them in random places elsewhere in the class file, notably under where it says ‘public:’, using a constructor, keeping the original declarations in place. This works, in the sense that the original compiler error message is gone and I now get a different one. This time the error occurred in the linking stage, with one of those horrific error message that looks like a long string of random letters and numbers that says absolutely nothing helpful, as opposed to the cryptic compilation error message I got earlier that at least told me what line the error was on. It’s complaining about doubly declared variables or something like that, and I notice that the end of the long string of meaningless characters, the names of the rogue const static ints are appended. This tells me that the problem remains with the initialization of these members, but nothing else. After multiple ham-handed attempts at working around the problems using all sorts of syntactical gymnastics, I finally declared them as global variables. Problem solved.
  • After finally getting everything to compile correctly, I tried a few things. Notice problems with the calculation of the rotation curve and density profiles of sample systems. The rotation curve problem is easy, because somehow the wrong calculation for the tangential velocity was included in the original source (long since fixed on my computer here). The density is slightly more difficult, because, for no clear reason, it calculates the density profile correctly, and then spits out the wrong value. More specifically, it gets the correct density values while in the for loop that loops through each bin (which I checked by outputting the density in every bin within the for loop), but outputs the wrong values immediately on exiting the for loop (which i checked by outputting the density values in a separate for loop). I still haven’t figured that one out. I could just output the density within the original for loop but I would like to know why it chokes after leaving the original loop.

So the analysis programs (more or less) work, or seem to. Compiling the N-body code itself goes smoothly and seems to work fine as well. Now to compile the galaxy-building programs.

  • There is no excuse in this day and age to be forced to submit to a 72-character per line limit (really 66 character because you can’t use the first six columns except for special cases) like fortran 77 requires. So I have to change a flag in the makefile. Fine. Except the Sun compiler’s flag can only handle up to 132 characters per line. I shouldn’t be making lines that long; fair enough.
  • After compilation, I compare two of the same models generated on HPCVL and on my computer here at Stirling. Slight discrepancy in the central potential. Hmm. Larry suggests doing some detective work to uncover the cause. After multiple write(*,*) statements, I discover that tiny roundoff errors on each machine can lead to big differences when they’re used as arguments to logs and are really close to zero. Shouldn’t be too big a deal.
  • I then try building a galaxy. Everything goes (more or less) well. Analyse the model, density still doesn’t work. Phooey. But I’ll manage. Everything else seems okay. So after getting a trial MCMC run started, I did a little more inspecting of the resulting models. I look at one model, generate the disk, bulge and halo, and work out the rotation curve. I notice a slight anomaly: the halo rotation is not being calculated correctly: it’s too low. Output the mass as a function of radius. Only then do I discover, to my horror, that the halo has giant hole in the centre. Everything else is correct – in particular, the total mass and tidal radius are correct. It’s not distributing the mass correctly. There’s just nothing that ends up in the centre. Nada. I try changing one halo parameter. This time, everything is fine. I try changing the whole parameter set. Everything is fine. Okay. So what’s going wrong? Larry thinks it’s a bug in the interpolation routines that would take a month to find. Since we can’t control where a MCMC chain goes in parameter space, there’s no telling if it’ll find one of these regions where the bug will manifest itself. What a mess. I’m toying with the idea of spending my Christmas holidays rewriting all this code in C. But for now we’ll do the MCMC run on my regular computer and use HPCVL for the simulations.

Now I’m trying to figure out how to account for the asymmetric drift in the galaxy we’re trying to model. Which, you may notice, is a problem that has actual physics involved. It’s nice to have one of those once in a while – cuz, ya know, I am technically an astrophysicist and not a programmer.

Read Full Post »