0
$\begingroup$

I am trying to understand the attached python 3X code.

The following dispersion relation is given

$$\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]$$

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.

To do so the code below should be checked in order to understand the plasma behavior either side of $L_{min}$ threshold, ie $L_{min} < L$ and $L < L_{min}$. To do so we should: 1) Adjust the number of simulation particles to grid points to improve the statistics, 2) Fix the number of particles per cell (i.e. npart/ngrid) and keep it $>> 1$, in order to reduce numerical noise and the runtime needed to observe the instability can be estimated from the maximum growth rate.

The theory is OK but I am struggling to actually find the minimum grid length $L_{min}$ once given the code. Might you please shed some light and guide me through? Thank you.

#! /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)
$\endgroup$
2
  • 2
    $\begingroup$ Unless you bump into someone who is an expert on plasma physics simulations it is very hard to help you. It is not clear what numerical method you use to simulate what particular equation / model. Furthermore, providing a bunch of code without any comments or explanations in the text is also not very helpful. Try to link your questions to concepts the community is more familiar with. $\endgroup$
    – Dan Doe
    Jan 2 at 18:16
  • $\begingroup$ @DanDoe you are absolutely right, I am thinking how to edit my question in the most efficient way to make it more to the point. $\endgroup$
    – JD_PM
    Jan 2 at 21:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.