From 14cac68c6742c450e01d7c228891946f72b33592 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Wed, 19 Jul 2023 11:30:46 +0200
Subject: [PATCH] First working draft
---
lbmpy/__init__.py | 2 +-
lbmpy/advanced_streaming/__init__.py | 5 +-
lbmpy/advanced_streaming/custom_code_nodes.py | 102 +++++++++++++++++
lbmpy/advanced_streaming/indexing.py | 104 +++---------------
lbmpy/boundaries/__init__.py | 5 +-
lbmpy/boundaries/boundaryconditions.py | 63 ++++++++++-
lbmpy/boundaries/boundaryhandling.py | 32 +-----
7 files changed, 186 insertions(+), 127 deletions(-)
create mode 100644 lbmpy/advanced_streaming/custom_code_nodes.py
diff --git a/lbmpy/__init__.py b/lbmpy/__init__.py
index 2ad617ec..26d7778f 100644
--- a/lbmpy/__init__.py
+++ b/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',
diff --git a/lbmpy/advanced_streaming/__init__.py b/lbmpy/advanced_streaming/__init__.py
index c88e3463..aeddb1f3 100644
--- a/lbmpy/advanced_streaming/__init__.py
+++ b/lbmpy/advanced_streaming/__init__.py
@@ -1,9 +1,10 @@
-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']
diff --git a/lbmpy/advanced_streaming/custom_code_nodes.py b/lbmpy/advanced_streaming/custom_code_nodes.py
new file mode 100644
index 00000000..0bebefc6
--- /dev/null
+++ b/lbmpy/advanced_streaming/custom_code_nodes.py
@@ -0,0 +1,102 @@
+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"
diff --git a/lbmpy/advanced_streaming/indexing.py b/lbmpy/advanced_streaming/indexing.py
index d953b682..af0af426 100644
--- a/lbmpy/advanced_streaming/indexing.py
+++ b/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})
diff --git a/lbmpy/boundaries/__init__.py b/lbmpy/boundaries/__init__.py
index 6dd746dd..cc9b28bb 100644
--- a/lbmpy/boundaries/__init__.py
+++ b/lbmpy/boundaries/__init__.py
@@ -1,8 +1,9 @@
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']
diff --git a/lbmpy/boundaries/boundaryconditions.py b/lbmpy/boundaries/boundaryconditions.py
index d1bb30d8..e894f1a1 100644
--- a/lbmpy/boundaries/boundaryconditions.py
+++ b/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.
diff --git a/lbmpy/boundaries/boundaryhandling.py b/lbmpy/boundaries/boundaryhandling.py
index 8a94802d..c6dffb27 100644
--- a/lbmpy/boundaries/boundaryhandling.py
+++ b/lbmpy/boundaries/boundaryhandling.py
@@ -1,13 +1,11 @@
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):
--
GitLab