Sand without grains - PART 1 - Abelian Sandpile Model

About a month ago I stumbled upon a project with an interesting FX : Moving Sand Dunes.
Yes ... like in the movie 'Dune' but without the worms.
And with 3 weeks dev time.

UPDATE:
The commercial is finally online so I can show you where I did apply this technique !

For some reason this video keeps being moved from YouTube and eventually has been deleted. Don't ask me why. I'm kinda sick of re-linking it every month so here's a frame that hopefully should show has article picture in Blogger.

This is the only version I found on the Net.

During the first day I kept watching and re-watching this YouTube video.

The part I found really fascinating is around second 16, where you can see the leading edge of the falling sand moving up, and the pile of sand underneath accumulating trying to chase the leading edge. Strangely enough this results in something visually moving up , opposed to what we would expect from sand falling down. 

What causes this behavior in real life is the result of millions and millions of tiny sand rocks. Each little rock has no choice than falling when there's nothing else supporting it, roll down colliding with other rocks and eventually stop as soon as the inclination is not strong enough to win the friction with the other rocks. 

All these little rocks will then eventually accumulate until the whole system becomes stable. 

If I had to recreate this in Houdini exactly as I just described, it would take a huge amount of CPU power, RAM and weeks of simulation. Not to mention disk space. In other words, a lot of resources. And the result would probably never be remotely close, because even with a terrific render farm, I could never have THAT many particles.

Actually, how many particles are we talking about ?

Well, given that the average medium-grained sized sand is 0.375 mm in diameter, how many grains of sand can I fit ...

• in a handful ? ~400.000
• in a bath tub ? ~9.000.000.000
• one sand dune ? probably this amount

Houdini Sand Solver is NOT an option in this case.

Let's make a few considerations:

• looking at the video, we cannot see the single sand rocks ! All we see is a large surface deforming. So maybe we can skip the particles approach entirely and focus on the behavior of this large surface.
• naively put, the behavior of the sand dune looks 'simple', hence there must be a 'simple' mechanism to reproduce it.
• we don't have to be physically correct ! (unless you're animating a cube : proof).

After a few minutes on YouTube I found this video :

Is really that simple ? Well ... seems too good to be true.
Good news is that since the algorithm is so simple, it will take no time to figure out if we're going in the right direction !

Before we start, let me show you the kind of result we are aiming for.

What you're looking at is a flipbook generated in Houdini and the surface is a HeightField.

The sand starts from a veeery unstable configuration and slowly rearranges in a way that is stable. Note how the little mountain in the foreground becomes stable soon, cause its height is lower. Note how the sand accumulates at the base, and the walls slowly become flat while growing and covering underlying elements.

Sand Pile Algorithm in 2D

After watching that video, I Googled "sandpile algorithm" and found this : Abelian Sandpile Model

I started reading and after a few lines I got to this point :

I stopped reading here, and I was already on Houdini.

NOTE:

Houdini has an HeightField Erode solver which is super mega cool, but in the long 5 minutes I tried I couldn't reproduce the sandpile behaviour described below. Regardless, I still wanted to experiment this idea, primarily because it's really quick to implement and once I do it myself I can control every single aspect of it.

Cool, let's get started.

As always, let's work in 2D first. It'll be easy to switch to 3D later.
Let's pretend this is our initial sand distribution:

Each column is a stack of sand grains.
For instance the stack E has 2 grains piled , the stack G has 3, and so on.

The idea is that if a stack of sand is too tall compared to the adjacent stacks (say by 2 grains), one grain of sand would fall down, towards one of the adjacent lower stacks.

When this happens, the taller stack becomes 'one grain' shorter , while the shorter stack will become 'one grain' taller ! Well, believe it or not, this is pretty much it.

Iteration - 1
Let's analyze the stacks one by one to identify the 'unstable' stacks:

stack C is 2 grains taller than B

stack D is 2 grains taller than E

stack G is 2 grains taller than F

stack H is 3 grains taller than I

So ...

C = C - 1 and B = B +1 (one grain of sand moves from C to B)

D = D - 1 and E = E +1 (one grain of sand moves from D to E)

G = G - 1 and F = F +1 (one grain of sand moves from G to F)

H = H - 1 and I = I +1  (one grain of sand moves from H to I)

Which brings us to the next state:

Iteration - 2
Let's analyze the stacks one by one again:

stack B is now 2 grains taller than A

So ...

B = B - 1 and A = A +1 (one grain of sand moves from B to A)

Iteration - 3
If we analyze the stacks again we notice there's nothing left to do. This means our sand distribution has reached a stable state and we're done !

What is interesting is that if you sum all the heights in each of the steps above , they'll return the same number.

0+1+3+4+2+1+3+3+0 = 17

0+2+2+3+3+2+2+2+1 = 17

1+1+2+3+3+2+2+2+1 = 17

This is good. It means that sand is not magically appearing or disappearing (the mass of the system is stable, like in reality). We just managed to re-arrange the stacks in a stable distribution using a simple iterative algorithm.

Assumptions about our initial sand distribution:

• if there is sand, there is sand underneath it.

This distribution would work with the algorithm ! Sand is supported by more sand underneath

This distribution would NOT work with the algorithm ! Sand in the red area is not supported.

• our sand grains are arranged into a grid.

Given the two assumptions above, we can simplify the system considering only the grains at the surface and ignore the sand below. In other words we'll only consider the visible surface and ignore whatever is happening in the depth.

The ideal data structure for this approach turns out to be a simple grid. Each cell of the grid correspond to one virtual stack of sand rocks (seen from above) and the number in it corresponds to the height of each stack.

Sand Pile Algorithm in 3D

Let's set this initial 'sand' distribution.

Remember that each square is a stack of sand, and the number is the height of that stack.

Before proceeding we need to set two 'global' variables:

• Mass transfer (M) : how much sand will fall from an unstable stack to one of the adjacent stacks.
• Threshold (T) : the minimum height difference between the current stack and its adjacent stacks that will mark that stack as unstable (and cause an amount M of sand to be transferred to some adjacent lower stack).

For this example let's set T = 2 and M = 0.5

Now, let's identify the stacks whose adjacent stacks are lower by T.

We can identify 3 ! For instance, both '2's are surrounded by many '0's (2-0 = T!), and '4' is obviously way taller than it's neighbors by a value bigger than T.
'1' instead is not selected because the difference between it's height and the adjacent stack's height is smaller than T.

The next step is to transfer some sand (M) from those 3 stacks to some adjacent lower stack.

Let's check our options:

I marked the possible 'falling' directions with arrows. For instance '4' can fall in all the directions cause it's literally surrounded by lower stacks. So ... what direction shall we choose ? Let's pick one randomly.
In the pic above I randomly choose one arrow (the red one), as the 'falling' direction.

Pic A

Cool. So the first '2' will fall left, '4' will fall down, and the last '2' will fall right.
What does 'fall' mean ?
It means that we'll have to remove M=0.5 from the tall stack (see pic below) ...

... and transfer it to the selected adjacent lower stack (see pic below).

More specifically, in this last step for each stack we need to add M as many times as incoming arrows.

At this point we have a new state, and we need to run the same steps again and again, until all the stacks are stable.

Remember that the orange stacks will become smaller (M will be subtracted) and the green stacks will become bigger (M will be added).

Step 2 (unstable)

  

Step 3 (unstable)

   

Step 4 (unstable)

   

Step 5 (stable !!)

Why is the last state stable ? Because no stack is surrounded by stacks lower than itself by T.

In other words, sand rocks have no reason to fall anywhere cause they're all supported by solid neighbors.

Cool.

There's an important consideration to make.

You might be tempted to do the following:

For each stack A ...

• get the 4 neighbors
• if any neighbor is lower than A by (say) 2 then
• reduce A by (say) 0.5
• increase one of those 'low' neighbors by 0.5

Well this would NOT work because during one iteration you're modifying one neighbor. But in the next iteration that neighbor might become the current stack !

So .. we need to do this in 2 steps, as showed above in Sandpile Algorithm 3d.

STEP 1 - identifying the unstable stacks only

For each stack A ...

• if it's surrounded by any stack lower than T, choose one randomly ... 

For instance if the stacks Up and Left of A are 'low', we can randomly choose one, say Up. So we'll mark A with 'U'. (so the choices are, U, D, L, R for an unstable A stack).

If A is NOT surrounded by 'low' stacks, then we mark A as 'S' (for stable).

Pic B: This corresponds to Pic A (scroll up a few pics).

STEP 2 - sand transfer

Now we can loop over each stack again and actually perform the transfer.

More specifically.

For every stack A...

• if A is NOT marked with S, reduce A by M (remember we set M = 0.5).
  Or, simply put, if A is unstable has to drop some sand.

• if the up stack is marked D (down) it means that the north stack wants to transfer some sand to the stack to his south, which is A. So increase A by M.
• if the down stack is marked U (up) it means that the south stack wants to transfer some sand to the stack to his north, which is A. So increase A by M.
• if the left stack is marked R (right) it means that the left stack wants to transfer some sand to the stack to his right, which is A. So increase A by M.
• if the right stack is marked L (left) it means that the right stack wants to transfer some sand to the stack to his left, which is A. So increase A by M.

Pic C

Note that in both Step 1 and 2 we're modifying only the stack currently looped on (A). The neighbors are inspected, never modified. Which happens to be exactly how Houdini VEX works ! What an incredible coincidence (it's not, I did it on purpose since the beginning ha!).

Cool.
Let's summarize the whole thing in pseudo-code.

STEP 1 (pseudo code)

1. set T = 2 (T is the height threshold)
2. for each stack A
1. h = height of A
2. hU = A's UP neighbor height
3. hD = A's DOWN neighbor height
4. hL = A's LEFT neighbor height
5. hR = A's RIGHT neighbor height
6. fU = 1 if hU is smaller than h by at least T, 0 otherwise.  (or ... h-hU>=T) 
7. fD = 1 if hD is smaller than h by at least T, 0 otherwise.  (or ... h-hD>=T)
8. fL = 1 if hL is smaller than h by at least T, 0 otherwise.  (or ... h-hL>=T)
9. fR = 1 if hR is smaller than h by at least T, 0 otherwise.  (or ... h-hR>=T)
10. if (any of fU or fD or fL or fR is 1) 
1. then pick one randomly and tag A with that label (example if we choose fU, A is marked with 'U').
2. else, mark A with 'S' (as stable)

At this stage all the stacks should have a tag that can be 'U', 'D', 'L', 'R', 'S' (like in Pic B).

STEP 2 (pseudo code)

1. set M = 0.5 (M is the mass transfer between sand stacks)
2. for each stack A
1. h = A's height
2. u = label of A's UP neighbor (could be 'U', 'D', 'L', 'R', 'S')
3. d = label of A's DOWN neighbor (could be 'U', 'D', 'L', 'R', 'S')
4. l = label of A's LEFT neighbor (could be 'U', 'D', 'L', 'R', 'S')
5. r = label of A's RIGHT neighbor (could be 'U', 'D', 'L', 'R', 'S')
6. if (u is 'D') then h = h + M
7. if (d is 'U') then h = h + M
8. if (l is 'R') then h = h + M
9. if (r is 'L') then h = h + M
10. if (A's label is different than 'S') then h = h - M

Cool ! This is pretty much the whole idea.

PART 2  - Implementing a HeightField Pandpile Solver in Houdini with VEX

PART 3 - Implementing a HeightField Sandpile Solver in Houdini with OpenCL (coming a little later than soon)

This algorithm is already fast cause it involves very simple math, nevertheless if you want to have production quality (don't we all ? ) you might want to use a very large grid with a ton of detail. In that case, why not speed up the calculation by say ... 10 or 20 times converting the code into OpenCL ?