PlasmaPy / PlasmaPy-PLEPs

A repository for PlasmaPy Enhancement Proposals
http://docs.plasmapy.org/
Other
8 stars 6 forks source link

Write PLEP on how to implement plasma simulation capabilities #17

Open namurphy opened 6 years ago

namurphy commented 6 years ago

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, https://github.com/PlasmaPy/PlasmaPy/pull/33, https://github.com/PlasmaPy/PlasmaPy/issues/46, https://github.com/PlasmaPy/PlasmaPy/issues/82, https://github.com/PlasmaPy/PlasmaPy/issues/167, https://github.com/PlasmaPy/PlasmaPy/issues/365, https://github.com/PlasmaPy/PlasmaPy/issues/457, https://github.com/PlasmaPy/PlasmaPy/pull/459, https://github.com/PlasmaPy/PlasmaPy/issues/69

namurphy commented 6 years ago

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:

Some considerations related to Julia:

Some additional considerations:

lemmatum commented 6 years ago

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 commented 6 years ago

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

namurphy commented 6 years ago

Update: Julia 1.0.0 has been officially released.

lemmatum commented 6 years ago

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 :smile:

lstagner commented 6 years ago

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.