- #1
JD_PM
- 1,125
- 156
Summary:: I am learning particle-in-cell (PIC) python 3X code. PIC currently represents one of the most important plasma simulation tools. It is particularly suited to the study of kinetic or non-Maxwellian effects.
I am learning particle-in-cell (PIC) python code. PIC currently represents one of the most important plasma simulation tools. It is particularly suited to the study of kinetic or non-Maxwellian effects.
Given the following dispersion relation
$$\epsilon (k, \omega) = 1 - \frac 1 2 \left[ \frac{\omega_p^2}{(\omega-kv_0)^2} + \frac{\omega_p^2}{(\omega+kv_0)^2}\right]$$
I found the range of wave numbers $k$ for which the oscillation frequency is imaginary to be ##-\left|\frac{w}{v_0}\right| \lt k \lt \left|\frac{w}{v_0}\right|##
What I am trying to understand is how to find the minimum grid length ##L_{min}## as a function of ##\frac{v_0}{w}##. ##L_{min}## indicates the needed minimum grid length to support such unstable modes.
I think we should be able to study the plasma behaviour for both ##L < L_{min}## and ##L > L_{min}##. I was told I should adjust the number of simulation particles to grid points to improve the statistics. Besides, the number of particles per cell (i.e. npart/ngrid) should be fixed and should be much greater than 1, in order to reduce numerical noise. The runtime needed (here in units of ##\omega_p^−1##) to observe the instability can be estimated from the maximum growth rate.
Here's the full python 3 code I am working with. Please note I have little experience with coding so I might ask lots of follow up questions.
Thank you! 
I am learning particle-in-cell (PIC) python code. PIC currently represents one of the most important plasma simulation tools. It is particularly suited to the study of kinetic or non-Maxwellian effects.
Given the following dispersion relation
$$\epsilon (k, \omega) = 1 - \frac 1 2 \left[ \frac{\omega_p^2}{(\omega-kv_0)^2} + \frac{\omega_p^2}{(\omega+kv_0)^2}\right]$$
I found the range of wave numbers $k$ for which the oscillation frequency is imaginary to be ##-\left|\frac{w}{v_0}\right| \lt k \lt \left|\frac{w}{v_0}\right|##
What I am trying to understand is how to find the minimum grid length ##L_{min}## as a function of ##\frac{v_0}{w}##. ##L_{min}## indicates the needed minimum grid length to support such unstable modes.
I think we should be able to study the plasma behaviour for both ##L < L_{min}## and ##L > L_{min}##. I was told I should adjust the number of simulation particles to grid points to improve the statistics. Besides, the number of particles per cell (i.e. npart/ngrid) should be fixed and should be much greater than 1, in order to reduce numerical noise. The runtime needed (here in units of ##\omega_p^−1##) to observe the instability can be estimated from the maximum growth rate.
Here's the full python 3 code I am working with. Please note I have little experience with coding so I might ask lots of follow up questions.
Python:
#! /usr/bin/python
#
# Python script for computing and plotting single charged particle
# trajectories in prescribed electric and magnetic fields.
# Roughly equivalent to boris.m MATLAB program
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.widgets import Slider, Button, RadioButtons
from mpl_toolkits.mplot3d import Axes3D
import os
import os.path
import sys
from sys import exit
from time import sleep
# ===================================
#
# Function to integrate particle trajectory
# in given E, B fields
#
# ===================================
def integrate(E0, B0, vz0):
global dt, v0, x0, xp, yp, zp, qom, larmor, nsteps
wc=qom*B0 # cyclotron frequency
larmor=vperp/wc
print ("Cyclotron frequency =",wc)
print ("Perpendicular velocity v_p=",vperp)
print ("Larmor radius=",larmor)
norm = 1. # choose whether to normalise plot axes dimensions to Larmor radius
trun=5*2*np.pi/wc # total runtime
dt=.1/wc # timestep - adjust to current B-field
nsteps=int(trun/dt) # timesteps
E=np.array([0.,E0,0.]) # initial E-field
B=np.array([0.,0.,B0]) # initial B-field
u=np.array([0.,0.,0.]) # intermediate velocity
h=np.array([0.,0.,0.]) # normalized B-field
xp[0]=x0[0]
yp[0]=x0[1]
zp[0]=x0[2]
v0[2]=vz0 # z-component
v=v0+.5*dt*qom*(E+np.cross(v0,B)) # shift initial velocity back 1/2 step
x=x0
for itime in range(1,nsteps):
x=x+dt*v
xp[itime]=x[0] /norm
yp[itime]=x[1] /norm
zp[itime]=x[2] /norm
tp[itime]=itime*dt
#
# Boris mover: solves dv/dt = q/m*(E + vxB) to 2nd order accuracy in dt
#
qomdt2 = dt*qom/2
h = qomdt2*B
s=2*h/(1+np.dot(h,h))
u = v + qomdt2*E
up=u+np.cross(u+np.cross(u,h),s)
v=up+qomdt2*E
# vxp[itime] = v[0]
# ===================================
# Make 2D plots of particle orbit
#
# ===================================
def plot_track2D():
global xp,yp,nsteps,ax1
fig = plt.figure(figsize=(8,8)) # initialize plot
xmin=np.min(xp)
xmax=np.max(xp)
ymin=np.min(yp)
ymax=np.max(yp)
fig.add_subplot(221) # 1st subplot in 2x2 arrangement
plt.cla()
plt.grid(True, which='both')
plt.xlim( (xmin, xmax) )
plt.ylim( (ymin, ymax) )
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.plot(xp[0:nsteps],yp[0:nsteps],c='b')
fig.add_subplot(222) # 2nd subplot
# fig.add_subplot(223) # 2nd subplot
# fig.add_subplot(224) # 2nd subplot
plt.draw()
plt.savefig('./particle_orbit.png') # Save plot to file
# ===================================
#
# Make 3D plot of particle orbit
#
# ===================================
def plot_track3D():
global xp,yp,zp,nsteps,ax1
xmin=np.min(xp)
xmax=np.max(xp)
ymin=np.min(yp)
ymax=np.max(yp)
zmin=np.min(zp)
zmax=np.max(zp)
ax1.cla()
plt.ion()
plt.grid(True, which='both')
ax1.set_xlim( (xmin, xmax) )
ax1.set_ylim( (ymin, ymax) )
ax1.set_zlim( (zmin, zmax) )
ax1.set_xlabel('$x $ [m]')
ax1.set_ylabel('$y $ [m]')
ax1.set_zlabel('$z $ [m]')
#ax1.set_aspect(1.)
ax1.scatter(xp,yp,zp,c=tp,marker='o') # tracks coloured by elapsed time since start
plt.draw()
# =============================================
#
# Main program
#
# =============================================
print ("Charged particle orbit solver")
plotboxsize = 8.
animated = True x0=np.array([0.,0.,0.]) # initial coords
vz0=0.
v0=np.array([-1e2,0.,vz0]) # initial velocity
vperp = np.sqrt(v0[0]**2+v0[2]**2)
E0=0.
B0=.1
e=1.602176e-19 # electron charge
m=9.109e-31 # electron mass
qom=e/m # charge/mass ratio
wc=qom*B0 # cyclotron frequency
larmor=vperp/wc
print (wc,vperp,larmor)
trun=5*2*np.pi/wc # total runtime
dt=.1/wc # timestep - adjust to current B-field
nsteps=int(trun/dt) # timesteps
B1=np.array([0.,0.,0.1]) # gradient B perturbation
#wc=qom*np.linalg.norm(B) # cyclotron frequency
#nsteps=2
tp = np.zeros(nsteps) # variables to store particle tracks
xp = np.zeros(nsteps)
yp = np.zeros(nsteps)
zp = np.zeros(nsteps)
vxp = np.zeros(nsteps)
vyp = np.zeros(nsteps)
vzp = np.zeros(nsteps)
# Compute orbit
integrate(E0, B0, vz0)
# 2D orbit plotter
plot_track2D()
exit(0) # Quit script before 3D plot - comment out to continue!
# Start 3D interactive mode with sliders for B, E and v0
plt.ion() # Turn on interactive plot display
fig = plt.figure(figsize=(8,8))
# Get instance of Axis3D
ax1 = fig.add_subplot(111, projection='3d')
# Get current rotation angle
print (ax1.azim)
# Set initial view to x-y plane
ax1.view_init(elev=90,azim=0)
ax1.set_xlabel('$x $[microns]')
ax1.set_ylabel('$y $[microns]')
ax1.set_zlabel('$z $[microns]')
plot_track3D()
#filename = 'a0_45/parts_p0000.%0*d'%(6, ts)
#plot_from_file(filename):
axcolor = 'lightgoldenrodyellow'
axe0 = fig.add_axes([0.1, 0.95, 0.3, 0.03])#, facecolor=axcolor) # box position, color & size
axb0 = fig.add_axes([0.5, 0.95, 0.3, 0.03])#, facecolor=axcolor)
axv0 = fig.add_axes([0.1, 0.9, 0.3, 0.03])#, facecolor=axcolor)
sefield = Slider(axe0, 'Ey [V/m]', -5.0,5.0, valinit=E0)
sbfield = Slider(axb0, 'Bz [T]', -1.0, 1.0, valinit=B0)
svz = Slider(axv0, 'vz [m/s]', 0.0, 1.0, valinit=0.)
def update(val):
E0 = sefield.val
B0 = sbfield.val
vz0 = svz.val
integrate(E0,B0,vz0)
plot_track3D()
plt.draw()
sefield.on_changed(update)
sbfield.on_changed(update)
svz.on_changed(update)
resetax = fig.add_axes([0.8, 0.025, 0.1, 0.04])
button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')
def reset(event):
global ax1
sefield.reset()
sbfield.reset()
svz.reset()
ax1.cla()
ax1.set_xlabel('$x $[microns]')
ax1.set_ylabel('$y $[microns]')
ax1.set_xlim( (0., 10.) )
# ax1.set_ylim( (-sigma, sigma) )
ax1.grid(True, which='both')
plt.draw()
button.on_clicked(reset)
#plt.show()
plt.show(block=False)
Last edited by a moderator:
