Write PLEP on how to implement plasma simulation capabilities

[namurphy]

One of the most important things that I want PlasmaPy to be able to do is make a high level interface for performing simulations - as high level as we can possibly make it.

Right now, setting up plasma simulations is often very tedious and Fortranic. To install a package, it's often necessary to sign a user agreement, download and install software libraries, try to compile the code, figure out why libraries aren't linking correctly, modify the initial and boundary conditions in the source code, run the simulation, debug the things you changed with semi-helpful error messages, and then trying to run the simulation again. There is a huge amount of overhead that goes into setting up and analyzing simulations. Often if we want to switch numerical methods (e.g., by going from finite difference to finite element) we then have to use an entirely different code, which makes benchmarks much more difficult. Many of the codes also prohibit redistribution. The data formats for different packages are typically different, which makes it harder to perform the same analysis on different simulation sets.

We should also keep in mind that there are some cool codes like BOUT++ that are open source, and some cool projects like PICKSE that are making codes more widely available and easier to use. There are also existing general frameworks with Python interfaces, for example FEniCS for finite element simulations.

It would be really helpful to write a PLEP that describes how we will implement our simulation framework for multiple types of methods. I believe our first focus should be on making it as easy as possible to set up a plasma simulation and to switch between different numerical methods and systems of equations. It should be possible to set up and run a simulation in ~10-20 lines of code, or possibly even less for pre-defined cases (like the Rayleigh-Taylor instability, which I again propose that we rename as the Swirly Mushroom Instability). We would probably do really well with an object-oriented framework as well, and this would very much relate to our Plasma abstract base class.

We will also need considerable geometric flexibility to model, e.g., tokamaks and stellarators. Finite element methods generally work really well for those. Plus there are other models like gyrokinetics and multi-fluid models (e.g., ions, electrons, and neutrals), so we want the equations to be easily switched in and out too, and we want to be able to turn certain terms on and off. We'll probably need a Grad-Shafranov (equilibrium) solver too. We also want to be able to perform PIC simulations, and compare to experimental results.

I also really appreciate the Dedalus spectral code which allow the equations to be written in plain text. I also very much appreciate the numerical method of the HiFi spectral element framework which allows for the solution of equations that can be written in a pretty general flux-source form.

We have two broad directions that we can take:

1. Use existing codes and write open source tools that couple them together.
2. Develop new functionality, while leveraging Python packages like FEniCS, SymPy, and SciPy that are available through conda-forge.

The first approach is what is being taken by OMFIT, which unfortunately has a license that prevents redistribution. Plus, many of the codes that have been coupled into OMFIT are also closed source. The strongest benefit of this approach is that it leverages existing codes. My apprehension about this approach is that these codes were written using different languages and coding conventions, and were not originally intended to be intercompatible.

I personally think we should take the second approach and develop a new simulation framework and work to make sure that it is intercompatible from the beginning. I believe this approach will in the long run lead to much more maintainable, consistent, and readable code. We could implement numerical methods from different packages as well. The drawback of this approach is that it would be a huge amount of work to get it set up. Another reason for not concentrating on the first approach is to not duplicate the great work being done by the OMFIT group.

We have a lot to figure out for this and a lot of different factors to weigh, so more thoughts are definitely welcome!

Related issues: #13, #14, #16, PlasmaPy/PlasmaPy#33, PlasmaPy/PlasmaPy#46, PlasmaPy/PlasmaPy#82, PlasmaPy/PlasmaPy#167, PlasmaPy/PlasmaPy#365, PlasmaPy/PlasmaPy#457, PlasmaPy/PlasmaPy#459, PlasmaPy/PlasmaPy#69

[namurphy]

And there's another question that we really need to consider: Is Python the right language to build the core of our simulation framework in? It may be more Pythonic to use Julia.

Some considerations related to Python:

• It's a very mature language that is still actively under development.
• We can use Cython, Numba, and Pythran to get compiled speeds.
• There are frameworks like Dask that can help with parallelization.
• There are packages like FEniCS that we can build upon for certain types of simulations.
• Developers and users will not need to learn a completely new language.
• Dealing with extensions like Cython, etc. can cause headaches at times.
• We can use packages like SymPy for, e.g., automatic code generation and other things.
• There are parts of Python that are going to be slow.
• Cythonizing functions often necessitates using a profiler to check which parts are Cythonized and which are performance bottlenecks.

Some considerations related to Julia:

• The purpose of Julia is to do exactly this sort of thing: large-scale scientific computing.
• We won't need to deal with bridging Python and Cython.
• Julia is a relatively new language. It is still evolving so there might be API changes, and we might end up having to either develop new functionality that is not yet available or use wrappers for packages like FEniCS.
• Developers would need to learn a new language.
• We would help move scientific computing away from Fortran.
• Julia can be installed using conda-forge.
• We would only need to decide between three units libraries!

Some additional considerations:

• This may be a both/and rather than an either/or decision. We could make use of both Python and Julia.
• The core computational package will be more maintainable in the long run if it is in a single language.
• We can create wrappers so that we can call Julia functionality from Python, though this would involve two interfaces. With wrappers, we could potentially make it so that users do not have to use Julia at all.
• We can design the simulation output so as to be easily compatible with our Plasma abstract base class.
• I like saying "abstract base class". It sounds really cool. 🐱
• I don't know Julia yet, so there are likely things I'm getting wrong.
• Who knows...maybe the future of Python will involve a just-in-time compiler?
• Whatever the case, we probably would want to package the simulation capabilities in a different package than PlasmaPy, but we could call that package from PlasmaPy.
• We should also think about designing wrappers to call PlasmaPy functionality from Julia.

[lemmatum]

I would note that there are a number of packages in Python which have been ported to Julia and I think they're avoiding using wrappers... Might want to ask them about their approach(es) and why they did it that way.

My personal preference would be to stick with Python until Julia is at least v1.0

[namurphy]

I just noticed that Julia 1.0.0.rc1 was released very recently!

[namurphy]

Update: Julia 1.0.0 has been officially released.

[lemmatum]

Just saw that as well! So yea, I'm up for learning some more Julia if we want to go that route for the high performance parts of our code 😄

[lstagner]

Hey if you guys are interested I started the JuliaFusion organization a while back to share the Julia code I wrote for my thesis. There isn't much in there at the moment just an a small library for reading EFIT geqdsk file, an Equilibrium library for representing the fields, and a code to calculate guiding center orbits. I have plans to write a package for representing flux functions, a package for getting atomic cross sections (probably using the Aladdin database), a collisional radiative model, and a package for simulating neutral beams and spectra. Basically all the things I need to make a Julia FIDASIM which I also develop.