Skip to content
Snippets Groups Projects
Commit b0b8e73c authored by Maja Warlich's avatar Maja Warlich
Browse files

Seperate functions, one code for inner and periodic.

parent 4a64a1a0
Branches
No related tags found
No related merge requests found
Pipeline #34268 failed
Stage: test
Expand all Hide whitespace changes
Inline Side-by-side
lbmpy/sparse/mapping.py
+ 85
− 213
View file @ b0b8e73c
......@@ -22,7 +22,7 @@ def pdf_index(cell_index, direction_index, domain_size):
def inverse_idx(stencil, idx):
return stencil.index(tuple(-d_i for d_i in stencil[idx]))
class SparseLbMapper:
"""Manages the mapping of cell coordinates to indices and back.
......@@ -129,8 +129,8 @@ class SparseLbMapper:
no_slip_flag = self.no_slip_flag
fluid_boundary_mask = self.ubb_flag | self.fluid_flag | self.density_flag
result = []
print("flag:", self.flag_array)
fucking_counter = 0
#print("flag:", self.flag_array)
num_border_cell = 0
for direction_idx, direction in enumerate(stencil):
print("dir:", direction_idx)
if all(d_i == 0 for d_i in direction):
......@@ -138,12 +138,13 @@ class SparseLbMapper:
continue
for own_cell_idx, cell in enumerate(self.fluid_coordinates):
if (cell[0] == 0 or cell[0] == self.domain_size[0]-1 or cell[1] == 0 or cell[1] == self.domain_size[1]-1):
fucking_counter += 1 #count skipped border cells for reshape function later
num_border_cell += 1 #count skipped border cells for reshape function later
continue #ignore fluid cells at the border ...
domain_size = (len(self.flag_array), len(self.flag_array[0]))
test = []
inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i in zip(cell, direction)])
inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)])
test.append(("Zelle", own_cell_idx))
print("write:", cell, "read:", inv_neighbor_cell)
if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask:
neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell))
result.append(pdf_index(neighbor_cell_idx, direction_idx, len(self)))
......@@ -172,16 +173,40 @@ class SparseLbMapper:
#print(test)
index_array = np.array(result, dtype=np.uint32)
print("number of fluid:", self.num_fluid_cells, "counter:", fucking_counter//8)
index_arr = index_array.reshape([len(stencil) - 1, self.num_fluid_cells-fucking_counter//8])
#print("number of fluid:", self.num_fluid_cells, "counter:", num_border_cell//8)
#index_arr = index_array.reshape([len(stencil) - 1, self.num_fluid_cells])
index_arr = index_array.reshape([len(stencil) - 1, self.num_fluid_cells-num_border_cell//8])
index_arr = index_arr.swapaxes(0, 1)
return index_arr
class SparseLbPeriodicityMapper:
DIR_SYMBOL = TypedSymbol("dir", np.int64)
def __init__(self, method, mapping: SparseLbMapper, dh):
index_field = Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=np.int64)
#if isinstance(boundary_eqs, Assignment):
# assignments = [boundary_eqs]
#if not isinstance(boundary_eqs, AssignmentCollection):
# assignments = AssignmentCollection(boundary_eqs)
#accessor_subs = dict()
#for fa in assignments.atoms(Field.Access):
# if fa.field == f_out:
# f_dir = 'out'
# elif fa.field == f_in:
# f_dir = 'in'
# else:
# continue
# i = 1 if f_dir == "in" else 2
# accessor_subs[fa] = pdf_field_sparse.absolute_access((index_field(i),), ())
#self.timestep_indexing = indexing
self.index_field = index_field
self.method = method
self.mapping = mapping
self.domain_size = dh.shape
......@@ -191,223 +216,70 @@ class SparseLbPeriodicityMapper:
self.no_slip_flag = mapping.no_slip_flag
self.density_flag = mapping.density_flag
self.ubb_flag = mapping.ubb_flag
def direction_index(self, direction):
direction = tuple(direction)
for direction_idx, dir_coordinates in enumerate(self.mapping.stencil):
if dir_coordinates == direction:
return direction_idx
raise ValueError("Could not find index for direction {}".format(direction))
def get_read_idx(self, read, write_idx, direction_idx):
if self.flag_arr[tuple(read)] & self.no_slip_flag: # Read cell is no-slip: flip PDF!
#read from write cell, inverse direction
return pdf_index(write_idx, inverse_idx(self.method.stencil, direction_idx), len(self.mapping))
#periodic_idx_array.append([direction_idx, inverse_idx(self.method.stencil, direction_idx), write, write]) #nur zu debug Zwecken
else:
return pdf_index(self.mapping.cell_idx(tuple(read)), direction_idx, len(self.mapping))
#periodic_idx_array.append([direction_idx, write, read]) #nur zu debug Zwecken
def get_assignment(self, direction_idx, direction, own_cell_idx, cell):
def get_assignment(self, direction, cell):
cell_idx = self.mapping.cell_idx(cell)
direction_idx = self.direction_index(direction)
inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
write_idx = pdf_index(own_cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, own_cell_idx, direction_idx)
print("write:", cell, "read:", inv_neighbor_cell)
write_idx = pdf_index(cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, cell_idx, direction_idx)
#print("write:", cell, "read:", inv_neighbor_cell)
return [direction_idx, write_idx, read_idx]
def create_inner_index_arr(self):
stencil = self.method.stencil
def fluid_border(self):
fluid_border_coord = []
for cell_idx, cell in enumerate(self.mapping.fluid_coordinates):
if (cell[0] == 0 or cell[0] == self.domain_size[0]-1 or cell[1] == 0 or cell[1] == self.domain_size[1]-1):
fluid_border_coord.append((cell_idx, cell))
result = []
fluid_boundary_mask = self.fluid_flag | self.ubb_flag | self.density_flag
for direction_idx, direction in enumerate(self.method.stencil):
if all(d_i == 0 for d_i in direction):# direction (0,0) irrelevant
continue
print("\n New direction:", direction, ", ", direction_idx)
naughty = [[int((ds_i-1)*(1-dir_i)/2)] if dir_i != 0 else [] for i, (dir_i, ds_i) in enumerate(zip(direction, self.domain_size))]
print(naughty)
for cell_description in fluid_border_coord:
cell = cell_description[1]
own_cell_idx = cell_description[0]
if cell[0] not in naughty[0] and cell[1] not in naughty[1]:
result.append(self.get_assignment(direction_idx, direction, own_cell_idx, cell))
return result
def create_index_arr(self): # erstellt index arrays für ALLE fluid Zellen, die sich am Rand der domain befinden.
# ein index array für alle Werte, die innerhalb des Blocks verschickt werden
# jeweils ein index array für Werte, die zu jeweils verschiedenen benachbarten Blocks geschickt werden (wenn verschiedene Kerne verschiedene Blöcke innerhalb einer Domain bearbeiten)
for cell in self.mapping.fluid_coordinates:
bool_border = [cell_i == 0 or cell_i == ds_i-1 for cell_i, ds_i in zip(cell, self.domain_size)]
if True in bool_border:
fluid_border_coord.append(cell)
return fluid_border_coord
def create_periodic_index_array(self):
stencil = self.method.stencil
print("domain_size:", self.domain_size)
fluid_border_coord = []
for cell_idx, cell in enumerate(self.mapping.fluid_coordinates):
if (cell[0] == 0 or cell[0] == self.domain_size[0]-1 or cell[1] == 0 or cell[1] == self.domain_size[1]-1):
fluid_border_coord.append((cell_idx, cell))
result = []
print(fluid_border_coord)
inner_idx_array = []
fluid_boundary_mask = self.fluid_flag | self.ubb_flag | self.density_flag
for direction_idx, direction in enumerate(self.method.stencil):
if all(d_i == 0 for d_i in direction):# direction (0,0) irrelevant
continue
print("\n New direction:", direction, ", ", direction_idx)
for pos in range(0,2):
print("new periodic_index_array")
periodic_idx_array = []
sop = (pos+1)%2
print("first iteration over cells")
for cell_description in fluid_border_coord:
cell = cell_description[1]
own_cell_idx = cell_description[0]
#print(cell)
if direction[pos] != 0:
slice_coord = int((self.domain_size[pos]-1)*(1-direction[pos])/2) #0 oder d_s-1
start = 1 if direction[sop] == 1 else 0
stop = self.domain_size[sop]-1 if direction[sop] == -1 else self.domain_size[sop]
if cell[pos] == slice_coord and cell[sop] >= start and cell[sop] < stop:
inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
write_idx = pdf_index(own_cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, own_cell_idx, direction_idx)
print("(p ) write:", cell)
periodic_idx_array.append([direction_idx, write_idx, read_idx])
slice_coord = [int((self.domain_size[sop]-1)*(1+direction[sop])/2)] if direction[sop] != 0 else [0, self.domain_size[sop]-1]
if cell[sop] in slice_coord and cell[pos] >= 1 and cell[pos] < self.domain_size[pos]-1:
inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
write_idx = pdf_index(own_cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, own_cell_idx, direction_idx)
print("(i_1) write:", cell)
inner_idx_array.append([direction_idx, write_idx, read_idx])
else: #if direction[pos] == 0
slice_coord = int((self.domain_size[sop]-1)*(1+direction[sop])/2)
if cell[sop] == slice_coord:
inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
write_idx = pdf_index(own_cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, own_cell_idx, direction_idx)
print("(i_2) write:", cell)
inner_idx_array.append([direction_idx, write_idx, read_idx])
print("feed result")
result.append(periodic_idx_array)
#Ecken
if (direction[0] == 0 or direction[1] == 0):
continue
print("second iteration over cells")
for cell_description in fluid_border_coord:
cell = cell_description[1]
own_cell_idx = cell_description[0]
corner = [int((ds_i-1)*(1-dir_i)/2) for dir_i, ds_i in zip(direction, self.domain_size)]
if cell[0] == corner[0] and cell[1] == corner[1]:
inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
write_idx = pdf_index(own_cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, own_cell_idx, direction_idx)
print("(c_p) write:", cell)
result.append([[direction_idx, write_idx, read_idx]])
corner = [int((ds_i-1)*(1+dir_i)/2) for dir_i, ds_i in zip(direction, self.domain_size)]
if cell[0] == corner[0] and cell[1] == corner[1]:
inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
write_idx = pdf_index(own_cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, own_cell_idx, direction_idx)
print("(c_i) write:", cell)
inner_idx_array.append([direction_idx, write_idx, read_idx])
print("End of Code")
for direction_idx, direction in enumerate(self.method.stencil):
if all(d_i == 0 for d_i in direction):# direction (0,0) irrelevant
fluid_border_coord = self.fluid_border()
#print(fluid_border_coord)
result = [[[] for j in range(0, len(stencil)-1)] for i in range(0, len(stencil)-1)]
inner_index_array = []
for direction_idx, direction in enumerate(stencil):
if all(d_i == 0 for d_i in direction): # Direction (0,0) irrelevant
continue
print("\n New direction:", direction, ", ", direction_idx)
for pos in range(0,2): # einmal für x, einmal für y Richtung ...
print("pos is ", pos)
sop = (pos+1)%2
if direction[pos] != 0:
# periodic/parallel: wird an anderen Block geschickt/periodisch gewrappt
print("(periodic:)")
periodic_idx_array = []
coord = int((self.domain_size[pos]-1)*(1-direction[pos])/2)
start = 1 if direction[sop] == 1 else 0
end = self.domain_size[sop]-1 if direction[sop] == -1 else self.domain_size[sop]
for i in range(start, end):
write = [0,0]
write[pos] = coord
write[sop] = i
if not (self.flag_arr[tuple(write)] & self.fluid_flag):
continue
read = [(write_i - dir_i)%ds_i for write_i, dir_i, ds_i in zip(write, direction, self.domain_size)]
print("write:", write)
write_idx = pdf_index(self.mapping.cell_idx(tuple(write)), direction_idx, len(self.mapping))
read_idx = self.get_read_idx(read, self.mapping.cell_idx(tuple(write)), direction_idx)
periodic_idx_array.append([direction_idx, write_idx, read_idx])
# "Die Zelle "write" bekommt ihren neuen Wert der jeweiligen direction von der Zelle "read"
result.append(tuple(periodic_idx_array))
# inner: wird zwischen benachbarten Zellen *im gleichen Block* geschickt
print("(inner1:)")
pos_bound = int((self.domain_size[pos]-1)*(1+direction[pos])/2)
pos_mid = pos_bound+direction[pos]*(-self.domain_size[pos]+1)
sop_position = [int((self.domain_size[sop]-1)*(1+direction[sop])/2)] if direction[sop] != 0 else [0, self.domain_size[sop]-1]
for b in sop_position:
for i in range(pos_bound, pos_mid, -direction[pos]):
write = [0,0]
write[pos] = i
write[sop] = b
print("write:", write)
if not (self.flag_arr[tuple(write)] & self.fluid_flag):
continue
read = [write_i - dir_i for write_i, dir_i in zip(write, direction)]
write_idx = pdf_index(self.mapping.cell_idx(tuple(write)), direction_idx, len(self.mapping))
read_idx = self.get_read_idx(read, write_idx, direction_idx)
inner_idx_array.append([direction_idx, write_idx, read_idx])
if direction[pos] == 0: #spricht directions 1, 2, 3 und 4 an
# inner: wird zwischen benachbarte Zellen *im gleichen Block* geschickt
print("(inner2:)")
pos_low = 1
pos_high = self.domain_size[pos]-1
sop_position = int((self.domain_size[sop]-1)*(1+direction[sop])/2)
for i in range(pos_low, pos_high):
write = [0,0]
write[pos] = i
write[sop] = sop_position
print("write:", write)
if not (self.flag_arr[tuple(write)] & self.fluid_flag):
continue
read = [write_i - dir_i for write_i, dir_i in zip(write, direction)]
write_idx = pdf_index(self.mapping.cell_idx(tuple(write)), direction_idx, len(self.mapping))
read_idx = self.get_read_idx(read, write_idx, direction_idx)
inner_idx_array.append([direction_idx, write, read]) #for debug
#Four corners: extra periodic_idx_array for each direction 5, 6, 7, 8
if (direction[0]*direction[1] != 0):
write = [int((self.domain_size[0]-1)*(1-direction[0])/2),int((self.domain_size[1]-1)*(1-direction[1])/2)]
print("write:", write)
if not (self.flag_arr[tuple(write)] & self.fluid_flag):
continue
read = [(write_i - dir_i)%ds_i for write_i, dir_i, ds_i in zip(write, direction, self.domain_size)]
write_idx = pdf_index(self.mapping.cell_idx(tuple(write)), direction_idx, len(self.mapping))
read_idx = self.get_read_idx(read, write_idx, direction_idx)
periodic_idx_array.append([direction_idx, write_idx, read_idx])
result.append(tuple(periodic_idx_array))
# result enthält *mehrere* index arrays
#result = list(dict.fromkeys(result)) # entferne doppelte index_arrays: speziell Ecken der Domain
# result ist eine liste von tuples von tuples --> [((...), (...)), ((...), (...), (...))]
#print(result)
# wandel result in list_result (liste von liste von listen) um: -->[[[...], [...]], [[...], [...], [...]]]
list_result = []
for index_array in result:
list_index_array = []
for write_read_pair in index_array:
list_index_array.append(list(write_read_pair))
list_result.append(list_index_array)
# zu den periodischen/parralel-orientierten index_arrays kommt noch der index array für die Werte, die nur innerhalb des Blocks verschickt werden:
list_result.append(inner_idx_array)
for index_array in list_result:
print(index_array)
return list_result
#print("\n New direction:", direction, ", ", direction_idx)
periodic_slice_coord = [[int((ds_i-1)*(1-dir_i)/2)] if dir_i != 0 else [] for i, (dir_i, ds_i) in enumerate(zip(direction, self.domain_size))]
for cell in fluid_border_coord:
bool_pos = [cell_i in slice_i for i, (cell_i, slice_i) in enumerate(zip(cell, periodic_slice_coord))]
#print(cell, ", ", periodic_slice_coord, ", ", bool_pos)
if True in bool_pos: # This cell & direction is periodical
block_direction = self.direction_index([int(bp_i)*dir_i for dir_i, bp_i in zip(direction, bool_pos)])
result[direction_idx-1][block_direction-1].append(self.get_assignment(direction, cell))
#print("Goes into", direction_idx-1, block_direction-1)
else: # This cell & direction is not periodical
inner_index_array.append(self.get_assignment(direction, cell))
#print("Inner")
result.append([inner_index_array])
#for a in result:
# print(a)
return result
def assignments(self):
return [Assignment(self.index_field(1), self.index_field(2))]
def get_kernel(self):
kernel = create_kernel(self.assignments(), ghost_layers=0)
return kernel
class SparseLbBoundaryMapper:
NEIGHBOR_IDX_NAME = 'nidx{}'
......@@ -467,7 +339,7 @@ class SparseLbBoundaryMapper:
self.index_field_dtype = index_field_dtype
self.neighbor_offsets = neighbor_offsets
self.index_field = index_field
self.timesptep_indexing = indexing
self.timestep_indexing = indexing
self.boundary_functor = boundary
def _build_substitutions(self):
......@@ -482,7 +354,7 @@ class SparseLbBoundaryMapper:
for dir_idx, subs_dict in enumerate(result):
inv_idx = stencil.index(tuple(-e for e in stencil[dir_idx]))
subs_dict[self.timesptep_indexing._index_array_symbol("in", True)] = inv_idx
subs_dict[self.timestep_indexing._index_array_symbol("in", True)] = inv_idx
return result
def create_index_arr(self, mapping: SparseLbMapper, boundary_mask, nr_of_ghost_layers=1):
......@@ -521,7 +393,7 @@ class SparseLbBoundaryMapper:
def get_kernel(self):
kernel = create_kernel(self.assignments(), ghost_layers=0)
index_arrs_node = self.timesptep_indexing.create_code_node()
index_arrs_node = self.timestep_indexing.create_code_node()
for node in self.boundary_functor.get_additional_code_nodes(self.method)[::-1]:
kernel.body.insert_front(node)
kernel.body.insert_front(index_arrs_node)
......
lbmpy_tests/simple_shear_flow.ipynb 0 → 100644
+ 303
− 0
View file @ b0b8e73c
This diff is collapsed.
lbmpy_tests/simple_sparse_shear_flow.ipynb 0 → 100644
+ 364
− 0
View file @ b0b8e73c
This diff is collapsed.
lbmpy_tests/test_sparse_lbm.ipynb
+ 312
− 180
View file @ b0b8e73c
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment