Discussion of PlasmaPy and NIMROD

Nick (as I understand that you are working on the design most closely related to this) and all,

I’ve been developing python 3 analysis code for NIMROD that I think it might be useful to review while considering the design of a plasmapy implementation. I’ve read the PLEPs and reviewed the plasmapy source and docs and I think I have some basic understanding of your approach (but not comprehensive). My code is far from perfect and I don’t mean for it to be used as the basis of any implementation but it is reasonably mature as a post-processing framework for MHD runs. We use it to compute flux-surface averages (annihilation operator for B.grad()), compute the trapped particle fraction, track particle drift orbits, a beam emission spectroscopy synthetic diagnostic, and plot general equations up to second order. It might be useful to discuss the design choices that I made (again, not perfect) and things that I think are good versus parts that I dislike.

Ultimately, I would like to abandon a NIMROD-centric postprocessing suite and integrate with plasmapy if possible. I like your test-driven development, sustainable code ethos and I think using a community solution would be good for all in the long term.

Let me know if you are interested.

Yes, we should definitely talk more about this! Right now we are in the architectural design phase for plasma simulation software. Within the next month or two, I want to get back to working on the PlasmaPy Enhancement Proposal (PLEP) on plasma simulation, so it’s a good time to have big picture discussions on how we are going to go about this. It would be good to have several of us work on the PLEP together in order to get a lot of different perspectives. I’d be interested to talk more with you about where the functionality that you described should live in the software ecosystem, and how to make it abstract enough that it can work with the outputs of a variety of simulation packages. This would be a good topic for future community meetings. Thanks for bringing this up!

I tagged this as “simulation” but perhaps a more appropriate tag would be analysis. It seems like you have a wide-ranging vision for not only pre/post processing but also creating PIC/fluid codes, right? This would be a discussion on the former.

Reviewing the plasmapy code my thinking is that I could create a new class that inherits from GenericPlasma (or maybe Plasma3D/PlasmaBlob). I have a number of questions on these classes: is it intended that the base classes will be extended as needed? I suspect that a way of querying the derivatives of fields will be needed (e.g. to track drift trajectories). Also, while the ability to handle multiple ionization states is included with the average concept, true multiple species seems to be precluded by the ion_temperature function (or the ability to track different ionization states with separate fluid equations).

I also have other ideas associated with this that I’d rather discuss to gauge interest first before writing anything detailed.

I would be happy to talk about this at a community meeting. Does Jitsi allow for screen sharing? That would be useful. I can also send out or post the code but I am concerned that without context it wouldn’t be that useful.

Definitely! They’re closer to drafts right now :smile: derivatives would probably be handy; I’m struggling to think of a generic solution to provide those, given the various kinds of coordinates we could have there, though… given Plasma3D as it stands consists of rectangular coordinates, perhaps we could have a diff method in GenericPlasma for which every subclass (e.g. kind of coordinates) would have to provide an implementation?

This got me thinking of xarray, and I checked how they’re handling it:

This feature is limited to simple cartesian geometry, i.e. coord must be one dimensional.

I guess it really isn’t a trivial problem…

I’d love to discuss those! Nick is supposed to be away for this Tuesday’s meeting, though, so if you’d like his input it may make sense to wait for the next one. But we can preliminarily discuss it with whoever shows up :slight_smile:

Sounds great. We can discuss tomorrow informally.

Here’s a written version of some talking points for today.

I’d like to advocate for the concept a of “view” as a separate base class from BasePlasma with the idea that instantiations of BasePlasma store the full data in an N-dimensional domain and can evaluate fields (and derivatives) at arbitrary points. A “view” is a concept that is useful for post-processing. I’ll call a proposed base class for a view BasePlasmaView to get started. An instantiation of BasePlasmaView could range from a point evaluation up an arbitrary N-dimensional grid and would also be specific to tensor rank (e.g. separate implementations for scalar, vector, and rank-2 tensor). The BasePlasmaView should contains standard mathematical operations such as curl, grad, dot and cross to enable user defined operations for post-processing.

It might be easier to consider what this interface would look like. For example consider the following as a proposed interface to plot the thermal particle drifts in 2D:

# I'm importing the base classes here but derived class would be needed
from plasma_base import *
from plasma_base_view import *

# best to get these from the formulary but let's not worry about that now
echrg=1.609e-19 # C
md=3.3435860e-27 # kg
kb=1.609e-19 # J/eV sorry not sorry

# plasma_base instance (consider this a 3D domain in this example)
simulation_data = BasePlasma()

# plasma_base_view function to generate a 25x25 2D grid defined by 3 points
grid = grid_2d_gen([0.9, -1.5, 0], [0.9, 1.5, 0], [2.5, -1.5, 0], 25, 25)

# get view of B on this grid with 1st order derivatives and T without derivatives
B = simulation_data.magnetic_field(grid, dord=1)
T = simulation_data.ion_temperature(grid, dord=0)

# compute the drifts (operations on plasmaBaseView are overloaded and return another PlasmaBaseView)
vti = sqrt(kb*T/md) 
Bmag = B.magnitude(dord=0) # returns a scalar from the vector B
omega = echrg/md*Bmag
bhat = B.hat()
grdb = (vti*vti/(2.0*omega*Bmag))*bhat.cross(grad(B.mag()))
curv = (vti*vti/omega)*bhat.cross(curl(bhat).cross(bhat))

# hypothetically we could also let plasmaBaseView plot
grdb.plot_quiver(title=r'v_{Ti}^2/2\Omega_{ci}|B| \mathbf{B}\times\nabla |B|',
                           units=r'm/s')
curv.plot_quiver(title=r'v_{Ti}^2/\Omega_{ci}|B| \mathbf{B}\times\mathbf{\kappa}',
                            units=r'm/s')

Here’s some example plots from tokamak data (these are from my APS-DPP poster which no one is obligated to look at but I’ll link anyway http://nucleus.txcorp.com/~jking/JKingDPP19.pdf):


If grid points in a view fall outside of the domain of a BasePlasma, field values values should be set to NaN as matplotlib does not plot a color for NaNs. Note that although the view is a 2D rectangle, points outside the computational domain are not present in the plot. There’s also some extra contours on this plot for the computational boundary, the last-closed flux surface and an approximate scrape-off layer boundary that aren’t present in the above pseudo-code.

What I like about this is that it lets the user define expressions with a syntax that mirrors how we would write this on paper. I first encountered this while using Firedrake (https://www.firedrakeproject.org/) I think it is a great idea.

Something implicit in the code is the optional variable dord (which should be differential_order for consistency with the plasmapy style). This sets the differential order of a field and can be reduced on derived quantities for optimization if we know that propagation of derivatives isn’t needed. If I didn’t care about optimization, I could just ignore it except in the magnetic_field call and propagate all derivatives. Without being explicitly specified each derived field would propagate as many derivatives as possible by applying the chain rule.

As another example we can plot the force balance in 1D along a radial cord. Typically tokamak folks assume that this is just JxB = grad§ but of course there are other terms. The sonic flow in edge might contribute through the advective force so let’s check.

from plasma_base import *
from plasma_base_view import *

# plasma_base instance (consider this a 3D domain in this example)
simulation_data = BasePlasma()

# plasma_base_view function to generate a 1000 point 1D grid defined by 2 points
grid = grid_1d_gen([1.8, -0.035, 0], [2.15, 0.74, 0], 1000)

B = simulation_data.magnetic_field(grid)
p = simulation_data.pressure(grid, dord=1)
nd = simulation_data.ion_density(grid)
v = simulation_data.center_of_mass_flow(grid, dord=1)

jxb=curl(B).cross(B)/mu0
grp=p.grad()
rhovdgrdv = (md*nd)*v.dot(grad(v))
force=jxb-grp
force.plot_lines(flabel=r'Force',
                f2=jxb.data, f2label=r'\mathbf{J}\times\mathbf{B}', 
                f3=grp.data, f3label=r'\nabla p', 
                f4=rhovdgrdv.data, f4label=r'v\cdot \nabla v',
                comp='perp', style='varied',
                legend_loc='lower left', ylabel=r'Force density N/m^3')

Here’s the result (JxB ~ grad( p )):


We can zoom in to see that the force from advection is on the order of the residual from our solve to satisfy the Grad-Shafranov equation:

I’ve run out of time before the meeting but there’s one final point: Why do we need point evaluations to be contained in BasePlasmaView? These enable things where a point quantity might be a derived quantity involving derivatives at points that are not know ahead of time like flux-surface averages , drift-orbit tracking and synthetic diagnostics that involve wave equations such as Doppler back scattering.

1 Like

There’s a lot to think about here, but for now I’ll just link the minutes from today’s telecon: https://docs.google.com/document/d/1XGcY3AjF8e95osuQt6_BEhTOjr1Z_PReRPYyn9BW49A/edit#heading=h.ke07y03icy21

I’ll get back to this once I’ve had some time to gather thoughts :slight_smile:

There was a request during the call today to post some code of a implementations. I don’t see how to do an attachment so I’ll just post a link.

Linked (nucleus.txcorp.com/~jking/field_class.py) is an example of an implementation that I used for the above plots. It uses numpy arrays and has a number of other NIMROD-specific things that wouldn’t appear in a plasmapy version.

Alternatively we discussed use of sympy which might be cleaner.