## Introduction

This article takes a look at Jos Stams Stable Fluids, the current standard for fluid simulation in games and beyond, and proposes a less complicated but also less scientific oriented approach.

## Background

What I was looking for was a paint stirring simulation. I came across Jos Stam's Stable Fluids. Its implementation does look like stirred fluids with high viscosity. It is a very popular method for fluid animation and is supposed to be scientifically sound as well, solving the Navier-Stokes equations. Alas, I do not understand it. Originally, it was presented as a method for fast fluid simulation for games, but many scientists make mention of it now. The first code was about injecting smoke into a liquid tank, the liquid an incompressable fluid. And I have three questions here.

- If I inject smoke into a liquid, it will bubble up, not spread like when I stir 2 colors of paint.
- If the tank is filled to the brim with an incompressible liquid, adding volume to it will make the tank flow over or explode.
- The implementation of advection. Advection generally means transport of mass, heat, momentum. The specific form here is a backward tracking.

What follows here is my simple, basic, experiments with fluid flow simulation and backward tracking advection. I do not have much knowledge of fluid mechanics. My implementation has no scientific value whatsoever. It is merely an effort to simulate fluid flow.

This is my simulation of the Kelvin-Helmholtz instability:

Private Sub InitiateHelm() pi16 = Pi / 48 For tX = 0 To 511 qy1 = Sin(tX * pi16) * tX * 0.01 For tY = 0 To 255 Vel(tX, tY, 1) = 1 ''velocity x Vel(tX, tY, 2) = qy1 ''velocity y For tX2 = 0 To 2 ImageData(tX2, tX, tY) = 255 ''white top half of the picture Next tX2 Intfield(tX, tY, 1) = 255 Next tY For tY = 256 To 511 Vel(tX, tY, 1) = -1 Vel(tX, tY, 2) = qy1 For tX2 = 0 To 2 ImageData(tX2, tX, tY) = 0 ''black bottom half of the picture Next tX2 Intfield(tX, tY, 1) = 0 Next tY Next tX Private Sub Draw() Dim a1 As Double, b1 As Double, qx1 As Double, qy1 As Double For tX = 0 To 510 For tY = 0 To 510 qx1 = tX + Vel(tX, tY, 1) qy1 = tY + Vel(tX, tY, 2) tX1 = Int(qx1) tY1 = Int(qy1) qx1 = qx1 - tX1 qy1 = qy1 - tY1 If tX1 < 0 Then tX1 = 0 If tX1 > 510 Then tX1 = 510 If tY1 < 0 Then tY1 = 0 If tY1 > 510 Then tY1 = 510 ''a1 = lerp(CDbl(Intfield(tX1, tY1, 1)), ''CDbl(Intfield(tX1 + 1, tY1, 1)), qx1) ''interpolation for smoother picture ''b1 = lerp(CDbl(Intfield(tX1, tY1 + 1, 1)), CDbl(Intfield(tX1 + 1, tY1 + 1, 1)), qx1) ''Intfield(tX, tY, 2) = lerp(a1, b1, qy1) Intfield(tX, tY, 2) = Intfield(tX1, tY1, 1) ''backward tracking advection ''Intfield(tX1, tY1, 2) = Intfield(tX, tY, 1) ''plain forward Next tY Next tX For tX = 0 To 511 For tY = 0 To 511 Intfield(tX, tY, 1) = (Intfield(tX , tY, 2) For tX2 = 0 to 2 ImageData(tX2, tX, tY) = Intfield(tX, tY, 1) Next tx2 Next tY Next tX End Sub

## Using the Code

It starts with a white field above a black, and I assign velocities. Not a surprise that the top velocity is opposite to the bottom one. But I also assume a periodical, sinusoid velocity perpendicular (y-direction) to these x-direction velocities. That is not supported by science. It was just something that I noticed in many pictures and videos of real fluid movement: Rayleigh-Taylor instability (looks like a variation to Kelvin-HelmHoltz), ink in water (looks like a form of Rayleigh-Taylor), and incense smoke columns. Some of these do not look like the Stable Fluid simulations. It looks as if there is a periodic velocity as well. A bit like low frequency sound waves, with a velocity straight ahead in horizontal direction accompanied by periodic amplitudes in vertical direction. It seems like something I can use in my fluid simulation, without explaining how it was caused, just assuming it.

I originally used this line to assign a value to every pixel:

Intfield(tX, tY, 2) = Intfield(tX1, tY1, 1) ''backward tracking advection

And the result is a curly wave. If I use the straightforward (literally) variation:

Intfield(tX1, tY1, 2) = Intfield(tX, tY, 1)

I get no curls! The curl in the picture is caused by the backward advection.

What happens is that a black pixel occasionally taps into the velocity of the white field, and vice versa.

It did not really come from there but the backward calculation makes it look that way. Backward advection causes all those nice curls. (I have seen many implementations where extra curl or vorticity is added. That does make even nicer pictures.)

As I said, I was originally looking for a simulation of stirring 2 colors of paint. It can be simulated by warping. That is the Liquify tool in Photoshop. Here is my implementation:

Private Sub InitWarp ReDim cSeg1(-rS To rS, -rS To rS) As Single ''rS is the radius of the warp brush For tX = 0 To rS dax = tX * tX For tY = 0 To rS dxy = Sqr(dax + tY * tY) If dxy <= rS Then dr = LineArr(Int(dxy)) ''LineArr stores the warp distances in the brush, ''larger distance of warp = closer to the middle of the brush cSeg1(tX, tY) = dr cSeg1(-tX, tY) = dr cSeg1(tX, -tY) = dr cSeg1(-tX, -tY) = dr End If Next tY Next tX End Sub Private Sub DoWarp dx = xM - x ''xM = middle(old mouse position x), x = new mouse position x dy = yM - y qX5 = pBox1.Width ''Actual picturebox qY5 = pBox1.Height For qY = -rS To rS qY1 = qY + yM If qY1 > -1 And qY1 < qY5 Then For qX = -rS To rS qX1 = qX + xM If qX1 > -1 And qX1 < qX5 Then ShiftArr(qX1, qY1, 1) = ShiftArr(qX1, qY1, 1) + cSeg1(qX, qY) * dx ''Shiftarr is a ''buffer that stores the warp distances of the brush, applied to the picture ShiftArr(qX1, qY1, 2) = ShiftArr(qX1, qY1, 2) + cSeg1(qX, qY) * dy End If Next qX End If Next qY For qY = -rS To rS qY1 = qY + yM If qY1 > -1 And qY1 < qY5 Then For qX = -rS To rS qX1 = qX + xM If qX1 > -1 And qX1 < qX5 Then qX3 = qX1 + ShiftArr(qX1, qY1, 1) qY3 = qY1 + ShiftArr(qX1, qY1, 2) If qX3 > -1 And qX3 < qX5 And qY3 > -1 And qY3 < qY5 Then For t = 0 To 2 ImageData(t, qX1, qY1) = ImageData2(t, qX3, qY3) ''Imagedata2 is the picture, ''Imagedata is a copy that will replace Imagedata2 once all changes are drawn Next t End If End If Next qX End If Next qY End Sub

Then when I looked at soap bubbles, the movements of the colored soap film, it seemed like an all-over warp. I implemented it like this:

Private Sub InitRand(b1 As Integer, b2 As Integer) Randomize Dim Randi(0 To 511) As Integer ''random warp distance across the picture tY1 = 255 tX1 = CInt(Rnd * 2) If tX1 < 1 Then tX1 = -1 Randi(0) = tX1 For tX = 1 To 511 ''y direction, random decide to add or distract 1 ''from the previous warp distance qy1 = Rnd If qy1 < 0.167 Then tX1 = Abs(tX1) - 1 If qy1 > 0.833 Then tX1 = -Abs(tX1) + 1 If Randi(tX - 1) < -b1 Then tX1 = 1 If Randi(tX - 1) > b1 Then tX1 = -1 Randi(tX) = Randi(tX - 1) + tX1 Next tX For tX = 3 To 508 ''smoothen the random values a bit Vel(tX, 256, 2) = (Randi(tX) + Randi(tX - 1) + Randi(tX + 1) + _ Randi(tX - 2) + Randi(tX + 2) + Randi(tX - 3) + Randi(tX + 3)) / 7 Next tX For tX = 0 To 511 qx1 = Vel(tX, 256, 2) / 256 For tY = 0 To 255 Vel(tX, tY, 2) = tY * qx1 ''warp distance at the edges = 0, '' calculate the warp (velocity) distance for all pixels, y direction Next tY qx1 = -qx1 For tY = 1 To 255 Vel(tX, tY + 256, 2) = Vel(tX, 256, 2) + tY * qx1 Next tY qy1 = Vel(tX, 511, 2) Next tX tY1 = 255 tX1 = CInt(Rnd * 2) If tX1 < 1 Then tX1 = -1 Randi(0) = tX1 For tX = 1 To 511 '' same for x direction qy1 = Rnd If qy1 < 0.167 Then tX1 = Abs(tX1) - 1 If qy1 > 0.833 Then tX1 = -Abs(tX1) + 1 If Randi(tX - 1) < -b2 Then tX1 = 1 If Randi(tX - 1) > b2 Then tX1 = -1 Randi(tX) = Randi(tX - 1) + tX1 Next tX For tX = 3 To 508 Vel(256, tX, 1) = (Randi(tX) + Randi(tX - 1) + Randi(tX + 1) + _ Randi(tX - 2) + Randi(tX + 2) + Randi(tX - 3) + Randi(tX + 3)) / 7 Next tX For tX = 0 To 511 qx1 = Vel(256, tX, 1) / 256 For tY = 0 To 255 Vel(tY, tX, 1) = tY * qx1 Next tY qx1 = -qx1 For tY = 1 To 255 Vel(tY + 256, tX, 1) = Vel(256, tX, 1) + tY * qx1 Next tY Next tX End Sub

The drawing routine is the same as above. Again, the backward advection creates curls and fingers.

## Further Reading

- Presentation by Jos Stam https://www.ljll.math.upmc.fr/~frey/cours/references/Stam%20J.,%20Stable%20fluids.pdf

Many more on the subject, including source codes for various languages if you type "Jos Stam Stable Fluids" in your search engine. - Kevin-Helmholtz instability: https://en.wikipedia.org/wiki/Kelvin%E2%80%93Helmholtz_instability
- Rayleigh-Taylor instability: https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Taylor_instability
- Ink in water: https://www.youtube.com/watch?v=t-Are2N7zUE
- Incense smoke: https://photodune.net/item/abstract-shapes-formed-by-scrolls-of-incense-smoke/21529834 (very visible are the periodic rings that start at the bottom of the smoke column as increase and decrease of the smoke density, and become curls higher up the smoke column)

## Conclusion and Points of Interest

Decide for yourself if this very simple algorithm creates realistic fluid simulations! It was coded in VB6, but by using C++ and a GPU implementation it can be much faster. I hope to inspire you to experiment more. As for myself, next I want to try to get better simulations of smoke columns and ink in water, as well as soap bubble films. Suggestions are welcome!

## History

- 10
^{th}November, 2020: Initial version