Abelian sandpile model

Abelian sandpile model
You are encouraged to solve this task according to the task description, using any language you may know.

This page uses content from Wikipedia. The original article was at Abelian sandpile model. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

Implement the Abelian sandpile model also known as Bak–Tang–Wiesenfeld model. It's history, mathematical definition and properties can be found under it's wikipedia article.

The task requires the creation of a 2D grid of arbitrary size on which "piles of sand" can be placed. Any "pile" that has 4 or more sand particles on it collapses, resulting in four particles being subtracted from the pile and distributed among it's neighbors.

It is recommended to display the output in some kind of image format, as terminal emulators are usually too small to display images larger than a few dozen characters tall. As an example of how to accomplish this, see the Bitmap/Write a PPM file task.

Examples:

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 4 0 0 -> 0 1 0 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 6 0 0 -> 0 1 2 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0  0 0  0  0    0 0 1 0 0
0  0 0  0  0    0 2 1 2 0
0  0 16 0  0 -> 1 1 0 1 1
0  0 0  0  0    0 2 1 2 0
0  0 0  0  0    0 0 1 0 0

Rust

<lang rust>// Set image size. const DIM: usize = 16;

// This function outputs the result to the console using UTF-8 block characters. fn draw_pile(pile: &Vec<Vec<usize>>) {

   for row in pile {
       let mut line = String::with_capacity(row.len());
       for elem in row {
           line.push(match elem {
               0 => ' ',
               1 => '░',
               2 => '▒',
               3 => '▓',
               _ => '█'
           });
       }

       println!("{}", line);
   }

}

// This function creates a file called "output.ppm" in the directory the program was run, which contains // a colored image of the pile. fn write_pile(pile: &Vec<Vec<usize>>) {

   use std::fs::File; // Used for opening the file.
   use std::io::Write; // Used for writing to the file.

   // Learn more about PPM here: http://netpbm.sourceforge.net/doc/ppm.html
   let mut file = File::create("./output.ppm").unwrap();

   // We write the signature, image dimensions and maximum color value to the file.
   let _ = write!(file, "P3\n {} {}\n255\n", DIM, DIM).unwrap();

   for row in pile {
       let mut line = String::with_capacity(row.len()*6);
       for elem in row {
           line.push_str(match elem {
               0 => "125 0 25 ", // Background color for cells that have no "sand" in them.

               // Depending on how many particles of sand is there in the cell we use a different shade of yellow.
               1 => "125 80 0 ",
               2 => "186 118 0 ",
               3 => "224 142 0 ",

               // It is impossible to have more than 3 particles of sand in one cell after the program has run,
               // however, Rust demands that all branches have to be considered in a match statement, so we
               // explicitly tell the compiler, that this is an unreachable branch.
               _ => unreachable!() 
           });
       }

       let _ = write!(file, "{}", line).unwrap();
   }

}

// This is the main part of the algorithm, a simple, recursive implementation of the model. fn handle_pile(x: usize, y: usize, pile: &mut Vec<Vec<usize>>) {

   if pile[y][x] >= 4 {
       pile[y][x] -= 4;

       // We check each neighbor, whether they have enough "sand" to collapse and if they do, 
       // we recursively call handle_pile on them.
       if y > 0 {
           pile[y-1][x] += 1;
           if pile[y-1][x] >= 4 {handle_pile(x, y-1, pile)}}

       if x > 0 {
           pile[y][x-1] += 1;
           if pile[y][x-1] >= 4 {handle_pile(x-1, y, pile)}}

       if y < DIM-1 {
           pile[y+1][x] += 1;
           if pile[y+1][x] >= 4 {handle_pile(x, y+1, pile)}}

       if x < DIM-1 {
           pile[y][x+1] += 1;
           if pile[y][x+1] >= 4 {handle_pile(x+1, y, pile)}}

       // Uncomment this line to show every iteration of the program. Not recommended with large input values.
       //draw_pile(&pile);

       // Finally we call the function on the current cell again, in case it had more than 4 particles.
       handle_pile(x,y,pile);
   }

}

fn main() {

   use std::thread::Builder; // Used to spawn a new thread.

   /* Rust by default uses a 2Mb stack, which gets quickly filled (resulting in a stack overflow) if we use any value larger than
    * about 30,000 as our input value. To circumvent this, we spawn a thread with 32Mbs of stack memory, which can easily handle
    * hundreds of thousands of sand particles. I tested the program using 256,000, but it should theoretically work with larger
    * values too. 
    */

   let _ = Builder::new().stack_size(33554432).spawn(|| {
       // This is our 2D grid. It's size can be set using the DIM constant found at the top of the code.
       let mut pile: Vec<Vec<usize>> = vec![vec![0;DIM]; DIM];

       // We place this much sand in the center of the grid.
       pile[DIM/2 - 1][DIM/2 - 1] = 16;

       // We start the algorithm on the pile we just created.
       handle_pile(DIM/2 - 1, DIM/2 - 1, &mut pile);

       draw_pile(&pile)
       
       // Uncomment this to save the image to a file after the recursive algorithm has ended.
       //write_pile(&pile)
   }).unwrap().join();

}</lang>

Output:

                
                
                
                
                
       ░        
      ▒░▒       
     ░░ ░░      
      ▒░▒       
       ░