I want to simulate fluids. There are two
main approaches: using particles which is called Lagrangian, or using fixed cells which is
called Eulerian in this video I want to focus on the fixed cell Eulerian approach. The particle
approach is fun though and in Unity is easy to do. I just create a tiny sphere with bounciness
zero and friction zero, and spawn a few thousand of them. I might experiment with particles in
a future journey. It's also possible to have a hybrid, where particles morph into cells, that
morph into particles, but I will leave that for future experimentation. There are so many places
I want to go that I can't. Like inside the earth, through the sun, see black holes colliding. I
can't do these things in reality. But maybe we can simulate them? In the Eulerian approach,
the fluid is represented by stationary cells. Each cell has a single vector which represents
the velocity of the fluid in that cell. Let's add some colored ink to the fluid. The ink
flows in the direction of the velocity. This is called ‘advection’. I'll talk about how we
can Implement advection a bit later. To simulate places I cannot go I feel a useful first step is:
can I simulate a simple fluid? No viscosity. No air. And I'm going to keep to two dimensions. So
that's my goal for this video. And I'm going to use GPU shaders. Because shaders are fun. Let's
make the color sources persistent and then I'll add a drop of each color. Let's increase the
resolution to 30p and turn off the arrows. I started this journey by watching some
videos by Matthius Müller who is physics research lead at Nvidia. I read through Jos
Stam's realtime fluid dynamics for games. I feel that Brian Trevethan’s thesis “Physically
Based Simulation for Computer Graphics” gives an excellent survey of fluid simulation techniques.
I couldn't understand his maths at first, so I spent several weeks watching videos on maths
by Commutant, Elucyder, and Dr Peyam. “In French ‘la place’ means ‘the place’. So you know you're
in the right place!” It feels a little static, so let's add some noise to the color. This gives
a little sensation of movement. We can actually increase the resolution more, to 480p. This is
pretty fun, but if I go back down to resolution 8p and turn the arrows back on we can see that
whilst the fluid is technically moving because we have velocity arrows, but the arrows themselves
are not moving. The fluid movement is not moving itself and look I drag these arrows so they're
all pointing to the same cell more or less. So that means we’ve got fluid flowing into the
cell and no fluid flowing out. So the fluid is either being destroyed or it's being infinitely
compressed, like a black hole. But generally it's not behaving like a typical fluid. There are two
things that we need to make the simplest fluid simulation that I'm aiming for in this video.
One of these things is conservation of mass, incompressibility. When fluid enters a
cell faster than it leaves this is called ‘negative divergence’ and when fluid leaves
a cell faster than it enters this is called ‘positive divergence’. So here we have white
fountains that represent matter being created, and black whirlpools represent matter being destroyed.
Now to make our simulation we need to clear this divergence. Jos Stam used some fancy maths for
this using Helmholtz Decomposition. But the maths was hard for me to understand and calculating
the curl-free field needs an iterative approach. Matthius Müller proposes an approach that does
still need iteration but I find much easier to understand. We do need to change representation
slightly though. Instead of storing the velocity at the center of each cell we're going to store
the velocity at the edge of each cell. It's a little less intuitive to look at but calculating
divergence is now really easy. All we have to do is add up the velocities flowing out of the cell
and subtract the velocities flowing in and that's our divergence. In my shader code I take the
velocity flowing out of the right hand side, add the velocity flowing out the top, and subtract
the velocity flowing into the left hand side, and into the bottom, and that's the divergence.
To clear this divergence what Matthius proposes is we take the divergence, divide it by four
- like the number of edges of the cell - and subtract that from each of the edge velocities.
This will immediately set the cell's divergence to zero. So I can click on these cells to clear
the divergence. When I clear the divergence for a cell I increase the divergence of its neighbors.
But that's okay: I can just keep clicking. It's a little like trying to put a screen protector on
your phone and there's a bubble of air trapped and you're pushing it around. But in our case
though the bubbles will go away eventually, if we click enough. In my shader code each
thread takes the divergence of the cell, divides it by the number of fluid neighbors to
form deltaV then subtracts deltaV from each of the edge velocities. The divergence of the cell
immediately becomes zero. Clearing divergence is also called ‘projection’ which is why this
function name starts with the word ‘project’. Let's go back to showing the center velocities
because I think they're easier to understand. Rather than clicking by hand I can get the
computer to click for me. And I can speed this up a bit. If I clear this, and draw a single
outwards velocity, the projection will gradually produce arrows flowing backwards, so that the
flow loops around. Of course I didn't spend 4 weeks writing shader code just to run projection
cell by cell. But there's a problem. Threads on a GPU are a little like rowers in a rowing team.
They move together in lock step. They all run the same instruction at the same time. This means that
if multiple threads try to update the same value, only one will succeed. On the left we have
some threads symbolized by spirals. Next to each thread we have an input color. I can
connect an input color to an output color using an arrow. When I click ‘step’ each thread
adds its input color to each connected output color. What happens if I connect two inputs
to the same output? We might expect the inputs to mix to form a blend. But actually only one
thread succeeds. The other thread’s update is ignored. Here is the shader code for this
demo. Each thread runs this function once, each with its own thread ID. Each thread loops
over the outputs. If there is a connection from the thread's input to an output then the
thread adds its own input to that output. When using a shader each output value should
be updated by only a single shader thread. Going back to clearing divergence, when I clear
the divergence for two adjacent cells each of them modifies the velocity of their shared edge.
So to clear the divergence of two neighbor cells we need to run the cell updates in sequence. On
a shader we want to run many threads in parallel, so this sounds impossible. There's a
solution though which Matthius shows which is called red-black Gauss Seidel. if
you look at all the cells I'm selecting, none of these cells have shared edges, so I can
update these cells together, in parallel, without write conflicts. And then we can select another
set of non-adjacent cells and clear those. There are four sets of such non-adjacent cells in total.
So, in four passes, we can update all cells. It's also possible to use only two passes but I get
slightly better frame rate using four passes. Let's run projection continuously. All right you
can't see it but it's running projection. So now, when I drag velocity around, the other flows
adjust immediately to keep divergence zero everywhere. What this means is that the
velocity loops around back to where it came from. It doesn't just vanish into
nowhere or be created from nowhere. If if I try to change the flow to point upwards
somewhere there's going to be downwards flow. Let's add some persistent upward flows. So now we
have some very strong back flows, like rip tides in surf. Let's bump up the simulation resolution
to 80p and use color to show velocity. I'm using HSV color to depict the velocity. The hue shows
the direction. There's a color wheel here. Black is stationary and then the brighter the color the
faster. These spots are yellow meaning they are upwards velocity. This blue color shows downwards
motion. This red color shows motion to the right, and the cyan shows motion to the left. Because
projection is running we get these gorgeous color gradients as the flows loop around. But we don't
really see movement so let's try switching to colored ink mode and add some colored ink and we
can drag the velocity around which is pretty fun. But we still don't really see movement
because the ink reaches a steady state so let's add some color noise. The noise shows
movement but it adds a bunch of extraneous color around the edges. So let's add some
decay to the ink. All right so we got rid of the extraneous color and we have some kind
of sensation of movement. But now our actual color kind of faded a lot. So what we can can
do is: we can increase the advection speed. So let's bump up the advection speed to about
140. And then I can drag the velocity around. All right, so it's not actually clear what is
different compared to not having projection on. let's reset, turn off projection, and draw a
single line of upwards persistent velocity. If I switch to the HSV view we can see that there is no
velocity except for that line. Now let's turn on projection. So now the whole simulation fills with
a beautiful color gradient effect. And this means that everywhere has some velocity. Now let's go
back to Ink color View and draw some ink sources. All right so we've definitely got movement now
and it seems plausible that we're enforcing incompressibility. But does this look right?
Does fluid normally spin around like this? Like maybe in a food blender. But what we have
here is a single source of directional velocity, like a hair dryer, or a fan. Wouldn't
we expect that fluid to have momentum and carry on in a straight line,
like water out of a tap? We would. To add conservation of momentum what we need is
to add advection of velocity. We saw advection before for the colored ink. ‘Advection’ means
moving things in the direction of the velocity. ‘Advection of velocity’ means: we're going to move
the velocity itself in the direction of itself. This seemed a bit mysterious to me. But if if we
switch to using particle representation it becomes a bit more obvious. In the case of particles
each particle has a velocity, and the velocity of each particle moves with that particle.
If I add Eulerian cell colors the cell colors also move with the particles, in the direction of
the particle velocity. This is especially obvious for the particles on the edges. So we need to add
this velocity advection to our Eulerian cell based simulation. Let's change back to low resolution
and turn the arrows back on. So if I move the arrows so the tail is in the cell centers then the
arrows point to where the velocity needs to move to next. So what we could do is take the velocity
from a cell and move that along the arrow to the head of the arrow, and somehow spread the velocity
over the nearby cells. But there's a problem. If you look at this cell here there's multiple
arrows pointing into it. In my shader I want to use one thread for each cell: one thread per
arrow. So two arrows pointing to the same cell means that we have two shader threads updating the
same value. From the demo earlier we saw that this is impossible: one of those updates will fail. But
what we can do is change from a push to a pull. Let's go back to the thread update demo. We
can assign each thread to an output instead of to an input. Each thread will pull colors
from the inputs. This way each output is only updated by a single thread. If we have two
inputs pointing to the same output that works just fine. switching from a push to a pull
can also help us run advection on a shader. Here is the shader code for this demo. Each
thread loops over the inputs. If there is a connection from the input to that thread's
output, then the thread adds that input value to its output variable. Finally the
thread writes its output variable to its output. And this works okay. So let's
use this pull approach for advection. I can slide the arrows backwards so that they
point into each cell rather than out of each cell. Then each thread can pull the velocity from
the tail of its arrow into its own cell, and this will run just fine on a shader. Unfortunately the
tail of each arrow won't usually lie exactly over a single cell, so we will need to interpolate
over the four nearest cells. Here is my shader code to interpolate between four cells given
the index of the bottom left cell and a lambda interpolation vector. The actual code has a bit
more bounds checking but it's otherwise similar. The advection is a bit fiddly because
we not actually storing the velocities at the center of each cell, but rather the
staggered velocities at the edges of each cell. Debugging this took me quite
a while but I got there in the end. Without projection, nothing happens when I
turn on advection and draw some velocity. This is because we are doing a pull. The
cells doing the pull have zero velocity, so the length of the arrow to the source cell
is also zero. Each cell pulls from itself. Since each cell has zero velocity itself,
then it pulls zero velocity from itself, and nothing changes. If we turn on projection
this problem goes away. The projection smears non-zero velocity everywhere, with this
beautiful gradient color effect. Even if the velocity is near zero it's still enough to
kick-start the advection process. It's not enough to only apply a projection though. With only
projection the fluid does not stream forward in a straight line. We just get small localized
loops of fluid flow. So there's no divergence but it's basically static. With both projection
and advection I feel it looks quite nice. Having to turn on projection in order to test
advection made debugging really hard. It was hard to know if bugs were in my projection code
or in my advection code. And let's be honest bugs were everywhere. Really I want to be able to
test projection and infection separately. This will make debugging easier. So let's think
about about trying to run advection without projection. The problem is without projection on
the advection doesn't pull from zero velocity, and nothing happens. So how to test? So one option
is to click somewhat and drag backwards in the opposite direction of the velocity. Since my mouse
position is supplying velocity this velocity flows forward from my mouse, and so we get a stable line
of velocity. When I release the mouse the velocity then advects away from where my mouse was. I
thought that adding noise to the velocity might help since then the velocity would be non-zero
everywhere. And it kind of does. But it turns out the direction of the velocity matters: and
it's wrong half the time. So let's turn the noise off. What does work is to pre-warm an area with
some velocity that's going in the same direction I want to run advection. And then I can advert
velocity over that. Still it's hard to be sure if it's really correct. In the end I wrote a bunch
of shader kernels just for unit testing. This was quite a lot of pain but at least I could be
somewhat confident that my code was working-ish. Let's turn projection back on, and then bump up
the resolution to 480p. It's a little slow so let's bump up the advection speed to 140. That
seems pretty fluid. I'm going to draw a line of persistent upwards velocity. All right so that
flows upwards and it hits the ceiling. But what's wrong with this picture? So we've got this huge
block of fluid moving upwards, and it's hitting the ceiling. The ceiling is solid. The fluid
cannot go through it. So the fluid's vanishing. Like, there's no fluid coming out the sides. So
this is violating conservation of mass. But we already have projection. Shouldn't projection fix
that? Well when we added projection we saw that we need to iterate the projection step and I'm
not. I'm just running the projection step once, in between each advection step. It turns out that
the number of projection iterations is important, and we need quite a few more than one,
in between each advection step. So let's try 3. and 5. 7. 10. So increasing the number of
projection iterations is definitely helping. What we need is for the information from the center
of the upflow to get to the edge of the upflow before we take the next advection step. Each
red black Gauss Seidell iteration communicates information about two pixels each time. Like no
more than two pixels. And I guess we've got like 50 to 100 pixels here. So we're going to need
at least 50 iterations on or so. So let's try 50. And let's try 100. All right, that's looking
a lot more convincing. And something I like is that now we see some nice turbulence effects. I
think turbulence is really fun and interesting. I think it looks much more convincingly like
a fluid with turbulence in. And I didn't have to do anything fancy to get the turbulence. I
just applied a few fundamental laws of physics. Conservation of mass. Conservation of momentum.
Incompressibility. I think that's pretty fun. Let's load in some velocities and colors I
drew earlier. And let's dab on a few spots of green ink. This is running at 50 frames
per second on my MacBook Air. I think this looks quite convincing as a fluid, and
I'm pretty happy with how this is looking. What I want to look into next is simulating
the inside of the Earth, the Sun, and black holes. If you're interested too please consider
subscribing. If you liked this video and you want me to make similar videos please consider
encouraging me by clicking like. Cool. See you!
