# Simple Particle In Cell Code in Matlab

### Introduction

This post includes the sample Particle-In-Cell (PIC) code that goes with our previous article on the electrostatic particle-in-cell method. It simulates the flow of uniform plasma past an obstruction – a charged plate in our case. This simple code is implemented in Matlab. If you don’t have Matlab, you can run the code using Octave, which is an open-source Matlab emulator. If you do end up getting Octave, as of the date of this article, the MinGW build available on Octave-Forge is out of date and contains a buggy contouring function. In order to use the latest version – and see the plots – you need to install Cygwin and add Octave as a Cygwin package…

### File Organization

The sample PIC code is contained in two files: flow_around_plate.m and eval_2dpot_GS.m. The first file is the actual PIC code, while the second file holds the potential solver. The solver computes electric potential using the Finite Difference method. The code simulates the flow of solar wind (low density plasma) over a charged plate. The code is not particularly fast (it’s in Matlab after all!), but shouldn’t take more than few minutes to run on a standard PC. If you run the code through Octave, make sure to first execute `more off` to avoid pagination of the output.

### Input Parameters

The simulation plasma environment is set on lines 24-28:

n0 = 1e12; %density in #/m^3 phi0 = 0; %reference potential Te = 1; %electron temperature in eV v_drift = 7000; %ion injection velocity, 7km/s phi_p = -7; %plate potential |

By changing the plate potential, you can alter the size of the wake that forms behind the plate. At sufficiently large values, ions are attracted back to the plate. This behavior is analogous to the plasma environment that forms around satellites utilizing electric propulsion. Slow moving ions are extracted out of the thruster plumes and can impact surfaces without a direct line of sight to the thruster.

The simulation domain size and the number of time steps is set starting with line 35:

nx = 16; %number of nodes in x direction ny = 10; %number of nodes in y direction ts = 200; %number of time steps |

At the default settings, 200 steps is sufficient to reach steady state. The simulation time step is such that a particle traveling at the drift velocity will cross 10% of the cell length per time step. Since particles speed up as they approach the plate, this conservative time step keeps particles from traveling more than a cell per time step. This is a necessary constraint in all PIC codes, but in general you want to keep the particles traveling a shorter distance (I generally use 25% in my simulations). The number of particles generated at each time step is specified on line 39,

np_insert = (ny-1)*15; %insert 15 particles per cell |

Plate dimension are set on lines 48 and 49:

box(1,:) = [floor(nx/3) floor(nx/3)+2]; %x range box(2,:) = [1 floor(ny/2)]; %y range |

For an explanation of the physics used in the code, please see the Particle In Cell page. Feel free to leave a comment if you have any questions.

### Download

Download the files here: flow_around_plate.m, eval_2dpot_GS.m

Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.
Please leave a comment if you have any questions. Thanks!

Hi Lubos;

in Matlab code your code the wall potential is given phi_p = -7V. Is there any special meaning under this or was it chosen randomly,

and, what is the result of choosing positive wall voltage?

Thank you..

Best regards,

Hi Ozgur,

Yes, this is just some arbitrary value that produced nice looking plot. The potential is negative since the potential in the bulk plasma is zero and hence this attracts the ions towards the plate. If you chose a positive value, the ions will be repelled away.

Our professor assigned us an assignment on Particle in Cell simulation for his Plasma Engineering graduate class. We are supposed to use your code as a basis. After altering your code a bit, we basically moved and extended the wall/barrier to the right edge and changed a few parameters.

Is there something in your code that generates a decrease in the density on the left hand side of your density plot? Since there’s a uniform injection of ions, shouldn’t there be a uniform density until it reaches the sheath region? We just can’t really understand why the left hand side is also going to zero since it’s technically not supposed to do that. In your density plot online, it seems that the density goes to zero, as well, so maybe we’re not wrong, I just don’t understand why.

Hi Ashley, that has to do with the way particles are injected into the domain. The particles are born at a random “x” position in the first cell and also with a random velocity. So you basically end up with a non-uniform density in this first cell. You can avoid this by loading particles at x=0 and also doing the density update prior to the push. I left it this way since I didn’t think it made much difference on the solution at the plate location.

Hi Lubos,

i try to simulate plasma flows by using Particle in Cell method. The situation is very similar to your matlab code.;(i.e. collisionless and quasi neutral plasma,and the model is hybrid). i have finished the code in non-rectangular domain. But i could not verified the results whether true or not. Can you help me about verification.

Thank you very much for your attention.

Sure Ozgur, what did you have in mind? Some ways you can also try to verify your code is by testing the potential solver with some analytical function and also checking for energy conservation.

A PIC code in Matlab could be adapted to run in a computer cluster? and/or adapted to run in parallel?, Thanks

Hi Sergio, if I am not mistaken, Matlab now supports parallel processing. However, I am not familiar with it – we do most of our actual development either in C/C++ or Java. Matlab is a great prototyping language, but if you are after speed, you’ll be better off using a non-interpreted language.

Hi…I am new in PIC..I want some help how to do assignment from particle position to grid position…how I can implement…

Also how I may know how many particle are there in one particular cell of the whole geometry?

I want to interpolate particle properties like velocities and charge to my grid (one dimensional only)

As in each cell there are so many particles so each has it own contribution on particular node (let say k) depending on weight factor so I want to know how I can sum all those properties on a node (k and k+1).

Why you have add 1.0 in the assignment part (fi = 1+part_x(p,1)/dh; )?

Please elaborate the assignment part.

Please help me…

Regards

Amitavo, you basically want to compute the average velocity. This is no different from computing the regular average, with the difference that you are averaging fractional contributions. When you scatter data to the grid, you place a fraction of the velocity on the left node, and the rest on the right node, assuming 1D problem. So as an example, let’s say the velocity is 100 m/s. Perhaps the particle is really close to the right node, so the right node may get 90 m/s, while only 10 m/s is deposited to the left node. Next imagine that this cell contains 5 particles, and that they are all located at the same spatial position. You will end up with 450 m/s total velocity on the right node and 50 m/s on the left node. If you were to divide the velocity sum by 5, you will get that the velocity varies from 10 m/s on the left end to 90 m/s on the right end. This is clearly incorrect since all particles in the cell are at 100 m/s. To get the correct answer you need to divide the sum by the fractional particle count. This is why you want to also keep a “count” variable on the grid nodes to which you scatter the value of 1 (if all particles have equal weight) for each particle. If particles have a variable weight than you scatter the respective specific weight. You then divide each node value by the count, i.e.

vel[i]=vel_sum[i]/count[i];. In this example, the counts would be 4.5 on right, and 0.5 on left, giving you 100 m/s on both the left and the right node.The second question about adding 1 is due to Matlab starting its indexing at 1 instead of 0. If your code is in C or Java, you don’t add this 1.

tHANKS FOR YOUR NICE REPLY. YES NOW IT IS UNDERSTAND THAT I AM WRONG IN AVERAGING.

BY “COUNT”..I UNDERSTAND THAT I HAVE TO MAKE COUNT OF NO. OF PARTICLES CONTRIBUTION OF PARTICULAR GRID…

BUT STILL “COUNT” STATEMENT IS NOT VERY CLEAR TO ME SO THAT I CAN IMPLEMENT IT….COULD YOU PLEASE MAKE CLEAR TO ME WITH AN SMALL EXAMPLE….IT HELPS ME TO CLEAR MY UNDERSTANDING…

It’s no different than scattering velocity, except that you are scattering particle weight. If all particles have the same weight, you simply scatter the value of 1.0 every time you scatter the velocity. If they have different weights, you scatter the particle weight, or the weight divided by a reference weight (to normalize, since weights are typically very large numbers).

Thanks a lot….

I am facing problems in assingment process after pushing the particles.

I am not able to implement the particle properties from the new position of the particles to the grid point which is shared by two cells.

After the localization process for the new particle positions how I can implement the assignment operation?

Please guide me as I am new in the field.

As an example: Suppose we have 10 particles in (K,K+1) cell and 5 particles in (K-1,K) cell after pushing then I want to assign velocities and charge to K node due to all these particles. Please suggest me/ help me how can I perform this operation.

Regards

Hi,

I have found your writings very useful. I am learning plasma out of interest. I have a doubt: if I have a 2d plasma inside 4 walls of 0 Volt each and I want to plot the potential (calculated as in your tutorial, electrons as Boltzmann particle, ions : moving macroparticles), why do I get an inverted Gaussian type curve? I thought potential should be uniform across a plasma! Thanks.

Hi John, that would be the case if the governing Poisson equation reduced to the zero forcing function Laplace form, . However, that is not the case. Even though plasma is quasi-neutral, this condition does not hold in the vicinity of walls. Near the walls you will get what’s known as a sheath, where (and hence ) since electrons are absorbed by the wall faster due to their higher thermal speed. This is somewhat analogous to the boundary layer in aerodynamics. Even if you have a uniform velocity flow, the flow will not be uniform near surfaces, since you need to have zero velocity at the surface boundary due to the no-slip condition. It is this non-zero charge density term that gives you the forcing function that gives you the Gaussian-like distribution. Hope this helped.

Am working on simulation of DBD using particle in cell method.PLS Sir how will i go about it.

That’s not a particularly easy problem. I must admit I do not have much background with modeling discharges (most of my work has been in plume modeling), however what I imagine you would need is a code that takes into account dielectric charging, surface deposition, as well as a model for breakdown. Have you seen the excellent book by Lieberman and Lichtenberg, Principles of Plasma Discharges and Materials Processing? It goes into details of various plasma discharges.

Hi,

Thanks for the explanation. If the electrons are also simulated as macroparticles like ions (with same specific weight but with opposite charge), what is to be done for the sheath formation in the simulation? Or, will it be inherent? How the time-step will be defined then?

Secondly, in the Matlab code for flow around a plate, why the factor 0.5 is missing for x-component of velocity in the Particle Generation section of the code? It is there for y-component of velocity.

The sheath will form self-consistently. It forms since the electrons are more mobile (lower mass -> higher velocity) and thus will be absorbed by the wall faster than the ions. This will in turn make the wall negative in respect to the bulk plasma (or the plasma potential will go up if the wall potential is held fixed). This potential gradient acts to reduce the electron flux. The wall potential will continue to go down until the electron current is balanced by the ion current. About your second question – good catch! I have no idea. It could very well be a typo. The way it’s written you will end up with anisotropic temperature. I wrote the code few years ago so I don’t quite remember what I was thinking back then.

Thanks for your reply. Also, thank you for maintaining such an active website.

If you had an unstructured triangular mesh instead of a mesh with quadrilateral elements, what would be a good algorithm for calculating the electric field at each of the mesh nodes. For the scatter and gather operations I was thinking you could use barycentric coordinates of a triangle.

http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)

But once you have the potential at each of the nodes on an unstructured mesh, what would be a good way to calculate the E field at each node.

This is really the foundation of the finite element method. Typically what you do is you define a basis element for each cell. These elements can be of different order. The derivative will be one order less. Quite typically, you may see linear elements used for potential which then give you zero-order (constant) electric field in the cell.

I tried implementing this PIC code using the finite element method solver from this blog post.

http://www.particleincell.com/2012/matlab-fem/

Here are some results I got.

https://dl.dropbox.com/u/5095342/Images/figure1.jpg

https://dl.dropbox.com/u/5095342/Images/figure2.jpg

https://dl.dropbox.com/u/5095342/Images/figure3.jpg

However, in order to get a stable result I needed to set

Te = 0.125; %electron temperature in eV

instead of Te = 1.0; . When I do this I get a result that looks similar to what the original PIC code from this blog post generates with Te set to 0.125 .

Do you know why the solver result might become unstable when Te is set higher than 0.125 ?

I used the following in the MIT FEM code for f(x,y).

% multiply Fe by f at centroid for load f(x,y): one-point quadrature!

Fe=Area/3*((mean(den(nodes,1))-n0*exp((mean(phi(nodes,1))-phi0)/Te))*QE/EPS0);

Does this look correct to you?

Hi Dr. Brieda

I have some question about ES-PIC code. Please help me.

In your code, was electron and ion isothermal together? otherwise What is your scientific reason?

How have been calculated time step (10^-12) for electrons?

Why the drift velocity has been considered only for ions? note that 10^5 (order of thermal velocity of electrons) isn’t much larger than 10^3 (order of drift velocity)

in calculation thermal velocity , vth=k*Te/m was used, in spite of temperature of electron was in ev not Kelvin!

best regards

Hi Mahzan, these assumptions are typically taken in the context of plume modeling. The reason is that including electrons as kinetic particles (or even a fluid) complicates the numerical model unnecessarily. If magnetic field is not present, the electron density follows the Boltzmann relationship.

Thanks for catching that vth mistake, I have corrected it. I have also specified a different temperature for ions to indicate that generally LEO plasma will not be isothermal.

In eval_2Dpot_GS.m you set the neumann boudary conditions around the border as follows.

%set boundaries

b(1:nx) = 0; %zero electric field on y=0;

b(nn-nx+1:nn) = 0; %zero electric field on y=L;

b(nx:nx:nn) = 0; %zero electric field on x=L;

b(1:nx:nn) = phi0; %fixed potential on x=0;

But in flow_around_plate.m you set the electric field around the border as follows.

efx(1,:) = 2*(phi(1,:) – phi(2,:)); %forward difference on x=0

efx(nx,:) = 2*(phi(nx-1,:) – phi(nx,:)); %backward difference on x=Lx

efy(:,1) = 2*(phi(:,1) – phi(:,2)); %forward difference on y=0

efy(:,ny) = 2*(phi(:,ny-1) – phi(:,ny)); %forward difference on y=Ly

Should the normal electric field on the borders y=0, y=L, x=L just be set to zero to match the boundary conditions?

efx(nx,:) = 0;

efy(:,1) = 0;

efy(:,ny) = 0;

Great example! Thank you.

I downloaded all files and ran them at the Octave. I received strange results: at every time step I have 0 particles and blank graphic.

Have you any suggestions about how I can solve my problem?

I’m sorry for dumb question..

I found where was a problem.

I downloaded all files and ran them at the Matlab. At every time step I have 0 particles and blank graphic.

Hi Sudhan, thanks for catching this. The equation for thermal velocity must have gotten screwed up at some point. I fixed it. Try it again and let me know if it still doesn’t work.

Thanks for your sharing. I was looking for MPM on the website to solve some normal fluid problems. Then I got here. I saw some interesting results from MPM, for example: http://grantkot.com/MPM/Liquid.html. Interesting, haha.

Thanks for the link!

Is it possible to incorporate magnetic field to this code??

Another thing, How I can plot the velocity distribution of the particles from this code?

Yes, you can have a static magnetic field. See the article on the Boris method for details of implementing the particle pusher.

You will generally need to write your own code to compute the VDF. Just create a number of bins and iterate over the particles and increment the bin counters based on their speed. Then you can normalize the bins if needed by the max value or the sum.

I have tried to plot the electric field using quiver command from the data efx and efy and got a result. Problem is that I am getting a swirl kind of structure in the electric field. which is not acceptable as i have also introduced two parallel plates instead of one. One I have kept at y=Ly and another at y=0. having width $ \delta $ x = 3

Hi, first of all I would like to say this is a great blog. I had been reading and testing some of the codes you posted here.

I am doing a thesis on an inertial electrostatic confinement reactor (basically have to model flow of plasma due to a static electric field). I have two concentric (mostly transparent) grids held at about 60KV potential. Gases are ionized at the outer grid and accelerate towards the central grid. In the center some undergo fusion. Ions will also undergo recirculations (if they don’t collide with the central grid).

So I am wondering if it will be possible to alter the geometries of this code (perhaps 2D circular geometry), and have grids rather than plates? Would it be viable to model an IEC device using this method? I would greatly appreciate any help you can provide.

Kind regards,

Jimmy Zhan

Hi Jimmy, you should check out our Starfish code (still in development). It is a general 2D solver that can be used for a number of different geometries.

Great blog,I have been working on the codes posted here. The link you posted here was very useful.I am learning plasma out of interest. Thank you so much.

How this particle in cell simulation may be used for “self-focusing of laser beam in plasmas”

Hi,

This is really a very fascinating site on PIC. With the above example of 2D PIC, is it possible to calculate current recorded by the plate as a function of time ? Can we use the formula (charge hitting the plate in time dT) / dT to calculate current at a particular time instant ? Since it is a 2D simulation, what will the current mean physically ? Is this the current collected by a plate of given x and y dimensions but with UNIT DIMENSION along Z (i.e. 1 meter) ? Your thought on this will be very helpful. Thanks!

Hi Prakash, yes, you would basically count the number of particles hitting the plate and average over some time. Since this is an ion-only simulation, you will also need to include the electron term, which for can be modeled analytically from the planar sheath equation. For a 2D XY simulation like this one, the current will be current per unit depth as you mentioned. If this were an axisymmetric simulation, then you could compute the actual surface area by considering the body of revolution.

Thanks a lot. This is the only elaborately-explained resource on the web on Particle-In-Cell. Please keep it up.

Hi,I would like to say this is a great blog. It helps me a lot in my final project.

I have tried to understand your code. About Debye length calculation,I could not understand lD = sqrt(EPS0*Te/(n0*QE)); %Debye length. According to formula, I thought is lD = sqrt(EPS0*Te*K/(n0*QE*QE)). Can you explain it, please? Thank you!

Hi Helen, this is because in the above example, electron temperature Te is in electron volts. The conversion from Kelvin to electron volts is k*Te(K)/QE = Te(eV).

Thanks,Lubos. Can i ask about u=(j-1)*nx+1? what is this for ?

hi

can i use this code for higher densities (about 10^20 #/m^3) ?

is this code suitable for tokamak plasma simulation?