Physics of Flow

What are we trying to do?

Let's program a graphics card to move useful information around in a way that is consistent with physical flow. We will use a grid of automata that repeatedly update their state by communicating with the automata around them. Here is an example of such a program.

What you need to know before we get started

1. The single most important fact about physics is that force is the negative gradient of energy.
2. F = ma isn't the whole truth.
The universe actually discretely exchanges momentum.
There is no such thing as a force - this is an abstract concept that makes sense to gigantic beings such as ourselves.

The universe is actually kind of like a network of automata that discretely exchange different forms of energy. I think of this as a bunch of simple merchants that act in their own self interest and buy and sell goods with the merchants around them. Nothing is created or destroyed - things can only be reordered.

But the interesting thing is, that even at scales where quantum mechanics don't apply, you can still arbitrary discretize space and time into cells and apply the same logic of the merchants trading information discretely to approximate continuous dynamics.

This is how we will simulate a fluid on a grid of automata.

The Manifold

According to loop quantum gravity, spacetime is like a foam of automata. Think of bubbles in your kitchen sink. Each bubble is directly touching some number of other bubbles. This is the bubble's neighborhood. I'm going to use the word neighborhood a lot here.

Moving forward it is very important that you have a picture in your head of an automata connected to an adjacent neighborhood of like automata. Each member of that neighborhood is connected to other neighborhoods and if you keep spreading out, you have the entire space (or manifold).

In our grid fluid, the neighborhood looks like a cross. For each cell, there are four neighbors. I like to call them :

n (north/up)
e (east/right)
s (south/down)
w (west/left)

and I like to call the cell that is being computed "me".

so "me" has "n", "e", "s", and "w" to worry about.

What does the automata need to know?

In our simulation, there are four pieces of relevant information.

x velocity
y velocity
pressure
spin/vorticity

Each cell contains four numbers that represent this information.
I like to store that information like this in the state vector of the cell :

me.xy = velocity
me.z = pressure
me.w = spin/vorticity

How do four numbers correspond to what's actually going on in water?

It's important to remember that the particles involved in water flow are way smaller than any flow features you can see and way smaller than the cell size of our simulation. So these are actually statistical values.

That means that the velocity is the average linear momentum of all the particles inside that cell.

The pressure is actually also the average linear momentum of each particle in the fluid... do you remember PV = nRT?
So temperature and pressure are macroscopic interpretations of the collisions of microscopic particles.
Here is a visualization of kinetic theory.

The vorticity is yet another average of the linear momentum of the particles in the cell... but this time it is the tendency of these momenta to go in a circle.

So all in all we have 3 different versions of kinetic energy :

xy : velocity : the tendency of the cell to move.
z : pressure : the tendency of the cell to bang into itself.
w : vorticity: the tendency of the cell to go in a circle.

Where is this information?

So the cool thing about fluids is that all of this state information is always moving and evolving. In the next section we'll talk about how it evolves. But for now let's talk about where it is.

In our fluid simulation we will need to find the state of the fluid at different locations.
Let's say we have the previous state of the fluid stored in a texture called "A".
Now let's say we want to find out what's going on at location "U" at this point in time.
If we just look up A(U) [the state of "A" at location "U"], we will not have the right information because fluids move.
Fortunately for us, fluids move predictably. So we can guess where "U" used to be in the last state of the fluid by using the velocity at point "U" in "A" [aka A(U).xy]

This is how I like to handle this in code :

vec4 sample_A (vec2 U) {
 // Make sure you're using linear texture filtering!
 return texture(A,U/resolution.xy);
}
vec4 sample_A_advected (vec2 U) {
 // here I subtract the velocity at A(U) from U to get a new position
 vec2 U_prev = U - sample_A (U).xy ;
 // U_prev is very likely to be where location "U" used to be
 return sample_A (U_prev);
}
So in order to get "me", "n", "e", "s", and "w" we would say :
me = sample_A_advected (U)
n = sample_A_advected (U + vec2(0,1))
e = sample_A_advected (U + vec2(1,0))
s = sample_A_advected (U - vec2(0,1))
w = sample_A_advected (U - vec2(1,0))

The Evolution of Flow

This is the most difficult section. Here we are going to talk about how our four values interact with their neighbors. This gets back to the part at the beginning where we talked about merchant automata. The neighborhood around our cell is going to interact with our cell. We have to make sure that these interactions are symmetrical so that nothing is created or destroyed. Different types of energy from the neighborhood will be traded with different types of energy from the cell.

The state of the fluid evolves in the following ways :
1. The velocity experiences a force against the gradient of the pressure. Pressure is energy stored inside the cell, so this is actually the very first thing I said in this article - force is the gradient of energy.

2. The velocity experiences a force perpendicular to the spin of the field. Imagine a baseball shooter : when two wheels with opposite spin come in contact with a baseball, they shoot it out perpendicular to the arrangement of the wheels.

3. Because pressure involves banging around - this value is very volatile. So the pressure in the cell is actually the average of the pressure of the neighborhood's previous state. That is to say, the pressure was completely traded with the neighborhood.

4. The pressure increases by the divergence of the field. That is to say, if the neighborhood is getting closer to "me", "me's" pressure is going to increase. If the neighborhood is moving away from "me", "me's" pressure is going to decrease.

5. The spin of the cell is similar to the pressure in that it is also completely traded with the neighborhood.

6. The spin trades with the curl of the fluid. The curl of the fluid is perpendicular to the divergence of the fluid. This is the amount that the neighborhood is going in a circle around "me"

The most important aspect of making the automata behave symmetrically is by dividing the interaction between each cell by the number of cells in the neighborhood. So with a grid there are four neighbors, so all of our equations will have a "0.25*( ... )".

The following are the discrete equations of spinor flow (that's a made up term to describe the algorithm I made up):

(order matters - don't compute z and w before x and y)
Pure algebra :
C.x -= 0.25*(e.z-w.z+(n.w*C.w-s.w*C.w));
C.y -= 0.25*(n.z-s.z+(e.w*C.w-w.w*C.w));
C.z = 0.25*((C.z,s.y-n.y+w.x-e.x)+(n.z+e.z+s.z+w.z));
C.w = 0.25*((n.x-s.x+w.y-e.y)-(n.w+e.w+s.w+w.w-0.*C.w));
English :
me velocity -= gradient(pressure) + spin interaction
me pressure = mean_pressure_of_neighborhood + divergence_of_velocity
me vorticity = curl_of_velocity - mean_vorticity_of_neighborhood

Digging deeper :

 Pressure gradient :
   gradient = vec2 (change in pressure across x, change in pressure across y);
   gradient = vec2 (e.z-w.z, n.z-e.z);

 Divergence :
  how much "n" is going up +
  how much "e" is going right +
  how much "s" is going down +
  how much "w" is going left
  divergence = n.y + e.x - s.y - w.x

 Curl :
  how much "n" is going right +
  how much "e" is going down +
  how much "s" is going left +
  how much "w" is going up
  curl = n.x - e.y - s.x + w.y

 Spin Interaction (this one is tricky)
  spin_interaction = (sign of the relativity of spins)*(spin1 * spin2)
   (in 3D, this would involve a cross product and the sign would come out of that)
  spin_force_x = n.w*me.w - s.w*me.w;
  spin_force_y = e.w*me.w - w.w*me.w;

Conclusion

The velocity of the neighborhood affects the spin and pressure of the cell.
The spin and pressure of the neighborhood affects the velocity of the cell.



In order to make a good fluid, you also have to define some boundary conditions. These might look something like this:
To make boundaries :
 If cell is on edge of grid, velocity perpendicular to the edge = zero
To make a jet :
 If cell is in this location, velocity = some velocity
...But if you want to make a good jet, make the edges smooth, and say:
 velocity = (0.1 or some small number) * (desired_velocity - actual_velocity)

Another important aspect of making a good fluid is having colors or particles move around.
To make color moving around, make a separate buffer that has access to the velocity field and use the sampling procedure in the "Where is this information" section to update the color at each cell. So that would look like this:
vec2 sample_velocity (vec2 U) {return the velocity at "U";}
vec4 sample_color (vec2 U) {
 prev_U = U - sample_velocity(U);
 return texture(color, prev_U/resolution.xy);
}
To make particles move around, you can use voronoi particle tracking.
This going to look kind of like color tracking, but instead you use nearest texture filtering and you store the location of the nearest particle in the state vector of each cell. Backwards advect to get a guess of where the particle just was, and then test all of the adjacent cells for the closest particle to the current cell. Once you found the best particle, add the velocity at its stored location to its location so that you will know where it is on the next iteration. Then you can go back and draw the particles by drawing a circle using the distance from the current cell to the location of the particle stored in it. The cells don't necessarily have to contain the particle they are tracking, it just has to be the closest particle to that cell.
Here is an example of voronoi tracking.

Please email me with any comments or questions! wyatthf@gmail.com