Point function sometimes fails to return values on cell edges

[zacjohnston]

Bug report

Bug summary

When using the point function to sample a FLASH checkpoint file, it will sometimes return an empty array as if the point is outside the domain. It seems to happen when the point lies on a cell edge, even when that edge is well inside the grid domain (though it doesn't happen on every edge).

Code for reproduction

data file: curl -JO http://use.yt/upload/ae71a17a

import yt

var = 'dens'
filename = 'sedov_hdf5_chk_0000'

data = yt.load(filename)
grids = data.index.grids

for b, grid in enumerate(grids):
    print(f'Grid left edge: {grid.LeftEdge}')
    print(f'Grid right edge: {grid.RightEdge}')

    nx = grid.ActiveDimensions[0]
    ny = grid.ActiveDimensions[1]

    grid_width = grid.RightEdge - grid.LeftEdge
    dx = grid_width[0] / nx
    dy = grid_width[1] / ny

    x_low = grid.LeftEdge[0]
    y_low = grid.LeftEdge[1]

    for i in range(nx):
        for j in range(ny):
            x = x_low + i*dx
            y = y_low + j*dy

            point = data.point([x, y, 0])

            if len(point['flash', var]) == 0:
                print(f'point returns null in grid {b} at x={x:.3f}, y={y:.3f}')

Actual outcome

Grid left edge: [0. 0. 0.] code_length
Grid right edge: [1. 1. 1.] code_length
point returns null in grid 0 at x=0.000, y=0.667
point returns null in grid 0 at x=0.167, y=0.667
point returns null in grid 0 at x=0.333, y=0.667
point returns null in grid 0 at x=0.500, y=0.667
point returns null in grid 0 at x=0.667, y=0.000
point returns null in grid 0 at x=0.667, y=0.167
point returns null in grid 0 at x=0.667, y=0.333
point returns null in grid 0 at x=0.667, y=0.500
point returns null in grid 0 at x=0.667, y=0.667
point returns null in grid 0 at x=0.667, y=0.833
point returns null in grid 0 at x=0.833, y=0.667
Grid left edge: [0. 0. 0.] code_length
Grid right edge: [0.5 0.5 1. ] code_length
point returns null in grid 1 at x=0.000, y=0.083
point returns null in grid 1 at x=0.083, y=0.000
point returns null in grid 1 at x=0.083, y=0.083
point returns null in grid 1 at x=0.083, y=0.167
[...]

Expected outcome

Always return a value at all points within the domain.

Version Information

• Operating System: Ubuntu 20.04
• Python Version: 3.7.7
• yt version: 3.4.1

yt was installed with conda install yt

[welcome]

Hi, and welcome to yt! Thanks for opening your first issue. We have an issue template that helps us to gather relevant information to help diagnosing and fixing the issue.

[neutrinoceros]

well it somehow seems to have gotten worse, so there may be more than one bug to solve here now
Here's what I get by running the embedded script with yt 4.0.1 + unyt 2.8.0

IterableUnitCoercionError: Received a list or tuple of quantities with nonuniform units: [unyt_quantity(0., 'code_length'), unyt_quantity(0., 'code_length'), 0]

Detailed traceback

╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ <ipython-input-1-17587cb3d73b>:28 in <module>                                                    │
│                                                                                                  │
│ /Users/robcleme/dev/yt-project/yt/yt/data_objects/selection_objects/point.py:52 in __init__      │
│                                                                                                  │
│   49 │   │   │   # we pass p through ds.arr to ensure code units are attached                    │
│   50 │   │   │   self.p = self.ds.arr(p)                                                         │
│   51 │   │   else:                                                                               │
│ ❱ 52 │   │   │   self.p = self.ds.arr(p, "code_length")                                          │
│   53                                                                                             │
│                                                                                                  │
│ /Users/robcleme/.pyenv/versions/yt-dev/lib/python3.9/site-packages/unyt/array.py:553 in __new__  │
│                                                                                                  │
│    550 │   │   │   pass                                                                          │
│    551 │   │   elif _iterable(input_array) and input_array:                                      │
│    552 │   │   │   if isinstance(input_array[0], unyt_array):                                    │
│ ❱  553 │   │   │   │   return _coerce_iterable_units(input_array, registry)                      │
│    554 │   │                                                                                     │
│    555 │   │   # Input array is an already formed ndarray instance                               │
│    556 │   │   # We first cast to be our class type                                              │
│                                                                                                  │
│ /Users/robcleme/.pyenv/versions/yt-dev/lib/python3.9/site-packages/unyt/array.py:248 in          │
│ _coerce_iterable_units                                                                           │
│                                                                                                  │
│    245 │   │   if any([isinstance(o, unyt_array) for o in input_object]):                        │
│    246 │   │   │   ff = getattr(input_object[0], "units", NULL_UNIT)                             │
│    247 │   │   │   if any([ff != getattr(_, "units", NULL_UNIT) for _ in input_object]):         │
│ ❱  248 │   │   │   │   raise IterableUnitCoercionError(input_object)                             │
│    249 │   │   │   # This will create a copy of the data in the iterable.                        │
│    250 │   │   │   return unyt_array(np.array(input_object), ff, registry=registry)              │
│    251 │   return np.asarray(input_object)                                                       │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯

edit: this is due to unyt becoming more strict about inputs, it is fixed with the following change in OP's script

-             point = data.point([x, y, 0])
+             point = data.point([x, y, data.quan(0, "code_length")])

with this patch, I report that I can reproduce the original problem.

[neutrinoceros]

I went down the rabbit hole and I am still not sure wether this is a frontend-specific issue.
What I know so far is that in yt.geometry.geometry_handler.Index._read_fluid_fields, we get some chunks (yt.geometry.geometry_handler.YTDataChunk) with a data_size of 0.
I think these chunks are retrieved via yt.geometry.grid_geometry_handler.GridIndex._identify_base_chunk, which we could override in FlashHierarchy if need be, but that's all I have so far.