DLA_numpy_diffusionCA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4597fa9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2e945ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sheet_cmap(reverse=False):\n",
    "    # a nice colormap\n",
    "    one = \"#f4f0e8\".lstrip('#')\n",
    "    two = \"#383b3e\".lstrip('#')\n",
    "    one = tuple(int(one[i:i+2], 16)/255 for i in (0, 2, 4))\n",
    "    two = tuple(int(two[i:i+2], 16)/255 for i in (0, 2, 4))\n",
    "    if reverse:\n",
    "        one, two = two, one\n",
    "    colors = [one, two]\n",
    "    cmap = mcolors.LinearSegmentedColormap.from_list('sheet_cmap', colors, N=256)\n",
    "    return cmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b3fc332",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configuration params\n",
    "size = (2000,2000)\n",
    "n_particles = 800000\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "# field: 1 for moving particles, 0 otherwise\n",
    "field = np.zeros(size[0]*size[1], dtype=np.int8).reshape(size)\n",
    "particles_pos = np.unravel_index(np.random.choice(np.arange(size[0]*size[1]), n_particles, replace=False), size)\n",
    "field[particles_pos] = 1\n",
    "\n",
    "# sticked_field: 1 for sticked particles, 0 otherwise\n",
    "sticked_field = np.zeros(size[0]*size[1], dtype=np.int8).reshape(size)\n",
    "sticked_field[size[0]//2, size[1]//2] = 1\n",
    "\n",
    "#array view for margolous blocking of field\n",
    "field_b_even = field.reshape((-1, 2, field.shape[1]//2, 2))  \n",
    "field_blocks_even = field_b_even.transpose((0,2,1,3))\n",
    "\n",
    "field_b_odd = field[1:-1,1:-1].reshape((-1, 2, field.shape[1]//2-1, 2))  \n",
    "field_blocks_odd = field_b_odd.transpose((0,2,1,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c671650",
   "metadata": {},
   "outputs": [],
   "source": [
    "#DLA simulation\n",
    "n_iter = 10000\n",
    "for i in range(n_iter):\n",
    "    #randomly move particles\n",
    "    if i%2==0:\n",
    "        rng.shuffle(field_blocks_even, axis=3)\n",
    "        rng.shuffle(field_blocks_even, axis=2)\n",
    "    else:\n",
    "        rng.shuffle(field_blocks_odd, axis=3)\n",
    "        rng.shuffle(field_blocks_odd, axis=2)\n",
    "    \n",
    "    #check for sticking\n",
    "    stick_1_0 = (field[2:,1:-1]+sticked_field[1:-1,1:-1]==2) & (sticked_field[2:,1:-1]==0)\n",
    "    field[2:,1:-1][stick_1_0]=0\n",
    "    sticked_field[2:,1:-1][stick_1_0]=1\n",
    "    \n",
    "    stick_m1_0 = (field[0:-2,1:-1]+sticked_field[1:-1,1:-1]==2) & (sticked_field[0:-2,1:-1]==0)\n",
    "    field[0:-2,1:-1][stick_m1_0]=0\n",
    "    sticked_field[0:-2,1:-1][stick_m1_0]=1\n",
    "    \n",
    "    stick_0_1 = (field[1:-1, 2:]+sticked_field[1:-1,1:-1]==2) & (sticked_field[1:-1, 2:]==0)\n",
    "    field[1:-1, 2:][stick_0_1]=0\n",
    "    sticked_field[1:-1, 2:][stick_0_1]=1\n",
    "    \n",
    "    stick_0_m1 = (field[1:-1, 0:-2]+sticked_field[1:-1,1:-1]==2) & (sticked_field[1:-1, 0:-2]==0)\n",
    "    field[1:-1, 0:-2][stick_0_m1]=0\n",
    "    sticked_field[1:-1, 0:-2][stick_0_m1]=1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbdcba57",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Sticked particles: {np.count_nonzero(sticked_field)} - total particles: {np.count_nonzero(sticked_field) + np.count_nonzero(field)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "320314b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)\n",
    "ax.spines[['top', 'right', 'bottom', 'left']].set_visible(False)\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.set_xlim([0,size[0]])\n",
    "ax.set_ylim([0,size[1]])\n",
    "#ax.set_axis_off()\n",
    "ax.imshow(sticked_field*2+field, cmap=sheet_cmap())\n",
    "plt.savefig(\"./dla_fast_2.png\", dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe22cacf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}