When I first stumbled upon fluid simulations in game dev I was amazed on how good the effect could be. I really wanted to learn how this works, but learning materials on this topic are suprisingly sparse - and those which I found were pretty difficult to understand. Nonetheless, I decided to give it a try; and - while I’m at it - why not create a blog post out of it to hopefully make it easier for the next guy?
Before we begin, I want to stress a few points:
- I’m not a mathematician. If you find errors in my explanations, please DM me on Bluesky or send an email and I’ll correct it
- This implementation is for learning purposes only. For that reason it is implemented in a way that is not the best performance-wise. Calculations are made solely on the CPU. We introduce way too many variables. All of that to make it easier to read and learn, not necessarily to squeeze the most FPS.
Learning materials I used:
- “Real-Time Fluid Dynamic for Games” by Jos Stam (PDF)
- “Fluid Simulation for Dummies” by Mike Ash
You can find all of the code in this repository: github.com/rskupnik/godot-fluid-simulation-demo
I used git commits to mark code checkpoints matching the chapters of this blog post, so if you don’t want to write the code alongside reading, you can make use of the commit view to follow along. For your convenience, I include a “project snapshot” and a “diff” link at each chapter, which lead to, respectively: the codebase at the discussed point and the commit diff view
AI disclosure: Every word of this blog post and every line of this codebase is written by me. All the diagrams and videos were created by me. AI was only used for research.
Finally, if you appreciate such work then consider buying me a coffee :)
Foundations
The algorithms we will use are based on physical equations for fluid flow - the Navier Stokes equations. Our use case is a game dev one - so we want to sacrifice precision of these calculations in favour of speed. The effect needs to be good enough but not overly expensive. We achieve this in three ways. First of all, we use a relatively small grid with large cells. Second - we advance the simulation in arbitrary time steps. Finally, we use approximation equations (such as Gauss-Seidel relaxation) to arrive at good-enough solutions to some equations.
Let me start with the mathematical description of what we will do in this blog post. This description might sound daunting, but don’t worry - we’ll explain everything as we go. Here goes: we will simulate fluid flow by moving a scalar density field through a vector velocity field. We’ll simulate velocity diffusion and advection as well as density diffusion and advection. Then we will add velocity projection with the goal of making the fluid obey the law of mass conservation - which will happen by balancing divergence with a pressure field. We will use bilinear interpolation and Gauss-Seidel relaxation for approximating values where needed.
Alright, with all of that out of the way, let’s begin!
The journey begins with a grid
Create a new Godot project and add a Node2D. I called mine “FluidGrid”. Attach a script to it. All the code goes into that script.
First, let’s define the grid:
@export var N := 16
@export var cell_size := 32
var size := 0
Here, N is the amount of cell in a row and column - basically the size of the grid - and cell_size is the size of a single cell in pixels. We also define size, which we will soon initialize. We will set it to N+2, because we want borders as well.
Next, let’s add arrays for storing the actual data.
# Density means "how much material does this cell contain"
var density: PackedFloat32Array
var density_prev: PackedFloat32Array
# "u" stores the horizontal velocity (x direction)
var u: PackedFloat32Array
var u_prev: PackedFloat32Array
# "v" stores the vertical velocity (y direction)
var v: PackedFloat32Array
var v_prev: PackedFloat32Array
For now, we need three arrays - density will store the density, u will store horizontal velocity and v - vertical velocity.
Density tells us how much material does a given cell contain - it will range between 0.0 and 1.0, where 0 is empty and 1 is full. Technically, it can go above 1.0 but we will only display up to 1.0. The way we display density is simple - with color. A cell full of density will be fully white and without density - fully transparent.
The velocity arrays also store floats and they describe velocity at a given cell. A velocity of 0 means no movement, then it can go positive or negative which corressponds to going right or left (for horizontal) and down or up (for vertical). Combined, they tell us the velocity at a given cell. We could just store the velocities as a single array of Vector2f, but it will be much easier for calculations if we separate them as two float arrays. When it comes to displaying this information - we will draw small blue arrows for each cell.
You might also wonder why each array has a _prev equivalent - that’s because for some of the calculations we will iterate over the real arrays and modify the data live - and in those cases we need to “snapshot” the data before we start iterating, so we can see what the values were before we started modifying them. That will be mostly used in the approximating equations. I follow the naming convention of Stam’s paper with this _prev name, although I believe _snapshot would be a more descriptive name.
Alright, let’s move on. Time to initialize all of these!
func _ready():
# Resize all the arrays properly
# We use a single-dimensional array to store the grid, which is why we need to multiply N
# The "+2" is added for borders, because there are two for each dimension (x, and y)
# For x dimension, there's a single cell border on the left and on the right, hence "+2". Same for the y direction
size = (N + 2) * (N + 2)
density.resize(size)
density_prev.resize(size)
u.resize(size)
v.resize(size)
u_prev.resize(size)
v_prev.resize(size)
queue_redraw()
We could use a two-dimensional array, but it will be simpler to work with a single-dimensional one. We just need one helper function to make it easier to index this array.
# This is a helper function that makes it easier to work with a grid when it is
# packed into a single-dimension array
# We can call it with the cell index (i and j) and it will translate it into
# an index in the single-dimension array
func IX(i: int, j: int) -> int:
return i + (N + 2) * j
With all of that, we can now implement the _draw() function to display the grid.
# This is the standard Godot function used for drawing
# We want to draw a simple grid of (N+2)*(N+2) rectangles of size cell_size
func _draw():
for j in range(0, N + 2):
for i in range(0, N + 2):
var x := i * cell_size # this translates the index into pixel position on screen
var y := j * cell_size
var rect := Rect2(x, y, cell_size, cell_size)
var is_boundary := i == 0 or j == 0 or i == N + 1 or j == N + 1
var fill := Color(0.16, 0.08, 0.08) if is_boundary else Color(0.08, 0.08, 0.08)
draw_rect(rect, fill, true)
draw_rect(rect, Color(0.35, 0.35, 0.35), false)
It’s pretty self-explanatory. We iterate over the grid and draw each cell as a simple Rect2, changing the color slightly for the border cells.
You can now run the project and you should see this:
Putting “fluid” in “fluid simulation”
Now that we have a grid, let’s add some fluid to it.
We will start very simple - we’ll make it possible to add density to a cell by clicking it with a mouse. Then we will make the density slowly fade away - this will be useful later so we can experiment without making the grid fill with fluid and requiring a restart.
# This helper function translates the position we clicked on with the mouse
# into the cell coordinates
# So if we click somewhere in the grid, it will return a Vector2i, where the
# first element is the index of the cell in that grid in x dimension
# and the other element is the index of the cell in the y dimension
func cell_from_mouse(pos: Vector2) -> Vector2i:
return Vector2i(floor(pos.x / cell_size), floor(pos.y / cell_size))
# This is the standard Godot function for processing input
# We want to detect a mouse click and inject density into the clicked cell
# Density is represented as a float number and is stored in the "density" array
func _input(event):
if event is InputEventMouseButton and event.pressed:
# figure out the cell that was clicked
var cell := cell_from_mouse(to_local(event.position))
var i := cell.x
var j := cell.y
if i >= 1 and i <= N and j >= 1 and j <= N:
density[IX(i, j)] += 1.0 # inject density into the cell
queue_redraw() # tell Godot to redraw the grid
And now we need to draw what’s inside the density array. Go to the _draw function and change the var fill := Color... line to this (see the diff linked above if confused).
var fill := Color(0.08, 0.08, 0.08)
if is_boundary:
fill = Color(0.16, 0.08, 0.08)
else:
# Even though density can go above 1.0, we need to clamp it to values between 0.0 and 1.0 for drawing
var d : float = clamp(density[IX(i, j)], 0.0, 1.0)
fill = Color(d, d, d)
We make use of the IX() function to index the density array and use that to determine the “intensity” of the colour in a cell. For now, each click of a mouse injects 1.0 density into a cell, so that cell should turn white.
This is the effect
Fade away, fade away, fade away
As mentioned, we want to add a simple fading effect so the density slowly disappears - to avoid clogging our grid later on.
Let’s begin with a simple variable controlling the intensity of this effect:
@export var density_fade_rate := 0.1
Now, time for the fade_density() function
# Fade density as the time passes
func fade_density(delta: float) -> void:
for j in range(1, N + 1):
for i in range(1, N + 1):
var idx := IX(i, j)
# We need to multiply the rate of density fade through delta
# to make it the same despite the framerate
density[idx] = max(0.0, density[idx] - density_fade_rate * delta)
In this function, we go over each cell and decrease the amount of density in that cell by the amount denoted by the rate of density fading multiplied with the time delta. If you come from the gamedev world, you are probably familiar with delta, but just for the sake of completion - this variable is provided by the Godot engine itself and it contains the amount of time that has passed since the last frame was drawn. It is meant to be used in various equations to either simulate the passage of time or bind some effect to the user’s framerate.
Finally, we need to call this fade_density() function every frame, which we will do using Godot’s standard _process() function
# This is the standard Godot function called every frame
# It's the heart of our simulation
# The "delta" variable holds the amount of time that passed since the last frame
# For now we use it to slowly fade the density
func _process(delta: float) -> void:
fade_density(delta)
queue_redraw()
Since we now modify the density array every frame (by fading it away), we also need to redraw the grid every frame, which is what queue_redraw() is for. Again - it’s a Godot built-in function.
At this point, you should be able to click the cells to add density and see how it slowly fades away
Time for some arrows
Alright, now that we can see the density, it’s time to also visualize velocity.
Start with a variable for controlling the scale of the arrows
@export var velocity_draw_scale := 20.0
This can be used to make the arrows prettier. It doesn’t affect the velocity itself, it only affects the scaling factor of the arrows we draw. Want bigger arrows? Increase the param. Arrows are way too big? Scale them down! Again, it’s important to note that this param does not affect the simulation in any way - it’s purely visual.
Our velocity field is currently initialized to 0.0 in every cell - so let’s hardcode some temporary velocities to make sure the drawing works properly. You should remove those at the end of this step.
Add those in _ready()
# TEMPORARY
# In this cell, we set the horizontal velocity to 1.0 and vertical to 0.0
# Result: horizontal arrow pointing right
u[IX(8, 8)] = 1.0
v[IX(8, 8)] = 0.0
# In this cell, we set the horizontal velocity to 0.0 and vertical to -1.0
# Result: vertical arrow pointing up
# Remember that in Godot, the y axis goes from top to bottom, hence why -1.0 points up
u[IX(9, 8)] = 0.0
v[IX(9, 8)] = -1.0
# In this cell, we set the horizontal velocity to 1.0 and vertical to -1.0
# Result: vertical arrow pointing up and right (diagonal)
# Remember that in Godot, the y axis goes from top to bottom, hence why -1.0 points up
u[IX(10, 8)] = 1.0
v[IX(10, 8)] = -1.0
Ok, now we need a function to draw the arrows. Before we do that - rename the current _draw() function to _draw_grid(). Then add the _draw_velocity_arrows() below.
# Draw velocity arrows in inner cells with tiny circles as arrow tips
func _draw_velocity_arrows():
for j in range(0, N + 2):
for i in range(0, N + 2):
var is_boundary := i == 0 or j == 0 or i == N + 1 or j == N + 1
if not is_boundary:
var idx := IX(i, j)
var center := Vector2(
i * cell_size + cell_size * 0.5,
j * cell_size + cell_size * 0.5
)
var velocity := Vector2(u[idx], v[idx])
var end := center + velocity * velocity_draw_scale
draw_line(center, end, Color(0.2, 0.8, 1.0), 2.0)
draw_circle(end, 2.5, Color(0.2, 0.8, 1.0))
In this function, we define a Vector2 which consists of the horizontal velocity in this cell (the u array) and the vertical velocity in this cell (the v array). We then define a center variable, which denotes the center of the cell and also the starting point of our velocity arrow; and an end variable which utilizes our velocity_draw_scale param. Then we pass both center and end to draw_line (which is a Godot function), which tells Godot to draw a line starting at center and ending at end. Since we scale the end with velocity_draw_scale, the arrow gets proportionally longer or shorter depending on how big we set velocity_draw_scale. Finally, we draw a tiny circle at the end to act as a sort of “arrow point”.
Alright, but we still need to call this function. Let’s restore _draw(), defined like this
func _draw():
_draw_grid()
_draw_velocity_arrows()
So now Godot’s standard _draw() function will draw both the grid itself (with density) and velocity arrows
Since we hardcoded a few velocities, you should now see this:
Let’s get this moving!
Alright, before moving further let’s just quickly revise our _input() so that we can inject both density and velocity by dragging the mouse.
First, some variables
@export var velocity_add_scale := 0.06
var is_dragging := false
var last_mouse_cell := Vector2i.ZERO
The velocity_add_scale will control how much velocity we add with a single mouse stroke. The rest are operational variables for tracking the mouse.
Now replace the entire _input() function with this
# This is the standard Godot function for processing input
# We want to detect when a mouse is dragged while clicked and the inject
# both density and velocity at the relevant cells
func _input(event):
if event is InputEventMouseButton:
is_dragging = event.pressed
last_mouse_cell = cell_from_mouse(to_local(event.position))
if event is InputEventMouseMotion and is_dragging:
var local_pos := to_local(event.position)
var cell := cell_from_mouse(local_pos)
var i := cell.x
var j := cell.y
if i >= 1 and i <= N and j >= 1 and j <= N:
# event.relative stores the relative difference between last time
# this function was called and this time
# in this case, it tells us how far the mouse has travelled
# we use that to decide how much velocity to add
var delta_velocity : Vector2 = event.relative * velocity_add_scale
var idx := IX(i, j)
# Inject the velocity into the cell
# u stores horizontal velocity, v stores vertical velocity
u[idx] += delta_velocity.x
v[idx] += delta_velocity.y
# Inject density
density[idx] += 1.0
queue_redraw()
If you run the project now, you can drag the mouse while holding left click to add both velocity and density. If you want to add density in a single place just keep the left mouse button clicked and do small drag movements in that cell.
Hold your horses
As a final step before me move to the Serious Stuff, let’s make velocity also slowly fade away. We will remove this effect at the end of this project, but for now it will come in handy.
First, as usual - a variable controlling the strength of this effect
@export var velocity_fade_rate := 0.2
Now a function that handles velocity fading
# Fade velocity as the time passes
func fade_velocity(delta: float) -> void:
for j in range(1, N + 1):
for i in range(1, N + 1):
var idx := IX(i, j)
# We need to multiply the rate of velocity fade through delta
# to make it the same despite the framerate
u[idx] = move_toward(u[idx], 0.0, velocity_fade_rate * delta)
v[idx] = move_toward(v[idx], 0.0, velocity_fade_rate * delta)
This is almost identical to how we fade density. The only real difference is we are using Godot’s move_toward function which smoothly interpolates the movement - just so it looks nice and smooth.
Finally, let’s make sure this function gets called from within _process()
func _process(delta: float) -> void:
fade_density(delta)
fade_velocity(delta)
queue_redraw()
With this, the velocity arrows should now slowly fade back to nothing
Let’s get the flow going
Alright, all the things we did so far were pretty basic. This chapter is where we start working on the mechanics that make the magic happen.
We start with density advection. You already know that density is basically our fluid - it’s a grid of floating point numbers that tells us how much material is currently present in a given cell. Advection is all about making that density move through our velocity field.
We have an array of density values and want to move them along the velocity vectors. How would you implement it? If you’re like me, your first instinct might be to do some variation of “iterate over the grid, grab the density at a given cell, then grab the velocity at that cell and make the density move to wherever the velocity arrow points to”. Apparently - this would be the hard way to approach this. I’ll be honest with you - I don’t fully understand why this would be difficult to calculate, but if people smarter than me claim so - I believe them.
According to Stam’s paper, a better approach is to track density backwards in time. It works like this - first, we snapshot the density array. Then, for each cell, we ask the question - where did the density value currently in this cell come from? To answer that question, we could look at the velocity arrow in that cell. Ok cool, but that only tells us where the density should go, not where it came from. To figure that out, we can simply reverse that velocity arrow and see where we land. Then we can interpolate the values between the four surrounding neighbours of that point to figure out how much density should we take from there.
Now, there might be a red lamp blinking in your head right now. Just because the velocity at the currently inspected cell points, for example, to the right - doesn’t mean that the density that came to this cell came from the left. It might have come from the top neighbour. Or the bottom one. I also asked myself the same question when first studying this algorithm - and apparently the answer is twofold. First of all - we are dealing with a velocity field. In such a field, it is very unlikely that the velocities would make a sudden turn like that. A situation where you have a velocity pointing straight down and then suddenly straight to the right is unlikely. This is a velocity field - they change gradually. Second of all - this is one of the places where we surrender perfect precision for good enough approximation.
These two arguments combined are why this solution makes much more sense - at least in my head :)
So what we are going to do is this (for each cell):
- Find the velocity, reverse it, multiply by the time delta - that will give us an approximation of where the density probably came from
- We will end up with a point that lays somewhere in the grid, not necessarily in the center of a cell. Let’s call that the “source cell”
- We will find all the neighbours of that source cell
- We will figure out how close the source point is to each fo the neighbours. Remember - we don’t necessary land on a cell center. The source point might be closer to the left and top (for example)
- With the source position, the source neighbouring tiles and the information of how close to each of the neighbours we are, we will use bilinear interpolation to figure out how much density we should take from each of those neighbours. The closer the point is to a neighbour - the more density should be taken from that neighbour (and less from the opposing neighbour)
Alright, let’s implement this. I will provide the rest of explanations as code comments to not break the continuity of the function
func advect_density(delta: float) -> void:
# Multiply delta by N to scale the time by grid size
var dt0 := delta * N
for j in range(1, N + 1):
for i in range(1, N + 1):
var idx := IX(i, j)
# Notice this tracks backwards through time
# First we subtract the velocity in this cell (multiplied by delta)
# from the cell index
# That tells us where the density in this cell came from (probably)
# Then we clamp to make sure we don't go outside the grid
# The result will be a position of where the density probably came from, lying somewhere
# between other grid cells (not necessarily in a cell center)
# For the sake of this example, let's assume we land at a point (7.3, 4.8)
# Which means it's the lower-left corner of the cell (7,4)
# (The center would be at (7.5, 4.5) )
# Remember - in Godot, the y axis points down
var x := i - dt0 * u[idx]
var y := j - dt0 * v[idx]
x = clamp(x, 0.5, N + 0.5)
y = clamp(y, 0.5, N + 0.5)
# This grabs the four cells surrounding the point we just calculated above
# Assuming we land at (7.3, 4.8)
# The four cells around this point (including the one it is part of)
# will be: (7,4) (8,4) (7,5) (8,5), because
# i0 = 7
# i1 = 8
# j0 = 4
# j1 = 5
# (it's a combination of all "i" with all "j")
var i0 := int(floor(x)) # floor will cut off the decimal part, rounding down, so 7.3 -> 7
var i1 := i0 + 1
var j0 := int(floor(y)) # 4.8 -> 4
var j1 := j0 + 1
# These are fractional weights
# They tell us how close the point is to each side of the cell
# Consider example where our point is at (7.3, 4.8)
# That gives us: s1 = 0.3, s0 = 0.7, t1 = 0.8, t0 = 0.2
# Basically we extract the decimal part into s1 and t1
# And then we put "the rest that is missing to 1.0" into s0 and t0
# ---
# What this tells us, is basically that:
# the point is at 30% distance from left side and 70% distance from right side
# which means it should be influenced stronger by the left side and less by the right side
# Same for y dimension, the point is at 80% distance towards bottom and 20% towards top
var s1 := x - i0
var s0 := 1.0 - s1
var t1 := y - j0
var t0 := 1.0 - t1
# With neighbouring tiles and the fractional weights calculated,
# we now need to decide how much density we should grab from each neighbouring tile
# based on how close we are to it (fractional weights tells us that)
# This function does exactly that. It's a standard mathematical operation called
# "bilinear interpolation". You might be familiar with "linear interpolation", often
# called "lerp" in the gamedev world, which interpolates between values in a single
# dimension (line). Bilinear interpolation does the same, except in two dimensions (lines)
# ---
# Notice also that we make use of the density_prev array here. That's because we
# keep updating the density array as we iterate, so we need to "snapshot" the densities
# before we start iterating - otherwise we would mess up results as they would keep changing
# So density_prev is pretty much a "snapshot of density before we start modifying it"
density[idx] = (
s0 * (t0 * density_prev[IX(i0, j0)] + t1 * density_prev[IX(i0, j1)]) +
s1 * (t0 * density_prev[IX(i1, j0)] + t1 * density_prev[IX(i1, j1)])
)
Now, let’s add a helper function to take the snapshot of density
func copy_density_to_prev() -> void:
for idx in range(size):
density_prev[idx] = density[idx]
And make sure our new code gets called. Replace _process()
func _process(delta: float) -> void:
copy_density_to_prev()
advect_density(delta)
fade_density(delta)
fade_velocity(delta)
queue_redraw()
The effect should look something like this
A drop in the ocean
You know how when you drop a bit of paint into water you can see it spread around nicely? Apparently that’s what fluids do - so it would be nice if our simulated fluid did the same! It’s called density diffusion
To solve this, we will use what is called the Gauss-Seidel relaxation. If you want a deep mathematical defintion - look around the internet, there’s plenty of stuff about it. I’m just gonna explain it with how I understand it.
What we want to do for diffusion is go for each cell, grab the average of the neighbours and add it to the density at this cell. Theoretically we could also decrease the neighbours density, but we can skip it and just let our fading effect take care of that - it will be believable enough.
Alright, but if we just do one pass - go through each cell, grab an average of the neighbours and add it to the cell’s density - it will look very choppy. What we want is to repeat this process several times to smoothen this out. This is basically what Gauss-Seidel relaxation is in this context - we first guess the proper value (in the first iteration) and then continue guessing (in the subsequent iterations) until our guess approaches a believable value. It’s an approximation. This works because we are working on the live density array - it is updated in each iteration, so the next iteration works on updated data, thus getting closer and closer to what the answer should be
Start with variables
@export var density_diffuse_rate := 0.006 # Strength of density diffusion effect
@export var density_diffuse_iterations := 20 # How many iterations when diffusing density
And now the implementation
# Diffusion simply means spreading the density to the neighbouring cells
# Think of it like putting a drop of paint in water - it will spread
# This uses what is called "Gauss-Seidel relaxation", which basically means
# we keep "guessing" (or rather approximating) the value density_diffuse_rate times
# With enough iterations, this is close enough for the simulation to be believable
func diffuse_density(delta: float) -> void:
# This is the strength of diffusion effect for this frame
# Delta is how much time passed
# Then we have the density diffusion strength parameter
# Finally we multiply by N squared so it square with cell size
var a := delta * density_diffuse_rate * N * N
for k in range(density_diffuse_iterations): # We iterate 20 times
for j in range(1, N + 1):
for i in range(1, N + 1):
var idx := IX(i, j)
density[idx] = (
density_prev[idx] + # This is the "anchor" - the original value at this cell, which stays constant throughout iterations (hence called "anchor")
a * ( # Here we multiply the strength of the effect with the sum of densities of all neighbours - this is the heart of "Gauss-Seidel" relaxation
density[IX(i - 1, j)] + # These are
density[IX(i + 1, j)] + # the four
density[IX(i, j - 1)] + # neighbours
density[IX(i, j + 1)] # :)
)
) / (1.0 + 4.0 * a) # This balances the math, since we add densities from 5 cells (this one + 4 neighbours). Without this, the density would grow too much (try it!)
Don’t forget to update _process()
func _process(delta: float) -> void:
copy_density_to_prev()
diffuse_density(delta)
copy_density_to_prev()
advect_density(delta)
fade_density(delta)
fade_velocity(delta)
queue_redraw()
And now we can see this nice diffusion effect. Remember you can hold left mouse button and move around in small circles to keep adding density in a single cell
Set some boundaries
It’s time to set some boundaries - and by that I mean making the walls of the simulation behave like walls. The effect won’t be much right now but it is necessary to implement.
The approach is rather simple. For each boundary cell, we want to copy the density value of the neighbour. That will make the regular cells around it behave properly in our advection and diffusion calculations. We also want to make those boundary cells reverse the velocity values of the nearest real cell - so a left boundary cell should look at the real cell to the right (which is the leftmost real cell) and reverse the horizontal velocity. That will make that boundary cell act as a wall which “reflects” the fluid. Same for top and bottom boundary cells and vertical velocities. For corner boundary cells we’ll simply average between their boundary neighbours.
Now, this equation described in Stam’s paper is condensed, but also difficult to understand. I have decided to include the original function to show how it could look like - but then we will implement it on our own in a way that is less condensed and probably less performant, but should be easier to understand what is going on.
So first - here is the original function as described by Stam. I only show it for reference - it is not called by anything in the codebase.
# This is the "original" function as defined in the Stam paper
# It is NOT USED by this program. I left it only for reference
# It is much more performant because it doesn't declare needless variables
# But that also make is harder to understand what is going on here
func set_bnd_original(b: int, grid: PackedFloat32Array) -> void:
for i in range(1, N + 1):
grid[IX(0, i)] = -grid[IX(1, i)] if b == 1 else grid[IX(1, i)]
grid[IX(N + 1, i)] = -grid[IX(N, i)] if b == 1 else grid[IX(N, i)]
grid[IX(i, 0)] = -grid[IX(i, 1)] if b == 2 else grid[IX(i, 1)]
grid[IX(i, N + 1)] = -grid[IX(i, N)] if b == 2 else grid[IX(i, N)]
grid[IX(0, 0)] = 0.5 * (grid[IX(1, 0)] + grid[IX(0, 1)])
grid[IX(0, N + 1)] = 0.5 * (grid[IX(1, N + 1)] + grid[IX(0, N)])
grid[IX(N + 1, 0)] = 0.5 * (grid[IX(N, 0)] + grid[IX(N + 1, 1)])
grid[IX(N + 1, N + 1)] = 0.5 * (grid[IX(N, N + 1)] + grid[IX(N + 1, N)])
Alright, we will now implement the same concept on our own without condensing it so much, so hopefully it should be easier to follow what is going on here
enum BoundaryType {
DENSITY,
VELOCITY_HORIZONTAL,
VELOCITY_VERTICAL
}
# This function handles boundary cells, which are not part of the simulation. It effectively
# establishes a "wall". We do this by making each boundary cell copy values of their simulated neighbour,
# modifying them accordingly if needed.
# For density, we simply copy the value
# This makes both density advection (moving along velocity field) and diffusion behave properly near borders
# For velocity, we also copy the value of a simulated neighbour but we invert part of it.
# For top/bottom walls, we invert the vertical values, for left/right we invert the horizontal values
# This makes the boundaries behave as "walls" for velocity - density will bounce off them due to reflected velocity
# Note: I realize declaring so many needless variables is bad for performance. This is done on purpose
# to make it easier to understand what is going on here. Original function included, see "set_bnd_original"
func set_bnd(boundary: BoundaryType, grid: PackedFloat32Array) -> void:
# We iterate from 1 to N+1 because we want to tackle boundaries, which are a single row/column
# With a single loop going from 1 to N+1 we handle each row and column at once, because each of them
# is exactly the same size
for i in range(1, N + 1):
# These define the indexes of boundary cells (the red ones)
# We will set values for those
var left_boundary_cell_index = IX(0, i)
var right_boundary_cell_index = IX(N + 1, i)
var bottom_boundary_cell_index = IX(i, 0)
var top_boundary_cell_index = IX(i, N + 1)
# These define the nearest cell neighbour of a boundary cell that is a real cell
# We will copy values from those
var leftmost_simulated_cell_index = IX(1, i)
var rightmost_simulated_cell_index = IX(N, i)
var bottommost_simulated_cell_index = IX(i, 1)
var topmost_simulated_cell_index = IX(i, N)
# For density, we simply copy the value of the nearest simulated cell
if boundary == BoundaryType.DENSITY:
grid[left_boundary_cell_index] = grid[leftmost_simulated_cell_index]
grid[right_boundary_cell_index] = grid[rightmost_simulated_cell_index]
grid[bottom_boundary_cell_index] = grid[bottommost_simulated_cell_index]
grid[top_boundary_cell_index] = grid[topmost_simulated_cell_index]
# For horizontal velocity, we flip the values on the horizontal axis (the boundary wall
# "reflects" the velocity) and just copy the vertical values
elif boundary == BoundaryType.VELOCITY_HORIZONTAL:
grid[left_boundary_cell_index] = -grid[leftmost_simulated_cell_index]
grid[right_boundary_cell_index] = -grid[rightmost_simulated_cell_index]
grid[bottom_boundary_cell_index] = grid[bottommost_simulated_cell_index]
grid[top_boundary_cell_index] = grid[topmost_simulated_cell_index]
# For vertical velocity, we flip the values on the vertical axis (the boundary wall
# "reflects" the velocity) and just copy the horizontal values
elif boundary == BoundaryType.VELOCITY_VERTICAL:
grid[left_boundary_cell_index] = grid[leftmost_simulated_cell_index]
grid[right_boundary_cell_index] = grid[rightmost_simulated_cell_index]
grid[bottom_boundary_cell_index] = -grid[bottommost_simulated_cell_index]
grid[top_boundary_cell_index] = -grid[topmost_simulated_cell_index]
# Corners
# Boundary corner cells are a special case because they have two neighbours (both of which are boundary cells as well)
# Since we already set the values for all the other boundary cells in the loop above,
# here we just set each corner to an average of their two neighbours
# For ease of understanding I labeled the vars with SW like South-West, etc.
var corner_SW = IX(0, 0)
var corner_SW_right_neighbour = IX(1, 0)
var corner_SW_top_neighbour = IX(0, 1)
var corner_NW = IX(0, N + 1)
var corner_NW_right_neighbour = IX(1, N + 1)
var corner_NW_bottom_neighbour = IX(0, N)
var corner_SE = IX(N + 1, 0)
var corner_SE_left_neighbour = IX(N, 0)
var corner_SE_top_neighbour = IX(N + 1, 1)
var corner_NE = IX(N + 1, N + 1)
var corner_NE_left_neighbour = IX(N, N + 1)
var corner_NE_bottom_neighbour = IX(N + 1, N)
grid[corner_SW] = 0.5 * (grid[corner_SW_right_neighbour] + grid[corner_SW_top_neighbour]) # SW corner = right neighbour + top neighbour
grid[corner_NW] = 0.5 * (grid[corner_NW_right_neighbour] + grid[corner_NW_bottom_neighbour]) # NW corner = right neighbour + bottom neighbour
grid[corner_SE] = 0.5 * (grid[corner_SE_left_neighbour] + grid[corner_SE_top_neighbour]) # SE corner = left neighbour + top neighbour
grid[corner_NE] = 0.5 * (grid[corner_NE_left_neighbour] + grid[corner_NE_bottom_neighbour]) # NE corner = left neighbour + bottom neighbour
Both the original function and mine accomplish the same thing.
With this function present we now need to call it. Add this at the end of advect_density and for each iterative loop in diffuse_density (refer to the git diff if confused)
set_bnd(BoundaryType.DENSITY, density)
The effect right now is small but this is what will allow are density to properly bounce off walls later on
Spread the love
Just as we make density diffuse over time, we now should do the same for velocities. We want them to spread around gently just like density does.
The way to accomplish this is exactly the same as with density diffusion - we will use Gauss-Seidel relaxation formula to iteratively approximate the values. If you need details, I refer you back to the chapter that discussed density diffusion.
As usual, let’s start with variables
@export var velocity_diffuse_rate := 0.001 # Strength of velocity diffusion effect
@export var velocity_diffuse_iterations := 20 # How many iterations when diffusing velocity
Now a helper function to snaphsot the velocity arrays
func copy_velocity_to_prev() -> void:
for idx in range(size):
u_prev[idx] = u[idx]
v_prev[idx] = v[idx]
And now the function itself, including the call to set_bnd to update the boundary cells
# Diffusion here can be understood as "spreading among neighbour cells"
# Just as we diffuse density, we now do the same for velocity
# This is pretty much exactly the same as density diffusion, except it does it for two arrays
# (because we store horizontal and vertical velocity separately)
# It uses the same "Gauss-Seidel relaxation" formula
func diffuse_velocity(delta: float) -> void:
# This is the strength of diffusion effect for this frame
# Delta is how much time passed
# Then we have the velocity diffusion strength parameter
# Finally we multiply by N squared so it scales with cell size
var a := delta * velocity_diffuse_rate * N * N
for k in range(velocity_diffuse_iterations):
for j in range(1, N + 1):
for i in range(1, N + 1):
var idx := IX(i, j)
# Those are exactly the same formulas as for density
# See density diffusion function for explanations :)
u[idx] = (
u_prev[idx] +
a * (
u[IX(i - 1, j)] +
u[IX(i + 1, j)] +
u[IX(i, j - 1)] +
u[IX(i, j + 1)]
)
) / (1.0 + 4.0 * a)
v[idx] = (
v_prev[idx] +
a * (
v[IX(i - 1, j)] +
v[IX(i + 1, j)] +
v[IX(i, j - 1)] +
v[IX(i, j + 1)]
)
) / (1.0 + 4.0 * a)
set_bnd(BoundaryType.VELOCITY_HORIZONTAL, u)
set_bnd(BoundaryType.VELOCITY_VERTICAL, v)
Finally, remember to update _process()
func _process(delta: float) -> void:
copy_velocity_to_prev() # Add this
diffuse_velocity(delta) # and this
copy_density_to_prev()
diffuse_density(delta)
copy_density_to_prev()
advect_density(delta)
fade_density(delta)
fade_velocity(delta)
queue_redraw()
Velocity should now spread just like density does
Note: I have modified the values of previously added modifiers: velocity_add_scale = 0.12 and velocity_diffuse_rate = 0.01. If you are following along and you effect is different from this recording - this might be the reason
The river flows
At this point we have:
- Velocity diffusion
- Density diffusion
- Density advection
What we are missing is velocity advection. That’s right - we want the velocity field to be able to move itself. As a reminder - advection means moving something along a vector field. In this case, we want the velocities themselves to follow the velocity field - this will simulate flow.
Luckily - we already know how to do that! It will be exactly the same approach as with density advection. We will go backwards in time, calculate where the velocity might have come from, find the source cells, calculate fractional weights to know how much each source cell should contribute and calculate the value with bilinear equations. There are two minor differences: one is that we did it on one array for density but for velocity we have two arrays. The other one is that for density we did not modify the velocity arrays so there was no need to snapshot them - but here we do. So we need a snapshot
I refer you back to density advection for details.
# When we advect density, we move density along the velocity field
# Here we do the same but with the velocity itself - so we move velocity along the velocity field
# This simulates sort of "momentum"
# It works exactly the same way as advecting density, with one caveat - we have to snapshot the
# velocity arrays and use the snapshot as source (because we modify the actual velocity array in this function,
# which doesn't happen for density)
# See "advect_density" function for details
func advect_velocity(delta: float) -> void:
var dt0 := delta * N
for j in range(1, N + 1):
for i in range(1, N + 1):
var idx := IX(i, j)
# Go backwards through time
# Notice we are using u_prev and v_prev here - we need stable values
# The rest of the function is the same as for advecting density
var x := i - dt0 * u_prev[idx]
var y := j - dt0 * v_prev[idx]
x = clamp(x, 0.5, N + 0.5)
y = clamp(y, 0.5, N + 0.5)
# Find cell neighborus for the point we arrive at
var i0 := int(floor(x))
var i1 := i0 + 1
var j0 := int(floor(y))
var j1 := j0 + 1
# Calculate fractional weights
var s1 := x - i0
var s0 := 1.0 - s1
var t1 := y - j0
var t0 := 1.0 - t1
# Approximate the values using bilinear interpolation
u[idx] = (
s0 * (t0 * u_prev[IX(i0, j0)] + t1 * u_prev[IX(i0, j1)]) +
s1 * (t0 * u_prev[IX(i1, j0)] + t1 * u_prev[IX(i1, j1)])
)
v[idx] = (
s0 * (t0 * v_prev[IX(i0, j0)] + t1 * v_prev[IX(i0, j1)]) +
s1 * (t0 * v_prev[IX(i1, j0)] + t1 * v_prev[IX(i1, j1)])
)
set_bnd(BoundaryType.VELOCITY_HORIZONTAL, u)
set_bnd(BoundaryType.VELOCITY_VERTICAL, v)
And update _process()
func _process(delta: float) -> void:
copy_velocity_to_prev()
diffuse_velocity(delta)
copy_velocity_to_prev() # Velocity advection -
advect_velocity(delta) # - put it here
copy_density_to_prev()
diffuse_density(delta)
copy_density_to_prev()
advect_density(delta)
fade_density(delta)
fade_velocity(delta)
queue_redraw()
For the effect, take a careful look at the recording below. Notice I only move my mouse in a small region, but the density spreads much further - that’s because the velocity also spreads in that direction now, making it possible for density to reach further.
Note: Some values were modified, velocity_add_scale = 0.2 and velocity_fade_rate = 0.02
Let’s get physical
Alright, this last step is what I believe to be the most difficult to understand part of this whole endeavour. It is also the part that makes our fluid simulation behave in a believable way.
First let’s discuss the problem we want to solve. We want to simulate a fluid that is incompressible. That means it should never compress (pile-up) or decompress (vanish). Unfortunately, both velocity diffusion and velocity advection introduce such artifacts. We need to get rid of them. So how do we do that?
We can fix this by introducing and calculating two helper scalar fields - divergence and pressure. I call them helper fields, because at no point will we display the values contained in them - they will just be used to fix the velocity array.
Because this piece is more complicated, I will introduce code in fragments and explain it along the way, instead of explaining first and then going through the entire function, like I did in other chapters.
So let’s start with just defining the function
func project_velocity() -> void:
var h := 1.0 / N # This is a scaling factor for our grid (basically the distance between one cell center and a neighbouring cell center)
Now let’s focus on divergence.
Divergence is the measure of how much density flows out and into the cell. Positive divergence means more density is flowing out of the cell than is flowing in - which means density is “appearing out of nowhere”. Negative divergence means the opposite - more density is flowing in than is flowing out - which means density is “disappearing”. What we want is divergence of 0, which means density is in balance.
So how do we calculate it? Using the velocity fields. It’s not overly complicated. Each cell has 4 neighbours - 2 on each axis. So for horizontal velocities, we need to look at the nearest left and right neighbour. For vertical - top and bottom. We can then take an average of those values to calculate divergence.
Ok so here’s how we do it in code
for j in range(1, N + 1):
for i in range(1, N + 1):
# We calculate divergence by observing the velocities of neghbouring cells
# Remember - we want to check if the fluid flowing in and out of this cell is in balance or not
# Velocity of the neighbours will tell us that. For vertical velocity we look one cell down and one cell up
# For horzontal velocity - one cell to the left and one to the right
divergence[IX(i, j)] = -0.5 * h * ( # Multiplying by 0.5 centers the difference because we are looking at two neighbours
# The negative sign is just for convenience - it matches the later equation of where pressure is added to the Poisson solve
u[IX(i + 1, j)] - u[IX(i - 1, j)] + # Horizontal velocity to the left and right
v[IX(i, j + 1)] - v[IX(i, j - 1)] # Vertical velocity to the up and down
)
pressure[IX(i, j)] = 0.0 # Initialize pressure to 0 at each cell, we will approximate it later
# As with any other field, we need to calculate the boundaries
# Both divergence and pressure are simple scalar fields, same as density, so we treat them the same
set_bnd(BoundaryType.SCALAR, divergence)
set_bnd(BoundaryType.SCALAR, pressure)
With divergence in hand we can now calcualte pressure. A pressure scalar field is simply a field that counteracts divergence. It answers the question “how much correction should be applied at a given point to counteract the divergence”. High pressure pushes velocity away, low pressure pulls it in.
How do we calculate it? With our old trick - approximation using Gauss-Seidel relaxation. We start with a guess (which is just 0) and continue iterating, averaging the pressure value by looking at the divergence in the current cell and the pressure of neighbouring cells. We do that for every cell and repeat 20 times, each repeat is working on new pressure values, thus coming nearer and nearer to the answer - until it is good enough.
Here’s the code
# This is the same familiar Gauss-Seidel relaxation equation that we have used before
# Calculating exact pressure is too expensive, so we approximate with Gauss-Seidel relaxation
# Iterating 20 times should be good enough
# Once we are finished, we get a pressure field to correct the effect of divergence
for k in range(20):
for j in range(1, N + 1):
for i in range(1, N + 1):
# New pressure at this cell = local divergence + average of neighbouring pressures
# Notice how the divergence stays static (this loop doesn't modify the divergence array)
# But the pressure array changes in each iteration
# After repeating that 20 times we arrive at a solution that is good enough - an approximation
pressure[IX(i, j)] = (
divergence[IX(i, j)] + # Local divergence
pressure[IX(i - 1, j)] + # Neighbouring pressures
pressure[IX(i + 1, j)] +
pressure[IX(i, j - 1)] +
pressure[IX(i, j + 1)]
) / 4.0 # Four neighbours so we divide to get an average
set_bnd(BoundaryType.SCALAR, pressure)
Now the pressure array contains the correction amount. We need to calculate a gradient of the pressure array though. A gradient of a scalar field is simply the information how the values of that field change when going through that field - it produces a vector field. Because we operate on a grid, we can cheat a little and simply calculate the gradient by looking at pressure values of neighbouring tiles. For any processed cell, simply look at pressure to the left and to the right to know how the horizontal velocity needs to be modified. Same for vertical velocity. Once we have the gradient, we subtract it from the velocity to remove the artifacts
Here’s the code
# Here we subtract the pressure gradient from velocity to eliminate the "unnatural" part
# The right hand side is the gradient calculation
# Technically, a gradient of a scalar field (which pressure is) is a vector field
# Such field contains vectors which point in the direction where pressure increases the fastest
# And the magnitude of those vectors tells us how steep that increase is
# We operate on a grid, so we can approximate this with neighbours - left and right for horizontal
# velocity and top and down for vertical velocity.
# We need to divide by 2h because the distance between the neighbours is 2h
# So we arrive at equation: (p[i+1] - p[i-1]) / (2h)
# After simplification it becomes: 0.5 * (p[i+1] - p[i-1]) / h
for j in range(1, N + 1):
for i in range(1, N + 1):
u[IX(i, j)] -= 0.5 * (pressure[IX(i + 1, j)] - pressure[IX(i - 1, j)]) / h
v[IX(i, j)] -= 0.5 * (pressure[IX(i, j + 1)] - pressure[IX(i, j - 1)]) / h
set_bnd(BoundaryType.VELOCITY_HORIZONTAL, u)
set_bnd(BoundaryType.VELOCITY_VERTICAL, v)
Remember to update _process(). The project_velocity() call should be inserted after velocity diffusion and advection (since both introduce divergence and we want to fix it each time). We should also disable velocity fading at this point as we shouldn’t be needing it anymore
func _process(delta: float) -> void:
copy_velocity_to_prev()
diffuse_velocity(delta)
project_velocity()
copy_velocity_to_prev()
advect_velocity(delta)
project_velocity()
copy_density_to_prev()
diffuse_density(delta)
copy_density_to_prev()
advect_density(delta)
fade_density(delta)
#fade_velocity(delta)
queue_redraw()
With this, our fluid simulation is complete. There are more things we could do, like add vorticity confinement, which would help keep the vortices from dissipating - but let’s end here.
(You might want to fiddle around with the parameters to get them to your liking)
So what now?
That’s it. You now understand how to do fluid simulations in a grid. How do you make use of it in gamedev? Sky is the limit! For example, I introduced a small spaceship pulled down by gravity and made it spew fluid-simulated fire from the engine
If you want to replicate it, checkout these two commits
- In this one I tried to make the ship use a flamethrower turret, but the effect wasn’t really that fun
- This one implements the rocket
This knowledge can also probably be used to get some nice fire effects with shaders. I’m not yet sure how to do it, but I will try to make everything we created here work entirely on the GPU by using a Compute Shader. After all, shaders work nicely with textures and what are textures if not arrays of floats? They should store our density, velocity, divergence and pressure fields nicely! Moreover, the current CPU-based simulation is not very performant, so getting this to work on a GPU might be a big step-up.
No idea how to do it yet, but I’ll try to figure it out! I’ll make a blog post about it if I manage to crack it :)
Thank you for reading. If you liked this and would like to share some support, consider buying me a coffee