The author claims that their code is optimized, yet still took 10 days to run. Unfortunately, they missed a key observation that would have cut their runtime by 1/4 or more.
When stablizing a sandpile that begins in a symmetric configuration, it can never enter a non-symmetric configuration. Therefore instead of computing the whole 10000x10000 pile, you can compute just one 5000x5000 quarter of it, and add a special case for the internal edges. They should not subtract 4 when overflowing, but instead subtract 3 (or 2 for the corner) because they will always trade grains with their mirror image.
Although it's a bit more complex, you can also improve upon that by tracking not a quarter of the grid, but instead 1/8th by using a triangle. When tracking a triangle the cells on the main diagonal of the square are not mirrored, but instead their neighbors are. This may not work quite as well with their parallelism setup though.
Well yes, using symmetry to reduce the amount of computation is a common trick.
For instance, programmers of reciprocal space density-functional codes are certainly familiar with symmetry groups and (ir)reducible k-point grids etc.
When stablizing a sandpile that begins in a symmetric configuration, it can never enter a non-symmetric configuration. Therefore instead of computing the whole 10000x10000 pile, you can compute just one 5000x5000 quarter of it, and add a special case for the internal edges. They should not subtract 4 when overflowing, but instead subtract 3 (or 2 for the corner) because they will always trade grains with their mirror image.
Although it's a bit more complex, you can also improve upon that by tracking not a quarter of the grid, but instead 1/8th by using a triangle. When tracking a triangle the cells on the main diagonal of the square are not mirrored, but instead their neighbors are. This may not work quite as well with their parallelism setup though.