mikeash.com pyblog/fluid-simulation-for-dummies.html comments

SPC - 2026-06-03 18:46:56

Wed, 03 Jun 2026 18:46:56 GMT

A revision of the code, including bug fixes and documentation, is available here:
pastebin[dot]com/p4kZHfuy

Fergal - 2026-06-03 18:07:20

Wed, 03 Jun 2026 18:07:20 GMT

For linsolve, it looks like you're solving the Laplace equation numerically:

You may find the Wikipedia article on
Laplace_equation_for_irrotational_flow

helpful

Catto - 2026-04-02 16:49:06

Thu, 02 Apr 2026 16:49:06 GMT

niggerlicious

Where U Elevate - 2026-01-08 09:18:50

Thu, 08 Jan 2026 09:18:50 GMT

Career guidance for freshers is crucial in today’s competitive market, and platforms like Where U Elevate help students discover the right career paths through startup exposure, internships, and real-world opportunities.

whereuelevate.com/blogs/career-guidance-for-freshers

Omar - 2024-07-06 17:37:01

Sat, 06 Jul 2024 17:37:01 GMT

Great to see an actual implementation after reading all the thoery, thank you for this. But seriously, what is going on with these variable names? Can anyone with a closer understanding of the code please explain what all the variables in the advect function do?

e - 2022-04-26 15:47:23

Tue, 26 Apr 2022 15:47:23 GMT

water

bentheone - 2022-02-15 09:04:28

Tue, 15 Feb 2022 09:04:28 GMT

@Anssi, yes ! Thank you, I can't believe nobody pointed that out before. The diffusion ratio being dependent on the grid size makes ZERO sense.

Anssi - 2022-01-10 20:22:55

Mon, 10 Jan 2022 20:22:55 GMT

I think there are many intances in the code (also in Stam's original article), where 'N', or 'N-2' should be replace by '1/SCALE'.

This because SCALE represents the spatial step size, not 1/(N-2). Right?

For example, the variable 'a' in the diffuse function should be defined as
a=dt*diff*SCALE^{-2}, instead of a=dt*diff*(N-2)^2.

Similar corrections seem to be in order in 'advect' and 'project' functions.

At least this works for me. Now the qualitative local behaviour does not depend on the size of the grid N, as it should not. You need to adjust the parameters again, if you had adjusted them according to the original code.

Btw. In diffuse function the 'a' is usually less than 1, and we can replace the 'lin_solve' there with the faster (and more straightforward)
x[IX(i, j)] = x0[IX(i, j)] + a*( x0[IX(i+1, j)]+x0[IX(i-1, j)] +x0[IX(i, j+1)]+x0[IX(i, j-1)]-4*x0[IX(i, j)])
which would be stable in this case. (This was considered in Stam's article but discarded unnecessarily in my mind.)

Aadhavan P - 2021-09-04 16:26:13

Sat, 04 Sep 2021 16:26:13 GMT

Hi!
Very good article. I'm saving this for later. Will try to implement it after my exams. I'm very interested in simulating stuff like these.

Polaris - 2021-03-03 20:27:03

Wed, 03 Mar 2021 20:27:03 GMT

Hey are you sure this is for dummies? As a dummy myself, I regret to inform that I was unable to understand.

rich engle - 2021-01-14 03:07:01

Thu, 14 Jan 2021 03:07:01 GMT

is there any way to have 2 fluids, like water and air?

Dimitri - 2020-10-26 20:05:42

Mon, 26 Oct 2020 20:05:42 GMT

Hey, great work!
Truly life-saving when this is not your field of work or studies, feels good to see the implementation beyond the theory
Thanks a lot for posting it

Eli Maynard - 2019-03-24 13:48:08

Sun, 24 Mar 2019 13:48:08 GMT

Hi, in



    x[IX(0, 0, N-1)]     = 0.33f * (x[IX(1, 0, N-1)]

                                  + x[IX(0, 1, N-1)]

                                  + x[IX(0, 0, N)]);

should the IX(0, 0, N) be IX(0, 0, N-2)?

Thanks for the article!

Sam - 2017-12-07 04:33:33

Thu, 07 Dec 2017 04:33:33 GMT

I haven't read the whole thing yet, but the pointers are confusing me a bit. Why initialize a pointer rather than a value?

Scott - 2015-08-23 23:47:11

Sun, 23 Aug 2015 23:47:11 GMT

Mike, Do you remember why you need to multiply 'a' by '(N-2)^2 ' in the diffuse function? For example Stam does this in 2D, so I'm expecting it to be '(N-2)^3' for 3D?

EWD - 2014-05-31 18:42:54

Sat, 31 May 2014 18:42:54 GMT

Can anybody explain to me why in 2006 people still use variable names such as "b", "N" or "p"??? Are they trying to save the compiler the work with reading long words?? Were the variables called 'shit' and 'fuck', it would still be more helpful!

But reading some of the comments here, some people deserve it.

bob: Yeah, N*N*N elements for a 3-dimensional array is a terrible waste. Try to allocate N elements for a 3-dimensional array.

bob - 2011-07-19 18:24:17

Tue, 19 Jul 2011 18:24:17 GMT

Why do you allocate N*N*N size for each dimension? This becomes impractical when work with large amounts of particles. For example,
sizeof(float) * 8 callocs * (1000 particles ^ 3) = 32 gb.

Martin - 2011-03-29 15:48:51

Tue, 29 Mar 2011 15:48:51 GMT

I had a deep look at the advect function and ask myself if it's correct to call it with the N-value. Shouldn't this be N-2 (or Nflaot be initialized by N-2)?

Anyway, thanks a lot for the code, it's really valuable for learning how this stuff works!

libero - 2011-02-04 21:45:29

Fri, 04 Feb 2011 21:45:29 GMT

does anyone have a working main() for this code because seems to be crashing and I cant find the reason for that? any help you be appreciated.

Akash Verma - 2010-10-23 21:31:08

Sat, 23 Oct 2010 21:31:08 GMT

Could someone please put up a main() for this to see it working?

Aaron Moffatt - 2010-06-23 09:01:50

Wed, 23 Jun 2010 09:01:50 GMT

Very, very useful document. Just to the point, outlining and giving example code of the key concepts. Really appreciate you writing and posting this.

uma sankar pradhan - 2010-04-02 22:13:04

Fri, 02 Apr 2010 22:13:04 GMT

the code is great.
But in the projection method,
during the calculation of divergence of velocity field,why your dividing it by N.I suppose ,it shoud be multiplied as dx=dy=dz=1/N
div[IX(i, j, k)] = -0.5f*(
                         velocX[IX(i+1, j , k )]
                        -velocX[IX(i-1, j , k )]
                        +velocY[IX(i , j+1, k )]
                        -velocY[IX(i , j-1, k )]
                        +velocZ[IX(i , j , k+1)]
                        -velocZ[IX(i , j , k-1)]
                    )*N;
instead of
div[IX(i, j, k)] = -0.5f*(
                         velocX[IX(i+1, j , k )]
                        -velocX[IX(i-1, j , k )]
                        +velocY[IX(i , j+1, k )]
                        -velocY[IX(i , j-1, k )]
                        +velocZ[IX(i , j , k+1)]
                        -velocZ[IX(i , j , k-1)]
                    )/N;
One more thing,
the Vx0,Vy0 need to be initialised to zero before each iteration.
I may be wrong.If i am wrong,kindly inform me.

fluid - 2009-10-09 10:33:03

Fri, 09 Oct 2009 10:33:03 GMT

really cool fluid simulation is here: http://www.escapemotions.com/experiments/fluid_fire_3/index.html

Kent - 2009-07-23 15:38:22

Thu, 23 Jul 2009 15:38:22 GMT

Thanks heaps for posting this code. I was just about to start moving Stams code from 2D to 3D myself. And I really like the way you have presented the code here. Nice work!

willem - 2009-05-14 06:34:43

Thu, 14 May 2009 06:34:43 GMT

great to see it written for a non-mathematician. thanks

student - 2009-05-10 06:36:10

Sun, 10 May 2009 06:36:10 GMT

Im around this same code the 3d code you provide was very helpfull thanks.
Im hava used several linear solvers gauss seidel relaxation (seq and parallel - Red Black), jacobi (seq and parallel) ....
I'm trying to implement the CG method also because I adjusted the stams code to allow internal boundaries (moving or not), among other things.
I have a pertinent question: I understood after considereble digging that we are solving a system of linear equations in the form Ax=b however what confuses me is what is A exactly?

Any comments or regards on this would be quitte helpfull anyone?

responses to: a14722@ubi.pt

mikeash - 2009-04-14 02:33:14

Tue, 14 Apr 2009 02:33:14 GMT

There is no main function because it's a discussion of techniques, not a full working example. You will have to incorporate what is provided into your own code. If you know little C then an application of this complexity is probably not a good place to start. You need to learn to walk before you can run.

Tom - 2009-04-14 02:15:33

Tue, 14 Apr 2009 02:15:33 GMT

Why is there no main function?
actually i know little about c, but i wanna learn this application and apply it into processing.

nabeel - 2009-04-03 16:20:54

Fri, 03 Apr 2009 16:20:54 GMT

Great stuff. Wonder why they never make sound this easy in texts. wonder if i can work out mixing in two phases with this???

Lars - 2008-12-17 05:12:52

Wed, 17 Dec 2008 05:12:52 GMT

Thanks very much for the 3D conversion of Jos Stam's 2D code -- it has saved me work and brainache. :-)
I was hoping to make a fun screensaver and this should give me a good boost.

Ivan DeWolf - 2008-09-17 10:25:18

Wed, 17 Sep 2008 10:25:18 GMT

here is another website with a really legible java implementation in 2d... be sure to look at the source code...
it also uses a gauss-seidell relaxation scheme for it's linear solver. I really wish I understood the linear solver part more.

http://www.multires.caltech.edu/teaching/demos/java/stablefluids.htm

Setis - 2008-06-16 20:40:59

Mon, 16 Jun 2008 20:40:59 GMT

About the rendering:

Since you know how to take a 2d slice, showing the density of ink, you can render it in 3d this way:

Take n 2d slices in a given time. Draw the slices in textures, so the more ink in a given cell, the less transparent and more colored is the corresponding pixel. Now draw a 3D cube of all the slices stacking. That way, you can see the ink moving in the cube. Just add some code for rotating the cube and so.

anonymous friend - 2008-01-25 03:07:26

Fri, 25 Jan 2008 03:07:26 GMT

The lin_solve function simply solve a system of linear equations. Since the matrix of this system is sparse (i.e contain a lot of zeros), this system can be solved efficiently (note that there exist more accurate algorithms) using a Gauss-Seidel algorithm.

What do we solve exactly ? For the diffusion, we try to find tthe densities which, when diffused backward in time, gives the density we started with. In the case of the projection, it is a little more complicated. We are solving a linear system called a poisson equation. In mathematics, Poisson's equation is a partial differential equation. To make sure that the fluid is incompressible, we need to satisfy the non-divergence condition.

The divergence of a vector field (the velicity of the fluid)at a particular location is a scalar value that essentially represents how much the field is ‘expanding’ at the point. A positive value represents a divergence and a negative value represents a convergence or compression.If we use the assumption that a fluid is incompressible, then that implies that the divergence of the velocity must be zero everywhere.

A useful splitting method for solving the incompressible flow equations is the so-called ‘projection method’. In this technique, one first computes a divergent vector field that results from the advection & viscosity terms, and then projects that entire vector field onto the space of non-divergent vector fields
This can be thought of as creating an instantaneous pressure field whose gradient compensates for any divergence caused by the advection & viscosity.

Getting back to the poisson equation, we create these beast to solve the unknown pressure field. The gradient of this pressure field is an additional force on the fluid and the divergence of this gradient compensates for the divergence resulting from the convection-viscosity equation. The result is a final velocity field that has no divergence.

mikeash - 2006-12-12 23:27:00

Tue, 12 Dec 2006 23:27:00 GMT

Brian, I’m glad you found it useful, and please do send me your results if you’d like.

Brian Silva - 2006-12-02 18:48:00

Sat, 02 Dec 2006 18:48:00 GMT

Nice work! I am doing exactly the same work with Stam’s implementation, but I was having trouble getting the simulation step to function correctly. I took your methodology from the FluidCubeStep and it solved my problem. So, thank you! 
 
I agree with you about the original sources: they are difficult to understand and don’t explain the math very well. I’m not entirely mathematically-inclined so I struggled for a long time to understand it. I’m setting out on a school project to explain the math from this particular implementation in a clear way. If you’re interested, I can send it to you when I finish. 
 
Again, great work!

neil - 2006-08-16 08:54:00

Wed, 16 Aug 2006 08:54:00 GMT

im a newbei and i want to learn more about simulation and animation, but my problem is i cant afford to enroll to schools, thats why im looking for free ebooks and tutorials which is available in this world. but my collections are not enough for my thirsty head. all i want is to learn more, if anyone could share their knowledge and books or sites i dont know yet, please enlighten me, by the way im making some simulation in physics principles or topics to help me in my carrier as future physics teacher, but i still need to learn lots of things, so guys please share your informations i assure you it wont go to waste. thanks. heres my add, <a href="mailto:kaerfmove@yahoo.com, thanks!" rel="nofollow">kaerfmove@yahoo.com, thanks!</a>

Artur W - 2006-07-12 14:23:00

Wed, 12 Jul 2006 14:23:00 GMT

Thank you very much! Finally something for a dummy like myself. Less equations and more detail on turning them into actual implementation. 
This helps!

N382AF - 2006-03-15 01:01:00

Wed, 15 Mar 2006 01:01:00 GMT

Nice work here, I saved it for reading later. A bit over my head, but I love learning this stuff! 
 
Thanks!