Skip to content
Snippets Groups Projects
Commit cb9cccfb authored by Michael Kuron's avatar Michael Kuron :mortar_board:
Browse files

Merge branch 'CreateIndexList' into 'master'

Use closest normal for boundary index list with single_link

See merge request pycodegen/pystencils!253
parents f06a6c30 2b88b63a
No related branches found
No related tags found
No related merge requests found
import itertools
import warnings import warnings
import numpy as np import numpy as np
...@@ -35,45 +34,67 @@ def numpy_data_type_for_boundary_object(boundary_object, dim): ...@@ -35,45 +34,67 @@ def numpy_data_type_for_boundary_object(boundary_object, dim):
+ [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data], align=True) + [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data], align=True)
def _create_boundary_neighbor_index_list_python(flag_field_arr, nr_of_ghost_layers, boundary_mask, def _create_index_list_python(flag_field_arr, boundary_mask,
fluid_mask, stencil, single_link): fluid_mask, stencil, single_link, inner_or_boundary=False, nr_of_ghost_layers=None):
if inner_or_boundary and nr_of_ghost_layers is None:
raise ValueError("If inner_or_boundary is set True the number of ghost layers "
"around the inner domain has to be specified")
if nr_of_ghost_layers is None:
nr_of_ghost_layers = 0
coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)] coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)]
index_arr_dtype = np.dtype([(name, np.int32) for name in coordinate_names] + [(direction_member_name, np.int32)]) index_arr_dtype = np.dtype([(name, np.int32) for name in coordinate_names] + [(direction_member_name, np.int32)])
result = [] # boundary cells are extracted via np.where. To ensure continous memory access in the compute kernel these cells
gl = nr_of_ghost_layers # have to be sorted.
for cell in itertools.product(*reversed([range(gl, i - gl) for i in flag_field_arr.shape])): boundary_cells = np.transpose(np.nonzero(flag_field_arr == boundary_mask))
cell = cell[::-1] for i in range(len(flag_field_arr.shape)):
if not flag_field_arr[cell] & fluid_mask: boundary_cells = boundary_cells[boundary_cells[:, i].argsort(kind='mergesort')]
continue
# First a set is created to save all fluid cells which are near boundary
fluid_cells = set()
for cell in boundary_cells:
cell = tuple(cell)
for dir_idx, direction in enumerate(stencil): for dir_idx, direction in enumerate(stencil):
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)]) neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
if flag_field_arr[neighbor_cell] & boundary_mask: # prevent out ouf bounds access. If boundary cell is at the border, some stencil directions would be out.
result.append(cell + (dir_idx,)) if any(not 0 + nr_of_ghost_layers <= e < upper - nr_of_ghost_layers
if single_link: for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
break continue
if flag_field_arr[neighbor_cell] & fluid_mask:
return np.array(result, dtype=index_arr_dtype) fluid_cells.add(neighbor_cell)
# then this is set is transformed to a list to make it sortable. This ensures continoous memory access later.
fluid_cells = list(fluid_cells)
if len(flag_field_arr.shape) == 3:
fluid_cells.sort(key=lambda tup: (tup[-1], tup[-2], tup[0]))
else:
fluid_cells.sort(key=lambda tup: (tup[-1], tup[0]))
def _create_boundary_cell_index_list_python(flag_field_arr, boundary_mask, cells_to_iterate = fluid_cells if inner_or_boundary else boundary_cells
fluid_mask, stencil, single_link): checkmask = boundary_mask if inner_or_boundary else fluid_mask
coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)]
index_arr_dtype = np.dtype([(name, np.int32) for name in coordinate_names] + [(direction_member_name, np.int32)])
result = [] result = []
for cell in itertools.product(*reversed([range(0, i) for i in flag_field_arr.shape])): for cell in cells_to_iterate:
cell = cell[::-1] cell = tuple(cell)
if not flag_field_arr[cell] & boundary_mask: sum_cells = np.zeros(len(cell))
continue
for dir_idx, direction in enumerate(stencil): for dir_idx, direction in enumerate(stencil):
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)]) neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
# prevent out ouf bounds access. If boundary cell is at the border, some stencil directions would be out.
if any(not 0 <= e < upper for e, upper in zip(neighbor_cell, flag_field_arr.shape)): if any(not 0 <= e < upper for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
continue continue
if flag_field_arr[neighbor_cell] & fluid_mask: if flag_field_arr[neighbor_cell] & checkmask:
result.append(cell + (dir_idx,))
if single_link: if single_link:
break sum_cells += np.array(direction)
else:
result.append(tuple(cell) + (dir_idx,))
# the discrete normal direction is the one which gives the maximum inner product to the stencil direction
if single_link and any(sum_cells != 0):
idx = np.argmax(np.inner(sum_cells, stencil))
result.append(tuple(cell) + (idx,))
return np.array(result, dtype=index_arr_dtype) return np.array(result, dtype=index_arr_dtype)
...@@ -91,7 +112,7 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask, ...@@ -91,7 +112,7 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
nr_of_ghost_layers: only relevant if neighbors is True nr_of_ghost_layers: only relevant if neighbors is True
inner_or_boundary: if true, the result contains the cell coordinates of the domain cells - inner_or_boundary: if true, the result contains the cell coordinates of the domain cells -
if false the boundary cells are listed if false the boundary cells are listed
single_link: if true only the first link is reported from this cell single_link: if true only the link in normal direction to this cell is reported
""" """
dim = len(flag_field.shape) dim = len(flag_field.shape)
...@@ -119,10 +140,8 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask, ...@@ -119,10 +140,8 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
else: else:
if flag_field.size > 1e6: if flag_field.size > 1e6:
warnings.warn("Boundary setup may take very long! Consider installing cython to speed it up") warnings.warn("Boundary setup may take very long! Consider installing cython to speed it up")
if inner_or_boundary: return _create_index_list_python(*args_no_gl, inner_or_boundary=inner_or_boundary,
return _create_boundary_neighbor_index_list_python(*args) nr_of_ghost_layers=nr_of_ghost_layers)
else:
return _create_boundary_cell_index_list_python(*args_no_gl)
def create_boundary_index_array(flag_field, stencil, boundary_mask, fluid_mask, boundary_object, def create_boundary_index_array(flag_field, stencil, boundary_mask, fluid_mask, boundary_object,
......
...@@ -22,20 +22,37 @@ def create_boundary_neighbor_index_list_2d(object[IntegerType, ndim=2] flag_fiel ...@@ -22,20 +22,37 @@ def create_boundary_neighbor_index_list_2d(object[IntegerType, ndim=2] flag_fiel
cdef int xs, ys, x, y cdef int xs, ys, x, y
cdef int dirIdx, num_directions, dx, dy cdef int dirIdx, num_directions, dx, dy
cdef int sum_x, sum_y
cdef float dot, maxn
cdef int calculated_idx
xs, ys = flag_field.shape xs, ys = flag_field.shape
boundary_index_list = [] boundary_index_list = []
num_directions = stencil.shape[0] num_directions = stencil.shape[0]
for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers): for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers):
for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers): for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers):
sum_x = 0; sum_y = 0;
if flag_field[x, y] & fluid_mask: if flag_field[x, y] & fluid_mask:
for dirIdx in range(num_directions): for dirIdx in range(num_directions):
dx = stencil[dirIdx,0] dx = stencil[dirIdx,0]; dy = stencil[dirIdx,1]
dy = stencil[dirIdx,1]
if flag_field[x + dx, y + dy] & boundary_mask: if flag_field[x + dx, y + dy] & boundary_mask:
boundary_index_list.append((x,y, dirIdx))
if single_link: if single_link:
break sum_x += dx; sum_y += dy;
else:
boundary_index_list.append((x, y, dirIdx))
dot = 0; maxn = 0; calculated_idx = 0
if single_link and (sum_x != 0 or sum_y != 0):
for dirIdx in range(num_directions):
dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1];
dot = dx * sum_x + dy * sum_y
if dot > maxn:
maxn = dot
calculated_idx = dirIdx
boundary_index_list.append((x, y, calculated_idx))
return boundary_index_list return boundary_index_list
...@@ -47,6 +64,10 @@ def create_boundary_neighbor_index_list_3d(object[IntegerType, ndim=3] flag_fiel ...@@ -47,6 +64,10 @@ def create_boundary_neighbor_index_list_3d(object[IntegerType, ndim=3] flag_fiel
cdef int xs, ys, zs, x, y, z cdef int xs, ys, zs, x, y, z
cdef int dirIdx, num_directions, dx, dy, dz cdef int dirIdx, num_directions, dx, dy, dz
cdef int sum_x, sum_y, sum_z
cdef float dot, maxn
cdef int calculated_idx
xs, ys, zs = flag_field.shape xs, ys, zs = flag_field.shape
boundary_index_list = [] boundary_index_list = []
num_directions = stencil.shape[0] num_directions = stencil.shape[0]
...@@ -54,15 +75,27 @@ def create_boundary_neighbor_index_list_3d(object[IntegerType, ndim=3] flag_fiel ...@@ -54,15 +75,27 @@ def create_boundary_neighbor_index_list_3d(object[IntegerType, ndim=3] flag_fiel
for z in range(nr_of_ghost_layers, zs - nr_of_ghost_layers): for z in range(nr_of_ghost_layers, zs - nr_of_ghost_layers):
for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers): for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers):
for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers): for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers):
sum_x = 0; sum_y = 0; sum_z = 0
if flag_field[x, y, z] & fluid_mask: if flag_field[x, y, z] & fluid_mask:
for dirIdx in range(num_directions): for dirIdx in range(num_directions):
dx = stencil[dirIdx,0] dx = stencil[dirIdx,0]; dy = stencil[dirIdx,1]; dz = stencil[dirIdx,2]
dy = stencil[dirIdx,1]
dz = stencil[dirIdx,2]
if flag_field[x + dx, y + dy, z + dz] & boundary_mask: if flag_field[x + dx, y + dy, z + dz] & boundary_mask:
boundary_index_list.append((x,y,z, dirIdx))
if single_link: if single_link:
break sum_x += dx; sum_y += dy; sum_z += dz
else:
boundary_index_list.append((x, y, z, dirIdx))
dot = 0; maxn = 0; calculated_idx = 0
if single_link and (sum_x != 0 or sum_y != 0 or sum_z != 0):
for dirIdx in range(num_directions):
dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]; dz = stencil[dirIdx, 2]
dot = dx * sum_x + dy * sum_y + dz * sum_z
if dot > maxn:
maxn = dot
calculated_idx = dirIdx
boundary_index_list.append((x, y, z, calculated_idx))
return boundary_index_list return boundary_index_list
...@@ -75,21 +108,39 @@ def create_boundary_cell_index_list_2d(object[IntegerType, ndim=2] flag_field, ...@@ -75,21 +108,39 @@ def create_boundary_cell_index_list_2d(object[IntegerType, ndim=2] flag_field,
cdef int xs, ys, x, y cdef int xs, ys, x, y
cdef int dirIdx, num_directions, dx, dy cdef int dirIdx, num_directions, dx, dy
cdef int sum_x, sum_y
cdef float dot, maxn
cdef int calculated_idx
xs, ys = flag_field.shape xs, ys = flag_field.shape
boundary_index_list = [] boundary_index_list = []
num_directions = stencil.shape[0] num_directions = stencil.shape[0]
for y in range(0, ys): for y in range(0, ys):
for x in range(0, xs): for x in range(0, xs):
sum_x = 0; sum_y = 0;
if flag_field[x, y] & boundary_mask: if flag_field[x, y] & boundary_mask:
for dirIdx in range(num_directions): for dirIdx in range(num_directions):
dx = stencil[dirIdx,0] dx = stencil[dirIdx,0]; dy = stencil[dirIdx,1]
dy = stencil[dirIdx,1]
if 0 <= x + dx < xs and 0 <= y + dy < ys: if 0 <= x + dx < xs and 0 <= y + dy < ys:
if flag_field[x + dx, y + dy] & fluid_mask: if flag_field[x + dx, y + dy] & fluid_mask:
boundary_index_list.append((x,y, dirIdx))
if single_link: if single_link:
break sum_x += dx; sum_y += dy
else:
boundary_index_list.append((x, y, dirIdx))
dot = 0; maxn = 0; calculated_idx = 0
if single_link and (sum_x != 0 or sum_y != 0):
for dirIdx in range(num_directions):
dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]
dot = dx * sum_x + dy * sum_y
if dot > maxn:
maxn = dot
calculated_idx = dirIdx
boundary_index_list.append((x, y, calculated_idx))
return boundary_index_list return boundary_index_list
...@@ -101,6 +152,10 @@ def create_boundary_cell_index_list_3d(object[IntegerType, ndim=3] flag_field, ...@@ -101,6 +152,10 @@ def create_boundary_cell_index_list_3d(object[IntegerType, ndim=3] flag_field,
cdef int xs, ys, zs, x, y, z cdef int xs, ys, zs, x, y, z
cdef int dirIdx, num_directions, dx, dy, dz cdef int dirIdx, num_directions, dx, dy, dz
cdef int sum_x, sum_y, sum_z
cdef float dot, maxn
cdef int calculated_idx
xs, ys, zs = flag_field.shape xs, ys, zs = flag_field.shape
boundary_index_list = [] boundary_index_list = []
num_directions = stencil.shape[0] num_directions = stencil.shape[0]
...@@ -108,14 +163,27 @@ def create_boundary_cell_index_list_3d(object[IntegerType, ndim=3] flag_field, ...@@ -108,14 +163,27 @@ def create_boundary_cell_index_list_3d(object[IntegerType, ndim=3] flag_field,
for z in range(0, zs): for z in range(0, zs):
for y in range(0, ys): for y in range(0, ys):
for x in range(0, xs): for x in range(0, xs):
sum_x = 0; sum_y = 0; sum_z = 0
if flag_field[x, y, z] & boundary_mask: if flag_field[x, y, z] & boundary_mask:
for dirIdx in range(num_directions): for dirIdx in range(num_directions):
dx = stencil[dirIdx,0] dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]; dz = stencil[dirIdx, 2]
dy = stencil[dirIdx,1]
dz = stencil[dirIdx,2]
if 0 <= x + dx < xs and 0 <= y + dy < ys and 0 <= z + dz < zs: if 0 <= x + dx < xs and 0 <= y + dy < ys and 0 <= z + dz < zs:
if flag_field[x + dx, y + dy, z + dz] & fluid_mask: if flag_field[x + dx, y + dy, z + dz] & fluid_mask:
boundary_index_list.append((x,y,z, dirIdx))
if single_link: if single_link:
break sum_x += dx; sum_y += dy; sum_z += dz
else:
boundary_index_list.append((x, y, z, dirIdx))
dot = 0; maxn = 0; calculated_idx=0
if single_link and (sum_x != 0 or sum_y !=0 or sum_z !=0):
for dirIdx in range(num_directions):
dx = stencil[dirIdx, 0]; dy = stencil[dirIdx, 1]; dz = stencil[dirIdx, 2]
dot = dx*sum_x + dy*sum_y + dz*sum_z
if dot > maxn:
maxn = dot
calculated_idx = dirIdx
boundary_index_list.append((x, y, z, calculated_idx))
return boundary_index_list return boundary_index_list
\ No newline at end of file
...@@ -26,11 +26,11 @@ def test_equivalence_cython_python_version(single_link): ...@@ -26,11 +26,11 @@ def test_equivalence_cython_python_version(single_link):
flag_field_3d[-1, :, :] = mask flag_field_3d[-1, :, :] = mask
flag_field_3d[7, 7, 7] = mask flag_field_3d[7, 7, 7] = mask
result_python_2d = cil._create_boundary_neighbor_index_list_python(flag_field_2d, 1, mask, fluid_mask, result_python_2d = cil._create_index_list_python(flag_field_2d, mask, fluid_mask,
stencil_2d, single_link) stencil_2d, single_link, True, 1)
result_python_3d = cil._create_boundary_neighbor_index_list_python(flag_field_3d, 1, mask, fluid_mask, result_python_3d = cil._create_index_list_python(flag_field_3d, mask, fluid_mask,
stencil_3d, single_link) stencil_3d, single_link, True, 1)
result_cython_2d = cil.create_boundary_index_list(flag_field_2d, stencil_2d, mask, result_cython_2d = cil.create_boundary_index_list(flag_field_2d, stencil_2d, mask,
fluid_mask, 1, True, single_link) fluid_mask, 1, True, single_link)
...@@ -62,11 +62,11 @@ def test_equivalence_cell_idx_list_cython_python_version(single_link): ...@@ -62,11 +62,11 @@ def test_equivalence_cell_idx_list_cython_python_version(single_link):
flag_field_3d[-1, :, :] = mask flag_field_3d[-1, :, :] = mask
flag_field_3d[7, 7, 7] = mask flag_field_3d[7, 7, 7] = mask
result_python_2d = cil._create_boundary_cell_index_list_python(flag_field_2d, mask, fluid_mask, result_python_2d = cil._create_index_list_python(flag_field_2d, mask, fluid_mask,
stencil_2d, single_link) stencil_2d, single_link, False)
result_python_3d = cil._create_boundary_cell_index_list_python(flag_field_3d, mask, fluid_mask, result_python_3d = cil._create_index_list_python(flag_field_3d, mask, fluid_mask,
stencil_3d, single_link) stencil_3d, single_link, False)
result_cython_2d = cil.create_boundary_index_list(flag_field_2d, stencil_2d, mask, fluid_mask, None, result_cython_2d = cil.create_boundary_index_list(flag_field_2d, stencil_2d, mask, fluid_mask, None,
False, single_link) False, single_link)
...@@ -75,3 +75,42 @@ def test_equivalence_cell_idx_list_cython_python_version(single_link): ...@@ -75,3 +75,42 @@ def test_equivalence_cell_idx_list_cython_python_version(single_link):
np.testing.assert_equal(result_python_2d, result_cython_2d) np.testing.assert_equal(result_python_2d, result_cython_2d)
np.testing.assert_equal(result_python_3d, result_cython_3d) np.testing.assert_equal(result_python_3d, result_cython_3d)
@pytest.mark.parametrize('inner_or_boundary', [False, True])
def test_normal_calculation(inner_or_boundary):
stencil = tuple((x,y) for x,y in product([-1, 0, 1], [-1, 0, 1]))
domain_size = (32, 32)
dtype = np.uint32
fluid_mask = dtype(1)
mask = dtype(2)
flag_field = np.ones([domain_size[0], domain_size[1]], dtype=dtype) * fluid_mask
radius_inner = domain_size[0] // 4
radius_outer = domain_size[0] // 2
y_mid = domain_size[1] / 2
x_mid = domain_size[0] / 2
for x in range(0, domain_size[0]):
for y in range(0, domain_size[1]):
if (y - y_mid) ** 2 + (x - x_mid) ** 2 < radius_inner ** 2:
flag_field[x, y] = mask
if (x - x_mid) ** 2 + (y - y_mid) ** 2 > radius_outer ** 2:
flag_field[x, y] = mask
args_no_gl = (flag_field, mask, fluid_mask, np.array(stencil, dtype=np.int32), True)
index_list = cil._create_index_list_python(*args_no_gl, inner_or_boundary=inner_or_boundary, nr_of_ghost_layers=1)
checkmask = mask if inner_or_boundary else fluid_mask
for cell in index_list:
idx = cell[2]
cell = tuple((cell[0], cell[1]))
sum_cells = np.zeros(len(cell))
for dir_idx, direction in enumerate(stencil):
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
if any(not 0 <= e < upper for e, upper in zip(neighbor_cell, flag_field.shape)):
continue
if flag_field[neighbor_cell] & checkmask:
sum_cells += np.array(direction)
assert np.argmax(np.inner(sum_cells, stencil)) == idx
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment