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

Call function for periodicity kernel.

parent a993a774
No related branches found
No related tags found
No related merge requests found
Pipeline #34381 failed
Stage: test
Expand all Hide whitespace changes
Inline Side-by-side
lbmpy/sparse/mapping.py
+ 28
− 20
View file @ 30c7ae00
......@@ -31,15 +31,14 @@ class SparseLbMapper:
flag_arr: integer array where each bit corresponds to a boundary or 'fluid'
"""
def __init__(self, stencil, domain_size, flag_arr, fluid_flag, no_slip_flag, ubb_flag):
def __init__(self, stencil, domain_size, flag_arr, fluid_flag, no_slip_flag, other_boundary_mask):
self._flag_arr = flag_arr
self._coordinate_arr = None
self._sorter = None # array of indices that sort _coordinate_arr
self._dirty = True
self.fluid_flag = fluid_flag
self.no_slip_flag = no_slip_flag
self.ubb_flag = ubb_flag
self.boundary_mask = ubb_flag
self.other_boundary_mask = other_boundary_mask
self._num_fluid_cells = None
self.stencil = stencil
self.domain_size = domain_size
......@@ -108,7 +107,7 @@ class SparseLbMapper:
# Add fluid cells
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.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]
total_cells = len(coordinates_fluid) + len(coordinates_boundary)
......@@ -126,7 +125,7 @@ class SparseLbMapper:
stencil = self.stencil
flag_arr = self.flag_array
no_slip_flag = self.no_slip_flag
fluid_boundary_mask = self.ubb_flag | self.fluid_flag
fluid_boundary_mask = self.other_boundary_mask | self.fluid_flag
result = []
for direction_idx, direction in enumerate(stencil):
#print("dir:", direction_idx)
......@@ -134,8 +133,7 @@ class SparseLbMapper:
assert direction_idx == 0
continue
for own_cell_idx, cell in enumerate(self.fluid_coordinates):
domain_size = (len(self.flag_array), len(self.flag_array[0]))
inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i, ds_i in zip(cell, direction, self.domain_size)])
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)]
if (True in at_border):
inv_neighbor_cell = list(cell)
......@@ -146,11 +144,10 @@ class SparseLbMapper:
elif flag_arr[tuple(inv_neighbor_cell)] & no_slip_flag: # no-slip before periodicity!
result.append(pdf_index(own_cell_idx, inverse_idx(stencil, direction_idx), len(self)))
else:
#periodicity handler
print("third")
#This is where we end up if inv_neighbor_cell is not fluid/no_slip/at the border;
at_border = False
for i, x_i in enumerate(inv_neighbor_cell):
if x_i == (ghost_layers - 1): #if inv_neighbor_cell NOT noslip/fuid & one coordinate 0
if x_i == (ghost_layers - 1):
inv_neighbor_cell[i] += flag_arr.shape[i] - (2 * ghost_layers)
at_border = True
elif x_i == flag_arr.shape[i] - ghost_layers:
......@@ -161,7 +158,7 @@ class SparseLbMapper:
neighbor_cell_idx = self.cell_idx(tuple(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))
raise ValueError("Could not find neighbor for {} direction {} ({})".format(cell, direction, inv_neighbor_cell))
index_array = np.array(result, dtype=np.uint32)
index_arr = index_array.reshape([len(stencil) - 1, self.num_fluid_cells])
......@@ -183,10 +180,21 @@ class SparseLbPeriodicityMapper:
self.flag_arr = mapping._flag_arr
self.fluid_flag = mapping.fluid_flag
self.no_slip_flag = mapping.no_slip_flag
self._index_arrays = None
self._kernel = None
self._dirty = True
def __call__(self, f):
if self._dirty:
self.create_periodic_index_array()
self.create_kernel()
for i_a in self._index_arrays:
self._kernel(f=f, idx=i_a)
def __call__(self, f, idx_arrays, kernel):
for i_a in idx_arrays:
kernel(f=f, idx=i_a)
@property
def get_kernel(self):
return self._kernel
def direction_index(self, direction):
direction = tuple(direction)
......@@ -219,7 +227,7 @@ class SparseLbPeriodicityMapper:
fluid_border_coord.append(cell)
return fluid_border_coord
def create_periodic_index_array(self):
def create_periodic_index_array(self):
stencil = self.mapping.stencil
fluid_border_coord = self.fluid_border()
#print(fluid_border_coord)
......@@ -237,10 +245,6 @@ class SparseLbPeriodicityMapper:
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])
#flatten result:
flattened_result = []
for block in result:
......@@ -248,6 +252,8 @@ class SparseLbPeriodicityMapper:
if (i_a):
flattened_result.append(np.array(i_a))
flattened_result = np.array(flattened_result)
self._index_arrays = flattened_result
self._dirty = False
return flattened_result
def assignments(self):
......@@ -255,8 +261,10 @@ class SparseLbPeriodicityMapper:
read = self.pdf_field_sparse.absolute_access((self.index_field(1),), ())
return [Assignment(self.DIR_SYMBOL, self.index_field(0)), Assignment(write, read)]
def get_kernel(self):
def create_kernel(self):
kernel = create_kernel(self.assignments(), ghost_layers=0)
_ast = kernel
self._kernel = _ast.compile()
return kernel
class SparseLbBoundaryMapper:
......
lbmpy_tests/simple_sparse_shear_flow.ipynb
+ 71
− 38
View file @ 30c7ae00
Source diff could not be displayed: it is too large. Options to address this: view the blob.
lbmpy_tests/test_sparse_lbm.ipynb
+ 172
− 67
View file @ 30c7ae00
This diff is collapsed.
lbmpy_tests/test_sparse_lbm_peridocity_only.ipynb
+ 11
− 14
View file @ 30c7ae00
%% Cell type:code id: tags:
``` python
import pytest
# pytest.importorskip('pycuda')
```
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.sparse import *
from pystencils.field import FieldType
import pystencils as ps
import warnings
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
from pystencils.data_types import TypedSymbol
```
%% Cell type:markdown id: tags:
# Sparse (list based) LBM
%% Cell type:code id: tags:
``` python
domain_size = (4, 3)
omega = 1.8
target = 'cpu'
stencil = get_stencil("D2Q9")
ghost_layers = 0
arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
```
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, True))
src = dh.add_array('src', values_per_cell=len(stencil))
dh.fill('src', 0.0, ghost_layers=True)
dst = dh.add_array('dst', values_per_cell=len(stencil))
dh.fill('dst', 0.0, ghost_layers=True)
velField = dh.add_array('velField', values_per_cell=dh.dim)
dh.fill('velField', 0.0, ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
ps.stencil.plot_2d(stencil)
```
%% Output
%% Cell type:code id: tags:
``` python
noslip = NoSlip()
flags = {
'fluid': 1,
noslip: 2,
}
flag_arr = np.zeros(arr_size, dtype=np.uint16)
flag_arr.fill(flags['fluid'])
```
%% Cell type:markdown id: tags:
### Mappings
%% Cell type:code id: tags:
``` python
mapping = SparseLbMapper(stencil, domain_size, flag_arr, flags['fluid'], flags[noslip], 4) #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)
# Arrays
#index_arr = index_arr_linear.reshape([len(method.stencil), mapping.num_fluid_cells])
#index_arr = index_arr.swapaxes(0, 1)
pdf_arr = np.empty((len(mapping), len(stencil)), order='f')
pdf_arr_tmp = np.empty_like(pdf_arr)
vel_arr = np.ones([mapping.num_fluid_cells, len(stencil[0])], order='f')
```
%% Cell type:code id: tags:
``` python
shape_pdf_arr = pdf_arr.shape
num = 0
for x in range(shape_pdf_arr[0]):
for y in range(shape_pdf_arr[1]):
pdf_arr[x, y] = num
num += 1
```
%% Cell type:code id: tags:
``` python
pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]")
pdf_field.field_type = FieldType.CUSTOM
pdf_field.pdf_field_tmp = FieldType.CUSTOM
```
%% Cell type:code id: tags:
``` python
periodicity_mapper = SparseLbPeriodicityMapper(mapping, dh, pdf_field)
periodicity_ast = periodicity_mapper.get_kernel()
periodicity_kernel = periodicity_ast.compile()
periodicity_index_arrays = periodicity_mapper.create_periodic_index_array()
def handle_periodicity():
periodicity_mapper(pdf_arr, periodicity_index_arrays, periodicity_kernel)
periodicity_mapper(pdf_arr)
```
%% Cell type:code id: tags:
``` python
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:
``` python
handle_periodicity()
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., 24., 34., 35.],
[ 36., 37., 38., 39., 40., 41., 42., 43., 44.],
[ 45., 46., 29., 48., 49., 50., 51., 61., 8.],
[ 54., 73., 56., 57., 58., 104., 51., 61., 62.],
[ 63., 64., 65., 66., 67., 68., 69., 70., 71.],
[ 72., 73., 56., 75., 76., 77., 78., 88., 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:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment