PIC-C, Particle In Cell Consulting, LLC. Plasma physics and rarefied gas simulation codes.

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.

plasma potential plasma density
Figure 1. Simulation results, plasma potential and ion density

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.

49 comments to “Simple Particle In Cell Code in Matlab”

  1. November 28, 2011 at 8:20 am

    Please leave a comment if you have any questions. Thanks!

  2. Ozgur Tumuklu
    December 1, 2011 at 2:37 am

    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,

    • December 1, 2011 at 7:00 am

      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.

  3. Ashley
    December 1, 2011 at 7:02 am

    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.

    • December 1, 2011 at 7:03 am

      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.

  4. Ozgur Tumuklu
    January 24, 2012 at 12:35 pm

    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.

    • January 24, 2012 at 2:09 pm

      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.

  5. Sergio Flores
    April 21, 2012 at 11:38 pm

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

    • April 22, 2012 at 8:51 pm

      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.

  6. Amitavo
    April 25, 2012 at 9:25 am

    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

    • April 27, 2012 at 5:28 am

      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.

  7. Amitavo
    April 27, 2012 at 7:16 am

    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…

    • April 27, 2012 at 8:09 am

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

      • Amitavo
        April 30, 2012 at 12:24 am

        Thanks a lot….

  8. Amitavo
    May 3, 2012 at 4:06 am

    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

  9. John Mathai
    May 28, 2012 at 3:28 am

    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.

    • May 28, 2012 at 7:40 am

      Hi John, that would be the case if the governing Poisson equation reduced to the zero forcing function Laplace form, \nabla^2\phi=0. 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 n_e<n_i (and hence \rho \ne 0) 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 \rho term that gives you the forcing function that gives you the Gaussian-like distribution. Hope this helped.

  10. olatunji
    May 30, 2012 at 11:03 am

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

    • June 1, 2012 at 6:51 am

      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.

  11. John Mathai
    June 1, 2012 at 1:19 am

    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.

    • June 1, 2012 at 6:58 am

      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.

      • John Mathai
        June 5, 2012 at 10:17 pm

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

  12. John
    June 18, 2012 at 9:59 am

    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.

    • June 18, 2012 at 10:03 am

      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.

  13. John
    June 24, 2012 at 12:08 pm

    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?

  14. mahnaz
    July 3, 2012 at 11:30 am

    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

    • July 13, 2012 at 5:18 am

      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.

  15. john
    July 6, 2012 at 9:14 am

    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;

  16. Maks
    October 18, 2012 at 7:10 pm

    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?

    • Maks
      October 18, 2012 at 7:23 pm

      I’m sorry for dumb question..
      I found where was a problem.

  17. sudhan
    November 17, 2012 at 6:47 pm

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

    • December 8, 2012 at 12:03 pm

      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.

  18. WEI Wei
    November 27, 2012 at 8:44 am

    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.

    • December 8, 2012 at 11:04 am

      Thanks for the link!

  19. Sayan
    February 4, 2013 at 4:07 am

    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?

    • February 4, 2013 at 4:30 pm

      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.

  20. Sayan
    February 7, 2013 at 3:55 am

    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

  21. February 11, 2013 at 11:40 pm

    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

    • April 10, 2013 at 4:12 am

      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.

  22. March 4, 2013 at 10:31 pm

    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.

  23. Sonu Sen
    February 8, 2014 at 2:54 am

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

  24. Prakash
    June 11, 2014 at 12:44 am

    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!

    • June 12, 2014 at 4:20 am

      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.

      • Prakash
        June 12, 2014 at 5:05 am

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

  25. Helen
    July 18, 2014 at 5:42 am

    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!

    • July 18, 2014 at 7:35 pm

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

      • Helen
        July 21, 2014 at 2:16 pm

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

  26. farid
    July 19, 2014 at 3:22 am

    hi
    can i use this code for higher densities (about 10^20 #/m^3) ?
    is this code suitable for tokamak plasma simulation?

  27. john
    August 11, 2014 at 12:09 am

    Hi Lubos,
    I’m interested to particle in cell LINAC (linear accelerator) application where it used protons charges. Is Your matlab code adapt for this application? However I have a series coupling resonant cavity where the proton beam through it. For example initial energy is 7 MeV and the cavity were designed to accelerate untill 10 MeV.
    Thank

Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> In addition, you can use$latex ...$ to include equations.