{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Self-organized criticality | Sandpile model \n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this exercise we simulate the sandpile model of Bak, Tang and Wiesenfeld (1987) discussed in the presentation. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import random\n",
"from scipy.stats import binned_statistic"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"############ No need to modify ###################\n",
"\n",
"def initialize_grid(N, orig_dist = \"zeros\", sink_value=-10000):\n",
" \n",
" \"\"\" Initializes a grid of size (N+2)x(N+2). Note that the grid is \"padded\"; the uttermost rows/columns represent\n",
" sinks on the edges of the surface; their value is set to -10000 in order to prevent them from reaching \n",
" the critical threshold.\n",
" \n",
" orig_dist - distribution of the grains at the beginning\n",
" \"zeros\" - empty grid \n",
" \"random\" - random values 0-3\n",
" sink_value - initial value of the sinks on the edges of the grid \"\"\"\n",
" \n",
" if orig_dist == \"zeros\":\n",
" grid = np.zeros((N+2,N+2))\n",
" \n",
" elif orig_dist == \"random\":\n",
" grid = np.random.randint(0,3,(N+2,N+2))\n",
" \n",
" \n",
" #set sink values\n",
" grid[0,:] = [sink_value]\n",
" grid[:,0] = [sink_value]\n",
" grid[N+1,:] = [sink_value]\n",
" grid[:,N+1] = [sink_value]\n",
" \n",
" return grid\n",
"\n",
"\n",
"def increment_by_one(grid, N, method = \"random\"):\n",
" \n",
" \"\"\" Adds one grain of sand to the grid.\n",
"\n",
" method - \"random\" or \"centre\". Decides whether the grains are dropped to random locations \n",
" or to the middle of the grid. \"\"\"\n",
" \n",
" if method == \"random\":\n",
" \n",
" x = random.randint(1,N)\n",
" y = random.randint(1,N)\n",
" grid[(x,y)] += 1\n",
" \n",
" elif method == \"centre\":\n",
" \n",
" x = int((N+2)/2)\n",
" y = int((N+2)/2)\n",
" grid[(x,y)] += 1\n",
" \n",
" return grid\n",
" \n",
"\n",
" \n",
"def logplot(sizes, xlabel, ylabel, num_bins):\n",
" \n",
" \"\"\" plots a log-log-pdf \"\"\"\n",
" \n",
" fig = plt.figure()\n",
" ax = fig.add_subplot(111)\n",
" bins = np.logspace(np.log10(1), np.log10(np.max(sizes)), num=num_bins)\n",
" pk, bin_edges = np.histogram(sizes, bins=bins, density=True)\n",
" bincenters, _, _ = binned_statistic(sizes, sizes, statistic='mean', bins=bins)\n",
" ax.loglog(bincenters, pk, marker = \"o\", markersize = 2)\n",
" ax.set_xlabel(xlabel)\n",
" ax.set_ylabel(ylabel)\n",
" plt.show()\n",
" \n",
" \n",
" \n",
"##############################\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We add grains one by one onto a NxN grid. A visualization of this process can be found in the lecture slides (page 8).\n",
"A pile becomes unstable once it reaches a critical threshold. This causes it to \"topple\", i.e. it gives one grain\n",
"to each of its four neighbors (diagonal neighbors are ignored), thus losing four grains. \n",
"This might cause the neighbors to become unstable, which leads to an avalanche. Note that one site can topple several times during one avalanche. A new grain is added only after the avalanche has ended, i.e. the height of all piles has once again fallen below the critical threshold.\n",
"\n",
"The functions for initializing the grid and dropping a new grain are given above; your task is to implement the loop in which the unstable piles redistribute the grains and the avalanches form. \n",
"\n",
"Keep track of the number of topplings (i.e. times a pile becomes unstable) after each grain is dropped. In addition, keep track of the duration of each avalanche. If, for example, a new grain causes one toppling, this toppling causes three topplings, and these three topplings cause altogether 5 topplings, the duration of this avalanche is three.\n",
" \n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n",
"def run_simulation(N, n_grains, threshold, orig_dist = \"zeros\", drop_method = \"random\"):\n",
" \n",
" \"\"\" Runs the simulation. \n",
" \n",
" N - length of the side of the grid, note, however, that initialize_grid() \n",
" returns a grid of size (N+2)x(N+2) to pad the surface with \"sinks\" where the grains can disappear\n",
" n_grains - number of grains dropped \n",
" threshold - critical threshold at which piles become unstable \n",
" orig_dist - distribution of grains on the grid at the beginning of the simulation;\n",
" \"zeros\" implies that there are no grains, \"random\" that each site has a random height between 0 and 3.\n",
" \n",
" drop_method - \"random\": grains dropped randomly, \"centre\": grains dropped to the centre of the grid\n",
" \n",
" Returns: \n",
" grid with the final heights\n",
" topplings - a list of the number of topplings after each new grain\n",
" durations - a list of the durations of the avalanches after each new grain\"\"\"\n",
" \n",
" grid = initialize_grid(N, orig_dist)\n",
" \n",
" topplings = []\n",
" durations = []\n",
" \n",
" for grain in range(n_grains):\n",
" \n",
" n_topplings = 0\n",
" duration = 0\n",
" \n",
" grid = increment_by_one(grid, N, drop_method) #a new grain is added\n",
" \n",
" #TODO: Implement a loop where the grains of a pile are redistributed if critical threshold is reached. \n",
" #Continue until the height of all piles falls under the threshold. Record the number of topplings and durations. \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" return grid, topplings, durations\n",
" \n",
" \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run the simulation and plot the pdfs of the topplings and the durations. Do the pdfs look like what you expected? Why/ why not? \n",
"\n",
"It is good to discard values from the beginning of the simulation before the critical state is reached. You can for example plot a time series of the variables to estimate the time needed to reach the critical state.\n",
"\n",
"Note that if you increase the number of added grains substantially, you might need to further \n",
"decrease the sink values in initialize_grid()."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"N = 30 #lenght of the side of the grid\n",
"n_grains = 10**4 #number of added grains\n",
"threshold = 4 #critical threshold\n",
"\n",
"grid, topplings, durations = run_simulation(N, n_grains, threshold, orig_dist = \"random\", drop_method = \"random\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"start_val = \n",
"\n",
"logplot(topplings[start_val:], \"number of topplings\", \"P(number of topplings)\", num_bins = 15)\n",
"logplot(durations[start_val:], \"duration\", \"P(duration)\", num_bins = 10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Repeat the simulation with dropping the grains to the centre. Is the end result affected? What about the time taken to reach the critical state? Why?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Optional: If you were to record the radius of each avalanche (i.e. the maximum distance from a toppled site to the site where the new grain was added), what kind of distribution do you think you would see? Why?\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}