Message Boards Message Boards


Abelian Sandpile Model Simulation with C++ and LTemplate

Posted 10 months ago
5 Replies
19 Total Likes


Abelian Sandpile Model

We call this Mathematica function to retrieve some background information from Wiki page: WikipediaData["Abelian sandpile model"]

The Abelian sandpile model, also known as the Bak[Dash]Tang[Dash]Wiesenfeld model, was the first discovered example of a dynamical system displaying self-organized criticality. It was introduced by Per Bak, Chao Tang and Kurt Wiesenfeld in a 1987 paper.The model is a cellular automaton.

Also in NKS there is discussion about self-organized criticality and some simple demo codes in Mathematica. You may also find this Numberhile video on this topic interesting, typically it talks about where "Abelian" comes from based on group theory point of view.


I found a GitHub repo with a nice implementation of the model in C++ and it works quite well for large number ( ~ 1 million initial sand particles at origin). I think this can be a good example to show how to use LTemplate package to call this C++ package from Mathematica if you do not want to start from scratch using the basic LibraryLink C API's.


Follow this link to install LTemplate package and refer to the document page coming with this package to find the basic example and some type conversion information (well, not like Mathematica you need to specify the type by yourself). In particular we need mma::IntMatrixRef as return type and to transfer a 2D array back to Mathematica. Meanwhile, we also need ra = mma::makeMatrix<mint> (ra stands for raw array) to store the grid information during computation. This ra also replaces the double pointer in the original code to handle the array for plot later. Here is the code

#include <LTemplate.h>
#include <numeric>
#include <iomanip>

class sandpile {

    const int dx[4] = { 1,0,-1,0 };
    const int dy[4] = { 0,1,0,-1 };
    const int N_size = 781; 
    const int N_size2 = 390;

    struct L_coord
       unsigned int x;
       unsigned int y;

    mma::IntMatrixRef CreateSandpile(int n ) {

        auto ra = mma::makeMatrix<mint>(N_size,N_size);
       bool** V_sites = new bool*[N_size];                    
       bool** Move_ = new bool*[N_size];                      
       unsigned int** Odometer = new unsigned int*[N_size];

       for (int k = 0; k < N_size; k++)
         V_sites[k] = new bool[N_size];
         Odometer[k] = new unsigned int[N_size];
         Move_[k] = new bool[N_size];

         for (int i = 0; i<N_size; i++)
          Move_[k][i] = false; V_sites[k][i] = false;

       L_coord* walking_ = new L_coord[N_size*N_size];        
       V_sites[N_size2][N_size2] = true;
       Move_[N_size2][N_size2] = true;
       ra(N_size2,N_size2) = n;
       walking_[0].x = N_size2;
       walking_[0].y = N_size2;  
       unsigned int x = 0, y = 0;
       int top_ = 0;
       unsigned int lx = 0, ly = 0;
       unsigned int d = 0; 
       unsigned int max_of_top = 0;
       unsigned long N_of_Moves = 0;

       while (top_ >= 0)
         N_of_Moves += 1;
         if (max_of_top < top_) { max_of_top = top_; }

         x = walking_[top_].x;  y = walking_[top_].y;

         top_--; Move_[x][y] = false;

         d = (ra(x,y) >> 2);
         ra(x,y) = ra(x,y) - (d << 2);

         for (int k = 0; k < 4; k++)
          lx = x + dx[k];
          ly = y + dy[k];

          V_sites[lx][ly] = true;
          ra(lx,ly) += d;
          if (Move_[lx][ly] == false && ra(lx,ly) >= 4)
              walking_[++top_].x = lx; 
              walking_[top_].y = ly; 
              Move_[lx][ly] = true;

        for (int k = 0; k<N_size; k++)
         delete[] V_sites[k]; delete[] Odometer[k]; delete[] Move_[k];
       delete[] V_sites; delete[] Odometer; delete[] Move_; delete[] walking_;
            return ra;

I keep most part intact, only made several essential modification according to the document of LTemplate:

  • Include LTemplate head file
  • Create a class called sandpile which wraps around the CreateSandpile function to be imported
  • Return mma::IntMatrixRef object instead of Clock object used to print computation time cost

Compile and Load

In a Mathematica notebook, I run the following lines. First I need to load the package

<< LTemplate`
SetDirectory[NotebookDirectory[]] (*put sandpile.h in the same directory as the working nb*)

we need to tell Mathematica which class to load and to pick the right member function at the same time

template =LClass["sandpile",{
    LFun["CreateSandpile", {Integer}, {Integer, 2}]

If our compiler is set up correctly, with verbose log, we should compile the class to dll just like in VS developer command or xcode or some clang/gcc environment on Linux.

"ShellCommandFunction" -> Print, "ShellOutputFunction"->Print


The last line is the dll in the $UserBaseDirectory. Then we can load the dll with


Import the class

sp = Make["sandpile"]

and call the member function:


I omitted one ArrayTrim function in C++ code to simplify the process. We actually need it to trim margins with all zeros in a row or column. I choose to implement it quickly in Mathematica with much shorter lines than the c++ version, even though we might transferred many null points back to kernel wasted. Nevertheless,

cleanRow[arr_] := DeleteCases[arr, _?(Total@# == 0 &)]

For small set, the

With[{arr = cleanRow@Transpose@cleanRow@(sp@"CreateSandpile"[5000])}, 
 ArrayPlot[arr, ColorRules -> {0 -> Black, 1 -> Magenta, 2 -> Red, 3 -> Blue}]]


Growing complexity

ArrayPlot[arr,ColorRules->{0->Black,1->Magenta,2-> Red,3-> Blue}]]&,


I tested the code on 10 million initial points and no problem at all to display it nicely with ArrayPlot. The animation at the beginning is only shows a quarter of the whole graph due to symmetry. The color scheme is the same as the original post.

  • 1 Million particles at center at t-zero

    With[{arr = cleanRow@Transpose@cleanRow@(sp2@"CreateSandpile"[1000000])}, 
    ArrayPlot[arr,ColorRules -> {0 -> Black, 1 -> Magenta, 2 -> Red, 3 -> Blue}]]


  • 10 Million particles at center at t-zero (only show quarter) and the result is attached in csv format.

      ColorRules -> {0 -> White, 1 -> Green, 2 -> Purple, 
        3 -> RGBColor["#D4AF37"]}]



  • Use LibraryLink based package LTemplate to call a C++ function with one integer input and 2D integer output.
  • Demonstrate how to pass 2D array in C++ back to Mathematica
  • Use ArrayPlot to visualize tensor computed in C++ environment
  • Safely wrap Manipulate or ListAnimate around external library function
  • Mathematica implementation can also be found here
  • Analytic research on this topic's complexity
5 Replies

Very nice! I don't fully understand how your code works, but the images are impressive! Keep up the good work!

Nice to see LTemplate in action :-)

You can also use

  ColorRules -> {0 -> Black, 1 -> Magenta, 2 -> Red, 3 -> Blue}

Since LTemplate works with classes/objects, one could store the state in the object, then have more member functions: one to propagate the model, one to retrieve the current state, and one to reset the state and choose N_size

enter image description here -- you have earned Featured Contributor Badge enter image description here

Your exceptional post has been selected for our editorial column Staff Picks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

Thanks so much for the demonstration Shengui!

Best regards, Eleazar

Of course. Your package is great and I just scratch the surface in the demo. I am working on a more comprehensive version from here with better OOP design. I learned many useful things from the documentations and tutorial examples coming with your package.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract