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[bot]]

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) ────────────────────────────────╮ │ :28 in │ │ │ │ /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.