doc/notebooks/00_tutorial_lbmpy_walberla_overview.ipynb

@@ -638,7 +638,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.5"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,

doc/notebooks/04_tutorial_cumulant_LBM.ipynb

@@ -673,7 +673,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -687,7 +687,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.2"
+   "version": "3.11.4"
   },
   "vscode": {
    "interpreter": {

doc/notebooks/demo_shallow_water_lbm.ipynb

@@ -541,7 +541,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -555,7 +555,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.2"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,

doc/sphinx/lbmpy.bib

@@ -266,4 +266,14 @@ journal = {Communications in Computational Physics}
   publisher = {Springer Link},
 }
 
+@article{BouzidiBC,
+    author = {Bouzidi, M’hamed and Firdaouss, Mouaouia and Lallemand, Pierre},
+    title = "{Momentum transfer of a Boltzmann-lattice fluid with boundaries}",
+    journal = {Physics of Fluids},
+    year = {2001},
+    month = {11},
+    doi = {10.1063/1.1399290},
+    url = {https://doi.org/10.1063/1.1399290},
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}

doc/sphinx/tutorials.rst

@@ -24,6 +24,7 @@ You can open the notebooks directly to play around with the code examples.
     /notebooks/demo_create_method_from_scratch.ipynb
     /notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
     /notebooks/demo_automatic_chapman_enskog_analysis.ipynb
+    /notebooks/demo_interpolation_boundary_conditions.ipynb
     /notebooks/demo_thermalized_lbm.ipynb
     /notebooks/demo_shallow_water_lbm.ipynb
     /notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb

lbmpy/__init__.py

@@ -12,7 +12,7 @@ from .stencils import LBStencil
 
 
 __all__ = ['create_lb_ast', 'create_lb_collision_rule', 'create_lb_function', 'create_lb_method',
-           'create_lb_method_from_existing', 'create_lb_update_rule', 'LBMConfig', 'LBMOptimisation',
+           'create_lb_update_rule', 'LBMConfig', 'LBMOptimisation',
            'Stencil', 'Method', 'ForceModel', 'CollisionSpace',
            'LatticeBoltzmannStep',
            'pdf_initialization_assignments', 'macroscopic_values_getter', 'compile_macroscopic_values_getter',

lbmpy/advanced_streaming/__init__.py

-from .indexing import BetweenTimestepsIndexing, NeighbourOffsetArrays
+from .indexing import BetweenTimestepsIndexing
 from .communication import get_communication_slices, LBMPeriodicityHandling
 from .utility import Timestep, get_accessor, is_inplace, get_timesteps, \
     numeric_index, numeric_offsets, inverse_dir_index, AccessPdfValues
 
-__all__ = ['BetweenTimestepsIndexing', 'NeighbourOffsetArrays',
+__all__ = ['BetweenTimestepsIndexing',
            'get_communication_slices', 'LBMPeriodicityHandling',
            'Timestep', 'get_accessor', 'is_inplace', 'get_timesteps',
            'numeric_index', 'numeric_offsets', 'inverse_dir_index', 'AccessPdfValues']

lbmpy/advanced_streaming/indexing.py

@@ -3,17 +3,12 @@ import sympy as sp
 import pystencils as ps
 
 from pystencils.typing 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 lbmpy.custom_code_nodes import TranslationArraysNode
 
 from itertools import product
 
 
-def _array_pattern(dtype, name, content):
-    return f"const {str(dtype)} {name} [] = {{ {','.join(str(c) for c in content)} }}; \n"
-
-
 class BetweenTimestepsIndexing:
 
     #   ==============================================
@@ -30,7 +25,7 @@ class BetweenTimestepsIndexing:
 
     @property
     def inverse_dir_symbol(self):
-        """Symbol denoting the inversion of a PDF field index. 
+        """Symbol denoting the inversion of a PDF field index.
         Use only at top-level of index to f_out or f_in, otherwise it can't be correctly replaced."""
         return sp.IndexedBase('invdir')
 
@@ -168,90 +163,21 @@ class BetweenTimestepsIndexing:
         return trivial_index_translations, trivial_offset_translations
 
     def create_code_node(self):
-        return BetweenTimestepsIndexing.TranslationArraysNode(self)
-
-    class TranslationArraysNode(CustomCodeNode):
-
-        def __init__(self, indexing):
-            code = ''
-            symbols_defined = set()
-
-            for f_dir, inv in indexing._required_index_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)
-                code += _array_pattern(indexing._index_dtype, index_array_symbol.name, indices)
-
-            for f_dir, inv in indexing._required_offset_arrays:
-                indices, offsets = indexing._get_translated_indices_and_offsets(f_dir, inv)
-                offset_array_symbols = indexing._offset_array_symbols(f_dir, inv)
-                symbols_defined |= set(offset_array_symbols)
-                for d, arrsymb in enumerate(offset_array_symbols):
-                    code += _array_pattern(indexing._offsets_dtype, arrsymb.name, offsets[d])
-
-            super(BetweenTimestepsIndexing.TranslationArraysNode, self).__init__(
-                code, symbols_read=set(), symbols_defined=symbols_defined)
+        array_content = list()
+        symbols_defined = set()
+        for f_dir, inv in self._required_index_arrays:
+            indices, offsets = self._get_translated_indices_and_offsets(f_dir, inv)
+            index_array_symbol = self._index_array_symbol(f_dir, inv)
+            symbols_defined.add(index_array_symbol)
+            array_content.append((self._index_dtype, index_array_symbol.name, indices))
 
-        def __str__(self):
-            return "Variable PDF Access Translation Arrays"
+        for f_dir, inv in self._required_offset_arrays:
+            indices, offsets = self._get_translated_indices_and_offsets(f_dir, inv)
+            offset_array_symbols = self._offset_array_symbols(f_dir, inv)
+            symbols_defined |= set(offset_array_symbols)
+            for d, arrsymb in enumerate(offset_array_symbols):
+                array_content.append((self._offsets_dtype, arrsymb.name, offsets[d]))
 
-        def __repr__(self):
-            return "Variable PDF Access Translation Arrays"
+        return TranslationArraysNode(array_content, symbols_defined)
 
 #   end class AdvancedStreamingIndexing
-
-
-class NeighbourOffsetArrays(CustomCodeNode):
-
-    @staticmethod
-    def neighbour_offset(dir_idx, stencil):
-        if isinstance(sp.sympify(dir_idx), sp.Integer):
-            return stencil[dir_idx]
-        else:
-            return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
-                         for symbol in NeighbourOffsetArrays._offset_symbols(len(stencil[0]))])
-
-    @staticmethod
-    def _offset_symbols(dim):
-        return [TypedSymbol(f"neighbour_offset_{d}", create_type(np.int32)) for d in ['x', 'y', 'z'][:dim]]
-
-    def __init__(self, stencil, offsets_dtype=np.int32):
-        offsets_dtype = create_type(offsets_dtype)
-        dim = len(stencil[0])
-
-        array_symbols = NeighbourOffsetArrays._offset_symbols(dim)
-        code = "\n"
-        for i, arrsymb in enumerate(array_symbols):
-            code += _array_pattern(offsets_dtype, arrsymb.name, (d[i] for d in stencil))
-
-        offset_symbols = NeighbourOffsetArrays._offset_symbols(dim)
-        super(NeighbourOffsetArrays, self).__init__(code, symbols_read=set(),
-                                                    symbols_defined=set(offset_symbols))
-
-
-class MirroredStencilDirections(CustomCodeNode):
-
-    @staticmethod
-    def mirror_stencil(direction, mirror_axis):
-        assert mirror_axis <= len(direction), f"only {len(direction)} axis available for mirage"
-        direction = list(direction)
-        direction[mirror_axis] = -direction[mirror_axis]
-
-        return tuple(direction)
-
-    @staticmethod
-    def _mirrored_symbol(mirror_axis):
-        axis = ['x', 'y', 'z']
-        return TypedSymbol(f"{axis[mirror_axis]}_axis_mirrored_stencil_dir", create_type(np.int32))
-
-    def __init__(self, stencil, mirror_axis, dtype=np.int32):
-        offsets_dtype = create_type(dtype)
-
-        mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(mirror_axis)
-        mirrored_directions = [stencil.index(MirroredStencilDirections.mirror_stencil(direction, mirror_axis))
-                               for direction in stencil]
-        code = "\n"
-        code += _array_pattern(offsets_dtype, mirrored_stencil_symbol.name, mirrored_directions)
-
-        super(MirroredStencilDirections, self).__init__(code, symbols_read=set(),
-                                                        symbols_defined={mirrored_stencil_symbol})

lbmpy/advanced_streaming/utility.py

@@ -58,24 +58,27 @@ odd_accessors = {
 }
 
 
+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', 'esopull', 'esopush']
+
+
 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 is_inplace(streaming_pattern) and (timestep == Timestep.BOTH):
+        raise ValueError(f"Invalid timestep for streaming pattern {streaming_pattern}: {str(timestep)}")
+
     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', 'esopull', 'esopush']
-
-
 def get_timesteps(streaming_pattern):
     return (Timestep.EVEN, Timestep.ODD) if is_inplace(streaming_pattern) else (Timestep.BOTH, )
 

lbmpy/boundaries/__init__.py

 from lbmpy.boundaries.boundaryconditions import (
     UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
-    ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant, FreeSlip)
+    ExtrapolationOutflow, NeumannByCopy, NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, StreamInConstant, FreeSlip)
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
 
-__all__ = ['NoSlip', 'FreeSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
-           'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
-           'LatticeBoltzmannBoundaryHandling', 'StreamInConstant']
+__all__ = ['NoSlip', 'NoSlipLinearBouzidi', 'QuadraticBounceBack', 'FreeSlip',
+           'UBB', 'FixedDensity',
+           'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
+           'DiffusionDirichlet', 'NeumannByCopy', 'StreamInConstant',
+           'LatticeBoltzmannBoundaryHandling']

lbmpy/boundaries/boundaries_in_kernel.py

 import sympy as sp
 
-from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
 from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
 from lbmpy.advanced_streaming.utility import Timestep, get_accessor
+from lbmpy.custom_code_nodes import LbmWeightInfo
 from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
 from pystencils.assignment import Assignment
 from pystencils.astnodes import Block, Conditional, LoopOverCoordinate, SympyAssignment

lbmpy/boundaries/boundaryconditions.py

 import abc
 from warnings import warn
 
-from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
-from pystencils.simp.assignment_collection import AssignmentCollection
 from pystencils import Assignment, Field
-from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
-from pystencils.typing import create_type
-from pystencils.sympyextensions import get_symmetric_part
-from lbmpy.simplificationfactory import create_simplification_strategy
-from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays, MirroredStencilDirections
+from pystencils.simp.assignment_collection import AssignmentCollection
 from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
+from pystencils.sympyextensions import get_symmetric_part, simplify_by_equality
+from pystencils.typing import create_type, TypedSymbol
+
+from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
+from lbmpy.custom_code_nodes import (NeighbourOffsetArrays, MirroredStencilDirections, LbmWeightInfo,
+                                     TranslationArraysNode)
+from lbmpy.maxwellian_equilibrium import discrete_equilibrium
+from lbmpy.simplificationfactory import create_simplification_strategy
 
 import sympy as sp
 
@@ -111,7 +113,207 @@ class NoSlip(LbBoundary):
         return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
 
 
-# end class NoSlip
+class NoSlipLinearBouzidi(LbBoundary):
+    """
+    No-Slip, (half-way) simple bounce back boundary condition with interpolation
+    to increase accuracy: :cite:`BouzidiBC`. In order to make the boundary condition work properly a
+    Python callback function needs to be provided to calculate the distance from the wall for each cell near to the
+    boundary. If this is not done the boundary condition will fall back to a simple NoSlip boundary.
+    Furthermore, for this boundary condition a second fluid cell away from the wall is needed. If the second fluid
+    cell is not available (e.g. because it is marked as boundary as well), the boundary condition should fall back to
+    a NoSlip boundary as well.
+
+    Args:
+        name: optional name of the boundary.
+        init_wall_distance: Python callback function to calculate the wall distance for each cell near to the  boundary
+        data_type: data type of the wall distance q
+    """
+
+    def __init__(self, name=None, init_wall_distance=None, data_type='double'):
+        self.data_type = data_type
+        self.init_wall_distance = init_wall_distance
+
+        super(NoSlipLinearBouzidi, self).__init__(name)
+
+    @property
+    def additional_data(self):
+        """Used internally only. For the NoSlipLinearBouzidi boundary the distance to the obstacle of every
+        direction is needed. This information is stored in the index vector."""
+        return [('q', create_type(self.data_type))]
+
+    @property
+    def additional_data_init_callback(self):
+        def default_callback(boundary_data, **_):
+            for cell in boundary_data.index_array:
+                cell['q'] = -1
+
+        if self.init_wall_distance:
+            return self.init_wall_distance
+        else:
+            warn("No callback function provided to initialise the wall distance for each cell "
+                 "(init_wall_distance=None). The boundary condition will fall back to a simple NoSlip BC")
+            return default_callback
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        f_xf = sp.Symbol("f_xf")
+        f_xf_inv = sp.Symbol("f_xf_inv")
+        d_x2f = sp.Symbol("d_x2f")
+        q = sp.Symbol("q")
+        one = sp.Float(1.0)
+        two = sp.Float(2.0)
+        half = sp.Rational(1, 2)
+
+        subexpressions = [Assignment(f_xf, f_out(dir_symbol)),
+                          Assignment(f_xf_inv, f_out(inv_dir[dir_symbol])),
+                          Assignment(d_x2f, f_in(dir_symbol)),
+                          Assignment(q, index_field[0]('q'))]
+
+        case_one = (half * (f_xf + f_xf_inv * (two * q - one))) / q
+        case_two = two * q * f_xf + (one - two * q) * d_x2f
+        case_three = f_xf
+
+        rhs = sp.Piecewise((case_one, sp.Ge(q, 0.5)),
+                           (case_two, sp.And(sp.Gt(q, 0), sp.Lt(q, 0.5))),
+                           (case_three, True))
+
+        boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), rhs)]
+
+        return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
+
+
+# end class NoSlipLinearBouzidi
+
+class QuadraticBounceBack(LbBoundary):
+    """
+    Second order accurate bounce back boundary condition. Implementation details are provided in a demo notebook here:
+    https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/demo_interpolation_boundary_conditions.html
+
+    Args:
+        relaxation_rate: relaxation rate to realise a BGK scheme for recovering the pre collision PDF value.
+        name: optional name of the boundary.
+        init_wall_distance: Python callback function to calculate the wall distance for each cell near to the  boundary
+        data_type: data type of the wall distance q
+    """
+
+    def __init__(self, relaxation_rate, name=None, init_wall_distance=None, data_type='double'):
+        self.relaxation_rate = relaxation_rate
+        self.data_type = data_type
+        self.init_wall_distance = init_wall_distance
+        self.equilibrium_values_name = "f_eq"
+        self.inv_dir_symbol = TypedSymbol("inv_dir", create_type("int32"))
+
+        super(QuadraticBounceBack, self).__init__(name)
+
+    @property
+    def additional_data(self):
+        """Used internally only. For the NoSlipLinearBouzidi boundary the distance to the obstacle of every
+        direction is needed. This information is stored in the index vector."""
+        return [('q', create_type(self.data_type))]
+
+    @property
+    def additional_data_init_callback(self):
+        def default_callback(boundary_data, **_):
+            for cell in boundary_data.index_array:
+                cell['q'] = 0.5
+
+        if self.init_wall_distance:
+            return self.init_wall_distance
+        else:
+            warn("No callback function provided to initialise the wall distance for each cell "
+                 "(init_wall_distance=None). The boundary condition will fall back to a simple NoSlip BC")
+            return default_callback
+
+    def get_additional_code_nodes(self, lb_method):
+        """Return a list of code nodes that will be added in the generated code before the index field loop.
+
+        Args:
+            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`
+
+        Returns:
+            list containing LbmWeightInfo
+        """
+        stencil = lb_method.stencil
+        inv_directions = [str(stencil.index(inverse_direction(direction))) for direction in stencil]
+        dtype = self.inv_dir_symbol.dtype
+        name = self.inv_dir_symbol.name
+        inverse_dir_node = TranslationArraysNode([(dtype, name, inv_directions), ], {self.inv_dir_symbol})
+        return [LbmWeightInfo(lb_method, self.data_type), inverse_dir_node, NeighbourOffsetArrays(lb_method.stencil)]
+
+    @staticmethod
+    def get_equilibrium(v, u, rho, drho, weight, compressible, deviation_only):
+        rho_background = sp.Integer(1)
+
+        result = discrete_equilibrium(v, u, rho, weight,
+                                      order=2, c_s_sq=sp.Rational(1, 3), compressible=compressible)
+        if deviation_only:
+            shift = discrete_equilibrium(v, [0] * len(u), rho_background, weight,
+                                         order=0, c_s_sq=sp.Rational(1, 3), compressible=False)
+            result = simplify_by_equality(result - shift, rho, drho, rho_background)
+
+        return result
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        inv = sp.IndexedBase(self.inv_dir_symbol, shape=(1,))[dir_symbol]
+        weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+        weight_of_direction = weight_info.weight_of_direction
+        pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
+
+        pdf_symbols = [sp.Symbol(f"pdf_{i}") for i in range(len(lb_method.stencil))]
+        f_xf = sp.Symbol("f_xf")
+        f_xf_inv = sp.Symbol("f_xf_inv")
+        q = sp.Symbol("q")
+        feq = sp.Symbol("f_eq")
+        weight = sp.Symbol("w")
+        weight_inv = sp.Symbol("w_inv")
+        v = [TypedSymbol(f"c_{i}", self.data_type) for i in range(lb_method.stencil.D)]
+        v_inv = [TypedSymbol(f"c_inv_{i}", self.data_type) for i in range(lb_method.stencil.D)]
+        one = sp.Float(1.0)
+        two = sp.Float(2.0)
+
+        subexpressions = [Assignment(pdf_symbols[i], pdf) for i, pdf in enumerate(pdf_field_accesses)]
+        subexpressions.append(Assignment(f_xf, f_out(dir_symbol)))
+        subexpressions.append(Assignment(f_xf_inv, f_out(inv_dir[dir_symbol])))
+        subexpressions.append(Assignment(q, index_field[0]('q')))
+        subexpressions.append(Assignment(weight, weight_of_direction(dir_symbol, lb_method)))
+        subexpressions.append(Assignment(weight_inv, weight_of_direction(inv, lb_method)))
+
+        for i in range(lb_method.stencil.D):
+            offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
+            subexpressions.append(Assignment(v[i], offset[i]))
+
+        for i in range(lb_method.stencil.D):
+            offset = NeighbourOffsetArrays.neighbour_offset(inv, lb_method.stencil)
+            subexpressions.append(Assignment(v_inv[i], offset[i]))
+
+        cqc = lb_method.conserved_quantity_computation
+        rho = cqc.density_symbol
+        drho = cqc.density_deviation_symbol
+        u = sp.Matrix(cqc.velocity_symbols)
+        compressible = cqc.compressible
+        deviation_only = lb_method.equilibrium_distribution.deviation_only
+
+        cqe = cqc.equilibrium_input_equations_from_pdfs(pdf_symbols, False)
+        subexpressions.append(cqe.all_assignments)
+
+        eq_dir = self.get_equilibrium(v, u, rho, drho, weight, compressible, deviation_only)
+        eq_inv = self.get_equilibrium(v_inv, u, rho, drho, weight_inv, compressible, deviation_only)
+
+        subexpressions.append(Assignment(feq, eq_dir + eq_inv))
+
+        # equation E.4 first summand
+        e41 = (f_xf_inv - f_xf) / two
+        # equation E.4 second summand
+        e42 = (f_xf_inv + f_xf - self.relaxation_rate * feq) / (two - two * self.relaxation_rate)
+        # equation E.3
+        fw = ((one - q) * (e41 + e42)) + q * f_xf_inv
+        # equation E.1
+        result = (one / (q + one)) * fw + (q / (q + one)) * f_xf
+
+        boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), result)]
+        return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
+
+
+# end class QuadraticBounceBack
 
 class FreeSlip(LbBoundary):
     """

lbmpy/boundaries/boundaryhandling.py

 import numpy as np
-import sympy as sp
-from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
-from lbmpy.advanced_streaming.utility import is_inplace, Timestep, AccessPdfValues
-from pystencils import Field, Assignment, TypedSymbol, create_kernel
-from pystencils.stencil import inverse_direction
-from pystencils import CreateKernelConfig, Target
+
+from pystencils import Assignment, CreateKernelConfig, create_kernel, Field, Target
 from pystencils.boundaries import BoundaryHandling
 from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
-from pystencils.backends.cbackend import CustomCodeNode
+from pystencils.simp import add_subexpressions_for_field_reads
+from pystencils.stencil import inverse_direction
+
+from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
+from lbmpy.advanced_streaming.utility import is_inplace, Timestep, AccessPdfValues
 
 
 class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
@@ -154,32 +154,6 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
 
 # end class LatticeBoltzmannBoundaryHandling
 
-
-class LbmWeightInfo(CustomCodeNode):
-    def __init__(self, lb_method, data_type='double'):
-        self.weights_symbol = TypedSymbol("weights", data_type)
-        data_type_string = "double" if self.weights_symbol.dtype.numpy_dtype == np.float64 else "float"
-
-        weights = [str(w.evalf(17)) for w in lb_method.weights]
-        if data_type_string == "float":
-            weights = "f, ".join(weights)
-            weights += "f"  # suffix for the last element
-        else:
-            weights = ", ".join(weights)
-        w_sym = self.weights_symbol
-        code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n"
-        super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
-
-    def weight_of_direction(self, dir_idx, lb_method=None):
-        if isinstance(sp.sympify(dir_idx), sp.Integer):
-            return lb_method.weights[dir_idx].evalf(17)
-        else:
-            return sp.IndexedBase(self.weights_symbol, shape=(1,))[dir_idx]
-
-
-# end class LbmWeightInfo
-
-
 def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method, boundary_functor,
                                              prev_timestep=Timestep.BOTH, streaming_pattern='pull',
                                              target=Target.CPU, **kernel_creation_args):
@@ -194,13 +168,16 @@ def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method,
     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
-
     config = CreateKernelConfig(index_fields=[index_field], target=target, default_number_int="int32",
                                 skip_independence_check=True, **kernel_creation_args)
 
+    default_data_type = config.data_type.default_factory()
+    if pdf_field.dtype != default_data_type:
+        boundary_assignments = add_subexpressions_for_field_reads(boundary_assignments, data_type=default_data_type)
+
+    elements = [Assignment(dir_symbol, index_field[0]('dir'))]
+    elements += boundary_assignments.all_assignments
+
     kernel = create_kernel(elements, config=config)
 
     #   Code Elements ahead of the loop

lbmpy/creationfunctions.py

@@ -112,13 +112,19 @@ class LBMConfig:
     """
     Sequence of relaxation rates, number depends on selected method. If you specify more rates than
     method needs, the additional rates are ignored.
+
+    If no relaxation rates are specified, the parameter `relaxation_rate` will be consulted.
     """
     relaxation_rate: Union[int, float, Type[sp.Symbol]] = None
     """
-    For SRT, TRT and polynomial cumulant models it is possible to define
-    a single ``relaxation_rate`` instead of a list (Internally this is converted to a list with a single entry).
-    The second rate for TRT is then determined via magic number. For the moment, central moment based and the
-    cumulant model, it sets only the relaxation rate corresponding to shear viscosity, setting all others to unity.
+    The method's primary relaxation rate. In most cases, this is the relaxation rate governing shear viscosity.
+    For SRT, this is the only relaxation rate.
+    For TRT, the second relaxation rate is then determined via magic number.
+    In the case of raw moment, central moment, and cumulant-based MRT methods, all other relaxation rates will be
+    set to unity.
+
+    If neither `relaxation_rate` nor `relaxation_rates` is specified, the behaviour is as if 
+    `relaxation_rate=sp.Symbol('omega')` was set.
     """
     compressible: bool = False
     """
@@ -332,9 +338,11 @@ class LBMConfig:
             self.stencil = LBStencil(self.stencil)
 
         if self.relaxation_rates is None:
-            self.relaxation_rates = [sp.Symbol("omega")] * self.stencil.Q
+            #   Fall back to regularized method
+            if self.relaxation_rate is None:
+                self.relaxation_rate = sp.Symbol("omega")
 
-        # if only a single relaxation rate is defined (which makes sense for SRT or TRT methods)
+        # if only a single relaxation rate is defined,
         # it is internally treated as a list with one element and just sets the relaxation_rates parameter
         if self.relaxation_rate is not None:
             if self.method in [Method.TRT, Method.TRT_KBC_N1, Method.TRT_KBC_N2, Method.TRT_KBC_N3, Method.TRT_KBC_N4]:

lbmpy/custom_code_nodes.py

+import numpy as np
+import sympy as sp
+
+from pystencils.typing import TypedSymbol, create_type
+from pystencils.backends.cbackend import CustomCodeNode
+
+
+class NeighbourOffsetArrays(CustomCodeNode):
+
+    @staticmethod
+    def neighbour_offset(dir_idx, stencil):
+        if isinstance(sp.sympify(dir_idx), sp.Integer):
+            return stencil[dir_idx]
+        else:
+            return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
+                         for symbol in NeighbourOffsetArrays._offset_symbols(len(stencil[0]))])
+
+    @staticmethod
+    def _offset_symbols(dim):
+        return [TypedSymbol(f"neighbour_offset_{d}", create_type('int32')) for d in ['x', 'y', 'z'][:dim]]
+
+    def __init__(self, stencil, offsets_dtype=np.int32):
+        offsets_dtype = create_type(offsets_dtype)
+        dim = len(stencil[0])
+
+        array_symbols = NeighbourOffsetArrays._offset_symbols(dim)
+        code = "\n"
+        for i, arrsymb in enumerate(array_symbols):
+            code += _array_pattern(offsets_dtype, arrsymb.name, (d[i] for d in stencil))
+
+        offset_symbols = NeighbourOffsetArrays._offset_symbols(dim)
+        super(NeighbourOffsetArrays, self).__init__(code, symbols_read=set(),
+                                                    symbols_defined=set(offset_symbols))
+
+
+class MirroredStencilDirections(CustomCodeNode):
+
+    @staticmethod
+    def mirror_stencil(direction, mirror_axis):
+        assert mirror_axis <= len(direction), f"only {len(direction)} axis available for mirage"
+        direction = list(direction)
+        direction[mirror_axis] = -direction[mirror_axis]
+
+        return tuple(direction)
+
+    @staticmethod
+    def _mirrored_symbol(mirror_axis):
+        axis = ['x', 'y', 'z']
+        return TypedSymbol(f"{axis[mirror_axis]}_axis_mirrored_stencil_dir", create_type('int32'))
+
+    def __init__(self, stencil, mirror_axis, dtype=np.int32):
+        offsets_dtype = create_type(dtype)
+
+        mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(mirror_axis)
+        mirrored_directions = [stencil.index(MirroredStencilDirections.mirror_stencil(direction, mirror_axis))
+                               for direction in stencil]
+        code = "\n"
+        code += _array_pattern(offsets_dtype, mirrored_stencil_symbol.name, mirrored_directions)
+
+        super(MirroredStencilDirections, self).__init__(code, symbols_read=set(),
+                                                        symbols_defined={mirrored_stencil_symbol})
+
+
+class LbmWeightInfo(CustomCodeNode):
+    def __init__(self, lb_method, data_type='double'):
+        self.weights_symbol = TypedSymbol("weights", data_type)
+        data_type_string = "double" if self.weights_symbol.dtype.numpy_dtype == np.float64 else "float"
+
+        weights = [str(w.evalf(17)) for w in lb_method.weights]
+        if data_type_string == "float":
+            weights = "f, ".join(weights)
+            weights += "f"  # suffix for the last element
+        else:
+            weights = ", ".join(weights)
+        w_sym = self.weights_symbol
+        code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n"
+        super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
+
+    def weight_of_direction(self, dir_idx, lb_method=None):
+        if isinstance(sp.sympify(dir_idx), sp.Integer):
+            return lb_method.weights[dir_idx].evalf(17)
+        else:
+            return sp.IndexedBase(self.weights_symbol, shape=(1,))[dir_idx]
+
+
+class TranslationArraysNode(CustomCodeNode):
+
+    def __init__(self, array_content, symbols_defined):
+        code = ''
+        for content in array_content:
+            code += _array_pattern(*content)
+        super(TranslationArraysNode, self).__init__(code, symbols_read=set(), symbols_defined=symbols_defined)
+
+    def __str__(self):
+        return "Variable PDF Access Translation Arrays"
+
+    def __repr__(self):
+        return "Variable PDF Access Translation Arrays"
+
+
+def _array_pattern(dtype, name, content):
+    return f"const {str(dtype)} {name} [] = {{ {','.join(str(c) for c in content)} }}; \n"

lbmpy/maxwellian_equilibrium.py

@@ -76,38 +76,49 @@ def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=sp.symbols(
     """
     weights = get_weights(stencil, c_s_sq)
     assert stencil.Q == len(weights)
-
     u = u[:stencil.D]
 
+    res = [discrete_equilibrium(e_q, u, rho, w_q, order, c_s_sq, compressible) for w_q, e_q in zip(weights, stencil)]
+    return tuple(res)
+
+
+def discrete_equilibrium(v=sp.symbols("v_:3"), u=sp.symbols("u_:3"), rho=sp.Symbol("rho"), weight=sp.Symbol("w"),
+                         order=2, c_s_sq=sp.Symbol("c_s") ** 2, compressible=True):
+    """
+    Returns the common discrete LBM equilibrium as a list of sympy expressions
+
+    Args:
+        v: symbols for mesoscopic velocity
+        u: symbols for macroscopic velocity
+        rho: sympy symbol for the density
+        weight: symbol for stencil weights
+        order: highest order of velocity terms (for hydrodynamics order 2 is sufficient)
+        c_s_sq: square of speed of sound
+        compressible: compressibility
+    """
     rho_outside = rho if compressible else sp.Rational(1, 1)
     rho_inside = rho if not compressible else sp.Rational(1, 1)
 
-    res = []
-    for w_q, e_q in zip(weights, stencil):
-        e_times_u = 0
-        for c_q_alpha, u_alpha in zip(e_q, u):
-            e_times_u += c_q_alpha * u_alpha
-
-        fq = rho_inside + e_times_u / c_s_sq
+    e_times_u = 0
+    for c_q_alpha, u_alpha in zip(v, u):
+        e_times_u += c_q_alpha * u_alpha
 
-        if order <= 1:
-            res.append(fq * rho_outside * w_q)
-            continue
+    fq = rho_inside + e_times_u / c_s_sq
 
-        u_times_u = 0
-        for u_alpha in u:
-            u_times_u += u_alpha * u_alpha
-        fq += sp.Rational(1, 2) / c_s_sq**2 * e_times_u ** 2 - sp.Rational(1, 2) / c_s_sq * u_times_u
+    if order <= 1:
+        return fq * rho_outside * weight
 
-        if order <= 2:
-            res.append(fq * rho_outside * w_q)
-            continue
+    u_times_u = 0
+    for u_alpha in u:
+        u_times_u += u_alpha * u_alpha
+    fq += sp.Rational(1, 2) / c_s_sq**2 * e_times_u ** 2 - sp.Rational(1, 2) / c_s_sq * u_times_u
 
-        fq += sp.Rational(1, 6) / c_s_sq**3 * e_times_u**3 - sp.Rational(1, 2) / c_s_sq**2 * u_times_u * e_times_u
+    if order <= 2:
+        return fq * rho_outside * weight
 
-        res.append(sp.expand(fq * rho_outside * w_q))
+    fq += sp.Rational(1, 6) / c_s_sq**3 * e_times_u**3 - sp.Rational(1, 2) / c_s_sq**2 * u_times_u * e_times_u
 
-    return tuple(res)
+    return sp.expand(fq * rho_outside * weight)
 
 
 @disk_cache

lbmpy/simplificationfactory.py

@@ -92,6 +92,7 @@ def _cumulant_space_simplification(split_inner_loop):
     s.add(expand_post_collision_central_moments)
     s.add(insert_aliases)
     s.add(insert_constants)
+    s.add(insert_aliases)
     s.add(add_subexpressions_for_divisions)
     s.add(add_subexpressions_for_constants)
     if split_inner_loop: