I was reading about the recent black hole mergers simulations performed by the people at Goddard, more thoroughly described here and in this forthcoming article. These are, undoubtedly, beautiful results, and a testament to the complexity of Einstein’s equations when it comes to obtain realistic results: according to the reports above, thousands of lines of code (plus an impressive array of supercomputers) were needed to obtain them. An impressive achievement, but still there’s something in it that makes me uneasy: if i’ve read it right, those thousands of lines of code are actually lines of Fortran (or, in the best case, C) code (more concretely, they’re using a library called Paramesh, written in Fortran 90). Now, if you ask anyone with a solid background in computer science, she will probably tell you that nobody (except physicists, that is) programs these days in Fortran. We know better languages, and have developed far better ways of writing computer programs in the 50 years since Fortran was invented. That is, we physicists are using obsolete technologies. Those newer languages (Scheme, Haskell, OCaml, and so on and so forth) are better in many ways, but specially in one that i am sure is close to any physicist’s heart: they provide far, far better means of abstraction. That is, you can write much shorter programs in a language that is conceptually closer to the problem at hand. And shorter may well mean something like a ten fold reduction in the number of lines of code; not to mention the benefits on clarity, maintanability and extensibility that greater abstraction entails. To use a metaphor, it’s like we were using Levi-Civita’s books and notation as our standard way of calculating in General Relativity, instead of modern differential geometry.
Of course, there’re perfectly understandable reasons for our using antics like Fortran, legacy code being probably the most prominent one; and physicists not having the needed expertise might well an important one too (but let me rush to say that efficiency of the code is not a good excuse these days). But i’m convinced that numerical physics would be vastly improved if we imported some expertise from the professional computing world. I’m told by friends in the field that some of the most ‘advanced’ guys are trying things like C++ and Java (instead of Fortran) these days: i’m sorry, but these languages were current some 20 years ago, and we’ve learnt since then how to avoid many of the pitfalls and unnecessary complexities they carry on. Much more interesting is to use interactive languages like Python (to be on the conservative side) or, if you ask me, functional languages like Scheme or Haskell. To give you a glimpse of what i’m talking about, here is how you’d write quicksort in Fortran 95; in Haskell, it’s a two-liner:
qsort  =  qsort (x:xs) = qsort (filter (< x) xs) ++ [x] ++ qsort (filter (>= x) xs)
Fortunately, not every one sticks to Fortran these days: Michele Vallisneri’s Synthetic LISA is a beautiful example of a step in the right direction, and i’m glad to see that numerical libraries like PETSC do in fact provide Python bindings. But, as i said, i think (after nine years or so of earning a living writing computer programs) that there are even better ways. As a matter of fact, i’m seriously considering the possibility of writing some LISA simulation code using Scheme. What deters me, besides lack of time, is the enormous weight of tradition: almost everybody out there in the physics community is using those C and Fortran libraries, and that means millions of lines of well-tested code and wondrous results like those black hole mergers. The easiest thing to do is to go with the flow, but still…
(By the way, these comments are by no means intended as a critique of Baker et al. work, which is impressive as it stands. Besides, for what i know, they may well be using far more sophisticated techniques than plain Fortran or C. My rants are more geared towards many other cases i’ve seen of physics programmers which were anything but sophisticated.)