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

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.