GSOC has finally come to a close. It seems like yesterday when I was blogging on here, young and foolish about the complexity of my project. What you see above are my modest results. Pretty as it might look, this fluid is most definitely code-name O(n^2), and that's not really a good thing! Let me explain the name in the three major aspects of my project.

Realism

It's one thing to implement a smoothed particle hydrodynamics paper in STEP, that was basically done by the midterm, but how can you make it... well, physical? In an SPH system, there are several parameters that must be adjusted to achieve that. Namely, the temperature dependent gas constant, the rest density, the smoothing kernel length, the surface tension coefficient, the viscosity constant, and the universal bounciness of object-object interactions. With the exception of the surface tension and bounciness, I was able to determine a number of these constants for a working fluid using trial and error.

Most significant for computation was the adjustment of the smoothing kernel length. This is the radius around each fluid particle at which the effects of the particle are spread out. According to the SPH algorithm I implemented, at each time step we compute the pressures at each fluid particle based on the density of all fluid particles within the smoothing kernel radius. That being said, at each time-step I made O(n^2) distance to fluid particle comparisons, and in the limit of an infinitely large smoothing kernel there would be O(n^2) arithmetic operations. Among some other issues I will discuss, I found this to be a limiting factor for the number of fluid particles I could use. 

Unfortunately, this prevented me from testing an important aspect of fluids, surface tension. By adding some additional fake surface tension forces outlined in Muller 2003, my fluid would be calmer and perhaps even support the weight of a boat at it's surface. Additionally, the addition of bounciness would allow gravity to act realistically on the fluid and have density collect at the bottom of a jar.

So, how can one improve on this brute-force enumeration of all pairs of particles? My expert mentor Vladimir suggested a new data structure might be the best approach. Since both fluids and gases benefit from a large number of particles, it would be advantageous to implement a Axis-Aligned Bounding Box data structure. Not only would this help for collision detection between any Step objects but this would provide a way to only calculate pressure forces for neighboring fluid particles. This works because we store a sorted list of boxes about both our X an Y axes and keep track of box overlaps. Unfortunately, my mentor and I did not have enough time to implement this complex rehaul to StepCore.

Measurement

Step is about education. For educational purposes, at the very least, it would be great to be able to know the average density and pressure within a region of the fluid. Despite my fluid being somewhat artificial and hacked together with strange internal values, I managed to achieve this. However, the computational complexity is, once again, not pretty. 

The calculation of pressure at a particle point in the World was already used for the dynamics of the fluid. I merely had to generalize it for discrete area elements of Step's "measurement rectangle". By moving the rectangle and resizing it within the Fluid, each time-step I am sliced the area into a grid, calculated the pressure and density at the center of each area, and then added it all together.

The trouble with that is probably quite clear. If you have a very large measurement rectangle you may want a widely spaced grid and for a small measurement rectangle you may want a closely spaced grid. A grid dependent on the size of the rectangle and the smoothing kernel radius is also a difficult value to fine tune.

In the future,  It would be ideal for the user to be able to select the precision of the pressure and density calculation at their own discretion, and have the variance values adjust accordingly. I was unable to implement variance calculations for these measurements, but that may also be an expensive calculation!

Visualization

Finally, the Pièce de résistance. The jaw-dropping visuals that all fluids deserve. As you can see above, there's lots of room for improvement. Visualization can be done in two ways. Muller 2003, the paper I was following for this SPH implementation, suggests a fancy isosurface calculation. That is, determine the normal to the fluid particle density field at every point. If this value is greater than some threshold, it's safe to say it's a surface particle and it could be used to draw a sexy looking surface. 

However, with a deadline looming, my mentor suggested a simpler approach. I would simply calculate density at each point in a grid and then paint a particle there with opacity based on the density. As seen above, these opacity calculations weren't quite calibrated, since we don't see much variation in blue except at the edges of the fluid. 

I had a problem though. Where should I "draw the line" for drawing the dots? I can't very well calculate the density at each point in the observable screen area. It would take forever to get a smooth fluid. The solution is to calculate the minimum bounding box for the fluid and use a cut-off value in-case one fluid particle flies away off the screen. This ended up being quite challenging for some reason, so I decided to only render the fluid within the "measurement rectangle".

I had an issue with converting between coordinate systems that still remains unsolved. I needed to draw my fluid about the origin of the WorldScene or else my fluid would render in the incorrect location. This is likely due to a simple mapping problem that I intend to figure out soon.

End

So, what's next for Fluids, Step, and Me? The realism, measurement, and visualization issues outlined above will continue to be worked on by Vladimir and I. After some of the fundamental problems are addressed, I forsee a very liquidy Step in the near future. Me? I'm starting a MSc. in Physics at University of Waterloo this Fall, but I had a lot of fun on this project and hope to continue making Step a richly featured open-source physics simulator!

It's about time I popped this bubble of silence.

I've been working hard on two programming projects this summer, namely my Google Summer of Code project and a fancy upgrade to some molecular dynamics software. In my previous posts, I looked at the back-end of Step and some of the mathematics of smoothed particle hydrodynamics. There are still plenty of outstanding problems in those areas but I'll address those gradually in the coming weeks. For this post, I'll give an overview of all the GUI/user interactivity stuff that i'll have to tackle on my quest to implement fast fluid simulation.

I'll just enumerate the ways the user can interact with a fluid and trace through exactly what is going on behind the scenes. I know, KDevelop/IDEs are pretty fancy, but I'm oldschool so I tend to just follow the code execution manually. It's kind of like a grep-based treasure hunt! So, I'll just give some running commentary with bonus screenshots from the current Gas classes. I doubt many readers of this blog will find this interesting, but my dream is that this post series may someday be useful for a new Step developer =P 

Oh, if you're a Qt newbie just remember that any class with a Q is a Qt class, and don't forget your slots and signals.

Read the rest of this entry »

The trials and tribulations of uncertainty run deep in the hearts of experimentalists. From our biased human minds, to the imprecision of our measurement devices, to the quantum uncertainty of our universe, errors have confounded scientists for eons. Even within the safe and deterministic world of computation, inexactness is rampant.

There are two main types of numerical errors. The first of which is result of computers having a limited number of bits (32 or 64) to represent numbers. All the coolest real numbers have an infinite number of digits so we're restricted to representing numbers that differ by machine epsilon.

The second source of numerical error pops up a lot in numerical algorithms, including those implemented in Step for solving differential equations. It would be nice if we could perform exact symbolic computation all the time, but physics can get messy, especially with errors! Almost any time we take a derivative on computer, we approximate it using a taylor series. This approximation works great if we keep an infinite number of terms but who has the time to calculate all that? Thus, our results differ from the exact answer based on where we truncate our taylor series. This truncation error combined with the floating point error mentioned above, has the potential to cause serious numerical instabilities and pain for computational scientists.

Numerical issues aside*, in true experimentalist fashion, Step allows users keep track of the propagation of user-inputted uncertainties over the course of a simulation. The mathematics behind this process is reviewed briefly here. As an example let's look at the calculation for the uncertainty in kinetic energy in the Step's ParticleErrors class,

 

double ParticleErrors::kineticEnergyVariance() const {
  return _velocityVariance.dot(particle()->velocity().cwise().square()) * square(particle()->mass())        + square(particle()->velocity().squaredNorm()/2) * _massVariance;
}

 

The kinetic energy (KE) of a particle is a function of both mass and velocity, where each of these variables have their respective uncertainties. In order to calculate the variance of kinetic energy we have to take the sum of d(KE)/dv multiplied by our velocityVariance and d(KE)/dm multiplied by our massVariance. The function only looks a bit confusing because we want a scalar value and we have to use a few Eigen functions to deal with 2D vectors. 

The trouble with smoothed particle hydrodynamics is that each "fluid particle" is inherently an approximation of many particles or a bulk region of fluid. Nonetheless, the calculated densities and pressures for each chunk of fluid are subject to uncertainty. As outlined by Muller et al., any scalar quantity A can be calculated by summing over all particles j with some smoothing kernel defined by W:

Given the user-inputted uncertainty in particle mass m and whatever calculated uncertainty we have for our quantites A and p, this formula can be used with the method above to calculate uncertainty in the newly calculated quantity A. These calculations will be coded in the FluidParticleError class outlined below. Note the addition of the member variables for densityVariance and pressureVariance. These values could be calculated on the fly, but in doing there would be a large amount of calculations due to the summation of all nearby particle errors. Later in this project this will also include Viscosity error calculations that depend on the velocity of nearby particles using the same equation above. 

It's a week into my GSOC project and I still have lots of work to do. Next post will be a bit more software development oriented, as I'll look more at the Qt GUI and how it connects with the numerical back-end of Step. Expect a post shortly!

*You can adjust the precision of your Solver in the properties dialog box of Step!

I went to the Sharcnet GPU and Cell Computing Symposium at my university last week. One of the keynote speakers, Dr. Paul Woodward, pretty much blew my weak and feeble mind. He used Roadrunner, a.k.a. the faster supercomputer in the world, to compute some impressive fluid dynamics. You know your research group is on the right track when you have a Wallpapers section on your web site.

Back on planet earth, us computational plebs don't have such luxurious resources at our finger tips. To start with, you don't have a lot of options when you feel like simulating fluids on a computer. Things get even more depressing when you want fluid flow computed in real time. In video games, or maybe just in old video games, accurate water simulation ends up being avoided all together with some elaborate ruse to keep computation costs down. Why do fluids have such a bad rep? To paraphrase Wikipedia: "We suck at solving systems of non-linear partial differential equations, here's a million dollars if you figure it out, good luck lol".

I find this whole water simulation extra-depressing because the universe computes the solution to these problems instantly and literally rubs it in our face whenever the wind blows (via turbulence). Just imagine with me what life would be like if we had an exact solution for fluid flow. Picture weather forecasts that actually forecast, reliable message-in-a-bottle delivery using ocean currents, and aerodynamic everything.

Wow, nice introduction. This post is actually about my fluid simulation project for Google Summer of Code. To achieve real-time fluid simulation in Step, I intend to follow the Smoothed Particle Hydrodynamics (SPH) algorithm outlined by Muller et al. in "Particle-Based Fluid Simulation for Interactive Applications". This method works by breaking up a large fluid body into discrete chunks and using these chunks to create smoothed out fields. In fluids, the stuff we care about are pressure, density, and velocity fields. SPH is just a fancy way of figuring out the forces that act on your fluid by taking these fields into account.

Without going into all the mathematical detail presented by Muller et al., we can begin to define our fluid particles class in StepCore, the mathematical back-end of Step. Note that all of this is pretty fuzzy right now since I haven't started coding anything. Many important properties of FluidParticle already exist the base class Particle in Step as depicted below. Nothing too special here, just get and set methods for the particle member variables. One thing to note is that FluidParticle has a specific radius that defines the core radius of each particle. However, at this point it is unlikely that the radius of each particle will vary. I'll talk about Errors in a separate post also.

Now that the FluidParticle is defined we can look at the collective group, the Fluid. All "bodies and forces" in Step are derived from the root class Item. A Fluid is simply a group of FluidParticle items in a handy "Item Group". The main purpose of this group is to extract macroscopic features of our fluid particles that exists within a rectangle drawn in the GUI.

This rectangular data inquiry box was designed for use with the Gas class where determining the number of particles, the kinetic energy, and the temperature are meaningful quantities for statistical mechanics. However, I'm not really sure how well the physics of statistical mechanics translates to the mathematics of a fluid. For instance, I may know the velocities of all my fluid chunks but is it correct to use these velocities to determine the so-called "temperature" of my fluid? I suppose so... but I have to do a bit more research.

Nonetheless, for now I'll just use the rectPressure(), rectVelocity() and rectDensity(). I'd like to be able to use these methods to not just calculate the average pressure, velocity, and density by looping over the values stored in each fluid chunk but averaging over the entire field created by the fluid chunks. I'm not sure how computationally demanding this calculation will be, but if I discretized space within the fluid into some sort of grid this may be possible. 

So now we have a Fluid that is composed of individual FluidParticles. Where does the physics come in? We need to look at the FluidForce class for that. The purpose of the FluidForce class is to do the work of applying the inter-fluid forces to each particle in our item group by means of the calcForce method that allows you to calculate the variances based on a boolean flag. In the process of calculating the forces on each fluid particle, I'll have to use the helper method calcPressureDensity() as well.

You'll notice that the only member variable for the fluidForce is a cutoff value and FluidErrors (which I will elaborate on later). This value represents a limiting distance for which the force contribution from a distant particle would be so small that it's not worth calculating. In the current Gas implementation in Step, calcForce requires O(n^2) comparisons, as we have to loop through all pairs of particles to determine if they have a distance less than the cutoff value. As mentioned briefly in Muller et al., a potential speed-up over this approach might be to place our fluid particles on a grid and only calculate the force contribution from particles attached to neighboring grid cells to our particle of interest. I'll address performance issues later in my project.

I've already overlooked a few critical issues that may require the modification of these classes. A fundamental characteristic of a fluid is viscosity, which is defined as resistance of a fluid to shear and stress forces. Obviously, my GSOC project will not be complete unless I can simulate some delicious sticky honey. At this stage, for simplicity, I'll probably just hard code some viscosity values. The second critical issue will be collisions. My mentor and I both agree that collision handling may get tricky. If you've ever been stuck in a wall or a tree when playing a video game you'll know what I mean.

I'll address that problem eventually but next I'll briefly discuss the Error classes associated with the fluid objects. The calculation of errors/variances is one of the cool features that makes Step unique!

After that, I'll explore the Qt and GUI side of Step. I'll connect the classes I've defined above to user interactions. Once I do that, things should get a lot clearer in terms of how the FluidForce, Fluid and FluidParticles work in conjunction with the World.

Step is a 2D open-source physics simulator that was recently added to the KDE Education project.

Vladimir, the original developer and former Google Summer of Code student, posted up a few ideas to improve Step for this years Google Summer of Code. He included one about improving fluid simulation and, being a computational physicist (TM), I decided to go for it. I wrote up a fancy proposal a few weeks ago and I'm happy to say that it was accepted.

A lot of GSOC projects are strictly geeking out behind-the-scenes, but this one is fun, visual, and physics-based. It's perfect for me.

Compared to Vladimir's original summer of code, spent designing Step in it's entirety, my proposal seems like a pretty modest addition. However, fluids are traditionally a very messy subject in the world of physics simulators. If you aren't convinced, read this thread on the Box2D forums from only a month ago. To quote ElectroDruid on page 9:

"I gave up with having the particles be handled by Box2D, and one of the main reasons was because of the amount of pairs generated. You can up the maximum limits, but eventually you just run out of storage space for the numbers of pairs that can be needed in certain cases. Box2D is great for rigid bodies, but for fluids where you need a lot of particles which interact with each other in close proximity all the time, it's easy to push the limits to breaking point."

I'm sure I will have to deal with similar problems encountered in that thread, but I'm optimistic that such issues can be resolved. Throughout the summer I'll be blogging my progress on the project here, so stay tuned.