Using the Task Parallel Library to multithread the fluid simulation.
Source: FluidPort_2.zip
Github: klutch/Box2DFluid
This is the second in a series of posts that will explain how to add a fluid simulation to your Box2D game. I will be using C#, XNA and the Farseer Physics Engine.
In this post, I’m going to show how to use the Task Parallel Library to run portions of the fluid simulation in parallel, which will boost its overall performance.
Before I get started, I feel like another disclaimer is in order: the code for my fluid simulation is the first multithreaded code I’ve ever written. So if you’re more experienced in this area than me, and find some errors, please post them in the comments below. Thanks!
The idea behind running the simulation in parallel is to take independent chunks of code and run them at the same time on multiple threads. The Task Parallel Library makes this really easy with its Parallel.For
method. We’re going to have to restructure some of the code we wrote in Part 1, but only a little bit.
Basically, we’ll be splitting the applyLiquidConstraints
method up into the following methods:
prepareSimulation
— Finds neighbors, resets pressure variables, etc…calculatePressure
— Determines a particle’s pressure based on its proximity to neighboring particlescalculateForce
— Uses pressure information to calculate a particle’s delta
moveParticle
— Updates the particle’s position (and cell if necessary)We’ll have to change the code in those methods to operate on a single particle, instead of looping through all of the active particles. We’ll let TPL’s For
method do the looping (with the exception of moveParticle
, explained later).
We’ll write the prepareSimulation
method first, but first remove the following code from applyLiquidConstraints
:
// Prepare simulation for (int i = 0; i < _numActiveParticles; i++) { int index = _activeParticles[i]; Particle particle = _liquid[index];
// Scale positions and velocities _scaledPositions[index] = particle.position * MULTIPLIER; _scaledVelocities[index] = particle.velocity * MULTIPLIER;
// Reset deltas _delta[index] = Vector2.Zero;
// Find neighbors findNeighbors(particle); }
Now create the following method:
// prepareSimulation private void prepareSimulation(int index) { Particle particle = _liquid[index];
// Find neighbors findNeighbors(particle);
// Scale positions and velocities _scaledPositions[index] = particle.position * MULTIPLIER; _scaledVelocities[index] = particle.velocity * MULTIPLIER;
// Reset deltas _delta[index] = Vector2.Zero; }
Since we’re splitting up the calculating of pressure and forces into two separate methods, lets move the pressure variables to the Particle
class. Add the following member variables to Particle
class in Particle.cs:
public float p;
public float pnear;
Since we’re storing the pressure in particles now, we need to reset those variables before every step to avoid accumulating them over time. So add the following to the end of prepareSimulation
:
// Reset pressures
_liquid[index].p = 0;
_liquid[index].pnear = 0;
Now, remove the following code from applyLiquidConstraints
:
// Calculate pressure float p = 0.0f; float pnear = 0.0f; for (int a = 0; a < particle.neighborCount; a++) { Vector2 relativePosition = _scaledPositions[particle.neighbors[a]] - _scaledPositions[index]; float distanceSq = relativePosition.LengthSquared();
//within idealRad check if (distanceSq < IDEAL_RADIUS_SQ) { particle.distances[a] = (float)Math.Sqrt(distanceSq); //if (particle.distances[a] < Settings.EPSILON) particle.distances[a] = IDEAL_RADIUS - .01f; float oneminusq = 1.0f - (particle.distances[a] / IDEAL_RADIUS); p = (p + oneminusq * oneminusq); pnear = (pnear + oneminusq * oneminusq * oneminusq); } else { particle.distances[a] = float.MaxValue; } }
And create the calculatePressure
method:
// calculatePressure private void calculatePressure(int index) { Particle particle = _liquid[index];
for (int a = 0; a < particle.neighborCount; a++) { Vector2 relativePosition = _scaledPositions[particle.neighbors[a]] - _scaledPositions[index]; float distanceSq = relativePosition.LengthSquared();
//within idealRad check if (distanceSq < IDEAL_RADIUS_SQ) { particle.distances[a] = (float)Math.Sqrt(distanceSq); //if (particle.distances[a] < Settings.EPSILON) particle.distances[a] = IDEAL_RADIUS - .01f; float oneminusq = 1.0f - (particle.distances[a] / IDEAL_RADIUS); particle.p = (particle.p + oneminusq * oneminusq); particle.pnear = (particle.pnear + oneminusq * oneminusq * oneminusq); } else { particle.distances[a] = float.MaxValue; } } }
Notice how the p
and pnear
local variables were removed, and replaced with references to the particle’s p
and pnear
.
Restructuring the code that calculates forces is a little bit more in-depth than what we’ve done so far. Up until now, we’ve been storing the forces in the _delta
array. We can’t do this with threaded code, because there’s no guarantee that multiple particles won’t be trying to update the same delta value at the same time. Luckily, TPL provides a solution to this problem. We create an array to hold the delta values, pass that into the calculateForce
method, store/modify the deltas, return the array, and then copy the accumulated deltas into our _delta
array. I’ll show how to do that once I get to describing the changes to the update
method.
Remove the following code from applyLiquidConstraints
:
// Calculate forces float pressure = (p - 5f) / 2.0f; //normal pressure term float presnear = pnear / 2.0f; //near particles term Vector2 change = Vector2.Zero; for (int a = 0; a < particle.neighborCount; a++) { Vector2 relativePosition = _scaledPositions[particle.neighbors[a]] - _scaledPositions[index];
if (particle.distances[a] < IDEAL_RADIUS) { float q = particle.distances[a] / IDEAL_RADIUS; float oneminusq = 1.0f - q; float factor = oneminusq * (pressure + presnear * oneminusq) / (2.0F * particle.distances[a]); Vector2 d = relativePosition * factor; Vector2 relativeVelocity = _scaledVelocities[particle.neighbors[a]] - _scaledVelocities[index];
factor = VISCOSITY * oneminusq * DT; d -= relativeVelocity * factor; _delta[particle.neighbors[a]] += d; change -= d; } } _delta[index] += change;
Create the calculateForce
method:
// calculateForce private Vector2[] calculateForce(int index, Vector2[] accumulatedDelta) { Particle particle = _liquid[index];
// Calculate forces float pressure = (particle.p - 5f) / 2.0f; //normal pressure term float presnear = particle.pnear / 2.0f; //near particles term Vector2 change = Vector2.Zero; for (int a = 0; a < particle.neighborCount; a++) { Vector2 relativePosition = _scaledPositions[particle.neighbors[a]] - _scaledPositions[index];
if (particle.distances[a] < IDEAL_RADIUS) { float q = particle.distances[a] / IDEAL_RADIUS; float oneminusq = 1.0f - q; float factor = oneminusq * (pressure + presnear * oneminusq) / (2.0F * particle.distances[a]); Vector2 d = relativePosition * factor; Vector2 relativeVelocity = _scaledVelocities[particle.neighbors[a]] - _scaledVelocities[index];
factor = VISCOSITY * oneminusq * DT; d -= relativeVelocity * factor; accumulatedDelta[particle.neighbors[a]] += d; change -= d; } } accumulatedDelta[index] += change;
return accumulatedDelta; }
p
and pnear
were replaced with particle.p
and particle.pnear
, and _delta
was replaced with accumulatedDelta
.
Remove the following code from applyLiquidConstraints
:
// Move particles for (int i = 0; i < _numActiveParticles; i++) { int index = _activeParticles[i]; Particle particle = _liquid[index];
particle.position += _delta[index] / MULTIPLIER; particle.velocity += _delta[index] / (MULTIPLIER * DT);
// Update particle cell int x = getGridX(particle.position.X); int y = getGridY(particle.position.Y);
if (particle.ci == x && particle.cj == y) continue; else { _grid[particle.ci][particle.cj].Remove(index);
if (_grid[particle.ci][particle.cj].Count == 0) { _grid[particle.ci].Remove(particle.cj);
if (_grid[particle.ci].Count == 0) { _grid.Remove(particle.ci); } }
if (!_grid.ContainsKey(x)) _grid[x] = new Dictionary<int, List<int>>(); if (!_grid[x].ContainsKey(y)) _grid[x][y] = new List<int>(20);
_grid[x][y].Add(index); particle.ci = x; particle.cj = y; } }
Now create the moveParticle
method:
// moveParticle private void moveParticle(int index) { Particle particle = _liquid[index]; int x = getGridX(particle.position.X); int y = getGridY(particle.position.Y);
// Move particle particle.position += _delta[index] / MULTIPLIER; particle.velocity += _delta[index] / (MULTIPLIER * DT);
// Update particle cell if (particle.ci == x && particle.cj == y) return; else { _grid[particle.ci][particle.cj].Remove(index);
if (_grid[particle.ci][particle.cj].Count == 0) { _grid[particle.ci].Remove(particle.cj);
if (_grid[particle.ci].Count == 0) { _grid.Remove(particle.ci); } }
if (!_grid.ContainsKey(x)) _grid[x] = new Dictionary<int, List<int>>(); if (!_grid[x].ContainsKey(y)) _grid[x][y] = new List<int>(20);
_grid[x][y].Add(index); particle.ci = x; particle.cj = y; } }
You can delete the applyLiquidConstraints
method now.
Add the following using
statement to the top of FluidSimulation.cs:
using System.Threading.Tasks;
Remove the call to applyLiquidConstraints
in update
.
First, we need to call prepareSimulation
using TPL’s Parallel.For
method. Add the following code to update
:
// Prepare simulation
Parallel.For(0, _numActiveParticles, i => { prepareSimulation(_activeParticles[i]); });
The structure of this method is similar to a normal for loop. Basically, we’re starting at index 0, ending at _numActiveParticles
, and calling prepareSimulation
with the active particle’s index as a parameter. (You can read Microsoft’s reference page for Parallel.For
here)
Now we do the same for the calculatePressure
method. Add the following to update
after the call to prepareSimulation
:
// Calculate pressures
Parallel.For(0, _numActiveParticles, i => { calculatePressure(_activeParticles[i]); });
Now for the tricky part. I mentioned above that the call to calculateForce
was a little bit more in-depth than the rest. Add the following code to update
after the call to calculatePressure
:
// Calculate forces
Parallel.For(
0,
_numActiveParticles,
() => new Vector2[MAX_PARTICLES],
(i, state, accumulatedDelta) => calculateForce(_activeParticles[i], accumulatedDelta),
(accumulatedDelta) =>
{
lock (_calculateForcesLock)
{
for (int i = _numActiveParticles - 1; i >= 0; i--)
{
int index = _activeParticles[i];
_delta[index] += accumulatedDelta[index];
}
}
}
);
The reference for this method can be viewed here
Basically, this code accumulates the delta values from the calls to calculateForce
and uses that to update _delta
. This is where my lack of knowledge on this subject really starts to show. From what I understand, the code that applies changes to the _delta
array is called after a “task” is complete, and it’s possible that multiple tasks will be ending at the same time. So we have to use a lock
to synchronize access to the _delta
array. To create the lock object, just add the following as a member variable to the FluidSimulation
class:
private object _calculateForcesLock = new object();
Finally, add the call to moveParticle
in update
after the call to calculateForce
:
// Update particle cells
for (int i = 0; i < _numActiveParticles; i++)
moveParticle(_activeParticles[i]);
We just use a normal loop for this, because the _grid
cannot be asynchronously updated. This could be split up into a moveParticle
and updateParticleCell
method, but in my experience the overhead of using TPL for something as simple as updating particle positions outweighed the benefit.
Everything should be ready to run now. Before running the simulation, you might want to increase the MAX_PARTICLES
value to 10,000 or so.
If you have multiple cores, you can see that the threaded code is now using more of your CPU by peeking at the Task Manager:
The performance increase from using TPL will vary depending on the number of cores available, but there’s other things you can do to increase the performance… I haven’t mentioned them until now, because they’re somewhat subjective. Increasing the radius and cell size will decrease the number of neighbors any given particle has, which means less time spent in calculatePressure
and calculateForce
. But it will also give the simulation a slightly less accurate appearance. I think it’s worth it, though.
Here are the settings I use:
public const float CELL_SIZE = 0.6f;
public const float RADIUS = 0.9f;
Another thing to keep in mind is that running the program attached to Visual Studio’s debugger slows the simulation down a lot.
Well, that’s enough for this post. Next post I’ll talk about adding gravity and resolving collisions with the Box2D geometry.