{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Bands\n\nThe classical Vicsek model in a square periodic domain `is known `_ to produce band-like structures in a very dilute regime. These structures also appears in a mean-field regime. To showcase the efficiency of the SiSyPHE library, we simulate a mean-field :class:`Vicsek ` model with the target :meth:`max_kappa() ` and $10^6$ particles. \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\n\nuse_cuda = torch.cuda.is_available()\ndtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the parameters and create an instance of the Vicsek model. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import sisyphe.models as models\n\nN = 1000000\nL = 1.\ndt = .01\n\nnu = 3\nsigma = 1.\nkappa = nu/sigma\n\nR = .01\nc = .1\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\nsimu=models.Vicsek(pos=pos,vel=vel,\n v=c,\n sigma=sigma,nu=nu,\n interaction_radius=R,\n box_size=L,\n boundary_conditions='periodic',\n variant = {\"name\" : \"max_kappa\", \"parameters\" : {\"kappa_max\" : 10.}},\n options = {},\n numerical_scheme='projection',\n dt=dt,\n block_sparse_reduction=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that we are in a mean field regime... \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Nneigh = simu.number_of_neighbours()\n\nprint(\"The most isolated particle has \" + str(Nneigh.min().item()) + \" neighbours.\")\nprint(\"The least isolated particle has \" + str(Nneigh.max().item()) + \" neighbours.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the block sparse parameters to their optimal value. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fastest, nb_cells, average_simu_time, simulation_time = simu.best_blocksparse_parameters(40,100)\n\nplt.plot(nb_cells,average_simu_time)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the function which compute the center of mass of the system (on the torus).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def center_of_mass(particles):\n cos_pos = torch.cos((2*math.pi / L) * particles.pos)\n sin_pos = torch.sin((2*math.pi / L) * particles.pos)\n average_cos = cos_pos.sum(0)\n average_sin = sin_pos.sum(0)\n center = torch.atan2(average_sin, average_cos)\n center = (L / (2*math.pi)) * torch.remainder(center, 2*math.pi)\n return center" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us save the positions and velocities of 100k particles and the center of mass of the system during 300 units of time. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sisyphe.display import save\n\nframes = [50., 100., 300.]\n\ns = time.time()\ndata = save(simu,frames,[\"pos\", \"vel\"],[center_of_mass], Nsaved=100000, save_file=False)\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": [ "At the end of the simulation, we plot the particles and the evolution of the center of mass. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sphinx_gallery_thumbnail_number = 2\nf = plt.figure(0, figsize=(12, 12))\nfor frame in range(len(data[\"frames\"])):\n x = data[\"pos\"][frame][:,0].cpu()\n y = data[\"pos\"][frame][:,1].cpu()\n u = data[\"vel\"][frame][:,0].cpu()\n v = data[\"vel\"][frame][:,1].cpu()\n ax = f.add_subplot(2,2,frame+1)\n plt.quiver(x,y,u,v)\n ax.set_xlim(xmin=0, xmax=simu.L[0].cpu())\n ax.set_ylim(ymin=0, ymax=simu.L[1].cpu())\n ax.set_title(\"time=\"+str(data[\"frames\"][frame]))\n\ncenter = data[\"center_of_mass\"]\n\ncenter_x = []\ncenter_y = []\n\nfor c in center:\n center_x.append(c[0])\n center_y.append(c[1])\n\nf = plt.figure(1)\nplt.plot(data[\"time\"],center_x)\nplt.ylabel(\"x-coordinate of the center of mass\")\nplt.xlabel(\"time\")\n\n\nf = plt.figure(2)\nplt.plot(data[\"time\"],center_y)\nplt.ylabel(\"y-coordinate of the center of mass\")\nplt.xlabel(\"time\")\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are still in a mean-field regime. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Nneigh = simu.number_of_neighbours()\n\nprint(\"The most isolated particle has \" + str(Nneigh.min().item()) + \" neighbours.\")\nprint(\"The least isolated particle has \" + str(Nneigh.max().item()) + \" neighbours.\")" ] } ], "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 }