{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Volume exclusion\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This model is introduced in\n\nS. Motsch, D. Peurichard, From short-range repulsion to Hele-Shaw problem in a model of tumor growth, *J. Math. Biology*, Vol. 76, No. 1, 2017.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First of all, some standard imports. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport sys\nimport time\nimport math\nimport torch\nimport numpy as np \nfrom matplotlib import pyplot as plt\nimport sisyphe.models as models\nfrom sisyphe.display import scatter_particles\n\n\nuse_cuda = torch.cuda.is_available()\ndtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Repulsion force\n\nEach particle is a disk with a (fixed) random radius. The particles repel each other when they overlap. The force exerted by a particle located at $x_j$ with radius $R_j$ on a particle located at $x_i$ with radius $R_i$ is\n\n\\begin{align}F = -\\frac{\\alpha}{R_i} \\nabla_{x_i} U\\left(\\frac{|x_i - x_j|^2}{(R_i + R_j)^2}\\right),\\end{align}\n\nwhere the potential is\n\n\\begin{align}U(s) = -\\log(s) + s - 1\\,\\,\\text{for}\\,\\, s<1 \\,\\,\\text{and}\\,\\, U(s) = 0\\,\\, \\text{for}\\,\\, s>1.\\end{align}\n\nInitially, the particles are clustered in a small region with a strong overlapping. \n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N = 10000\nrmin = .1\nrmax = 1.\nR = (rmax-rmin)*torch.rand(N).type(dtype)+rmin\nL = 100.\nD0 = 20.\npos = (D0*torch.rand((N,2)).type(dtype)-D0/2)+torch.tensor([L/2,L/2]).type(dtype)\n\ndt = .1\n\nsimu = models.VolumeExclusion(pos=pos,\n                 interaction_radius=R,\n                 box_size=L,\n                 alpha=2.5,\n                 division_rate=0., \n                 death_rate=0.,                    \n                 dt=dt)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the simulation over 200 units of time using an adaptive time-step which ensures that the energy :attr:`E` of the system decreases.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 13\n\nframes = [0,1,2,3,4,5,10,30,50,75,100,150,200]\n\ns = time.time()\nscatter_particles(simu,frames)\ne = time.time()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Print the total simulation time and the average time per iteration. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print('Total time: '+str(e-s)+' seconds')\nprint('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note this funny behaviour: the particles are clustered by size! \n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Repulsion force, random births and random deaths\n\nSame system but this time, particles die at a constant rate and give birth to new particles at the same rate. A new particle is added next to its parent and has the same radius. \n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N = 10000\nrmin = .1\nrmax = 1.\nR = (rmax-rmin)*torch.rand(N).type(dtype)+rmin\nL = 100.\nD0 = 20.\npos = (D0*torch.rand((N,2)).type(dtype)-D0/2)+torch.tensor([L/2,L/2]).type(dtype)\n\ndt = .1\n\nsimu = models.VolumeExclusion(pos=pos,\n                 interaction_radius=R,\n                 box_size=L,\n                 alpha=2.5,\n                 division_rate=.3, \n                 death_rate=.3,                    \n                 dt=dt,\n                 Nmax = 20000)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the simulation over 200 units of time using an adaptive time-step which ensures that the energy :attr:`E <sisyphe.models.VolumeExclusion.E>` of the system decreases.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frames = [0,1,2,3,4,5,10,30,50,75,100,150,200]\n\ns = time.time()\nscatter_particles(simu,frames)\ne = time.time()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Print the total simulation time and the average time per iteration. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print('Total time: '+str(e-s)+' seconds')\nprint('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.8.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}