First working draft

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 .custom_code_nodes import NeighbourOffsetArrays, MirroredStencilDirections
+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', 'NeighbourOffsetArrays', 'MirroredStencilDirections',
            'get_communication_slices', 'LBMPeriodicityHandling',
            'Timestep', 'get_accessor', 'is_inplace', 'get_timesteps',
            'numeric_index', 'numeric_offsets', 'inverse_dir_index', 'AccessPdfValues']

lbmpy/advanced_streaming/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/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 .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:
 
     #   ==============================================
@@ -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/boundaries/__init__.py

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

lbmpy/boundaries/boundaryconditions.py

@@ -4,11 +4,10 @@ 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 lbmpy.advanced_streaming.custom_code_nodes import NeighbourOffsetArrays, MirroredStencilDirections, LbmWeightInfo
 from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
 
 import sympy as sp
@@ -604,6 +603,66 @@ class FixedDensity(LbBoundary):
 
 # end class FixedDensity
 
+class NoslipEquilibriumReconstruction(LbBoundary):
+
+    def __init__(self, normal_direction, velocity=None, name=None):
+        if isinstance(normal_direction, str):
+            normal_direction = direction_string_to_offset(normal_direction, dim=len(normal_direction))
+
+        if name is None:
+            name = f"Noslip ER: {offset_to_direction_string(normal_direction)}"
+
+        if velocity is None:
+            self.velocity = tuple([0] * len(normal_direction))
+        else:
+            assert len(normal_direction) == len(velocity), "normal_direction and velocity must have same size"
+
+        self.normal_direction = tuple([int(n) for n in normal_direction])
+        self.velocity = velocity
+        assert all([n in [-1, 0, 1] for n in self.normal_direction]), \
+            "Only -1, 0 and 1 allowed for defining the normal direction"
+        super(NoslipEquilibriumReconstruction, self).__init__(name)
+
+    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 NeighbourOffsetArrays
+
+        """
+        return [NeighbourOffsetArrays(lb_method.stencil)]
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
+        tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction))
+
+        cqc = lb_method.conserved_quantity_computation
+        density_symbol = sp.Symbol("rho")
+        pdf_field_accesses = [f_out[tangential_offset](i) for i in range(len(lb_method.stencil))]
+        density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
+
+        substitutions = {u_symbolic: u for u_symbolic, u in zip(cqc.velocity_symbols, self.velocity)}
+
+        conditions = list()
+        equilibrium_terms = lb_method.get_equilibrium_terms().subs(substitutions)
+        for i, term in enumerate(equilibrium_terms[1:]):
+            conditions.append((term, sp.Eq(dir_symbol, sp.Integer(i + 1))))
+
+        expr = sp.Piecewise(*conditions, (equilibrium_terms[0], True))
+
+        subexpressions = density_equations.all_assignments
+        main_assignments = [Assignment(f_in(inv_dir[dir_symbol]), expr)]
+
+        ac = AssignmentCollection(main_assignments, subexpressions=subexpressions)
+        ac = ac.new_without_unused_subexpressions()
+        ac.topological_sort()
+        return ac
+
+# end class NoslipEquilibriumReconstruction
+
 
 class DiffusionDirichlet(LbBoundary):
     """Boundary condition for advection-diffusion problems that fixes the concentration at the obstacle.

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 import Field, Assignment, create_kernel
 from pystencils.stencil import inverse_direction
 from pystencils import CreateKernelConfig, Target
 from pystencils.boundaries import BoundaryHandling
 from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
-from pystencils.backends.cbackend import CustomCodeNode
 
 
 class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
@@ -152,34 +150,6 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
         return dh.reduce_float_sequence(list(result), 'sum')
 
 
-# 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):