Somehow, I've always been fascinated by fluidsimulations. They live at the point where physics, programming and creativity meet. Their characteristics are so weirdly mesmerising: acting as a whole when applying a force, only to break up into smaller blobs when the force becomes too large. Fiercely splashing around before settling down into a calm puddle again, waiting for the next interaction.
https://twitter.com/peeke__/status/1090363535299891200 (Re)published onI've tried to simulate fluids many times before, but never really nailed it. The basic concepts are deceivingly simple, but the implementation so incredibly complex.
Exploring the lesser known corners of the web, I stumbled upon a paper titled Particlebased Viscoelastic Fluid Simulation, by Simon Clavet, Philippe Beaudoin, and Pierre Poulin. It took me quite some time to figure it out (partially anyways) but I finally got my simulation running. I wanted to share my implementation with you (the gist of it). If you're really interested you can always check out the source for this article's header graphic.
The implementation is done in JavaScript using Three.js for rendering.
The basic theory
If we look at the behaviour of a fluid more closely, is possesses a few important characteristics:
 A fluid attracts other fluid of the same type. This makes the fluid into a coherent whole and gives rise to surface tension effects. If there were no other forces, a fluid would form a perfect sphere.
 A fluid has a maximum density. Fluid flows from highdensity regions to lowerdensity regions (this is also called a pressure gradient).
 A fluid is usually (mostly) incompressible: squishing it at some point should make the fluid expand in other places.
If we imagine a fluid consisting of a whole bunch of particles (molecules) we can say that particles attract or repel other particles, based on some rules. We'll be looking into how the double density relaxation algorithm implements these rules.
The algorithm
First off, all particles have a bunch of forces acting on them. These forces accelerate or decelerate their velocities. All particles are moved to their new position according to their velocities (even when this would move them into other particles or outside of the container).
Then two densities are calculated for each particle. The density is based on the amount of neighbouring particles within a certain distance (the interaction radius).
The first density is used to calculate an overall attracting and repelling of neighbouring particles within the interaction radius. The second density (near density) is used solely to push particles that are too close to each other away. This attracting and repelling is called relaxation.
Concretely, the algorithm consists of three passes, in which we loop over all of the simulated particles. The passes consist of the following steps:
Pass 1
 Apply forces to the particles and advance their position based on their velocity
 Store the particles' position in a spatial hash map for a fast lookup in later steps
Pass 2
 Find interesting (close enough) neighbours
 Calculate pressure and nearpressure for the particles' location
 Apply relaxation
Pass 3
 Contain the particles within a boundary
 Calculate their new velocity
Performance considerations
Before we dive in, a performance consideration: since we'll be looping over ~1000–2000 particles, it's absolutely necessary that we use Typed Arrays. Typed arrays can only contain numbers of a predefined type. Because of this constraint the browser can heavily optimise access to them.
We'll track the following properties for each particle:
const state = {
x: new Float32Array(PARTICLE_COUNT), // x location
y: new Float32Array(PARTICLE_COUNT), // y location
oldX: new Float32Array(PARTICLE_COUNT), // previous x location
oldY: new Float32Array(PARTICLE_COUNT), // previous y location
vx: new Float32Array(PARTICLE_COUNT), // horizontal velocity
vy: new Float32Array(PARTICLE_COUNT), // vertical velocity
p: new Float32Array(PARTICLE_COUNT), // pressure
pNear: new Float32Array(PARTICLE_COUNT), // pressure near
g: new Float32Array(PARTICLE_COUNT), // 'nearness' to neighbour
mesh: [] // Three.js mesh for rendering
};
The mesh property has to be a regular array since it will store Three.js Mesh objects, which aren't numbers (of course it's in no way necessary to use Three.js to implement this algorithm).
## Pass 1
Apply forces and advance particle based on velocity
The first step in our algorithm is to update the particle velocities (vx
and vy
), based on the forces acting on the particles and advance their positions based on their new velocity. Before we do this though, it's important that we store the starting positions (oldX
and oldY
), because we'll need these values later on.
// Pass 1
for (let i = 0; i < PARTICLE_COUNT; i++) {
// Update old position
state.oldX[i] = state.x[i];
state.oldY[i] = state.y[i];
applyGlobalForces(i, dt);
// Update positions
state.x[i] += state.vx[i] * dt;
state.y[i] += state.vy[i] * dt;
// Update hashmap
const gridX = (state.x[i] / canvasRect.w + 0.5) * GRID_CELLS;
const gridY = (state.y[i] / canvasRect.h + 0.5) * GRID_CELLS;
hashMap.add(gridX, gridY, i);
}
The dt
parameter in this code is our time step in seconds. This should be around .0166
seconds in order to reach 60 frames per second. Since speed is measured in distance per second (distance / t
) we increase our particles position by adding the velocity, multiplied by the time step dt
.
The applyGlobalForces
method applies external forces to our particles, accelerating their velocity.
const applyGlobalForces = (i, dt) => {
const force = GRAVITY;
state.vx[i] += force[0] * dt;
state.vy[i] += force[1] * dt;
};
Acceleration is expressed as the velocity increase per second (v / t
). If we have calculated our acceleration, we can multiply it by the time step dt
to find the increase in velocity. Our velocity is therefor increased by acceleration * dt
.
You can calculate the acceleration by a force with the following formula: acceleration = force / mass
. Assuming all of our particles have a mass of exactly 1
, we can simply say the velocity increase equals force / 1 * dt
, or force * dt
.
Store the particles' position in a spatial hashmap
After updating the positions of our particles, we store their locations in a spatial hashmap.
A spatial hashmap is an effective way to index (and access) a large set of points. We divide our canvas into a grid (an array), where each grid cell is a bucket (also an array). In our case the grid is 54 by 54 cells. Once a frame we clear the hash map and distribute all of our particles over the buckets.
To calculate the bucket index under which to store a particle, we calculate:
const index = Math.round(cellX) + Math.round(cellY) * gridCellsInRow
Querying the hashmap uses the same index to return the bucket.
The spatial hashmap expects the 0, 0
coordinate to be in the top left corner and uses cells as units. Our simulation renders the 0, 0
coordinate in the center of the canvas and uses pixels. This means we have to do a conversion before add to or quering the hashmap. In our case this is done in the following way:
const gridX = (state.x[i] / canvasRect.w + 0.5) * GRID_CELLS;
const gridY = (state.y[i] / canvasRect.h + 0.5) * GRID_CELLS;
As an important performance optimisation we reuse our bucket arrays. Creating an array is an expensive operation if you're doing it about 175.000 times per second. That's why clear our existing bucket arrays using splice
instead.
You can find the specific implementation of the spatial hashmap on github.
Pass 2
Calculate densities and neardensities
We now have a bunch of particles, easily searchable through the spatial hash map. Their positions have been updated, but we haven't taken into account any of the behaviour which makes them behave like a fluid yet. That's exactly what we'll do in pass two.
// Pass 2
for (let i = 0; i < PARTICLE_COUNT; i++) {
const neighbours = getNeighboursWithGradients(i);
updateDensities(i, neighbours);
// perform double density relaxation
relax(i, neighbours, dt);
}
Find interesting (close enough) neighbours
To calculate the two densities, we first have to find all the neighbours within the set interaction radius. We'll calculate a gradient (g
) for each neighbour, signifying it's distance towards the particle.
In this context, a gradient means a value between 0
and 1
. If a neighbour has the exact same position as our particle i
, the gradient value equals 1
(very close). The value gradually approaches 0
while moving away from the particle until finally the gradient equals 0
(not close at all) at the distance equal to the interaction radius.
This all happens in the getNeighboursWithGradients
function:
const getNeighboursWithGradients = i => {
const gridX = (state.x[i] / canvasRect.w + 0.5) * GRID_CELLS;
const gridY = (state.y[i] / canvasRect.h + 0.5) * GRID_CELLS;
const radius = (INTERACTION_RADIUS / canvasRect.w) * GRID_CELLS;
const results = hashMap.query(gridX, gridY, radius);
const neighbours = [];
for (let k = 0; k < results.length; k++) {
const n = results[k];
if (i === n) continue; // Skip itself
const g = gradient(i, n);
if (g === 0) continue
state.g[n] = g; // Store the gradient
neighbours.push(n); // Push the neighbour to neighbours
}
return neighbours;
};
We query the the hashmap for all the buckets within a certain radius around the particle's position.
For each found neighbour, we calculate and store the gradient. Since the array with found neighbours will also contain the particle itself, we'll have to check for that case (if (i === n) continue
). Neighbours with a gradient of 0
are also discarded, since they don't have any effect. The gradients are stored and then the neighbours are returned.
Calculating the gradient is done with the following function:
const gradient = (i, n) => {
const particle = [state.x[i], state.y[i]]; // position of i
const neighbour = [state.x[n], state.y[n]]; // position of n
const lsq = lengthSq(subtract(particle, neighbour));
if (lsq > INTERACTION_RADIUS_SQ) return 0;
const distance = Math.sqrt(lsq)
return 1  distance / INTERACTION_RADIUS;
};
We first calculate the squared distance (lsq
, length squared) between the particle and its neighbour. From here we can easily find the distance by taking the square root of lsq
, but since that's an expensive operation we first check wether lsq
is smaller than the squared interaction radius (squared, since the distance is also squared). If lsq
is larger, the neighbour falls outside of the interaction radius so we just return 0
.
We then calculate the actual distance between the particles, which we will convert into the gradient (g
) by dividing the distance by the interaction radius. This will return 0
when the positions are the same, and 1
when their distance is equal to the interaction radius. You might have noticed this is the exact reverse from what we wanted our gradient value to look like, which is why we calculate 1  distance / INTERACTION_RADIUS
to invert the value before returning.
It's probably smart to memoize the gradient function's result somehow, since it will calculate an expensive square root many, many times per second.
Calculate pressure and nearpressure
Next up, calculating pressure (p
) and nearpressure (pNear
) to use in our relaxation function. We already collected all interesting neighbours and their gradient values. We are going to weigh them and sum them together, to find the density and neardensity.
The weighing of the densities is an essential part of the algorithm. We'll use two different kernels for each density (a kernel is a weighing function): density = g * g
and nearDensity = g * g * g
. This means particles with gradient values close to 0
, will barely count towards these densities, but values close to 1
will. For the neardensity this effect is even more pronounced, effectively acting at very close distances only.
The following updatePressure
function calculates and stores both pressures:
const updatePressure = (i, neighbours) => {
let density = 0;
let nearDensity = 0;
for (let k = 0; k < neighbours.length; k++) {
const g = state.g[neighbours[k]]; // Get g for this neighbour
density += g * g;
nearDensity += g * g * g;
}
state.p[i] = STIFFNESS * (density  REST_DENSITY);
state.pNear[i] = STIFFNESS_NEAR * nearDensity;
};
In this code, STIFFNESS
and STIFFNESS_NEAR
are two constant that determine how pronounced the pressure effect influences the fluid. Pressurenear is calculated by multiplying STIFFNESS_NEAR
and nearDensity
. pressureNear
will always be positive, which means it will always exert a repelling force. This enforces the incompressibility characteristic of the fluid.
For regular density, this works a little bit differently. Before multiplying with STIFFNESS
, we subtract REST_DENSITY
from the density. This means the pressure will be positive when the density is above the rest density, resulting in a repelling force. Below the rest density the pressure will be negative, resulting in an attractive force. Effectively, this will result in the movement of particles from highpressure to lowerpressure regions. Also, this part is responsible for all the surface tension effects you can see.
Apply relaxation
Finally, the heart of the algorithm, the relaxation step.
const relax = (i, neighbours, dt) => {
const pos = [state.x[i], state.y[i]];
for (let k = 0; k < neighbours.length; k++) {
const n = neighbours[k];
const g = state.g[n];
const nPos = [state.x[n], state.y[n]];
const magnitude = state.p[i] * g + state.pNear[i] * g * g;
const direction = unitApprox(subtract(nPos, pos))
const force = multiplyScalar(direction, magnitude);
const d = multiplyScalar(force, dt * dt);
state.x[i] += d[0] * .5;
state.y[i] += d[1] * .5;
state.x[n] += d[0] * .5;
state.y[n] += d[1] * .5;
}
};
In this step we loop over all neighbours and apply a force based on the pressures (p
and pNear
) we calculated in the previous step.
To calculate the force, we first calculate the unit vector pointing from our particle towards the neighbour particle. A unit vector always has a length of exactly 1
. If you multiply it with the magnitude
value, you end up with a force pointing from our particle towards the neighbour particle with a strength equal to magnitude
. In our calculation, the magnitude is calculated by weighing the pressures again: p * g + pNear * g * g
.
Since forces act on the acceleration of particles, but we're interested in the actual displacement, we multiple the force once by dt
to convert the acceleration to velocity increase, then multiply by dt
again to find the actual displacement (d
) of the particles position.
Finally we subtract half of this displacement from particle i
's position and add the other half to the neighbour particle. Should we have had a positive pressure, this means the particles are moved away from each other. If the pressure was negative, they would have attracted each other.
The reason we apply half of the displacement to both particles, is to ensure that we do not violate Newton's third law:
Failing to preserve this law will cause the fluid to become unbalanced, causing it to spontaneously start moving by itself.
To optimise the performance of this function I've used a unitApprox function, which is copied from Nick Vogt. It's a good enough approximation of the unit vector, avoiding an expensive square root calculation.
Pass 3
Contain the particles within a boundary
With that pass out of the way, all that's left to do is to make sure our particles are contained within our container, calculate new velocities for them and to update the actual Three.js mesh positions, to reflect the simulation.
// Pass 3
for (let i = 0; i < PARTICLE_COUNT; i++) {
// Constrain the particles to a container
contain(i, dt);
// Calculate new velocities
calculateVelocity(i, dt);
// Update
state.mesh[i].position.set(state.x[i], state.y[i], 0);
}
Containing particles is simply the act of resetting their position to the boundary of the container, should they have crossed it. Our container has the shape of a circle, for which this boundary check is relatively easy:
const contain = i => {
const pos = [state.x[i], state.y[i]];
if (lengthSq(pos) > canvasRect.radiusSq) {
const unitPos = unit(pos)
const newPos = multiplyScalar(clone(unitPos), canvasRect.radius)
state.x[i] = newPos[0]
state.y[i] = newPos[1]
}
}
First we check if the particle passed the boundary, using the squared lengths so we don't have to calculate a square root. If they did cross the boundary, we take the unit vector of the particle's position and multiply it by the containing circle's radius. This will effectively place it exactly at the edge of the circle.
Calculate their new velocity
Next, we calculate the new velocity for the particle:
const calculateVelocity = (i, dt) => {
const pos = [state.x[i], state.y[i]];
const old = [state.oldX[i], state.oldY[i]];
const v = multiplyScalar(subtract(pos, old), 1 / dt);
state.vx[i] = v[0];
state.vy[i] = v[1];
};
Part of what makes this algorithm fast, is that it doesn't have to keep track of the particles' velocity during each single interaction. With ~10002000 particles interacting with one another, you can imagine that would yield a lot of interactions.
Instead, we subtract the position we stored at the beginning of our simulation step from the current position, to get the actual distance the particle moved. We divide this by the time step dt
we used to calculate the new velocity. It might not be exactly as accurate as the real thing, but it looks pretty good!
Bonus step: unsticking particles from the boundary of the container
So if you would now go ahead and implement these steps, you'd notice something strange: particles tend to stack up against the boundaries of their container. After breaking my head wondering about why this happens, I think I figured it out:
Once a particle remains at the same position (e.g., the boundary) for two time steps, this effectively reset its velocity to zero. And because there are no particles outside of the container, there is no pressure pushing the particle away from the boundary, causing it to get stuck there. While particles eventually will be pulled away from the border by the attraction of other particles inside the container, this effect is noticeable and doesn't look very nice.
To fix this, I've modified the contain function a little bit:
const contain = (i, dt) => {
const pos = [state.x[i], state.y[i]];
if (lengthSq(pos) > canvasRect.radiusSq) {
const unitPos = unit(pos)
const newPos = multiplyScalar(clone(unitPos), canvasRect.radius)
state.x[i] = newPos[0]
state.y[i] = newPos[1]
const antiStick = multiplyScalar(
unitPos,
INTERACTION_RADIUS * dt
)
state.oldX[i] += antiStick[0]
state.oldY[i] += antiStick[1]
}
};
After constraining the position to the container, the old position (oldX
and oldY
) is moved towards the outside of the container, perpendicular to the boundary. After the velocity is updated, this results in a net velocity away from the boundary.

That's it, the fluid simulation. Please let me know what you think of it on Twitter! If you mention the article link, your comment will be included below.
Happy coding!