Page 1 of 2 12 Last
  • Jump to page:
    #1
  1. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0

    Optimizing PIC (Particle in Cell) Code


    Hello all.

    I'm a terrible programmer who is trying to learn.

    Code:
    #!/usr/bin/env python
    #
    # Electrostatic PIC code in a 1D cyclic domain
    
    from numpy import arange, concatenate, zeros, linspace, floor, array, pi
    from numpy import diff, sign, sin, cos, sqrt, random, histogram
    from numpy import *
    
    import matplotlib.pyplot as plt # Matplotlib plotting library
    
    try:
        import matplotlib.gridspec as gridspec  # For plot layout grid
        got_gridspec = True
    except:
        got_gridspec = False
    
    # Need an FFT routine, either from SciPy or NumPy
    try:
        from scipy.fftpack import fft, ifft
    except:
        # No SciPy FFT routine. Import NumPy routine instead
        from numpy.fft import fft, ifft
    
    def rk4step(f, y0, dt, args=()):
        """ Takes a single step using RK4 method """
        k1 = f(y0, *args)
        k2 = f(y0 + 0.5*dt*k1, *args)
        k3 = f(y0 + 0.5*dt*k2, *args)
        k4 = f(y0 + dt*k3, *args)
    
        return y0 + (k1 + 2.*k2 + 2.*k3 + k4)*dt / 6.
    
    def calc_density(position, ncells, L):
        """ Calculate charge density given particle positions
        
        Input
          position  - Array of positions, one for each particle
                      assumed to be between 0 and L
          ncells    - Number of cells
          L         - Length of the domain
    
        Output
          density   - contains 1 if evenly distributed
        """
        # This is a crude method and could be made more efficient
        
        density = zeros([ncells])
        nparticles = len(position)
        
        dx = L / ncells       # Uniform cell spacing
        for p in position / dx:    # Loop over all the particles, converting position into a cell number
            plower = int(p)        # Cell to the left (rounding down)
            offset = p - plower    # Offset from the left
            density[plower] += 1. - offset
            density[(plower + 1) % ncells] += offset
        # nparticles now distributed amongst ncells
        density *= float(ncells) / float(nparticles)  # Make average density equal to 1
        return density
    
    def periodic_interp(y, x):
        """
        Linear interpolation of a periodic array y at index x
        
        Input
    
        y - Array of values to be interpolated
        x - Index where result required. Can be an array of values
        
        Output
        
        y[x] with non-integer x
        """
        ny = len(y)
        if len(x) > 1:
            y = array(y) # Make sure it's a NumPy array for array indexing
        xl = floor(x).astype(int) # Left index
        dx = x - xl
        xl = ((xl % ny) + ny) % ny  # Ensures between 0 and ny-1 inclusive
        return y[xl]*(1. - dx) + y[(xl+1)%ny]*dx
    
    def fft_integrate(y):
        """ Integrate a periodic function using FFTs
        """
        n = len(y) # Get the length of y
        
        f = fft(y) # Take FFT
        # Result is in standard layout with positive frequencies first then negative
        # n even: [ f(0), f(1), ... f(n/2), f(1-n/2) ... f(-1) ]
        # n odd:  [ f(0), f(1), ... f((n-1)/2), f(-(n-1)/2) ... f(-1) ]
        
        if n % 2 == 0: # If an even number of points
            k = concatenate( (arange(0, n/2+1), arange(1-n/2, 0)) )
        else:
            k = concatenate( (arange(0, (n-1)/2+1), arange( -(n-1)/2, 0)) )
        k = 2.*pi*k/n
        
        # Modify frequencies by dividing by ik
        f[1:] /= (1j * k[1:]) 
        f[0] = 0. # Set the arbitrary zero-frequency term to zero
        
        return ifft(f).real # Reverse Fourier Transform
       
    
    def pic(f, ncells, L):
        """ f contains the position and velocity of all particles
        """
        nparticles = len(f)/2     # Two values for each particle
        pos = f[0:nparticles] # Position of each particle
        vel = f[nparticles:]      # Velocity of each particle
    
        dx = L / float(ncells)    # Cell spacing
    
        # Ensure that pos is between 0 and L
        pos = ((pos % L) + L) % L
        
        # Calculate number density, normalised so 1 when uniform
        density = calc_density(pos, ncells, L)
        
        # Subtract ion density to get total charge density
        rho = density - 1.
        
        # Calculate electric field
        E = -fft_integrate(rho)*dx
        
        # Interpolate E field at particle locations
        accel = -periodic_interp(E, pos/dx)
    
        # Put back into a single array
        return concatenate( (vel, accel) )
    
    ####################################################################
    
    def run(pos, vel, L, ncells=None, out=[], output_times=linspace(0,20,100), cfl=0.5):
        
        if ncells == None:
            ncells = int(sqrt(len(pos))) # A sensible default
    
        dx = L / float(ncells)
        
        f = concatenate( (pos, vel) )   # Starting state
        nparticles = len(pos)
        
        time = 0.0
        for tnext in output_times:
            # Advance to tnext
            stepping = True
            while stepping:
                # Maximum distance a particle can move is one cell
                dt = cfl * dx / max(abs(vel))
                if time + dt >= tnext:
                    # Next time will hit or exceed required output time
                    stepping = False
                    dt = tnext - time
                #print "Time: ", time, dt
                f = rk4step(pic, f, dt, args=(ncells, L))
                time += dt
                
            # Extract position and velocities
            pos = ((f[0:nparticles] % L) + L) % L
            vel = f[nparticles:]
            
            # Send to output functions
            for func in out:
                func(pos, vel, ncells, L, time)
            
        return pos, vel
    
    ####################################################################
    # 
    # Output functions and classes
    #
    
    class Plot:
        """
        Displays three plots: phase space, charge density, and velocity distribution
        """
        def __init__(self, pos, vel, ncells, L):
            
            d = calc_density(pos, ncells, L)
            vhist, bins  = histogram(vel, int(sqrt(len(vel))))
            vbins = 0.5*(bins[1:]+bins[:-1])
            
            # Plot initial positions
            if got_gridspec:
                self.fig = plt.figure()
                self.gs = gridspec.GridSpec(4, 4)
                ax = self.fig.add_subplot(self.gs[0:3,0:3])
                self.phase_plot = ax.plot(pos, vel, '.')[0]
                ax.set_title("Phase space")
                
                ax = self.fig.add_subplot(self.gs[3,0:3])
                self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]
                
                ax = self.fig.add_subplot(self.gs[0:3,3])
                self.vel_plot = ax.plot(vhist, vbins)[0]
            else:
                self.fig = plt.figure()
                self.phase_plot = plt.plot(pos, vel, '.')[0]
                
                self.fig = plt.figure()
                self.density_plot = plt.plot(linspace(0, L, ncells), d)[0]
                
                self.fig = plt.figure()
                self.vel_plot = plt.plot(vhist, vbins)[0]
            plt.ion()
            plt.show()
            
        def __call__(self, pos, vel, ncells, L, t):
            d = calc_density(pos, ncells, L)
            vhist, bins  = histogram(vel, int(sqrt(len(vel))))
            vbins = 0.5*(bins[1:]+bins[:-1])
            
            self.phase_plot.set_data(pos, vel) # Update the plot
            self.density_plot.set_data(linspace(0, L, ncells), d)
            self.vel_plot.set_data(vhist, vbins)
            plt.draw()
            
    
    class Summary:
        def __init__(self):
            self.t = []
            self.firstharmonic = []
            
        def __call__(self, pos, vel, ncells, L, t):
            # Calculate the charge density
            d = calc_density(pos, ncells, L)
            
            # Amplitude of the first harmonic
            fh = 2.*abs(fft(d)[1]) / float(ncells)
            
            print "Time:", t, "First:", fh
            
            self.t.append(t)
            self.firstharmonic.append(fh)
    
    ####################################################################
    # 
    # Functions to create the initial conditions
    #
    
    def landau(npart, L, alpha=0.2):
        """
        Creates the initial conditions for Landau damping
        
        """
        # Start with a uniform distribution of positions
        pos = random.uniform(0., L, npart)
        pos0 = pos.copy()
        k = 2.*pi / L
        for i in range(10): # Adjust distribution using Newton iterations
            pos -= ( pos + alpha*sin(k*pos)/k - pos0 ) / ( 1. + alpha*cos(k*pos) )
            
        # Normal velocity distribution
        vel = random.normal(0.0, 1.0, npart)
        
        return pos, vel
    
    def twostream(npart, L, vbeam=2):
        # Start with a uniform distribution of positions
        pos = random.uniform(0., L, npart)
        # Normal velocity distribution
        vel = random.normal(0.0, 1.0, npart)
        
        np2 = int(npart / 2)
        vel[:np2] += vbeam  # Half the particles moving one way
        vel[np2:] -= vbeam  # and half the other
        
        return pos,vel
    
    ####################################################################
    
    if __name__ == "__main__":
        # Generate initial condition
        #
        if False:
            # 2-stream instability
            L = 100
            ncells = 20
            pos, vel = twostream(10000, L, 3.)
        else:
            # Landau damping
            L = 4.*pi
            ncells = 20
            npart = 20000
            pos, vel = landau(npart, L)
        
        # Create some output classes
        #p = Plot(pos, vel, ncells, L) # This displays an animated figure
        s = Summary()                 # Calculates, stores and prints summary info
       
    	
        # Run the simulation
        pos, vel = run(pos, vel, L, ncells, 
                       out=[s],    #Animation off
    				   #out=[p, s],            #Animation On          # These are called each output
                       output_times=linspace(0.,40,100)) # The times to output
      
        # Summary stores an array of the first-harmonic amplitude
        # Make a semilog plot to see exponential damping
    	
        x = s.t
        data = s.firstharmonic
    	
    	#Local minimum and local maximum
        b = (diff(sign(diff(data))) > 0).nonzero()[0] + 1 # local min
        c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max
    	
        timemin = []
        minh = []
    	
        for j in b:
    	    timemin.append(x[j])
    	    minh.append(data[j])
    		
        derivListmin = [] #Finding the change in gradients between the minimum peaks of the Harmonics
    	
        for j in range(len(timemin)-1):
            derivmin = (minh[j]-minh[j+1])/(timemin[j]-timemin[j+1])
            derivListmin.append(derivmin)
    		
        print derivListmin
    	
        for indexmin, valuemin in enumerate(derivListmin):
    	    if valuemin > 0: #j refers to each element in derivlist
    		   print indexmin, valuemin
    		   variablemin = indexmin 
    		   break
    		   
        print "timemin =", timemin[variablemin]
        print "minh =", minh[variablemin]
    		   
        higherindexminh = []
        for j in range(len(minh)):
    		if j > variablemin:
    		   higherindexminh.append(minh[j])
    	
        print "higherindexminh = ", higherindexminh
    	
        noiselevelmin = mean(higherindexminh)
    	
        print noiselevelmin
    	
    	
        timemax = []
        maxh = []
        for i in c:
            timemax.append(x[i])
            maxh.append(data[i])
    		
    	
        derivListmax = [] #Finding the change in gradients between the maximum peaks of the Harmonics
    	
        for i in range(len(timemax)-1):
            derivmax = (maxh[i]-maxh[i+1])/(timemax[i]-timemax[i+1])
            derivListmax.append(derivmax) #This is the gradient
    		
        print derivListmax
    	
    	#Index is the position of the gradient in the list of derivList
        for indexmax, valuemax in enumerate(derivListmax):
    	    if valuemax > 0: #For when the derivative between the peaks first becomes postive, print all postive values after it
    		   print indexmax, valuemax
    		   variablemax = indexmax
    		   break
    		   
        print "timemax =", timemax[variablemax]
        print "maxh =", maxh[variablemax]
    		   
        higherindexmaxh = [] #creating a new list that includes every peak that follows when the derivative first becomes positive
    	
        for i in range(len(maxh)):
    		if i > variablemax: #An index/position of higher value than when the derivative first becomes positive
    		   higherindexmaxh.append(maxh[i]) #Fill the list with such values
    	
        print "higherindexmaxh = ", higherindexmaxh
    	
        noiselevelmax = mean(higherindexmaxh) #Average all the values in the list to imitate a constant noise.
    	
        print noiselevelmax
    	
        noiselevelavg = (noiselevelmax + noiselevelmin)/2
    	
    	
        plt.figure()
        plt.plot(s.t, s.firstharmonic)
        plt.xlabel("Time [Normalised]")
        plt.ylabel("First harmonic amplitude [Normalised]")
        plt.yscale('log')
        plt.plot(timemax, maxh, 'or')
        plt.plot(timemin, minh, 'oy')
        plt.plot([0,40],[noiselevelavg, noiselevelavg])
        #plt.plot([0,40],[noiselevelmin, noiselevelmin])
        
    	
    	#show()
    
        plt.ioff() # This so that the windows stay open
        plt.show()
    All the parts highlighted in Red are what I struggled to write with my limited knowledge.

    It works and gives me an output, but is there a way to make this code run faster and more efficiently? Its a computational project and we get extra marks for making the WHOLE code (even the parts I didn't write) more efficient.

    Apparently a sure fire way is to rewrite the whole thing in Fortran, but that would be impossible for me so Python is ok.

    What this code does, it that it generates a distribution of particles in a plasma and applies a randomized damping (Landau Damping) in order to imitate experimental conditions. I can change the number of particles, the magnitude of the damping and so forth.

    My code at the end crudely finds the noise level of the damping (where noise is indistinguishable from Landau Damping) by using derivatives of the maximal and minimal values of Landau Harmonics.

    Please help me guys
  2. #2
  3. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    Please, for whom are you working?
    [code]Code tags[/code] are essential for python code and Makefiles!
  4. #3
  5. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0
    I don't work for anyone..

    I'm doing a Postgrad in Plasma Physics & Fusion.. The code simply imitates conditions in an experimental Tokamak.

    You can run it in your Python GUI to see how it works.

    Is there a way I can make the code run faster or more efficiently.
  6. #4
  7. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    I removed the final statement of program

    #plt.show() # removed by dwl for profiling.

    PROFILING COMMAND

    Then profiled the code
    Code:
    $ ( cd /tmp && python -m cProfile p.py > /tmp/pic.profile.txt & ) &
    Then sorted the interesting lines of pic.profile.txt (piping an emacs region to sort)
    SORT THE PROFILE TO MAKE SENSE OF IT
    Code:
    sort --numeric-sort --key=2,2
    or
    $ sort --numeric-sort --key=2,2 /tmp/pic.profile.txt
    the tail of which, with titles, looks like
    Code:
             200749 function calls (198054 primitive calls) in 231.770 seconds
    
       Ordered by: standard name (reordered by sort)
    
       ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    
    ...
    
            1    0.141    0.141  231.045  231.045 p.py:131(run)
         4805    0.185    0.000    0.185    0.000 {numpy.core.multiarray.concatenate}
          595    0.987    0.002  219.888    0.370 p.py:22(rk4step)
          813    2.186    0.003    2.186    0.003 {max}
         2380    2.753    0.001  218.901    0.092 p.py:102(pic)
         2380    4.352    0.002    4.481    0.002 p.py:58(periodic_interp)
         2480  219.846    0.089  219.862    0.089 p.py:31(calc_density)
    showing that calc_density completely dominates the computation time. Therefor, without actually having done it, I advise replace this loop
    Code:
        for p in position / dx:    # Loop over all the particles, converting position into a cell number
            plower = int(p)        # Cell to the left (rounding down)
            offset = p - plower    # Offset from the left
            density[plower] += 1. - offset
            density[(plower + 1) % ncells] += offset
    with numpy array operations.

    It looks like you distributed charge to nearest nodes using position consideration only. This will cause a net force on the particle, yet particles should not have a "self-force" driving them in some direction.
    Last edited by b49P23TIvg; January 29th, 2013 at 12:48 PM.
    [code]Code tags[/code] are essential for python code and Makefiles!
  8. #5
  9. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0
    Can you explain how you profiled the code like that? It looks really interesting - it would be great to replicate that!

    I will get right on the numpy array - expect me to ask for help again though.

    I'm not too sure about the "net force".. I will study the code again.

    I'm assuming I will need to manipulate this:

    astype(int) ???
  10. #6
  11. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    This command profiles the code
    Code:
    $ ( cd /tmp && python -m cProfile p.py > /tmp/pic.profile.txt & ) &
    You are using a unix. You don't need all the parenthesis and background. I was using my shell for other purpose.

    Convince yourself that distributing charge based on distance-weighted nearest two grid points causes self-force. (Other distributions can have other problems, of course.) If I recall the electrostatics correctly,
    Code:
    0 . . q 1
    
    Node 0 at x==0, node 1 at x==1
    Charge at 3/4
    If you distribute the charge as
    3/4 at node 1 and 1/4 at grid 0 then the field
    at 3/4 is proportional to
    
       (-&((%*~)/) |.)1r4 3r4  NB. www.jsoftware.com
    _104r9
       
    Magnitudewise, that's bigger than 10.
    I'll send you an article on probepic if you send your email address. The document is too large for me to attach here.
    Last edited by b49P23TIvg; January 29th, 2013 at 01:18 PM.
    [code]Code tags[/code] are essential for python code and Makefiles!
  12. #7
  13. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0
    I'll PM you my email thanks

    Here is what I've done so far:

    Old code:

    Code:
    for p in position / dx:    # Loop over all the particles, converting position into a cell number
            plower = int(p) #Changing values in the array to integers
            #plower = int(p)        # Cell to the left (rounding down)
            offset = p - plower    # Offset from the left
            density[plower] += 1. - offset
            density[(plower + 1) % ncells] += offset
    New code:

    Code:
    density, _ = histogram(position, ncells, (0, L))
    I figured I should keep it simple since I'm not very good. The program runs (at an exceptionally fast speed I might add) but it gives silly values for the density.

    I think it is binning the particles in a random fashion, so all my plots and fittings don't look good.

    A comparison between the old density and new histogram density shows that the new density is slightly out of place from the old one.

    I think it is to do with not including this section of code:

    Code:
    density[plower] += 1. - offset
            density[(plower + 1) % ncells] += offset
    This part of the code makes sure that particles not in integer positions are placed into their nearest cells. But it looks extremely tricky (as in I have no idea) trying to include it into a histogram version.

    Maybe the histogram is too crude?

    I have updated my code since I last replied with some fittings, take a look and run it if you want (if I didn't explain things properly):

    Code:
    #!/usr/bin/env python
    #
    # Electrostatic PIC code in a 1D cyclic domain
    
    from numpy import arange, concatenate, zeros, linspace, floor, array, pi
    from numpy import diff, sign, sin, cos, sqrt, random, histogram
    from numpy import *
    
    from numpy import histogram
    
    from scipy.optimize import curve_fit
    
    import matplotlib.pyplot as plt # Matplotlib plotting library
    
    try:
        import matplotlib.gridspec as gridspec  # For plot layout grid
        got_gridspec = True
    except:
        got_gridspec = False
    
    # Need an FFT routine, either from SciPy or NumPy
    try:
        from scipy.fftpack import fft, ifft
    except:
        # No SciPy FFT routine. Import NumPy routine instead
        from numpy.fft import fft, ifft
    
    def rk4step(f, y0, dt, args=()):
        """ Takes a single step using RK4 method """
        k1 = f(y0, *args)
        k2 = f(y0 + 0.5*dt*k1, *args)
        k3 = f(y0 + 0.5*dt*k2, *args)
        k4 = f(y0 + dt*k3, *args)
    
        return y0 + (k1 + 2.*k2 + 2.*k3 + k4)*dt / 6.
    
    def calc_density(position, ncells, L):
        """ Calculate charge density given particle positions
        
        Input
          position  - Array of positions, one for each particle
                      assumed to be between 0 and L
          ncells    - Number of cells
          L         - Length of the domain
    
        Output
          density   - contains 1 if evenly distributed
        """
        # This is a crude method and could be made more efficient
        
        density = zeros([ncells])
        nparticles = len(position)
    	
        '''Can replace the following with Numpy array operations'''
    
        dx = L / ncells       # Uniform cell spacing, L and ncells are values
    	
        '''p = position / dx
        plower = p.astype(int) #always rounds down
        offset = p - plower    # Offset from the left
    	
        density2 = zeros([ncells])
    	
        density2, _ = histogram(position, ncells, (0, L), normed = True, weights )
        density2[plower] += 1. - offset
        density2[(plower + 1) % ncells] += offset
    	'''
    	
    	
        
        for p in position / dx:    # Loop over all the particles, converting position into a cell number
            plower = int(p) #Changing values in the array to integers
            #plower = int(p)        # Cell to the left (rounding down)
            offset = p - plower    # Offset from the left
            density[plower] += 1. - offset
            density[(plower + 1) % ncells] += offset
    	
        '''plt.figure()
        plt.plot(density)
        plt.plot(density2)
        plt.ioff()
        plt.show()
        '''
    	
        # nparticles now distributed amongst ncells
        density *= float(ncells) / float(nparticles)  # Make average density equal to 1
        return density
    
    def periodic_interp(y, x):
        """
        Linear interpolation of a periodic array y at index x
        
        Input
    
        y - Array of values to be interpolated
        x - Index where result required. Can be an array of values
        
        Output
        
        y[x] with non-integer x
        """
        ny = len(y)
        if len(x) > 1:
            y = array(y) # Make sure it's a NumPy array for array indexing
        xl = floor(x).astype(int) # Left index
        dx = x - xl
        xl = ((xl % ny) + ny) % ny  # Ensures between 0 and ny-1 inclusive
        return y[xl]*(1. - dx) + y[(xl+1)%ny]*dx
    
    def fft_integrate(y):
        """ Integrate a periodic function using FFTs
        """
        n = len(y) # Get the length of y
        
        f = fft(y) # Take FFT
        # Result is in standard layout with positive frequencies first then negative
        # n even: [ f(0), f(1), ... f(n/2), f(1-n/2) ... f(-1) ]
        # n odd:  [ f(0), f(1), ... f((n-1)/2), f(-(n-1)/2) ... f(-1) ]
        
        if n % 2 == 0: # If an even number of points
            k = concatenate( (arange(0, n/2+1), arange(1-n/2, 0)) )
        else:
            k = concatenate( (arange(0, (n-1)/2+1), arange( -(n-1)/2, 0)) )
        k = 2.*pi*k/n
        
        # Modify frequencies by dividing by ik
        f[1:] /= (1j * k[1:]) 
        f[0] = 0. # Set the arbitrary zero-frequency term to zero
        
        return ifft(f).real # Reverse Fourier Transform
       
    
    def pic(f, ncells, L):
        """ f contains the position and velocity of all particles
        """
        nparticles = len(f)/2     # Two values for each particle
        pos = f[0:nparticles] # Position of each particle
        vel = f[nparticles:]      # Velocity of each particle
    
        dx = L / float(ncells)    # Cell spacing
    
        # Ensure that pos is between 0 and L
        pos = ((pos % L) + L) % L
        
        # Calculate number density, normalised so 1 when uniform
        density = calc_density(pos, ncells, L)
        
        # Subtract ion density to get total charge density
        rho = density - 1.
        
        # Calculate electric field
        E = -fft_integrate(rho)*dx
        
        # Interpolate E field at particle locations
        accel = -periodic_interp(E, pos/dx)
    
        # Put back into a single array
        return concatenate( (vel, accel) )
    
    ####################################################################
    
    def run(pos, vel, L, ncells=None, out=[], output_times=linspace(0,20,100), cfl=0.5):
        
        if ncells == None:
            ncells = int(sqrt(len(pos))) # A sensible default
    
        dx = L / float(ncells)
        
        f = concatenate( (pos, vel) )   # Starting state
        nparticles = len(pos)
        
        time = 0.0
        for tnext in output_times:
            # Advance to tnext
            stepping = True
            while stepping:
                # Maximum distance a particle can move is one cell
                dt = cfl * dx / max(abs(vel))
                if time + dt >= tnext:
                    # Next time will hit or exceed required output time
                    stepping = False
                    dt = tnext - time
                #print "Time: ", time, dt
                f = rk4step(pic, f, dt, args=(ncells, L))
                time += dt
                
            # Extract position and velocities
            pos = ((f[0:nparticles] % L) + L) % L
            vel = f[nparticles:]
            
            # Send to output functions
            for func in out:
                func(pos, vel, ncells, L, time)
            
        return pos, vel
    
    ####################################################################
    # 
    # Output functions and classes
    #
    
    class Plot:
        """
        Displays three plots: phase space, charge density, and velocity distribution
        """
        def __init__(self, pos, vel, ncells, L):
            
            d = calc_density(pos, ncells, L)
            vhist, bins  = histogram(vel, int(sqrt(len(vel))))
            vbins = 0.5*(bins[1:]+bins[:-1])
            
            # Plot initial positions
            if got_gridspec:
                self.fig = plt.figure()
                self.gs = gridspec.GridSpec(4, 4)
                ax = self.fig.add_subplot(self.gs[0:3,0:3])
                self.phase_plot = ax.plot(pos, vel, '.')[0]
                ax.set_title("Phase space")
                
                ax = self.fig.add_subplot(self.gs[3,0:3])
                self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]
                
                ax = self.fig.add_subplot(self.gs[0:3,3])
                self.vel_plot = ax.plot(vhist, vbins)[0]
            else:
                self.fig = plt.figure()
                self.phase_plot = plt.plot(pos, vel, '.')[0]
                
                self.fig = plt.figure()
                self.density_plot = plt.plot(linspace(0, L, ncells), d)[0]
                
                self.fig = plt.figure()
                self.vel_plot = plt.plot(vhist, vbins)[0]
            plt.ion()
            plt.show()
            
        def __call__(self, pos, vel, ncells, L, t):
            d = calc_density(pos, ncells, L)
            vhist, bins  = histogram(vel, int(sqrt(len(vel))))
            vbins = 0.5*(bins[1:]+bins[:-1])
            
            self.phase_plot.set_data(pos, vel) # Update the plot
            self.density_plot.set_data(linspace(0, L, ncells), d)
            self.vel_plot.set_data(vhist, vbins)
            plt.draw()
            
    
    class Summary:
        def __init__(self):
            self.t = []
            self.firstharmonic = []
            
        def __call__(self, pos, vel, ncells, L, t):
            # Calculate the charge density
            d = calc_density(pos, ncells, L)
            
            # Amplitude of the first harmonic
            fh = 2.*abs(fft(d)[1]) / float(ncells)
            
            print "Time:", t, "First:", fh
            
            self.t.append(t)
            self.firstharmonic.append(fh)
    
    ####################################################################
    # 
    # Functions to create the initial conditions
    #
    
    def landau(npart, L, alpha=0.2):
        """
        Creates the initial conditions for Landau damping
        
        """
        # Start with a uniform distribution of positions
        pos = random.uniform(0., L, npart)
        pos0 = pos.copy()
        k = 2.*pi / L
        for i in range(10): # Adjust distribution using Newton iterations
            pos -= ( pos + alpha*sin(k*pos)/k - pos0 ) / ( 1. + alpha*cos(k*pos) )
            
        # Normal velocity distribution
        vel = random.normal(0.0, 1.0, npart)
        
        return pos, vel
    
    def twostream(npart, L, vbeam=2):
        # Start with a uniform distribution of positions
        pos = random.uniform(0., L, npart)
        # Normal velocity distribution
        vel = random.normal(0.0, 1.0, npart)
        
        np2 = int(npart / 2)
        vel[:np2] += vbeam  # Half the particles moving one way
        vel[np2:] -= vbeam  # and half the other
        
        return pos,vel
    
    ####################################################################
    
    if __name__ == "__main__":
        # Generate initial condition
        #
        if False:
            # 2-stream instability
            L = 100
            ncells = 20
            pos, vel = twostream(10000, L, 3.)
        else:
            # Landau damping
            L = 4.*pi
            ncells = 20
            npart = input("How many Particles? ") #Requires the users input in order to aid repeat runs.
            pos, vel = landau(npart, L)
    		
    		
        
        # Create some output classes
        #p = Plot(pos, vel, ncells, L) # This displays an animated figure
        s = Summary()                 # Calculates, stores and prints summary info
       
    	
        # Run the simulation
        pos, vel = run(pos, vel, L, ncells, 
                       out=[s],    #Animation off
    				   #out=[p, s],            #Animation On          # These are called each output
                       output_times=linspace(0.,40,100)) # The times to output
      
        # Summary stores an array of the first-harmonic amplitude
        # Make a semilog plot to see exponential damping
    	
        x = s.t
        data = s.firstharmonic
    	
    	#Local minimum and local maximum
        b = (diff(sign(diff(data))) > 0).nonzero()[0] + 1 # local min
        c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max
    	
        timemin = []
        minh = []
    	
        for j in b:
    	    timemin.append(x[j])
    	    minh.append(data[j])
    		
        derivListmin = [] #Finding the change in gradients between the minimum peaks of the Harmonics
    	
        for j in range(len(timemin)-1):
            derivmin = (minh[j]-minh[j+1])/(timemin[j]-timemin[j+1])
            derivListmin.append(derivmin)
    		
        #print derivListmin
    	
        for indexmin, valuemin in enumerate(derivListmin):
    	    if valuemin > 0: #j refers to each element in derivlist
    		   #print indexmin, valuemin
    		   variablemin = indexmin 
    		   break
    		   
        #print "timemin =", timemin[variablemin]
        #print "minh =", minh[variablemin]
    		   
        higherindexminh = []
        for j in range(len(minh)):
    		if j > variablemin:
    		   higherindexminh.append(minh[j])
    	
        #print "higherindexminh = ", higherindexminh
    	
        noiselevelmin = mean(higherindexminh)
    	
        #print noiselevelmin
    	
    	
        timemax = []
        maxh = []
        for i in c:
            timemax.append(x[i])
            maxh.append(data[i])
    		
    	
        derivListmax = [] #Finding the change in gradients between the maximum peaks of the Harmonics
    	
        for i in range(len(timemax)-1):
            derivmax = (maxh[i]-maxh[i+1])/(timemax[i]-timemax[i+1])
            derivListmax.append(derivmax) #This is the gradient
    		
        #print derivListmax
    	
    	#Index is the position of the gradient in the list of derivList
        for indexmax, valuemax in enumerate(derivListmax):
    	    if valuemax > 0: #For when the derivative between the peaks first becomes postive, print all postive values after it
    		   #print indexmax, valuemax
    		   variablemax = indexmax
    		   break
    		   
        #print "timemax =", timemax[variablemax]
        #print "maxh =", maxh[variablemax]
    		   
        higherindexmaxh = [] #creating a new list that includes every peak that follows when the derivative first becomes positive
    	
        for i in range(len(maxh)):
    		if i > variablemax: #An index/position of higher value than when the derivative first becomes positive
    		   higherindexmaxh.append(maxh[i]) #Fill the list with such values
    	
        #print "higherindexmaxh = ", higherindexmaxh
    	
        noiselevelmax = mean(higherindexmaxh) #Average all the values in the list to imitate a constant noise.
    	
        #print noiselevelmax
    	
        noiselevelavg = (noiselevelmax + noiselevelmin)/2
    
        file = open('c:/Users/Pinky/Programming/PIC.txt','a')
        file.write((str(noiselevelavg) + '   ' + str(npart)) + "\n")
        file.close()
    	
    	# file = open('c:/Users/Pinky/Programming/PIC.txt','r')
    	#file.read()
    	#file.close()
    	
        noiselist = []
        particlelist = []
    	
        for line in open('c:/Users/Pinky/Programming/PIC.txt', 'r'):
            words = line.split() #This will split the rows in the file into 2 list by spaces
            noiselist.append(words[0])
            particlelist.append(words[1])
    		
        file.close()
    		
        print noiselist
    	
        '''
    	Investigating a way of calculating the accuracy of the noise level in conjuction
    	with finding the period of the Laundau damping.  Will superimpose a sinusoidal
    	function which will exponentially damp itself over time.  Subtracting this function away
    	from the original Laundau function will give a residue - take this as noise
        '''
        def func(x, b, k, a):
    	    return (b*abs(cos(k*x))*exp(-x*a))
    	
        x = linspace(0, 40, 100)
        popt, pcov = curve_fit(func, x, data, [0.2, 1.4, 0.1])
        print "params = ", popt
        fitdata = func(x, popt[0], popt[1], popt[2])
    	
        residue = ((s.firstharmonic - fitdata))
    	
        residueavg = mean(residue)
    	
        plt.plot(particlelist, noiselist, 'or')
        plt.xlabel("Noise")
        plt.ylabel("Number of Particles")
    	
        plt.figure()
        plt.plot(s.t, s.firstharmonic)
        plt.plot([0,40],[residueavg,residueavg],"-o")
        plt.xlabel("Time [Normalised]")
        plt.ylabel("First harmonic amplitude [Normalised]")
        plt.yscale('log')
        plt.plot(timemax, maxh, 'or')
        plt.plot(timemin, minh, 'oy')
        plt.plot([0,40],[noiselevelavg, noiselevelavg])
        plt.plot(x, fitdata)
    	
    	#plt.plot([0,40],[noiselevelmax, noiselevelmax])
        #plt.plot([0,40],[noiselevelmin, noiselevelmin])
        
    	
    	#show()
    
        plt.ioff() # This so that the windows stay open
        plt.show()
    The program will prompt you for the number of particles, 5000 should do it.

    Any pointers on how to tackle this problem?

    Many thanks for your help so far, it is much appreciated.
  14. #8
  15. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    I haven't run the code with histogram.

    You could use a more sophisticated ODE solver,
    scipy interface to odepack to get possibly lower error with fewer function evaluations.---
    never mind. Your time step is controlled by 1 cell traversal by fastest particle.
    Last edited by b49P23TIvg; January 29th, 2013 at 09:30 PM.
    [code]Code tags[/code] are essential for python code and Makefiles!
  16. #9
  17. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    Rewritten calc_density function improves run time to 1/5 that of post 1. Do the results agree? I think so! I attempted to retain your position-weighting particle distribution. I reduced the number of computations that take place in that loop over particles by moving the array computations outside the loop. Replacing that loop with c code would strongly drive down the run time.
    Code:
    # Electrostatic PIC code in a 1D cyclic domain
    
    import numpy
    from numpy import arange, concatenate, zeros, linspace, floor, array, pi
    from numpy import diff, sign, sin, cos, sqrt, random, histogram
    from numpy import *
    
    import matplotlib.pyplot as plt # Matplotlib plotting library
    
    try:
        import matplotlib.gridspec as gridspec  # For plot layout grid
        got_gridspec = True
    except:
        got_gridspec = False
    
    # Need an FFT routine, either from SciPy or NumPy
    try:
        from scipy.fftpack import fft, ifft,integrate
    except:
        # No SciPy FFT routine. Import NumPy routine instead
        from numpy.fft import fft, ifft
    
    def rk4step(f, y0, dt, args=()):
        """ Takes a single step using RK4 method """
        k1 = f(y0, *args)
        k2 = f(y0 + 0.5*dt*k1, *args)
        k3 = f(y0 + 0.5*dt*k2, *args)
        k4 = f(y0 + dt*k3, *args)
        return y0 + (k1 + 2.*k2 + 2.*k3 + k4)*dt / 6.
    
    def calc_density(position, ncells, L):
        """ Calculate charge density given particle positions
    
        Input
          position  - Array of positions, one for each particle
                      assumed to be between 0 and L
          ncells    - Number of cells
          L         - Length of the domain
    
        Output
          density   - contains 1 if evenly distributed
        """
        # This is a crude method and could be made more efficient
    
        #density = zeros([ncells])
        nparticles = len(position)
    
        dx = L / ncells       # Uniform cell spacing
        #for p in position / dx:    # Loop over all the particles, converting position into a cell number
        #    plower = int(p)        # Cell to the left (rounding down)
        #    offset = p - plower    # Offset from the left
        #    density[plower] += 1. - offset
        #    density[(plower + 1) % ncells] += offset
        p = position/dx
        plower = numpy.array(p,dtype=numpy.integer)
        offset = p - plower
        density = zeros([ncells+1])
        ld = 1-offset
        for (low,ld,ud,) in zip(plower,ld,offset,):
            density[low] += ld
            density[low+1] += ud
        density[0] += density[-1]
        density = density[:-1]
        # nparticles now distributed amongst ncells
        density *= float(ncells) / float(nparticles)  # Make average density equal to 1
        return density
    
    def periodic_interp(y, x):
        """
        Linear interpolation of a periodic array y at index x
    
        Input
    
        y - Array of values to be interpolated
        x - Index where result required. Can be an array of values
    
        Output
    
        y[x] with non-integer x
        """
        ny = len(y)
        if len(x) > 1:
            y = array(y) # Make sure it's a NumPy array for array indexing
        xl = floor(x).astype(int) # Left index
        dx = x - xl
        xl = ((xl % ny) + ny) % ny  # Ensures between 0 and ny-1 inclusive
        return y[xl]*(1. - dx) + y[(xl+1)%ny]*dx
    
    def fft_integrate(y):
        """ Integrate a periodic function using FFTs
        """
        n = len(y) # Get the length of y
    
        f = fft(y) # Take FFT
        # Result is in standard layout with positive frequencies first then negative
        # n even: [ f(0), f(1), ... f(n/2), f(1-n/2) ... f(-1) ]
        # n odd:  [ f(0), f(1), ... f((n-1)/2), f(-(n-1)/2) ... f(-1) ]
    
        if n % 2 == 0: # If an even number of points
            k = concatenate( (arange(0, n/2+1), arange(1-n/2, 0)) )
        else:
            k = concatenate( (arange(0, (n-1)/2+1), arange( -(n-1)/2, 0)) )
        k = 2.*pi*k/n
    
        # Modify frequencies by dividing by ik
        f[1:] /= (1j * k[1:])
        f[0] = 0. # Set the arbitrary zero-frequency term to zero
    
        return ifft(f).real # Reverse Fourier Transform
    
    
    def pic(f, ncells, L):
        """ f contains the position and velocity of all particles
        """
        nparticles = len(f)/2     # Two values for each particle
        pos = f[0:nparticles] # Position of each particle
        vel = f[nparticles:]      # Velocity of each particle
    
        dx = L / float(ncells)    # Cell spacing
    
        # Ensure that pos is between 0 and L
        pos = ((pos % L) + L) % L
    
        # Calculate number density, normalised so 1 when uniform
        density = calc_density(pos, ncells, L)
    
        # Subtract ion density to get total charge density
        rho = density - 1.
    
        # Calculate electric field
        E = -fft_integrate(rho)*dx
    
        # Interpolate E field at particle locations
        accel = -periodic_interp(E, pos/dx)
    
        # Put back into a single array
        return concatenate( (vel, accel) )
    
    ####################################################################
    
    def run(pos, vel, L, ncells=None, out=[], output_times=linspace(0,20,100), cfl=0.5):
    
        if ncells == None:
            ncells = int(sqrt(len(pos))) # A sensible default
    
        dx = L / float(ncells)
    
        f = concatenate( (pos, vel) )   # Starting state
        nparticles = len(pos)
    
        time = 0.0
        for tnext in output_times:
            # Advance to tnext
            stepping = True
            while stepping:
                # Maximum distance a particle can move is one cell
                dt = cfl * dx / max(abs(vel))
                if time + dt >= tnext:
                    # Next time will hit or exceed required output time
                    stepping = False
                    dt = tnext - time
                #print "Time: ", time, dt
                #scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)[source]
                f = rk4step(pic, f, dt, args=(ncells, L))
                time += dt
    
            # Extract position and velocities
            pos = ((f[0:nparticles] % L) + L) % L
            vel = f[nparticles:]
    
            # Send to output functions
            for func in out:
                func(pos, vel, ncells, L, time)
    
        return pos, vel
    
    ####################################################################
    #
    # Output functions and classes
    #
    
    class Plot:
        """
        Displays three plots: phase space, charge density, and velocity distribution
        """
        def __init__(self, pos, vel, ncells, L):
    
            d = calc_density(pos, ncells, L)
            vhist, bins  = histogram(vel, int(sqrt(len(vel))))
            vbins = 0.5*(bins[1:]+bins[:-1])
    
            # Plot initial positions
            if got_gridspec:
                self.fig = plt.figure()
                self.gs = gridspec.GridSpec(4, 4)
                ax = self.fig.add_subplot(self.gs[0:3,0:3])
                self.phase_plot = ax.plot(pos, vel, '.')[0]
                ax.set_title("Phase space")
    
                ax = self.fig.add_subplot(self.gs[3,0:3])
                self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]
    
                ax = self.fig.add_subplot(self.gs[0:3,3])
                self.vel_plot = ax.plot(vhist, vbins)[0]
            else:
                self.fig = plt.figure()
                self.phase_plot = plt.plot(pos, vel, '.')[0]
    
                self.fig = plt.figure()
                self.density_plot = plt.plot(linspace(0, L, ncells), d)[0]
    
                self.fig = plt.figure()
                self.vel_plot = plt.plot(vhist, vbins)[0]
            plt.ion()
            plt.show()
    
        def __call__(self, pos, vel, ncells, L, t):
            d = calc_density(pos, ncells, L)
            vhist, bins  = histogram(vel, int(sqrt(len(vel))))
            vbins = 0.5*(bins[1:]+bins[:-1])
    
            self.phase_plot.set_data(pos, vel) # Update the plot
            self.density_plot.set_data(linspace(0, L, ncells), d)
            self.vel_plot.set_data(vhist, vbins)
            plt.draw()
    
    
    class Summary:
        def __init__(self):
            self.t = []
            self.firstharmonic = []
    
        def __call__(self, pos, vel, ncells, L, t):
            # Calculate the charge density
            d = calc_density(pos, ncells, L)
    
            # Amplitude of the first harmonic
            fh = 2.*abs(fft(d)[1]) / float(ncells)
    
            print "Time:", t, "First:", fh
    
            self.t.append(t)
            self.firstharmonic.append(fh)
    
    ####################################################################
    #
    # Functions to create the initial conditions
    #
    
    def landau(npart, L, alpha=0.2):
        """
        Creates the initial conditions for Landau damping
    
        """
        # Start with a uniform distribution of positions
        pos = random.uniform(0., L, npart)
        pos0 = pos.copy()
        k = 2.*pi / L
        for i in range(10): # Adjust distribution using Newton iterations
            pos -= ( pos + alpha*sin(k*pos)/k - pos0 ) / ( 1. + alpha*cos(k*pos) )
    
        # Normal velocity distribution
        vel = random.normal(0.0, 1.0, npart)
    
        return pos, vel
    
    def twostream(npart, L, vbeam=2):
        # Start with a uniform distribution of positions
        pos = random.uniform(0., L, npart)
        # Normal velocity distribution
        vel = random.normal(0.0, 1.0, npart)
    
        np2 = int(npart / 2)
        vel[:np2] += vbeam  # Half the particles moving one way
        vel[np2:] -= vbeam  # and half the other
    
        return pos,vel
    
    ####################################################################
    
    if __name__ == "__main__":
        # Generate initial condition
        #
        if True:
            # 2-stream instability
            L = 100
            ncells = 20
            pos, vel = twostream(10000, L, 3.)
        else:
            # Landau damping
            L = 4.*pi
            ncells = 20
            npart = 20000
            pos, vel = landau(npart, L)
    
        # Create some output classes
        #p = Plot(pos, vel, ncells, L) # This displays an animated figure
        s = Summary()                 # Calculates, stores and prints summary info
    
    
        # Run the simulation
        pos, vel = run(pos, vel, L, ncells,
                       out=[s],    #Animation off
    				   #out=[p, s],            #Animation On          # These are called each output
                       output_times=linspace(0.,40,100)) # The times to output
    
        # Summary stores an array of the first-harmonic amplitude
        # Make a semilog plot to see exponential damping
    
        x = s.t
        data = s.firstharmonic
    
    	#Local minimum and local maximum
        b = (diff(sign(diff(data))) > 0).nonzero()[0] + 1 # local min
        c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max
    
        timemin = []
        minh = []
    
        for j in b:
    	    timemin.append(x[j])
    	    minh.append(data[j])
    
        derivListmin = [] #Finding the change in gradients between the minimum peaks of the Harmonics
    
        for j in range(len(timemin)-1):
            derivmin = (minh[j]-minh[j+1])/(timemin[j]-timemin[j+1])
            derivListmin.append(derivmin)
    
        print derivListmin
    
        for indexmin, valuemin in enumerate(derivListmin):
    	    if valuemin > 0: #j refers to each element in derivlist
    		   print indexmin, valuemin
    		   variablemin = indexmin
    		   break
    
        print "timemin =", timemin[variablemin]
        print "minh =", minh[variablemin]
    
        higherindexminh = []
        for j in range(len(minh)):
    		if j > variablemin:
    		   higherindexminh.append(minh[j])
    
        print "higherindexminh = ", higherindexminh
    
        noiselevelmin = mean(higherindexminh)
    
        print noiselevelmin
    
    
        timemax = []
        maxh = []
        for i in c:
            timemax.append(x[i])
            maxh.append(data[i])
    
    
        derivListmax = [] #Finding the change in gradients between the maximum peaks of the Harmonics
    
        for i in range(len(timemax)-1):
            derivmax = (maxh[i]-maxh[i+1])/(timemax[i]-timemax[i+1])
            derivListmax.append(derivmax) #This is the gradient
    
        print derivListmax
    
    	#Index is the position of the gradient in the list of derivList
        for indexmax, valuemax in enumerate(derivListmax):
    	    if valuemax > 0: #For when the derivative between the peaks first becomes postive, print all postive values after it
    		   print indexmax, valuemax
    		   variablemax = indexmax
    		   break
    
        print "timemax =", timemax[variablemax]
        print "maxh =", maxh[variablemax]
    
        higherindexmaxh = [] #creating a new list that includes every peak that follows when the derivative first becomes positive
    
        for i in range(len(maxh)):
    		if i > variablemax: #An index/position of higher value than when the derivative first becomes positive
    		   higherindexmaxh.append(maxh[i]) #Fill the list with such values
    
        print "higherindexmaxh = ", higherindexmaxh
    
        noiselevelmax = mean(higherindexmaxh) #Average all the values in the list to imitate a constant noise.
    
        print noiselevelmax
    
        noiselevelavg = (noiselevelmax + noiselevelmin)/2
    
    
        plt.figure()
        plt.plot(s.t, s.firstharmonic)
        plt.xlabel("Time [Normalised]")
        plt.ylabel("First harmonic amplitude [Normalised]")
        plt.yscale('log')
        plt.plot(timemax, maxh, 'or')
        plt.plot(timemin, minh, 'oy')
        plt.plot([0,40],[noiselevelavg, noiselevelavg])
        #plt.plot([0,40],[noiselevelmin, noiselevelmin])
    
    
    	#show()
    
        plt.ioff() # This so that the windows stay open
        plt.show() # removed by dwl for profiling.


    Code:
             151233 function calls (148538 primitive calls) in 46.396 seconds
    
       Ordered by: reordered by tottime
    
       ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    
    ...
          417    0.348    0.001    0.348    0.001 {max}
          796    0.461    0.001   40.268    0.051 p.py:112(pic)
          796    0.761    0.001    0.785    0.001 p.py:68(periodic_interp)
         1085    3.717    0.003    3.717    0.003 {zip}
          896   39.986    0.045   43.725    0.049 p.py:31(calc_density)
    [code]Code tags[/code] are essential for python code and Makefiles!
  18. #10
  19. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481

    Brief introduction to PIC method.


    To compute force on a charged particle due to all particles (10 thousand, perhaps) is a daunting task. The particle in cell method distributes the charges onto a coarse grid (20 points in this linear model) and then computes the force on each particle from the grid, making the computation tractable. Using sum of external forces = mass * acceleration the code can integrate to update the particle positions.

    Other notes: each "particle" in the model would represent many particles. In this plasma simulation the ions having mass millions of times heavier than the electrons are fixed and evenly distributed. The "-1." accounts for their charge:
    Code:
        # Calculate number density, normalised so 1 when uniform
        density = calc_density(pos, ncells, L)
    
        # Subtract ion density to get total charge density
        rho = density - 1.
    [code]Code tags[/code] are essential for python code and Makefiles!
  20. #11
  21. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0
    This is absolutely insane.. Seriously I am in your debt. I'm taking a good look at the code now, so that I will understand it and replicate it.

    I will report back with my musings!
  22. #12
  23. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    I don't understand how to use the devshed mail system. If you sent me an address where I can send the report I've got, well, in any case, I didn't receive it. You can send me mail to an account I read by request,

    b49p23tivg@yahoo.com

    As far as the changes I made, I'd think you understand the offset and other computations.

    The mathematical equivalent of zip is approximately transpose.

    I studied plasma physics (30 years ago), so it's not completely mystical dumb luck that I just read descriptions of PROBEPIC. But almost.
    [code]Code tags[/code] are essential for python code and Makefiles!
  24. #13
  25. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0
    Originally Posted by b49P23TIvg
    I don't understand how to use the devshed mail system. If you sent me an address where I can send the report I've got, well, in any case, I didn't receive it. You can send me mail to an account I read by request,

    b49p23tivg@yahoo.com

    As far as the changes I made, I'd think you understand the offset and other computations.

    The mathematical equivalent of zip is approximately transpose.

    I studied plasma physics (30 years ago), so it's not completely mystical dumb luck that I just read descriptions of PROBEPIC. But almost.
    The code all works fine.

    The reason why the Histogram version didn't work (although much faster) was because it did not smoothly allocate the particles to cells depending on their offset from the integer cells - it just binned them irrespective of whether or not they were slowly varying across the cells (requiring the evolution of the offset).

    This method is definitely more robust, thanks again.

    My next personal attempt to improve the overall running of the code was to compile the density calculation function into Cython and then call it out as a function.

    Apparently this saves alot of work and time.

    However I am getting an error when trying to compile the Cython file. This is what I get in my command line:

    Code:
    density_calc.pyx:23:18: undeclared name not builtin: numpy
    building 'density_calc' extension
    error: Unable to find vcvarsall.bat
  26. #14
  27. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,995
    Rep Power
    481
    I'd try scipy.weave . I've never used it and won't be able to try it until tomorrow.
    [code]Code tags[/code] are essential for python code and Makefiles!
  28. #15
  29. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jan 2013
    Posts
    22
    Rep Power
    0
    Hmm, in-lining C/C++ within Python code. Powerful for speeding up and optimization. Attempting it now, will let you know what I've done.
Page 1 of 2 12 Last
  • Jump to page:

IMN logo majestic logo threadwatch logo seochat tools logo