# Box2D Fluid (Part 2)

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:

• `prepareSimulation` — Finds neighbors, resets pressure variables, etc…
• `calculatePressure` — Determines a particle’s pressure based on its proximity to neighboring particles
• `calculateForce` — 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).

#### 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
{
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
{
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:

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