Skip to content
Snippets Groups Projects

Refactor packaging, part I

3 files
+ 109
55
Compare changes
  • Side-by-side
  • Inline
Files
3
@@ -2,26 +2,22 @@ import warnings
@@ -2,26 +2,22 @@ import warnings
import numpy as np
import numpy as np
 
try:
try:
# Try to import right away - assume compiled code is available
import pyximport
# compile with: python setup.py build_ext --inplace --use-cython
from pystencils.boundaries.createindexlistcython import create_boundary_neighbor_index_list_2d, \
create_boundary_neighbor_index_list_3d, create_boundary_cell_index_list_2d, create_boundary_cell_index_list_3d
 
pyximport.install(language_level=3)
cython_funcs_available = True
cython_funcs_available = True
except ImportError:
except ImportError:
try:
cython_funcs_available = False
# If not, try development mode and import via pyximport
import pyximport
if cython_funcs_available:
from pystencils.boundaries.createindexlistcython import (
pyximport.install(language_level=3)
create_boundary_neighbor_index_list_2d,
cython_funcs_available = True
create_boundary_neighbor_index_list_3d,
except ImportError:
create_boundary_cell_index_list_2d,
cython_funcs_available = False
create_boundary_cell_index_list_3d,
if cython_funcs_available:
)
from pystencils.boundaries.createindexlistcython import create_boundary_neighbor_index_list_2d, \
create_boundary_neighbor_index_list_3d, create_boundary_cell_index_list_2d, \
create_boundary_cell_index_list_3d
boundary_index_array_coordinate_names = ["x", "y", "z"]
boundary_index_array_coordinate_names = ["x", "y", "z"]
direction_member_name = "dir"
direction_member_name = "dir"
@@ -30,40 +26,59 @@ default_index_array_dtype = np.int32
@@ -30,40 +26,59 @@ default_index_array_dtype = np.int32
def numpy_data_type_for_boundary_object(boundary_object, dim):
def numpy_data_type_for_boundary_object(boundary_object, dim):
coordinate_names = boundary_index_array_coordinate_names[:dim]
coordinate_names = boundary_index_array_coordinate_names[:dim]
return np.dtype([(name, default_index_array_dtype) for name in coordinate_names]
return np.dtype(
+ [(direction_member_name, default_index_array_dtype)]
[(name, default_index_array_dtype) for name in coordinate_names]
+ [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data], align=True)
+ [(direction_member_name, default_index_array_dtype)]
+ [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data],
align=True,
def _create_index_list_python(flag_field_arr, boundary_mask,
)
fluid_mask, stencil, single_link, inner_or_boundary=False, nr_of_ghost_layers=None):
 
def _create_index_list_python(
 
flag_field_arr,
 
boundary_mask,
 
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:
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 "
raise ValueError(
"around the inner domain has to be specified")
"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:
if nr_of_ghost_layers is None:
nr_of_ghost_layers = 0
nr_of_ghost_layers = 0
coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)]
coordinate_names = boundary_index_array_coordinate_names[
index_arr_dtype = np.dtype([(name, default_index_array_dtype) for name in coordinate_names]
: len(flag_field_arr.shape)
+ [(direction_member_name, default_index_array_dtype)])
]
 
index_arr_dtype = np.dtype(
 
[(name, default_index_array_dtype) for name in coordinate_names]
 
+ [(direction_member_name, default_index_array_dtype)]
 
)
# boundary cells are extracted via np.where. To ensure continous memory access in the compute kernel these cells
# boundary cells are extracted via np.where. To ensure continous memory access in the compute kernel these cells
# have to be sorted.
# have to be sorted.
boundary_cells = np.transpose(np.nonzero(flag_field_arr == boundary_mask))
boundary_cells = np.transpose(np.nonzero(flag_field_arr == boundary_mask))
for i in range(len(flag_field_arr.shape)):
for i in range(len(flag_field_arr.shape)):
boundary_cells = boundary_cells[boundary_cells[:, i].argsort(kind='mergesort')]
boundary_cells = boundary_cells[boundary_cells[:, i].argsort(kind="mergesort")]
# First a set is created to save all fluid cells which are near boundary
# First a set is created to save all fluid cells which are near boundary
fluid_cells = set()
fluid_cells = set()
for cell in boundary_cells:
for cell in boundary_cells:
cell = tuple(cell)
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)]
 
)
# prevent out ouf bounds access. If boundary cell is at the border, some stencil directions would be out.
# 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
if any(
for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
not 0 + nr_of_ghost_layers <= e < upper - nr_of_ghost_layers
 
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] & fluid_mask:
fluid_cells.add(neighbor_cell)
fluid_cells.add(neighbor_cell)
@@ -83,9 +98,14 @@ def _create_index_list_python(flag_field_arr, boundary_mask,
@@ -83,9 +98,14 @@ def _create_index_list_python(flag_field_arr, boundary_mask,
cell = tuple(cell)
cell = tuple(cell)
sum_cells = np.zeros(len(cell))
sum_cells = np.zeros(len(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)]
 
)
# prevent out ouf bounds access. If boundary cell is at the border, some stencil directions would be out.
# 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] & checkmask:
if flag_field_arr[neighbor_cell] & checkmask:
if single_link:
if single_link:
@@ -101,8 +121,15 @@ def _create_index_list_python(flag_field_arr, boundary_mask,
@@ -101,8 +121,15 @@ def _create_index_list_python(flag_field_arr, boundary_mask,
return np.array(result, dtype=index_arr_dtype)
return np.array(result, dtype=index_arr_dtype)
def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
def create_boundary_index_list(
nr_of_ghost_layers=1, inner_or_boundary=True, single_link=False):
flag_field,
 
stencil,
 
boundary_mask,
 
fluid_mask,
 
nr_of_ghost_layers=1,
 
inner_or_boundary=True,
 
single_link=False,
 
):
"""Creates a numpy array storing links (connections) between domain cells and boundary cells.
"""Creates a numpy array storing links (connections) between domain cells and boundary cells.
Args:
Args:
@@ -119,11 +146,20 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
@@ -119,11 +146,20 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
"""
"""
dim = len(flag_field.shape)
dim = len(flag_field.shape)
coordinate_names = boundary_index_array_coordinate_names[:dim]
coordinate_names = boundary_index_array_coordinate_names[:dim]
index_arr_dtype = np.dtype([(name, default_index_array_dtype) for name in coordinate_names]
index_arr_dtype = np.dtype(
+ [(direction_member_name, default_index_array_dtype)])
[(name, default_index_array_dtype) for name in coordinate_names]
 
+ [(direction_member_name, default_index_array_dtype)]
 
)
stencil = np.array(stencil, dtype=default_index_array_dtype)
stencil = np.array(stencil, dtype=default_index_array_dtype)
args = (flag_field, nr_of_ghost_layers, boundary_mask, fluid_mask, stencil, single_link)
args = (
 
flag_field,
 
nr_of_ghost_layers,
 
boundary_mask,
 
fluid_mask,
 
stencil,
 
single_link,
 
)
args_no_gl = (flag_field, boundary_mask, fluid_mask, stencil, single_link)
args_no_gl = (flag_field, boundary_mask, fluid_mask, stencil, single_link)
if cython_funcs_available:
if cython_funcs_available:
@@ -142,22 +178,42 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
@@ -142,22 +178,42 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
return np.array(idx_list, dtype=index_arr_dtype)
return np.array(idx_list, dtype=index_arr_dtype)
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(
return _create_index_list_python(*args_no_gl, inner_or_boundary=inner_or_boundary,
"Boundary setup may take very long! Consider installing cython to speed it up"
nr_of_ghost_layers=nr_of_ghost_layers)
)
return _create_index_list_python(
*args_no_gl,
def create_boundary_index_array(flag_field, stencil, boundary_mask, fluid_mask, boundary_object,
inner_or_boundary=inner_or_boundary,
nr_of_ghost_layers=1, inner_or_boundary=True, single_link=False):
nr_of_ghost_layers=nr_of_ghost_layers,
idx_array = create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
)
nr_of_ghost_layers, inner_or_boundary, single_link)
 
 
def create_boundary_index_array(
 
flag_field,
 
stencil,
 
boundary_mask,
 
fluid_mask,
 
boundary_object,
 
nr_of_ghost_layers=1,
 
inner_or_boundary=True,
 
single_link=False,
 
):
 
idx_array = create_boundary_index_list(
 
flag_field,
 
stencil,
 
boundary_mask,
 
fluid_mask,
 
nr_of_ghost_layers,
 
inner_or_boundary,
 
single_link,
 
)
dim = len(flag_field.shape)
dim = len(flag_field.shape)
if boundary_object.additional_data:
if boundary_object.additional_data:
coordinate_names = boundary_index_array_coordinate_names[:dim]
coordinate_names = boundary_index_array_coordinate_names[:dim]
index_arr_dtype = numpy_data_type_for_boundary_object(boundary_object, dim)
index_arr_dtype = numpy_data_type_for_boundary_object(boundary_object, dim)
extended_idx_field = np.empty(len(idx_array), dtype=index_arr_dtype)
extended_idx_field = np.empty(len(idx_array), dtype=index_arr_dtype)
for prop in coordinate_names + ['dir']:
for prop in coordinate_names + ["dir"]:
extended_idx_field[prop] = idx_array[prop]
extended_idx_field[prop] = idx_array[prop]
idx_array = extended_idx_field
idx_array = extended_idx_field
Loading