Hello all.

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

All the parts highlighted in Red are what I struggled to write with my limited knowledge.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 librarytry: 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()

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

Tweet This+ 1 thisPost To Linkedin