Sandpile Model Prototype

As promised; I recently whipped up a Non-Abelian Sandpile model written in C. It's a fairly simple code where there is a grid of cells and if the sand grains in a cell is greater than the threshold, then the sand in that cell topples to the neighbors. In an abelian model this sand is distributed equally to its surrounding neighbors. In the model I have chosen: this is distributed randomly to neighboring cells. There's still much to be done. For one I need to test the power law for this code. And play around to see what kind of patterns need to be generated. This will be done in Part 2 if ever.     

#include <stdio.h>

#include <math.h>

#include <stdbool.h>

#include <stdlib.h>

const int size1 = 50;

const int size2 = 50;

int amount = 300;

int xi = 25;

int yi = 25;

int threshold = 4;

int timeFinal = 500000;

int grid[size1][size2];

void generategrid(){

int i = 0, j = 0;

for(int i = 0; i < size1; i++){

for(int j = 0; j < size2; j++){

grid[i][j] = 0;

}

}

}

void printmatrix(){

int i = 0,j = 0;

for(int i = 0; i < size1; i++){

for(int j = 0; j < size2; j++){

printf(" %4d ", grid[i][j]);

}

printf("\n");

}

}

void Topple(int i, int j) {

int total;

if( (i < 0 || i >= size1) || (j < 0 && j >= size2) ) {// check for invalid grid coorinates

return;

}

if(grid[i][j] < threshold){ // do we really need to topple?

return;

}

total = threshold;

if(i==0 || j==0 || i==size1-1 || j==size2-1){

//on borders, sand may drop off the edges/corners

if(i==0 && j==0 ){ //bottom left corner

int r1 = rand()%(total);

grid[0][1] += r1;

total -= r1;

grid[1][0] += rand()%total;

}

else if(i==size1-1 && j==size2-1 ){ //top right corner

int r1 = rand()%(total);

grid[size1-2][size2-1] += r1;

total -= r1;

grid[size1-1][size2-2] +=rand()% total;

}

else if(i==size1-1 && j==0 ){ // bottom right corner

int r1 = rand()%(total);

grid[size1-2][0] += r1;

total -= r1;

grid[size1-1][1] += rand()%total;

}

else if(i==0 && j==size2-1){ // top left corner

int r1 =rand()%(total);

grid[0][size2-2] +=r1;

total -= r1;

grid[1][size2-1] += rand()%total;

}

else if(i==0 && (j!=0 || j!=size2-1)){ //left edge but not corners

int r1 = rand()%(total);

grid[i][j+1] += r1;

total -= r1;

int r2 = rand()%(total);

grid[i][j-1] += r2;

total -= r2;

grid[i+1][j] += rand()%total;

}

else if(j==0 && (i!=0 || i!=size1-1)){ //bottom edge but not corners

int r1 = rand()%(total);

grid[i+1][j] += r1;

total -=r1;

int r2 = rand()%(total);

grid[i-1][j] += r2;

total -= r2;

grid[i][j+1] += rand()%total;

}

else if(i==size1-1 && (j!=0 || j!=size2-1)){ //right edge but not corners

int r1 = rand()%(total);

grid[i][j+1] += r1;

total -= r1;

int r2 = rand()%(total);

grid[i][j-1] += r2;

total -= r2;

grid[i-1][j] += rand()%total;

}

else if(j== size2-1 && (i!=0 || i!=size2-1)){ //top edge but not corners

int r1 = rand()%(total);

grid[i+1][j] += r1;

total -= r1;

int r2 = rand()%(total);

grid[i-1][j] += r2;

total -= r2;

grid[i][j-1] += rand()%total;

}

}

else {

int r1 = rand()%total;

grid[i+1][j] += r1;

total -= r1;

int r2 = rand()%total;

grid[i-1][j] += r2;

total -= r2;

int r3 = rand()%total;

grid[i][j+1] += r3;

total -= r3;

grid[i][j-1] += total; // no sand lost in this case

}

grid[i][j] = 0;

//debuggers:

//printmatrix();

// printf("Press any key to continue\n");

//getchar();

return;

}

void CheckThreshold(){

int i, j = 0;

    int gridfinal;

for( i = 0; i < size1; i++){

for( j = 0; j < size2; j++){

if(grid[i][j]>=threshold){

                gridfinal = grid[i][j]-threshold;

do{

Topple(i,j);

Topple(i+1,j);

Topple(i,j+1);

Topple(i,j-1);

}while(grid[i][j]>gridfinal);

}

}

}

}

void addgrain(){

int count;

int t;

do{

grid[xi][yi] += amount;

CheckThreshold();

t += 1;

}while(t < timeFinal);

}

int main(){

generategrid();

addgrain();

printmatrix();

}

Update July 13: I'll test these when I get the chance. But I plan on improving the random addition of sand grains as well. Here are the updates to the code bellow. An important limitation to the new add grain feature is that is assumes uniform distribution.  

                             int rx = 4; // new var for boundary of sand addition

                int ry = 4; // new var for boundary of sand addition

grid[xi+rand()%(rx)][yi+rand()%(ry)] += amount; //mod for add grain

Update Dec 3. I realized there was an issue with the way I defined the toppling conditions. In alignment with sandpile I modded the code (but beware its not been tested yet).

I caught this when reading the rules for the basic abelian sandpile model with the assumed threshold of 4 grains of sand.

Any site ${\displaystyle (x,y)}$ with

${\displaystyle z(x,y)\ge 4}$

is unstable and can topple (or fire), sending one of its chips to each of its four neighbors:

${\displaystyle z(x,y)\to z(x,y)-4,}$

${\displaystyle z(x\pm 1,y)\to z(x\pm 1,y)+1,}$

${\displaystyle z(x,y\pm 1)\to z(x,y\pm 1)+1.}$

Obviously 4 was the threshold for the tile so the while condition should operate until the grid value is gridfinal= grid[i][j]-threshold.

and then the threshold value is distributed to its neighbors as seen with the fact that the topple function's total (which is the max number of grains it can distribute) = to the global variable: threshold.