Skip to content

Commit

Permalink
Merge pull request #807 from bp/combine_tmult_optimise
Browse files Browse the repository at this point in the history 
optimising the combining of tr mults from multiple grid connection sets
  • Loading branch information
@andy-beer
andy-beer committed Jul 22, 2024
2 parents 4394e14 + 7232a99 commit 8394bb0
Show file tree
Hide file tree
Showing 4 changed files with 106 additions and 87 deletions.
12 changes: 6 additions & 6 deletions resqpy/fault/_gcs_functions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -439,29 +439,29 @@ def combined_tr_mult_from_gcs_mults(model,

# merge in each of the three directional face arrays for this gcs with combined arrays
for (combo_trm, gcs_trm) in [(combo_trm_k, gcs_trm_k), (combo_trm_j, gcs_trm_j), (combo_trm_i, gcs_trm_i)]:
mask = np.logical_not(np.isnan(gcs_trm)) # true where this tr mult is present
mask = np.logical_not(np.isnan(gcs_trm)).astype(bool) # true where this tr mult is present
clash_mask = np.logical_and(mask, np.logical_not(np.isnan(combo_trm))) # true where combined value clashes
if np.any(clash_mask):
if merge_mode == 'exception':
raise ValueError('gcs transmissibility multiplier conflict when merging')
if merge_mode == 'minimum':
combo_trm[:] = np.where(clash_mask, np.minimum(combo_trm, gcs_trm), combo_trm)
combo_trm[clash_mask] = np.minimum(combo_trm, gcs_trm)[clash_mask]
elif merge_mode == 'maximum':
combo_trm[:] = np.where(clash_mask, np.maximum(combo_trm, gcs_trm), combo_trm)
combo_trm[clash_mask] = np.maximum(combo_trm, gcs_trm)[clash_mask]
elif merge_mode == 'multiply':
combo_trm[:] = np.where(clash_mask, combo_trm * gcs_trm, combo_trm)
combo_trm[clash_mask] = (combo_trm * gcs_trm)[clash_mask]
else:
raise Exception(f'code failure with unrecognised merge mode {merge_mode}')
mask = np.logical_and(mask,
np.logical_not(clash_mask)) # remove clash faces from mask (already handled)
if np.any(mask):
combo_trm[:] = np.where(mask, gcs_trm, combo_trm) # update combined array from individual array
combo_trm[mask] = gcs_trm[mask] # update combined array from individual array

# for each of the 3 combined tr mult arrays, replace unused values with the default fill value
# also check that any set values are non-negative
for combo_trm in (combo_trm_k, combo_trm_j, combo_trm_i):
if fill_value is not None and not np.isnan(fill_value):
combo_trm[:] = np.where(np.isnan(combo_trm), fill_value, combo_trm)
combo_trm[:][np.isnan(combo_trm).astype(bool)] = fill_value
assert np.all(combo_trm >= 0.0)
else:
assert np.all(np.logical_or(np.isnan(combo_trm), combo_trm >= 0.0))
Expand Down
67 changes: 29 additions & 38 deletions resqpy/fault/_grid_connection_set.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2100,9 +2100,6 @@ def grid_face_arrays(self,
ak = np.full((nk + 1, nj, ni), default_value, dtype = dtype)
aj = np.full((nk, nj + 1, ni), default_value, dtype = dtype)
ai = np.full((nk, nj, ni + 1), default_value, dtype = dtype)
# mk = np.zeros((nk + 1, nj, ni), dtype = bool)
# mj = np.zeros((nk, nj + 1, ni), dtype = bool)
# mi = np.zeros((nk, nj, ni + 1), dtype = bool)

# populate arrays from faces of gcs, optionally filtered by feature index
cip, fip = self.list_of_cell_face_pairs_for_feature_index(None)
Expand All @@ -2113,42 +2110,36 @@ def grid_face_arrays(self,
else:
indices = self.indices_for_feature_index(feature_index)

# opposing_count = 0
side_list = ([0] if lazy else [0, 1])
for fi in indices:
# fi = int(i)
if active_mask is not None and not active_mask[fi]:
continue
value = gcs_prop_array[fi]
if baffle_mask is not None and baffle_mask[fi]:
value = 0 # will be cast to float (or bool) if needed when assigned below
for side in side_list:
cell_kji0 = cip[fi, side].copy()
# opposing = cell_kji0.copy()
axis, polarity = fip[fi, side]
assert 0 <= axis <= 2 and 0 <= polarity <= 1
cell_kji0[axis] += polarity
# opposing[axis] += (1 - polarity)
if axis == 0:
ak[tuple(cell_kji0)] = value
# mk[tuple(cell_kji0)] = True
# if mk[tuple(opposing)]:
# opposing_count += 1
elif axis == 1:
aj[tuple(cell_kji0)] = value
# mj[tuple(cell_kji0)] = True
# if mj[tuple(opposing)]:
# opposing_count += 1
else:
ai[tuple(cell_kji0)] = value
# mi[tuple(cell_kji0)] = True
# if mi[tuple(opposing)]:
# opposing_count += 1

# if opposing_count:
# log.warning(f'{opposing_count} suspicious opposing faces of {len(indices)} detected in gcs: {self.title}')
# else:
# log.debug(f'no suspicious opposing faces detected in gcs: {self.title}')

value_array = gcs_prop_array.copy()
if baffle_mask is not None:
value_array[baffle_mask] = 0 # will be cast to float (or bool) if needed
if active_mask is not None:
cip = cip[active_mask, :, :]
value_array = value_array[active_mask]

for side in side_list:
cell_kji0 = cip[:, side].copy() # shape (N, 3)
axis = fip[:, side, 0] # shape (N,)
polarity = fip[:, side, 1] # shape (N,)
assert 0 <= np.min(axis) and np.max(axis) <= 2
assert 0 <= np.min(polarity) and np.max(polarity) <= 1

axis_mask = (axis == 0).astype(bool)
ak_kji0 = cell_kji0[axis_mask, :]
ak_kji0[:, 0] += polarity[axis_mask]
ak[ak_kji0[:, 0], ak_kji0[:, 1], ak_kji0[:, 2]] = value_array[axis_mask]

axis_mask = (axis == 1).astype(bool)
aj_kji0 = cell_kji0[axis_mask, :]
aj_kji0[:, 1] += polarity[axis_mask]
aj[aj_kji0[:, 0], aj_kji0[:, 1], aj_kji0[:, 2]] = value_array[axis_mask]

axis_mask = (axis == 2).astype(bool)
ai_kji0 = cell_kji0[axis_mask, :]
ai_kji0[:, 2] += polarity[axis_mask]
ai[ai_kji0[:, 0], ai_kji0[:, 1], ai_kji0[:, 2]] = value_array[axis_mask]

return (ak, aj, ai)

Expand Down
6 changes: 3 additions & 3 deletions resqpy/grid/_grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2255,9 +2255,9 @@ def combined_tr_mult_properties_from_gcs_mults(self,
pc = self.extract_property_collection()

if baffle_triplet is not None:
trm_k[1:-1] = np.where(baffle_triplet[0], 0.0, trm_k[1:-1])
trm_j[:, 1:-1] = np.where(baffle_triplet[1], 0.0, trm_j[:, 1:-1])
trm_i[:, :, 1:-1] = np.where(baffle_triplet[2], 0.0, trm_i[:, :, 1:-1])
trm_k[1:-1][baffle_triplet[0]] = 0.0
trm_j[:, 1:-1][baffle_triplet[1]] = 0.0
trm_i[:, :, 1:-1][baffle_triplet[2]] = 0.0

if composite_property:
tr_composite = np.concatenate((trm_k.flat, trm_j.flat, trm_i.flat))
Expand Down
108 changes: 68 additions & 40 deletions tests/unit_tests/fault/test_fault_connection_set.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,9 +120,11 @@ def test_fault_connection_set(tmp_path):
# model.store_epc()

g2_fcs, g2_fa = rqtr.fault_connection_set(g2)
g2_fcs_again, _ = rqtr.fault_connection_set(g2)

assert g2_fcs is not None
assert g2_fa is not None
assert g2_fcs_again is not None

# show_fa(g2, g2_fcs, g2_fa)

Expand All @@ -132,6 +134,8 @@ def test_fault_connection_set(tmp_path):
# check grid faces property array generation
g2_fcs.write_hdf5()
g2_fcs.create_xml()
g2_fcs_again.write_hdf5()
g2_fcs_again.create_xml()
gcs_p_a = np.array((37, 51), dtype = int)
p = rqp.Property.from_array(model,
gcs_p_a,
Expand All @@ -143,6 +147,7 @@ def test_fault_connection_set(tmp_path):
discrete = True,
null_value = -1)
gcs_trm_a = np.array((0.1, 0.5), dtype = float)
gcs_trm_b = np.array((0.2, 0.4), dtype = float)
trm_p = rqp.Property.from_array(model,
gcs_trm_a,
'test',
Expand All @@ -152,6 +157,15 @@ def test_fault_connection_set(tmp_path):
indexable_element = 'faces',
discrete = False,
uom = 'Euc')
trm_again_p = rqp.Property.from_array(model,
gcs_trm_b,
'test',
'trmult',
g2_fcs_again.uuid,
property_kind = 'transmissibility multiplier',
indexable_element = 'faces',
discrete = False,
uom = 'Euc')
for lazy in [True, False]:
pk, pj, pi = g2_fcs.grid_face_arrays(property_uuid = p.uuid,
default_value = p.null_value(),
Expand All @@ -163,46 +177,60 @@ def test_fault_connection_set(tmp_path):
assert np.all(pi == -1)
assert np.count_nonzero(pj > 0) == 2 if lazy else 4
assert tuple(np.unique(pj)) == (-1, 37, 51)
for merge_mode in ['minimum', 'maximum', 'multiply']:
for sided in ([False] if merge_mode == 'multiply' else [False, True]):
c_trm_k, c_trm_j, c_trm_i = \
rqf.combined_tr_mult_from_gcs_mults(model,
[trm_p.uuid],
merge_mode = merge_mode,
sided = sided,
fill_value = 1.0)
assert all([a is not None for a in (c_trm_k, c_trm_j, c_trm_i)])
assert c_trm_k.shape == (3, 2, 2) and c_trm_j.shape == (2, 3, 2) and c_trm_i.shape == (2, 2, 3)
assert_array_almost_equal(c_trm_k, 1.0)
assert_array_almost_equal(c_trm_i, 1.0)
assert np.count_nonzero(c_trm_j < 0.9) == 2 if lazy else 4
assert not np.any(np.isnan(c_trm_j))
g_trm_k_uuid, g_trm_j_uuid, g_trm_i_uuid = \
g2.combined_tr_mult_properties_from_gcs_mults(gcs_uuid_list = [g2_fcs.uuid],
merge_mode = merge_mode,
sided = sided,
fill_value = 1.0,
composite_property = False)
assert g_trm_k_uuid is not None and g_trm_j_uuid is not None and g_trm_i_uuid is not None
g_trm_k = rqp.Property(model, uuid = g_trm_k_uuid).array_ref()
g_trm_j = rqp.Property(model, uuid = g_trm_j_uuid).array_ref()
g_trm_i = rqp.Property(model, uuid = g_trm_i_uuid).array_ref()
assert g_trm_k is not None and g_trm_j is not None and g_trm_i is not None
assert np.all(g_trm_k == c_trm_k)
assert np.all(g_trm_j == c_trm_j)
assert np.all(g_trm_i == c_trm_i)
g_trm_list = g2.combined_tr_mult_properties_from_gcs_mults(gcs_uuid_list = [g2_fcs.uuid],
merge_mode = merge_mode,
sided = sided,
fill_value = 1.0,
composite_property = True)
assert len(g_trm_list) == 1
g_trm = rqp.Property(model, uuid = g_trm_list[0]).array_ref()
assert g_trm.ndim == 1
assert g_trm.size == g_trm_k.size + g_trm_j.size + g_trm_i.size
assert np.all(g_trm[:g_trm_k.size] == g_trm_k.flat)
assert np.all(g_trm[g_trm_k.size:g_trm_k.size + g_trm_j.size] == g_trm_j.flat)
assert np.all(g_trm[g_trm_k.size + g_trm_j.size:] == g_trm_i.flat)
for pair, gcs_uuid_list, tr_uuid_list in [(False, [g2_fcs.uuid], [trm_p.uuid]),
(True, [g2_fcs.uuid, g2_fcs_again.uuid], [trm_p.uuid, trm_again_p.uuid])]:
for merge_mode in ['minimum', 'maximum', 'multiply']:
for sided in ([False] if merge_mode == 'multiply' else [False, True]):
c_trm_k, c_trm_j, c_trm_i = \
rqf.combined_tr_mult_from_gcs_mults(model,
tr_uuid_list,
merge_mode = merge_mode,
sided = sided,
fill_value = 1.0)
assert all([a is not None for a in (c_trm_k, c_trm_j, c_trm_i)])
assert c_trm_k.shape == (3, 2, 2) and c_trm_j.shape == (2, 3, 2) and c_trm_i.shape == (2, 2, 3)
assert_array_almost_equal(c_trm_k, 1.0)
assert_array_almost_equal(c_trm_i, 1.0)
assert np.count_nonzero(c_trm_j < 0.9) == 2 if lazy else 4
assert not np.any(np.isnan(c_trm_j))
unique_c_trm_j = np.unique(c_trm_j)
assert len(unique_c_trm_j) == 3
if merge_mode == 'minimum':
assert np.isclose(np.min(c_trm_j), 0.1)
assert np.isclose(unique_c_trm_j[1], 0.4 if pair else 0.5)
elif merge_mode == 'maximum':
assert np.isclose(np.min(c_trm_j), 0.2 if pair else 0.1)
assert np.isclose(unique_c_trm_j[1], 0.5)
else: # multiply
assert np.isclose(np.min(c_trm_j), 0.02 if pair else 0.1)
assert np.isclose(unique_c_trm_j[1], 0.2 if pair else 0.5)
assert np.isclose(unique_c_trm_j[2], 1.0)
g_trm_k_uuid, g_trm_j_uuid, g_trm_i_uuid = \
g2.combined_tr_mult_properties_from_gcs_mults(gcs_uuid_list = gcs_uuid_list,
merge_mode = merge_mode,
sided = sided,
fill_value = 1.0,
composite_property = False)
assert g_trm_k_uuid is not None and g_trm_j_uuid is not None and g_trm_i_uuid is not None
g_trm_k = rqp.Property(model, uuid = g_trm_k_uuid).array_ref()
g_trm_j = rqp.Property(model, uuid = g_trm_j_uuid).array_ref()
g_trm_i = rqp.Property(model, uuid = g_trm_i_uuid).array_ref()
assert g_trm_k is not None and g_trm_j is not None and g_trm_i is not None
assert np.all(g_trm_k == c_trm_k)
assert np.all(g_trm_j == c_trm_j)
assert np.all(g_trm_i == c_trm_i)
g_trm_list = g2.combined_tr_mult_properties_from_gcs_mults(gcs_uuid_list = gcs_uuid_list,
merge_mode = merge_mode,
sided = sided,
fill_value = 1.0,
composite_property = True)
assert len(g_trm_list) == 1
g_trm = rqp.Property(model, uuid = g_trm_list[0]).array_ref()
assert g_trm.ndim == 1
assert g_trm.size == g_trm_k.size + g_trm_j.size + g_trm_i.size
assert np.all(g_trm[:g_trm_k.size] == g_trm_k.flat)
assert np.all(g_trm[g_trm_k.size:g_trm_k.size + g_trm_j.size] == g_trm_j.flat)
assert np.all(g_trm[g_trm_k.size + g_trm_j.size:] == g_trm_i.flat)

# I face split with full juxtaposition of kji0 (1, *, 0) with (0, *, 1)
# pattern 4, 4 (or 3, 3) diagram 1
Expand Down

0 comments on commit 8394bb0

Please sign in to comment.