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:

resdata/lib/resdata/rd_grid.cpp

Lines 917 to 984 in 55d970d

	/**
	The problem is that some grids have invalid cells. Typically
	the cells accomodating numerical aquifers are located at an utm
	position (0,0).
	Cells which have some pillars located in (0,0) and some cells
	located among the rest of the grid become completely warped - with
	insane volumes, parts of the reservoir volume doubly covered, and
	so on.
	To keep these cells out of the real-world (i.e. involving utm
	coordinates) computations they are marked as 'tainted' in this
	function. The tainting procedure is completely heuristic, and
	probably wrong.
	--------------
	There is second heuristic which marks cells as invalid. In some
	cases (which ??) cells outside the area of interest are just set to
	have all four corner at the same arbitrary depth; these cells are
	inactive and do not affect flow simulations - however the arbitrary
	location of the cells warps visualisation of normal inactive cells
	completely. We therefore try to invalidate such cells here. The
	algorithm used is the same as used for RMS; however RMS will mark
	the cells as inactive - whereas we mark already inactive cells as
	invalid.
	It is very important that this mechanism is NOT used to set cells
	to inactive; as that will break completely when it comes to
	e.g. data like PRESSURE and SWAT which are defined as a number of
	active elements.
	*/
	static void rd_cell_taint_cell(rd_cell_type *cell) {
	int c;
	for (c = 0; c < 8; c++) {
	const point_type p = cell->corner_list[c];
	if ((p.x == 0) && (p.y == 0)) {
	SET_CELL_FLAG(cell, CELL_FLAG_TAINTED);
	break;
	}
	}
	/*
	Second heuristic to invalidate cells.
	*/
	if (cell->active == CELL_NOT_ACTIVE) {
	if (!GET_CELL_FLAG(cell, CELL_FLAG_TAINTED)) {
	const point_type p0 = cell->corner_list[0];
	int cell_index = 1;
	while (true) {
	const point_type pi = cell->corner_list[cell_index];
	if (pi.z != p0.z)
	// There is a difference - the cell is certainly valid.
	break;
	else {
	cell_index++;
	if (cell_index == 8) {
	// They have all been at the same height up until now;
	// the cell is marked as invalid.
	SET_CELL_FLAG(cell, CELL_FLAG_TAINTED);
	break;
	}
	}
	}
	}
	}
	}

which then is used to make find_cell return False:

resdata/lib/resdata/rd_grid.cpp

Lines 3576 to 3577 in 55d970d

	if (GET_CELL_FLAG(cell, CELL_FLAG_TAINTED))
	return false;

[oyvindeide]

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