Nature, in Code

The Abelian Sandpile Model

Yesterday, I came across an interesting article via Paul Smaldino on Twitter (@psmaldino) that talked about a simple model called the "Abelian Sandpile Model". The model is a classical model of "Self-organized criticality".

The model is simple enough to explain. Imagine a rectangular grid, where each cell can hold maximally three grains of sand. As soon as you try to put a fourth grain of sand on the pile in the cell, it will topple. Four grains will fall to the cell's adjacent cells - one to the top, one to the bottom, and to the left, and one to the right.

If you start adding sand grains to the same cell over and over again, you get very interesting dynamics, as the article describes very well.

For those of you who've read my book Nature, In Code, this model is very simple to implement. To see interesting dynamics, however, we have to implement the model on a rather large grid. Below, I chose a 300 x 300 grid - if you run it locally I suggest 600 x 600. The model calculations are still very fast with 90,000 cells, but the visual update is rather slow. To get around this visual performance update, I decided to update the visual only every 1,000 sand grains.

Take a look below and watch the sand pile get to 100,000 sand grains (or bigger, if you are patient enough).

(I'm pretty sure the way I draw the visuals using D3 is not the smartest way to do it, performance-wise. If you are a D3 expert and can help me turn this into a smooth animation, I'd love to hear from you!). Thanks to Todd Bodnar for pointing out that canvas drawing should be much faster for this. Also thanks to Boris Conforty for spotting an error in an earlier version of the code.

RELOAD SIMULATION STOP SIMULATION

Code

//////////////////////////////////////
// START DRAWING CODE
//////////////////////////////////////

var _data; // data from previous update

function draw_grid(data, colors) {
	var color_obj = {};
	for (var i = 0; i < colors.length; i+=2) {
		color_obj[colors[i]] = colors[i+1];	
	}
	var width = 600;
	var height = 600;
	var grid_length = data.length;
	var width_cell = width/grid_length;
	var height_cell = height/grid_length;

	var canvas = document.getElementById("grid")
	if (canvas == null) {
		canvas = document.createElement('canvas');
		canvas.id = "grid";
		canvas.width = width;
		canvas.height = height;
		document.getElementsByTagName('body')[0].appendChild(canvas);
	}

	var context = canvas.getContext("2d");

	function draw_cells(){
		for (var i = 0; i < grid_length; i++) {
			for (var ii = 0; ii < grid_length; ii++) {
				if (_data && _data[i][ii] === data[i][ii]) {
					continue;
				} 
				context.fillStyle = color_obj[data[i][ii]];
        		context.fillRect(i*width_cell, ii*height_cell, width_cell, height_cell);
			}
		}
	}
	draw_cells();
	if (!_data) {
		_data = [];
	}
	for (var i = 0; i < grid_length; i++) {
		_data[i] = data[i].slice();
	}
}

function update_grid(data,colors) {
	draw_grid(data, colors);
}


//////////////////////////////////////
// END DRAWING CODE
//////////////////////////////////////			
			
var grid_length = 300;
var grid = [];
var temp_grid = [];
var colors = ["0","#ffffff","1","#00bfff","2","#ffd700","3","#b03060"];
var population = [];

var start_i;
var start_ii;

var continue_drawing = true;
var grain_counter = 0;



function init_grid() {
	for (var i = 0; i < grid_length; i = i + 1) {
		grid[i] = [];
		for (var ii = 0; ii < grid_length; ii = ii + 1) {
			grid[i][ii] = 0;
		}
	}
	start_i = Math.round(grid_length / 2);
	start_ii = Math.round(grid_length / 2);
}

init_grid();

draw_grid(grid,colors);

run_and_draw();

function run_and_draw() {
	var cd = continue_drawing;
	// the global variable continue_drawing can be set via iframe
	// to stop the simulation
	while (cd) {
		run_time_step();
		if (grain_counter % 1000 == 0) {
			cd = false;
			update_grid(grid,colors);
			setTimeout( function(){
				run_and_draw();
			},1000);
		}	
	}
}


function run_time_step() {
    add_sand(start_i, start_ii);
    grain_counter++;
}
			
function add_sand(i,ii) {
    var grains = grid[i][ii];
    if (grains < 3) {
        grid[i][ii]++;
    }
    else {
        grid[i][ii] = grains - 3;
        if (i > 0) {
            add_sand(i - 1, ii);
        }
        if (i < grid_length - 1) {
            add_sand(i + 1, ii);
        }
        if (ii > 0) {
            add_sand(i, ii - 1);
        }
        if (ii < grid_length - 1) {
            add_sand(i, ii + 1);
        }
    }
}