.gitignore

@@ -10,4 +10,5 @@ __pycache__
 _build
 /.idea
 .cache
-_local_tmp
\ No newline at end of file
+_local_tmp
+**/.vscode
\ No newline at end of file

lbmpy/advanced_streaming/__init__.py

+from .boundaryindexing import AdvancedStreamingBoundaryIndexing
+from .boundaryconditions import FlexibleBoundary, FlexibleNoSlip
+from .boundaryhandling import FlexibleLBMBoundaryHandling, \
+    create_advanced_streaming_boundary_kernel
+from .communication import get_communication_slices, PeriodicityHandling
+from .utility import Timestep, get_accessor
+
+__all__ = ['AdvancedStreamingBoundaryIndexing', 'FlexibleBoundary',
+           'FlexibleNoSlip', 'create_advanced_streaming_boundary_kernel',
+           'FlexibleLBMBoundaryHandling',
+           'get_communication_slices', 'PeriodicityHandling',
+           'Timestep', 'get_accessor']

lbmpy/advanced_streaming/boundaryconditions.py

+import pystencils as ps
+
+
+class FlexibleBoundary:
+    """Base class for all boundary classes that use the AdvancedStreamingBoundaryIndexing"""
+
+    inner_or_boundary = True
+    single_link = False
+
+    def __init__(self, name=None):
+        self._name = name
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        """
+        This function defines the boundary behavior and must therefore be implemented by all boundaries.
+        The boundary is defined through a list of sympy equations from which a boundary kernel is generated.
+
+
+        :param f_out: a pystencils field acting as a proxy to access the populations streaming out of the current
+                cell, i.e. the post-collision PDFs of the previous LBM step
+        :param f_in: a pystencils field acting as a proxy to access the populations streaming into the current
+                cell, i.e. the pre-collision PDFs for the next LBM step
+        :param dir_symbol: a sympy symbol that can be used as an index to f_out and f_in. It describes the direction
+                    pointing from the fluid to the boundary cell.
+        :param inv_dir: an indexed sympy symbol which describes the inversion of a direction index. It can be used in
+                the indices of f_out and f_in for retrieving the PDF of the inverse direction.
+        :param lb_method: an instance of the LB method used. Use this to adapt the boundary to the method
+                    (e.g. compressibility)
+        :param index_field: the boundary index field that can be used to retrieve and update boundary data
+
+        :return: list of pystencils assignments, or pystencils.AssignmentCollection
+        """
+        raise NotImplementedError("Boundary class has to overwrite __call__")
+
+    @property
+    def additional_data(self):
+        """Return a list of (name, type) tuples for additional data items required in this boundary
+        These data items can either be initialized in separate kernel see additional_data_kernel_init or by
+        Python callbacks - see additional_data_callback """
+        return []
+
+    @property
+    def additional_data_init_callback(self):
+        """Return a callback function called with a boundary data setter object and returning a dict of
+        data-name to data for each element that should be initialized"""
+        return None
+
+    @property
+    def additional_code_nodes(self):
+        """Return a list of code nodes that will be added in the generated code before the index field loop."""
+        return []
+
+    @property
+    def name(self):
+        if self._name:
+            return self._name
+        else:
+            return type(self).__name__
+
+    @name.setter
+    def name(self, new_value):
+        self._name = new_value
+
+# end class FlexibleBoundary
+
+
+class FlexibleNoSlip(FlexibleBoundary):
+
+    def __init__(self, name=None):
+        """Set an optional name here, to mark boundaries, for example for force evaluations"""
+        super(FlexibleNoSlip, self).__init__(name)
+
+    """
+        No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
+        Extended for use with any streaming pattern.
+    """
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        return ps.Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
+
+    def __hash__(self):
+        return hash(self.name)
+
+    def __eq__(self, other):
+        if not isinstance(other, FlexibleNoSlip):
+            return False
+        return self.name == other.name
+
+# end class FlexibleNoSlip

lbmpy/advanced_streaming/boundaryhandling.py

+from lbmpy.advanced_streaming.boundaryindexing import AdvancedStreamingBoundaryIndexing
+from lbmpy.advanced_streaming.utility import Timestep
+from pystencils import Field, Assignment, create_indexed_kernel
+from pystencils.boundaries import BoundaryHandling
+from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
+
+
+class FlexibleLBMBoundaryHandling(BoundaryHandling):
+    """
+    Enables boundary handling for LBM simulations with advanced streaming patterns. 
+    For the in-place patterns AA and EsoTwist, two kernels are generated for a boundary 
+    object and the right one selected depending on the time step.
+    """
+
+    def __init__(self, lb_method, data_handling, pdf_field_name, streaming_pattern='pull',
+                 name="boundary_handling", flag_interface=None, target='cpu', openmp=True):
+        self._lb_method = lb_method
+        self._streaming_pattern = streaming_pattern
+        self._two_fields_kernel = streaming_pattern in ['pull', 'push']
+        self._prev_timestep = None
+        super(FlexibleLBMBoundaryHandling, self).__init__(data_handling, pdf_field_name, lb_method.stencil,
+                                                          name, flag_interface, target, openmp)
+
+    #   TODO: Force On Boundary
+
+    #   ------------------------- Overridden methods of pystencils.BoundaryHandling -------------------------
+
+    @property
+    def prev_timestep(self):
+        return self._prev_timestep
+
+    def __call__(self, prev_timestep=Timestep.BOTH, **kwargs):
+        self._prev_timestep = prev_timestep
+        super(FlexibleLBMBoundaryHandling, self).__call__(**kwargs)
+        self._prev_timestep = None
+
+    def add_fixed_steps(self, fixed_loop, **kwargs):
+        raise NotImplementedError("Adding to fixed loop is not supported by FlexibleLBMBoundaryHandling")
+
+    def _add_boundary(self, boundary_obj, flag=None):
+        if self._two_fields_kernel:
+            return super(FlexibleLBMBoundaryHandling, self)._add_boundary(boundary_obj, flag)
+        else:
+            return self._add_flexible_boundary(boundary_obj, flag)
+
+    def _add_flexible_boundary(self, boundary_obj, flag=None):
+        if boundary_obj not in self._boundary_object_to_boundary_info:
+            sym_index_field = Field.create_generic('indexField', spatial_dimensions=1,
+                                                   dtype=numpy_data_type_for_boundary_object(boundary_obj, self.dim))
+            kernels = [self._create_boundary_kernel(
+                self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.EVEN).compile(),
+                self._create_boundary_kernel(
+                self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.ODD).compile()]
+            if flag is None:
+                flag = self.flag_interface.reserve_next_flag()
+            boundary_info = self.FlexibleBoundaryInfo(self, boundary_obj, flag, kernels)
+            self._boundary_object_to_boundary_info[boundary_obj] = boundary_info
+        return self._boundary_object_to_boundary_info[boundary_obj].flag
+
+    def _create_boundary_kernel(self, symbolic_field, symbolic_index_field, boundary_obj, prev_timestep=Timestep.BOTH):
+        return create_advanced_streaming_boundary_kernel(
+            symbolic_field, symbolic_index_field, self._lb_method, boundary_obj,
+            prev_timestep=prev_timestep, streaming_pattern=self._streaming_pattern,
+            target=self._target, openmp=self._openmp)
+
+    class FlexibleBoundaryInfo(object):
+
+        @property
+        def kernel(self):
+            prev_timestep = self._boundary_handling.prev_timestep
+            if prev_timestep is None:
+                raise Exception(
+                    "The flexible boundary kernel property was accessed while "
+                    + "there was no boundary handling in progress.")
+            return self._kernels[prev_timestep]
+
+        def __init__(self, boundary_handling, boundary_obj, flag, kernels):
+            self._boundary_handling = boundary_handling
+            self.boundary_object = boundary_obj
+            self.flag = flag
+            self._kernels = kernels
+
+# end class FlexibleLBMBoundaryHandling
+
+
+def create_advanced_streaming_boundary_kernel(pdf_field, index_field, lb_method, boundary_functor,
+                                              prev_timestep=Timestep.BOTH, streaming_pattern='pull',
+                                              target='cpu', openmp=True):
+    index_dtype = index_field.dtype.numpy_dtype.fields['dir'][0]
+    offsets_dtype = index_field.dtype.numpy_dtype.fields['x'][0]
+    indexing = AdvancedStreamingBoundaryIndexing(
+        pdf_field, lb_method.stencil, prev_timestep, streaming_pattern, index_dtype, offsets_dtype)
+
+    f_out, f_in = indexing.proxy_fields
+    dir_symbol = indexing.dir_symbol
+    inv_dir = indexing.inverse_dir_symbol
+
+    boundary_assignments = boundary_functor(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field)
+    boundary_assignments = indexing.substitute_proxies(boundary_assignments)
+
+    #   Code Elements inside the loop
+    elements = [Assignment(dir_symbol, index_field[0]('dir'))]
+    elements += boundary_assignments.all_assignments
+
+    kernel = create_indexed_kernel(elements, [index_field], target=target, cpu_openmp=openmp)
+
+    #   Code Elements ahead of the loop
+    index_arrs_node = indexing.create_code_node()
+    for node in boundary_functor.additional_code_nodes[::-1]:
+        kernel.body.insert_front(node)
+    kernel.body.insert_front(index_arrs_node)
+    return kernel

lbmpy/advanced_streaming/boundaryindexing.py

+import numpy as np
+import sympy as sp
+import pystencils as ps
+
+from pystencils.data_types import TypedSymbol, create_type
+from pystencils.backends.cbackend import CustomCodeNode
+
+from lbmpy.advanced_streaming.utility import get_accessor, inverse_dir_index, is_inplace, Timestep
+
+from itertools import product
+
+
+class AdvancedStreamingBoundaryIndexing:
+
+    #   =======================================
+    #       Symbols for usage in boundaries
+    #   =======================================
+
+    @property
+    def proxy_fields(self):
+        return ps.fields(f"f_out({self._q}), f_in({self._q}): [{self._dim}D]")
+
+    @property
+    def dir_symbol(self):
+        return TypedSymbol('dir', create_type(self._index_dtype))
+
+    @property
+    def inverse_dir_symbol(self):
+        return sp.IndexedBase('invdir')
+
+    #   =============================
+    #       Constructor and State
+    #   =============================
+
+    def __init__(self, pdf_field, stencil, prev_timestep=Timestep.BOTH, streaming_pattern='pull',
+                 index_dtype=np.int32, offsets_dtype=np.int32):
+        if prev_timestep == Timestep.BOTH and is_inplace(streaming_pattern):
+            raise ValueError('Cannot create index arrays for both kinds of timesteps for inplace streaming pattern '
+                             + streaming_pattern)
+
+        prev_accessor = get_accessor(streaming_pattern, prev_timestep)
+        next_accessor = get_accessor(streaming_pattern, prev_timestep.next())
+
+        outward_accesses = prev_accessor.write(pdf_field, stencil)
+        inward_accesses = next_accessor.read(pdf_field, stencil)
+
+        self._accesses = {'out': outward_accesses, 'in': inward_accesses}
+
+        self._pdf_field = pdf_field
+        self._stencil = stencil
+        self._dim = len(stencil[0])
+        self._q = len(stencil)
+        self._coordinate_names = ['x', 'y', 'z'][:self._dim]
+
+        self._index_dtype = create_type(index_dtype)
+        self._offsets_dtype = create_type(offsets_dtype)
+
+        self._required_arrays = set()
+        self._trivial_translations = self._collect_trivial_translations()
+
+    def _index_array_symbol(self, f_dir, inverse):
+        assert f_dir in ['in', 'out']
+        inv = '_inv' if inverse else ''
+        name = f"f_{f_dir}{inv}_dir_idx"
+        return TypedSymbol(name, self._index_dtype)
+
+    def _offset_array_symbols(self, f_dir, inverse):
+        assert f_dir in ['in', 'out']
+        inv = '_inv' if inverse else ''
+        name_base = f"f_{f_dir}{inv}_offsets_"
+        symbols = [TypedSymbol(name_base + d, self._index_dtype) for d in self._coordinate_names]
+        return symbols
+
+    def _array_symbols(self, f_dir, inverse, index):
+        if (f_dir, inverse) in self._trivial_translations:
+            return {'index': index, 'offsets': (0, ) * self._dim}
+
+        index_array_symbol = self._index_array_symbol(f_dir, inverse)
+        offset_array_symbols = self._offset_array_symbols(f_dir, inverse)
+        self._required_arrays.add((f_dir, inverse))
+        return {
+            'index': sp.IndexedBase(index_array_symbol, shape=(1,))[index],
+            'offsets': tuple(sp.IndexedBase(s, shape=(1,))[index]
+                             for s in offset_array_symbols)
+        }
+
+    #   =================================
+    #       Proxy fields substitution
+    #   =================================
+
+    def substitute_proxies(self, assignments):
+        if isinstance(assignments, ps.Assignment):
+            assignments = [assignments]
+
+        if not isinstance(assignments, ps.AssignmentCollection):
+            assignments = ps.AssignmentCollection(assignments)
+
+        accesses = self._accesses
+        f_out, f_in = self.proxy_fields
+        inv_dir = self.inverse_dir_symbol
+
+        accessor_subs = dict()
+
+        for fa in assignments.atoms(ps.Field.Access):
+            if fa.field == f_out:
+                f_dir = 'out'
+            elif fa.field == f_in:
+                f_dir = 'in'
+            else:
+                continue
+
+            inv = False
+            idx = fa.index[0]
+            if isinstance(idx, sp.Indexed) and idx.base == inv_dir:
+                idx = idx.indices[0]
+                if isinstance(sp.sympify(idx), sp.Integer):
+                    idx = inverse_dir_index(self._stencil, idx)
+                inv = True
+
+            if isinstance(sp.sympify(idx), sp.Integer):
+                accessor_subs[fa] = accesses[f_dir][idx].get_shifted(*(fa.offsets))
+            else:
+                arr = self._array_symbols(f_dir, inv, idx)
+                accessor_subs[fa] = self._pdf_field[arr['offsets']](arr['index']).get_shifted(*(fa.offsets))
+
+        return assignments.new_with_substitutions(accessor_subs)
+
+    #   =================
+    #       Internals
+    #   =================
+
+    def _get_translated_indices_and_offsets(self, f_dir, inv):
+        accesses = self._accesses[f_dir]
+
+        if inv:
+            inverse_indices = [inverse_dir_index(self._stencil, i)
+                               for i in range(len(self._stencil))]
+            accesses = [accesses[idx] for idx in inverse_indices]
+
+        indices = [a.index[0] for a in accesses]
+        offsets = []
+        for d in range(self._dim):
+            offsets.append([a.offsets[d] for a in accesses])
+        return indices, offsets
+
+    def _collect_trivial_translations(self):
+        trivials = set()
+        trivial_indices = list(range(self._q))
+        trivial_offsets = [[0] * self._q] * self._dim
+        for f_dir, inv in product(['in', 'out'], [False, True]):
+            indices, offsets = self._get_translated_indices_and_offsets(f_dir, inv)
+            if indices == trivial_indices and offsets == trivial_offsets:
+                trivials.add((f_dir, inv))
+        return trivials
+
+    def create_code_node(self):
+        return AdvancedStreamingBoundaryIndexing.TranslationArraysNode(self)
+
+    class TranslationArraysNode(CustomCodeNode):
+
+        def __init__(self, indexing):
+            code = ''
+            symbols_defined = set()
+
+            for f_dir, inv in indexing._required_arrays:
+                indices, offsets = indexing._get_translated_indices_and_offsets(f_dir, inv)
+
+                index_array_symbol = indexing._index_array_symbol(f_dir, inv)
+                symbols_defined.add(index_array_symbol)
+                acc_indices = ", ".join([str(i) for i in indices])
+                code += self._array_pattern(indexing._index_dtype, index_array_symbol.name, acc_indices)
+
+                offset_array_symbols = indexing._offset_array_symbols(f_dir, inv)
+                symbols_defined |= set(offset_array_symbols)
+                for d, arrsymb in enumerate(offset_array_symbols):
+                    acc_offsets = ", ".join([str(o) for o in offsets[d]])
+                    code += self._array_pattern(indexing._offsets_dtype, arrsymb.name, acc_offsets)
+
+            super(AdvancedStreamingBoundaryIndexing.TranslationArraysNode, self).__init__(
+                code, symbols_read=set(), symbols_defined=symbols_defined)
+
+        def _array_pattern(self, dtype, name, content):
+            return f"const {str(dtype)} {name} [] = {{ {content} }}; \n"
+
+        def __str__(self):
+            return "Boundary Access Translation Arrays"
+
+        def __repr__(self):
+            return "Boundary Access Translation Arrays"
+
+#   end class AdvancedStreamingIndexing

lbmpy/advanced_streaming/communication.py

+from pystencils import Field, Assignment
+from pystencils.slicing import shift_slice, get_slice_before_ghost_layer, normalize_slice
+from lbmpy.advanced_streaming.utility import is_inplace, get_accessor, numeric_index, \
+    numeric_offsets, Timestep, get_timesteps
+from pystencils.datahandling import SerialDataHandling
+from itertools import chain
+
+
+def trim_slice_in_direction(slices, direction):
+    assert len(slices) == len(direction)
+
+    result = []
+    for s, d in zip(slices, direction):
+        if isinstance(s, int):
+            result.append(s)
+            continue
+        start = s.start + 1 if d == -1 else s.start
+        stop = s.stop - 1 if d == 1 else s.stop
+        result.append(slice(start, stop, s.step))
+
+    return tuple(result)
+
+
+def extend_dir(direction):
+    if len(direction) == 0:
+        yield tuple()
+    elif direction[0] == 0:
+        for d in [-1, 0, 1]:
+            for rest in extend_dir(direction[1:]):
+                yield (d, ) + rest
+    else:
+        for rest in extend_dir(direction[1:]):
+            yield (direction[0], ) + rest
+
+
+def _get_neighbour_transform(direction, ghost_layers):
+    return tuple(d * (ghost_layers + 1) for d in direction)
+
+
+def _fix_length_one_slices(slices):
+    """Slices of length one are replaced by their start value for correct periodic shifting"""
+    if isinstance(slices, int):
+        return slices
+    elif isinstance(slices, slice):
+        if slices.stop is not None and abs(slices.start - slices.stop) == 1:
+            return slices.start
+        elif slices.stop is None and slices.start == -1:
+            return -1  # [-1:] also has length one
+        else:
+            return slices
+    else:
+        return tuple(_fix_length_one_slices(s) for s in slices)
+
+
+def get_communication_slices(
+        stencil, comm_stencil=None, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1):
+    """
+    Return the source and destination slices for periodicity handling or communication between blocks.
+
+    :param stencil: The stencil used by the LB method.
+    :param comm_stencil: The stencil defining the communication directions. If None, it will be set to stencil.
+    :param streaming_pattern: The streaming pattern.
+    :param prev_timestep: Timestep after which communication is run.
+    :param ghost_layers: Number of ghost layers in each direction.
+
+    """
+    if comm_stencil is None:
+        comm_stencil = stencil
+
+    pdfs = Field.create_generic('pdfs', spatial_dimensions=len(stencil[0]), index_shape=(len(stencil),))
+    write_accesses = get_accessor(streaming_pattern, prev_timestep).write(pdfs, stencil)
+    slices_per_comm_direction = dict()
+
+    for comm_dir in comm_stencil:
+        if all(d == 0 for d in comm_dir):
+            continue
+
+        slices_for_dir = []
+
+        for streaming_dir in set(extend_dir(comm_dir)) & set(stencil):
+            d = stencil.index(streaming_dir)
+            write_offsets = numeric_offsets(write_accesses[d])
+            write_index = numeric_index(write_accesses[d])[0]
+
+            tangential_dir = tuple(s - c for s, c in zip(streaming_dir, comm_dir))
+            origin_slice = get_slice_before_ghost_layer(comm_dir, ghost_layers=ghost_layers, thickness=1)
+            origin_slice = _fix_length_one_slices(origin_slice)
+            src_slice = shift_slice(trim_slice_in_direction(origin_slice, tangential_dir), write_offsets)
+
+            neighbour_transform = _get_neighbour_transform(comm_dir, ghost_layers)
+            dst_slice = shift_slice(src_slice, neighbour_transform)
+
+            src_slice = src_slice + (write_index, )
+            dst_slice = dst_slice + (write_index, )
+
+            slices_for_dir.append((src_slice, dst_slice))
+
+        slices_per_comm_direction[comm_dir] = slices_for_dir
+    return slices_per_comm_direction
+
+
+def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
+                             domain_size=None, target='gpu',
+                             opencl_queue=None, opencl_ctx=None):
+    """Copies a rectangular array slice onto another non-overlapping array slice"""
+    from pystencils.gpucuda.kernelcreation import create_cuda_kernel
+
+    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]
+
+    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):
+        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_cuda_kernel([copy_eq], iteration_slice=dst_slice, skip_independence_check=True)
+    if target == 'gpu':
+        from pystencils.gpucuda import make_python_function
+        return make_python_function(ast)
+    elif target == 'opencl':
+        from pystencils.opencl import make_python_function
+        return make_python_function(ast, opencl_queue, opencl_ctx)
+    else:
+        raise ValueError('Invalid target:', target)
+
+
+class PeriodicityHandling:
+
+    def __init__(self, stencil, data_handling, pdf_field_name,
+                 streaming_pattern='pull', ghost_layers=1, target='cpu',
+                 opencl_queue=None, opencl_ctx=None,
+                 pycuda_direct_copy=True):
+        """
+            Periodicity Handling for Lattice Boltzmann Streaming.
+
+            **On the usage with cuda/opencl:** 
+            - pycuda allows the copying of sliced arrays within device memory using the numpy syntax,
+            e.g. `dst[:,0] = src[:,-1]`. In this implementation, this is the default for periodicity
+            handling. Alternatively, if you set `pycuda_direct_copy=False`, GPU kernels are generated and
+            compiled. The compiled kernels are almost twice as fast in execution as pycuda array copying,
+            but especially for large stencils like D3Q27, their compilation can take up to 20 seconds. 
+            Choose your weapon depending on your use case.
+
+            - pyopencl does not support copying of non-contiguous sliced arrays, so the usage of compiled
+            copy kernels is forced on us. On the positive side, compilation of the OpenCL kernels appears
+            to be about four times faster.
+        """
+        if not isinstance(data_handling, SerialDataHandling):
+            raise ValueError('Only single node data handling is supported!')
+
+        assert target in ['cpu', 'gpu', 'opencl']
+
+        self.stencil = stencil
+        self.dh = data_handling
+        self.pdf_field_name = pdf_field_name
+        self.ghost_layers = ghost_layers
+        periodicity = data_handling.periodicity
+        self.inplace_pattern = is_inplace(streaming_pattern)
+        self.target = target
+        self.cpu = target == 'cpu'
+        self.opencl_queue = opencl_queue
+        self.opencl_ctx = opencl_ctx
+        self.pycuda_direct_copy = target == 'gpu' and pycuda_direct_copy
+
+        def is_copy_direction(direction):
+            for d, p in zip(direction, periodicity):
+                if d != 0 and not p:
+                    return False
+
+            return True
+
+        copy_directions = tuple(filter(is_copy_direction, stencil[1:]))
+        self.comm_slices = []
+        timesteps = get_timesteps(streaming_pattern)
+        for timestep in timesteps:
+            slices_per_comm_dir = get_communication_slices(stencil=stencil,
+                                                           comm_stencil=copy_directions,
+                                                           streaming_pattern=streaming_pattern,
+                                                           prev_timestep=timestep,
+                                                           ghost_layers=ghost_layers)
+            self.comm_slices.append(list(chain.from_iterable(v for k, v in slices_per_comm_dir.items())))
+
+        if target == 'opencl' or (target == 'gpu' and not pycuda_direct_copy):
+            self.device_copy_kernels = []
+            for timestep in timesteps:
+                self.device_copy_kernels.append(self._compile_copy_kernels(timestep))
+
+    def __call__(self, prev_timestep=Timestep.BOTH):
+        if self.cpu:
+            self._periodicity_handling_cpu(prev_timestep)
+        else:
+            self._periodicity_handling_gpu(prev_timestep)
+
+    def _periodicity_handling_cpu(self, prev_timestep):
+        arr = self.dh.cpu_arrays[self.pdf_field_name]
+        comm_slices = self.comm_slices[prev_timestep.idx]
+        for src, dst in comm_slices:
+            arr[dst] = arr[src]
+
+    def _compile_copy_kernels(self, timestep):
+        pdf_field = self.dh.fields[self.pdf_field_name]
+        kernels = []
+        for src, dst in self.comm_slices[timestep.idx]:
+            kernels.append(
+                periodic_pdf_copy_kernel(
+                    pdf_field, src, dst, target=self.target, 
+                    opencl_queue=self.opencl_queue, opencl_ctx=self.opencl_ctx))
+        return kernels
+
+    def _periodicity_handling_gpu(self, prev_timestep):
+        arr = self.dh.gpu_arrays[self.pdf_field_name]
+        if self.pycuda_direct_copy:
+            for src, dst in self.comm_slices[prev_timestep.idx]:
+                arr[dst] = arr[src]
+        else:
+            kernel_args = {self.pdf_field_name: arr}
+            for kernel in self.device_copy_kernels[prev_timestep.idx]:
+                kernel(**kernel_args)

lbmpy/advanced_streaming/utility.py

+from lbmpy.fieldaccess import PdfFieldAccessor, \
+    StreamPullTwoFieldsAccessor, \
+    StreamPushTwoFieldsAccessor, \
+    AAEvenTimeStepAccessor, \
+    AAOddTimeStepAccessor, \
+    EsoTwistEvenTimeStepAccessor, \
+    EsoTwistOddTimeStepAccessor
+
+import pystencils as ps
+from enum import IntEnum
+
+
+class Timestep(IntEnum):
+    EVEN = 0
+    ODD = 1
+    BOTH = 2
+
+    def next(self):
+        return self if self == Timestep.BOTH else Timestep((self + 1) % 2)
+
+    @property
+    def idx(self):
+        """To use this timestep as an array index"""
+        return self % 2
+
+
+streaming_patterns = ['push', 'pull', 'aa', 'esotwist']
+
+even_accessors = {
+    'pull': StreamPullTwoFieldsAccessor,
+    'push': StreamPushTwoFieldsAccessor,
+    'aa': AAEvenTimeStepAccessor,
+    'esotwist': EsoTwistEvenTimeStepAccessor
+}
+
+odd_accessors = {
+    'pull': StreamPullTwoFieldsAccessor,
+    'push': StreamPushTwoFieldsAccessor,
+    'aa': AAOddTimeStepAccessor,
+    'esotwist': EsoTwistOddTimeStepAccessor
+}
+
+
+def get_accessor(streaming_pattern: str, timestep: Timestep) -> PdfFieldAccessor:
+    if streaming_pattern not in streaming_patterns:
+        raise ValueError(
+            "Invalid value of parameter 'streaming_pattern'.", streaming_pattern)
+
+    if timestep == Timestep.EVEN:
+        return even_accessors[streaming_pattern]
+    else:
+        return odd_accessors[streaming_pattern]
+
+
+def is_inplace(streaming_pattern):
+    if streaming_pattern not in streaming_patterns:
+        raise ValueError('Invalid streaming pattern', streaming_pattern)
+
+    return streaming_pattern in ['aa', 'esotwist']
+
+
+def get_timesteps(streaming_pattern):
+    return (Timestep.EVEN, Timestep.ODD) if is_inplace(streaming_pattern) else (Timestep.BOTH, )
+
+
+def numeric_offsets(field_access: ps.Field.Access):
+    return tuple(int(o) for o in field_access.offsets)
+
+
+def numeric_index(field_access: ps.Field.Access):
+    return tuple(int(i) for i in field_access.index)
+
+
+def inverse_dir_index(stencil, dir):
+    return stencil.index(tuple(-d for d in stencil[dir]))
+
+
+class AccessPdfValues:
+    """Allows to access values from a PDF array correctly depending on 
+    the streaming pattern."""
+
+    def __init__(self, pdf_field, stencil, streaming_pattern='pull', timestep=Timestep.BOTH, accessor=None):
+        if accessor is None:
+            accessor = get_accessor(streaming_pattern, timestep)
+        self.read_accs = accessor.read(pdf_field, stencil)
+        self.write_accs = accessor.write(pdf_field, stencil)
+
+    def write_pdf(self, pdf_arr, pos, d, value):
+        offsets = self.write_accs[d].offsets
+        pos = tuple(p + o for p, o in zip(pos, offsets))
+        i = self.write_accs[d].index[0]
+        pdf_arr[pos + (i,)] = value
+
+    def read_pdf(self, pdf_arr, pos, d):
+        offsets = self.read_accs[d].offsets
+        pos = tuple(p + o for p, o in zip(pos, offsets))
+        i = self.read_accs[d].index[0]
+        return pdf_arr[pos + (i,)]

lbmpy/macroscopic_value_kernels.py

@@ -4,6 +4,8 @@ from copy import deepcopy
 from lbmpy.simplificationfactory import create_simplification_strategy
 from pystencils.field import Field, get_layout_of_array
 
+from lbmpy.advanced_streaming.utility import get_accessor, Timestep
+
 
 def pdf_initialization_assignments(lb_method, density, velocity, pdfs):
     """Assignments to initialize the pdf field with equilibrium"""
@@ -28,8 +30,27 @@ def macroscopic_values_getter(lb_method, density, velocity, pdfs):
 
 macroscopic_values_setter = pdf_initialization_assignments
 
+#   Extensions for any streaming patterns
+
+
+def flexible_macroscopic_values_setter(lb_method, density, velocity,
+                                       pdf_field, streaming_pattern='pull', previous_timestep=Timestep.BOTH):
+    previous_step_accessor = get_accessor(streaming_pattern, previous_timestep)
+    write_accesses = previous_step_accessor.write(pdf_field, lb_method.stencil)
+    return pdf_initialization_assignments(lb_method, density, velocity, write_accesses)
+
 
-def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None, field_layout='numpy', target='cpu'):
+def flexible_macroscopic_values_getter(lb_method, density, velocity,
+                                       pdf_field, streaming_pattern='pull', previous_timestep=Timestep.BOTH):
+    previous_step_accessor = get_accessor(streaming_pattern, previous_timestep)
+    write_accesses = previous_step_accessor.write(pdf_field, lb_method.stencil)
+    return macroscopic_values_getter(lb_method, density, velocity, write_accesses)
+
+
+def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None,
+                                      ghost_layers=1, iteration_slice=None,
+                                      field_layout='numpy', target='cpu',
+                                      streaming_pattern='pull', previous_timestep=Timestep.BOTH):
     """
     Create kernel to compute macroscopic value(s) from a pdf field (e.g. density or velocity)
 
@@ -37,8 +58,13 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
         lb_method: instance of :class:`lbmpy.methods.AbstractLbMethod`
         output_quantities: sequence of quantities to compute e.g. ['density', 'velocity']
         pdf_arr: optional numpy array for pdf field - used to get optimal loop structure for kernel
+        ghost_layers: a sequence of pairs for each coordinate with lower and upper nr of ghost layers
+                      that should be excluded from the iteration. If None, the number of ghost layers 
+                      is determined automatically and assumed to be equal for all dimensions.        
+        iteration_slice: if not None, iteration is done only over this slice of the field
         field_layout: layout for output field, also used for pdf field if pdf_arr is not given
         target: 'cpu' or 'gpu'
+        previous_step_accessor: The accessor used by the streaming pattern of the previous timestep
 
     Returns:
         a function to compute macroscopic values:
@@ -83,15 +109,19 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
             output_mapping[output_quantity] = output_mapping[output_quantity][0]
 
     stencil = lb_method.stencil
-    pdf_symbols = [pdf_field(i) for i in range(len(stencil))]
+    previous_step_accessor = get_accessor(streaming_pattern, previous_timestep)
+    pdf_symbols = previous_step_accessor.write(pdf_field, stencil)
+
     eqs = cqc.output_equations_from_pdfs(pdf_symbols, output_mapping).all_assignments
 
     if target == 'cpu':
         import pystencils.cpu as cpu
-        kernel = cpu.make_python_function(cpu.create_kernel(eqs))
+        kernel = cpu.make_python_function(cpu.create_kernel(
+            eqs, ghost_layers=ghost_layers, iteration_slice=iteration_slice))
     elif target == 'gpu':
         import pystencils.gpucuda as gpu
-        kernel = gpu.make_python_function(gpu.create_cuda_kernel(eqs))
+        kernel = gpu.make_python_function(gpu.create_cuda_kernel(
+            eqs, ghost_layers=ghost_layers, iteration_slice=iteration_slice))
     else:
         raise ValueError("Unknown target '%s'. Possible targets are 'cpu' and 'gpu'" % (target,))
 
@@ -107,7 +137,10 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
     return getter
 
 
-def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None, field_layout='numpy', target='cpu'):
+def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None, 
+                                      ghost_layers=1, iteration_slice=None,
+                                      field_layout='numpy', target='cpu', 
+                                      streaming_pattern='pull', previous_timestep=Timestep.BOTH):
     """
     Creates a function that sets a pdf field to specified macroscopic quantities
     The returned function can be called with the pdf field to set as single argument
@@ -116,8 +149,13 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
         lb_method: instance of :class:`lbmpy.methods.AbstractLbMethod`
         quantities_to_set: map from conserved quantity name to fixed value or numpy array
         pdf_arr: optional numpy array for pdf field - used to get optimal loop structure for kernel
+        ghost_layers: a sequence of pairs for each coordinate with lower and upper nr of ghost layers
+                      that should be excluded from the iteration. If None, the number of ghost layers 
+                      is determined automatically and assumed to be equal for all dimensions.        
+        iteration_slice: if not None, iteration is done only over this slice of the field
         field_layout: layout of the pdf field if pdf_arr was not given
         target: 'cpu' or 'gpu'
+        previous_step_accessor: The accessor used by the streaming pattern of the previous timestep
 
     Returns:
         function taking pdf array as single argument and which sets the field to the given values
@@ -155,7 +193,10 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
     else:
         eq = eq.new_without_subexpressions()
 
-    substitutions = {sym: pdf_field(i) for i, sym in enumerate(lb_method.post_collision_pdf_symbols)}
+    previous_step_accessor = get_accessor(streaming_pattern, previous_timestep)
+    write_accesses = previous_step_accessor.write(pdf_field, lb_method.stencil)
+
+    substitutions = {sym: write_accesses[i] for i, sym in enumerate(lb_method.post_collision_pdf_symbols)}
     eq = eq.new_with_substitutions(substitutions).all_assignments
 
     if target == 'cpu':

lbmpy_tests/advanced_streaming/test_advanced_streaming_noslip.py

+import numpy as np
+import sympy as sp
+
+from lbmpy.advanced_streaming import FlexibleNoSlip, create_advanced_streaming_boundary_kernel, Timestep
+from lbmpy.advanced_streaming.utility import even_accessors, odd_accessors, streaming_patterns, inverse_dir_index, AccessPdfValues
+from lbmpy.creationfunctions import create_lb_method
+from lbmpy.stencils import get_stencil
+
+import pystencils as ps
+
+from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
+from pystencils.data_types import TypedSymbol, create_type
+from pystencils.field import Field, FieldType
+
+from itertools import product
+
+import pytest
+
+@pytest.mark.parametrize("stencil", [ 'D2Q9', 'D3Q19', 'D3Q27'])
+@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
+@pytest.mark.parametrize("prev_timestep", [Timestep.EVEN, Timestep.ODD])
+def test_advanced_streaming_noslip_single_cell(stencil, streaming_pattern, prev_timestep):
+    """
+    Flexible NoSlip Test
+    """
+
+    stencil = get_stencil(stencil)
+    q = len(stencil)
+    dim = len(stencil[0])
+    pdf_field = ps.fields(f'pdfs({q}): [{dim}D]')
+
+    prev_pdf_access = AccessPdfValues(pdf_field, stencil, streaming_pattern, prev_timestep)
+    next_pdf_access = AccessPdfValues(pdf_field, stencil, streaming_pattern, prev_timestep.next())
+
+    pdfs = np.zeros((3,) * dim + (q,))
+    pos = (1,) * dim
+    for d in range(q):
+        prev_pdf_access.write_pdf(pdfs, pos, d, d)
+
+    lb_method = create_lb_method(stencil=stencil, method='srt')
+    noslip = FlexibleNoSlip()
+
+    index_struct_dtype = numpy_data_type_for_boundary_object(noslip, dim)
+
+    index_field = Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0],
+                    shape=(TypedSymbol("indexVectorSize", create_type(np.int64)), 1), strides=(1, 1))
+    index_vector = np.array([ pos + (d,) for d in range(q) ], dtype=index_struct_dtype)
+
+    flex_ast = create_advanced_streaming_boundary_kernel(pdf_field, 
+                    index_field, lb_method, noslip,
+                    prev_timestep=prev_timestep,
+                    streaming_pattern=streaming_pattern)
+
+    flex_kernel = flex_ast.compile()
+    
+    flex_kernel(pdfs=pdfs, indexVector=index_vector, indexVectorSize=len(index_vector))
+
+    reflected_pdfs = [ next_pdf_access.read_pdf(pdfs, pos, d) for d in range(q)]
+    inverse_pdfs = [ inverse_dir_index(stencil, d) for d in range(q) ] 
+    assert reflected_pdfs == inverse_pdfs
\ No newline at end of file

lbmpy_tests/advanced_streaming/test_communication.py

+import numpy as np
+from lbmpy.stencils import get_stencil
+from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice
+from lbmpy.advanced_streaming.communication import get_communication_slices, _fix_length_one_slices
+from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
+
+import pytest
+
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
+def test_slices_not_empty(stencil, streaming_pattern, timestep):
+    stencil = get_stencil(stencil)
+    dim = len(stencil[0])
+    q = len(stencil)
+    arr = np.zeros( (4,) * dim + (q,) )
+    slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep, ghost_layers=1)
+    for _, slices_list in slices.items():
+        for src, dst in slices_list:
+            assert all(s != 0 for s in arr[src].shape)
+            assert all(s != 0 for s in arr[dst].shape)
+
+
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
+def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
+    stencil = get_stencil(stencil)
+    dim = len(stencil[0])
+    q = len(stencil)
+    arr = np.zeros( (4,) * dim + (q,) )
+    slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep, ghost_layers=1)
+    for _, slices_list in slices.items():
+        for src, dst in slices_list:
+            src_shape = arr[src].shape
+            dst_shape = arr[dst].shape
+            assert src_shape == dst_shape
+
+
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
+def test_pull_communication_slices(stencil):
+    stencil = get_stencil(stencil)
+
+    slices = get_communication_slices(
+        stencil, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1)
+    for i, d in enumerate(stencil):
+        if i == 0:
+            continue
+
+        for s in slices[d]:
+            if s[0][-1] == i:
+                src = s[0][:-1]
+                dst = s[1][:-1]
+                break
+        
+        inner_slice = _fix_length_one_slices(get_slice_before_ghost_layer(d, ghost_layers=1))
+        inv_dir = (-e for e in d)
+        gl_slice = _fix_length_one_slices(get_ghost_region_slice(inv_dir, ghost_layers=1))
+        assert src == inner_slice
+        assert dst == gl_slice
\ No newline at end of file

lbmpy_tests/advanced_streaming/test_fully_periodic_flow.py

+import numpy as np
+import sympy as sp
+
+from pystencils.datahandling import create_data_handling
+from pystencils import create_kernel
+from pystencils.plot import scalar_field, vector_field, vector_field_magnitude
+
+from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_function
+from lbmpy.macroscopic_value_kernels import flexible_macroscopic_values_getter, flexible_macroscopic_values_setter
+from lbmpy.stencils import get_stencil
+
+from lbmpy.advanced_streaming import PeriodicityHandling
+from lbmpy.advanced_streaming.utility import is_inplace, streaming_patterns, get_timesteps
+
+import pytest
+from numpy.testing import assert_allclose, assert_array_equal
+
+all_results = dict()
+
+targets = ['cpu']
+
+try:
+    import pycuda.autoinit
+    targets += ['gpu']
+except ImportError:
+    pass
+
+try:
+    import pystencils.opencl.autoinit
+    targets += ['opencl']
+except ImportError:
+    pass
+
+
+@pytest.mark.parametrize('target', targets)
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+def test_fully_periodic_flow(target, stencil, streaming_pattern):
+
+    if target == 'opencl':
+        from pystencils.opencl.opencljit import get_global_cl_queue
+        opencl_queue = get_global_cl_queue()
+    else:
+        opencl_queue = None
+
+    gpu = target in ['gpu', 'opencl']
+
+    #   Stencil
+    stencil = get_stencil(stencil)
+    q = len(stencil)
+    dim = len(stencil[0])
+
+    #   Streaming
+    inplace = is_inplace(streaming_pattern)
+    timesteps = get_timesteps(streaming_pattern)
+    zeroth_timestep = timesteps[0]
+
+    #   Data Handling and PDF fields
+    domain_size = (30,) * dim
+    periodicity = (True,) * dim
+
+    dh = create_data_handling(domain_size=domain_size, periodicity=periodicity,
+                              default_target=target, opencl_queue=opencl_queue)
+
+    pdfs = dh.add_array('pdfs', q)
+    if not inplace:
+        pdfs_tmp = dh.add_array_like('pdfs_tmp', pdfs.name)
+
+    #   LBM Streaming and Collision
+    method_params = {
+        'stencil': stencil,
+        'method': 'srt',
+        'relaxation_rate': 1.0
+    }
+
+    optimization = {
+        'symbolic_field': pdfs,
+        'target': target
+    }
+
+    if not inplace:
+        optimization['symbolic_temporary_field'] = pdfs_tmp
+
+    lb_collision = create_lb_collision_rule(optimization=optimization, **method_params)
+    lb_method = lb_collision.method
+
+    lb_kernels = []
+    if inplace:
+        lb_kernels.append(create_lb_function(collision_rule=lb_collision,
+                                             optimization=optimization,
+                                             kernel_type=streaming_pattern + '_even',
+                                             **method_params))
+        lb_kernels.append(create_lb_function(collision_rule=lb_collision,
+                                             optimization=optimization,
+                                             kernel_type=streaming_pattern + '_odd',
+                                             **method_params))
+    else:
+        if streaming_pattern == 'pull':
+            kernel_type = 'stream_pull_collide'
+        elif streaming_pattern == 'push':
+            kernel_type = 'collide_stream_push'
+        lb_kernels.append(create_lb_function(collision_rule=lb_collision,
+                                             optimization=optimization,
+                                             kernel_type=kernel_type,
+                                             **method_params))
+
+    #   Macroscopic Values
+    density = 1.0
+    density_field = dh.add_array('rho', 1)
+    u_x = 0.01
+    velocity = (u_x,) * dim
+    velocity_field = dh.add_array('u', dim)
+
+    u_ref = np.full(domain_size + (dim,), u_x)
+
+    setter = flexible_macroscopic_values_setter(
+        lb_method, density, velocity, pdfs,
+        streaming_pattern=streaming_pattern, previous_timestep=zeroth_timestep)
+    setter_kernel = create_kernel(setter, ghost_layers=1, target=target).compile()
+
+    getter_kernels = []
+    for t in timesteps:
+        getter = flexible_macroscopic_values_getter(
+            lb_method, density_field, velocity_field, pdfs,
+            streaming_pattern=streaming_pattern, previous_timestep=t)
+        getter_kernels.append(create_kernel(getter, ghost_layers=1, target=target).compile())
+
+    #   Periodicity
+    periodicity_handler = PeriodicityHandling(stencil, dh, pdfs.name, streaming_pattern=streaming_pattern, target=target)
+
+    # Initialization and Timestep
+    current_timestep = zeroth_timestep
+
+    def init():
+        global current_timestep
+        current_timestep = zeroth_timestep
+        dh.run_kernel(setter_kernel)
+
+    def one_step():
+        global current_timestep
+
+        # Periodicty
+        periodicity_handler(current_timestep)
+
+        # Here, the next time step begins
+        current_timestep = current_timestep.next()
+
+        # LBM Step
+        dh.run_kernel(lb_kernels[current_timestep.idx])
+
+        # Field Swaps
+        if not inplace:
+            dh.swap(pdfs.name, pdfs_tmp.name)
+
+        # Macroscopic Values
+        dh.run_kernel(getter_kernels[current_timestep.idx])
+
+    #   Run the simulation
+    init()
+
+    for _ in range(100):
+        one_step()
+
+    #   Evaluation
+    if gpu:
+        dh.to_cpu(velocity_field.name)
+    u = dh.gather_array(velocity_field.name)
+
+    #   Equal to the steady-state velocity field up to numerical errors
+    assert_allclose(u, u_ref)
+
+    #   Flow must be equal up to numerical error for all streaming patterns
+    global all_results
+    for key, prev_u in all_results.items():
+        if key[0] == stencil:
+            prev_pattern = key[1]
+            assert_allclose(
+                u, prev_u, err_msg=f'Velocity field for {streaming_pattern} differed from {prev_pattern}!')
+    all_results[(stencil, streaming_pattern)] = u

lbmpy_tests/advanced_streaming/test_periodic_pipe_with_force.py

+import numpy as np
+import sympy as sp
+
+from pystencils.datahandling import create_data_handling
+from pystencils import create_kernel
+from pystencils.slicing import make_slice
+
+from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_function
+from lbmpy.macroscopic_value_kernels import flexible_macroscopic_values_getter, flexible_macroscopic_values_setter
+from lbmpy.stencils import get_stencil
+
+from lbmpy.advanced_streaming import PeriodicityHandling, FlexibleNoSlip, FlexibleLBMBoundaryHandling
+from lbmpy.advanced_streaming.utility import is_inplace, streaming_patterns, Timestep, get_timesteps
+
+import pytest
+from numpy.testing import assert_allclose
+
+all_results = dict()
+
+
+class PeriodicPipeFlow:
+    def __init__(self, stencil, streaming_pattern):
+
+        #   Stencil
+        self.stencil = stencil
+        self.q = len(self.stencil)
+        self.dim = len(self.stencil[0])
+
+        #   Streaming
+        self.streaming_pattern = streaming_pattern
+        self.inplace = is_inplace(self.streaming_pattern)
+        self.timesteps = get_timesteps(streaming_pattern)
+        self.zeroth_timestep = self.timesteps[0]
+
+        #   Domain, Data Handling and PDF fields
+        self.pipe_length = 60
+        self.pipe_radius = 15
+        self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
+        self.periodicity = (True, ) + (False, ) * (self.dim - 1)
+        self.force = (0.0001, ) + (0.0,) * (self.dim - 1)
+
+        self.dh = create_data_handling(domain_size=self.domain_size, periodicity=self.periodicity)
+
+        self.pdfs = self.dh.add_array('pdfs', self.q)
+        if not self.inplace:
+            self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
+
+        #   LBM Streaming and Collision
+        method_params = {
+            'stencil': stencil,
+            'method': 'srt',
+            'relaxation_rate': 1.0,
+            'force_model': 'guo',
+            'force': self.force
+        }
+
+        optimization = {
+            'symbolic_field': self.pdfs,
+            'target': 'cpu'
+        }
+
+        if not self.inplace:
+            optimization['symbolic_temporary_field'] = self.pdfs_tmp
+
+        self.lb_collision = create_lb_collision_rule(optimization=optimization, **method_params)
+        self.lb_method = self.lb_collision.method
+
+        self.lb_kernels = []
+        if self.inplace:
+            self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
+                                                      optimization=optimization,
+                                                      kernel_type=self.streaming_pattern + '_even',
+                                                      **method_params))
+            self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
+                                                      optimization=optimization,
+                                                      kernel_type=self.streaming_pattern + '_odd',
+                                                      **method_params))
+        else:
+            if self.streaming_pattern == 'pull':
+                kernel_type = 'stream_pull_collide'
+            elif self.streaming_pattern == 'push':
+                kernel_type = 'collide_stream_push'
+            self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
+                                                      optimization=optimization,
+                                                      kernel_type=kernel_type,
+                                                      **method_params))
+
+        #   Macroscopic Values
+        self.density = 1.0
+        self.density_field = self.dh.add_array('rho', 1)
+        u_x = 0.0
+        self.velocity = (u_x,) * self.dim
+        self.velocity_field = self.dh.add_array('u', self.dim)
+
+        setter = flexible_macroscopic_values_setter(
+            self.lb_method, self.density, self.velocity, self.pdfs,
+            streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
+        self.init_kernel = create_kernel(setter, ghost_layers=1).compile()
+
+        self.getter_kernels = []
+        for t in self.timesteps:
+            getter = flexible_macroscopic_values_getter(
+                self.lb_method, self.density_field, self.velocity_field, self.pdfs,
+                streaming_pattern=self.streaming_pattern, previous_timestep=t)
+            self.getter_kernels.append(create_kernel(getter, ghost_layers=1).compile())
+
+        #   Periodicity
+        self.periodicity_handler = PeriodicityHandling(self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
+
+        #   Boundary Handling
+        self.noslip = FlexibleNoSlip()
+        self.bh = FlexibleLBMBoundaryHandling(self.lb_method, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
+
+        self.bh.set_boundary(boundary_obj=self.noslip, mask_callback=self.mask_callback)
+
+        self.current_timestep = self.zeroth_timestep
+
+    def mask_callback(self, x, y, z=None):
+        y = y - self.pipe_radius
+        z = z - self.pipe_radius if z is not None else 0
+        return np.sqrt(y**2 + z**2) >= self.pipe_radius
+
+    def init(self):
+        self.current_timestep = self.zeroth_timestep
+        self.dh.run_kernel(self.init_kernel)
+
+    def step(self):
+        #   Order matters! First communicate, then boundaries, otherwise
+        #   periodicity handling overwrites reflected populations
+        # Periodicty
+        self.periodicity_handler(self.current_timestep)
+
+        # Boundaries
+        self.bh(prev_timestep=self.current_timestep)
+
+        # Here, the next time step begins
+        self.current_timestep = self.current_timestep.next()
+
+        # LBM Step
+        self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
+
+        # Field Swaps
+        if not self.inplace:
+            self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
+
+        # Macroscopic Values
+        self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
+
+    def run(self, iterations):
+        for _ in range(iterations):
+            self.step()
+
+    @property
+    def velocity_array(self):
+        return self.dh.gather_array(self.velocity_field.name)
+
+    def get_trimmed_velocity_array(self):
+        u = np.copy(self.dh.gather_array(self.velocity_field.name))
+        mask = self.bh.get_mask(None, self.noslip)
+        for idx in np.ndindex(u.shape[:-1]):
+            if mask[idx] != 0:
+                u[idx] = np.full((self.dim, ), np.nan)
+        return u
+
+
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+def test_periodic_pipe(stencil, streaming_pattern):
+    stencil = get_stencil(stencil)
+    pipeflow = PeriodicPipeFlow(stencil, streaming_pattern)
+    pipeflow.init()
+    pipeflow.run(100)
+    u = pipeflow.get_trimmed_velocity_array()
+
+    #   Flow must be equal up to numerical error for all streaming patterns
+    global all_results
+    for key, prev_u in all_results.items():
+        if key[0] == stencil:
+            prev_pattern = key[1]
+            assert_allclose(
+                u, prev_u,
+                rtol=1, atol=1e-16,
+                err_msg=f'Velocity field for {streaming_pattern} differed from {prev_pattern}!')
+    all_results[(stencil, streaming_pattern)] = u

lbmpy_tests/test_macroscopic_value_kernels.py

@@ -3,30 +3,45 @@ import numpy as np
 from lbmpy.creationfunctions import create_lb_method
 from lbmpy.macroscopic_value_kernels import (
     compile_macroscopic_values_getter, compile_macroscopic_values_setter)
+from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
 
+import pytest
 
-def test_set_get_density_velocity_with_fields():
-    for stencil in ['D2Q9', 'D3Q19']:
-        for force_model in ['guo', 'luo', 'none']:
-            for compressible in [True, False]:
-                force = (0.1, 0.12, -0.17)
-                method = create_lb_method(stencil=stencil, force_model=force_model, method='trt',
-                                          compressible=compressible, force=force)
-                size = (3, 7, 4)[:method.dim]
-                pdf_arr = np.zeros(size + (len(method.stencil),))
-                density_input_field = 1 + 0.2 * (np.random.random_sample(size) - 0.5)
-                velocity_input_field = 0.1 * (np.random.random_sample(size + (method.dim, )) - 0.5)
-                setter = compile_macroscopic_values_setter(method, pdf_arr=pdf_arr,
-                                                           quantities_to_set={'density': density_input_field,
-                                                                              'velocity': velocity_input_field}, )
-                setter(pdf_arr)
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19'])
+@pytest.mark.parametrize('force_model', ['guo', 'luo', None])
+@pytest.mark.parametrize('compressible', ['compressible', False])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+@pytest.mark.parametrize('prev_timestep', [Timestep.EVEN, Timestep.ODD])
+def test_set_get_density_velocity_with_fields(stencil, force_model, compressible, streaming_pattern, prev_timestep):
+    force = (0.1, 0.12, -0.17)
+    method = create_lb_method(stencil=stencil, force_model=force_model, method='trt',
+                              compressible=compressible, force=force)
+    ghost_layers = 1
+    inner_size = (3, 7, 4)[:method.dim]
+    total_size = tuple(s + 2 * ghost_layers for s in inner_size)
+    pdf_arr = np.zeros(total_size + (len(method.stencil),))
+
+    inner_slice = (slice(ghost_layers, -ghost_layers, 1), ) * method.dim
+    density_input_field = np.zeros(total_size)
+    velocity_input_field = np.zeros(total_size + (method.dim, ))
+    density_input_field[inner_slice] = 1 + 0.2 * (np.random.random_sample(inner_size) - 0.5)
+    velocity_input_field[inner_slice] = 0.1 * \
+        (np.random.random_sample(inner_size + (method.dim, )) - 0.5)
+
+    quantities_to_set = {'density': density_input_field, 'velocity': velocity_input_field}
+
+    setter = compile_macroscopic_values_setter(method, pdf_arr=pdf_arr, quantities_to_set=quantities_to_set,
+                                               ghost_layers=ghost_layers, streaming_pattern=streaming_pattern, previous_timestep=prev_timestep)
+    setter(pdf_arr)
+
+    getter = compile_macroscopic_values_getter(method, ['density', 'velocity'], pdf_arr=pdf_arr, 
+                                               ghost_layers=ghost_layers, streaming_pattern=streaming_pattern, previous_timestep=prev_timestep)
 
-                getter = compile_macroscopic_values_getter(method, ['density', 'velocity'], pdf_arr=pdf_arr)
-                density_output_field = np.empty_like(density_input_field)
-                velocity_output_field = np.empty_like(velocity_input_field)
-                getter(pdfs=pdf_arr, density=density_output_field, velocity=velocity_output_field)
-                np.testing.assert_almost_equal(density_input_field, density_output_field)
-                np.testing.assert_almost_equal(velocity_input_field, velocity_output_field)
+    density_output_field = np.zeros_like(density_input_field)
+    velocity_output_field = np.zeros_like(velocity_input_field)
+    getter(pdfs=pdf_arr, density=density_output_field, velocity=velocity_output_field)
+    np.testing.assert_almost_equal(density_input_field, density_output_field)
+    np.testing.assert_almost_equal(velocity_input_field, velocity_output_field)
 
 
 def test_set_get_constant_velocity():
@@ -40,11 +55,11 @@ def test_set_get_constant_velocity():
                                           compressible=compressible, force=force)
                 size = (1, 1, 1)[:method.dim]
                 pdf_arr = np.zeros(size + (len(method.stencil),))
-                setter = compile_macroscopic_values_setter(method, pdf_arr=pdf_arr,
+                setter = compile_macroscopic_values_setter(method, pdf_arr=pdf_arr, ghost_layers=0,
                                                            quantities_to_set={'velocity': ref_velocity[:method.dim]}, )
                 setter(pdf_arr)
 
-                getter = compile_macroscopic_values_getter(method, ['velocity'], pdf_arr=pdf_arr)
+                getter = compile_macroscopic_values_getter(method, ['velocity'], pdf_arr=pdf_arr, ghost_layers=0)
                 velocity_output_field = np.zeros(size + (method.dim, ))
                 getter(pdfs=pdf_arr, velocity=velocity_output_field)
                 if method.dim == 2: