{ "metadata": { "name": "", "signature": "sha256:16301910a63c2a257ba7348839cb4784828f6370eeeaf705c517582e1b7a83ca" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Particle-In-Cell Simulations" ] }, { "cell_type": "heading", "level": 5, "metadata": {}, "source": [ "Created by Cameron Parvini for MAE 6286" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import numpy as np\n", "import sympy\n", "import pylab\n", "import random\n", "sympy.init_printing()\n", "\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d.axes3d import Axes3D\n", "from matplotlib import animation\n", "from JSAnimation.IPython_display import display_animation" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "def create_grid(L,d,x_i,n,nt):\n", " \n", " \"\"\" Creates a discrete grid based on our desired length and step size,\n", " then also initializes our location storage matrix and calculates the number\n", " of grid points in the directional arrays.\n", " \n", " Parameters\n", " ----------\n", " L : array of floats\n", " Desired distances from origin to the edge of the cube in each direction\n", " d : float\n", " Discrete grid step size\n", " x_i : float\n", " Initial position of particles in cartesian coordinates\n", " n : int\n", " Number of particles\n", " nt : int\n", " Number of time steps\n", " \n", " Returns\n", " -------\n", " x,y,z : arrays of floats\n", " Arrays containing all of the discrete grid locations and their\n", " true cartesian value\n", " loc : array of floats\n", " Storage matrix for the (x,y,z) location (0th dimension) for each \n", " particle (1st dimension)\n", " nx,ny,nz : ints\n", " Number of grid points in each direction\n", " X,Y,Z : Coordinate Matrices of floats \n", " 3D matrices (size nx,ny,nz) containing all of the discrete grid \n", " points and their respective x, y, or z value. Allows for conversion\n", " from index number to grid value.\n", " \"\"\"\n", " \n", " x = np.arange(-L[0],L[0],d)\n", " y = np.arange(-L[1],L[1],d)\n", " z = np.arange(-L[2],L[2],d)\n", " X,Y,Z = np.meshgrid(x,y,z)\n", " loc = np.zeros((3,n,nt))\n", " loc[:,:,0] = x_i[:,:]\n", " nx = int(2*L[0]/d)\n", " ny = int(2*L[1]/d)\n", " nz = int(2*L[2]/d)\n", " \n", " return x, y, z, loc, nx, ny, nz, X, Y, Z" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "def PICsimAll(v_i,dt,nt,m,Zi):\n", " \n", " \"\"\" Takes in the intial conditions of the problem to provide the velocity and \n", " location of every particle for each time step.\n", " \n", " Parameters\n", " ----------\n", " v_i : Matrix of floats \n", " Matrix that stores the cartesian initial velocity components \n", " (0th dimension) for each particle (1st dimension)\n", " dt : Array of float\n", " Time step for each particle type\n", " nt : int\n", " Number of time steps desired\n", " m : Array of float\n", " Mass of each particle species in the simulation\n", " Zi : float\n", " Average charge number for the ion species\n", " \n", " Returns\n", " -------\n", " vn : Matrix of floats \n", " Matrix that stores the cartesian velocity components (0th dimension) for\n", " each particle (1st dimension) throughout all of time (2nd dimension)\n", " loca : Matrix of floats\n", " Matrix that stores the cartesian location (0th dimension) for each \n", " particle (1st dimension) throughout all of time (2nd dimension)\n", " \"\"\"\n", " \n", " # Create our local holder copies\n", " \n", " vn = np.zeros((3,n,nt))\n", " loca = np.zeros_like(vn)\n", " rho = np.zeros_like(vn) # This will use the number density to hold rho\n", " vn[:,:,0] = v_i[:,:]\n", " pot = np.zeros((nx, ny, nz))\n", " N = np.zeros_like(X) # This will hold the number of particles in each index \n", " pd = np.zeros_like(pot) # Our holder matrix for recalculating phi\n", " Exi = np.ones_like(X)*1000 # These are the initial cartesian field holders\n", " Eyi = np.zeros_like(Y)\n", " Ezi = np.zeros_like(Z)\n", " Bxi = np.zeros_like(X)\n", " Byi = np.zeros_like(Y)\n", " Bzi = np.ones_like(Z)*0.1\n", " \n", " # To simplify the math we're about to do...\n", " den = 2*(dy**2*dz**2+dx**2*dz**2+dx**2*dy**2)\n", " \n", " for i in range(1,nt):\n", " \n", " # Get the new velocities\n", " vn[:,:,i] = intEOMAll(vn[:,:,i-1],m,i,dt,Exi,Eyi,Ezi,Bxi,Byi,Bzi,loca[:,:,i])\n", " \n", " # Use the old velocities to get the current location\n", " loca[:,:n/2-1,i] = vn[:,:n/2-1,i-1]*dt[0]\n", " loca[:,n/2:,i] = vn[:,n/2:,i-1]*dt[1]\n", " \n", " # Use the current locations, find their nearest index\n", " locInd = loca[:,:,i].astype(int)\n", " \n", " # Fancy slicing to increment by the effective charge\n", " N[locInd[:,:n/2-1]] += -1.\n", " N[locInd[:,n/2:]] += Zi*1.\n", " \n", " # Find rho!\n", " rho = N * e / d\n", " \n", " # Calculate the potential for all interior points given the current rho,\n", " # so that you can set the stage for the next time iteration\n", " pot[1:nx-1,1:ny-1,1:nz-1] = ( dx**2*dy**2*dz**2/(eps*den)*\\\n", " (rho[1:nx-1,1:ny-1,1:nz-1]) +\\\n", " dy**2*dz**2/(den)*(pd[2:nx,1:ny-1,1:nz-1]+pd[0:nx-2,1:ny-1,1:nz-1])+\\\n", " dx**2*dz**2/(den)*(pd[1:nx-1,2:ny,1:nz-1]+pd[1:nx-1,0:ny-2,1:nz-1])+\\\n", " dx**2*dy**2/(den)*(pd[1:nx-1,1:ny-1,2:nz]+pd[1:nx-1,1:ny-1,0:nz-2]) )\n", " \n", " # No boundary conditions specified, eh? Why not Neumann?\n", " pot[0,:,:] = pot[1,:,:]\n", " pot[nx-1,:,:] = pot[nx-2,:,:]\n", " pot[:,0,:] = pot[:,1,:]\n", " pot[:,ny-1,:] = pot[:,ny-2,:]\n", " pot[:,:,0] = pot[0,:,1]\n", " pot[:,:,nz-1] = pot[:,:,nz-2]\n", " \n", " # Finally, recalculate the Electric Field for the next time step\n", " Exi[1:nx-1,1:ny-1,1:nz-1] = -(pot[2:nx,1:ny-1,1:nz-1]-\\\n", " pot[0:nx-2,1:ny-1,1:nz-1])/(2*dx)\n", " Eyi[1:nx-1,1:ny-1,1:nz-1] = -(pot[1:nx-1,2:ny,1:nz-1]-\\\n", " pot[1:nx-1,0:ny-2,1:nz-1])/(2*dy)\n", " Ezi[1:nx-1,1:ny-1,1:nz-1] = -(pot[1:nx-1,1:ny-1,2:nz]-\\\n", " pot[1:nx-1,1:ny-1,0:nz-2])/(2*dz)\n", " \n", " return vn, loca\n", "\n", "def intEOMAll(v,m,t,dt,Exi,Eyi,Ezi,Bxi,Byi,Bzi,loca):\n", " \n", " \"\"\" Takes in the intial conditions of the problem to provide the velocity at the \n", " next time step\n", " \n", " Parameters\n", " ----------\n", " v : Matrix of floats \n", " Matrix that stores the cartesian velocity components from the previous \n", " step (0th dimension) for each particle (1st dimension)\n", " m : Array of float\n", " Mass of each particle species in the simulation\n", " t : int\n", " Current time step\n", " dt : Array of float\n", " Time step for each particle type\n", " Exi, \n", " Eyi, \n", " Ezi : 3D Matrix of floats\n", " Cartesian electric field components for each grid point at the current\n", " time step\n", " Bxi, \n", " Byi, \n", " Bzi : 3D Matrix of floats\n", " Cartesian magnetic field components for each grid point at the current\n", " time step\n", " loca : Matrix of floats\n", " Matrix that contains the cartesian location (0th dimension) for each \n", " particle (1st dimension) at the current time step\n", " \n", " Returns\n", " -------\n", " v2 : Matrix of floats \n", " Matrix that contains the cartesian velocity components (0th dimension) \n", " for each particle (1st dimension) at the next time step\n", " \"\"\"\n", " \n", " v2 = np.zeros_like(v)\n", " \n", " for i in range(0,n):\n", " \n", " # We have to grab the right m, q, and dt values for each species\n", " if i < n/2:\n", " q = -e\n", " mass = m[0]\n", " dthold = dt[0]\n", " else:\n", " q = e\n", " mass = m[1]\n", " dthold = dt[0]\n", "\n", " # Find the current E and B field values according to the particle's location\n", " # from the cartesian component matrices for each particle individually. \n", " # Notice that there is no need to save these values between steps, so we\n", " # need only use holder arrays!\n", " \n", " Ehold = np.array([Exi[loca[0,i],loca[1,i],loca[2,i]],\\\n", " Eyi[loca[0,i],loca[1,i],loca[2,i]],\\\n", " Ezi[loca[0,i],loca[1,i],loca[2,i]]])\n", " Bhold = np.array([Bxi[loca[0,i],loca[1,i],loca[2,i]],\\\n", " Byi[loca[0,i],loca[1,i],loca[2,i]],\\\n", " Bzi[loca[0,i],loca[1,i],loca[2,i]]])\n", " \n", " # This is all to build up to the final line, equations 1 & 4 discretized \n", " # and solved for the next velocity:\n", " v2[:,i] = (Ehold + np.cross(v[:,i],Bhold))*dthold*q/mass + v[:,i]\n", " \n", " return v2" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "def animateAll(frame):\n", " x = (loc_all[0,:n/2,frame]-L[0])\n", " y = (loc_all[1,:n/2,frame]-L[1])\n", " z = (loc_all[2,:n/2,frame]-L[2])\n", " ax.scatter(x,y,z,c='b',marker='o')\n", " x = (loc_all[0,n/2:,frame]-L[0])\n", " y = (loc_all[1,n/2:,frame]-L[1])\n", " z = (loc_all[2,n/2:,frame]-L[2])\n", " ax.scatter(x,y,z,c='r',marker='o')\n", " \n", " ax.set_xlim((np.ndarray.min(view[0,:,frame]), np.ndarray.max(view[0,:,frame])))\n", " ax.set_ylim((np.ndarray.min(view[0,:,frame]), np.ndarray.max(view[0,:,frame])))\n", " ax.set_zlim((np.ndarray.min(view[0,:,frame]), np.ndarray.max(view[0,:,frame])))\n", " \n", "def initAll():\n", " ax.scatter([],[],[],c='b',marker='o')\n", " ax.scatter([],[],[],c='r',marker='o')\n", " " ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "def initial_v(n):\n", " \n", " \"\"\" Creates randomized initial velocities for ions and electrons separately \n", " based on the desired velocity ranges.\n", " \n", " Parameters\n", " ----------\n", " n : int\n", " Number of particles\n", " \n", " Returns\n", " -------\n", " v_i : Matrix of floats \n", " Matrix that stores the cartesian velocity components (0th dimension) for\n", " each particle (1st dimension)\n", " \"\"\"\n", " \n", " v_i = np.zeros((3,n))\n", " \n", " for i in range(0,n/2-1):\n", " v_i[0,i] = random.randrange(1e5,5e5)*random.choice([-1,1]) # [m/s]\n", " v_i[1,i] = random.randrange(1e5,5e5)*random.choice([-1,1]) # [m/s]\n", " v_i[2,i] = random.randrange(1e5,5e5)*random.choice([-1,1]) # [m/s]\n", " \n", " for i in range(n/2,n-1):\n", " v_i[0,i] = random.randrange(1e3,5e3)*random.choice([-1,1]) # [m/s]\n", " v_i[1,i] = random.randrange(1e3,5e3)*random.choice([-1,1]) # [m/s]\n", " v_i[2,i] = random.randrange(1e3,5e3)*random.choice([-1,1]) # [m/s]\n", " \n", " return v_i" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "# Initial Conditions\n", "n = 20 # [particles]\n", "dt_e = 2.50e-11 # [s]\n", "dt_i = 5.0e-5 # [s]\n", "dt = np.array([dt_e,dt_i])\n", "nt = 200 # [steps]\n", "m_e = 9.1094e-31 # [kg]\n", "m_He = 1.6605e-27 # [kg]\n", "Zi = 1.0 # [Average ion charge number, He]\n", "m_tot = np.array([m_e,m_He])\n", "e = 1.609e-19 # [C]\n", "eps = 8.854e-12 # [F/m]\n", "\n", "# Grid Info\n", "dx = dy = dz = d = 5.0e-5 # [m]\n", "L = 75.*np.array([d, d, d]) # [m]\n", "x_i = np.zeros((3,n)) # [m]\n", "x, y, z, loc, nx, ny, nz, X, Y, Z = create_grid(L,d,x_i,n,nt)\n", "\n", "# Begin Math\n", "v_i_tot = initial_v(n)\n", "vt_all,loc_all = PICsimAll(v_i_tot,dt,nt,m_tot,Zi)" ], "language": "python", "metadata": {}, "outputs": [ { "ename": "IndexError", "evalue": "index 156 is out of bounds for axis 0 with size 150", "output_type": "pyerr", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mIndexError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;31m# Begin Math\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 21\u001b[0m \u001b[0mv_i_tot\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0minitial_v\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 22\u001b[1;33m \u001b[0mvt_all\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mloc_all\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mPICsimAll\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv_i_tot\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdt\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnt\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mm_tot\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mZi\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;32m\u001b[0m in \u001b[0;36mPICsimAll\u001b[1;34m(v_i, dt, nt, m, Zi)\u001b[0m\n\u001b[0;32m 60\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 61\u001b[0m \u001b[1;31m# Fancy slicing to increment by the effective charge\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 62\u001b[1;33m \u001b[0mN\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mlocInd\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 63\u001b[0m \u001b[0mN\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mlocInd\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mZi\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mIndexError\u001b[0m: index 156 is out of bounds for axis 0 with size 150" ] } ], "prompt_number": 24 }, { "cell_type": "code", "collapsed": false, "input": [ "# Animate All\n", "fig = pylab.figure(figsize=(8,6));\n", "ax = Axes3D(fig)\n", "ax.view_init(45,45)\n", "\n", "view = loc_all\n", "ax.set_xlim((np.ndarray.min(view[0,:,0]), np.ndarray.max(view[0,:,0])))\n", "ax.set_ylim((np.ndarray.min(view[1,:,0]), np.ndarray.max(view[1,:,0])))\n", "ax.set_zlim((np.ndarray.min(view[2,:,0]), np.ndarray.max(view[2,:,0])))\n", "ax.set_xlabel(r'X Position [m]')\n", "ax.set_ylabel(r'Y Position [m]')\n", "ax.set_zlabel(r'Z Position [m]')\n", "ax.set_title(r'20 Particle Simulation for Problem 3')\n", "\n", "anim = animation.FuncAnimation(fig, animateAll, nt,\\\n", " init_func = initAll, interval=50)\n", "display_animation(anim, default_mode='loop')" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "" ] } ], "prompt_number": 21 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 } ], "metadata": {} } ] }