Screenshots and Code
Fluid simulation looks pretty cool, but there are no web pages (well, I couldn't find any) that describe in basic terms how it's done. So, after writing a mostly working implementation of Stable Fluids, I thought I'd make one.
An enclosed area, or volume, of fluid can be represented as two fields: the velocity field and the density field.
Each element (voxel) of the velocity field has an x and y velocity stored at the right and bottom face, respectively, and the density field is composed of cells with velocities defined at the center.
In each simulation step, the density is advected based on the corresponding areas in the velocity grid. The change in the velocity grid is described by the equations above.
The first part of (2) is advection. This involves two grids: the old grid and the new grid, which is initially a copy of the old one. For each voxel in the new grid, one advects backwards by dt * v to an old voxel, where v is the velocity at the voxel in the old grid. The interpolated velocity at this voxel is assigned to the voxel in the new grid.
Note: One can advect backwards n times each step, moving dt/n time each small step, for a more accurate representation of the flow. I believe this is called a Runge-Kutta method.
The second component is viscosity.
The equation for solving, implicitly, for the effects of viscosity is:
where nu and t are constants, I is the identity matrix, ∇2 is the vector laplacian, and w3(x) and w2(x) are the velocity fields (in column form) after and before the viscosity calculation, respectively.
The vector laplacian, when multiplied by an element in our velocity, gives a vector of the second derivative of each component in each direction. Since we don't have a field defined by an analytic function, finite differencing must be used.
The final expression for the second derivative at b is: f''(b)=f(a)-2f(b)+f(c). It is clear that these values "line up"; i.e. the differences in total velocity falls at a face.
The same technique is used to find the second derivatives at points in the field; 1,-2, and 1 end up being the coefficients in the rows of the vector laplacian matrix. Then any number of methods can be used to solve the linear system; I chose conjugate gradient, but I don't have particularly strong reasons for it. If you know of a better way, please email me.
Here it is most useful to realize that the values line up. In the first equation, the divergence of the velocity field should fall at the center of a voxel; similarly, the second derivative of a scalar field is also defined at the center. In the second equation, the gradient of the scalar field approximates best at the face of the voxels; coincidentally, that's where the velocities are defined.
The pressure solver works much the same way. The laplacian here is a scalar one; it is used to calculate an intermediate field, q, whose gradient is subtracted from the old velocity field. The ∇ · represents the divergence at each velocity field voxel.
Update: So, the simulator doesn't completely work; it loses mass as the simulation progresses. After plotting divergence and density, I realized that mass is lost at the edges. I tried simulating only the inner edges, but it still failed :(. If you have any ideas about what else to try, please email me.