find_cell fails for a lot of simple input including SPE1 grid generated from flow

[eivindjahren]

The following fails:

from hypothesis import given
import hypothesis.strategies as st

from resdata.grid import Grid
from resdata.grid import GridGenerator

coordinates = st.integers(min_value=1, max_value=30)
increments = st.integers().filter(lambda x: x != 0)

grids = st.builds(
    GridGen.create_rectangular,
    st.tuples(*([coordinates] * 3)),
    st.tuples(*([increments] * 3)),
)

@given(grids)
def test_that_find_cell_finds_the_center_point_of_each_cell(grid):
    for cell in grid.cells():
        assert cell.valid # succeeds
        assert cell.ijk == grid.find_cell(*cell.coordinate) # fails

failing example is GridGenerator.create_rectangular((1, 1, 1), (1, 1, -1)) and test succeeds for

coordinates = st.integers(min_value=1, max_value=30)
increments = st.integers(min_value=1, max_value=10000)

so this definitely has to do with negative increments.

More worryingly, both the check for valid geometry and find_cell does not work with the SPE1 grid generated from flow:

def test_that_flow_grid_works():
    g = Grid(
        ".../ert/test-data/ert/flow_example/spe1_out/realization-0/iter-0/SPE1.EGRID"
    )
    for cell in g.cells():
        # assert cell.valid # fails
        center_point = cell.coordinate
        assert g.find_cell(*center_point) == cell.ijk  # fails

[eivindjahren]

Not sure what the difference is. The following only fails for the grid generated by flow:

cell = list(GridGen.create_rectangular((10, 20, 30), (1, 1, 1)).cells())[0]
print(cell.corners)
# > [(0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (0.0, 0.0, 1.0), (1.0, 0.0, 1.0), (0.0, 1.0, 1.0), (1.0, 1.0, 1.0)]
assert cell.valid # success

cell = list(Grid(
    ".../ert/test-data/ert/flow_example/spe1_out/realization-0/iter-0/SPE1.EGRID"
).cells())[0]
print(cell.corners)
# > [(0.0, 0.0, 8325.0), (1000.0, 0.0, 8325.0), (0.0, 1000.0, 8325.0), (1000.0, 1000.0, 8325.0), (0.0, 0.0, 8345.0), (1000.0, 0.0, 8345.0), (0.0, 1000.0, 8345.0), (1000.0, 1000.0, 8345.0)]
assert cell.valid # fail

[eivindjahren]

With hypothesis it seems that rectangular grids generated by resdata generally works:

The following fails:

from hypothesis import given
import hypothesis.strategies as st

from resdata.grid import Grid
from resdata.grid import GridGenerator

coordinates = st.integers(min_value=1, max_value=30)
increments = st.integers(min_value=1, max_value=10000)

grids = st.builds(
    GridGen.create_rectangular,
    st.tuples(*([coordinates] * 3)),
    st.tuples(*([increments] * 3)),
)

@given(grids)
def test_that_find_cell_finds_the_center_point_of_each_cell(grid):
    for cell in grid.cells():
        assert cell.valid # succeeds
        assert cell.ijk == grid.find_cell(*cell.coordinate) # fails

failing example is GridGenerator.create_rectangular((1, 1, 1), (1, 1, -1)) and test succeeds for

coordinates = st.integers(min_value=1, max_value=30)
increments = st.integers(min_value=1, max_value=10000)

so this definitely has to do with negative increments. However, it does not work again if floating point increments are used:

coordinates = st.integers(min_value=1, max_value=30)
increments = st.floats(min_value=1, max_value=10000)

which fails with the counterexample:

>           assert cell.ijk == grid.find_cell(*cell.coordinate)
E           assert (0, 1, 0) == None
E            +  where (0, 1, 0) = Cell(0, 1, 0, active, (0.500, 1228.892, 3474.750), grid=).ijk

where grid=GridGenerator.create_rectangular((1, 2, 1), (1.0, 819.2613253644572, 6949.5)).

[eivindjahren]

This is due to the heuristic to set cells as invalid:

https://github.com/equinor/resdata/blob/55d970da2aa1e2724f330a109941a1bfafa9a643/lib/resdata/rd_grid.cpp#L917-L984

which then is used to make find_cell return False:

https://github.com/equinor/resdata/blob/55d970da2aa1e2724f330a109941a1bfafa9a643/lib/resdata/rd_grid.cpp#L3576-L3577

[oyvindeide]

Not great, but will not be prioritized now. If there are user reports of this, we can revisit.