Over the past couple of weeks or so, I’ve been going over SPH (smoothed particle hydrodynamics).
I once again went through the beginning references noted in my earlier post, here [^]. However, instead of rewriting my notes (which I lost in the last HDD crash), this time round, I went straight to programming. … In this post, let me
recap recall what all I did.
First, I went through the great post “Why my fluids don’t flow” [^] by Tom Madams. … His blog has the title: “I am doing it wrong,” with the sub-text: “most of the time.” [Me too, Tom, me too!] This post gave a listing of what looked like a fully working C++ code. Seeing this code listing (even if the videos are no longer accessible), I had to try it.
So, I installed the Eclipse CDT. [Remember my HDD crash?—the OS on the new HDD had no IDEs/C++ compilers installed till now; I had thus far installed only Python and PyCharm]. I also installed MinGW, freeglut, Java JRE, but not all of them in the correct order. [Remember, I too do it wrong, most of the time.] I then created a “Hello World” project, and copy-pasted Tom’s code.
The program compiled well. [If you are going to try Tom’s code in Eclipse CDT + MinGW on Windows, the only issue you would now (in 2016) run into would be in the Project Settings, both in the compiler and linker settings parts, and both for OpenGL and freeglut. The latest versions of Eclipse & MinGW have undergone changes and don’t work precisely as explained even in the most helpful Web resources about this combination. … It’s not a big deal, but it’s not exactly what the slightly out-of-date online resources on this topic continue telling you either. … The options for the linker are a bit trickier to get than those for the compiler; the options for freeglut certainly are trickier to get than those for OpenGL. … If you have problems with this combination (Eclipse + MinGW on Windows 7 64-bit, with OpenGL and freeglut), then drop me a line and I will help you out.]
Tom’s program not only compiled well, but it also worked beautifully. Quite naturally, I had to change something about it.
So I removed his call to glDrawArrays(), and replaced the related code with the even older glBegin(GL_POINTS), glVertex2d(), glEnd() sort of a code. As I had anticipated, there indeed was no noticeable performance difference. If the fluid in the original code required something like a minute (of computer’s physical time) to settle down to a certain quiescent state, then so did the one with the oldest-style usage of OpenGL. The FPS in the two cases were identical in almost all of the release runs, and they differed by less than 5–7% for the debug runs as well, when the trials were conducted on absolutely fresh cold-starts (i.e. with no ready-to-access memory pages in either physical or virtual memory).
Happy with my change, I then came to study Tom’s SPH code proper. I didn’t like the emitters. True to my core engineering background, what I wanted to simulate was the dam break. That means, all the 3000 particles would be present in the system right from the word go, thereby also having a slower performance throughout, including in the beginning. But Tom’s code was too tied up with the emitters. True to my software engineering background, rather than search and remove the emitters-related portion and thus waste my time fixing the resulting compiler errors, I felt like writing my own code. [Which true programmer doesn’t?]
So I did that, writing only stubs for the functions involving the calculations of the kernels and the accelerations. … I, however, did implement the grid-based nearest-neighbor search. Due to laziness, I simply reused the STL lists, rather than implementing the more basic (and perhaps slightly more efficient) “p->next” idiom.
Then I once again came back to Tom’s code, and began looking more carefully at his SPH-specific computations.
What I now didn’t like was the variables defined for the near-density and the near-pressure. These quantities didn’t fit very well into my preconceived notions of how a decent SPH code ought to look like.
So, I decided to deprove [which word is defined as an antonym of “improve”] this part, by taking this 2010 code from its 2007 (Becker et al.) theoretical basis, to a 2003 basis (Muller et al., Eurographics).
Further following my preconceived notions, I also decided to keep the values of the physical constants (density, gas stiffness, viscosity, surface tension) the same as those for the actual water.
The code, of course, wouldn’t work. The fluid would explode as if it were a gas, not water.
I then turned my learner’s attention to David Bindel’s code (see the “Resources” section at the bottom of his page here [^]).
Visiting Bindel’s pages once again, this time round, I noticed that he had apparently written this code only as a background material for a (mere) course-assignment! It was not even an MS thesis! And here I was, still struggling with SPH, even after having spent something like two weeks of full-time effort on it! [The difference was caused by the use of the realistic physical constants, of course. But I didn’t want to simply copy-paste Tom’s or Bindel’s parameter values; I wanted to understand where they came from—what kind of physical and computational contexts made those specific values give reasonable results.]
I of course liked some of the aspects of Bindel’s code better—e.g. kernels—and so, I happily changed my code here and there to incorporate them.
But I didn’t want to follow Bindel’s normalize_mass routine. Two reasons: (i) Once again according to my preconceived notions, I wanted to first set aside a sub-region of the overall domain for the fluid; then decide with how many particles to populate it, and what lattice arrangement to follow (square? body centered-cubic? hexagonal close-packed?); based on that, calculate each particle’s radius; then compute the volume of each particle; and only then set its mass using the gross physical density of the material from which it is composed (using the volume the particle would occupy if it were to be isolated from all others, as an intermediate step). The mass of a particle, thus computed (and not assumed) would remain fixed for all the time-steps in the program. (ii) I eventually wanted a multi-phase dam-break, and so wasn’t going to assume a global constant for the mass. Naturally, my code wouldn’t be able to blindly follow Bindel on all counts.
I also didn’t like the version of the leapfrog he has implemented. His version requires you to maintain additional quantities of the velocities at the half time-steps (I didn’t mind that), and also to implement a separate leapfrog_start() function (which I did mind—an additional sequence of very similar-looking function calls becomes tricky to modify and maintain). So, I implemented the other version of the leapfrog, viz., the “velocity Verlet.” It has exactly the same computational properties (of being symplectic and time-reversible), the same error/convergence properties (it too is second-order accurate), but it comes with the advantage that the quantities are defined only at the integer time-steps—no half-time business, and no tricky initialization sequence to be maintained.
My code, of course, still didn’t work. The fluid would still explode. The reason, still, was: the parameter values. But the rest of the code now was satisfactory. How do I know this last part? Simple. Because, I commented out the calls to all the functions involving all other accelerations, and retained only the acceleration due to gravity. I could then see the balls experiencing the correct free-fall under gravity, with the correct bouncing-back from the floor of the domain. Both the time for the ball to hit the floor as well as the height reached after bouncing were in line with what physics predicts. Thus I knew that my time integration routines would be bug-free. Using some debug tracings, I also checked that the nearest-neighbour routines were working correctly.
I then wrote a couple of Python scripts to understand the different kernels better; I even plotted them using MatPlotLib. I felt better. A program I wrote was finally producing some output that I could in principle show someone else (rather than having just randomly exploding fluid particles). Even if it was doing only kernel calculations and not the actual SPH simulation. I had to feel [slightly] better, and I did.
At this stage, I stopped writing programs. I began thinking. [Yes, I do that, too.]
To cut a long story short, I ended up formulating two main research ideas concerning SPH. Both these ideas are unlike my usual ones.
Usually, when I formulate some new research idea, it is way too conceptual—at least as compared to the typical research reported in the engineering journals. Typically, at that stage (of my formulation of a new research idea), I am totally unable to see even an outline of what kind of a sequence of journal papers could possibly follow from it.
For instance, in the case of my diffusion equation-related result, it took me years before an outline for a good conference paper—
almost like really speaking, at par with a journal paper—could at all evolve. I did have the essential argument ready. But I didn’t know what all context—the specifically mathematical context—would be expected in a paper based on that idea. I (and all the mathematicians I contacted) also had no idea as to how (or where) to go hunting for that context. And I certainly didn’t have any concrete idea as to how I would pull it all together to build a concrete and sufficiently rigorous argument. I knew nothing of that; I only knew that the instantaneous action-at-a-distance (IAD) was now dead; summarily dead. Similarly, in case of QM, I do have some new ideas, but I am still light-years away from deciding on a specific sequence of what kind of papers could be written about it, let alone have a good, detailed idea for the outline of the next journal paper to write on the topic.
However, in this case—this research on SPH—my ideas happen to be more like
what [other] people typically use when they write papers for [even very high impact] journals those which lie behind the routine journal papers. So, papers should follow easily, once I work on these ideas.
Indeed, papers must follow those ideas. …There is another reason to it, too.
… Recently, I’ve come to develop an appreciation, a very deep kind of an appreciation, of the idea of having one’s own Google Scholar page, complete with a [fairly] recent photo, a verified email account at an educational institution (preferably with a .edu, or a .ac.in (.in for India) domain, rather than a .org or a .com domain), and a listing of one’s own h-index. [Yes, my own Google Scholar page, even if the h-Index be zero, initially. [Time heals all wounds.] I have come to develop that—an appreciation of this idea of having a Google Scholar page. … One could provide a link to it from one’s personal Web site, one could even cite the page in one’s CV, it could impress UGC/NBA/funding folks…. There are many uses to having a Google Scholar page.
…That is another reason why [journal] papers must come out, at least now.
And I expect that the couple of ideas regarding SPH should lead to at least a couple of journal papers.
Since these ideas are more like the usual/routine research, it would be possible to even plan for their development execution. Accordingly, let me say (as of today) that I should be able to finish both these papers within the next 4–5 months. [That would be the time-frame even if I have no student assistant. [Having a student assistant—even a very brilliant student studying at an IIT, say at IIT Bombay—would still not shorten the time to submission, neither would it reduce my own work-load any more than by about 10–20% or so. That’s the reason I am not planning on a student assistant on these ideas.]
But, yes, with all this activity in the recent past, and with all the planned activity, it is inevitable that papers would fall out. Papers must, in fact, fall out. …. Journal papers! [Remember Google Scholar?]
Of course, when it comes to execution, it’s a different story that even before I begin any serious work on them, I still have to first complete writing my CFD notes, and also have to write a few FDM, FVM and VoF/LevelSet programs scripts or OpenFOAM cases. Whatever I had written in the past, most of it was lost in my last HDD crash. I thus have a lot of territory to recover first.
Of course, rewriting notes/codes is fast. I could so rapidly progress on SPH this year—a full working C++ code in barely 2–3 weeks flat—only because I had implemented some MD (molecular dynamics) code in 2014, no matter how simple MD it was. The algorithms for collision detection and reflections at boundaries remain the same for all particles approaches: MD with hard disks, MD with LJ potential, and SPH. Even if I don’t have the previously written code, the algorithms are no longer completely new to me. As I begin to write code, the detailed considerations and all come back very easily, making the progress very swift, as far as programming is concerned.
When it comes to notes, I somehow find that writing them down once again takes almost the same length of time—just because you had taken out your notes earlier, it doesn’t make writing them down afresh takes any less time.
Thus, overall, recovering the lost territory would still take quite some effort and time.
My blogging would therefore continue to remain sparse even in the near future; expect at the most one more post this month (May 2016).
The work on the journal papers itself should begin in the late-June/early-July, and it should end by mid-December. It could. Nay, it must. … Yes, it must!
Papers must come out of all these activities, else it’s no research at all—it’s nothing. It’s a zero, a naught, a nothing, if there are no papers to show that you did research.
… Papers must fall out! … Journal papers!!
A Song I Like:
(Western, Instrumental) “The rain must fall”
May be one quick editing pass later today, and I will be done with this post. Done on 12th May 2016.]