Thanks for sharing! Unit directions can also be created using:
UnitVector[3, #] & /@ Range[3]
But using permutations is also quite neat actually... I'm not sure if '4' should be the toppling value for all dimensions; I think it should be equal to the number of neighbours. So for squares: 2 in 1D, 4 in 2D, 6 in 3D et cetera. There are examples with hexagons (in 2D), where they take 6 as the toppling value...
I'm not quite sure why it is so slow though... maybe because sand is created because of the '4' thing.
Maybe you can explain me the difference between definition of the forms:
Sandpile /: Sandpile[a_List] + Sandpile[n_Integer] := ...
Sandpile[a_List] + Sandpile[n_Integer] ^:= ...
I'm not very cognisant on this...