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

Modify sparseLbMapper to ghost_layers.

parent f372cad7
No related branches found
No related tags found
No related merge requests found
...@@ -40,8 +40,10 @@ class SparseLbMapper: ...@@ -40,8 +40,10 @@ class SparseLbMapper:
self.no_slip_flag = no_slip_flag self.no_slip_flag = no_slip_flag
self.other_boundary_mask = other_boundary_mask self.other_boundary_mask = other_boundary_mask
self._num_fluid_cells = None self._num_fluid_cells = None
self._num_ghost_cells = None
self.stencil = stencil self.stencil = stencil
self.domain_size = domain_size self.domain_size = domain_size
self.ghost_layers = None
@property @property
...@@ -55,7 +57,10 @@ class SparseLbMapper: ...@@ -55,7 +57,10 @@ class SparseLbMapper:
def fluid_coordinates(self): def fluid_coordinates(self):
if self._dirty: if self._dirty:
self._assemble() self._assemble()
return self._coordinate_arr[:self.num_fluid_cells] all_fluid_coordinates = self._coordinate_arr[:self.num_fluid_cells]
non_ghost_fluid_coord = [f for f in self._coordinate_arr[:self.num_fluid_cells] if True not in [(x_i == 0) or (x_i == self._flag_arr.shape[i]-1) for i, x_i in enumerate(f)]]
#return self._coordinate_arr[:self.num_fluid_cells]
return non_ghost_fluid_coord
@property @property
def num_fluid_cells(self): def num_fluid_cells(self):
...@@ -109,7 +114,8 @@ class SparseLbMapper: ...@@ -109,7 +114,8 @@ class SparseLbMapper:
coordinates_fluid = np.argwhere(np.bitwise_and(self._flag_arr, self.fluid_flag)).astype(np.uint32) coordinates_fluid = np.argwhere(np.bitwise_and(self._flag_arr, self.fluid_flag)).astype(np.uint32)
coordinates_boundary = np.argwhere(np.bitwise_and(self._flag_arr, self.other_boundary_mask)).astype(np.uint32) coordinates_boundary = np.argwhere(np.bitwise_and(self._flag_arr, self.other_boundary_mask)).astype(np.uint32)
self._num_fluid_cells = coordinates_fluid.shape[0] self._num_fluid_cells = coordinates_fluid.shape[0]
self._num_ghost_cells = np.prod(np.array(tuple(ds + 2 * self.ghost_layers for ds in self.domain_size)))\
- np.prod(np.array(self.domain_size))
total_cells = len(coordinates_fluid) + len(coordinates_boundary) total_cells = len(coordinates_fluid) + len(coordinates_boundary)
self._coordinate_arr = np.empty((total_cells,), dtype=struct_type) self._coordinate_arr = np.empty((total_cells,), dtype=struct_type)
...@@ -122,6 +128,7 @@ class SparseLbMapper: ...@@ -122,6 +128,7 @@ class SparseLbMapper:
def create_index_array(self, ghost_layers=1): def create_index_array(self, ghost_layers=1):
# TODO support different layouts here # TODO support different layouts here
self.ghost_layers = ghost_layers
stencil = self.stencil stencil = self.stencil
flag_arr = self.flag_array flag_arr = self.flag_array
no_slip_flag = self.no_slip_flag no_slip_flag = self.no_slip_flag
...@@ -134,10 +141,10 @@ class SparseLbMapper: ...@@ -134,10 +141,10 @@ class SparseLbMapper:
continue continue
for own_cell_idx, cell in enumerate(self.fluid_coordinates): for own_cell_idx, cell in enumerate(self.fluid_coordinates):
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 in zip(cell, direction)])
at_border = [n_cell_i == -1 or n_cell_i == ds_i for n_cell_i, ds_i in zip(inv_neighbor_cell, self.domain_size)] #at_border = [n_cell_i == 0 or n_cell_i == ds_i+1 for n_cell_i, ds_i in zip(cell, self.domain_size)]
if (True in at_border):
inv_neighbor_cell = list(cell)
#print("write:", cell, "read:", inv_neighbor_cell) #print("write:", cell, "read:", inv_neighbor_cell)
#if (True in at_border):
# continue
if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask: if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask:
neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell)) neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell))
result.append(pdf_index(neighbor_cell_idx, direction_idx, len(self))) result.append(pdf_index(neighbor_cell_idx, direction_idx, len(self)))
...@@ -147,31 +154,33 @@ class SparseLbMapper: ...@@ -147,31 +154,33 @@ class SparseLbMapper:
result.append(pdf_index(own_cell_idx, inverse_idx(stencil, direction_idx), len(self))) result.append(pdf_index(own_cell_idx, inverse_idx(stencil, direction_idx), len(self)))
#print("write:", cell, "read:", cell, "direction:", inverse_idx(stencil, direction_idx)) #print("write:", cell, "read:", cell, "direction:", inverse_idx(stencil, direction_idx))
else: else:
#This is where we end up if inv_neighbor_cell is not fluid/no_slip/at the border; #This is where we end up if inv_neighbor_cell is not fluid/no_slip/ghost cell;
at_border = False #legacy periodicity handler
for i, x_i in enumerate(inv_neighbor_cell): #print("N is not normal")
if x_i == (ghost_layers - 1): #at_border = False
inv_neighbor_cell[i] += flag_arr.shape[i] - (2 * ghost_layers) #for i, x_i in enumerate(inv_neighbor_cell):
at_border = True # if x_i == (ghost_layers - 1):
elif x_i == flag_arr.shape[i] - ghost_layers: # inv_neighbor_cell[i] += flag_arr.shape[i] - (2 * ghost_layers)
inv_neighbor_cell[i] -= flag_arr.shape[i] - (2 * ghost_layers) # at_border = True
at_border = True # elif x_i == flag_arr.shape[i] - ghost_layers:
if at_border: # inv_neighbor_cell[i] -= flag_arr.shape[i] - (2 * ghost_layers)
assert flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask # at_border = True
neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell)) #if at_border:
result.append(pdf_index(neighbor_cell_idx, direction_idx, len(self))) # assert flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask
else: # neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell))
raise ValueError("Could not find neighbor for {} direction {} ({})".format(cell, direction, inv_neighbor_cell)) # result.append(pdf_index(neighbor_cell_idx, direction_idx, len(self)))
#else:
raise ValueError("Could not find neighbor for {} direction {} ({})".format(cell, direction, inv_neighbor_cell))
index_array = np.array(result, dtype=np.uint32) index_array = np.array(result, dtype=np.uint32)
index_arr = index_array.reshape([len(stencil) - 1, self.num_fluid_cells]) index_arr = index_array.reshape([len(stencil) - 1, len(self.fluid_coordinates)])
index_arr = index_arr.swapaxes(0, 1) index_arr = index_arr.swapaxes(0, 1)
return index_arr return index_arr
def set_velocity(self, full_vel_arr): #macht aus full_vel_arr einen sparse vel_arr def set_velocity(self, full_vel_arr): #macht aus full_vel_arr einen sparse vel_arr
# Erstmal ghost layers von data_handling herausfinden und wegschneiden... # Erstmal ghost layers von data_handling herausfinden und wegschneiden...
ghost_layers = [(full_vel_arr.shape[i] - self.domain_size[i])//2 for i in range(0, full_vel_arr.shape[-1])] #ghost_layers = [(full_vel_arr.shape[i] - self.domain_size[i])//2 for i in range(0, full_vel_arr.shape[-1])]
full_vel_arr = full_vel_arr[ghost_layers[0]:-ghost_layers[0], ghost_layers[1]:-ghost_layers[1]] #full_vel_arr = full_vel_arr[ghost_layers[0]:-ghost_layers[0], ghost_layers[1]:-ghost_layers[1]]
# Dann initial velocities in sparsen vector schreiben # Dann initial velocities in sparsen vector schreiben
fluid_cell_arr = self.coordinates fluid_cell_arr = self.coordinates
arr = fluid_cell_arr[:self.num_fluid_cells] arr = fluid_cell_arr[:self.num_fluid_cells]
...@@ -188,9 +197,9 @@ class SparseLbPeriodicityMapper: ...@@ -188,9 +197,9 @@ class SparseLbPeriodicityMapper:
self.pdf_field_sparse = pdf_field_sparse self.pdf_field_sparse = pdf_field_sparse
self.index_field = index_field self.index_field = index_field
self.mapping = mapping self.mapping = mapping
self.domain_size = dh.shape
self.periodicity = dh.periodicity self.periodicity = dh.periodicity
self.flag_arr = mapping._flag_arr self.flag_arr = mapping._flag_arr
self.domain_size = self.flag_arr.shape
self.fluid_flag = mapping.fluid_flag self.fluid_flag = mapping.fluid_flag
self.no_slip_flag = mapping.no_slip_flag self.no_slip_flag = mapping.no_slip_flag
self._index_arrays = None self._index_arrays = None
...@@ -219,44 +228,47 @@ class SparseLbPeriodicityMapper: ...@@ -219,44 +228,47 @@ class SparseLbPeriodicityMapper:
def get_read_idx(self, read, write_idx, direction_idx): 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! if self.flag_arr[tuple(read)] & self.no_slip_flag: # Read cell is no-slip: flip PDF!
#read from write cell, inverse direction #read from write cell, inverse direction
#print("Neighbor is solid")
return pdf_index(write_idx, inverse_idx(self.mapping.stencil, direction_idx), len(self.mapping)) return pdf_index(write_idx, inverse_idx(self.mapping.stencil, direction_idx), len(self.mapping))
else: else:
return pdf_index(self.mapping.cell_idx(tuple(read)), direction_idx, len(self.mapping)) return pdf_index(self.mapping.cell_idx(tuple(read)), direction_idx, len(self.mapping))
def get_assignment(self, direction, cell): def get_assignment(self, direction, cell, bool_pos):
cell_idx = self.mapping.cell_idx(cell) cell_idx = self.mapping.cell_idx(cell)
direction_idx = self.direction_index(direction) 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)] #inv_neighbor_cell = [(cell_i - dir_i)%ds_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)]
inv_neighbor_cell = [if b_p: (cell_i + 2*dir_i)%ds_i else: cell_i for cell_i, dir_i, ds_i, b_p in zip(cell, direction, self.domain_size, bool_pos)]
write_idx = pdf_index(cell_idx, direction_idx, len(self.mapping)) write_idx = pdf_index(cell_idx, direction_idx, len(self.mapping))
read_idx = self.get_read_idx(inv_neighbor_cell, cell_idx, direction_idx) read_idx = self.get_read_idx(inv_neighbor_cell, cell_idx, direction_idx)
#print("write:", cell, "read:", inv_neighbor_cell) print("write:", cell, "read:", inv_neighbor_cell)
return [write_idx, read_idx] return [write_idx, read_idx]
def fluid_border(self): def fluid_border(self):
fluid_border_coord = [] fluid_border_coord = []
for cell in self.mapping.fluid_coordinates: 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)] bool_border = [cell_i == 1 or cell_i == ds_i-2 for cell_i, ds_i in zip(cell, self.domain_size)]
if True in bool_border: if True in bool_border:
fluid_border_coord.append(cell) fluid_border_coord.append(cell)
return fluid_border_coord return fluid_border_coord
def create_periodic_index_array(self): def create_periodic_index_array(self):
print(self.domain_size)
stencil = self.mapping.stencil stencil = self.mapping.stencil
fluid_border_coord = self.fluid_border() fluid_border_coord = self.fluid_border()
#print(fluid_border_coord) print(fluid_border_coord)
result = [[[] for j in range(0, len(stencil)-1)] for i in range(0, len(stencil)-1)] 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): for direction_idx, direction in enumerate(stencil):
if all(d_i == 0 for d_i in direction): # Direction (0,0) irrelevant if all(d_i == 0 for d_i in direction): # Direction (0,0) irrelevant
continue continue
#print("\n New direction:", direction, ", ", direction_idx) 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))] periodic_slice_coord = [[int((ds_i-1)*(1+dir_i)/2-dir_i)] if dir_i != 0 else [] for i, (dir_i, ds_i) in enumerate(zip(direction, self.domain_size))]
for cell in fluid_border_coord: 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))] 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) #print(cell, ", ", periodic_slice_coord, ", ", bool_pos)
if True in bool_pos: # This cell & direction is periodical 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)]) # block_direction: where to put this part of index_array (which block to send to...)
result[direction_idx-1][block_direction-1].append(self.get_assignment(direction, cell)) 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, bool_pos))
#print("Goes into", direction_idx-1, block_direction-1) #print("Goes into", direction_idx-1, block_direction-1)
#flatten result: #flatten result:
flattened_result = [] flattened_result = []
......
Source diff could not be displayed: it is too large. Options to address this: view the blob.
Source diff could not be displayed: it is too large. Options to address this: view the blob.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import pytest import pytest
# pytest.importorskip('pycuda') # pytest.importorskip('pycuda')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from lbmpy.session import * from lbmpy.session import *
from lbmpy.sparse import * from lbmpy.sparse import *
from pystencils.field import FieldType from pystencils.field import FieldType
import pystencils as ps import pystencils as ps
import warnings import warnings
from IPython.core.interactiveshell import InteractiveShell from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all" InteractiveShell.ast_node_interactivity = "all"
from pystencils.data_types import TypedSymbol from pystencils.data_types import TypedSymbol
``` ```
%% Output
Traceback (most recent call last):
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3331, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "<ipython-input-2-3785c7fcd9fd>", line 2, in <module>
from lbmpy.sparse import *
File "/home/maja/Studium/Ba/lbmpy/lbmpy/sparse/__init__.py", line 1, in <module>
from .mapping import SparseLbBoundaryMapper, SparseLbMapper, SparseLbPeriodicityMapper
File "/home/maja/Studium/Ba/lbmpy/lbmpy/sparse/mapping.py", line 240
inv_neighbor_cell = [if b_p: (cell_i + 2*dir_i)%ds_i else: cell_i for cell_i, dir_i, ds_i, b_p in zip(cell, direction, self.domain_size, bool_pos)]
^
SyntaxError: invalid syntax
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Sparse (list based) LBM # Sparse (list based) LBM
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
domain_size = (4, 3) domain_size = (3, 2)
omega = 1.8 omega = 1.8
target = 'cpu' target = 'cpu'
stencil = get_stencil("D2Q9") stencil = get_stencil("D2Q9")
ghost_layers = 0 ghost_layers = 1
arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size) arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, True)) dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, True))
src = dh.add_array('src', values_per_cell=len(stencil)) src = dh.add_array('src', values_per_cell=len(stencil))
dh.fill('src', 0.0, ghost_layers=True) dh.fill('src', 0.0, ghost_layers=True)
dst = dh.add_array('dst', values_per_cell=len(stencil)) dst = dh.add_array('dst', values_per_cell=len(stencil))
dh.fill('dst', 0.0, ghost_layers=True) dh.fill('dst', 0.0, ghost_layers=True)
velField = dh.add_array('velField', values_per_cell=dh.dim) velField = dh.add_array('velField', values_per_cell=dh.dim)
dh.fill('velField', 0.0, ghost_layers=True) dh.fill('velField', 0.0, ghost_layers=True)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
ps.stencil.plot_2d(stencil) ps.stencil.plot_2d(stencil)
``` ```
%% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
noslip = NoSlip() noslip = NoSlip()
flags = { flags = {
'fluid': 1, 'fluid': 1,
noslip: 2, noslip: 2,
} }
flag_arr = np.zeros(arr_size, dtype=np.uint16) flag_arr = np.zeros(arr_size, dtype=np.uint16)
flag_arr.fill(flags['fluid']) flag_arr.fill(flags['fluid'])
plt.scalar_field(flag_arr) plt.scalar_field(flag_arr)
plt.colorbar(); plt.colorbar();
``` ```
%% Output
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Mappings ### Mappings
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mapping = SparseLbMapper(stencil, domain_size, flag_arr, flags['fluid'], flags[noslip], 0) #Warum müssen (dürfen!) hier Nullen stehen? mapping = SparseLbMapper(stencil, domain_size, flag_arr, flags['fluid'], flags[noslip], 0) #Warum müssen (dürfen!) hier Nullen stehen?
index_arr = mapping.create_index_array(ghost_layers) index_arr = mapping.create_index_array(ghost_layers)
# Arrays # Arrays
#index_arr = index_arr_linear.reshape([len(method.stencil), mapping.num_fluid_cells]) #index_arr = index_arr_linear.reshape([len(method.stencil), mapping.num_fluid_cells])
#index_arr = index_arr.swapaxes(0, 1) #index_arr = index_arr.swapaxes(0, 1)
pdf_arr = np.empty((len(mapping), len(stencil)), order='f') pdf_arr = np.empty((len(mapping), len(stencil)), order='f')
pdf_arr_tmp = np.empty_like(pdf_arr) pdf_arr_tmp = np.empty_like(pdf_arr)
vel_arr = np.ones([mapping.num_fluid_cells, len(stencil[0])], order='f') vel_arr = np.ones([mapping.num_fluid_cells, len(stencil[0])], order='f')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
shape_pdf_arr = pdf_arr.shape shape_pdf_arr = pdf_arr.shape
num = 0 num = 0
for x in range(shape_pdf_arr[0]): for x in range(shape_pdf_arr[0]):
for y in range(shape_pdf_arr[1]): for y in range(shape_pdf_arr[1]):
pdf_arr[x, y] = num pdf_arr[x, y] = num
num += 1 num += 1
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]") pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]")
pdf_field.field_type = FieldType.CUSTOM pdf_field.field_type = FieldType.CUSTOM
pdf_field.pdf_field_tmp = FieldType.CUSTOM pdf_field.pdf_field_tmp = FieldType.CUSTOM
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
periodicity_mapper = SparseLbPeriodicityMapper(mapping, dh, pdf_field) periodicity_mapper = SparseLbPeriodicityMapper(mapping, dh, pdf_field)
periodicity_mapper.create_periodic_index_array()
def handle_periodicity(): def handle_periodicity():
periodicity_mapper(pdf_arr) periodicity_mapper(pdf_arr)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
pdf_arr pdf_arr
``` ```
%% Output
array([[ 0., 1., 2., 3., 4., 5., 6., 7., 8.],
[ 9., 10., 11., 12., 13., 14., 15., 16., 17.],
[ 18., 19., 20., 21., 22., 23., 24., 25., 26.],
[ 27., 28., 29., 30., 31., 32., 33., 34., 35.],
[ 36., 37., 38., 39., 40., 41., 42., 43., 44.],
[ 45., 46., 47., 48., 49., 50., 51., 52., 53.],
[ 54., 55., 56., 57., 58., 59., 60., 61., 62.],
[ 63., 64., 65., 66., 67., 68., 69., 70., 71.],
[ 72., 73., 74., 75., 76., 77., 78., 79., 80.],
[ 81., 82., 83., 84., 85., 86., 87., 88., 89.],
[ 90., 91., 92., 93., 94., 95., 96., 97., 98.],
[ 99., 100., 101., 102., 103., 104., 105., 106., 107.]])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
handle_periodicity() handle_periodicity()
pdf_arr pdf_arr
``` ```
%% Output
array([[ 0., 19., 2., 3., 85., 50., 105., 7., 98.],
[ 9., 10., 11., 12., 94., 14., 78., 16., 62.],
[ 18., 19., 2., 21., 103., 23., 96., 34., 89.],
[ 27., 46., 29., 30., 31., 77., 96., 34., 35.],
[ 36., 37., 38., 39., 40., 41., 42., 43., 44.],
[ 45., 46., 29., 48., 49., 50., 51., 61., 98.],
[ 54., 73., 56., 57., 58., 14., 51., 61., 62.],
[ 63., 64., 65., 66., 67., 68., 69., 70., 71.],
[ 72., 73., 56., 75., 76., 77., 78., 16., 35.],
[ 81., 100., 83., 3., 85., 23., 78., 16., 89.],
[ 90., 91., 92., 12., 94., 50., 96., 34., 98.],
[ 99., 100., 83., 21., 103., 14., 105., 7., 62.]])
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment