# How to determine the minimum grid length | Numerical Plasma Physics

• JD_PM
In summary: X (cm)") ylabel("Y (cm)") title("Track of particle (X, Y)") fig = plt.figure() fig.subplots_adjust(left = 0.5, right = 1.5, bottom = 0.5, top = 1.5) ax = fig.add_subplot(211, projection='3d') # Axes for x- and y-ax
JD_PM
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.

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)

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')

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

#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

# 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)
Thank you!

Last edited by a moderator:
But that code is nothing to do with plasma simulation it is a
JD_PM said:
Python:
 #  Python script for computing and plotting single charged particle
#  trajectories in prescribed electric and magnetic fields.
It doesn't use particle-in-cell or any other mesh algorithm, it appears to use a leap-frog integrator for the ODEs for a single particle.

Hi pbuk.

pbuk said:
It doesn't use particle-in-cell or any other mesh algorithm, it appears to use a leap-frog integrator for the ODEs for a single particle.

Honestly, I do not have the knowledge to argue that back.

Let me summarize my question, given that it might be uncomfortable to be presented with such a long post.

I am interested in finding the minimum grid length ##L_{min}## needed to support unstable modes. To do so the code should be checked in order to understand the plasma behavior either side of ##L_{min}## threshold, ie ##L_{min} < L## and ##L < L_{min}##.

I am given the following hints

• Adjust the number of simulation particles to grid points to improve the statistics.
• Fix the number of particles per cell (i.e. npart/ngrid) and keep it ##>> 1##, in order to reduce numerical noise.
• The runtime needed to observe the instability can be estimated from the maximum growth rate.
If you guys could guide me through I would really appreciate it

You might point out I should deal with some issues first, such as where we expect an instability to occur given the dispersion relation (in other words, finding the range of wave numbers ##k## for which the oscillation frequency is imaginary). If you need me to address that first let me know and I will deal with such smaller question in a different forum :)

Where in the code you posted do you enter ## L ##?

JD_PM
pbuk said:
Where in the code you posted do you enter ## L ##?

A general grid length ##L## is not introduced. I guess we have to do so in the "main program", but I am stuck in how to proceed.

JD_PM said:
I am given the following hints

• Adjust the number of simulation particles to grid points to improve the statistics.
• Fix the number of particles per cell (i.e. npart/ngrid) and keep it ##>> 1##, in order to reduce numerical noise.
• The runtime needed to observe the instability can be estimated from the maximum growth rate.
I do not see how the hints can help.

Here's a summary of what the code contains at the moment:
• A function to integrate particle trajectory in given ##E, B## fields.
• A Lorentz force law solver up to 2nd order accuracy in dt.
• How to make 2D plots of particle orbit.
• How to make 3D plot of particle orbit.
• "Main program"

Let me share the original statement

So question 2 should be solved first; let me share that question's statement too (together with 1's, given that 1 and 2 are related)

I have been dealing with question 2 for days and this is what I got.

If you see it is really unavoidable to properly solve 2. first let me know and I will ask it in PF this time

JD_PM said:
I have been dealing with question 2 for days.
Then you have been wasting your time, as well as that of stackexchange members.

Didn't you notice the ## = 0 ## at the end of ref ## (1) ##, or did you think it was unimportant? And in any case why are you trying to use the dispersion function ## \varepsilon (k, \omega) = 0 ## given in question 1 to find unstable values of ## k ##? Question 2 states that you should be using the dispersion relation derived as the answer in question 1.

OK, we are dealing with #2 here.

pbuk said:
And in any case why are you trying to use the dispersion function ## \varepsilon (k, \omega) = 0 ## given in question 1 to find unstable values of ## k ##? Question 2 states that you should be using the dispersion relation derived as the answer in question 1.

You are right, thanks for pointing it out! I used the dispersion relation (as can be checked in the other thread) and obtained the following range

$$-\left| \frac{\omega_p}{v_0}\right| < k < \left| \frac{\omega_p}{v_0}\right|$$

The corresponding range for ##\alpha## is ##-1 < \alpha < 1##.

We want to relate the wave number ##k## to a (minimum grid) length. We know that the wavelength ##\lambda## relates to the wave number ##k##

$$\lambda = \frac{2 \pi}{k}$$

So assuming that with "grid length" they meant "wavelength associated to the wave number ##k##", then the minimum grid length should be

$$L_{min} > 2 \pi \left| \frac{v_0}{\omega_p}\right|$$

So that we can actually see the grid.

OK but now the question is how to implement this into the code.

## 1. How do you calculate the minimum grid length for numerical plasma physics simulations?

The minimum grid length for numerical plasma physics simulations can be calculated using the Debye length, which is a fundamental parameter in plasma physics. It is given by the formula λD = √(ε0kBTe2/ne), where ε0 is the permittivity of free space, kB is the Boltzmann constant, Te is the electron temperature, and ne is the electron density. This length scale determines the spatial resolution required to accurately capture the plasma dynamics.

## 2. What factors affect the minimum grid length in numerical plasma physics simulations?

The minimum grid length in numerical plasma physics simulations is affected by several factors, including the plasma parameters (such as temperature, density, and magnetic field strength), the type of simulation (e.g. fluid or particle-in-cell), the desired accuracy, and the computational resources available.

## 3. How does the minimum grid length affect the accuracy of numerical plasma physics simulations?

The minimum grid length plays a crucial role in determining the accuracy of numerical plasma physics simulations. A smaller grid length allows for a more detailed representation of the plasma dynamics, resulting in a more accurate simulation. However, using a smaller grid length also requires more computational resources and may lead to longer simulation times.

## 4. Can the minimum grid length be adjusted during a simulation?

In most cases, the minimum grid length cannot be adjusted during a simulation. Once the simulation parameters are set, including the grid length, the simulation must be run until completion. However, there are some adaptive mesh refinement techniques that can dynamically adjust the grid resolution in certain regions of the simulation domain.

## 5. How can the minimum grid length be validated in numerical plasma physics simulations?

The minimum grid length can be validated by comparing the simulation results with analytical solutions or experimental data. If the simulation results closely match the expected behavior, then the chosen minimum grid length is likely appropriate. Additionally, convergence studies can be performed by varying the grid length and observing how the results change. A smaller grid length should result in more accurate and consistent results, up to a certain point where further refinement does not significantly improve the simulation.