Box2D Fluid (Part 2)

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.

Introduction

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.

Modifying applyLiquidConstraints

Basically, we’ll be splitting the applyLiquidConstraints method up into the following methods:

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).

Preparing the Simulation

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;
}

Calculating Pressures

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.

Calculating Forces

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.

Moving the Particles

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.

Modifying update

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.

Running the Simulation

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:

Single threaded (2.6k particles)

Multithreaded (5.2k particles)

Performance

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.