.gitignore

@@ -16,4 +16,5 @@ _local_tmp
 doc/bibtex.json
 /db
 /lbmpy/phasefield/simplex_projection.*.so
-/lbmpy/phasefield/simplex_projection.c
\ No newline at end of file
+/lbmpy/phasefield/simplex_projection.c
+*.swp

lbmpy/boundaries/boundaryconditions.py

@@ -173,12 +173,11 @@ class UBB(LbBoundary):
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         vel_from_idx_field = callable(self._velocity)
         vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
-        direction = dir_symbol
 
         assert self.dim == lb_method.dim, \
             f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"
 
-        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(direction, lb_method.stencil)
+        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
 
         velocity = tuple(v_i.get_shifted(*neighbor_offset)
                          if isinstance(v_i, Field.Access) and not vel_from_idx_field
@@ -193,7 +192,7 @@ class UBB(LbBoundary):
         c_s_sq = sp.Rational(1, 3)
         weight_of_direction = LbmWeightInfo.weight_of_direction
         vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
-            direction, lb_method)
+            dir_symbol, lb_method)
 
         # Better alternative: in conserved value computation
         # rename what is currently called density to "virtual_density"
@@ -205,12 +204,12 @@ class UBB(LbBoundary):
             density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
             density_symbol = lb_method.conserved_quantity_computation.defined_symbols()['density']
             result = density_equations.all_assignments
-            result += [Assignment(f_in(inv_dir[direction]),
-                                  f_out(direction) - vel_term * density_symbol)]
+            result += [Assignment(f_in(inv_dir[dir_symbol]),
+                                  f_out(dir_symbol) - vel_term * density_symbol)]
             return result
         else:
-            return [Assignment(f_in(inv_dir[direction]),
-                               f_out(direction) - vel_term)]
+            return [Assignment(f_in(inv_dir[dir_symbol]),
+                               f_out(dir_symbol) - vel_term)]
 
 # end class UBB
 

lbmpy/sparse/__init__.py

-from .mapping import SparseLbBoundaryMapper, SparseLbMapper
+from .mapping import SparseLbBoundaryMapper, SparseLbMapper, SparseLbPeriodicityMapper, SparseLbCommunicationMapper
 from .update_rule_sparse import (
     create_lb_update_rule_sparse, create_macroscopic_value_getter_sparse,
     create_macroscopic_value_setter_sparse, create_symbolic_list)
 
-__all__ = ['SparseLbBoundaryMapper', 'SparseLbMapper', 'create_lb_update_rule_sparse',
+__all__ = ['SparseLbBoundaryMapper', 'SparseLbPeriodicityMapper', 'SparseLbMapper', 'SparseLbCommunicationMapper', 'create_lb_update_rule_sparse',
            'create_macroscopic_value_setter_sparse', 'create_macroscopic_value_getter_sparse',
            'create_symbolic_list']

lbmpy/sparse/mapping.py

 from typing import Tuple
+import operator
 
 import numpy as np
 import sympy as sp
 
+from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing, NeighbourOffsetArrays
+from lbmpy.advanced_streaming.utility import Timestep
 from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
-from pystencils import Assignment, Field, TypedSymbol
+from pystencils import Assignment, Field, TypedSymbol, create_kernel
 from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
 from pystencils.boundaries.createindexlist import (
     boundary_index_array_coordinate_names, create_boundary_index_list, direction_member_name)
 
+from pystencils.field import FieldType
+from pystencils import Assignment, AssignmentCollection
 
+
+def pdf_index(cell_index, direction_index, domain_size):
+    return 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.
 
     Args:
           flag_arr: integer array where each bit corresponds to a boundary or 'fluid'
     """
-    def __init__(self, stencil, flag_arr, fluid_flag, no_slip_flag, other_boundary_mask):
+
+    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
@@ -25,7 +40,11 @@ class SparseLbMapper:
         self.no_slip_flag = no_slip_flag
         self.other_boundary_mask = other_boundary_mask
         self._num_fluid_cells = None
+        self._num_ghost_cells = None
         self.stencil = stencil
+        self.domain_size = domain_size
+        self.ghost_layers = None
+        
 
     @property
     def coordinates(self):
@@ -34,20 +53,49 @@ class SparseLbMapper:
 
         return self._coordinate_arr
 
+    @property
+    def flag_array(self):
+        if self._dirty:
+            self._assemble()
+        return self._flag_arr
+    
+    
     @property
     def fluid_coordinates(self):
         if self._dirty:
             self._assemble()
+        #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
     def num_fluid_cells(self):
         return self._num_fluid_cells
-
-    @property
-    def flag_array(self):
-        return self._flag_arr
-
+    
+    def _at_border(self,cell):
+        dim = self._flag_arr.shape
+        if True in [cell_i == 0 or cell_i == ds_i-1 for ds_i, cell_i in zip(dim, cell)]:
+            return [0]
+        border = [cell_i == 1 or cell_i == ds_i-2 for ds_i, cell_i in zip(dim, cell)]
+        return border    
+    
+    def _amend_flag_arr(self):
+        # Fix Flag field jetzt neu: für n-dim zum Sparpreis!
+        coordinates_no_slip = np.argwhere(np.bitwise_and(self._flag_arr, self.no_slip_flag)).astype(np.uint32)
+        res = []
+        dim = self._flag_arr.shape
+        for cell in coordinates_no_slip:
+            border = self._at_border(cell)
+            if True in border:
+                direction = [-1 if cell_i == 1 else 1 if cell_i == ds_i-2 else 0 for ds_i, cell_i in zip(dim, cell)]
+                l = [(cell_i + sl_i*2)%ds_i for cell_i, bb_i, ds_i, sl_i in zip(cell, border, dim, direction)]
+                pos = np.array([cell, l])
+                for i in range(pow(2, len(l))):
+                    i_bin = (format(i, '0'+str(len(l))+'b'))
+                    opp = [pos[int(i_bin[arr]), arr] for arr in range(len(l))]
+                    self._flag_arr[tuple(opp)] = self.no_slip_flag
+    
     def cell_idx(self, coordinate: Tuple[int, ...]) -> np.uint32:
         """Maps from coordinates (x,y,z) or (x,y) tuple to the list index. Raises ValueError if coordinate not found."""
         if self._dirty:
@@ -87,12 +135,15 @@ class SparseLbMapper:
     def _assemble(self):
         dim = len(self._flag_arr.shape)
         struct_type = np.dtype([(name, np.uint32) for name in boundary_index_array_coordinate_names[:dim]])
+        
+        self._amend_flag_arr()
 
         # 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.other_boundary_mask)).astype(np.uint32)
         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)
 
         self._coordinate_arr = np.empty((total_cells,), dtype=struct_type)
@@ -105,103 +156,330 @@ class SparseLbMapper:
 
     def create_index_array(self, ghost_layers=1):
         # TODO support different layouts here
+        self.ghost_layers = ghost_layers
         stencil = self.stencil
-
-        def pdf_index(cell_index, direction_index):
-            return cell_index + direction_index * len(self)
-
-        def inverse_idx(idx):
-            return stencil.index(tuple(-d_i for d_i in stencil[idx]))
-
         flag_arr = self.flag_array
+        #print(flag_arr)
         no_slip_flag = self.no_slip_flag
         fluid_boundary_mask = self.other_boundary_mask | self.fluid_flag
-
         result = []
         for direction_idx, direction in enumerate(stencil):
+            #print("dir:", direction_idx)
             if all(d_i == 0 for d_i in direction):
                 assert direction_idx == 0
                 continue
             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)])
-                if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask:
+                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)]
+                #print("write:", cell, "read:", inv_neighbor_cell)
+                if (True in at_border):
+                    #print("at border")
+                    result.append(pdf_index(own_cell_idx, direction_idx, len(self)))
+                elif 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))
+                    result.append(pdf_index(neighbor_cell_idx, direction_idx, len(self)))
+                    #print("write:", cell, "read:", inv_neighbor_cell, "direction:", direction_idx)
                 elif flag_arr[tuple(inv_neighbor_cell)] & no_slip_flag:  # no-slip before periodicity!
-                    result.append(pdf_index(own_cell_idx, inverse_idx(direction_idx)))
+                    #print("no_slip")
+                    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))
                 else:
-                    at_border = False
-                    for i, x_i in enumerate(inv_neighbor_cell):
-                        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:
-                            inv_neighbor_cell[i] -= flag_arr.shape[i] - (2 * ghost_layers)
-                            at_border = True
-                    if at_border:
-                        assert 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))
-                    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])
-        index_arr = index_arr.swapaxes(0, 1)
+        index_arr = index_array.reshape([len(stencil) - 1, len(self.fluid_coordinates)])
+        index_arr = index_arr.swapaxes(0, 1) 
         return index_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...
+        #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]]
+        # Dann initial velocities in sparsen vector schreiben
+        fluid_cell_arr = self.coordinates
+        arr = fluid_cell_arr[:self.num_fluid_cells]
+        if (len(self._flag_arr.shape) == 3):
+            sparse_vel_arr = full_vel_arr[arr['x'], arr['y'], arr['z']]
+        elif (len(self._flag_arr.shape) == 2):
+            sparse_vel_arr = full_vel_arr[arr['x'], arr['y']]
+        return sparse_vel_arr
+
+    
+class SparseLbPeriodicityMapper:
+    DIR_SYMBOL = TypedSymbol("dir", np.uint32)
 
-
+    def __init__(self, mapping: SparseLbMapper, dh, pdf_field_sparse):
+        
+        index_field = Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=np.uint32)
+        self.pdf_field_sparse = pdf_field_sparse
+        self.index_field = index_field
+        self.mapping = mapping
+        self.periodicity = dh.periodicity
+        self.flag_arr = mapping._flag_arr
+        self.domain_size = self.flag_arr.shape
+        self.no_slip_flag = mapping.no_slip_flag
+        self._ghost_cells = None
+        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)
+        self._kernel(f=f, idx=self._index_arrays)
+    
+    @property
+    def get_kernel(self):
+        return self._kernel
+    
+    def at_border(self, cell):
+        return True in [cell_i == 0 or cell_i  == ds_i-1 for cell_i, ds_i in zip(cell, self.domain_size)]
+    
+    def ghost_cells(self):
+        border_bool = [self.at_border(cell) for cell in self.mapping._coordinate_arr[:self.mapping.num_fluid_cells]]
+        self._ghost_cells = self.mapping._coordinate_arr[:self.mapping.num_fluid_cells][border_bool]
+        return
+    
+    def get_assignment(self, direction, cell, bool_slice):
+        cell_idx = self.mapping.cell_idx(cell)
+        direction_idx = self.mapping.stencil.index(direction)
+        pdf_read_idx = pdf_index(self.mapping.cell_idx(tuple(read)), direction_idx, len(self.mapping))
+        pdf_write_idx = pdf_index(cell_idx, direction_idx, len(self.mapping))
+        return [pdf_write_idx, pdf_read_idx]
+    
+    def create_periodic_index_array(self):
+        self.ghost_cells()
+        stencil = self.mapping.stencil
+        result = [[[] for j in range(0, len(stencil)-1)] for i in range(0, len(stencil)-1)]
+        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)
+            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 i, cell in enumerate(self._ghost_cells): # Only fluid ghost cells
+                bool_slice = [cell_i in slice_i for i, (cell_i, slice_i) in enumerate(zip(cell, periodic_slice_coord))]
+                periodic_coord_range = [[1 - dir_i, ds_i - 2 - dir_i] for dir_i, ds_i in zip(direction, self.domain_size)]
+                bool_range = [pcr_i[0] <= cell_i <= pcr_i[1] for pcr_i, cell_i in zip(periodic_coord_range, cell)]
+                #print(cell, ", range", periodic_coord_range, "bool", bool_range)
+                if all(bool_range): # This ghost cell needs this direction value
+                    # block_direction: where to put this part of index_array (which block to send to...)
+                    block_direction = [int(bp_i)*dir_i for dir_i, bp_i in zip(direction, bool_slice)]
+                    block_index = self.mapping.stencil.index(tuple(block_direction))
+                    result[block_index-1][direction_idx-1].append(self.get_assignment(direction, cell, bool_slice))
+                    #print("Goes into", direction_idx-1, block_index-1)
+        # Flatten result array:
+        flattened_result = []
+        #print(result)
+        for block in result:
+            for i_a in block:
+                if (i_a):
+                    for pair in i_a:
+                        flattened_result.append(np.array(pair, dtype=np.uint32))
+        flattened_result = np.array(flattened_result, dtype=np.uint32)
+        self._index_arrays = flattened_result
+        self._dirty = False
+        #print(flattened_result)
+        return flattened_result
+    
+    def assignments(self):
+        write = self.pdf_field_sparse.absolute_access((self.index_field(0),), ())
+        read = self.pdf_field_sparse.absolute_access((self.index_field(1),), ())
+        return [Assignment(self.DIR_SYMBOL, self.index_field(0)), Assignment(write, read)]
+
+    def create_kernel(self):
+        kernel = create_kernel(self.assignments(), ghost_layers=0)
+        _ast = kernel
+        self._kernel = _ast.compile()
+        return kernel
+
+class SparseLbCommunicationMapper:
+    
+    def __init__(self, mapping: SparseLbMapper, dh, pdf_field_sparse):
+        index_field = Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=np.int64)
+        self.pdf_field_sparse = pdf_field_sparse
+        self.index_field = index_field
+        self.mapping = mapping
+        self.flag_arr = mapping._flag_arr
+        self.domain_size = self.flag_arr.shape
+        self.no_slip_flag = mapping.no_slip_flag
+        #self._ghost_cells = None
+        #self._border_cells = None
+        self._parallel_index_array = None
+        self._bounce_back_index_array = None
+        self._send_packages = None
+        self._receive_here = None
+        self._received_packages = None
+        self._kernel = None
+        self._dirty = True
+        
+    def _ghost(self, cell):
+        return True in [cell_i == 0 or cell_i  == ds_i-1 for cell_i, ds_i in zip(cell, self.domain_size)]
+    
+    #def _ghost_cells(self):
+    #    return [cell if self._ghost(cell) for cell in self._coordinates_all]
+    
+    @property
+    def ghost_cells(self):
+        ghost_bool = [self._ghost(cell) for cell in self.mapping._coordinate_arr[:self.mapping.num_fluid_cells]]
+        #self._ghost_cells = 
+        return self.mapping._coordinate_arr[:self.mapping.num_fluid_cells][ghost_bool]
+    
+    def _at_border_not_ghost(self, cell):
+        if self._ghost(cell):
+            return False
+        return True in [cell_i == 1 or cell_i  == ds_i-2 for cell_i, ds_i in zip(cell, self.domain_size)]
+    @property
+    def border_cells(self):
+        border_bool = [self._at_border_not_ghost(cell) for cell in self.mapping._coordinate_arr[:self.mapping.num_fluid_cells]]
+        #self._border_cells 
+        return self.mapping._coordinate_arr[:self.mapping.num_fluid_cells][border_bool]
+    
+    #def _get_assignment(self, direction, cell, neighbor, bool_slice): #Falls man auf solide Zellen achten muss
+    #    cell = tuple(cell)
+    #    #print("pack:", cell, "neighbor:", tuple(neighbor))
+    #    if (self.flag_arr[cell] | self.flag_arr[tuple(neighbor)]) & self.no_slip_flag :
+    #        #print("is solid")
+    #        return -1
+    #    cell_idx = self.mapping.cell_idx(cell)
+    #    direction_idx = self.mapping.stencil.index(direction)
+    #    pdf_cell_idx = pdf_index(cell_idx, direction_idx, len(self.mapping))
+    #    return pdf_cell_idx
+    
+    def create_packages(self):
+        stencil = self.mapping.stencil
+        result = [[[] for j in range(0, len(stencil)-1)] for i in range(0, len(stencil)-1)]
+        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)
+            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 self.border_cells:
+                bool_slice = [cell_i in slice_i for (cell_i, slice_i) in zip(cell, periodic_slice_coord)]
+                if True in bool_slice: # This ghost cell needs this direction value
+                    # block_direction: where to put this part of index_array (which block to send to...)
+                    block_direction = [int(bp_i)*dir_i for dir_i, bp_i in zip(direction, bool_slice)]
+                    block_index = self.mapping.stencil.index(tuple(block_direction))
+                    #neighbor_ghost_cell = [(cell_i + dir_i*int(bs_i)) for cell_i, dir_i, bs_i in zip(cell, direction, bool_slice)]
+                    #result[block_index-1][direction_idx-1].append(self._get_assignment(direction, cell, neighbor_ghost_cell, bool_slice))
+                    cell_idx = self.mapping.cell_idx(cell)
+                    #print(cell, cell_idx)
+                    result[block_index-1][direction_idx-1].append(pdf_index(cell_idx, direction_idx, len(self.mapping)))
+                    #print("Goes into", block_index-1, direction_idx-1)
+        self._send_packages = result
+        return result
+    
+    def receive_here(self):
+        self._received_packages = self._send_packages #tmp
+        stencil = self.mapping.stencil
+        result = [[[] for j in range(0, len(stencil)-1)] for i in range(0, len(stencil)-1)]
+        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)
+            periodic_slice_coord = [[int((ds_i-1)*(1-dir_i)/2)] if dir_i != 0 else [] for (dir_i, ds_i) in zip(direction, self.domain_size)]          
+            for cell in self.ghost_cells: # Fluid and solid ghost cells
+                bool_slice = [cell_i in slice_i for (cell_i, slice_i) in zip(cell, periodic_slice_coord)]
+                periodic_coord_range = [[1 - dir_i, ds_i - 2 - dir_i] for dir_i, ds_i in zip(direction, self.domain_size)]
+                bool_range = [pcr_i[0] <= cell_i <= pcr_i[1] for pcr_i, cell_i in zip(periodic_coord_range, cell)]
+                #print(cell, ", range", periodic_coord_range, "bool", bool_range)
+                if all(bool_range): # This ghost cell needs this direction value
+                    # block_direction: where to put this part of index_array (which block to send to...)
+                    block_direction = [int(bp_i)*dir_i for dir_i, bp_i in zip(direction, bool_slice)]
+                    block_index = self.mapping.stencil.index(tuple(block_direction))
+                    #future_pull_cell = [(cell_i + dir_i) for cell_i, dir_i, bs_i in zip(cell, direction, bool_slice)]
+                    cell_idx = self.mapping.cell_idx(cell)
+                    cell_idx = self.mapping.cell_idx(cell)
+                    #print(cell, cell_idx)
+                    result[block_index-1][direction_idx-1].append(pdf_index(cell_idx, direction_idx, len(self.mapping)))
+                    #result[block_index-1][direction_idx-1].append(self._get_assignment_rec(direction, cell, future_pull_cell, bool_slice))
+                    #print("Goes into", block_index-1, direction_idx-1)
+        self._receive_here = result
+        return result
+    
+    def create_index_array(self):
+        #During initialization, somehow receive packages ... 
+        self._received_packages = self._send_packages #tmp
+        self._parallel_index_array = []
+         
+            
 class SparseLbBoundaryMapper:
     NEIGHBOR_IDX_NAME = 'nidx{}'
-    DIR_SYMBOL = TypedSymbol("dir", np.uint32)
+    DIR_SYMBOL = TypedSymbol("dir", np.int64)
 
     def __init__(self, boundary, method, pdf_field_sparse):
         full_pdf_field = Field.create_generic('pdfFull', spatial_dimensions=method.dim, index_dimensions=1)
+        full_pdf_field.field_type = FieldType.CUSTOM
+        prev_timestep = Timestep.BOTH
+        streaming_pattern = 'pull'
 
         additional_data_field = Field.create_generic('additionalData', spatial_dimensions=1,
                                                      dtype=boundary.additional_data)
-        boundary_eqs = boundary(full_pdf_field, self.DIR_SYMBOL, method, additional_data_field)
+
+        indexing = BetweenTimestepsIndexing(
+            full_pdf_field, method.stencil, prev_timestep, streaming_pattern, np.int64, np.int64)
+
+        f_out, f_in = indexing.proxy_fields
+        dir_symbol = indexing.dir_symbol
+        inv_dir = indexing.inverse_dir_symbol
+
+        boundary_eqs = boundary(f_out, f_in, dir_symbol, inv_dir, method, additional_data_field)
+        #print("Boundary:", boundary_eqs)
+
         neighbor_offsets = {fa.offsets for eq in boundary_eqs for fa in eq.atoms(Field.Access)}
         neighbor_offsets = list(neighbor_offsets)
+        #print("Neighbor_offsets", neighbor_offsets)
 
-        neighbor_offsets_dtype = [(self.NEIGHBOR_IDX_NAME.format(i), np.uint32)
-                                  for i in range(len(neighbor_offsets))]
+        neighbor_offsets_dtype = [(self.NEIGHBOR_IDX_NAME.format(i), np.int64)
+                                  for i in range(method.dim)]
 
-        index_field_dtype = np.dtype([('dir', np.uint32),
+        index_field_dtype = np.dtype([('dir', np.int64),
                                       *neighbor_offsets_dtype,
                                       *boundary.additional_data])
-        index_field = Field.create_generic('indexField', spatial_dimensions=1, dtype=index_field_dtype)
-        boundary_eqs = boundary(full_pdf_field, self.DIR_SYMBOL, method, index_field)
-
-        offset_subs = {off: sp.Symbol(self.NEIGHBOR_IDX_NAME.format(i)) for i, off in enumerate(neighbor_offsets)}
-
-        new_boundary_eqs = []
-        for eq in boundary_eqs:
-            substitutions = {
-                fa: pdf_field_sparse.absolute_access([index_field(offset_subs[fa.offsets].name)], fa.index)
-                for fa in eq.atoms(Field.Access)
-                if fa.field == full_pdf_field
-            }
-            new_boundary_eqs.append(eq.subs(substitutions))
-
-        self.boundary_eqs = new_boundary_eqs
-        self.boundary_eqs_orig = boundary_eqs
+        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.boundary_eqs = assignments.new_with_substitutions(accessor_subs)
+        self.boundary_eqs_orig = assignments
         self.method = method
         self.index_field_dtype = index_field_dtype
         self.neighbor_offsets = neighbor_offsets
         self.index_field = index_field
+        self.timestep_indexing = indexing
+        self.boundary_functor = boundary
 
     def _build_substitutions(self):
         dim = self.method.dim
         stencil = self.method.stencil
+        neighbour_offset = NeighbourOffsetArrays(stencil)
 
-        result = [{BoundaryOffsetInfo.offset_from_dir(self.DIR_SYMBOL, dim)[current_dim]: offset[current_dim]
+        result = [{self.neighbor_offsets[0][current_dim]: offset[current_dim]
                    for current_dim in range(dim)}
                   for i, offset in enumerate(self.method.stencil)
                   ]
+
         for dir_idx, subs_dict in enumerate(result):
             inv_idx = stencil.index(tuple(-e for e in stencil[dir_idx]))
-            subs_dict[BoundaryOffsetInfo.inv_dir(self.DIR_SYMBOL)] = 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):
@@ -211,28 +489,37 @@ class SparseLbBoundaryMapper:
                                              flag_dtype(boundary_mask), flag_dtype(mapping.fluid_flag),
                                              nr_of_ghost_layers)
 
-        result = np.empty(idx_arr.shape, dtype=self.index_field_dtype)
+        result = []
 
         dim = self.method.dim
 
         coord_names = boundary_index_array_coordinate_names[:dim]
         center_coordinates = idx_arr[coord_names]
 
-        substitutions = self._build_substitutions()
-        for j, neighbor_offset in enumerate(self.neighbor_offsets):
-            neighbor_coordinates = center_coordinates.copy()
-            offsets = np.array([tuple(int(sp.sympify(e).subs(substitution)) for e in neighbor_offset)
-                                for substitution in substitutions])
-            for i, coord_name in enumerate(coord_names):
-                neighbor_coordinates[coord_name] += offsets[:, i][idx_arr['dir']]
-            result[self.NEIGHBOR_IDX_NAME.format(j)] = mapping.cell_idx_bulk(neighbor_coordinates)
-
-        result[direction_member_name] = idx_arr[direction_member_name]
-        # Find neighbor indices
-        return result
+        result = []
+        for direction_idx, direction in enumerate(stencil):
+            if all(d_i == 0 for d_i in direction):
+                assert direction_idx == 0
+                continue
+            for own_cell_idx, cell in enumerate(mapping.fluid_coordinates):
+                #inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i in zip(cell, direction)])
+                inv_neighbor_cell = [(write_i - dir_i)%ds_i for write_i, dir_i, ds_i in zip(cell, direction, mapping.domain_size)]
+                if mapping.flag_array[tuple(inv_neighbor_cell)] & boundary_mask: 
+                    write = pdf_index(own_cell_idx, direction_idx, len(mapping))
+                    read = pdf_index(own_cell_idx, inverse_idx(stencil, direction_idx), len(mapping))
+                    result.append([direction_idx, write, read])
+
+        return np.array(result, dtype=np.int64)
 
     def assignments(self):
-        return [BoundaryOffsetInfo(self.method.stencil),
-                LbmWeightInfo(self.method),
-                Assignment(self.DIR_SYMBOL, self.index_field(self.DIR_SYMBOL.name)),
+        return [Assignment(self.DIR_SYMBOL, self.index_field(0)),
                 *self.boundary_eqs]
+
+    def get_kernel(self):
+        kernel = create_kernel(self.assignments(), ghost_layers=0)
+
+        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)
+        return kernel

lbmpy_tests/Simple_LB_kernel.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from lbmpy.session import *
+
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+
+from IPython.core.interactiveshell import InteractiveShell
+InteractiveShell.ast_node_interactivity = "all"
+```
+
+%% Cell type:code id: tags:
+
+``` python
+gpu = False
+target = 'cpu'
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reference_length = 15
+maximal_velocity = 0.005
+reynolds_number = 1000
+kinematic_vicosity = (reference_length * maximal_velocity) / reynolds_number
+initial_velocity=(maximal_velocity, 0)
+
+omega = relaxation_rate_from_lattice_viscosity(kinematic_vicosity)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = get_stencil('D2Q9')
+#domain_size = (reference_length * 12, reference_length * 4)
+domain_size = (400, 100)
+dim = len(domain_size)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dh = ps.create_data_handling(domain_size=domain_size, periodicity=(False, False))
+
+src = dh.add_array('src', values_per_cell=len(stencil), alignment=True)
+dh.fill('src', 0.0, ghost_layers=True)
+dst = dh.add_array('dst', values_per_cell=len(stencil), alignment=True)
+dh.fill('dst', 0.0, ghost_layers=True)
+
+velField = dh.add_array('velField', values_per_cell=dh.dim, alignment=True)
+dh.fill('velField', 0.0, ghost_layers=True)
+
+domain_size
+type(dh.cpu_arrays[src.name])
+dh.cpu_arrays[src.name].nbytes
+dh.cpu_arrays[src.name].size
+```
+
+%% Output
+
+
+    $\displaystyle \left( 400, \  100\right)$
+    (400, 100)
+
+    numpy.ndarray
+
+
+    $\displaystyle 2952288$
+    2952288
+
+
+    $\displaystyle 369036$
+    369036
+
+%% Cell type:code id: tags:
+
+``` python
+method = create_lb_method(stencil='D2Q9', method='cumulant', relaxation_rate=omega)
+method
+```
+
+%% Output
+
+    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f2774701460>
+
+%% Cell type:markdown id: tags:
+
+# Initialisation
+
+%% Cell type:code id: tags:
+
+``` python
+opt =  None # {'instruction_set': 'avx', 'assume_aligned': False, 'nontemporal': False, 'assume_inner_stride_one': False}
+```
+
+%% Cell type:code id: tags:
+
+``` python
+init = pdf_initialization_assignments(method, 1.0, initial_velocity, src.center_vector)
+
+ast_init = ps.create_kernel(init, target=dh.default_target)
+kernel_init = ast_init.compile()
+
+dh.run_kernel(kernel_init)
+type(src.center_vector)
+```
+
+%% Output
+
+    sympy.matrices.dense.MutableDenseMatrix
+
+%% Cell type:markdown id: tags:
+
+# Update Rules
+
+%% Cell type:code id: tags:
+
+``` python
+update = create_lb_update_rule(lb_method=method,
+                               output={'velocity': velField},
+                               optimization={"symbolic_field": src,
+                                             "symbolic_temporary_field": dst},
+                               kernel_type='stream_pull_collide')
+
+ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=True)
+kernel = ast_kernel.compile()
+```
+
+%% Cell type:markdown id: tags:
+
+# Boundary Handling
+
+%% Cell type:code id: tags:
+
+``` python
+def set_obstacle(x, y, *_):
+    size = 5
+    left = domain_size[0] // size
+    right = domain_size[0] - domain_size[0] // size
+    lower = domain_size[1] // size
+    upper = domain_size[1] - domain_size[1] // size
+    left = 25
+    right = 373
+    lower = 7
+    upper = 93
+    print("obstacle", left, right, lower, upper)
+    return np.logical_not(np.logical_or(np.logical_or(x < left, x > right), np.logical_or(y < lower, y > upper)))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name="bh")
+
+inflow = UBB(initial_velocity)
+outflow = ExtrapolationOutflow(stencil[4], method)
+wall = NoSlip("wall")
+
+bh.set_boundary(inflow, slice_from_direction('W', dim))
+bh.set_boundary(outflow, slice_from_direction('E', dim))
+for direction in ('N', 'S'):
+    bh.set_boundary(wall, slice_from_direction(direction, dim))
+
+bh.set_boundary(NoSlip("obstacle"), mask_callback=set_obstacle)
+
+plt.figure(dpi=200)
+plt.boundary_handling(bh)
+```
+
+%% Output
+
+    2
+
+    4
+
+    8
+
+    8
+
+    obstacle 25 373 7 93
+
+    16
+
+    <Figure size 3200x1200 with 0 Axes>
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+obstacle_flag = bh.get_flag(NoSlip("obstacle"))
+fluid_flag = 1
+
+number_of_fluid_cells = len(np.where(dh.cpu_arrays['bhFlags'] == 1)[0])
+number_of_obstacle_cells = len(np.where(dh.cpu_arrays['bhFlags'] == 16)[0])
+
+
+ratio = number_of_fluid_cells / number_of_obstacle_cells
+print(ratio)
+```
+
+%% Output
+
+    0.33654103180967654
+
+%% Cell type:code id: tags:
+
+``` python
+def timeloop(timeSteps):
+    for i in range(timeSteps):
+        bh()
+        dh.run_kernel(kernel)
+        dh.swap("src", "dst")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+mask = np.fromfunction(set_obstacle, (domain_size[0], domain_size[1], len(domain_size)))
+
+def run():
+    timeloop(100)
+    #print(np.ma.array(dh.gather_array('velField'), mask=mask))
+    return np.ma.array(dh.gather_array('velField'), mask=mask)
+
+#animation = plt.vector_field_magnitude_animation(run, frames=1000, rescale=True)
+#set_display_mode('video')
+#res = display_animation(animation)
+
+#res
+```
+
+%% Output
+
+    obstacle 25 373 7 93
+
+%% Cell type:code id: tags:
+
+``` python
+```

lbmpy_tests/Simple_LB_kernel_comm.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from lbmpy.session import *
+
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+
+from lbmpy.advanced_streaming.utility import get_timesteps
+```
+
+%% Cell type:code id: tags:
+
+``` python
+gpu = False
+target = 'cpu'
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reference_length = 30
+maximal_velocity = 0.05
+reynolds_number = 100000
+kinematic_vicosity = (reference_length * maximal_velocity) / reynolds_number
+initial_velocity=(maximal_velocity, 0)
+
+omega = relaxation_rate_from_lattice_viscosity(kinematic_vicosity)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+omega
+```
+
+%% Output
+
+
+    $\displaystyle 1.9998200161985422$
+    1.9998200161985422
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = get_stencil('D2Q9')
+#domain_size = (reference_length * 12, reference_length * 4)
+domain_size = (3,3)
+dim = len(domain_size)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ps.stencil.plot_2d(stencil)
+```
+
+%% Output
+
+
+
+%% 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), alignment=True)
+dh.fill('src', 0.0, ghost_layers=True)
+dst = dh.add_array('dst', values_per_cell=len(stencil), alignment=True)
+dh.fill('dst', 0.0, ghost_layers=True)
+
+velField = dh.add_array('velField', values_per_cell=dh.dim, alignment=True)
+dh.fill('velField', maximal_velocity, 0, ghost_layers=True)
+dh.fill('velField', 0.0, 1, ghost_layers=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+method = create_lb_method(stencil='D2Q9', method='cumulant', relaxation_rate=omega)
+method
+```
+
+%% Output
+
+    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f13c73de370>
+
+%% Cell type:markdown id: tags:
+
+# Initialisation
+
+%% Cell type:code id: tags:
+
+``` python
+opt =  None # {'instruction_set': 'avx', 'assume_aligned': False, 'nontemporal': False, 'assume_inner_stride_one': False}
+```
+
+%% Cell type:code id: tags:
+
+``` python
+init = pdf_initialization_assignments(method, 1.0, initial_velocity, src.center_vector)
+
+ast_init = ps.create_kernel(init, target=dh.default_target)
+kernel_init = ast_init.compile()
+
+dh.run_kernel(kernel_init)
+```
+
+%% Cell type:markdown id: tags:
+
+# Update Rules
+
+%% Cell type:code id: tags:
+
+``` python
+update = create_lb_update_rule(lb_method=method,
+                               output={'velocity': velField},
+                               optimization={"symbolic_field": src,
+                                             "symbolic_temporary_field": dst},
+                               kernel_type='stream_pull_collide')
+
+ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=True)
+kernel = ast_kernel.compile()
+```
+
+%% Cell type:markdown id: tags:
+
+# Boundary Handling
+
+%% Cell type:code id: tags:
+
+``` python
+def set_obstacle(x, y, *_):
+    size = 2
+    left = domain_size[0] // size
+    right = domain_size[0] - domain_size[0] // size
+    lower = domain_size[1] // size
+    upper = domain_size[1] - domain_size[1] // size
+    return np.logical_not(np.logical_or(np.logical_or(x < left, x > right), np.logical_or(y < lower, y > upper)))
+
+
+def set_sphere(x, y, *_):
+    mid = (domain_size[0] // 3, domain_size[1] // 2)
+    radius = reference_length // 2
+    return (x-mid[0])**2 + (y-mid[1])**2 < radius**2
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name="bh")
+
+inflow = UBB(initial_velocity)
+outflow = ExtrapolationOutflow(stencil[4], method)
+wall = NoSlip("wall")
+#freeslip = FreeSlip(stencil)
+
+# bh.set_boundary(inflow, slice_from_direction('W', dim))
+# bh.set_boundary(outflow, slice_from_direction('E', dim))
+#bh.set_boundary(wall, slice_from_direction('S', dim))
+#bh.set_boundary(wall, slice_from_direction('N', dim))
+
+#bh.set_boundary(NoSlip("obstacle"), mask_callback=set_sphere)
+
+plt.figure(dpi=200)
+plt.boundary_handling(bh)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['bhFlags'][0, :]
+```
+
+%% Output
+
+    array([1, 1, 1, 1, 1], dtype=uint32)
+
+%% Cell type:code id: tags:
+
+``` python
+dh.periodicity
+```
+
+%% Output
+
+    (True, True)
+
+%% Cell type:code id: tags:
+
+``` python
+periodic = LBMPeriodicityHandling(stencil=stencil, data_handling=dh,
+                                  pdf_field_name=src.name, streaming_pattern='pull')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# obstacle_flag = bh.get_flag(NoSlip("obstacle"))
+# fluid_flag = 1
+
+# number_of_fluid_cells = len(np.where(dh.cpu_arrays['bhFlags'] == 1)[0])
+# number_of_obstacle_cells = len(np.where(dh.cpu_arrays['bhFlags'] == 16)[0])
+
+
+# ratio = number_of_fluid_cells / number_of_obstacle_cells
+# print(ratio)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def timeloop(timeSteps):
+    for i in range(timeSteps):
+        periodic()
+        bh()
+        dh.run_kernel(kernel)
+        dh.swap("src", "dst")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+mask = np.fromfunction(set_sphere, (domain_size[0], domain_size[1], len(domain_size)))
+
+def run():
+    timeloop(100)
+    return np.ma.array(dh.gather_array('velField'), mask=mask)
+
+#animation = plt.vector_field_magnitude_animation(run, frames=5, rescale=True)
+#set_display_mode('video')
+#res = display_animation(animation)
+
+#res
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src'][0, :, 3]
+```
+
+%% Output
+
+    array([0.09527778, 0.09527778, 0.09527778, 0.09527778, 0.09527778])
+
+%% Cell type:code id: tags:
+
+``` python
+pdf_field = dh.fields[src.name]
+kernels = []
+for src_slice, dst_slice in periodic.comm_slices[0]:
+
+    pdf_idx = src_slice[-1]
+    assert isinstance(pdf_idx, int), "PDF index needs to be an integer constant"
+    assert pdf_idx == dst_slice[-1], "Source and Destination PDF indices must be equal"
+    print("idx", pdf_idx)
+    print("src", src_slice)
+    print("dst", dst_slice)
+    src_slice = src_slice[:-1]
+    dst_slice = dst_slice[:-1]
+
+    if domain_size is None:
+        domain_size = pdf_field.spatial_shape
+
+    normalized_from_slice = normalize_slice(src_slice, domain_size)
+    normalized_to_slice = normalize_slice(dst_slice, domain_size)
+
+    def _start(s):
+        print("start", s)
+        return s.start if isinstance(s, slice) else s
+
+    def _stop(s):
+        print("stop", s)
+        return s.stop if isinstance(s, slice) else s
+
+    offset = [_start(s1) - _start(s2) for s1, s2 in zip(normalized_from_slice, normalized_to_slice)]
+    print("offset", offset)
+    assert offset == [_stop(s1) - _stop(s2) for s1, s2 in zip(normalized_from_slice, normalized_to_slice)], \
+        "Slices have to have same size"
+
+    copy_eq = Assignment(pdf_field(pdf_idx), pdf_field[tuple(offset)](pdf_idx))
+    print(copy_eq)
+    print("src", src_slice)
+    print("dst", dst_slice)
+
+    print("------")
+```
+
+%% Output
+
+    idx 7
+    src (1, 1, 7)
+    dst (-1, -1, 7)
+    start 1
+    start 2
+    start 1
+    start 2
+    offset [-1, -1]
+    stop 1
+    stop 2
+    stop 1
+    stop 2
+    src[0,0](7) ← src[-1,-1](7)
+    src (1, 1)
+    dst (-1, -1)
+    ------
+    idx 3
+    src (1, slice(1, -1, None), 3)
+    dst (-1, slice(1, -1, None), 3)
+    start 1
+    start 2
+    start slice(1, 2, 1)
+    start slice(1, 2, 1)
+    offset [-1, 0]
+    stop 1
+    stop 2
+    stop slice(1, 2, 1)
+    stop slice(1, 2, 1)
+    src[0,0](3) ← src[-1,0](3)
+    src (1, slice(1, -1, None))
+    dst (-1, slice(1, -1, None))
+    ------
+    idx 5
+    src (1, slice(1, -2, None), 5)
+    dst (-1, slice(1, -2, None), 5)
+    start 1
+    start 2
+    start slice(1, 1, 1)
+    start slice(1, 1, 1)
+    offset [-1, 0]
+    stop 1
+    stop 2
+    stop slice(1, 1, 1)
+    stop slice(1, 1, 1)
+    src[0,0](5) ← src[-1,0](5)
+    src (1, slice(1, -2, None))
+    dst (-1, slice(1, -2, None))
+    ------
+    idx 7
+    src (1, slice(2, -1, None), 7)
+    dst (-1, slice(2, -1, None), 7)
+    start 1
+    start 2
+    start slice(2, 2, 1)
+    start slice(2, 2, 1)
+    offset [-1, 0]
+    stop 1
+    stop 2
+    stop slice(2, 2, 1)
+    stop slice(2, 2, 1)
+    src[0,0](7) ← src[-1,0](7)
+    src (1, slice(2, -1, None))
+    dst (-1, slice(2, -1, None))
+    ------
+    idx 5
+    src (1, -2, 5)
+    dst (-1, 0, 5)
+    start 1
+    start 2
+    start 1
+    start 0
+    offset [-1, 1]
+    stop 1
+    stop 2
+    stop 1
+    stop 0
+    src[0,0](5) ← src[-1,1](5)
+    src (1, -2)
+    dst (-1, 0)
+    ------
+    idx 7
+    src (slice(2, -1, None), 1, 7)
+    dst (slice(2, -1, None), -1, 7)
+    start slice(2, 2, 1)
+    start slice(2, 2, 1)
+    start 1
+    start 2
+    offset [0, -1]
+    stop slice(2, 2, 1)
+    stop slice(2, 2, 1)
+    stop 1
+    stop 2
+    src[0,0](7) ← src[0,-1](7)
+    src (slice(2, -1, None), 1)
+    dst (slice(2, -1, None), -1)
+    ------
+    idx 2
+    src (slice(1, -1, None), 1, 2)
+    dst (slice(1, -1, None), -1, 2)
+    start slice(1, 2, 1)
+    start slice(1, 2, 1)
+    start 1
+    start 2
+    offset [0, -1]
+    stop slice(1, 2, 1)
+    stop slice(1, 2, 1)
+    stop 1
+    stop 2
+    src[0,0](2) ← src[0,-1](2)
+    src (slice(1, -1, None), 1)
+    dst (slice(1, -1, None), -1)
+    ------
+    idx 8
+    src (slice(1, -2, None), 1, 8)
+    dst (slice(1, -2, None), -1, 8)
+    start slice(1, 1, 1)
+    start slice(1, 1, 1)
+    start 1
+    start 2
+    offset [0, -1]
+    stop slice(1, 1, 1)
+    stop slice(1, 1, 1)
+    stop 1
+    stop 2
+    src[0,0](8) ← src[0,-1](8)
+    src (slice(1, -2, None), 1)
+    dst (slice(1, -2, None), -1)
+    ------
+    idx 1
+    src (slice(1, -1, None), -2, 1)
+    dst (slice(1, -1, None), 0, 1)
+    start slice(1, 2, 1)
+    start slice(1, 2, 1)
+    start 1
+    start 0
+    offset [0, 1]
+    stop slice(1, 2, 1)
+    stop slice(1, 2, 1)
+    stop 1
+    stop 0
+    src[0,0](1) ← src[0,1](1)
+    src (slice(1, -1, None), -2)
+    dst (slice(1, -1, None), 0)
+    ------
+    idx 5
+    src (slice(2, -1, None), -2, 5)
+    dst (slice(2, -1, None), 0, 5)
+    start slice(2, 2, 1)
+    start slice(2, 2, 1)
+    start 1
+    start 0
+    offset [0, 1]
+    stop slice(2, 2, 1)
+    stop slice(2, 2, 1)
+    stop 1
+    stop 0
+    src[0,0](5) ← src[0,1](5)
+    src (slice(2, -1, None), -2)
+    dst (slice(2, -1, None), 0)
+    ------
+    idx 6
+    src (slice(1, -2, None), -2, 6)
+    dst (slice(1, -2, None), 0, 6)
+    start slice(1, 1, 1)
+    start slice(1, 1, 1)
+    start 1
+    start 0
+    offset [0, 1]
+    stop slice(1, 1, 1)
+    stop slice(1, 1, 1)
+    stop 1
+    stop 0
+    src[0,0](6) ← src[0,1](6)
+    src (slice(1, -2, None), -2)
+    dst (slice(1, -2, None), 0)
+    ------
+    idx 8
+    src (-2, 1, 8)
+    dst (0, -1, 8)
+    start 1
+    start 0
+    start 1
+    start 2
+    offset [1, -1]
+    stop 1
+    stop 0
+    stop 1
+    stop 2
+    src[0,0](8) ← src[1,-1](8)
+    src (-2, 1)
+    dst (0, -1)
+    ------
+    idx 4
+    src (-2, slice(1, -1, None), 4)
+    dst (0, slice(1, -1, None), 4)
+    start 1
+    start 0
+    start slice(1, 2, 1)
+    start slice(1, 2, 1)
+    offset [1, 0]
+    stop 1
+    stop 0
+    stop slice(1, 2, 1)
+    stop slice(1, 2, 1)
+    src[0,0](4) ← src[1,0](4)
+    src (-2, slice(1, -1, None))
+    dst (0, slice(1, -1, None))
+    ------
+    idx 6
+    src (-2, slice(1, -2, None), 6)
+    dst (0, slice(1, -2, None), 6)
+    start 1
+    start 0
+    start slice(1, 1, 1)
+    start slice(1, 1, 1)
+    offset [1, 0]
+    stop 1
+    stop 0
+    stop slice(1, 1, 1)
+    stop slice(1, 1, 1)
+    src[0,0](6) ← src[1,0](6)
+    src (-2, slice(1, -2, None))
+    dst (0, slice(1, -2, None))
+    ------
+    idx 8
+    src (-2, slice(2, -1, None), 8)
+    dst (0, slice(2, -1, None), 8)
+    start 1
+    start 0
+    start slice(2, 2, 1)
+    start slice(2, 2, 1)
+    offset [1, 0]
+    stop 1
+    stop 0
+    stop slice(2, 2, 1)
+    stop slice(2, 2, 1)
+    src[0,0](8) ← src[1,0](8)
+    src (-2, slice(2, -1, None))
+    dst (0, slice(2, -1, None))
+    ------
+    idx 6
+    src (-2, -2, 6)
+    dst (0, 0, 6)
+    start 1
+    start 0
+    start 1
+    start 0
+    offset [1, 1]
+    stop 1
+    stop 0
+    stop 1
+    stop 0
+    src[0,0](6) ← src[1,1](6)
+    src (-2, -2)
+    dst (0, 0)
+    ------
+
+%% Cell type:code id: tags:
+
+``` python
+copy_eq
+ps.stencil.plot_2d(stencil)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+periodic.comm_slices
+```
+
+%% Output
+
+    [[((1, 1, 7), (-1, -1, 7)),
+      ((1, slice(1, -1, None), 3), (-1, slice(1, -1, None), 3)),
+      ((1, slice(1, -2, None), 5), (-1, slice(1, -2, None), 5)),
+      ((1, slice(2, -1, None), 7), (-1, slice(2, -1, None), 7)),
+      ((1, -2, 5), (-1, 0, 5)),
+      ((slice(2, -1, None), 1, 7), (slice(2, -1, None), -1, 7)),
+      ((slice(1, -1, None), 1, 2), (slice(1, -1, None), -1, 2)),
+      ((slice(1, -2, None), 1, 8), (slice(1, -2, None), -1, 8)),
+      ((slice(1, -1, None), -2, 1), (slice(1, -1, None), 0, 1)),
+      ((slice(2, -1, None), -2, 5), (slice(2, -1, None), 0, 5)),
+      ((slice(1, -2, None), -2, 6), (slice(1, -2, None), 0, 6)),
+      ((-2, 1, 8), (0, -1, 8)),
+      ((-2, slice(1, -1, None), 4), (0, slice(1, -1, None), 4)),
+      ((-2, slice(1, -2, None), 6), (0, slice(1, -2, None), 6)),
+      ((-2, slice(2, -1, None), 8), (0, slice(2, -1, None), 8)),
+      ((-2, -2, 6), (0, 0, 6))]]
+
+%% Cell type:code id: tags:
+
+``` python
+get_timesteps('pull')[0]
+```
+
+%% Output
+
+
+    $\displaystyle Both$
+    <Timestep.BOTH: 2>
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src'][0, :, 3]
+```
+
+%% Output
+
+    array([0.09527778, 0.09527778, 0.09527778, 0.09527778, 0.09527778])
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src'][0, slice(0, -1, None), 3]
+```
+
+%% Output
+
+    array([0.09527778, 0.09527778, 0.09527778, 0.09527778])
+
+%% Cell type:code id: tags:
+
+``` python
+periodic_kernel = periodic._compile_copy_kernel()
+```
+
+%% Output
+
+    ---------------------------------------------------------------------------
+    AttributeError                            Traceback (most recent call last)
+    <ipython-input-27-db85c5babc91> in <module>
+    ----> 1 periodic_kernel = periodic._compile_copy_kernel()
+
+    AttributeError: 'LBMPeriodicityHandling' object has no attribute '_compile_copy_kernel'

lbmpy_tests/Simple_LB_kernel_comm_2.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from lbmpy.session import *
+
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+
+from lbmpy.advanced_streaming.utility import get_timesteps
+```
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = get_stencil('D2Q9')
+domain_size = (3, 3)
+dim = len(domain_size)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ps.stencil.plot_2d(stencil)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, True), default_target="cpu")
+
+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)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+shape = dh.cpu_arrays['src'].shape
+
+num = 0
+for i in range(shape[0]):
+    for j in range(shape[1]):
+        for k in range(shape[2]):
+            dh.cpu_arrays['src'][i, j, k] = num
+            num += 1
+```
+
+%% Cell type:code id: tags:
+
+``` python
+periodic = LBMPeriodicityHandling(stencil=stencil, data_handling=dh,
+                                  pdf_field_name=src.name, streaming_pattern='pull')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# periodic()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# dh.cpu_arrays['src'][4, 4]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# dh.cpu_arrays['src'][1, 1]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# obstacle_flag = bh.get_flag(NoSlip("obstacle"))
+# fluid_flag = 1
+
+# number_of_fluid_cells = len(np.where(dh.cpu_arrays['bhFlags'] == 1)[0])
+# number_of_obstacle_cells = len(np.where(dh.cpu_arrays['bhFlags'] == 16)[0])
+
+
+# ratio = number_of_fluid_cells / number_of_obstacle_cells
+# print(ratio)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# from pystencils.gpucuda.kernelcreation import create_cuda_kernel
+
+pdf_field = dh.fields[src.name]
+kernels = []
+
+num = 0
+for src_slice, dst_slice in periodic.comm_slices[0]:
+
+    pdf_idx = src_slice[-1]
+    assert isinstance(pdf_idx, int), "PDF index needs to be an integer constant"
+    assert pdf_idx == dst_slice[-1], "Source and Destination PDF indices must be equal"
+    src_slice = src_slice[:-1]
+    dst_slice = dst_slice[:-1]
+
+    normalized_from_slice = normalize_slice(src_slice, pdf_field.spatial_shape)
+    normalized_to_slice = normalize_slice(dst_slice, pdf_field.spatial_shape)
+
+    def _start(s):
+        return s.start if isinstance(s, slice) else s
+
+    def _stop(s):
+        return s.stop if isinstance(s, slice) else s
+
+    offset = [_start(s1) - _start(s2) for s1, s2 in zip(normalized_from_slice, normalized_to_slice)]
+    assert offset == [_stop(s1) - _stop(s2) for s1, s2 in zip(normalized_from_slice, normalized_to_slice)], \
+        "Slices have to have same size"
+
+    copy_eq = Assignment(pdf_field(pdf_idx), pdf_field[tuple(offset)](pdf_idx))
+    ast = create_kernel([copy_eq], iteration_slice=dst_slice, skip_independence_check=True)
+    kernel = ast.compile()
+
+    dh.run_kernel(kernel)
+    print(copy_eq)
+    print(dst_slice)
+
+    if num == 1:
+        break
+    num += 1
+```
+
+%% Output
+
+    src[0,0](7) ← src[-3,-3](7)
+    (-1, -1)
+    src[0,0](3) ← src[-3,0](3)
+    (-1, slice(1, -1, None))
+
+%% Cell type:code id: tags:
+
+``` python
+ps.show_code(ast)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+kernel = ast.compile()
+
+dh.run_kernel(kernel)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src'][2, 4]
+```
+
+%% Output
+
+    array([126., 127., 101., 129., 130., 131., 132., 106., 107.])
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src'][3, 3]
+```
+
+%% Output
+
+    array([162., 163., 164., 165., 166., 167., 168., 169., 170.])
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src']
+```
+
+%% 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.],
+            [108., 109., 110., 111., 112., 113., 114., 115., 116.],
+            [117., 118., 119., 120., 121., 122., 123., 124., 125.],
+            [126., 127., 128., 129., 130., 131., 132., 133., 134.]],
+    
+           [[135., 136., 137., 138., 139., 140., 141., 142., 143.],
+            [144., 145., 146., 147., 148., 149., 150., 151., 152.],
+            [153., 154., 155., 156., 157., 158., 159., 160., 161.],
+            [162., 163., 164., 165., 166., 167., 168., 169., 170.],
+            [171., 172., 173., 174., 175., 176., 177., 178., 179.]],
+    
+           [[180., 181., 182., 183., 184., 185., 186., 187., 188.],
+            [189., 190., 191., 192., 193., 194., 195., 196., 197.],
+            [198., 199., 200., 201., 202., 203., 204., 205., 206.],
+            [207., 208., 209., 210., 211., 212., 213., 214., 215.],
+            [216., 217., 218., 219., 220., 221., 222.,  61., 224.]]])
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays['src'][0, slice(0, -1, None), 3]
+```
+
+%% Output
+
+    array([ 3., 12., 21., 30.])

lbmpy_tests/test_esotwist_visualization.ipynb

@@ -127,7 +127,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.5"
+   "version": "3.8.2"
   }
  },
  "nbformat": 4,

lbmpy_tests/test_single_boundary_node.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from lbmpy.session import *
+
+from lbmpy.stencils import get_stencil
+from lbmpy.updatekernels import create_stream_only_kernel
+from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
+```
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = get_stencil("D2Q9")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+domain_size = (3, 2)
+dim = len(domain_size)
+omega = 1.98
+```
+
+%% Cell type:code id: tags:
+
+``` python
+lb_method = create_lb_method(stencil=stencil, method='mrt', relaxation_rates=[omega, 1, 1, 1])
+```
+
+%% 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)
+denField = dh.add_array('denField', values_per_cell=1)
+dh.fill('denField', 1.0, ghost_layers=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# init = pdf_initialization_assignments(lb_method, density=denField.center,
+#                                      velocity=velField.center_vector, pdfs=src.center_vector)
+
+# ast_init = ps.create_kernel(init, target=dh.default_target)
+# kernel_init = ast_init.compile()
+
+# dh.run_kernel(kernel_init)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+update = create_stream_only_kernel(stencil, src, dst_field=dst, accessor=StreamPullTwoFieldsAccessor())
+
+ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=True)
+kernel = ast_kernel.compile()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+#update
+```
+
+%% Cell type:code id: tags:
+
+``` python
+periodic = dh.synchronization_function(src.name, target=dh.default_target)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, src.name, name="bh")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def set_sphere(x, y, *_):
+    result = np.zeros(x.shape, dtype=bool)
+    result[3, 1] = True
+    result[0,3] = True
+    #for i in range(1, x.shape[1] - 1):
+    #    if i == 2:
+    #        result[1, i] = True
+    return result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bh.set_boundary(NoSlip("obstacle"), mask_callback=set_sphere)
+
+plt.figure(dpi=200)
+plt.boundary_handling(bh)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+i = 0
+for ii in range(dh.cpu_arrays[src.name].shape[0]):
+    j = 0
+    for jj in range(dh.cpu_arrays[src.name].shape[1]):
+        k = 0
+        for kk in range(dh.cpu_arrays[src.name].shape[2]):
+            dh.cpu_arrays[src.name][ii,jj,kk] = i+j+k
+            k = k+1
+        j = j+10
+    i = i+100
+dh.cpu_arrays[src.name]
+```
+
+%% Output
+
+    array([[[  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.],
+            [ 10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.],
+            [ 20.,  21.,  22.,  23.,  24.,  25.,  26.,  27.,  28.],
+            [ 30.,  31.,  32.,  33.,  34.,  35.,  36.,  37.,  38.]],
+    
+           [[100., 101., 102., 103., 104., 105., 106., 107., 108.],
+            [110., 111., 112., 113., 114., 115., 116., 117., 118.],
+            [120., 121., 122., 123., 124., 125., 126., 127., 128.],
+            [130., 131., 132., 133., 134., 135., 136., 137., 138.]],
+    
+           [[200., 201., 202., 203., 204., 205., 206., 207., 208.],
+            [210., 211., 212., 213., 214., 215., 216., 217., 218.],
+            [220., 221., 222., 223., 224., 225., 226., 227., 228.],
+            [230., 231., 232., 233., 234., 235., 236., 237., 238.]],
+    
+           [[300., 301., 302., 303., 304., 305., 306., 307., 308.],
+            [310., 311., 312., 313., 314., 315., 316., 317., 318.],
+            [320., 321., 322., 323., 324., 325., 326., 327., 328.],
+            [330., 331., 332., 333., 334., 335., 336., 337., 338.]],
+    
+           [[400., 401., 402., 403., 404., 405., 406., 407., 408.],
+            [410., 411., 412., 413., 414., 415., 416., 417., 418.],
+            [420., 421., 422., 423., 424., 425., 426., 427., 428.],
+            [430., 431., 432., 433., 434., 435., 436., 437., 438.]]])
+
+%% Cell type:code id: tags:
+
+``` python
+periodic()
+bh()
+dh.run_kernel(kernel)
+# dh.swap("src", "dst")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays[src.name]
+```
+
+%% Output
+
+    array([[[320., 321., 322., 323., 324., 325., 326., 327., 328.],
+            [310., 311., 312., 313., 314., 315., 316., 317., 318.],
+            [320., 321., 322., 323., 324., 325., 326., 327., 328.],
+            [310., 311., 312., 313., 314., 315., 316., 317., 125.]],
+    
+           [[120., 121., 122., 123., 124., 125., 126., 127., 128.],
+            [110., 111., 112., 113., 114., 115., 116., 117., 118.],
+            [120., 121., 122., 123., 124., 125., 126., 127., 128.],
+            [110., 111., 112., 113., 114., 115., 116., 117., 118.]],
+    
+           [[220., 221., 222., 223., 224., 225., 226., 227., 228.],
+            [210., 211., 212., 213., 214., 215., 216., 217., 218.],
+            [220., 221., 222., 223., 224., 225., 226., 227., 228.],
+            [210., 211., 212., 213., 214., 215., 216., 217., 218.]],
+    
+           [[320., 321., 322., 323., 324., 325., 326., 327., 328.],
+            [310., 322., 312., 214., 314., 228., 316., 317., 318.],
+            [320., 321., 322., 323., 324., 325., 326., 327., 328.],
+            [310., 311., 312., 313., 314., 315., 316., 317., 318.]],
+    
+           [[120., 121., 122., 123., 124., 125., 126., 127., 128.],
+            [110., 111., 112., 113., 114., 115., 116., 117., 118.],
+            [120., 121., 122., 123., 124., 125., 126., 127., 128.],
+            [110., 111., 112., 113., 114., 115., 116., 117., 118.]]])
+
+%% Cell type:code id: tags:
+
+``` python
+dh.cpu_arrays[dst.name]
+```
+
+%% Output
+
+    array([[[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.]],
+    
+           [[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [110., 121., 122., 213., 314., 225., 326., 227., 328.],
+            [120., 111., 112., 223., 324., 215., 316., 217., 125.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.]],
+    
+           [[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [210., 221., 222., 214., 114., 325., 126., 327., 128.],
+            [220., 211., 212., 323., 124., 228., 116., 317., 118.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.]],
+    
+           [[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [310., 321., 322., 113., 214., 125., 226., 127., 228.],
+            [320., 322., 312., 123., 224., 115., 216., 117., 218.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.]],
+    
+           [[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+            [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.]]])
+
+%% Cell type:code id: tags:
+
+``` python
+```

lbmpy_tests/test_sparse_lbm.ipynb

-%% Cell type:code id: tags:
-
-``` python
-import pytest
-pytest.importorskip('pycuda')
-```
-
-%% Output
-
-    <module 'pycuda' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/pycuda/__init__.py'>
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.sparse import *
-from pystencils.field import FieldType
-import pycuda.gpuarray as gpuarray
-
-import pystencils as ps
-```
-
-%% Cell type:markdown id: tags:
-
-# Sparse (list based) LBM
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (20, 20)
-omega = 1.8
-target = 'cpu'
-
-ghost_layers = 1
-arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
-lid_velocity = 0.01
-force = 1e-6
-
-channel = True
-if channel:
-    kwargs={'force': (force, 0)}
-else:
-    kwargs = {}
-
-method = create_lb_method(stencil='D2Q9', relaxation_rate=omega, compressible=False, **kwargs)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-ubb = UBB( (lid_velocity, 0) )
-noslip = NoSlip()
-flags = {
-    'fluid': 1,
-    noslip: 2,
-    ubb: 4,
-}
-flag_arr = np.zeros(arr_size, dtype=np.uint16)
-flag_arr.fill(flags['fluid'])
-
-if channel:
-    flag_arr[0, :] = 0
-    flag_arr[-1, :] = 0
-    flag_arr[:, 0] = flags[noslip]
-    flag_arr[:, -1] = flags[noslip]
-else:
-    flag_arr[:, -1] = flags[ubb]
-    flag_arr[:, 0] = flags[noslip]
-    flag_arr[0, :] = flags[noslip]
-    flag_arr[-1, :] = flags[noslip]
-
-plt.scalar_field(flag_arr)
-plt.colorbar();
-```
-
-%% Output
-
-
-
-%% Cell type:markdown id: tags:
-
-### Mappings
-
-%% Cell type:code id: tags:
-
-``` python
-mapping = SparseLbMapper(method.stencil, flag_arr, flags['fluid'], flags[noslip], flags[ubb])
-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(method.stencil)), order='f')
-pdf_arr_tmp = np.empty_like(pdf_arr)
-
-vel_arr = np.ones([mapping.num_fluid_cells, method.dim], order='f')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]",
-                                                #f=pdf_arr[:mapping.num_fluid_cells],
-                                                #d=pdf_arr_tmp[:mapping.num_fluid_cells],
-                                                #u=vel_arr
-                                               )
-pdf_field.field_type = FieldType.CUSTOM
-pdf_field.pdf_field_tmp = FieldType.CUSTOM
-```
-
-%% Cell type:markdown id: tags:
-
-### Macroscopic quantities
-
-%% Cell type:code id: tags:
-
-``` python
-cqc = method.conserved_quantity_computation
-inp_eqs = cqc.equilibrium_input_equations_from_init_values()
-setter_eqs = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
-setter_eqs = setter_eqs.new_with_substitutions({sym: pdf_field(i)
-                                                for i, sym in enumerate(method.post_collision_pdf_symbols)})
-kernel_initialize = ps.create_kernel(setter_eqs, ghost_layers=((0, 0),), ).compile()
-
-def init():
-    kernel_initialize(f=pdf_arr[:mapping.num_fluid_cells])
-init()
-
-
-getter_eqs = cqc.output_equations_from_pdfs(pdf_field.center_vector,
-                                            {'velocity': vel_field})
-kernel_compute_u = ps.create_kernel(getter_eqs, ghost_layers=((0,0),) ).compile()
-
-def get_velocity(arr=pdf_arr):
-    fluid_cell_arr = mapping.coordinates
-    kernel_compute_u(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
-    full_velocity_arr = np.zeros(flag_arr.shape + (2,), dtype=np.float64)
-    full_velocity_arr.fill(np.nan)
-    arr = fluid_cell_arr[:mapping.num_fluid_cells]
-    full_velocity_arr[arr['x'], arr['y']] = vel_arr
-    return full_velocity_arr[1:-1, 1:-1]
-```
-
-%% Cell type:markdown id: tags:
-
-### Stream Collide Kernel
-
-%% Cell type:code id: tags:
-
-``` python
-#index_field = ps.Field.create_from_numpy_array("idx", index_arr, index_dimensions=1)
-index_field = ps.Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=index_arr.dtype)
-
-collision_rule = method.get_collision_rule()
-
-Q = len(method.stencil)
-symbol_subs = {sym: pdf_field.absolute_access((index_field(i-1),),())
-               for i, sym in enumerate(method.pre_collision_pdf_symbols)}
-
-symbol_subs.update({sym: pdf_field_tmp(i) for i, sym in enumerate(method.post_collision_pdf_symbols)})
-
-symbol_subs[method.pre_collision_pdf_symbols[0]] = pdf_field(0) # special case for center
-symbol_subs
-```
-
-%% Output
-
-
-    $\displaystyle \left\{ d_{0} : {{d}_{0}^{0}}, \  d_{1} : {{d}_{0}^{1}}, \  d_{2} : {{d}_{0}^{2}}, \  d_{3} : {{d}_{0}^{3}}, \  d_{4} : {{d}_{0}^{4}}, \  d_{5} : {{d}_{0}^{5}}, \  d_{6} : {{d}_{0}^{6}}, \  d_{7} : {{d}_{0}^{7}}, \  d_{8} : {{d}_{0}^{8}}, \  f_{0} : {{f}_{0}^{0}}, \  f_{1} : {{f}_{\mathbf{{idx}_{0}^{0}}}^{0}}, \  f_{2} : {{f}_{\mathbf{{idx}_{0}^{1}}}^{0}}, \  f_{3} : {{f}_{\mathbf{{idx}_{0}^{2}}}^{0}}, \  f_{4} : {{f}_{\mathbf{{idx}_{0}^{3}}}^{0}}, \  f_{5} : {{f}_{\mathbf{{idx}_{0}^{4}}}^{0}}, \  f_{6} : {{f}_{\mathbf{{idx}_{0}^{5}}}^{0}}, \  f_{7} : {{f}_{\mathbf{{idx}_{0}^{6}}}^{0}}, \  f_{8} : {{f}_{\mathbf{{idx}_{0}^{7}}}^{0}}\right\}$
-    {d₀: d_C__0, d₁: d_C__1, d₂: d_C__2, d₃: d_C__3, d₄: d_C__4, d₅: d_C__5, d₆: d
-    _C__6, d₇: d_C__7, d₈: d_C__8, f₀: f_C__0, f₁: f_b59661980ed8, f₂: f_c4b733d80
-    293, f₃: f_46a6e6d3a3b5, f₄: f_31387ae7897f, f₅: f_24d70d6a047a, f₆: f_c9de861
-    6062f, f₇: f_c16814a33bbc, f₈: f_52ab73425808}
-
-%% Cell type:code id: tags:
-
-``` python
-collision_rule = method.get_collision_rule()
-update_rule = collision_rule.new_with_substitutions(symbol_subs)
-kernel_stream_collide = ps.create_kernel(update_rule, ghost_layers=[(0,0)], target=target).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-### Boundary Kernels
-
-%% Cell type:code id: tags:
-
-``` python
-if not channel:
-    if target == 'gpu':
-        raise NotImplementedError("UBB on GPU not working yet")
-
-    ubb_mapper = SparseLbBoundaryMapper(ubb, method, pdf_field)
-    #TODO the following line is wrong: kernel contains accesses to index_field and pdf_field which have
-    #different size: a correct kernel comes out when by change the shape is taken from index field,
-    # when taken from pdf field, a wrong kernel is generated
-    ubb_kernel = ps.create_kernel(ubb_mapper.assignments(), ghost_layers=0).compile()
-    ubb_idx_arr = ubb_mapper.create_index_arr(mapping, flags[ubb])
-    ps.show_code(ubb_kernel.ast)
-    def handle_ubb():
-        ubb_kernel(indexField=ubb_idx_arr, f=pdf_arr[:mapping.num_fluid_cells])
-else:
-    def handle_ubb():
-        pass
-```
-
-%% Cell type:markdown id: tags:
-
-### Time Loop
-
-%% Cell type:code id: tags:
-
-``` python
-def time_step():
-    global pdf_arr, pdf_arr_tmp, index_arr
-    handle_ubb()
-    if target == 'gpu':
-        gpu_pdf_arr = gpuarray.to_gpu(pdf_arr)
-        gpu_pdf_arr_tmp = gpuarray.to_gpu(pdf_arr_tmp)
-        gpu_index_arr = gpuarray.to_gpu(index_arr)
-
-        kernel_stream_collide(f=gpu_pdf_arr[:mapping.num_fluid_cells],
-                          d=gpu_pdf_arr_tmp[:mapping.num_fluid_cells],
-                          idx=gpu_index_arr)
-
-        pdf_arr = gpu_pdf_arr.get()
-        pdf_arr_tmp = gpu_pdf_arr_tmp.get()
-        index_arr = gpu_index_arr.get()
-    else:
-        kernel_stream_collide(f=pdf_arr[:mapping.num_fluid_cells],
-                              d=pdf_arr_tmp[:mapping.num_fluid_cells],
-                              idx=index_arr)
-    pdf_arr_tmp, pdf_arr = pdf_arr, pdf_arr_tmp
-
-def run(steps=100):
-    for t in range(steps):
-        time_step()
-    return get_velocity()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-result = run(100)
-plt.vector_field(result, step=1);
-```
-
-%% Output
-
-
-
-%% Cell type:markdown id: tags:
-
-### Check against reference
-
-%% Cell type:code id: tags:
-
-``` python
-if channel:
-    reference = create_channel(domain_size, force=force, lb_method=method)
-else:
-    reference = create_lid_driven_cavity(domain_size, relaxation_rate=omega, lid_velocity=lid_velocity,
-                                         compressible=False)
-reference.run(100)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-np.testing.assert_almost_equal(reference.velocity[:, :], result)
-```

lbmpy_tests/test_sparse_lbm_parallel_only.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import pytest
+import sys
+# 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 = 1
+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'])
+#flag_arr[3,1] = 2
+#flag_arr[0,3] = 2
+#flag_arr[3,3] = 2
+#flag_arr[0,1] = 2
+#plt.scalar_field(flag_arr)
+#plt.colorbar();
+```
+
+%% Cell type:markdown id: tags:
+
+### Mappings
+
+%% Cell type:code id: tags:
+
+``` python
+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
+parallel_mapper = SparseLbCommunicationMapper(mapping, dh, pdf_field)
+parallel_mapper.create_packages()
+parallel_mapper.receive_here()
+parallel_mapper.create_index_array()
+```
+
+%% Output
+
+    
+     New direction: (0, 1) ,  1
+    (1, 3) 8
+    (2, 3) 13
+    (3, 3) 18
+    (4, 3) 23
+    
+     New direction: (0, -1) ,  2
+    (1, 1) 6
+    (2, 1) 11
+    (3, 1) 16
+    (4, 1) 21
+    
+     New direction: (-1, 0) ,  3
+    (1, 1) 6
+    (1, 2) 7
+    (1, 3) 8
+    
+     New direction: (1, 0) ,  4
+    (4, 1) 21
+    (4, 2) 22
+    (4, 3) 23
+    
+     New direction: (-1, 1) ,  5
+    (1, 1) 6
+    (1, 2) 7
+    (1, 3) 8
+    (2, 3) 13
+    (3, 3) 18
+    (4, 3) 23
+    
+     New direction: (1, 1) ,  6
+    (1, 3) 8
+    (2, 3) 13
+    (3, 3) 18
+    (4, 1) 21
+    (4, 2) 22
+    (4, 3) 23
+    
+     New direction: (-1, -1) ,  7
+    (1, 1) 6
+    (1, 2) 7
+    (1, 3) 8
+    (2, 1) 11
+    (3, 1) 16
+    (4, 1) 21
+    
+     New direction: (1, -1) ,  8
+    (1, 1) 6
+    (2, 1) 11
+    (3, 1) 16
+    (4, 1) 21
+    (4, 2) 22
+    (4, 3) 23
+
+    [[[38, 43, 48, 53], [], [], [], [163, 168, 173], [188, 193, 198], [], []],
+     [[], [66, 71, 76, 81], [], [], [], [], [221, 226, 231], [246, 251, 256]],
+     [[], [], [96, 97, 98], [], [156, 157], [], [217, 218], []],
+     [[], [], [], [141, 142, 143], [], [201, 202], [], [262, 263]],
+     [[], [], [], [], [158], [], [], []],
+     [[], [], [], [], [], [203], [], []],
+     [[], [], [], [], [], [], [216], []],
+     [[], [], [], [], [], [], [], [261]]]
+
+    
+     New direction: (0, 1) ,  1
+    (1, 0) 5
+    (2, 0) 10
+    (3, 0) 15
+    (4, 0) 20
+    
+     New direction: (0, -1) ,  2
+    (1, 4) 9
+    (2, 4) 14
+    (3, 4) 19
+    (4, 4) 24
+    
+     New direction: (-1, 0) ,  3
+    (5, 1) 26
+    (5, 2) 27
+    (5, 3) 28
+    
+     New direction: (1, 0) ,  4
+    (0, 1) 1
+    (0, 2) 2
+    (0, 3) 3
+    
+     New direction: (-1, 1) ,  5
+    (2, 0) 10
+    (3, 0) 15
+    (4, 0) 20
+    (5, 0) 25
+    (5, 1) 26
+    (5, 2) 27
+    
+     New direction: (1, 1) ,  6
+    (0, 0) 0
+    (0, 1) 1
+    (0, 2) 2
+    (1, 0) 5
+    (2, 0) 10
+    (3, 0) 15
+    
+     New direction: (-1, -1) ,  7
+    (2, 4) 14
+    (3, 4) 19
+    (4, 4) 24
+    (5, 2) 27
+    (5, 3) 28
+    (5, 4) 29
+    
+     New direction: (1, -1) ,  8
+    (0, 2) 2
+    (0, 3) 3
+    (0, 4) 4
+    (1, 4) 9
+    (2, 4) 14
+    (3, 4) 19
+
+    [[[35, 40, 45, 50], [], [], [], [160, 165, 170], [185, 190, 195], [], []],
+     [[], [69, 74, 79, 84], [], [], [], [], [224, 229, 234], [249, 254, 259]],
+     [[], [], [116, 117, 118], [], [176, 177], [], [237, 238], []],
+     [[], [], [], [121, 122, 123], [], [181, 182], [], [242, 243]],
+     [[], [], [], [], [175], [], [], []],
+     [[], [], [], [], [], [180], [], []],
+     [[], [], [], [], [], [], [239], []],
+     [[], [], [], [], [], [], [], [244]]]
+
+%% Cell type:code id: tags:
+
+``` python
+periodicity_mapper = SparseLbPeriodicityMapper(mapping, dh, pdf_field)
+#periodicity_mapper.create_periodic_index_array()
+
+def handle_periodicity():
+    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.],
+           [108., 109., 110., 111., 112., 113., 114., 115., 116.],
+           [117., 118., 119., 120., 121., 122., 123., 124., 125.],
+           [126., 127., 128., 129., 130., 131., 132., 133., 134.],
+           [135., 136., 137., 138., 139., 140., 141., 142., 143.],
+           [144., 145., 146., 147., 148., 149., 150., 151., 152.],
+           [153., 154., 155., 156., 157., 158., 159., 160., 161.],
+           [162., 163., 164., 165., 166., 167., 168., 169., 170.],
+           [171., 172., 173., 174., 175., 176., 177., 178., 179.],
+           [180., 181., 182., 183., 184., 185., 186., 187., 188.],
+           [189., 190., 191., 192., 193., 194., 195., 196., 197.],
+           [198., 199., 200., 201., 202., 203., 204., 205., 206.],
+           [207., 208., 209., 210., 211., 212., 213., 214., 215.],
+           [216., 217., 218., 219., 220., 221., 222., 223., 224.],
+           [225., 226., 227., 228., 229., 230., 231., 232., 233.],
+           [234., 235., 236., 237., 238., 239., 240., 241., 242.],
+           [243., 244., 245., 246., 247., 248., 249., 250., 251.],
+           [252., 253., 254., 255., 256., 257., 258., 259., 260.],
+           [261., 262., 263., 264., 265., 266., 267., 268., 269.]])
+
+%% Cell type:code id: tags:
+
+``` python
+handle_periodicity()
+pdf_arr
+```
+
+%% Output
+
+    array([[  0.,   1.,   2.,   3.,   4.,   5., 213.,   7.,   8.],
+           [  9.,  10.,  11.,  12., 193.,  14., 195.,  16.,  17.],
+           [ 18.,  19.,  20.,  21., 202.,  23., 204.,  25., 206.],
+           [ 27.,  28.,  29.,  30., 211.,  32.,  33.,  34., 215.],
+           [ 36.,  37.,  38.,  39.,  40.,  41.,  42.,  43., 197.],
+           [ 45.,  73.,  47.,  48.,  49.,  50.,  78.,  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.,  56.,  84.,  85.,  86.,  87.,  88.,  62.],
+           [ 90., 118.,  92.,  93.,  94., 122., 123.,  97.,  98.],
+           [ 99., 100., 101., 102., 103., 104., 105., 106., 107.],
+           [108., 109., 110., 111., 112., 113., 114., 115., 116.],
+           [117., 118., 119., 120., 121., 122., 123., 124., 125.],
+           [126., 127., 101., 129., 130., 131., 132., 106., 107.],
+           [135., 163., 137., 138., 139., 167., 168., 142., 143.],
+           [144., 145., 146., 147., 148., 149., 150., 151., 152.],
+           [153., 154., 155., 156., 157., 158., 159., 160., 161.],
+           [162., 163., 164., 165., 166., 167., 168., 169., 170.],
+           [171., 172., 146., 174., 175., 176., 177., 151., 152.],
+           [180., 208., 182., 183., 184., 212., 186., 187., 188.],
+           [189., 190., 191., 192., 193., 194., 195., 196., 197.],
+           [198., 199., 200., 201., 202., 203., 204., 205., 206.],
+           [207., 208., 209., 210., 211., 212., 213., 214., 215.],
+           [216., 217., 191., 219., 220., 221., 222., 196., 224.],
+           [225., 226., 227., 228., 229.,  77., 231., 232., 233.],
+           [234., 235., 236.,  57., 238.,  59., 240., 241., 242.],
+           [243., 244., 245.,  66., 247.,  68., 249.,  70., 251.],
+           [252., 253., 254.,  75., 256., 257., 258.,  79., 260.],
+           [261., 262., 263., 264., 265., 266., 267.,  61., 269.]])
+
+%% Cell type:code id: tags:
+
+``` python
+print("sparse pdf arrays:", pdf_arr.shape)
+pdf_arr.nbytes
+pdf_arr_tmp.nbytes
+
+print("fluid index array:", index_arr.shape)
+index_arr.nbytes
+
+print("periodicity index array:")
+pia = 0
+i = 0
+for i_a in periodicity_mapper._index_arrays:
+    i+=1
+    pia += i_a.nbytes
+print("Pia:", pia, i)
+pia += periodicity_mapper._index_arrays.nbytes
+periodicity_mapper._index_arrays.nbytes
+
+sparse = index_arr.nbytes + pdf_arr_tmp.nbytes + pdf_arr.nbytes + periodicity_mapper._index_arrays.nbytes + pia
+
+print("full pdf arrays:", dh.cpu_arrays[src.name].shape)
+dh.cpu_arrays[src.name].nbytes
+dh.cpu_arrays[dst.name].nbytes
+
+full = dh.cpu_arrays[dst.name].nbytes + dh.cpu_arrays[src.name].nbytes
+
+print("Compare: full", full, "sparse", sparse, "porosity 1")
+```
+
+%% Output
+
+    sparse pdf arrays: (30, 9)
+
+
+    $\displaystyle 2160$
+    2160
+
+
+    $\displaystyle 2160$
+    2160
+
+    fluid index array: (30, 8)
+
+
+    $\displaystyle 960$
+    960
+
+    periodicity index array:
+    Pia: 608 16
+
+
+    $\displaystyle 128$
+    128
+
+    full pdf arrays: (6, 5, 9)
+
+
+    $\displaystyle 2160$
+    2160
+
+
+    $\displaystyle 2160$
+    2160
+
+    Compare: full 4320 sparse 6144 porosity 1
+
+%% Cell type:code id: tags:
+
+``` python
+```

lbmpy_tests/test_sparse_lbm_peridocity_only.ipynb

+%% 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 = (2, 2, 2)
+omega = 1.8
+target = 'cpu'
+stencil = get_stencil("D3Q27")
+
+ghost_layers = 1
+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_3d(stencil)
+stencil
+```
+
+%% Output
+
+
+    $\displaystyle \left( \left( 0, \  0, \  0\right), \  \left( 0, \  1, \  0\right), \  \left( 0, \  -1, \  0\right), \  \left( -1, \  0, \  0\right), \  \left( 1, \  0, \  0\right), \  \left( 0, \  0, \  1\right), \  \left( 0, \  0, \  -1\right), \  \left( -1, \  1, \  0\right), \  \left( 1, \  1, \  0\right), \  \left( -1, \  -1, \  0\right), \  \left( 1, \  -1, \  0\right), \  \left( 0, \  1, \  1\right), \  \left( 0, \  -1, \  1\right), \  \left( -1, \  0, \  1\right), \  \left( 1, \  0, \  1\right), \  \left( 0, \  1, \  -1\right), \  \left( 0, \  -1, \  -1\right), \  \left( -1, \  0, \  -1\right), \  \left( 1, \  0, \  -1\right), \  \left( 1, \  1, \  1\right), \  \left( -1, \  1, \  1\right), \  \left( 1, \  -1, \  1\right), \  \left( -1, \  -1, \  1\right), \  \left( 1, \  1, \  -1\right), \  \left( -1, \  1, \  -1\right), \  \left( 1, \  -1, \  -1\right), \  \left( -1, \  -1, \  -1\right)\right)$
+    ((0, 0, 0), (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1
+    ), (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0), (0, 1, 1), (0, -1, 1), (-1,
+     0, 1), (1, 0, 1), (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1), (1, 1, 1)
+    , (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1), (
+    -1, -1, -1))
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+# Fix Flag field # nur für 2D :( Wird in mapping.py benutzt
+def at_border(cell):
+    if True in [cell_i == 0 or cell_i == ds_i-1 for ds_i, cell_i in zip(flag_arr.shape, cell)]:
+        return [0]
+    bool_border = [cell_i == 1 or cell_i == ds_i-2 for ds_i, cell_i in zip(flag_arr.shape, cell)]
+    return bool_border
+
+def amend_flag_arr():
+    coords = [[x, y] for y in range(flag_arr.shape[1]) for x in range(flag_arr.shape[0])]
+    res = []
+    for cell in coords:
+        bool_border = at_border(cell)
+        if flag_arr[tuple(cell)] == flags[noslip] and True in bool_border:
+            sl = [-1 if cell_i == 1 else 1 if cell_i == ds_i-2 else 0 for ds_i, cell_i in zip(flag_arr.shape, cell)]
+            l = [(cell_i + sl_i*2)%ds_i for cell_i, bb_i, ds_i, sl_i in zip(cell, bool_border, flag_arr.shape, sl)]
+            pos = np.array([cell, l])
+            for i in range(pow(2, len(l))):
+                i_bin = (format(i, '0'+str(len(l))+'b'))
+                opp = [pos[int(i_bin[arr]), arr] for arr in range(len(l))]
+    for cell in res:
+        flag_arr[cell[0], cell[1]] = flags[noslip]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Fix Flag field # nur für 2D :( Wird in mapping.py benutzt
+#flag_arr = np.zeros((5,5,5), dtype=np.uint16)
+#flag_arr.fill(1)
+#flag_arr[1,1,1] = 2
+#print(flag_arr)
+
+def at_(cell):
+    if True in [cell_i == 0 or cell_i == ds_i-1 for ds_i, cell_i in zip(flag_arr.shape, cell)]:
+        return [0]
+    bool_border = [cell_i == 1 or cell_i == ds_i-2 for ds_i, cell_i in zip(flag_arr.shape, cell)]
+    return bool_border
+
+def amend():
+    #coords = [[x, y] for y in range(flag_arr.shape[1]) for x in range(flag_arr.shape[0])]
+    #coords = [[i] for i in range(flag_arr.shape[0])]
+    coords = [[x, y, z] for z in range(flag_arr.shape[2]) for y in range(flag_arr.shape[1]) for x in range(flag_arr.shape[0])]
+    res = []
+    for cell in coords:
+        bool_border = at_(cell)
+        #print(flag_arr[tuple(cell)])
+        if flag_arr[tuple(cell)] == 2 and True in bool_border:
+            sl = [-1 if cell_i == 1 else 1 if cell_i == ds_i-2 else 0 for ds_i, cell_i in zip(flag_arr.shape, cell)]
+            l = [(cell_i + sl_i*2)%ds_i for cell_i, bb_i, ds_i, sl_i in zip(cell, bool_border, flag_arr.shape, sl)]
+            pos = np.array([cell, l])
+            for i in range(pow(2, len(l))):
+                i_bin = (format(i, '0'+str(len(l))+'b'))
+                opp = [pos[int(i_bin[arr]), arr] for arr in range(len(l))]
+                flag_arr[tuple(opp)] = 2
+
+#flag_arr.shape
+#amend()
+#print("2", flag_arr)
+```
+
+%% 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'])
+#flag_arr[3,1] = flags[noslip]
+#flag_arr[1,2,2] = flags[noslip]
+#amend_flag_arr()
+#plt.scalar_field(flag_arr)
+#plt.colorbar();
+```
+
+%% Cell type:markdown id: tags:
+
+### Mappings
+
+%% Cell type:code id: tags:
+
+``` python
+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')
+```
+
+%% Output
+
+    [[[1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]]
+    
+     [[1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]]
+    
+     [[1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]]
+    
+     [[1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]
+      [1 1 1 1]]]
+
+%% Cell type:code id: tags:
+
+``` python
+shape_pdf_arr = pdf_arr.shape
+print(shape_pdf_arr)
+
+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
+```
+
+%% Output
+
+    (64, 27)
+
+%% 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_mapper.create_periodic_index_array()
+
+def handle_periodicity():
+    periodicity_mapper(pdf_arr)
+
+```
+
+%% Cell type:code id: tags:
+
+``` python
+pdf_arr
+```
+
+%% Output
+
+    array([[0.000e+00, 1.000e+00, 2.000e+00, ..., 2.400e+01, 2.500e+01,
+            2.600e+01],
+           [2.700e+01, 2.800e+01, 2.900e+01, ..., 5.100e+01, 5.200e+01,
+            5.300e+01],
+           [5.400e+01, 5.500e+01, 5.600e+01, ..., 7.800e+01, 7.900e+01,
+            8.000e+01],
+           ...,
+           [1.647e+03, 1.648e+03, 1.649e+03, ..., 1.671e+03, 1.672e+03,
+            1.673e+03],
+           [1.674e+03, 1.675e+03, 1.676e+03, ..., 1.698e+03, 1.699e+03,
+            1.700e+03],
+           [1.701e+03, 1.702e+03, 1.703e+03, ..., 1.725e+03, 1.726e+03,
+            1.727e+03]])
+
+%% Cell type:code id: tags:
+
+``` python
+handle_periodicity()
+pdf_arr
+```
+
+%% Output
+
+    
+     New direction: (0, 1, 0) ,  1
+    (1, 0, 1) [0, 1, 0] [False, True, False]
+    (1, 0, 2) [0, 1, 0] [False, True, False]
+    (2, 0, 1) [0, 1, 0] [False, True, False]
+    (2, 0, 2) [0, 1, 0] [False, True, False]
+    
+     New direction: (0, -1, 0) ,  2
+    (1, 3, 1) [0, -1, 0] [False, True, False]
+    (1, 3, 2) [0, -1, 0] [False, True, False]
+    (2, 3, 1) [0, -1, 0] [False, True, False]
+    (2, 3, 2) [0, -1, 0] [False, True, False]
+    
+     New direction: (-1, 0, 0) ,  3
+    (3, 1, 1) [-1, 0, 0] [True, False, False]
+    (3, 1, 2) [-1, 0, 0] [True, False, False]
+    (3, 2, 1) [-1, 0, 0] [True, False, False]
+    (3, 2, 2) [-1, 0, 0] [True, False, False]
+    
+     New direction: (1, 0, 0) ,  4
+    (0, 1, 1) [1, 0, 0] [True, False, False]
+    (0, 1, 2) [1, 0, 0] [True, False, False]
+    (0, 2, 1) [1, 0, 0] [True, False, False]
+    (0, 2, 2) [1, 0, 0] [True, False, False]
+    
+     New direction: (0, 0, 1) ,  5
+    (1, 1, 0) [0, 0, 1] [False, False, True]
+    (1, 2, 0) [0, 0, 1] [False, False, True]
+    (2, 1, 0) [0, 0, 1] [False, False, True]
+    (2, 2, 0) [0, 0, 1] [False, False, True]
+    
+     New direction: (0, 0, -1) ,  6
+    (1, 1, 3) [0, 0, -1] [False, False, True]
+    (1, 2, 3) [0, 0, -1] [False, False, True]
+    (2, 1, 3) [0, 0, -1] [False, False, True]
+    (2, 2, 3) [0, 0, -1] [False, False, True]
+    
+     New direction: (-1, 1, 0) ,  7
+    (2, 0, 1) [0, 1, 0] [False, True, False]
+    (2, 0, 2) [0, 1, 0] [False, True, False]
+    (3, 0, 1) [-1, 1, 0] [True, True, False]
+    (3, 0, 2) [-1, 1, 0] [True, True, False]
+    (3, 1, 1) [-1, 0, 0] [True, False, False]
+    (3, 1, 2) [-1, 0, 0] [True, False, False]
+    
+     New direction: (1, 1, 0) ,  8
+    (0, 0, 1) [1, 1, 0] [True, True, False]
+    (0, 0, 2) [1, 1, 0] [True, True, False]
+    (0, 1, 1) [1, 0, 0] [True, False, False]
+    (0, 1, 2) [1, 0, 0] [True, False, False]
+    (1, 0, 1) [0, 1, 0] [False, True, False]
+    (1, 0, 2) [0, 1, 0] [False, True, False]
+    
+     New direction: (-1, -1, 0) ,  9
+    (2, 3, 1) [0, -1, 0] [False, True, False]
+    (2, 3, 2) [0, -1, 0] [False, True, False]
+    (3, 2, 1) [-1, 0, 0] [True, False, False]
+    (3, 2, 2) [-1, 0, 0] [True, False, False]
+    (3, 3, 1) [-1, -1, 0] [True, True, False]
+    (3, 3, 2) [-1, -1, 0] [True, True, False]
+    
+     New direction: (1, -1, 0) ,  10
+    (0, 2, 1) [1, 0, 0] [True, False, False]
+    (0, 2, 2) [1, 0, 0] [True, False, False]
+    (0, 3, 1) [1, -1, 0] [True, True, False]
+    (0, 3, 2) [1, -1, 0] [True, True, False]
+    (1, 3, 1) [0, -1, 0] [False, True, False]
+    (1, 3, 2) [0, -1, 0] [False, True, False]
+    
+     New direction: (0, 1, 1) ,  11
+    (1, 0, 0) [0, 1, 1] [False, True, True]
+    (1, 0, 1) [0, 1, 0] [False, True, False]
+    (1, 1, 0) [0, 0, 1] [False, False, True]
+    (2, 0, 0) [0, 1, 1] [False, True, True]
+    (2, 0, 1) [0, 1, 0] [False, True, False]
+    (2, 1, 0) [0, 0, 1] [False, False, True]
+    
+     New direction: (0, -1, 1) ,  12
+    (1, 2, 0) [0, 0, 1] [False, False, True]
+    (1, 3, 0) [0, -1, 1] [False, True, True]
+    (1, 3, 1) [0, -1, 0] [False, True, False]
+    (2, 2, 0) [0, 0, 1] [False, False, True]
+    (2, 3, 0) [0, -1, 1] [False, True, True]
+    (2, 3, 1) [0, -1, 0] [False, True, False]
+    
+     New direction: (-1, 0, 1) ,  13
+    (2, 1, 0) [0, 0, 1] [False, False, True]
+    (2, 2, 0) [0, 0, 1] [False, False, True]
+    (3, 1, 0) [-1, 0, 1] [True, False, True]
+    (3, 1, 1) [-1, 0, 0] [True, False, False]
+    (3, 2, 0) [-1, 0, 1] [True, False, True]
+    (3, 2, 1) [-1, 0, 0] [True, False, False]
+    
+     New direction: (1, 0, 1) ,  14
+    (0, 1, 0) [1, 0, 1] [True, False, True]
+    (0, 1, 1) [1, 0, 0] [True, False, False]
+    (0, 2, 0) [1, 0, 1] [True, False, True]
+    (0, 2, 1) [1, 0, 0] [True, False, False]
+    (1, 1, 0) [0, 0, 1] [False, False, True]
+    (1, 2, 0) [0, 0, 1] [False, False, True]
+    
+     New direction: (0, 1, -1) ,  15
+    (1, 0, 2) [0, 1, 0] [False, True, False]
+    (1, 0, 3) [0, 1, -1] [False, True, True]
+    (1, 1, 3) [0, 0, -1] [False, False, True]
+    (2, 0, 2) [0, 1, 0] [False, True, False]
+    (2, 0, 3) [0, 1, -1] [False, True, True]
+    (2, 1, 3) [0, 0, -1] [False, False, True]
+    
+     New direction: (0, -1, -1) ,  16
+    (1, 2, 3) [0, 0, -1] [False, False, True]
+    (1, 3, 2) [0, -1, 0] [False, True, False]
+    (1, 3, 3) [0, -1, -1] [False, True, True]
+    (2, 2, 3) [0, 0, -1] [False, False, True]
+    (2, 3, 2) [0, -1, 0] [False, True, False]
+    (2, 3, 3) [0, -1, -1] [False, True, True]
+    
+     New direction: (-1, 0, -1) ,  17
+    (2, 1, 3) [0, 0, -1] [False, False, True]
+    (2, 2, 3) [0, 0, -1] [False, False, True]
+    (3, 1, 2) [-1, 0, 0] [True, False, False]
+    (3, 1, 3) [-1, 0, -1] [True, False, True]
+    (3, 2, 2) [-1, 0, 0] [True, False, False]
+    (3, 2, 3) [-1, 0, -1] [True, False, True]
+    
+     New direction: (1, 0, -1) ,  18
+    (0, 1, 2) [1, 0, 0] [True, False, False]
+    (0, 1, 3) [1, 0, -1] [True, False, True]
+    (0, 2, 2) [1, 0, 0] [True, False, False]
+    (0, 2, 3) [1, 0, -1] [True, False, True]
+    (1, 1, 3) [0, 0, -1] [False, False, True]
+    (1, 2, 3) [0, 0, -1] [False, False, True]
+    
+     New direction: (1, 1, 1) ,  19
+    (0, 0, 0) [1, 1, 1] [True, True, True]
+    (0, 0, 1) [1, 1, 0] [True, True, False]
+    (0, 1, 0) [1, 0, 1] [True, False, True]
+    (0, 1, 1) [1, 0, 0] [True, False, False]
+    (1, 0, 0) [0, 1, 1] [False, True, True]
+    (1, 0, 1) [0, 1, 0] [False, True, False]
+    (1, 1, 0) [0, 0, 1] [False, False, True]
+    
+     New direction: (-1, 1, 1) ,  20
+    (2, 0, 0) [0, 1, 1] [False, True, True]
+    (2, 0, 1) [0, 1, 0] [False, True, False]
+    (2, 1, 0) [0, 0, 1] [False, False, True]
+    (3, 0, 0) [-1, 1, 1] [True, True, True]
+    (3, 0, 1) [-1, 1, 0] [True, True, False]
+    (3, 1, 0) [-1, 0, 1] [True, False, True]
+    (3, 1, 1) [-1, 0, 0] [True, False, False]
+    
+     New direction: (1, -1, 1) ,  21
+    (0, 2, 0) [1, 0, 1] [True, False, True]
+    (0, 2, 1) [1, 0, 0] [True, False, False]
+    (0, 3, 0) [1, -1, 1] [True, True, True]
+    (0, 3, 1) [1, -1, 0] [True, True, False]
+    (1, 2, 0) [0, 0, 1] [False, False, True]
+    (1, 3, 0) [0, -1, 1] [False, True, True]
+    (1, 3, 1) [0, -1, 0] [False, True, False]
+    
+     New direction: (-1, -1, 1) ,  22
+    (2, 2, 0) [0, 0, 1] [False, False, True]
+    (2, 3, 0) [0, -1, 1] [False, True, True]
+    (2, 3, 1) [0, -1, 0] [False, True, False]
+    (3, 2, 0) [-1, 0, 1] [True, False, True]
+    (3, 2, 1) [-1, 0, 0] [True, False, False]
+    (3, 3, 0) [-1, -1, 1] [True, True, True]
+    (3, 3, 1) [-1, -1, 0] [True, True, False]
+    
+     New direction: (1, 1, -1) ,  23
+    (0, 0, 2) [1, 1, 0] [True, True, False]
+    (0, 0, 3) [1, 1, -1] [True, True, True]
+    (0, 1, 2) [1, 0, 0] [True, False, False]
+    (0, 1, 3) [1, 0, -1] [True, False, True]
+    (1, 0, 2) [0, 1, 0] [False, True, False]
+    (1, 0, 3) [0, 1, -1] [False, True, True]
+    (1, 1, 3) [0, 0, -1] [False, False, True]
+    
+     New direction: (-1, 1, -1) ,  24
+    (2, 0, 2) [0, 1, 0] [False, True, False]
+    (2, 0, 3) [0, 1, -1] [False, True, True]
+    (2, 1, 3) [0, 0, -1] [False, False, True]
+    (3, 0, 2) [-1, 1, 0] [True, True, False]
+    (3, 0, 3) [-1, 1, -1] [True, True, True]
+    (3, 1, 2) [-1, 0, 0] [True, False, False]
+    (3, 1, 3) [-1, 0, -1] [True, False, True]
+    
+     New direction: (1, -1, -1) ,  25
+    (0, 2, 2) [1, 0, 0] [True, False, False]
+    (0, 2, 3) [1, 0, -1] [True, False, True]
+    (0, 3, 2) [1, -1, 0] [True, True, False]
+    (0, 3, 3) [1, -1, -1] [True, True, True]
+    (1, 2, 3) [0, 0, -1] [False, False, True]
+    (1, 3, 2) [0, -1, 0] [False, True, False]
+    (1, 3, 3) [0, -1, -1] [False, True, True]
+    
+     New direction: (-1, -1, -1) ,  26
+    (2, 2, 3) [0, 0, -1] [False, False, True]
+    (2, 3, 2) [0, -1, 0] [False, True, False]
+    (2, 3, 3) [0, -1, -1] [False, True, True]
+    (3, 2, 2) [-1, 0, 0] [True, False, False]
+    (3, 2, 3) [-1, 0, -1] [True, False, True]
+    (3, 3, 2) [-1, -1, 0] [True, True, False]
+    (3, 3, 3) [-1, -1, -1] [True, True, True]
+
+    array([[0.000e+00, 1.000e+00, 2.000e+00, ..., 2.400e+01, 2.500e+01,
+            2.600e+01],
+           [2.700e+01, 2.800e+01, 2.900e+01, ..., 5.100e+01, 5.200e+01,
+            5.300e+01],
+           [5.400e+01, 5.500e+01, 5.600e+01, ..., 7.800e+01, 7.900e+01,
+            8.000e+01],
+           ...,
+           [1.647e+03, 1.648e+03, 1.649e+03, ..., 1.671e+03, 1.672e+03,
+            1.673e+03],
+           [1.674e+03, 1.675e+03, 1.676e+03, ..., 1.698e+03, 1.699e+03,
+            6.200e+02],
+           [1.701e+03, 1.702e+03, 1.703e+03, ..., 1.725e+03, 1.726e+03,
+            5.930e+02]])
+
+%% Cell type:code id: tags:
+
+``` python
+```

lbmpy_tests/test_sparse_with_obstacles_3d.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import pytest
+pytest.importorskip('pycuda')
+```
+
+%% Output
+
+    <module 'pycuda' from '/usr/lib/python3/dist-packages/pycuda/__init__.py'>
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.sparse import *
+from pystencils.field import FieldType
+import pycuda.gpuarray as gpuarray
+
+from IPython.core.interactiveshell import InteractiveShell
+InteractiveShell.ast_node_interactivity = "all"
+import pystencils as ps
+
+import math
+```
+
+%% Cell type:markdown id: tags:
+
+# Sparse (list based) LBM
+
+%% Cell type:code id: tags:
+
+``` python
+domain_size = (200, 80, 80)
+omega = 1.99
+target = 'cpu'
+
+ghost_layers = 1
+arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
+lid_velocity = 0.01
+force = 1e-6
+
+channel = False
+stencil='D3Q27'
+
+method = create_lb_method(stencil=stencil, method='mrt', relaxation_rates=[omega, 1, 1, 1, 1, 1, 1, 1],
+                             force=(5e-6, 0, 0))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = get_stencil(stencil)
+dh = ps.create_data_handling(domain_size=domain_size, periodicity=(True, False, False))
+
+src = dh.add_array('src', values_per_cell=len(stencil), alignment=True)
+dh.fill('src', 0.0, ghost_layers=True)
+dst = dh.add_array('dst', values_per_cell=len(stencil), alignment=True)
+dh.fill('dst', 0.0, ghost_layers=True)
+
+velField = dh.add_array('velField', values_per_cell=dh.dim, alignment=True)
+dh.fill('velField', 0.0, ghost_layers=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ps.stencil.plot_3d(stencil)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+def set_sphere_xxx(x, y, z, *_):
+    c = -1
+    mid = list()
+    for i in [20, 30, 40, 50, 60, 70]:
+        c = c * -1
+        for j in [9, 21, 33]:
+            for k in [9, 21, 33]:
+                mid.append([i, j + c*8, k + c*8])
+    radius = 5
+    result = np.zeros(x.shape, dtype=bool)
+    for m in mid:
+        result += (x-m[0])**2 + (y-m[1])**2 + (z-m[1])**2 < radius**2
+    return result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+obstacles_factor = 3
+def set_sphere(x, y, z, *_):
+    c_y = 0
+    c_z = 0
+    mid = list()
+    factor = obstacles_factor
+    radius = math.ceil(domain_size[1]*0.07)
+    d = math.ceil(radius*0.8)
+    start = math.ceil(math.ceil(radius+1))
+    end = math.ceil(domain_size[0]-math.ceil(radius+1))
+    dist_x = math.ceil(1.5*radius+d+factor)
+    dist_y = math.ceil(1.5*radius+d+factor)
+    mid_y = math.ceil(domain_size[1]/2)
+    dist_z = math.ceil(1.5*radius+d+factor)
+    mid_z = math.ceil(domain_size[2]/2)
+    for i in range(start, end+1, dist_x):
+        c_y = 0 if c_y else 1
+        for k, j in enumerate(range(0, mid_y, dist_y)):
+            y_1 = (mid_y+(k+0.5)*dist_y-c_y*dist_y/2)
+            y_2 = (mid_y-(k+0.5)*dist_y-c_y*dist_y/2)
+            c_z = 0 if c_z else 1
+            for l, n in enumerate(range(0, mid_z, dist_z)):
+                z_1 = (mid_z+(l+0.5)*dist_z-c_z*dist_z/2)
+                z_2 = (mid_z-(l+0.5)*dist_z-c_z*dist_z/2)
+                mid.append([i, y_1, z_1])
+                mid.append([i, y_1, z_2])
+                mid.append([i, y_2, z_1])
+                mid.append([i, y_2, z_2])
+            #if (c == 1):
+            #    print([i, (mid_y+(k+0.5)*dist_y-c*dist_y/2)%domain_size[1]])
+            #    print([i, (mid_y-(k+0.5)*dist_y-c*dist_y/2)%domain_size[1]])
+        #for j in range(radius+d, domain_size[1]-radius-d+1, dist_y):
+            #mid.append([i, math.ceil(j + c*radius*0.8)])
+    result = np.zeros(x.shape, dtype=bool)
+    for m in mid:
+        result += (x-m[0])**2 + (y-m[1])**2 + (z-m[2])**2 < radius**2
+    return result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+x = np.arange(0, arr_size[0], 1)
+y = np.arange(0, arr_size[1], 1)
+z = np.arange(0, arr_size[2], 1)
+xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")
+
+bool_index = set_sphere(xx, yy, zz)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ubb = UBB( (lid_velocity, 0, 0) )
+noslip = NoSlip()
+flags = {
+    'fluid': 1,
+    noslip: 2,
+}
+flag_arr = np.zeros(arr_size, dtype=np.uint16)
+flag_arr.fill(flags['fluid'])
+flag_arr[bool_index] = flags[noslip]
+
+flag_arr[:, -1, :] = flags[noslip]
+flag_arr[:, 0, :] = flags[noslip]
+flag_arr[:, :, -1] = flags[noslip]
+flag_arr[:, :, 0] = flags[noslip]
+flag_arr.shape
+fl = 20
+ft = 20
+plt.xlabel("x", fontsize=fl)
+plt.ylabel("y", fontsize=fl)
+plt.xticks(fontsize=ft)
+plt.yticks(fontsize=ft)
+#plt.pcolormesh(flag_arr, cmap=plt.cm.get_cmap('RdBu'))
+#nipy_spectral_r
+plt.scalar_field(flag_arr[:, :, 35], cmap=plt.cm.get_cmap('nipy_spectral_r'))
+#plt.colorbar();
+plt.savefig("yyy.png")
+```
+
+%% Output
+
+
+    $\displaystyle \left( 202, \  82, \  82\right)$
+    (202, 82, 82)
+
+    Text(0.5, 0, 'x')
+
+    Text(0, 0.5, 'y')
+
+    (array([0. , 0.2, 0.4, 0.6, 0.8, 1. ]), <a list of 6 Text xticklabel objects>)
+
+    (array([0. , 0.2, 0.4, 0.6, 0.8, 1. ]), <a list of 6 Text yticklabel objects>)
+
+    <matplotlib.image.AxesImage at 0x7f1bb20fc8e0>
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+#Get porosity
+def porosity():
+    print("All", flag_arr.size)
+    print("Fluid", flag_arr[flag_arr == flags['fluid']].size)
+    print("Solid", flag_arr[flag_arr == flags[noslip]].size)
+    return flag_arr[flag_arr == flags['fluid']].size/flag_arr.size
+
+porosity()
+```
+
+%% Output
+
+    All 1358248
+    Fluid 1073969
+    Solid 284279
+
+
+    $\displaystyle 0.7907016980698665$
+    0.7907016980698665
+
+%% Cell type:markdown id: tags:
+
+### Mappings
+
+%% Cell type:code id: tags:
+
+``` python
+mapping = SparseLbMapper(method.stencil, domain_size, flag_arr, flags['fluid'], flags[noslip], 0)
+index_arr = mapping.create_index_array(ghost_layers)
+index_arr.shape
+# 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(method.stencil)), order='f')
+pdf_arr_tmp = np.empty_like(pdf_arr)
+
+vel_arr = np.ones((len(mapping), len(domain_size)))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+pdf_field, pdf_field_tmp, vel_field = ps.fields("f(27), d(27), u(3): [1D]",
+                                                #f=pdf_arr[:mapping.num_fluid_cells],
+                                                #d=pdf_arr_tmp[:mapping.num_fluid_cells],
+                                                #u=vel_arr
+                                               )
+pdf_field.field_type = FieldType.CUSTOM
+pdf_field.pdf_field_tmp = FieldType.CUSTOM
+```
+
+%% Cell type:code id: tags:
+
+``` python
+method.first_order_equilibrium_moment_symbols
+```
+
+%% Cell type:markdown id: tags:
+
+### Macroscopic quantities
+
+%% Cell type:code id: tags:
+
+``` python
+cqc = method.conserved_quantity_computation
+
+inp_eqs = cqc.equilibrium_input_equations_from_init_values()
+setter_eqs = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
+setter_eqs = setter_eqs.new_with_substitutions({sym: pdf_field(i)
+                                                for i, sym in enumerate(method.post_collision_pdf_symbols)})
+ur = setter_eqs.all_assignments
+ur[1] = ps.Assignment(method.first_order_equilibrium_moment_symbols[0], vel_field(0))
+ur[2] = ps.Assignment(method.first_order_equilibrium_moment_symbols[1], vel_field(1))
+ur[3] = ps.Assignment(method.first_order_equilibrium_moment_symbols[2], vel_field(2))
+kernel_initialize = ps.create_kernel(ur, ghost_layers=((0, 0, 0),), ).compile()
+
+def init():
+    vel_arr = mapping.set_velocity(dh.cpu_arrays[velField.name])
+    kernel_initialize(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
+init()
+
+vel_arr.shape # Enthält auch Ghost Zellen. (sparse)
+
+getter_eqs = cqc.output_equations_from_pdfs(pdf_field.center_vector,
+                                            {'velocity': vel_field})
+kernel_compute_u = ps.create_kernel(getter_eqs, ghost_layers=((0,0, 0),) ).compile()
+
+def get_velocity(arr=pdf_arr, ghost_layers=True): #macht aus sparse vel_arr einen full_velocity_arr
+    fluid_cell_arr = mapping.coordinates
+    kernel_compute_u(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
+    full_velocity_arr = np.zeros(flag_arr.shape + (3,), dtype=np.float64)
+    arr = fluid_cell_arr[:mapping.num_fluid_cells]
+    full_velocity_arr[arr['x'], arr['y'], arr['z']] = vel_arr
+    #np.nan_to_num(full_velocity_arr, copy=False, nan=0.0, posinf=None, neginf=None) #PROBLEM: bei no-slip kommt nan raus...
+    if ghost_layers:
+        return full_velocity_arr[1:-1, 1:-1, domain_size[2]//2]
+    else:
+        return full_velocity_arr[1:-1, 1:-1, domain_size[2]//2]
+```
+
+%% Cell type:markdown id: tags:
+
+### Stream Collide Kernel
+
+%% Cell type:code id: tags:
+
+``` python
+#index_field = ps.Field.create_from_numpy_array("idx", index_arr, index_dimensions=1)
+index_field = ps.Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=index_arr.dtype)
+
+collision_rule = method.get_collision_rule()
+
+Q = len(method.stencil)
+symbol_subs = {sym: pdf_field.absolute_access((index_field(i-1),),())
+               for i, sym in enumerate(method.pre_collision_pdf_symbols)}
+
+symbol_subs.update({sym: pdf_field_tmp(i) for i, sym in enumerate(method.post_collision_pdf_symbols)})
+
+symbol_subs[method.pre_collision_pdf_symbols[0]] = pdf_field(0) # special case for center
+```
+
+%% Cell type:code id: tags:
+
+``` python
+collision_rule = method.get_collision_rule()
+update_rule = collision_rule.new_with_substitutions(symbol_subs)
+kernel_stream_collide = ps.create_kernel(update_rule, ghost_layers=[(0,0,0)], target=target).compile()
+```
+
+%% Cell type:markdown id: tags:
+
+### Boundary Kernels
+
+%% Cell type:code id: tags:
+
+``` python
+if not channel:
+    if target == 'gpu':
+        raise NotImplementedError("UBB on GPU not working yet")
+
+    #ubb_mapper = SparseLbBoundaryMapper(ubb, method, pdf_field)
+    #TODO the following line is wrong: kernel contains accesses to index_field and pdf_field which have
+    #different size: a correct kernel comes out when by change the shape is taken from index field,
+    # when taken from pdf field, a wrong kernel is generated
+    #ubb_ast = ubb_mapper.get_kernel()
+
+    #ps.show_code(ubb_ast)
+
+    #ubb_kernel = ubb_ast.compile()
+    #ubb_idx_arr = ubb_mapper.create_index_arr(mapping, flags[ubb])
+    #def handle_ubb():
+    #    ubb_kernel(f=pdf_arr[:mapping.num_fluid_cells], idx=ubb_idx_arr)
+
+    periodicity_mapper = SparseLbPeriodicityMapper(mapping, dh, pdf_field)
+
+    def handle_periodicity():
+        periodicity_mapper(pdf_arr)
+
+else:
+    def handle_ubb():
+        pass
+```
+
+%% Cell type:markdown id: tags:
+
+### Time Loop
+
+%% Cell type:code id: tags:
+
+``` python
+def time_step():
+    global pdf_arr, pdf_arr_tmp, index_arr
+    handle_periodicity()
+    if target == 'gpu':
+        gpu_pdf_arr = gpuarray.to_gpu(pdf_arr)
+        gpu_pdf_arr_tmp = gpuarray.to_gpu(pdf_arr_tmp)
+        gpu_index_arr = gpuarray.to_gpu(index_arr)
+
+        kernel_stream_collide(f=gpu_pdf_arr[:mapping.num_fluid_cells],
+                          d=gpu_pdf_arr_tmp[:mapping.num_fluid_cells],
+                          idx=gpu_index_arr)
+
+        pdf_arr = gpu_pdf_arr.get()
+        pdf_arr_tmp = gpu_pdf_arr_tmp.get()
+        index_arr = gpu_index_arr.get()
+    else:
+        kernel_stream_collide(f=pdf_arr[:mapping.num_fluid_cells],
+                              d=pdf_arr_tmp[:mapping.num_fluid_cells],
+                              idx=index_arr)
+    pdf_arr_tmp, pdf_arr = pdf_arr, pdf_arr_tmp
+    #handle_periodicity()
+    #handle_ubb()
+
+def run(steps=100):
+    for t in range(steps):
+        time_step()
+    return get_velocity()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+init()
+
+mask = np.fromfunction(set_sphere, (domain_size[0], domain_size[1], domain_size[2], 3))
+print("mask", mask[:, :, 35, :].shape, )
+def create_animation():
+    result = run(50)
+    result = np.ma.array(result, mask=mask[:, :, 35, :])
+    return result
+
+animation = plt.vector_field_magnitude_animation(create_animation, frames=100, rescale=True)
+set_display_mode('video')
+res = display_animation(animation)
+res
+```
+
+%% Cell type:code id: tags:
+
+``` python
+#Get porosity
+def porosity():
+    print("All", mapping._flag_arr.size)
+    print("Fluid", mapping._flag_arr[flag_arr == flags['fluid']].size)
+    print("Solid", mapping._flag_arr[flag_arr == flags[noslip]].size)
+    return mapping._flag_arr[flag_arr == flags['fluid']].size/mapping._flag_arr.size
+
+porosity()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+fluid = str(flag_arr[flag_arr == flags['fluid']].size)
+solid = str(flag_arr[flag_arr == flags[noslip]].size)
+full = str(flag_arr.size)
+phi = str(porosity())
+```
+
+%% Cell type:code id: tags:
+
+``` python
+print("sparse pdf arrays:", pdf_arr.shape, pdf_arr.nbytes, pdf_arr_tmp.nbytes)
+m_pdf = str(pdf_arr.nbytes)
+
+print("fluid index array:", index_arr.shape, index_arr.nbytes)
+m_fia = str(index_arr.nbytes)
+
+print("periodicity index array:", periodicity_mapper._index_arrays.nbytes)
+m_pia = str(periodicity_mapper._index_arrays.nbytes)
+
+sparse = index_arr.nbytes + pdf_arr_tmp.nbytes + pdf_arr.nbytes + periodicity_mapper._index_arrays.nbytes
+m_i = str(sparse)
+
+print("full pdf arrays:", dh.cpu_arrays[src.name].shape, dh.cpu_arrays[src.name].nbytes)
+full = dh.cpu_arrays[dst.name].nbytes + dh.cpu_arrays[src.name].nbytes
+m_d = str(full)
+
+print("Compare: full", full, "sparse", sparse, "porosity", porosity())
+ratio = str(sparse/full)
+#f = open("../../tests/data.txt", "a")
+#f.write(phi+"\t\t"+fluid+"\t"+solid+"\t"+str(obstacles_factor)+"\t\t"+m_d+"\t\t"+m_i+"\t\t"+m_pdf+"\t\t"+m_fia+"\t\t"+m_pia+"\t"+ratio+"\n")
+#f.close()
+```
+
+%% Cell type:markdown id: tags:
+
+### Check against reference
+
+%% Cell type:code id: tags:
+
+``` python
+if channel:
+    reference = create_channel(domain_size, force=force, lb_method=method)
+else:
+    reference = create_lid_driven_cavity(domain_size, relaxation_rate=omega, lid_velocity=lid_velocity,
+                                         compressible=False)
+reference.run(timesteps)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+np.testing.assert_almost_equal(reference.velocity[:, :], result)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.vector_field(reference.velocity[:, :], step=1)
+```

lbmpy_tests/test_sparse_with_obstacles_lbm.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import pytest
+pytest.importorskip('pycuda')
+```
+
+%% Output
+
+    <module 'pycuda' from '/usr/lib/python3/dist-packages/pycuda/__init__.py'>
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.sparse import *
+from pystencils.field import FieldType
+import pycuda.gpuarray as gpuarray
+
+from IPython.core.interactiveshell import InteractiveShell
+InteractiveShell.ast_node_interactivity = "all"
+import pystencils as ps
+
+import math
+```
+
+%% Cell type:markdown id: tags:
+
+# Sparse (list based) LBM
+
+%% Cell type:code id: tags:
+
+``` python
+domain_size = (400,140)
+omega = 1.99
+target = 'cpu'
+
+ghost_layers = 1
+arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
+lid_velocity = 0.01
+force = 1e-6
+
+channel = False
+stencil='D2Q9'
+
+method = create_lb_method(stencil=stencil, method='mrt', relaxation_rates=[omega, 1, 1, 1],
+                             force=(5e-6, 0))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = get_stencil("D2Q9")
+dh = ps.create_data_handling(domain_size=domain_size, periodicity=(False, False))
+
+src = dh.add_array('src', values_per_cell=len(stencil), alignment=True)
+dh.fill('src', 0.0, ghost_layers=True)
+dst = dh.add_array('dst', values_per_cell=len(stencil), alignment=True)
+dh.fill('dst', 0.0, ghost_layers=True)
+
+velField = dh.add_array('velField', values_per_cell=dh.dim, alignment=True)
+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
+def set_sphere_exact(x, y, *_):
+    c = -1
+    mid = list()
+    for i in [80, 110, 140, 170, 200, 230, 260]:
+        c = c * -1
+        for j in [18, 46, 74, 102]:
+            mid.append([i, j + c*8])
+    radius = 10
+    result = np.zeros(x.shape, dtype=bool)
+    for m in mid:
+        result += (x-m[0])**2 + (y-m[1])**2 < radius**2
+    return result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+obstacles_factor = 10
+def set_sphere(x, y, *_):
+    c = 0
+    mid = list()
+    factor = obstacles_factor
+    radius = math.ceil(domain_size[1]*1/14*0.7)
+    d = math.ceil(radius*0.8)
+    start = math.ceil(math.ceil(radius+1))
+    end = math.ceil(domain_size[0]-math.ceil(radius+1))
+    dist_x = math.ceil(1.5*radius+d+factor)
+    dist_y = math.ceil(1.5*radius+d+factor)
+    mid_y = math.ceil(domain_size[1]/2)
+    for i in range(start, end+1, dist_x):
+        c = 0 if c else 1
+        for k, j in enumerate(range(0, mid_y, dist_y)):
+            a = (mid_y+(k+0.5)*dist_y-c*dist_y/2)
+            b = (mid_y-(k+0.5)*dist_y-c*dist_y/2)
+            #if (a > domain_size[1]-1):
+            #    mid.append([i, a%domain_size[1]])
+            #if (b < 0 or b > domain_size[1]-1):
+            #    mid.append([i, b%domain_size[1]])
+            mid.append([i, a])
+            mid.append([i, b])
+            #if (c == 1):
+            #    print([i, (mid_y+(k+0.5)*dist_y-c*dist_y/2)%domain_size[1]])
+            #    print([i, (mid_y-(k+0.5)*dist_y-c*dist_y/2)%domain_size[1]])
+        #for j in range(radius+d, domain_size[1]-radius-d+1, dist_y):
+            #mid.append([i, math.ceil(j + c*radius*0.8)])
+    result = np.zeros(x.shape, dtype=bool)
+    for m in mid:
+        result += (x-m[0])**2 + (y-m[1])**2 < radius**2
+    return result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+x = np.arange(0, arr_size[0], 1)
+y = np.arange(0, arr_size[1], 1)
+xx, yy = np.meshgrid(x, y, indexing="ij")
+
+bool_index = set_sphere(xx, yy)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ubb = UBB( (lid_velocity, 0) )
+noslip = NoSlip()
+flags = {
+    'fluid': 1,
+    noslip: 2,
+}
+flag_arr = np.zeros(arr_size, dtype=np.uint16)
+flag_arr.fill(flags['fluid'])
+flag_arr[bool_index] = flags[noslip]
+
+flag_arr[:, -1] = flags[noslip]
+flag_arr[:, 0] = flags[noslip]
+fl = 20
+ft = 20
+plt.xlabel("x", fontsize=fl)
+plt.ylabel("y", fontsize=fl)
+plt.xticks(fontsize=ft)
+plt.yticks(fontsize=ft)
+#plt.pcolormesh(flag_arr, cmap=plt.cm.get_cmap('RdBu'))
+#nipy_spectral_r
+plt.scalar_field(flag_arr, cmap=plt.cm.get_cmap('nipy_spectral_r'))
+#plt.colorbar();
+plt.savefig("../../my_pngs/flag_field.png")
+```
+
+%% Output
+
+    Text(0.5, 0, 'x')
+
+    Text(0, 0.5, 'y')
+
+    (array([0. , 0.2, 0.4, 0.6, 0.8, 1. ]), <a list of 6 Text xticklabel objects>)
+
+    (array([0. , 0.2, 0.4, 0.6, 0.8, 1. ]), <a list of 6 Text yticklabel objects>)
+
+    <matplotlib.image.AxesImage at 0x7f69e1e2b1f0>
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+#Get porosity
+def porosity():
+    print("All", flag_arr.size)
+    print("Fluid", flag_arr[flag_arr == flags['fluid']].size)
+    print("Solid", flag_arr[flag_arr == flags[noslip]].size)
+    return flag_arr[flag_arr == flags['fluid']].size/flag_arr.size
+
+porosity()
+```
+
+%% Output
+
+    All 57084
+    Fluid 44635
+    Solid 12449
+
+
+    $\displaystyle 0.781917875411674$
+    0.781917875411674
+
+%% Cell type:markdown id: tags:
+
+### Mappings
+
+%% Cell type:code id: tags:
+
+``` python
+mapping = SparseLbMapper(method.stencil, domain_size, flag_arr, flags['fluid'], flags[noslip], 0)
+index_arr = mapping.create_index_array(ghost_layers)
+index_arr.shape
+# 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(method.stencil)), order='f')
+pdf_arr_tmp = np.empty_like(pdf_arr)
+
+vel_arr = np.ones((len(mapping), len(domain_size)))
+```
+
+%% Output
+
+
+    $\displaystyle \left( 44635, \  8\right)$
+    (44635, 8)
+
+%% Cell type:code id: tags:
+
+``` python
+pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]",
+                                                #f=pdf_arr[:mapping.num_fluid_cells],
+                                                #d=pdf_arr_tmp[:mapping.num_fluid_cells],
+                                                #u=vel_arr
+                                               )
+pdf_field.field_type = FieldType.CUSTOM
+pdf_field.pdf_field_tmp = FieldType.CUSTOM
+```
+
+%% Cell type:code id: tags:
+
+``` python
+method.first_order_equilibrium_moment_symbols
+```
+
+%% Output
+
+
+    $\displaystyle \left( u_{0}, \  u_{1}\right)$
+    (u₀, u₁)
+
+%% Cell type:markdown id: tags:
+
+### Macroscopic quantities
+
+%% Cell type:code id: tags:
+
+``` python
+cqc = method.conserved_quantity_computation
+
+inp_eqs = cqc.equilibrium_input_equations_from_init_values()
+setter_eqs = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
+setter_eqs = setter_eqs.new_with_substitutions({sym: pdf_field(i)
+                                                for i, sym in enumerate(method.post_collision_pdf_symbols)})
+ur = setter_eqs.all_assignments
+ur[1] = ps.Assignment(method.first_order_equilibrium_moment_symbols[0], vel_field(0))
+ur[2] = ps.Assignment(method.first_order_equilibrium_moment_symbols[1], vel_field(1))
+kernel_initialize = ps.create_kernel(ur, ghost_layers=((0, 0),), ).compile()
+
+def init():
+    vel_arr = mapping.set_velocity(dh.cpu_arrays[velField.name])
+    kernel_initialize(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
+init()
+
+vel_arr.shape # Enthält auch Ghost Zellen. (sparse)
+
+getter_eqs = cqc.output_equations_from_pdfs(pdf_field.center_vector,
+                                            {'velocity': vel_field})
+kernel_compute_u = ps.create_kernel(getter_eqs, ghost_layers=((0,0),) ).compile()
+
+def get_velocity(arr=pdf_arr, ghost_layers=True): #macht aus sparse vel_arr einen full_velocity_arr
+    fluid_cell_arr = mapping.coordinates
+    kernel_compute_u(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
+    full_velocity_arr = np.zeros(flag_arr.shape + (2,), dtype=np.float64)
+    arr = fluid_cell_arr[:mapping.num_fluid_cells]
+    full_velocity_arr[arr['x'], arr['y']] = vel_arr
+    #np.nan_to_num(full_velocity_arr, copy=False, nan=0.0, posinf=None, neginf=None) #PROBLEM: bei no-slip kommt nan raus...
+    if ghost_layers:
+        return full_velocity_arr[1:-1, 1:-1]
+    else:
+        return full_velocity_arr[1:-1, 1:-1]
+```
+
+%% Output
+
+
+    $\displaystyle \left( 44635, \  2\right)$
+    (44635, 2)
+
+%% Cell type:markdown id: tags:
+
+### Stream Collide Kernel
+
+%% Cell type:code id: tags:
+
+``` python
+#index_field = ps.Field.create_from_numpy_array("idx", index_arr, index_dimensions=1)
+index_field = ps.Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=index_arr.dtype)
+
+collision_rule = method.get_collision_rule()
+
+Q = len(method.stencil)
+symbol_subs = {sym: pdf_field.absolute_access((index_field(i-1),),())
+               for i, sym in enumerate(method.pre_collision_pdf_symbols)}
+
+symbol_subs.update({sym: pdf_field_tmp(i) for i, sym in enumerate(method.post_collision_pdf_symbols)})
+
+symbol_subs[method.pre_collision_pdf_symbols[0]] = pdf_field(0) # special case for center
+```
+
+%% Cell type:code id: tags:
+
+``` python
+collision_rule = method.get_collision_rule()
+update_rule = collision_rule.new_with_substitutions(symbol_subs)
+kernel_stream_collide = ps.create_kernel(update_rule, ghost_layers=[(0,0)], target=target).compile()
+```
+
+%% Cell type:markdown id: tags:
+
+### Boundary Kernels
+
+%% Cell type:code id: tags:
+
+``` python
+if not channel:
+    if target == 'gpu':
+        raise NotImplementedError("UBB on GPU not working yet")
+
+    #ubb_mapper = SparseLbBoundaryMapper(ubb, method, pdf_field)
+    #TODO the following line is wrong: kernel contains accesses to index_field and pdf_field which have
+    #different size: a correct kernel comes out when by change the shape is taken from index field,
+    # when taken from pdf field, a wrong kernel is generated
+    #ubb_ast = ubb_mapper.get_kernel()
+
+    #ps.show_code(ubb_ast)
+
+    #ubb_kernel = ubb_ast.compile()
+    #ubb_idx_arr = ubb_mapper.create_index_arr(mapping, flags[ubb])
+    #def handle_ubb():
+    #    ubb_kernel(f=pdf_arr[:mapping.num_fluid_cells], idx=ubb_idx_arr)
+
+    periodicity_mapper = SparseLbPeriodicityMapper(mapping, dh, pdf_field)
+
+    def handle_periodicity():
+        periodicity_mapper(pdf_arr)
+
+else:
+    def handle_ubb():
+        pass
+```
+
+%% Cell type:markdown id: tags:
+
+### Time Loop
+
+%% Cell type:code id: tags:
+
+``` python
+def time_step():
+    global pdf_arr, pdf_arr_tmp, index_arr
+    handle_periodicity()
+    if target == 'gpu':
+        gpu_pdf_arr = gpuarray.to_gpu(pdf_arr)
+        gpu_pdf_arr_tmp = gpuarray.to_gpu(pdf_arr_tmp)
+        gpu_index_arr = gpuarray.to_gpu(index_arr)
+
+        kernel_stream_collide(f=gpu_pdf_arr[:mapping.num_fluid_cells],
+                          d=gpu_pdf_arr_tmp[:mapping.num_fluid_cells],
+                          idx=gpu_index_arr)
+
+        pdf_arr = gpu_pdf_arr.get()
+        pdf_arr_tmp = gpu_pdf_arr_tmp.get()
+        index_arr = gpu_index_arr.get()
+    else:
+        kernel_stream_collide(f=pdf_arr[:mapping.num_fluid_cells],
+                              d=pdf_arr_tmp[:mapping.num_fluid_cells],
+                              idx=index_arr)
+    pdf_arr_tmp, pdf_arr = pdf_arr, pdf_arr_tmp
+    #handle_periodicity()
+    #handle_ubb()
+
+def run(steps=100):
+    for t in range(steps):
+        time_step()
+    return get_velocity()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+init()
+
+mask = np.fromfunction(set_sphere, (domain_size[0], domain_size[1], 2))
+mask.shape
+
+def create_animation():
+    result = run(2)
+    result = np.ma.array(result, mask=mask)
+    return result
+
+animation = plt.vector_field_magnitude_animation(create_animation, frames=1, rescale=True)
+set_display_mode('video')
+res = display_animation(animation)
+res
+```
+
+%% Output
+
+
+    $\displaystyle \left( 400, \  140, \  2\right)$
+    (400, 140, 2)
+
+    <IPython.core.display.HTML object>
+
+%% Cell type:code id: tags:
+
+``` python
+#Get porosity
+def porosity():
+    print("All", mapping._flag_arr.size)
+    print("Fluid", mapping._flag_arr[flag_arr == flags['fluid']].size)
+    print("Solid", mapping._flag_arr[flag_arr == flags[noslip]].size)
+    return mapping._flag_arr[flag_arr == flags['fluid']].size/mapping._flag_arr.size
+
+porosity()
+```
+
+%% Output
+
+    All 57084
+    Fluid 44635
+    Solid 12449
+
+
+    $\displaystyle 0.781917875411674$
+    0.781917875411674
+
+%% Cell type:code id: tags:
+
+``` python
+fluid = str(flag_arr[flag_arr == flags['fluid']].size)
+solid = str(flag_arr[flag_arr == flags[noslip]].size)
+full = str(flag_arr.size)
+phi = str(porosity())
+```
+
+%% Output
+
+    All 57084
+    Fluid 44635
+    Solid 12449
+
+%% Cell type:code id: tags:
+
+``` python
+print("sparse pdf arrays:", pdf_arr.shape, pdf_arr.nbytes, pdf_arr_tmp.nbytes)
+m_pdf = str(pdf_arr.nbytes)
+
+print("fluid index array:", index_arr.shape, index_arr.nbytes)
+m_fia = str(index_arr.nbytes)
+
+print("periodicity index array:", periodicity_mapper._index_arrays.nbytes)
+m_pia = str(periodicity_mapper._index_arrays.nbytes)
+
+sparse = index_arr.nbytes + pdf_arr_tmp.nbytes + pdf_arr.nbytes + periodicity_mapper._index_arrays.nbytes
+m_i = str(sparse)
+
+print("full pdf arrays:", dh.cpu_arrays[src.name].shape, dh.cpu_arrays[src.name].nbytes)
+full = dh.cpu_arrays[dst.name].nbytes + dh.cpu_arrays[src.name].nbytes
+m_d = str(full)
+
+print("Compare: full", full, "sparse", sparse, "porosity", porosity())
+ratio = str(sparse/full)
+f = open("../../tests/data.txt", "a")
+f.write(phi+"\t\t"+fluid+"\t"+solid+"\t"+str(obstacles_factor)+"\t\t"+m_d+"\t\t"+m_i+"\t\t"+m_pdf+"\t\t"+m_fia+"\t\t"+m_pia+"\t"+ratio+"\n")
+f.close()
+```
+
+%% Output
+
+    sparse pdf arrays: (44635, 9) 3213720 3213720
+    fluid index array: (44635, 8) 1428320
+    periodicity index array: 6688
+    full pdf arrays: (402, 142, 9) 4110048
+    All 57084
+    Fluid 44635
+    Solid 12449
+    Compare: full 8220096 sparse 7862448 porosity 0.781917875411674
+
+
+    $\displaystyle 95$
+    95
+
+%% Cell type:markdown id: tags:
+
+### Check against reference
+
+%% Cell type:code id: tags:
+
+``` python
+if channel:
+    reference = create_channel(domain_size, force=force, lb_method=method)
+else:
+    reference = create_lid_driven_cavity(domain_size, relaxation_rate=omega, lid_velocity=lid_velocity,
+                                         compressible=False)
+reference.run(timesteps)
+```
+
+%% Output
+
+    ---------------------------------------------------------------------------
+    NameError                                 Traceback (most recent call last)
+    <ipython-input-23-580a42c05a1e> in <module>
+          4     reference = create_lid_driven_cavity(domain_size, relaxation_rate=omega, lid_velocity=lid_velocity,
+          5                                          compressible=False)
+    ----> 6 reference.run(timesteps)
+
+    NameError: name 'timesteps' is not defined
+
+%% Cell type:code id: tags:
+
+``` python
+np.testing.assert_almost_equal(reference.velocity[:, :], result)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.vector_field(reference.velocity[:, :], step=1)
+```