# Simulating fluwids

Nerd sniping. It doesn't apply only to mathematical puzzles, and the victim doesn't have to be another person. It's not uncommon to nerd snipe oneself with a nifty programming idea. Especially one that looks simple enough to turn into a working demo in a matter of hours, if not minutes.

This is a story of how I decided to have a working fluid simulation later in the evening.

This article uses MathML for rendering equations. If you don't see any, consider using a browser that can display maths, like Firefox.

## Divergence

The idea entered my consciousness as I was daydreaming about cells floating in liquids, to give Breedmatic a biological special effect. As I considered what makes things flow, I was struck with a realization: a simple flowing fluid can be modelled with one equation!

If we simulate the fluid flow as a vector field *F*, then any volume *V* inside it satisfies our first equation:

\iint_{S} F n̂ \cdot d S = 0

where *n̂* is the normal vector to the surface *S*.

In simpler words, we can slice the fluid into any volume *V*, but the outflow of fluid is always going to be 0.

This itself is derived from the **divergence formula**:

d i v (F) = lim_{V \to 0} \frac{\iint_{S} F n̂ \cdot d S}{V}

Divergence is essentially the measure of the amount of "fluid" coming out from a given volume. Except in this case, the volume is infinitely small, and divergence is defined at any *point* inside the vector field.

Can you see how simple and powerful the concept is? Fluid into = fluid out of. We could easily create a grid representation of a 2D fluid that upholds this property. Start with a trivial grid, then subdivide it cleverly… until a large and pretty flowing vector field is formed. Piece of cake, that shouldn't take more than 1 hour of work.

## The grid

I took inspiration for the computer representation of the map from wind forecasts. Like this example from meteo.pl:

Baltic Sea and neighboring lands overlaid with wind barbs between 5 and 25 knots. Most of them point East, but on the Eastern shore they turn North This map uses wind barbs to indicate wind strength and direction in each cell of the underlying grid.

If we use square cells, it becomes a straightforward computer analogue of a vector field, with one important difference: the cell size inside a vector field is infinitely small. Computers can't do that (easily), but it's good enough for me.

A grid with an arrow in every cell The grid keeps the direction of flow in every cell.

Our first equation still needs to be upheld here, and not just in an abstract form. The easiest way to do that is to use one of the many formulas approximating divergence on a grid that we can find online. They have the advantage of being simple to calculate, operating on just 4 neighboring cells. Instead of writing them down, I will show what they calculate:

A grid with 4 cells in a square. with distances between the start and the end of each arrow marked with a letter. a to d Caption reads "div=a+b+c+d" Divergence is approximated by 4 components, one for each field. Each says how much the flow points outward of the middle of the group of 4 cells, and then they are summed up.

Our goal is to keep the inwardness of the arrows equal to their outwardness. Simple, right?

Writing this project, I made a couple of mathematical assertions that I can't be bothered or am not able to prove. One of those is wrong, can you spot it before the grand reveal?

**Assertion 1** enter stage right.

That if divergence = 0 in every possible group of 4 neighboring cells on the grid, then it's not possible to find any set of cells where divergence ≠ 0.

I think it's reasonable, because groups of 4 neighboring cells are overlapping, and cover the entire fluid together with their divergence = 0 guarantee. If it's false, then divergence could appear on a larger scale without appearing on a smaller scale first.

*Correction:* There is a hidden mistake here, that nevertheless doesn't affect algorithms as described in this post. Imagine a grid of cells where those in odd rows flow left, and those in even rows flow right. Taking each pair, the divergence is 0, but take 3 of them, and there will be 2 of one kind and 1 of the other. Unbalanced, and diverging! A better formulation should state:

**Assertion 1 (corrected)**

If divergence = 0 in every possible group of 4 neighboring cells on the grid, then it's not possible to find any set of cells where divergence ≠ 0, if that set is made of non-overlapping groups of 4 neighboring cells.

If we consider only pieces with 0 divergence, then the collection of them obviously has 0 divergence too.

Thanks to Roman for pointing it out.

## Building the grid

We have an equation which is defined on a 4-group. How do we create a grid of size, say, 16×16 cells while maintaining the equation? Obviously, divide and conquer.

Start with a single 4-group. Insert some flow between the cells. Split it while maintaining the equation, and then insert some extra flow. Rinse, repeat.

A 2×2 grid with zero flows gains some flows, and then turns into a 4×4 grid with 4 big flows, one in each quadrant. Finally, the 4×4 grid gets a separate flow in each cell.

Unrolling this step by step, we start with a group of 4 undisturbed cells of fluid flow.

4 cells in rows, each with a zero symbol inside.

There is no flow, so divergence = 0 and our requirement is maintained. So far so good. But how do we introduce some flow into it? A fluid that doesn't flow makes for a boring simulation indeed.

Introducing stir. Notice that divergence will not be changed if the fluid keeps flowing in circles, no matter how fast, so let's do that:

4 cells in arrows chasing each other, marking that no arrow points outwards. Each arrow has a letter "a" to "d" next to it, and all are =0. Caption says "div=a+b+c+d=0" Flow in a stirred 4-group.

Great, now let's split it and cross the boundary from a 2×2 grid into a 4×4 world. We need new assertions here:

**Assertion 2**: if a single cell is split into 4 cells all of the same flow, the resulting 4-group has 0 divergence.

This seems straightforward to me. First, a single cell in our model cannot have any divergence. Fluid comes in, fluid goes out. Making it bigger doesn't change anything, the flow is still straight through.

**Assertion 3**: when a grid of 0 divergence can be split into 4-cells by turning each cell into 4 identical ones, the resulting divergence stays 0.

This one reaches a bit farther, but it's essentially the same principle: making a small thing bigger doesn't change the basic properties of the whole.

A grid of 4 arrows turns into a grid of 16 arrows, with each quadrant having only 1 kind of arrow. This is what the grid looks like after splitting.

With assertions stated, we can do the stirring again:

A 4×4 grid with curly arrows embracing each intersection point between grid lines. We stir around each middle point, each with a different force.

But now we hit another assertion!

**Assertion 4**: Stirring a 4-group does not affect the divergence of any overlapping 4-groups.

This one sounds less reasonable. Let's quickly check the basic cases with stir magnitude = a.

First, is overlapping on the side enough to disturb the neighbors? The divergence before is indexed with 0, and after applying neighbor's stir with 1.

Shows a grid made of two 4-groups overlapping with 2 cells, and marks the contribution of new arrows from one 4-group to the other 4-group as "a" and "-a". The green 4-group and the one to the left share 2 cells. The one to the left was stirred with force *a*.

d i v_{1} = d i v_{0} + a - (- a) = d i v_{0}

The extra flow comes in, but it leaves again. What about overlapping on the corner?

Shows a grid made of two 4-groups overlapping with 1 corner cell, and marks the contribution of new arrows from one 4-group to the other 4-group as 0.

d i v_{1} = d i v_{0} + 0 = d i v_{0}

Here, the flow is perpendicular to the center between the 4 cells, so it never even enters the equation. Great!

Assertion upheld, we have all the pieces to divide and conquer. Both splitting and introducing stir is ours to have. It's time to implement the algorithm and enjoy the results.

## The reckoning

After all this thinking, armed with the above knowledge and confidence, writing the code seemed like a formality.

Not so. As I progressed through stages of stirring and splitting, a red light furiously engaged. CODE RED! CODE RED! INTEGRATION TESTS FAILING. It was the 4x4 grid. And the 8×8 grid. The divergence scanner reported anomalies never seen on smaller scales. Manual checks revealed that the checks were not just due to inaccuracies that computer calculations always make. Divergence was in places almost as strong as the average cell's flow!

But how? Where did I go wrong?

After some head-wall interactions, I extracted the offending 4-group. Half of it was the edge, which I added with an unchangeable flow of 0, the other half came from splitting a cell with a positive flow equal *f*. This is what it looked like:

A diagram of 2 cells: one with an arrow, other gray and with a null symbol. This diagram turns into a bigger diagram where the left 4-group has the same arrow in each cell, and all cells to the right have zeros inside. Marked are diverge components on the edge 4-group, all of which except one are 0. Edge of the grid can never have any flow other than 0, so it's marked in gray.

d i v = a + b + c + d = - f \neq 0

As you can see, the divergence comes out to *-f*. This proves that *Assertion 3* is wrong. Splitting can affect flow after all, on the boundary between cells.

But maybe I could salvage it? Find a better way to split? Come up with necessary equations? So I have tried, but never found any set of equations I could write down. I had to give up hope and search for another way.

At that point, the clock struck midnight.

## Try again

Burning with even more nerd-snipedness than before, my attention turned to overanalysis. Accurate splitting got me in trouble. And the splitting comes from the way my grid looks. What are other representations of flows inside a vessel? Is there one where splitting flows is easy?

There is one. But it does not hold directions like a weather map. It resembles a series of barriers instead, tracking flows *between*, not *inside* cells. And divergence is easily calculated, too, for each region that would have been a cell in the previous model:

4 lines arranged in a square, an arrow pointing outwards crosses each. Arrows marked a-b.

d i v = a + b + c + d

Can it be stirred?

4 cells, with walls they share highlighted. The shared walls have arrows of equal size crossing them. Following the arrows brings us back to starting cell. When the stirring flow crosses into a cell, the same amount leaves the cell, leaving its divergence unchanged.

And what does splitting look like?

A square with sides marked a-d gets split with 2 lighter crossing lines in the middle. All lines are split into 2 as well, resulting in 4 squares sharing sides. Here, the added walls are marked in lighter gray. Arrows have been replaced by their values.

Almost good. The values are not defined, but they are also not horribly floating. Let's try limiting the unknowns by simply forcing the boundary values to be halved.

The same diagram as above: 4 cells sharing 4 walls. 8 walls which are not shared are each marked as 1/2 times a-d, with the same letter on both walls on one side.

The flows at 4 barriers inside still need to be calculated, but now they don't depend on anything outside this small piece. But now another piece of picture is revealed. We can calculate how much surplus fluid each subcell receives from the outside by comparing its 2 boundary edges:

The same 4 cell arrangement, except this time each external corner has an arrow from one wall to the other. Arrows are marked s₁-s₄. s₁ = (b - a)/2 and so on.

If we distribute the surplus across the subcells, we can calculate the flows between them! But it's not so easy either. Notice that we can stir inside our 4-group, and that won't affect how much surplus fluid each cell has. That gives us some freedom we don't necessarily need.

## Waterfall

One way to estimate flows between subcells is to look at the surplus, and consider the cells to be a waterfall. Where would the surplus flow? Downwards, from the highest point.

4 labels s₁-s₄ arranged on neighboring columns of different heights.

Calculating this is quite straightforward if we know which cell has the highest surplus. The trick to do that is arranging cells in a circle, where consecutive ones share a barrier: *abcd*, and *a* again. If we create an array of *abcdabc*, we can always find the index of the highest value and start a new array there.

Different slices of length 4 from the abcdabc array.

The actual algorithm practically writes itself, and is left as an exercise for the reader. It can be somewhat freeform, but if it eliminates divergence, it's good enough. My version is available in fluw.

## Simulation

The new approach knows the same tricks as the old one: stirring and splitting. And it can do them without mistakes, too. But the representation is different. It's a network of barriers:

4 squares made up of lines. Each shares 2 sides with its neighbors.

which is not so easy to use in simulations. Thankfully, we can transform it into a grid by summing up barrier flows as components of the flow vector.

v⃗ = a - c, b - d

And we can also run our previous divergence detector, which… detects no divergence!

Video of the fluid simulation.

Mission complete! After several evenings of work, I ended up with fluw, which you can run on your computer too.

*Author's note:* The events didn't happen exactly as I described them here, but I prefer to spare the readers from reality which is way more boring and less dramatic than that.

Comments