Include a full particle in cell code

[StanczakDominik]

With the particle stepper pretty much being finished and live I believe we could begin thinking about including a full pic simulation within the code.

I don't think we need to focus on beating commercial and research codes at scalability, efficiency or speed, it'd be nice if that turned out to be the case, but it'd take a lot of effort and it's probably better to start small.

I think what would work better within open source is to have a simulation code that's feature complete, reliable, well documented, with ample test cases, visualization, something people could study and learn from, something you can compare results with. I know that part has been a pain in my experience.

It'd also be somewhat neat to have the same interface as the MHD simulation.

Generality is also nice, being able to easily replace and compare the numerical methods used seems like a good idea.

Of course, that doesn't mean we're going slow on it either, I'd like to start with numpy, numba and possibly add pycuda later on (that'd need a separate issue).

I've got a code of mine that I'd like to adapt here.

What does everyone think about this? I'd very much appreciate input and help with this one. Literature suggestions beyond Birdsall and Landon would also be nice!

[antoinetavant]

I'd like to contribute. I totally agree with the objectives you said, especially about speed and well documented and tested.

However, I think we have to discuss about few points to begin with :

1. For what simulation? Space plasma (collision less, no borders, maybe relativistic), laser-plasma interaction (electro-magnetic ), cold plasma (collisionnal, with wall interaction) ? We could try to make for all of them, but that will be a pain.
2. To which complexity ? 1D ? 2D ? how many particles ? and chemistry ? I think the simpler the better.
3. Should we provide a complete PIC function, with user-friendly interface and visualization, or a toolbox of modules (particle stepper, collisions, Poisson solver, ...) with a simple framework, that users would assemble to do the simulation they want ?

I would answer the 3. by the second proposition : let's provide a toolbox and a framework. So that every users can simulate what they want, the way they want, and we won't have to think of all of the possibilities ourselves. And of course we would provide examples. I guess this make reference to #73 .

I'm not sure I'm clear. Anyway I'm supper down for that !

PS : there is some Python PIC codes https://github.com/fbpic/fbpic https://github.com/skuschel/postpic I've got plenty of paper about PIC validation and benshmarks, but Birdsall is still my bible ! :)

[StanczakDominik]

@antoinelpp re: your 1, I know that's a pain, but that's kind of a given for a generic package for plasma physics, probably one of the broadest fields in physics.

I think the framework idea is best: letting people pick and choose the kind of simulation they want to run helps the most people, though it is indeed going to be difficult. Besides, if they built something on top of it, that can always be pulled back into the main package.

If I understand it right, this would mean making a few presets: Cold plasma, Laser-plasma, Space... maybe picking one that's relatively simple and inviting collaborators for other kinds of plasmas?

Do you think it's too crazy to want to, at least eventually, support full 3D PIC simulation in PlasmaPy, or is that potentially realistic?

[antoinetavant]

@StanczakDominik 3D pic is very time consuming. Even in the HPC sector, people don't try it too much. However, we can think of it, starting from low but keeping in mind that 3D can eventually supported.

Building a generic framework for the different presets sound nice ! I like that.

[StanczakDominik]

Awesome, let's start with 1D and 2D first then :) I made a project for this one. I'll try to make a bunch of issues for things to do to make this stepper a full fledged PIC code, but I probably won't be able to do that today. If you want to try, go right ahead!

[lemmatum]

Is anyone working on this? Seems like we don't even have a PR for the PIC stuff. I'd be interested in contributing, especially as I've been meaning to learn about PIC codes for my Langmuir probe experiments.

[StanczakDominik]

Remember my slow thesis code? That was what I was gonna contribute for this, and then it came out slow as heck and I got demotivated.

I can into Cython now, though, so I may have to revisit that!

[antoinetavant]

Hi there, I tried a few months ago to contribute on this issue and I arrived to the conclusion that :

1. I wasn't yet sur of the structure of the pic code that we would need
2. I had to build one outside of PlasmaPy So I now have a 1D Pic code working nicely, I can show it to you and I think I can try to do a PR with it.

However it is not perfect yet, especially question coverage and versatility

There is it. it is not user-friendly, I'll quickly update the Readme! :grinning:

[StanczakDominik]

@antoinelpp so we essentially have two 1D codes. Wanna compare and combine? :D

[lemmatum]

@antoinelpp excellent, thanks for linking it! You might be able to just start a PR as is and then use our feedback to help with coverage and versatility (if you need it).

By the way, do we want a module for electromagnetic equations and such? I assume we will need things like Maxwell's equations and Lorentz force for a much of our simulation needs.

[antoinetavant]

@lemmatum Allright I'll do that in a few ! :smile: Concerning improvements the Lorentz force is quite easy to take into account, but for the electro-magnetic equation I really don't know how to do it...

@StanczakDominik I took a look at yours, and it is really good ! Definitely we should combine !

[lemmatum]

Perhaps some of these resources could be of use? I think the two PDFs mention EM fields.

• http://plasma.kulgun.net/simulation/
• http://cips.colorado.edu/mathcad/13_PIC_Beam_Plasma_5.pdf
• http://hedpschool.lle.rochester.edu/2011SummerSchool/lectures/Ren.pdf
• https://www.particleincell.com/blog/

[tulasinandan]

Guys, I am back after a long hiauts. I'll join forces. We also have https://github.com/UCLA-Plasma-Simulation-Group/PIC-skeleton-codes to get some guidance from. These are written in Fortran but have python interfaces. The coding style is not my preferred style but the codes are very well tested and should be good for benchmarking.

[antoinetavant]

Hi @tulasinandan thanks for reviving the PIC part ! :smile: I didn't know about the UCLA-Plasma-simulation Group, i'll take a close look at it ! Do you think of using C/Fortran libraries for our package ? And what do you guys think of parallel code ?

[tulasinandan]

I wish I oculd have put in more effort earlier but job applications and proposals are time consuming projects by themselves!

I would vote against using C/Fortran libraries, for portability issues. Although we could create small interfaces for codes like UCLA codes. This should be easy as they already have example python codes. All we need to do it to write an interface that lets PlasmaPy communicate with their fortran libraries. That should be part of the simulations interface class.

As far as parallelism is concerned, I would vote against that as well in the beginning with a single caveat. We should create our data structures such that we do not have to modify them when we decide to parallelize things. The reason why I am against parallelism at the first cut is because particle exchange and the related memory management is a very tricky business which would take up a lot of time.

https://picsar.net/features/particle-sorting/

This one actually does heavy optimization but even with minimal optimization, adding/removing particles in a processor's memory needs a lot of thought about dynamic memory management.

[tulasinandan]

Another code we might want to keep checking https://github.com/fbpic/fbpic

I don't know how flexible the code is to be used as a framework but they have a lot of things implemented already.

[lemmatum]

Easy solution: make the data structure all arrays where the elements don't depend on each other, so that you can just map functions over them in an embarrassingly parallel way (in Cython).. can't have race conditions if there if there are no clashes in memory :stuck_out_tongue_closed_eyes:

More serious solution: Write as much of the code in a purely functional way as possible, and then those functions which are "impure" will be where you have to be careful about optimization and memory management. (Have any of you programmed in Haskell before?)

[tulasinandan]

@lemmatum I have not programmed in Haskell but I usually like to follow the fundamental principle of functional programming, that functions should be as independent of everything else as possible.

Now, whatever we do, I am not sure how we would take care of the memory management issues for particle arrays in an easy way (for parallel applications). For serial, it's trivial. You assign memory and you're done as long as you work in the classical limit and don't have particle being created or destroyed.

In parallel PIC simulations, the number of particles on a given processor is never fixed in time. At this instance processor 23 might have 10234 particles, and a in a few time steps 439 particles move outside its domain (and are all from different parts of its memory), and maybe 675 particles move into its domain. These can be easily added at the end. However, this means dynamically sorting and reallocating memory positions for particles.

Another point: If we do not want to dynamically change total memory allocated for particles, we would have to assign a big chunk with a buffer initially, and then hope that number of incoming particles is never more than available memory. This works in a number of nearly incompressible situations. However, some times in very compressive cases, one does hit the wall and the simulation has to crash.

[lemmatum]

@tulasinandan good to hear! I know too many physicists that build code with tons of side-effects which then makes any sort of optimization impossible. :frowning_face:

On memory management: Presumably each particle will occupy the same amount of memory. Could we write a memory allocator that:

1. deletes particles which have left the domain
2. keeps a cache of address spaces for particles which are still in the domain and ones which have been deleted
3. chooses the earliest spot in the cache of empty particles

One way to do this may be to just keep two arrays, one of existing particles, and one of deleted particles. Write a class with methods that allow you to transfer the "pointer" to a particle's memory from one cache to the other. And for added safety, fill all attributes of deleted particles with None or the like (that way you can still keep the particle model even in the deleted list). Then if you want to get fancy with parallelization, you can make one or two more cache arrays that keep track of which existing particles and deleted particles are "locked" (in use by a processor) at any point in time. There may be some asynchronous libraries we could use to this end.

[tulasinandan]

@lemmatum Interesting idea! I think we could improve this a bit more. However given my primitive programming style I am not sure how this would be implemented in practice. I'll rephrase your version to make sure that I understand it (and maybe add a bit more?):

Is it possible to have an array of pointers only? If yes, then for parallel implementation, one could allocate 2*n pointer arrays (n=parallelization-dimensions, 2n to have one for each neighboring processor) to send particles and a temporary buffer to recieve particles. The particles in this temporary buffer would be moved to the sent particle slots after the "send" is completed via the pointer arrays. This way we would simply be using one smaller scratch array for recieving particles.

I am not sure if something like this is implementable in python. Also, all this is completely redundant in a serial code.

[lemmatum]

@tulasinandan Python has the global interpreter lock (GIL) which prevents true parallelization. To get low-level parallelization and memory management we would have to implement this in Cython (which is pretty close to Python - it's like Python + static types + some different defaults for memory allocation to numerical types). Looks like Cython has a better abstraction than just having an array of pointers. We should dig a bit and see what sort of abstractions Cython has for locking memory while parallelizing - the model you suggests sounds general enough that someone might've already implemented it under the hood in Cython.

As for serial code, we could have a flag or something that just totally omits this memory management and data locking mechanism.

[tulasinandan]

Oh I so wish that I had learned C at some point in my life. I tried a couple of times but never pushed any project hard enough to learn and use C on a regular basis. I guess this will eventually be my first dig at something that is close to C. At the first dig we should probably stick to the serial versions that @StanczakDominik and @antoinelpp have and increment on top of that base.

[lemmatum]

@StanczakDominik you've read the Cython book, do you remember anything about memory allocation for parallelization? It may just be that we need to slice up the PIC simulation into subregions and then let each core handle a region.

@tulasinandan Here's some stuff on the GIL. Looks like Numba is another way of dealing with GIL/parallelization.

[tulasinandan]

It may just be that we need to slice up the PIC simulation into subregions and then let each core handle a region.

@lemmatum Yes, that is the standard strategy: slice up the domain into subregions handled by individual cores. For fields it's pretty straightforward, the cores simply exchange boundaries with their neighbors. For particles we have to think about all the mess mentioned above.

[StanczakDominik]

As for the Cython book, that only skimmed on parallelism. It touched on disabling the GIL, showed a couple of examples and mentioned OpenMP, if I remember correctly, but there wasn't much about the subject. I'll keep digging into this...

I'm loving the discussion in this thread, by the way!

[lemmatum]

@tulasinandan I see! I didn't initially understand that this was about dealing with particles being passed off between subregions where each subregion is assigned to a processor/thread.

@StanczakDominik yea, I remember parallelism being a bit scant. I guess we'll have to dig into OpenMP and Cython's online documentation.

Also Cython has a structure called Memoryviews which may be applicable. Here is some stuff on Cython memory allocation as well.