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.