# PIC (particle-in-cell), space charge tracking simulation?

**URL:** https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077
**Category:** Numerics
**Tags:** package
**Created:** [July 6, 2019, 8:34pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077 "2019-07-06T20:34:03Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [July 6, 2019, 8:34pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/1 "2019-07-06T20:34:03Z")

</div>

There’s a lot of demand in plasma/accelerator science to simulate space charge through so-called particle-in-cell model (which is somewhat like finite element analysis in structural dynamics etc)

Many of them, for example [fbpic](https://github.com/fbpic/fbpic) uses Python for the apparent advantage for scientists users and for performance reason has to use `Numba`. Also, there are more legacy code written in Fortan (e.g [Using ASTRA at DESY-Hamburg](http://www.desy.de/~mpyflo/))

Is there any related project going on in Julia community? If people were to write something in Julia, would there be a real advantage over Numpy?

---

<div class="post-metadata">

### Author: ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)
#### Post date: [July 6, 2019, 9:38pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/2 "2019-07-06T21:38:00Z")

</div>

I am not an expert of Python/Numpy/numba… (I have developed // and HPC physics solvers in C++ for many years) so my answer is only a strong opinion.

1. Pros of Julia for HPC simulations

- Solve the problem of two languages (productivity of a dynamic language and performances of C++ or Fortran codes)).
- Very expressive, complete and rich (too rich ?) language that allows you to express accurately and efficiently an amazing variety of algorithms.

1. Cons.

- People form classical HPC (e.g. CFD) have not yet massively made the switch to Julia and so they are (IMHO) a bit late to react compared to other scientific domains (e.g. mathematical optimization, data science,…). As a consequence, the ecosystem may appear a bit sparse and you may have to develop more stuff by yourself. For example, I am presently coding parallel geometric multigrid solvers (and it is a real pleasure to do it in Julia).
- The multi-threading API is not yet totally completed.

With this in mind, I estimate that the benefits (generic, modular, composable algorithms, functional programming, meta-programming, performance tools…) to start new development in Julia largely exceed the drawbacks of the relative youth of the ecosystem: the productivity in Julia allows you to focus on algorithms and deliver more performant and maintainable simulation tools.

P.S. It takes some time to learn Julia but I think it is a good investment 😉

---

<div class="post-metadata">

### Author: ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)
#### Post date: [July 7, 2019, 9:52am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/3 "2019-07-07T09:52:23Z")

</div>

@LaurentPlagne You make a very good reply. It would be really interesting to hear more about the features of the Julia language help in your new work.

@jling I cannot help you directly. I have a background in physics, but not in accelerator physics. I would say though - if you are starting a new project from scratch use Julia. Don’t code any fundamental algorithms by yourself - if you ask on this forum you will probably discover that there is a Julia package which has that function.  
However if you need to use the Astra package from DESY you can use the Fortran functions in that package quite easily. SO you use Julia as a ‘wrapper’ around existing functions.  
From your post I do not think you are using the Astra package, rather using it as an example.

---

<div class="post-metadata">

### Author: ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [July 7, 2019, 1:42pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/4 "2019-07-07T13:42:21Z")

</div>

I am using ASTRA at the moment, the problem with ASTRA though, is that it’s design to run like:

```julia
$ ASTRA xxx.in

```

and spits out some files record the results of the simulation, There’s not really any `interface` can be used. And I think this is a common pattern, that performant software just meant to be ran like a binary executable if it’s impossible to abstract any _common_ core set out (for example what openblas does)

---

<div class="post-metadata">

### Author: ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)
#### Post date: [July 10, 2019, 11:40am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/5 "2019-07-10T11:40:27Z")

</div>

First, I apologize for my general feelings about the Julia suitability for HPC physics simulation: probably not what the OP expected.

I have already mentioned (in the Pros paragraph) the main features that makes Julia so convenient for my work. Functions that are generic by default (duck typing), abstract types for implementation swiches, nice and light syntax for parametric types, first class citizenship for functions, introspection, macros… Each time I have to switch back to C++, I feel like putting on an heavy dragon’s class armor.

I should also mention integrated features like @code\_native and @btime that accelerate the tuning of computationally intensive functions. Obviously there are very good tools (vtune, valgrind…) for C++ for profiling or assembly investigation but they are external/commercial tools. The Julia’s ecosystem is also becoming very competitive. For example Makie that allows for rather sophisticated and performant visualizations 😉  
 ![test](https://global.discourse-cdn.com/julialang/original/3X/8/d/8df9559ccbe6d26e7a33b7e8c6e79fa8b5ea64ac.gif)

Package system of modern languages like Julia (or Rust or Python) does not exist for C++ or Fortran.

Of course they are still some missing features/ perf issues like:

- allocating array views
- OpenCl (Vulkan ?) backends competitive to CUDA’s Julia handling
- Polished multi-threading API

but the Julia’s evolution is so fast that I am pretty confident for the future.

---

<div class="post-metadata">

### Author: ![laurence-feng](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurence-feng/32/11221_2.png) [@laurence-feng](https://discourse.julialang.org/u/laurence-feng)
#### Post date: [July 10, 2019, 12:03pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/6 "2019-07-10T12:03:33Z")

</div>

I had the same question too…

---

<div class="post-metadata">

### Author: ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [July 10, 2019, 3:09pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/7 "2019-07-10T15:09:56Z")

</div>

you made great point and to me they all sound like relevant things to PIC or any simulation really; but all these algorithms are quite complicated so I don’t know how to feel about technical details and how doing it in Julia would feel.

Me posting this was more like: google ‘julia particle in cell’ resulted nothing, let me put up something so people could discuss or even better, if a team is writing a PIC program and can make some comments.

---

<div class="post-metadata">

### Author: ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)
#### Post date: [July 11, 2019, 12:42am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/8 "2019-07-11T00:42:04Z")

</div>

Hey everyone, since the publication of Anthony Peratt’s 2015 book I have been very curious about all this. There is a guy named Donald Scott who made a [presentation](https://www.youtube.com/watch?v=t8tqgntbjyE) about Anthony Peratt’s particle-in-cell simulations of spiral galaxies using plasma discharge. I have taken some excerpts and made it into a song:

[![](https://global.discourse-cdn.com/julialang/original/3X/0/d/0d042d6eedebe71013ce87231d9ccae6297eca3a.jpeg "Escaped Electron (feat. JazzJet)") ](https://www.youtube.com/watch?v=DuUCEJPvebs)

if someone wants to do some spiral galactic PIC of a plasma discharge in Julia, do it while listening to this!

---

<div class="post-metadata">

### Author: ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)
#### Post date: [July 12, 2019, 5:45am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/9 "2019-07-12T05:45:39Z")

</div>

I’d be interested in participating in this project. IMO Julia will excel in speeding development of the peripheral components of a PIC code, since the core of these codes is a small fraction of all LOC.

There are umpteen different PIC codes written in a variety of languages. What should make the Julia version stand apart is it’s interoperability with FE grids rather than just Cartesian boxes. Otherwise, why recreate so much effort?

It should also allow the user to apply their own Ohm’s law for hybrid equations, and use the Darwin approach. Also implicit timestepping would set it apart from other free implementations.

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [July 12, 2019, 6:08am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/10 "2019-07-12T06:08:55Z")

</div>

Check this out [Programming the material point method in Julia - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S0965997816302769)

---

<div class="post-metadata">

### Author: ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)
#### Post date: [July 12, 2019, 7:45pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/11 "2019-07-12T19:45:42Z")

</div>

It’s so interesting that the same technique has so many different names! The particle shape function in plasmas needs to be higher order than that described in that paper, but it’s a simple extension.

It’s all quite a lot of work though, and I don’t have the bandwidth. Any other takers?

---

<div class="post-metadata">

### Author: ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [July 12, 2019, 9:55pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/12 "2019-07-12T21:55:52Z")

</div>

Problem is I don’t have experience writing performance-critial code like this, it would be nice if there’s an article talking about toy model / code (even in python or whatever)

---

<div class="post-metadata">

### Author: ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)
#### Post date: [July 12, 2019, 10:05pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/13 "2019-07-12T22:05:30Z")

</div>

This particle / grid approach is exactly what PIC is about.

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [July 13, 2019, 5:13am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/14 "2019-07-13T05:13:11Z")

</div>

I know, that’s why I posted the link. 🙂

---

<div class="post-metadata">

### Author: ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)
#### Post date: [July 13, 2019, 8:33am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/15 "2019-07-13T08:33:56Z")

</div>

This paper has some nice pictures [https://arxiv.org/pdf/1310.7866.pdf](https://arxiv.org/pdf/1310.7866.pdf), and it’s got Hartmut Ruhl and Amitava Bhattacharjee on the author list. They know their stuff. It covers the basics of particle shapes and mapping quantities from the grid to particle and vice-versa.

All the gory details are in the Particle Simulation Code manual - I can’t seem to google a copy of it.

On the Maxwell solver side there is [https://github.com/wsshin/MaxwellFDM.jl](https://github.com/wsshin/MaxwellFDM.jl), which appears to be under development. I know that the Julia FEM package is crying out for developers, and really feel as though that should be the ultimate goal.

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [July 13, 2019, 10:40am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/16 "2019-07-13T10:40:31Z")

</div>

> [@jcook](#):
>
> under development. I know that the Julia FEM package is crying out for developers, and really feel as though that should be the ultimate goal.

I’m not sure what you mean here. FDTD is not really finite elements, and PIC isn’t finite elements either.

---

<div class="post-metadata">

### Author: ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)
#### Post date: [July 13, 2019, 11:30am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/17 "2019-07-13T11:30:41Z")

</div>

I should be clearer. In PIC moments of the distribution function gathered from particles are projected onto any grid. Then any grid will do. FDTD is the easiest to implement but there shouldn’t be any restriction to generalise to FEM, except some extra coding trickiness for dealing with shape functions and element shapes.

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [July 13, 2019, 11:37am UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/18 "2019-07-13T11:37:33Z")

</div>

I believe that FDTD relies on the “nested” character of the Cartesian grid.  
It gives it the conservation properties that make it a success.

I think you’re right, PIC should be able to work with any grid at all. I’m not sure though that it will necessarily lead to a more efficient code. Part of the appeal of PIC is that it can deal with “fronts” moving through space. In that case static graded grids don’t necessarily help very much. And if the graded grid needs to be redone many times over the simulation one could just as well work with a Lagrangean method without the complication of advection.

---

<div class="post-metadata">

### Author: ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)
#### Post date: [July 13, 2019, 1:08pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/19 "2019-07-13T13:08:32Z")

</div>

Could you explain what you mean by nested and graded grids? Do you mean adaptive mesh refinement?

Regarding conservation, explicit Cartesian FDTD plasma PIC codes can be formulated to exactly conserve either momentum or energy but not both. I’m most familiar with the momentum conserving kind. With FEM formulations there are more ways of getting it wrong, but it would still work.

A PIC code that works with arbitrary grids would open up opportunities for gyro kinetic particles for which Cartesian geometries fall short.

Are you also refering to the level set method for front tracking? I’ve never worked with that so would be interested to hear your thoughts.

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [July 13, 2019, 1:14pm UTC](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077/20 "2019-07-13T13:14:49Z")

</div>

I am possibly making an unjustified assumption that this is a topic of general interest, as it is only tangentially related to Julia. Perhaps this should be discussed through private messages?

The FDTD method relies on the “staggered” (or nested) grids for the dual quantities. That works well when the grid is Cartesian.

I have to state for the record that my familiarity with PIC codes is for the solid mechanics equations of fast dynamic deformation. The fronts refers then to shock waves or shear bands. I’m not really familiar with plasma physics or the application of these codes plasma physics. Sorry.

[Next page](https://discourse.julialang.org/t/pic-particle-in-cell-space-charge-tracking-simulation/26077.md?page=2)
