Use closest normal for boundary index list with single_link

For creating the index list just the first stencil entry was taken which is a neighbour of the investigated cell if single_link=True. With this MR the discrete normal is calculated and the neighbouring cell in the normal direction is taken to build up the index array.

Furthermore, the computational cost of the python versions for create_boundary_index_list is reduced drastically because the iteration space is now restricted to the boundary cells and not the entire domain anymore.

pystencils/boundaries/createindexlist.py

@@ -47,21 +47,27 @@ def _create_index_list_python(flag_field_arr, boundary_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)])
 
+    # boundary cells are extracted via np.where. To ensure continous memory access in the compute kernel these cells
+    # have to be sorted.
     boundary_cells = np.transpose(np.where(flag_field_arr == boundary_mask))
     for i in range(len(flag_field_arr.shape)):
         boundary_cells = boundary_cells[boundary_cells[:, i].argsort(kind='mergesort')]
 
-    fluid_cells = []
+    # 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):
             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 + nr_of_ghost_layers <= e < upper - nr_of_ghost_layers
                    for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
                 continue
-            if flag_field_arr[neighbor_cell] & fluid_mask and fluid_cells.count(neighbor_cell) == 0:
-                fluid_cells.append(neighbor_cell)
+            if flag_field_arr[neighbor_cell] & fluid_mask:
+                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:
@@ -76,6 +82,7 @@ def _create_index_list_python(flag_field_arr, boundary_mask,
         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)])
+            # 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)):
                 continue
             if flag_field_arr[neighbor_cell] & checkmask:
@@ -84,6 +91,7 @@ def _create_index_list_python(flag_field_arr, boundary_mask,
                 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,))
@@ -104,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
         inner_or_boundary: if true, the result contains the cell coordinates of the domain cells -
                     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)