Sandpile Animations

Roman E. Maeder, Wolfram Research

A simple simulation consists of a square array of cells, each holding a number of grains of sand. Any pile containing more than 3 grains becomes unstable and sheds one grain to each of its four immediate neighbours. These cells may then themselves contain too many grains and shed grains to their neighbours, and so on. The evolution of this system eventually terminates, and the result is independent of the order in which the redistribution process is carried out.
We look at the fascinating toppling sequences that arise from various initial conditions, leading to fractal images, and discuss efficient code for generating animations and images with millions of grains of sand.

Sandpile Algebra

A simple simulation consists of a square array of cells, each holding a number of grains of sand. Any pile containing more than 3 grains becomes unstable and sheds one grain to each of its four immediate neighbours. These cells may then themselves contain more than 3 grains and shed grains to their neighbours, and so on. The evolution of this system eventually terminates (with each cell holding at most 3 grains), and the result is independent of the order in which the redistribution process is carried out.

This model was popularized on Numberphile [https://www.youtube.com/watch?v=1MtEUErz7Gg], and there is a good contribution about sandpile addition in the Wolfram community [http://community.wolfram.com/groups/-/m/t/1058615].

Here we want to explore what happens if we add 4 grains to the central cell, let it settle, and then repeat.

The color coding used is

colors={0->Black,1->Brown,2->Darker[Yellow],3->LightBrown,_->Orange}

Efficient Grain Distribution

We need to write fairly efficient code if we want to explore millions of grains added. The cells become a two-dimensional array, and we aim to compute the new cell contents after a toppling operation in one step.

If a is our 2-dimensional array of cells, the grains of sand remaining on a cell are whatever is left after distributing any multiples of four present (Mod[a,4]) plus the contributions from the adjacent cells. Each cell contributes q=Quotient[a,4] grains to each of its neighbours. We can add them all together by shifting q to the left, right, up, and down, and then adding corresponding cells. This works best if there is a border of zeroes at the edges of the array, because we can the perform the shifts via a rotation of the array in the first or second dimension.

With the center cell that receives the original grains in the middle of the array (having and odd dimension), the array is symmetric and we can simplify the rotation of the inner dimension, RotateLeft[q,{0,1}], to Transpose[RotateLeft[Transpose[q]]]. This is indeed a simplification because Transpose[q] is equal to q, and we need the value RotateLeft[q] anyway. Here, then is our procedure that works for such symmetric layouts, to compute the effects of all cells with more than 3 grains spilling over into their neighbours.

topple[a_]:=Module[{aa},
	With[{q=BitShiftRight[a,2]},
		 {qq=RotateLeft[q]+RotateRight[q]},
		aa=BitAnd[a,3]+qq+Transpose[qq];
	];
	If[Max[aa[[-1]]]>0,
		aa=ArrayPad[aa,{{1,1},{1,1}}]
	]; 
aa]

We use a bit shift in place of a quotient, andBitAnd in place of the mod, for some extra speed. The test at the end extends the array by a new border of zeroes, in case the sand has spilled into the border.

The use of quotient and remainder has the additional benefit of toppling all multiples of four present at the same time, and not just one.

Let us start with 1000 grains in the center, surrounded by a border of zeroes, ...

a=CenterArray[1000,{3,3}]

... and evolve it until all cells contain less than 4 grains.

n=0;
While[Max[a]>3, a=topple[a];n++ ];
{n,Total[a,2],Max[a],Dimensions[a]}
ArrayPlot[a,PixelConstrained->Min[16,Ceiling[400/Length[a]]],
ColorRules->colors]

If we want 20000 grains, we can simply multiply the previous result by 20 and continue.
The temporary cell keeps us entertained during the possibly lengthy computations:

a=20a;
PrintTemporary[Dynamic[{n,Max[a],Dimensions[a]}]];
While[Max[a]>3, a=topple[a];n++; ];
{n,Total[a,2],Max[a],Dimensions[a]}
ArrayPlot[a,PixelConstrained->Min[16,Ceiling[600/Length[a]]],
ColorRules->colors]

Exploiting Symmetries

4-fold Symmetry (plus Reflections)

All arrays starting from a single pile have a four-fold rotational symmetry. For a large gain in efficiency we can restrict the computation to one quarter of all cells, such as the lower-right one. To avoid having to deal with the borders, we can keep one row/column left and above the center row/column, preceded by the usual border of zeroes.

The center is, therefore, always at position {3,3}, and the computation needs to be slightly altered to add the new border of zeros only at the bottom and the right, and restore the second row/column and zeros in the first row/column after each step.

toppleq[a_]:=Module[{aa},
With[{q=BitShiftRight[a,2]},{qq=RotateLeft[q]+RotateRight[q]},
aa=BitAnd[a,3]+qq+Transpose[qq];
];
If[Max[aa[[-1]]]>0,aa=ArrayPad[aa,{{0,1},{0,1}}]]; 
aa[[All,2]]=aa[[2]]=aa[[4]];
aa[[All,1]]=aa[[1]]=aa[[-1]];
aa]

To start, we need at least 5 rows, with the center at {3,3}, and 10^5 grains in this example:

a=CenterArray[10^5,{5,5}];

This loop restricts the temporary progress display to an update once every second, rather than an update for each iteration, which may cost a significant portion of the total computation time.

m=0;
PrintTemporary[Dynamic[{m,Max[a],Dimensions[a]},UpdateInterval->1,TrackedSymbols:>{}]];
While[Max[a]>3, a=toppleq[a];m++; ];
{m,Total[a,2],Max[a],Dimensions[a]}
ArrayPlot[a,PixelConstrained->Min[16,Ceiling[400/Length[a]]],
ColorRules->colors]

Recovering the full image requires reflections along the horizontal and vertical axes and some stitching.

restore[q_]:=
With[{v=Transpose[Join[Reverse[Take[q,{4,-1}]],Take[q,{3,-1}]]]},
Join[Reverse[Take[v,{4,-1}]],Take[v,{3,-1}]]
]

afull=restore[a];{Total[afull,2],Dimensions[afull]}

ArrayPlot[afull,PixelConstrained->Min[16,Ceiling[600/Length[afull]]],
ColorRules->colors]

Millions of Grains

With this method we can deal with 1 million grains in about an hour. This color scheme uses black for empty cells, brown for one grain, red for two grains, and yellow for 3 grains. Here is the resulting distribution of grains after all cells have settled.

Animations

Toppling Sequences

Why does the computation take so long? The area of the sand pile grows much quicker than its border, where extra sand dumped into the center needs to be deposited eventually. With most cells full, they will just push grains back and forth between them, and not much progress happens.

Here is an animation of the amazing action going on when a single grain is added to a settled pile of 32599 grains originally dumped onto the center square. A bright orange is used to highlight the unstable cells about to topple. Watching this reminds me a lot of the game of Life. It takes 2319 steps until the array settles.

init=;
DynamicModule[{array,m,run},
Dynamic[
If[Length[array]==0,array=ArrayPad[add[init],{{2,2},{2,2}}]];
If[run&&Max[array]>3,array=topple[array];m++;];
ArrayPlot[array,PixelConstrained->5,ColorRules->colors,
Epilog->{Inset[Toggler[Dynamic[run],{True->Button[Style[m,12],Appearance->"Pressed"],False->Button[Style[m,12]]}],{Left,Top},{-1.2,1.1}]}],SynchronousUpdating->False],
SaveDefinitions->True,
Initialization:>(m=0;array={};run=False;)
]

Preloaded Sandpiles

Much more stunning animation sequences can be found if we do not start with an empty array, but rather with having a large area of almost toppling cells, with 3 grains each. One single grain dropped into the middle cell then triggers quite a spectacular sequence of events.

This sequence of frames show an initial setup with a square of cells with 3 grains each, and 4 grains in the center.

Now we can add another 4 grains, and so on.

Here is an animation of what happens if we repeat the process by adding 4 more grains whenever the toppling sequence settles.

ctr=41;pad=5;
colors={0->Black,1->Brown,2->Darker[Yellow], 3->LightBrown, _->Orange};
DynamicModule[{array={},m=0,n,run=False},
Dynamic[
If[Length[array]==0,array=ArrayPad[ConstantArray[3,{ctr,ctr}],pad*{{1,1},{1,1}}];
n=Total[array,2];
];
If[run,
If[Max[array]>3,
array=topple[array];m++,
array=add[array,4];n=Total[array,2];m=0;Pause[0.1]]
];
ArrayPlot[array,PixelConstrained->Min[12,Ceiling[600/Length[array]]],ColorRules->colors,
Epilog->{Inset[Toggler[Dynamic[run],{True->Button[Style[n,12],Appearance->"Pressed"],False->Button[Style[n,12]]}],{Left,Top},{-1.2,1.1}]}],SynchronousUpdating->False,TrackedSymbols:>{array,run}],
SaveDefinitions->True,
Initialization:>(m=0;n=4;array={};run=False;)
]

I' ve produced a video, available here [https://youtu.be/4-UosqDf9Ck], with the evolution of a 60x60 square over the first 18 minutes.
Another good choice for an initial condition is a Manhattan disk, because the topplings occur only horizontally and vertically, so the natural metric is the Manhattan metric. Here is the video: [https://youtu.be/E4WbG6Ymzhw].
For video production, I output individual frames in PNG format using Export[], then use ffmpeg to assemble them into an MP4 movie suitable for upload to Youtube.

Hexagonal Sandpiles

The square arrangement is not the only way to have symmetric cells. Another pretty setup uses hexagons in a honeycomb arrangement. Here each cell has six neighbours, and so we topple multiples of 6, not 4, as before. Here is the result of repeatedly dumping 6 grains into the middle hexagon, then letting it settle.

Watch out for some videos of these hexagonal topplings on my Youtube channel soon: https://www.youtube.com/c/RomanMäder.

Style[URLShorten["https://www.youtube.com/c/RomanMäder"],Large]

Bonus Animation

DynamicModule[{array={},m,n,run=False,c=25,pix,chc={0->Darker[Red],1->Lighter[Red],2->Yellow, 3->White, _->Orange}},
Dynamic[
If[Length[array]==0,array=ArrayPad[Table[3Boole[Abs[i]<c||Abs[j]<c],{i,-3c,3c},{j,-3c,3c}],5];
pix=Min[8,Ceiling[600/Length[array]]]];
n=Total[array,2];m=0;
If[run,
If[Max[array]>3,
array=topple[array];m++,
array=Nest[add,array,4];n=Total[array,2];m=0;Pause[0.2]]
];
ArrayPlot[array,PixelConstrained->pix,ColorRules->chc,
Epilog->{Inset[Toggler[Dynamic[run],{True->Button[Style[n,12],Appearance->"Pressed"],False->Button[Style[n,12]]}],{Left,Top},{-1.2,1.1}]}],SynchronousUpdating->False,TrackedSymbols:>{array,run}],
SaveDefinitions->True,
Initialization:>(array={};run=False;)
]

Extra Code

topple[a_]:=Module[{aa},
With[{q=BitShiftRight[a,2]},{qq=RotateLeft[q]+RotateRight[q]},
aa=BitAnd[a,3]+qq+Transpose[qq];
];
If[Max[aa[[-1]]]>0,aa=ArrayPad[aa,{{1,1},{1,1}}]]; 
aa]

center[a_]:=Extract[a,(Dimensions[a]+1)/2]

add[a_,k_:1]:=a+CenterArray[k,Dimensions[a]]

colors={0->Black,1->Brown,2->Darker[Yellow],3->LightBrown,_->Orange};