Stable Fluids

For my final project in Scientific Visualization, a class taught by my advisor Han-Wei Shen in Spring 2022, I experiment with an implementation of the paper Stable Fluids by Jos Stam. My code is available here, and can be ran with "python Code/main.py" after downloading. For this project, I accomodate an existing Stable Fluids solution in Python by adding functionality for the experiments performed, such as boundaries, custom ink setups, and moving sources. From the other Stable Fluids implementation, I use the fluid solver (fluid.py), which is responsible for a single timestep advection, and numerical.py, which creates a few math operations used by the fluid solver. However, adjustments are made to the to the initialization and step in fluid.py to accomodate boundaries, which I added. Main.py is written by me for the examples shown at the end, each of which are not able to be done with the original code.

Fluid like behaviour is a very fascinating topic, especially because we still haven't found exact equations to model the fluids themselves. However, we do have approximations or estimates of equations that determine fluid motion. The Navier-Stokes equations are one such model for fluid motion that do well, for which we do not have an analytical solution yet. If you have a solution to the Navier-Stokes equations, you would win $1,000,000 since it is a Millenium Prize Problem! Therefore, to solve nonlinear differential equation, different approximations and discritized solutions exist, including the one this project focuses on called "Stable Fluids".

This implementation is not as rigorous as other solutions to the Navier Stokes equations, but gives a good enough visualization, especially for cinematic reasons within games and movies. The solutions using this paper's implementation are more "diffused" than other Eulerian solutions. As a benefit, the simulation speed is extremely fast, and it is guaranteed to not "blow up" as other solutions will under certain conditions. On the engineering side, it is also much easier to both understand and implement, and lends itself to parallelization nicely.

Above is the form of the Naiver-Stokes equations that are used in this paper. Stable Fluids uses a semi-Lagrangian and implicit approach, broken into four steps. For this project, only three of the steps are used, with the "diffusion" step being skipped. Diffusion is skipped because viscosity tends to dampen the fluid simulation, and inviscous fluids look cooler. Therefore, in our simulations, we have 0 viscosity, and can remove this diffusion term altogether.

Above shows the 4 steps within a single update to the vector field. "Add force" is the f term in the Navier Stokes equation, which represents external forces. This can be used to create artificial boundaries, sources, and sinks, as well as gravity.

The next step is the advection, which is extremely intuitive in this paper. All that is needed is a particle tracer, which in this class, we did for Lab 3! In the advection step, the velocity values themselves are advected to the next step using the current velocity field. However, if a particle is dropped at each grid point and advected forward, the interpolation becomes challenging because final velocity values may land off of grid points, shown below.

Instead, the author advects the velocity field backward by the same amount of time, and then uses linear interpolation to get the final velocity value for the grid point advection started at. The process is shown below.

Since we are skipping over diffusion, the last step is to project. In the previous steps, we may have broken our incompressible rule. This step is responsible for fixing those problems. The author uses two solutions given two boundary conditions: a general solution and an elegant solution for periodic boundaries.

The general solution for the projection step is solving a Poisson problem, requiring a Poisson solver. In essence, the Poisson equation comes down to solving a large, sparse, linear system. The original paper uses a multigrid method for the solution, which is linear time, just like the advection step above.

However, when the boundaries are periodic, the solution is more elegant. The periodicity allows the velocity field to be transformed into the Fourier domain. In the Fourier domain, the vector field can easily be decomposed into its divergence free and irrotational components. The author creates an operator P which is responsible for taking only the divergence free component of the vector field in the Fourier domain. After taking that component, it can be inverse transformed to the original velocity space, and the update is complete.

Above on the left shows the operator that takes only the divergence free component of the vector field w. k is the wavenumber, and by multiplying k by w(k), the result is orthogonal, which gives a divergence free result. For more information, I suggest a conference talk here.

Above on the right shows the entire algorithm. We didn't explain the diffusion step, as we assume the fluid is inviscous, so v=0. These steps are very efficient and do not require timestep refinement, and are unconditionally stable.

To visualize results, we use the same advection calculation we use for the velocity field advection step. We can inject dye into the fluid at arbitrary points with arbitrary color, and keep track of them as they're advected in the fluid over time. Below are some visualizations made.

All experiments are done using Python with NumPy and SciPy. The boundaries are not periodic in my examples, so a laplacian is solved using SciPy instead of the Fourier domain solution. Grid sizes are 512^2 for video 1 and 2, and 256^2 for video 3. Video 1 has 500 timesteps, and video 2 and 3 have 1500. Time per frame is between 10ms and 400 ms, which can be realtime if lower framerates are allowed.

In videos 1 and 2, I add boundaries such that fluid cannot escape the domain. However in the third video, dye is permitted to leave the domain since results were odd when constrained.

In the first video, the external forces are used to stop ink from entering the letters "SCIVIS" inside the domain. The ink is colored based on vorticity, and is injected constantly from the source spots, where the vector field is set to point inward.

In the second video, the ink is initialized as a previous photo taken from a scivis group dinner in December 2021. The ink is mixed throughout the domain from three sources at the bottom of the image. Toward the end, it just looks like coffee being mixed!

The last video shows when two colors are mixed together. To make the visualization, I "overwrite" the output image with a white circle where my "mixing" stick is. I have coded it to mix clockwise by adding external forces at the positions folowing the movement of the mixing stick.

Overall quality is good, but the slow spread of ink is apparent, especially in the beginning of video 2, where the ink for the image quickly blurs.

In the future, I may implement this in Javascript on my website so a realtime tool could be used, kind of like this tool. That implemenation uses WebGL, which would be more performant than javascript, however.