### Thread: Optimizing PIC (Particle in Cell) Code

Page 1 of 2 12 Last
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:
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)
self.phase_plot = ax.plot(pos, vel, '.')[0]
ax.set_title("Phase space")

self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]

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

"""
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):
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])

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.

2. Please, for whom are you working?
3. 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.
4. 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.
5. 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) ???
6. 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.```
Last edited by b49P23TIvg; January 29th, 2013 at 01:18 PM.
7. 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:
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)
self.phase_plot = ax.plot(pos, vel, '.')[0]
ax.set_title("Phase space")

self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]

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

"""
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):
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])

#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.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.
8. 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.
9. 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:
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)
self.phase_plot = ax.plot(pos, vel, '.')[0]
ax.set_title("Phase space")

self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]

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

"""
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):
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])

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)```
10. #### 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.```
11. 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!
12. 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.
13. 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```
14. I'd try scipy.weave . I've never used it and won't be able to try it until tomorrow.
15. 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