{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Mills\n\nExamples of milling behaviours. \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 display_kinetic_particles\n\nuse_cuda = torch.cuda.is_available()\ndtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Milling in the D'Orsogna et al. model\n\nLet us create an instance of the attraction-repulsion model introduced in\n\nM. R. D\u2019Orsogna, Y. L. Chuang, A. L. Bertozzi, L. S. Chayes, Self-Propelled Particles with Soft-Core Interactions: Patterns, Stability, and Collapse, *Phys. Rev. Lett.*, Vol. 96, No. 10 (2006).\n\nThe particle system satisfies the ODE:\n\n\\begin{align}\\frac{\\mathrm{d}X^i_t}{\\mathrm{d}t} = V^i_t\\end{align}\n\n\\begin{align}\\frac{\\mathrm{d}V^i_t}{\\mathrm{d}t} = (\\alpha-\\beta|V^i_t|^2)V^i_t - \\frac{m}{N}\\nabla_{x^i}\\sum_{j\\ne i} U(|X^i_t-X^j_t|)\\end{align}\n\nwhere $U$ is the Morse potential \n\n\\begin{align}U(r) := -C_a\\mathrm{e}^{-r/\\ell_a}+C_r\\mathrm{e}^{-r/\\ell_r}\\end{align}\n\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N = 10000\nmass = 1000.\nL = 10. \n\nCa = .5\nla = 2.\nCr = 1.\nlr = .5\n\nalpha = 1.6\nbeta = .5\nv0 = math.sqrt(alpha/beta)\n\npos = L*torch.rand((N,2)).type(dtype)\nvel = torch.randn(N,2).type(dtype)\nvel = vel/torch.norm(vel,dim=1).reshape((N,1))\nvel = v0*vel\n\ndt = .01\n\nsimu = models.AttractionRepulsion(pos=pos,\n                 vel=vel,\n                 interaction_radius=math.sqrt(mass),\n                 box_size=L,\n                 propulsion = alpha,\n                 friction = beta,                        \n                 Ca = Ca,\n                 la = la,\n                 Cr = Cr,\n                 lr = lr,                        \n                 dt=dt,\n                 p=1,                        \n                 isaverage=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the simulation over 100 units of time and plot 10 frames. The ODE system is solved using the Runge-Kutta 4 numerical scheme. \n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frames = [0,1,2,3,4,5,10,40,70,100]\n\ns = time.time()\nit, op = display_kinetic_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": [
        "## Milling in the Vicsek model \n\nLet us create an instance of the Asynchronuos Vicsek model with a bounded cone of vision and a bounded angular velocity, as introduced in: \n\nA. Costanzo, C. K. Hemelrijk, Spontaneous emergence of milling (vortex state) in a Vicsek-like model, *J. Phys. D: Appl. Phys.*, 51, 134004\n\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N = 10000\nR = 1.\nL = 20.\n\nnu = 1\nsigma = .02\nkappa = nu/sigma\n\nc = .175\nangvel_max = .175/nu\n\npos = L*torch.rand((N,2)).type(dtype)\nvel = torch.randn(N,2).type(dtype)\nvel = vel/torch.norm(vel,dim=1).reshape((N,1))\n\ndt = .01"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We add an option to the target\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "target = {\"name\" : \"normalised\", \"parameters\" : {}}\noption = {\"bounded_angular_velocity\" : {\"angvel_max\" : angvel_max, \"dt\" : 1./nu}}\n\nsimu=models.AsynchronousVicsek(pos=pos,vel=vel,\n                 v=c,\n                 jump_rate=nu,kappa=kappa,\n                 interaction_radius=R,\n                 box_size=L,\n                 vision_angle=math.pi, axis = None,\n                 boundary_conditions='periodic',\n                 variant=target,\n                 options=option,\n                 sampling_method='projected_normal')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the simulation over 200 units of time and plot 10 frames. \n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = -1\n\nframes = [0, 10, 30, 50, 75, 100, 125, 150, 175, 200]\n\ns = time.time()\nit, op = display_kinetic_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
}