The Sandpile Model - Part 1

Intro

In this short workshop, we implement a simple sandpile model on a 2D grid. The sandpile model is a cellular automaton that simulates the behavior of granular materials. Each cell in the grid can hold a certain number of “grains of sand”. When the number of grains in a cell exceeds a certain threshold, the cell “topples”, distributing its grains to neighboring cells. This process can lead to complex patterns and self-organized criticality.

Let’s define the sandpile model more concretely: - We have a 2D grid of size N x N, where each cell can hold an integer number of grains. - We start with an initial configuration of grains in the grid. - We repeatedly apply the following rules until no cell has more than 3 grains: - If a cell has 4 or more grains, it topples: - The cell loses 4 grains. - Each of the four orthogonal neighbors (up, down, left, right) receives 1 grain. - If a neighbor is outside the grid, the grain is lost (it falls off the edge). - The process continues until all cells have 3 or fewer grains.

We implement this model in Python and visualize the final state of the sandpile after it stabilizes. Let’s get started!

Simulation

A simple example

Before we write the simulation code, let’s illustrate the toppling process with a simple example. Consider a 3x3 grid where the center cell has 21 grains and all other cells are empty.

We can apply the toppling rules step by step to see how the grains are distributed. The center cell topples multiple times, distributing grains to its neighbors, which may also topple if they exceed the threshold. This process continues until all cells have 3 or fewer grains. By following this example, we can understand how the toppling mechanism works and how it leads to a stable configuration.

In the code, we model the grid as a 2D numpy array, where each element represents the number of grains in that cell. The toppling process is implemented as a loop that checks for cells with 4 or more grains, topples them 4 grains at a time, and updates the grid according to the distribution rules until it stabilizes.

This 4-grain at a time toppling algorithm is straightforward, but can be slow for large grids and many grains. Since we are not interested in the dynamics but only in the final state, we can optimize the process by toppling multiple of 4 grains at a time.

For example, if a cell has 21 grains, we can topple it 5 times in one step, the grid loses 4x5 grains and has 1 grain remaining, and distributes 5 grains to each neighbor. This optimization can significantly speed up the simulation.

We use the divmod() function in the numpy library to calculate how many times each cell topples and the remaining grains after toppling in a single step. The divmod() function takes three arguments: (1) the 2D grid array, (2) the divisor (4 in our case), and (3) an out argument that specifies where to store the results of the quotient (number of topples) and remainder (grains left after toppling). This allows us to efficiently update the grid in each iteration until it stabilizes.1

1 the out argument allows us to store the results directly into pre-allocated arrays rather than creating new ones. This prevents Python from repeatedly allocating and freeing memory, which drastically speeds up code that runs in large loops with intensive data processing.

We see that comparing with the 4-grain at a time toppling, the optimized algorithm can reduce the number of iterations significantly, especially when there are cells with a large number of grains. This allows us to efficiently simulate larger grids and more complex configurations.

Full Implementation

We are now ready to implement the sandpile simulation. Our simulation uses a large grid with a large number of grains in the center. This configuration leads to a complex pattern after stabilization.

We first create three functions:

  1. stabilize_sandpile(): This function takes the grid as input and apply the toppling rules until the sandpile is stable (no cell has more than 3 grains).
  2. add_grains_at_center(): This function adds a specified number of grains to the center cell of the grid.
  3. simulate_sandpile(): This function initializes the grid, adds grains to the center, and calls the stabilization function to get the final stable state of the sandpile.

We are now ready to run the simulation. This may take a few moments as our grid size and the number of grains can be quite large. We measure the time taken to complete the simulation to get an idea of its performance.

Visualization

We visualize the final state of the sandpile using a heatmap, where the color intensity represents the number of grains in each cell. This allows us to see the complex patterns that emerge from the toppling process.

We use the plotly library for the visualization. The plotly library provides interactive plotting capabilities, which allows us to explore the final configuration of the sandpile in more detail. For example, we can hover over cells to see the exact number of grains, or zoom in on specific areas of the grid to see finer details of the pattern.

Exercise

  • Try different initial configurations, such as adding grains to multiple cells or starting with a random distribution of grains. Observe how the final patterns change with different initial conditions.
  • Try hexagonal or triangular grids instead of square grids. Change the toppling rule accordingly. What kind of patterns emerge?

References

  1. Sandpile Model - Wikipedia
  2. Searching for a fast Abelian Sandpile implementation - Julien Palard: A repository with many implementations of the sandpile model in Python and C, including optimizations and visualizations.
  3. Abelian sandpile model - Rosetta Code: Sandpile model implemented in many different programming languages.