# Stable Fluids - Jos Stam - GDC03 code

This topic is 2768 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi,

I try to fully understand the Stable Fluids method of Jos Stam. I am experimenting with the code of Stam's GDC03 implementation - paper .

In the process of understanding the solver it's vital to learn the behaviour of one step only - density diffusion, density advection and so on.
So, I would like to just play with the density solver only (no velocity update and no velocities whatsoever).
static void idle_func ( void ){	get_from_UI ( dens_prev, u_prev, v_prev );	//vel_step ( N, u, v, u_prev, v_prev, visc, dt );	dens_step ( N, dens, dens_prev, u, v, diff, dt );	glutSetWindow ( win_id );	glutPostRedisplay ();}

With no active velocities I guess this is now just a diffusion solver. And its working as it should in my mind.
Note: default value for density diffusion is zero so one have to modify it to make it work - my diffusion coefficient is 0.00075f.

But with no velocities I think it's ok to neglect the advection term in the density step.
void dens_step ( int N, float * x, float * x0, float * u, float * v, float diff, float dt ){	add_source ( N, x, x0, dt );	SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt );	//SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, dt );}

But this seems to be incorrect. The density values are frozen as they are added from the user input, no diffusion at all - and I don't know why.

I would really appreciate if someone could point out what I miss.

##### Share on other sites

It might be that it's because diffusion causes self-advection (search for the term in the paper). Diffusion automatically implies that different parts of the cloud do/will have non-zero relative speed with respect to each other.

As you've already perceived, this self-generated velocity due to cloud expansion is different than velocity caused by mouse/external-driven forces.

In the end, it seems to me that it's all just velocity no matter what the source is, and so you've always got to use the advect() function.

You might want to email Stam for further information. He might reply, you never know.

If you're not familiar with entropy, it might be an interesting subject for you. Diffusion causes randomly directed (e.g., disordered) velocity vectors, where an external force causes coherently directed (e.g., ordered) velocity vectors.

Randomness seems natural: The molecules of a book falling in a gravitational field at a very high speed have very coherent velocity vectors. When the book hits the ground all of that kinetic energy is turned to heat and sound waves, which disperse in random directions. As well, the velocity vectors of the books molecules also randomize at this time. You will (practically) never find the reverse of this -- heat and sound waves from the environment don't magically come together to force the book's molecules to take up a coherent velocity and send the book shooting straight up into the sky. Why? Because it seems that there's no such thing as anti-diffusion.

This is what Eddington called the "arrow of time"... Disorder naturally tends to increase over time in an isolated system such as the Universe, not decrease.

[Edited by - taby on July 20, 2010 2:36:06 PM]

##### Share on other sites
Thanks for your thoughts.

Self-advection can be interpreted as “velocity should move along itself” - as Stam mentioned in his papers. This self-advection term can be found in the velocity update.

In the density update advection moves the density (dye) through the known velocity field. I can't identify any term that could generate velocity in the density update (dye diffusion + dye advection) - from theoretical and from coding point of view as well.

But I won't give in so easily, hopefully I can post the answer in the near future.

##### Share on other sites
Quote:
 Original post by TooropThanks for your thoughts.Self-advection can be interpreted as “velocity should move along itself” - as Stam mentioned in his papers. This self-advection term can be found in the velocity update.In the density update advection moves the density (dye) through the known velocity field. I can't identify any term that could generate velocity in the density update (dye diffusion + dye advection) - from theoretical and from coding point of view as well.But I won't give in so easily, hopefully I can post the answer in the near future.

Yes, sorry, I think I meant to say self-generated self-advection. If you keep the advect call in, do the vectors of the velocity field gradually become non-zero in length in the region surrounding the place where you put the initial bit of gas, even when there are zero external forces? If so, then it must be diffusion that's responsible for this apparently spontaneous creation of velocity on the macroscopic scale -- this is how it works in the real world anyway. Commenting out the diffuse call should confirm or deny that notion.

##### Share on other sites
Maybe I can't describe the situation in an adequate way or I miss your point.

There are no forces at initialization and I won't give any external forces to the simulation while on the run - I can check this because no velocity vectors can be exhibited in the grid.

I add simply dye density sources and with a positive dye diffusion coefficient everything is just fine - they are forming nice white puffs, gradually fading to gray. No spontaneous creation of velocity.

I can turn off the velocity update, it's still working fine since there's no velocity in the system and the dye diffusion step can't produce any velocity. My further guess is with no active velocity there's no point to advect the dye properties further in time through the velocity field (it's zero), so I turn off the advection part in the density update. But this backfires on the simulation, the dye won't diffuse at all, it's frozen - and I can't figure it out why is that.

Thank you for your input on the issue!

##### Share on other sites
Quote:
 Original post by TooropMaybe I can't describe the situation in an adequate way or I miss your point.

It's not your fault. It's been a long time since I went through this paper, and I'm having trouble explaining it. Thanks for your patience.

Still, my hunch is that the diffusion process creates some velocity.

If you look at Figure 1 in Stam's paper, you'll see on the right hand side that the equations have a wavelike (second-order derivative) term:

e.g.,
$+ \nu \nabla^2 \rho$

Now see the second equation in this wikipedia subsection called "Incompressible flow of Newtonian fluids":
http://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Incompressible_flow_of_Newtonian_fluids

You'll find that the wavelike term is related to viscosity. Wikipedia uses the mu symbol for the viscosity coefficient. Stam uses a = N*N*dt*diff;

Diffusion implies wave behaviour. A wave implies that something is travelling outward in a radially symmetric fashion. Travel implies velocity.

Basically, my hunch is that any differences in the density field between steps that is caused by the diffuse call is self-generating some velocity, if you don't comment out the advect step. Familiarity with viscosity (the diffusion of momentum) is why.

Sorry I can't be of much more help.

##### Share on other sites
I think you get the diffusion terms (kinematic viscosity and dye (density) diffusion) wrong.

velocity update - The Navier-Stokes equation describes how we should update the velocity field in the fluid - GDC03.pdf Figure 1 (top) with nu (=kinematic viscosity) + the projection step.

dye (density) update - Advection diffusion type equation describes how to move substances (dye) through the fluid - GDC03.pdf Figure 1 (bottom) with kappa (=dye diffusion).
In my opinion this step cannot generate any velocity. On the other hand a temperature scalar field would generate buoyant forces.

I would like to create a case with zero velocity field. So I skip the velocity update entirely. Since there's no velocity I thought I can leave the advect term in the density update as well. But as I mentioned earlier something is wrong and I can't figure it out from theoretical viewpoint why can't I neglect this advection term.

Thank you for not giving up!

##### Share on other sites
Quote:
 Original post by TooropI think you get the diffusion terms (kinematic viscosity and dye (density) diffusion) wrong.

That's kind of true, though diffusion of momentum and diffusion of ink are both wavelike phenomena following identical rules.

A fast moving ball loses speed in honey faster than water because the viscosity coefficient is higher.

A high ink diffusion rate equates to fast dispersal of ink, which is equivalent to high viscosity in the previous scenario of ball/honey.

The first scenario involves an expanding cloud of kinetic energy (the ball's momentum). The second scenario involves an expanding cloud of mass-energy (ink).

If you figure it out within the next couple of days, please let me know. I'll have some time soon to look into the code more thoroughly. This is now officially driving me nuts. :)

##### Share on other sites
I've been tracking this thread for three days but did not have the time to write. A month ago I tried to implement the simulation in C# actually using Stam's paper and code. I had similar errors and "gave up" then. Yesterday I looked into it again.

I'm glad Stam provided a good to follow explanation and working code. The latter is a rarity with papers. But the code is so concise it's hard to follow. So my reverse engineering and using a different representation and catchier naming is not complete, so buggy.

I agree with Turop that you can turn on and off the various stages, but to keep it stable you have to make sure the conservation laws still work (at least, probably the boundary conditions are important, too). Can't say what's going wrong in your implementation, Toorop, but I narrowed down my bug(s) by explicitly dumping everything after each step: In my case something is going wrong with the Swap() logic of the grid values, some states freeze, no matter how I add density or forces.

@taby: I can more or less follow your explanation and I agree that diffusion is an effect of moving masses: But this happens on a microscopic level, and in this simulation is "abstracted away" with a time-varying density field. If you turn off advection there won't be velocities at all.

I hope this helps a bit. Maybe I find time to get it running, I promise to report back. And thanks for the link to the other paper. It would be so cool to get it stable and O(n). ("Only" have to learn yet another math trick. Never heard of method of characteristics before)

unbird

##### Share on other sites
Quote:
 Original post by unbird@taby: I can more or less follow your explanation and I agree that diffusion is an effect of moving masses: But this happens on a microscopic level, and in this simulation is "abstracted away" with a time-varying density field. If you turn off advection there won't be velocities at all.

I agree entirely. I'm just about done my project and then I'll look through the code to confirm this. :)

##### Share on other sites
(... several facepalms later...)

I started over and now use an almost 1:1 copy for my C# implementation (it really was a bad idea to switch language and refactor at the same time). Even then I encountered similar problems - until I found my stupid bugs. Now it seems to work, so on your side, Toorop, I definitively suspect either a quirk with the SWAP or something going wrong with handling your input (sources) rather than some physics going wrong. I can turn off/on the advection like you did - and/or even the whole velocity step - I get the same results as long as there are no forces added: Here I "drop" density at a random point every 4th frame (read left-to-right, top-to-bottom):

Can you post your whole (relevant) code, please ?

##### Share on other sites
Quote:
 Original post by unbird(... several facepalms later...)

Thanks for letting us know. I definitely need to refresh my memory! I should be facepalming, not you.

##### Share on other sites
Finally, I have my own working version of Stam's code. I have messed up the pointers in my implementation of SWAP - unbird was right! :)

The code is nothing special, just some minor name convection modifications and stuff.

I haven't implemented the velocity update part, so I'm solving just the advection diffusion equation for the non-reactive substance (dye). The problem is still here - with zero velocity field I cannot skip the advection for dye. It's a misery because advection won't do anything with zero velocities.

main loop
void onIdle( ){	updateFromUI(dyeDensity_prev, velU_prev, velV_prev);	dyeUpdate(dyeDensity, dyeDensity_prev, velU, velV, dyeDiffRate);	glutPostRedisplay();}

advection diffusion equation solver
void dyeUpdate(float* x, float* x0, float* u, float* v, float diff){	addSource(x, x0);	swapBuffers(x0, x); diffuse(0, x, x0, diff);	swapBuffers(x0, x); advect(0, x, x0, u, v);	// won't work without this but IMO it should be fine}

In the diffuse() and advect() functions I haven't changed anything - it's Stam's code.

##### Share on other sites
Glad to hear it's working, sort of. Since my facepalming wasn't over yet here some pitfall you might want to look into. The naming of the globals with a _prev suffix suggest they hold the previous state, but this is only true after the update (maybe not even that). Before calling vel_step/dens_step they have to be initialized with delta values through get_from_UI - they will be added through an add_source-call at the very beginning of both parts. If no interaction is desired, the _prevs have to be zeroed, and this is where I was failing: I inserted my density/forces punctually, forgot to zero. So, the get_from_UI call is a rather important part of the simulation loop.

##### Share on other sites
Little update here!
I have extended the code to support different dye colors in order to create some nice mixing effects. The below picture depicts the evolution of four smoke puffs with different colors. I add additional forces with the mouse. Pretty nice looking! ;)

I implemented a basic dissipate routine to avoid the case of filling the whole simulation area.

But unfortunately, the initial question still persists and it's quite frustrating ...

##### Share on other sites
Niiiiiice [cool]

I can imagine your frustration, that's what I went through yesterday. You can still post your code (or PM me), maybe I spot the bug.

Thought about dissipation, too, actually had once such a bug where my density exploded without adding mass. So I provided a function which tells me the current total mass to check for preservation.

Though I now understand the maths better, I probably won't go much further with refactoring or optimizations now, just did some minor C#-obvious changes. I mainly wanted to have a CPU-reference: This thing has to go shaders!

Thanks again to make me look into this!

unbird

##### Share on other sites
I'm thankful for all your efforts - I'll send the code via a PM.

One can find a nice follow-up to the Stable Fluids method here.

##### Share on other sites
Thanks to the helpful remarks from unbird I managed to get the dyeUpdate work with zero velocity field without the advection step.

In the original implementation the dyeUpdate function gets the buffers by value. The double swapBuffer call refreshes the original buffers correctly. But if only just one swapBuffer call is present the data isn't propagated back correctly. Using pass by reference in the dyeUpdate solves the whole misery!

It was a pretty lame mistake, thanks for all your comments!

##### Share on other sites
Hi, at the end of that Stam paper, he states he's implemented the algorithm in 16-bit fixed point. Is this code available someplace? Anyone know?

##### Share on other sites
WTH - marching cubism, Batman, this little thing is a patented algorithm?! Then why did he give everyone the algo? Is this Stam guy a sadist? Why should I much care if I can't use it? mostly waste of time then. I like experimenting and working on algorithms, but code is meant to be code, not theory.