The keep_yboundaries argument to open_boutdataset is implemented in the _trim function and only removes the y-boundary points at the edge of the individual files loaded. If the data has been squashed into a single file that includes a second (upper) target, those points will not be removed.

Currently open_boutdataset() will trim off all ghost cells, even at the outside edges of the domain. There should be an option to retain the guard cells at the edge of the domain.

This is not completely trivial because currently the trimming is done by a preprocess argument to open_boutdataset(), so the same function is applied to every dataset. In order to know when to leave the guard cells untouched, the preprocess function has to know the position of the dataset it is currently operating on within the grid. Right now this isn't possible, but there are two ways it could be done:

1. Use the filename of the dataset to determine it's position, which would require a change upstream in xarray as discussed here .

2. Get BOUT to write out the global position indexes to the dump files, which can then be used to work out the dataset's position in the grid.

Another option is to not use the preprocess function at all, and instead concatenate all the untrimmed data then cut out the unwanted chunks from the result. I don't think this is the best solution however because then you aren't doing the trimming in parallel, it's not as clean conceptually, and it could break a lot of the tests.

I think the best solution would be to use global indexes once BOUT actually writes them out.

To test the operation of the collect function in #59 I am using the create_bout_ds_list function from test_load, which was used to generate test data for test_animate.

The problem lies with the keep_xboundaries and keep_yboundaries kwargs of _auto_open_mfboutdataset. If True, these should increase the x and y dimensions, respectively, but it has the opposite effect on the data from create_bout_ds_list.

I have tried it with actual BOUT++ data, and the parameters behave as expected. It seems to be specific to create_bout_ds_list data, but I can't ascertain why this is.

Sometimes BOUT crashes whilst writing out the dump netcdf files. If this happens it can cause opaque errors, but we should be able to catch at least some of those and supply an additional error message suggesting that this might be the problem.

We probably want to use a try... except... statement around open_mfdataset(), and catch any MergeError, NetcdfError, or any ValueErrors matching certain regexes (depending on the type of error thrown inside xarray.concat()).

When restarting a simulation, BOUT includes the initial timestep in the output, so when you open from multiple restarts (e.g. open_boutdataset('./run*/BOUT.dmp.*.nc')) the returned dataset has duplicated timesteps. These should be automatically dropped.

Should we setup Travis-CI so that we get automated testing?

xBOUT should contain a new collect() function whose purpose is to maintain backwards compatibility with the old collect().

These would make it easier to phase out usage of the old collect.

It should have the same signature, so needs to return a numpy array. Would look something like

def collect(var, path, prefix='BOUT.dmp', info=True):
    datapath = path + prefix + '*.nc' 
    ds = open_boutdataset(datapath=datapath, info=info)
    return ds[var].values

There is no way that it would be able to print out the same info, but it might be able to throw some of the same errors.

Combining a large number of dump files (one per processor) into a single file is a common, useful operation. Working with a single file speeds up xBOUT operation considerably, and by dropping guard cells a single file minimises network traffic when downloading data.

This operation is simple to implement with xBOUT (although re-loading the resulting file still needs work #18), but is much slower than using the squashoutput function provided by boutdata. I have timed 4m23 for a 2.1GB dataset with squashoutput, while opening with xBOUT and calling ds.bout.save() took 21m38 (interactive jobs on a Marconi login node, but my experience from trying a couple of times before is that this is typical).

Don't know what the cause of this is - some profiling would probably be useful. Presumably any fix would be in xarray rather than here in xBOUT, unless there are some magic options to xarray.open_mfdataset() or xarray.Dataset.to_netcdf() that would speed things up.

We want to fully integrate a class structure storing all input file options.

We can't just use boutdata.BoutOptionsFile though, because it uses boututils.DataFile in places, so it will need at least a partial rewrite to work with xBOUT.

I actually tried to see if there was a way to rewrite BoutOptionsFile completely (aiming for a much shorter, less home-grown implementation), but the unique structure of BOUT.inp files foiled me. BOUT.inp files are similar to .ini files and YAML files, but have distinguishing properties which make handling them with variations of standard python library config classes difficult:

1. Missing initial section header. This is annoying but solvable - you can still read and write .inp files with configparser.ConfigParser but you have to special-case the first line.

2. Nested sections. This rules out easily subclassing configparser.ConfigParser - I thought maybe you could override ConfigParser.sections() with a recursive lookup method but as ConfigParser inherits its structure from collections.MutableMapping I couldn't see any easy way to make this work.

3. Specifying nesting in header names. The [section:subsection] syntax isn't shared by any commonly-used config file structure I could find. I think this rules out using configobj.ConfigObj, which is a shame because that one supports nesting (and nested substitution).

While trying out these things I did write some tests though, so we could use them (boutdata.BoutOptionsFile currently has no direct tests).

Therefore I suggest what we should do is:

1. Write some tests defining the expected behaviour of the options file class (I have some already, which I will submit later),

2. Copy over the basic parts of the boutdata.BoutOptionsFile code to xBOUT,

3. Reimplement the methods which rely on boututils.DataFile.

This should hopefully allow us to improve the boutdata.BoutOptionsFile code if needed but without breaking backwards compatibility.

Should the examples be saved as python scripts or as Jupyter notebooks?

Presumably we should change the licence on this repo (from Apache) to the same GNU one that BOUT-dev has?

There are many types of analysis functions which BOUT users might want to apply to their data, but which are fairly specific to plasma physics. Examples could include:

• Mode growth rates
• Flux-surface averages
• Conditional averaging
• Calculating ExB velocity
• Probability Distribution Functions of fluctuations
• Bicoherence
• Fourier analysis (using xrft) and spectrograms
• Spectral energy transfer diagrams
• Finding blobs? (maybe too specific to Storm)

These should live in xBOUT, either as functions or as accessor methods.

Should also find out if there is any crossover with PlasmaPy here.

If units are made mandatory attrs in BOUT, then these functions could be written to adjust the units of the data they operate on (best to do this with a decorator).

See the description at boutproject/SD1D#3

The following example doesn't work with matplotlib 3.1.1 and never. It works with matplotlib 3.0.1

array = np.random.randint(0, 100, size=(3, 4, 5))
xarray = xr.DataArray(array, dims=('t','x','y'), coords={'t': [10, 20,30], 'x': [1, 2,3,4],'y': [100, 200,300,400,500]})
dask_da = xarray.chunk(chunks={'t':3})
dask_da.bout.animate2D( animate_over='t', x='x', y='y')

my other package versions are:
animatplot 0.5.0.dev1
matplotlib 3.1.1
xarray 0.12.3+25.g86799b78
xBOUT 0.1.0
xstorm 0.1

The error message is the following:

Traceback (most recent call last):
File "/home/gregor/anaconda3/lib/python3.7/site-packages/matplotlib/animation.py", line 230, in saving
yield self
File "/home/gregor/anaconda3/lib/python3.7/site-packages/matplotlib/animation.py", line 1156, in save
writer.grab_frame(**savefig_kwargs)
File "/home/gregor/anaconda3/lib/python3.7/site-packages/matplotlib/animation.py", line 572, in grab_frame
buf.getvalue()))
File "/home/gregor/anaconda3/lib/python3.7/site-packages/PIL/Image.py", line 2542, in frombytes
im.frombytes(data, decoder_name, args)
File "/home/gregor/anaconda3/lib/python3.7/site-packages/PIL/Image.py", line 829, in frombytes
raise ValueError("not enough image data")
ValueError: not enough image data

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "test_blob_label.py", line 32, in
dask_da.bout.animate2D( animate_over='t', x='x', y='y')
File "/home/gregor/Documents/xBOUT/xbout/boutdataarray.py", line 83, in animate2D
save_as=save_as, ax=ax, **kwargs)
File "/home/gregor/Documents/xBOUT/xbout/plotting/animate.py", line 96, in animate_imshow
anim.save(save_as + '.gif', writer=PillowWriter(fps=fps))
File "/home/gregor/Documents/animatplot/animatplot/animation.py", line 187, in save
self.animation.save(*args, **kwargs)
File "/home/gregor/anaconda3/lib/python3.7/site-packages/matplotlib/animation.py", line 1156, in save
writer.grab_frame(**savefig_kwargs)
File "/home/gregor/anaconda3/lib/python3.7/contextlib.py", line 130, in exit
self.gen.throw(type, value, traceback)
File "/home/gregor/anaconda3/lib/python3.7/site-packages/matplotlib/animation.py", line 232, in saving
self.finish()
File "/home/gregor/anaconda3/lib/python3.7/site-packages/matplotlib/animation.py", line 575, in finish
self._frames[0].save(
IndexError: list index out of range

The few functions and methods which have been defined should have documentation. These could be auto-generated from their docstrings using something like Sphinx.

Can this be integrated with BOUT's documentation?

Once BOUT enforces that normalisations should be saved in the attrs of each variable, then xBOUT can automatically read them and (optionally) apply them to get back to physical units.

Should read through the CF conventions to see if there is a standard for this which we should follow.

Need to fix the tests which check that the old version of collect gives the same results.

Once I have a backwards compatible API function then this should be even more explicit.

xarray is deprecating the Dataset.drop() method in favour of Dataset.drop_vars(). However, Dataset.drop_vars() was only introduced in xarray-0.14.0, so if we change it now we have to increase the minimum supported version of xarray. Suggest leaving as-is for now, until our oldest supported version of xarray has to increase to 0.14.0 for some other reason, or until the newest xarray does not support Dataset.drop().

While trying to find the docs I realised there is currently no link to them on the readme. Are they okay to be added @ZedThree ?

The homepage of the docs is also out of date with the requirements - it should say that xarray v0.12.2 or above is required.

If you try to open only a single dump file with open_boutdataset():

ds = open_boutdataset('BOUT.dmp.0.nc')

then it will currently throw an error:

"Each run directory does not contain an equal number
of output files. If the parallelization scheme of 
your simulation changed partway-through, then please
load each directory separately and concatenate them
along the time dimension with xarray.concat()."

This is because it reads the processor parallelisation, then assumes that the the total number of processors divides the total number of files perfectly:

nprocs = nxpe * nype
n_runs = int(len(filepaths) / nprocs)
if len(filepaths) % nprocs != 0:
    raise ValueError("Each run directory... " 

We can either:

a) Remove this check, which makes the behaviour a bit less explicit.

b) Have a separate logical path for opening only one file. This then excludes the possibly of opening some other subset of all the files (like a single x row), but I don't know why you would want that. Also if you want to only open one file then you can always just use xarray.open_dataset() instead.

It may be possible to apply the BoutCore compiled C++ operators which David Schworer wrote to xBOUT datasets using xr.apply_ufunc().

This could be very powerful as then we could have a set of methods (perhaps accessed as da.boutcore.method()?) which apply the BoutCore operations directly.

This should definitely be @dschwoerer's job, because I have no idea how BoutCore works!

I tried making a gif, but it didn't loop.

This problem was with matplotlib calling Pillow Writer (matplotlib/matplotlib#11787), and the fix was merged in matplotlib 3.1.

We should upgrade our requirements for matplotlib to 3.1.1 (the most recent version) to avoid this odd bug.

Is it possible to test that the installed files are working as expected?

I tried to run pytest with the installed files in PYTHONPATH - but that fails:

collected 9 items / 7 errors / 1 skipped / 1 selected

==================================== ERRORS ====================================
_________________ ERROR collecting xbout/tests/test_animate.py _________________
import file mismatch:
imported module 'xbout.tests.test_animate' has this __file__ attribute:
  /builddir/build/BUILDROOT/python-xBOUT-0.0.1-1.fc33.x86_64/usr/lib/python3.8/site-packages/xbout/tests/test_animate.py
which is not the same as the test file we want to collect:
  /builddir/build/BUILD/xBOUT-master/xbout/tests/test_animate.py
HINT: remove __pycache__ / .pyc files and/or use a unique basename for your test file modules

If that is not possible, that is fine as well, but I would prefer testing the installed files, rather the source.

It shouldn't be a huge issue, so 'Won't fix' is fine to.

There are also some warnings:

/usr/lib/python3.8/site-packages/xbout/load.py:171: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/load.py:556: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/load.py:558: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/tests/test_load.py:290: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/tests/test_load.py:294: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/tests/test_load.py:311: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/load.py:171: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/load.py:556: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/load.py:558: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/tests/test_load.py:290: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/tests/test_load.py:294: SyntaxWarning: "is" with a literal. Did you mean "=="?
/usr/lib/python3.8/site-packages/xbout/tests/test_load.py:311: SyntaxWarning: "is" with a literal. Did you mean "=="?

Full output is here

The coordinates argument to apply_geometry is unnecessary, because its functionality can just be implemented as a type of geometry instead.

This would be simpler, as then the default would just be not to change the dimensions of the dump files at all and not add any coordinates. If you want something different then you have to specify a geometry (and we can provide a default set of choices).

@johnomotani we should just make the default toroidal geometry be optionless. What would be a sensible choice of dimension names then (to not conflict with likely variable names)? 'psi', 'theta', 'zeta'?

The routine to plot the limiter positions attempts to show solid surfaces using hatching, but sometimes ends up looking bad, in particular for unusual configurations such as the TORPEX simulations. In particular the hatching ticks aren't scaled to match the overall size of the plot, sometimes aren't all equally-spaced, and it's not clear whether they should behave differently in the case of including/not including x boundaries.

This is purely cosmetic, but it would be nice to improve this, because then these plots would be much closer to publication-ready.

xarray and dask are obviously hard dependencies, but animatplot could be made optional by only importing it if you try to plot a gif, and natsort could possibly be eliminated by writing our own natural sorting function (instead of the call here).

xBOUT also currently depends on a specific fork of xarray, and this should be merged as soon as possible so that we can just have a dependency on the master version of xarray.

Missing installation instructions from README -- add to top

Me and @johnomotani thought it would make sense to refactor so that open_boutdataset can also open a grid file, and discard open_grid. That's because if you have dump files, you should never actually need to open the grid file in addition. So as it's either grid OR dump files, we should just have one argument to open_boutdataset: filepath.

There are some .nc files in the repo that were originally for testing, but are no longer needed as we create temporary files for tests. They make the repo ~70MB to download, so we should remove them. This requires re-writing the repo history, as described here
https://help.github.com/en/github/authenticating-to-github/removing-sensitive-data-from-a-repository
It is recommended to do this with all PRs closed. We're aiming to do this before the BOUT++ workshop next week, before getting new users, as existing users will need to re-clone or reset/rebase their copies of xBOUT.

BOUT writes out the git commit hash of the compiled executable used into the log files. We could read this and store in the in the dataset attrs.

Ideally BOUT would be changed upstream so that it writes the commit hash into the dump files instead, so that xBOUT doesn't have to be given another filepath argument.

Should we delete the old branches brought over in this fork? They're all behind master with no new commits ahead so we wouldn't lose anything.

xBOUT should be installable with pip.

I don't know how that would work with the current dependency on a specific branch of a forked version of xarray (once it's merged that won't be an issue).

@ZedThree you said this was pretty simple?

If BOUT is changed upstream to write out the global grid point indexes to each file, then xBOUT can use them to infer the order in which the datasets should be concatenated.

This would also require implementing an infer_order_from_coords option in xarray.auto_combine(), as discussed here.

This would move even more of the concatenation logic upstream to xarray, as then we wouldn't need to create a nested list-of-lists anymore, we could just pass the unordered glob of all dump files.

On the other hand this might make the open_mfdataset() call considerably slower, as it will have to read the coordinate variables in order to order the datasets (related discussion here).

I had to add nx, ny, nz to the benchmark file, but after that, I get the following:

      Checking variable flag0 ... Traceback (most recent call last):
  File "./runtest", line 48, in <module>
    result = collect(v, path="data", info=False)
  File "/home/peter/.local/lib/python3.7/site-packages/xbout/load.py", line 207, in collect
    keep_yboundaries=yguards, info=info)
  File "/home/peter/.local/lib/python3.7/site-packages/xbout/load.py", line 284, in _auto_open_mfboutdataset
    chunks=chunks)
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/backends/api.py", line 945, in open_mfdataset
    join=join,
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/combine.py", line 307, in _nested_combine
    join=join,
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/combine.py", line 198, in _combine_nd
    join=join,
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/combine.py", line 220, in _combine_all_along_first_dim
    datasets, dim, compat, data_vars, coords, fill_value, join
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/combine.py", line 248, in _combine_1d
    join=join,
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/concat.py", line 133, in concat
    return f(objs, dim, data_vars, coords, compat, positions, fill_value, join)
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/concat.py", line 340, in _dataset_concat
    var, to_merge[var], compat=compat, equals=equals.get(var, None)
  File "/home/peter/.local/lib/python3.7/site-packages/xarray/core/merge.py", line 143, in unique_variable
    f"conflicting values for variable {name!r} on objects to be combined. "
xarray.core.merge.MergeError: conflicting values for variable 'flag3ad' on objects to be combined. You can skip this check by specifying compat='override'.

Would be nice to format the whole repo with black, and add black as a PR-check. Need to do this at the point where it causes least pain with merging...

The animated 2D plots are useful, and it would be nice to also have these for other types of plots (line, histogram).

These could also use animatplot behind the scenes.

The regions introduced in #107 are based on slicing the Dataset with global indices. The indices stored in Region objects will not be correct if the Dataset is sliced for some reason, e.g. ds = ds.isel(x=slice(20, None).

One way to resolve this might be to add global index coordinates, e.g. for a dataset with dimensions {'x', 'theta', 'zeta'} have some coordinates {'x_ind', 'theta_ind', 'zeta_ind'} that give global integer indices. Then we could use this in Region.getSlices() (which might need an extra Dataset or DataArray argument to get the coordinates from) to get offset. E.g. if the passed-in DataArray has x_ind that does not start from 0, can subtract the da[x_ind][0] from Region.xinner_ind, and check that Region.xinner_ind is within the bounds of x_ind.

Once real-space coordinates are included, we can implement tokamak/BOUT-specific plotting methods.

For example we could make a method
da.bout.plot_poloidal_slice().

Not sure whether xarray's built-in handling of plotting with non-uniform coordinates would be sufficient, or whether it would be harder than that.

Someone has almost certainly done this already, it might be in BOUT's pylib (boutdata.pol_slice() and boututils.plotpolslice() would certainly be useful) or maybe even in OMFIT? Could again have crossover with PlasmaPy?

If we end up reimplementing this then might want to look at the internals of cartopy.

It would be useful to have a geometry='slab' option, as well as the current 'toroidal' and 's-alpha'

Prompted by the discussion here #108 (comment).

'toroidal' should actually work for a slab, but it adds stuff that doesn't necessarily make sense. A simple slab probably just wants coordinates x = dx.isel(y=0).cumsum() - dx.isel(y=0)/2., y = dy.isel(x=0).cumsum() - dy.isel(x=0)/2., z = np.arange(0.5, nz)*dz. It can call _create_toroidal_regions() to create regions in case those are useful (although they should be fairly trivial).

To be implemented after #107 and #108 are merged.

Once BOUT writes out the real space coordinates to the dump files, then xBOUT can read them.

Then xBOUT can return a dataset with coords representing real space coordinates, and dims representing the logical indexes (as it has now).

Help wanted because I don't know much about how the real-space coordinates, metric, Jacobian, etc in BOUT are related to each other.

Currently the unit tests create fake BOUT data and check that xBOUT reads them as expected.

There should also be integrated tests which read some examples of real BOUT data. These should be from a variety of BOUT modules and be of varying dimensions. (SD1D, Storm2D and Hermes would be good candidates.)

The data should be as small as possible.

Help wanted because it would be better if other people gave me the data to perform the tests on.

open_boutdataset() will currently assume you want to use the NETCDF4 engine to read the dump files if they end with .nc, and will throw an error if you try to open a file which has any other extension.

We should identify which types of files this will and won't work with (e. g. hdf5?), and support any others which BOUT can produce and xarray supports.

Should start by looking through boututils.datafile.py.

The bout.save() accessor method will save the whole dataset to either a single file, or one file for each time-evolving variable.

open_boutdataset() needs to have a second path to detect this and load them correctly.

So I'm making some plots where I have different-sized parallel domains all of which I want to centre on the same point in real space (the mid-plane). Opening the data with xstorm returns the parallel coordinate in real space as starting at 0 on the lower target every time. It would be nice if xBOUT/xstorm automatically centred the results on the mid-plane before returning them to me, but I'm not actually sure whether there is enough information in the dump files and input file to know this? BOUT only returns dy right, so we don't know the integration constant? Also how would xBOUT know if I was doing a target-to-target or a mid-plane-to-target simulation?

It would be sweet to be able to make 3D plots of BOUT data. Unfortunately matplotlib3D is not great, and mayavi is a bit outdated (not sure if it's still actively developed?).

But I recently found this library, which runs in a jupyter notebook and uses WebGL as the plotting backend so it's really fast. The examples are insanely cool and the one of a chest cavity in particular is worth a look.

Here's me using K3D-jupyter to plot some filaments from a 3D slab turbulence simulation in like 3 lines of code (after preloading the values DataArray using xbout):

@johnomotani

create_test_data appears to be an undocumented xarray internal function and importing it sets the matplotlib backend. This gives me about 60 lines of warnings when running the tests as I have a specific backend set in my matplotlibrc file.

We should probably reimplement this function in xBOUT to avoid relying on internals of other libraries.

I've come across a problem where Python would crash when trying to interpolate. I think the problem is when the input contains NaN (for me these were coming from corner boundary cells). The help (scipy-1.4.1) for scipy.interpolate.interp1d says "Note that calling interp1d with NaNs present in input values results in undefined behaviour", which sounds like this is a bug that scipy does not intend to fix. That probably makes it an xarray issue really, but I don't have a minimal working example to report it at the moment.

A workaround for my issue is to set keep_xboundaries=False when calling open_boutdataset, so that the corners are not included but we can still have y-boundary cells that are useful for parallel interpolation.

error message says: 'QuadMesh' object has no property 'titles'

I think we should look to rewrite how the subplots are are made in animate_list. #80 is a decent temporary fix for when blank subplots get made. But won't work in all cases.

This fix will only work properly if ncells is a square number and only if ncells is up to 3 greater than nvars. I think this could be generalised better.

The BOUT grid file is just another netCDF file, but we should read it and merge it into the dataset.

Ideally BOUT would be changed upstream so that it writes the grid data into the dump files, so that open_boutdataset() doesn't need another filepath argument.