Getting dusty…

I have been getting dusty for some time now.

… No, by “dusty” I don’t mean that dust of the “Heat and Dust” kind, even though it’s been quite the regular kind of an “unusually hot” summer this year, too.

[In case you don’t know, “Heat and Dust” was a neat movie that I vaguely recall I had liked when it had come on the scene some 2–3 decades ago. Guess I was an undergrad student at COEP or so, back then or so. (Google-devataa now informs me that the movie was released in 1983, the same year that I graduated from COEP.)]

Anyway, about the title of this post: By getting dusty, I mean that I have been trying to get some definite and very concrete beginning, on the software development side, on modelling things using “dust.” That is, particles. I mean to say, methods like molecular dynamics (MD), smoothed particle hydrodynamics (SPH), lattice Boltzmann method (LBM), etc. … I kept on postpoing writing a blog post here with the anticipation that I would succeed in tossing a neat toy code for a post.

… However, I soon came face-to-face with the  sobering truth that since becoming a professor, my programming skills have taken a (real) sharp dip.

Can you believe that I had trouble simply getting wxWidgets to work on Ubuntu/Win64? Or to get OpenGL to work on Ubuntu? It took almost two weeks for me to do that! (And, I haven’t yet got OpenGL to work with wxWidgets on Ubuntu!) … So, finally, I (once again) gave up the idea of doing some neat platform-neutral C++ work, and, instead (once again) went back to VC++. Then there was a problem awaiting me regarding VC++ too.

Actually, I had written a blog post against the modern VC++ at iMechanica and all (link to be inserted) but that was quite some time back (may be a couple of years or so). In the meanwhile, I had forgotten how bad VC++ has really become over the years, and I had to rediscover that discouraging fact once again!

So, I then tried installing VC++ 6 on Win7 (64-bit)—and, yet another ugly surprise. VC++ 6 won’t run on Win7. It’s possible to do that using some round-about ways, but it all is now a deprecated technology.

Finally, I resigned myself to using VC++ 10 on Win7. Three weeks of precious vacation time already down!

That‘s what I meant when I said how bad a programmer I have turned out to be, these days.

Anyway, that’s when I finally could begin writing some real though decidedly toy code for some simple MD, just so that I could play around with it a bit.

Though writing MD code seems such a simple, straight-forward game (what’s the Verlet algorithm if not plain and simple FDM?), I soon realized that there are some surprises in it, too.

All of the toy MD code to be found on the Internet (and some of the serious code, too) assumes only a neat rectangular or a brick-like domain. If you try to accommodate an arbitrary shaped domain (even if only with straight-lines for boundaries), you immediately run into the run-time efficiency issues. The matter gets worse if you try to accommodate holes in the domain—e.g., a square hole in a square domain like what my Gravatar icon shows. (It was my initial intention to quickly do an MD simulation for this flow through square cavity having a square hole.)

Next, once you are able to handle arbitrarily shaped domains with arbitrarily shaped holes, small bugs begin to show up. Sometimes, despite no normal flux condition, my particles were able to slip past the domain boundaries, esp. near the points where two bounding edges meet. However, since the occurrence was rare, and hard to debug for (what do you do if it happens only in the 67,238th iteration? hard-code a break after 67,237, recompile, run, go off for a coffee?) I decided to downscale the complexity of the basic algorithm.

So, instead of using the Lennard-Jones potential (or some other potential), I decided to switch off the potential field completely, and have some simple perfectly hard and elastically colliding disks. (I did implement separate shear and normal frictions at the time of collisions, though. (Of course, the fact that frictions don’t let the particles be attracted, ever, is a different story, but, remember, we are trying to be simple here.)) The particles were still managing to escape the domain, at rare times. But at least, modelling with the hard elastic disks allowed me to locate the bugs better. Turned out to be a very, very stupid-some matter: my algorithm had to take care for a finite-sized particle interacting with the boundary edges.

But, yes, I could then get something like more than 1000 particles happily go colliding with each other for more than 10^8 collisions. (No, I haven’t parallelized the code. I merely let it run on my laptop while I was attending in a longish academics-related meeting.)

Another thing: Some of the code on the ‘net, I found, simply won’t work for even a modestly longer simulation run.  For instance, Amit Kumar’s code here [^]. (Sorry Amit, I should have written you first. Will drop you a line as soon as this over-over-overdue post is out the door.) The trouble with such codes is with the time-step, I guess. … I don’t know for sure, yet; I am just loud-thinking about the reason. And, I know for a fact that if you use Amit’s parameter values, a gas explosion is going to occur rather soon, may be right after something like 10^5 collisions or so. Java runs too slowly and so Amit couldn’t have noticed it, but that’s what happens with those parameter values in my C++ code.

I haven’t yet fixed all my bugs, and in fact, haven’t yet implemented the Lennard-Jones model for the arbitrarily shaped domains (with (multiple) holes). I thought I should first factor out the common code well, and then proceed. … And that’s when other matters of research took over.

Well, why did I get into MD, in the first place? Why didn’t I do something useful starting off with the LBM?

Well, the thing is like this. I know from my own experience that this idea of a stationary control volume and a moving control volume is difficult for students to get. I thought it would be easy to implement an MD fluid, and then, once I build in the feature of “selecting” (or highlighting) a small group of molecules close to each other, with these molecules together forming a “continuous” fluid parcel, I could then directly show my students how this parcel evolves with the fluid motion—how the location, shape, momentum and density of the parcel undergoes changes. They could visually see the difference between the Eulerian and Lagrangian descriptions. That, really speaking, was my motivation.

But then, as I told you, I discovered that I have regressed low enough to have become a very bad programmer by now.

Anyway, playing around this way also gave me some new ideas. If you have been following this blog for some time, you would know that I have been writing in favor of homeopathy. While implementing my code, I thought that it might be a good idea to implement not just LJ, but also the dipole nature of water molecules, and see how the virtual water behaves: does it show the hypothesized character of persistence of structure or not. (Yes, you read it here first.) But, what the hell, I have far too many other things lined up for me to pursue this thread right away. But, sure, it’s very interesting to me, and so, I will do something in that direction some time in future.

Once I get my toy MD code together (for both hard-disk and LJ/other potentials models, via some refactorability/extensibility thrown in), then I guess I would move towards a toy SPH code. … Or at least that’s what I guess would be the natural progression in all this toys-building activity. This way, I could reuse many of my existing MD classes. And so, LBM should logically follow after SPH—what do you think? (And, yes, I will have to take the whole thing from the current 2D to 3D, sometime in future, too.)

And, all of that is, of course, assuming that I manage to become a better programmer.

Before closing, just one more note: This vacation, I tried Python too, but found that (i) to accomplish the same given functional specs, my speed of development using C++ is almost the same (if not better) than my speed with Python, and (ii) I keep missing the nice old C/C++ braces for the local scope. Call it the quirky way in which an old programmer’s dusty mind works, but at least for the time being, I am back to C++ from that brief detour to Python.

Ok, more, later.

* * * * *   * * * * *   * * * * *

A Song I Like:
(Marathi) “waaT sampataa sampenaa, kuNi waaTet bhetenaa…”
Music: Datta Dawajekar
Lyrics: Devakinandan Saaraswat
Singer: Jayawant Kulkarni

[PS: A revision is perhaps due for streamlining the writing. May be I will come back and do that.]