Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Plotting grids using open_boutdataset seems broken #306

Open
mrhardman opened this issue Nov 12, 2024 · 1 comment · May be fixed by #307
Open

Plotting grids using open_boutdataset seems broken #306

mrhardman opened this issue Nov 12, 2024 · 1 comment · May be fixed by #307
Labels
bug Something isn't working

Comments

@mrhardman
Copy link
Collaborator

There is a suggested example for plotting the fields of a grid.nc file https://github.com/boutproject/xBOUT/blob/master/examples/plot_grid.py. I have tried to use this code to plot the output from a simple Hyponotoad-generated grid. The result is that I received error messages which suggest that the advertised functionality is broken.

To generate the grid I use (thanks to J. Omotani)

from hypnotoad.cases.circular import CircularEquilibrium
from hypnotoad.core.mesh import BoutMesh
import numpy as np

R0 = 2.3
B0 = 3.2
r_inner = 0.1
r_outer = 1.4
q = 4.1
settings = {
    "R0": R0,
    "B0": B0,
    "poloidal_spacing_method": "linear",
    "q_coefficients": [q],
    "r_inner": r_inner,
    "r_outer": r_outer,
    "refine_methods": "line",
    "refine_width": 1.0e-3,
}

equilib = CircularEquilibrium(settings, nonorthogonal_settings=settings)

mesh = BoutMesh(equilib, settings)
mesh.geometry()
mesh.writeGridfile("grid.nc")

I then run the following script to try to plot

from matplotlib import pyplot as plt
from xbout import open_boutdataset
#from xbout.load import _open_grid

gridfilepath = "grid.nc"
grid = open_boutdataset(gridfilepath, geometry="toroidal")
#grid = _open_grid(gridfilepath, chunks={},keep_xboundaries=True,keep_yboundaries=True)

grid["psi_poloidal"].bout.contourf()
grid["psi_poloidal"].bout.contour()
grid["psi_poloidal"].bout.pcolormesh()
grid["psi_poloidal"].bout.pcolormesh(shading="gouraud")
grid["psi_poloidal"].bout.regions()
plt.show()

With this script, I get the following error

$ python3 plot_test.py
NXPE not found, setting to 1
NYPE not found, setting to 1
MXG not found, setting to 2
MYG not found, setting to 0
MXSUB not found, setting to 2
MYSUB not found, setting to 8
Traceback (most recent call last):
  File "*/lib/python3.13/site-packages/xarray/core/dataset.py", line 1317, in _construct_dataarray
    variable = self._variables[name]
               ~~~~~~~~~~~~~~~^^^^^^
KeyError: 't_array'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "*/plot_test.py", line 11, in <module>
    grid = open_boutdataset(gridfilepath, geometry="toroidal")
  File "*/lib/python3.13/site-packages/xbout/load.py", line 279, in open_boutdataset
    ds, remove_yboundaries = _auto_open_mfboutdataset(
                             ~~~~~~~~~~~~~~~~~~~~~~~~^
        datapath=datapath,
        ^^^^^^^^^^^^^^^^^^
    ...<4 lines>...
        **kwargs,
        ^^^^^^^^^
    )
    ^
  File "*/lib/python3.13/site-packages/xbout/load.py", line 736, in _auto_open_mfboutdataset
    _, unique_indices = unique(ds["t_array"], return_index=True)
                               ~~^^^^^^^^^^^
  File "*/lib/python3.13/site-packages/xarray/core/dataset.py", line 1410, in __getitem__
    return self._construct_dataarray(key)
           ~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^
  File "*/lib/python3.13/site-packages/xarray/core/dataset.py", line 1319, in _construct_dataarray
    _, name, variable = _get_virtual_variable(self._variables, name, self.dims)
                        ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "*/lib/python3.13/site-packages/xarray/core/dataset.py", line 175, in _get_virtual_variable
    raise KeyError(key)
KeyError: 't_array'

Trying to use _open_grid directly also failed. Fixing this would be very helpful for plotting Hypnotoad output.

@mrhardman mrhardman added the bug Something isn't working label Nov 12, 2024
@johnomotani
Copy link
Collaborator

This problem was mostly caused by a change to hypnotoad which had resulted in a 't' dimension being created in grid files, which confused the function that detects whether the file being opened is 'dump', 'restart' or 'grid'. Also a name-conflict error.

@johnomotani johnomotani linked a pull request Nov 13, 2024 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants