Short Straming pattern demo

doc/conf.py

@@ -27,7 +27,7 @@ source_suffix = '.rst'
 master_doc = 'index'
 
 copyright = f'{datetime.datetime.now().year}, Martin Bauer, Markus Holzer'
-author = 'Martin Bauer, Markus Holzer'
+author = 'Martin Bauer, Markus Holzer, Frederik Hennig'
 # The short X.Y version (including .devXXXX, rcX, b1 suffixes if present)
 version = re.sub(r'(\d+\.\d+)\.\d+(.*)', r'\1\2', lbmpy.__version__)
 version = re.sub(r'(\.dev\d+).*?$', r'\1', version)

doc/sphinx/tutorials.rst

@@ -20,6 +20,7 @@ You can open the notebooks directly to play around with the code examples.
     /notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
     /notebooks/11_tutorial_Non_Newtonian_Flow.ipynb
     /notebooks/demo_stencils.ipynb
+    /notebooks/demo_streaming_patterns.ipynb
     /notebooks/demo_create_method_from_scratch.ipynb
     /notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
     /notebooks/demo_automatic_chapman_enskog_analysis.ipynb

lbmpy/advanced_streaming/utility.py

@@ -4,7 +4,11 @@ from lbmpy.fieldaccess import PdfFieldAccessor, \
     AAEvenTimeStepAccessor, \
     AAOddTimeStepAccessor, \
     EsoTwistEvenTimeStepAccessor, \
-    EsoTwistOddTimeStepAccessor
+    EsoTwistOddTimeStepAccessor, \
+    EsoPullEvenTimeStepAccessor, \
+    EsoPullOddTimeStepAccessor, \
+    EsoPushEvenTimeStepAccessor, \
+    EsoPushOddTimeStepAccessor
 
 import numpy as np
 import pystencils as ps
@@ -33,20 +37,24 @@ class Timestep(IntEnum):
             return 'Both'
 
 
-streaming_patterns = ['push', 'pull', 'aa', 'esotwist']
+streaming_patterns = ['push', 'pull', 'aa', 'esotwist', 'esopull', 'esopush']
 
 even_accessors = {
     'pull': StreamPullTwoFieldsAccessor,
     'push': StreamPushTwoFieldsAccessor,
     'aa': AAEvenTimeStepAccessor,
-    'esotwist': EsoTwistEvenTimeStepAccessor
+    'esotwist': EsoTwistEvenTimeStepAccessor,
+    'esopull': EsoPullEvenTimeStepAccessor,
+    'esopush': EsoPushEvenTimeStepAccessor
 }
 
 odd_accessors = {
     'pull': StreamPullTwoFieldsAccessor,
     'push': StreamPushTwoFieldsAccessor,
     'aa': AAOddTimeStepAccessor,
-    'esotwist': EsoTwistOddTimeStepAccessor
+    'esotwist': EsoTwistOddTimeStepAccessor,
+    'esopull': EsoPullOddTimeStepAccessor,
+    'esopush': EsoPushOddTimeStepAccessor
 }
 
 
@@ -65,7 +73,7 @@ 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']
+    return streaming_pattern in ['aa', 'esotwist', 'esopull', 'esopush']
 
 
 def get_timesteps(streaming_pattern):

lbmpy/fieldaccess.py

@@ -14,6 +14,8 @@ __all__ = ['PdfFieldAccessor', 'CollideOnlyInplaceAccessor', 'StreamPullTwoField
            'AAEvenTimeStepAccessor', 'AAOddTimeStepAccessor',
            'PeriodicTwoFieldsAccessor', 'StreamPushTwoFieldsAccessor',
            'EsoTwistEvenTimeStepAccessor', 'EsoTwistOddTimeStepAccessor',
+           'EsoPullEvenTimeStepAccessor', 'EsoPullOddTimeStepAccessor',
+           'EsoPushEvenTimeStepAccessor', 'EsoPushOddTimeStepAccessor',
            'visualize_pdf_field_accessor', 'visualize_field_mapping']
 
 
@@ -140,7 +142,7 @@ class AAEvenTimeStepAccessor(PdfFieldAccessor):
 
     @staticmethod
     def write(field, stencil):
-        return [field(stencil.index(inverse_direction(d))) for d in stencil]
+        return [field(stencil.inverse_index(d)) for d in stencil]
 
 
 class AAOddTimeStepAccessor(PdfFieldAccessor):
@@ -148,57 +150,169 @@ class AAOddTimeStepAccessor(PdfFieldAccessor):
 
     @staticmethod
     def read(field, stencil):
-        res = []
-        for i, d in enumerate(stencil):
-            inv_dir = inverse_direction(d)
-            field_access = field[inv_dir](stencil.index(inv_dir))
-            res.append(field_access)
-        return res
+        return [field[inverse_direction(d)](stencil.inverse_index(d)) for i, d in enumerate(stencil)]
 
     @staticmethod
     def write(field, stencil):
         return [field[d](i) for i, d in enumerate(stencil)]
 
 
+class EsoTwistEvenTimeStepAccessor(PdfFieldAccessor):
+    is_inplace = True
+
+    @staticmethod
+    def read(field, stencil):
+        return [field[tuple(max(-e, 0) for e in d)](i) for i, d in enumerate(stencil)]
+
+    @staticmethod
+    def write(field, stencil):
+        return [field[tuple(max(e, 0) for e in d)](stencil.inverse_index(d)) for d in stencil]
+
+
 class EsoTwistOddTimeStepAccessor(PdfFieldAccessor):
     is_inplace = True
 
     @staticmethod
     def read(field, stencil):
-        result = []
-        for direction in stencil:
-            inv_dir = inverse_direction(direction)
-            spatial_offset = tuple(max(e, 0) for e in inv_dir)
-            result.append(field[spatial_offset](stencil.index(inv_dir)))
+        return [field[tuple(max(e, 0) for e in inverse_direction(d))](stencil.inverse_index(d)) for d in stencil]
+
+    @staticmethod
+    def write(field, stencil):
+        return [field[tuple(max(e, 0) for e in d)](i) for i, d in enumerate(stencil)]
+
+
+class EsoPullEvenTimeStepAccessor(PdfFieldAccessor):
+    is_inplace = True
+
+    @staticmethod
+    def read(field, stencil):
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[inverse_direction(d)](i))
+            else:
+                result.append(field[center_cell](i))
         return result
 
     @staticmethod
     def write(field, stencil):
-        result = []
-        for i, direction in enumerate(stencil):
-            spatial_offset = tuple(max(e, 0) for e in direction)
-            result.append(field[spatial_offset](i))
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[center_cell](stencil.inverse_index(d)))
+            else:
+                result.append(field[d](stencil.inverse_index(d)))
+
         return result
 
 
-class EsoTwistEvenTimeStepAccessor(PdfFieldAccessor):
+class EsoPullOddTimeStepAccessor(PdfFieldAccessor):
     is_inplace = True
 
     @staticmethod
     def read(field, stencil):
-        result = []
-        for i, direction in enumerate(stencil):
-            spatial_offset = tuple(max(-e, 0) for e in direction)
-            result.append(field[spatial_offset](i))
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[inverse_direction(d)](stencil.inverse_index(d)))
+            else:
+                result.append(field[center_cell](stencil.inverse_index(d)))
         return result
 
     @staticmethod
     def write(field, stencil):
-        result = []
-        for direction in stencil:
-            inv_dir = inverse_direction(direction)
-            spatial_offset = tuple(max(e, 0) for e in direction)
-            result.append(field[spatial_offset](stencil.index(inv_dir)))
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[center_cell](i))
+            else:
+                result.append(field[d](i))
+
+        return result
+
+
+class EsoPushEvenTimeStepAccessor(PdfFieldAccessor):
+    is_inplace = True
+
+    @staticmethod
+    def read(field, stencil):
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[center_cell](stencil.inverse_index(d)))
+            else:
+                result.append(field[inverse_direction(d)](stencil.inverse_index(d)))
+
+        return result
+
+    @staticmethod
+    def write(field, stencil):
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[d](i))
+            else:
+                result.append(field[center_cell](i))
+
+        return result
+
+
+class EsoPushOddTimeStepAccessor(PdfFieldAccessor):
+    is_inplace = True
+
+    @staticmethod
+    def read(field, stencil):
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            inv_dir = inverse_direction(d)
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[center_cell](i))
+            else:
+                result.append(field[inv_dir](i))
+
+        return result
+
+    @staticmethod
+    def write(field, stencil):
+        lehmann_stencil = _get_lehmann_stencil(stencil)
+        center_cell = tuple([0] * stencil.D)
+        result = [field.center]
+        for i, d in enumerate(stencil):
+            if i == 0:
+                continue
+            if lehmann_stencil.index(d) % 2 == 0:
+                result.append(field[d](stencil.inverse_index(d)))
+            else:
+                result.append(field[center_cell](stencil.inverse_index(d)))
+
         return result
 
 
@@ -250,3 +364,25 @@ def visualize_pdf_field_accessor(pdf_field_accessor, title=True, read_plot_param
     if title:
         ax_left.set_title("Read")
         ax_right.set_title("Write")
+
+# -------------------------------------------- Helpers -----------------------------------------------------------
+
+
+def _get_lehmann_stencil(stencil):
+    """
+    EsoPull and EsoPush streaming is only simple to implement with a specific stencil ordering, that comes from
+    "High Performance Free Surface LBM on GPUs" by moritz lehmann
+
+    Args:
+        stencil: lattice Boltzmann stencil
+    """
+    if stencil.Q == 9:
+        return LBStencil(Stencil.D2Q9, ordering="lehmann")
+    elif stencil.Q == 15:
+        return LBStencil(Stencil.D3Q15, ordering="lehmann")
+    elif stencil.Q == 19:
+        return LBStencil(Stencil.D3Q19, ordering="lehmann")
+    elif stencil.Q == 27:
+        return LBStencil(Stencil.D3Q27, ordering="lehmann")
+    else:
+        ValueError("EsoPull or EsoPush is only available for D2Q9, D3Q15, D3Q19 and D3Q27 stencil")

lbmpy/maxwellian_equilibrium.py

@@ -12,6 +12,7 @@ In particular, the continuous and discrete Maxwellians are now represented by
 import warnings
 
 import sympy as sp
+from sympy.core.numbers import Zero
 
 from pystencils.cache import disk_cache
 from pystencils.sympyextensions import remove_higher_order_terms
@@ -37,7 +38,7 @@ get_weights.weights = {
         2: sp.Rational(1, 36),
     },
     7: {
-        0: sp.simplify(0.0),
+        0: Zero(),
         1: sp.Rational(1, 6),
     },
     15: {

lbmpy/stencils.py

@@ -56,6 +56,10 @@ class LBStencil:
             raise ValueError("The stencil you have created is not valid. "
                              "It probably contains elements with different lengths")
 
+        if len(set(self._stencil_entries)) < len(self._stencil_entries):
+            raise ValueError("The stencil you have created is not valid. "
+                             "It contains duplicated elements")
+
         self._ordering = ordering
         self._dim = len(self._stencil_entries[0])
         self._q = len(self._stencil_entries)
@@ -87,8 +91,14 @@ class LBStencil:
     def plot(self, slice=False, **kwargs):
         ps.stencil.plot(stencil=self._stencil_entries, slice=slice, **kwargs)
 
-    def index(self, index):
-        return self._stencil_entries.index(index)
+    def index(self, direction):
+        assert len(direction) == self.D, "direction must match stencil.D"
+        return self._stencil_entries.index(direction)
+
+    def inverse_index(self, direction):
+        assert len(direction) == self.D, "direction must match stencil.D"
+        direction = ps.stencil.inverse_direction(direction)
+        return self._stencil_entries.index(direction)
 
     def __getitem__(self, index):
         return self._stencil_entries[index]
@@ -112,6 +122,7 @@ class LBStencil:
         table = """
         <table style="border:none; width: 100%">
             <tr {nb}>
+                <th {nb} >Nr.</th>
                 <th {nb} >Direction Name</th>
                 <th {nb} >Direction </th>
             </tr>
@@ -119,13 +130,15 @@ class LBStencil:
         </table>
         """
         content = ""
-        for direction in self._stencil_entries:
+        for i, direction in enumerate(self._stencil_entries):
             vals = {
+                'nr': sp.latex(i),
                 'name': sp.latex(ps.stencil.offset_to_direction_string(direction)),
                 'entry': sp.latex(direction),
                 'nb': 'style="border:none"'
             }
             content += """<tr {nb}>
+                            <td {nb}>${nr}$</td>
                             <td {nb}>${name}$</td>
                             <td {nb}>${entry}$</td>
                          </tr>\n""".format(**vals)
@@ -147,7 +160,11 @@ def _predefined_stencils(stencil: str, ordering: str):
             'uk': ((0, 0),
                    (1, 0), (-1, 0), (0, 1), (0, -1),
                    (1, 1), (-1, -1), (-1, 1), (1, -1),
-                   )
+                   ),
+            'lehmann': ((0, 0),
+                        (1, 0), (-1, 0), (0, 1), (0, -1),
+                        (1, 1), (-1, -1), (1, -1), (-1, 1),
+                        )
         },
         'D2V17': {
             'walberla': (
@@ -170,6 +187,7 @@ def _predefined_stencils(stencil: str, ordering: str):
                          (0, 1, 0), (0, -1, 0),
                          (-1, 0, 0), (1, 0, 0),
                          (0, 0, 1), (0, 0, -1))
+
         },
         'D3Q15': {
             'walberla':
@@ -181,6 +199,10 @@ def _predefined_stencils(stencil: str, ordering: str):
                          (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
                          (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
                          (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
+            'lehmann': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 1), (-1, -1, -1), (1, 1, -1), (-1, -1, 1),
+                        (1, -1, 1), (-1, 1, -1), (-1, 1, 1), (1, -1, -1)),
         },
         'D3Q19': {
             'walberla': ((0, 0, 0),
@@ -188,21 +210,26 @@ def _predefined_stencils(stencil: str, ordering: str):
                          (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
                          (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
                          (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1)),
+            'counterclockwise': ((0, 0, 0),
+                                 (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                                 (1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
+                                 (0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
+                                 (1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1)),
             'braunschweig': ((0, 0, 0),
-                             (1, 0, 0), (-1, 0, 0),
-                             (0, 1, 0), (0, -1, 0),
-                             (0, 0, 1), (0, 0, -1),
-                             (1, 1, 0), (-1, -1, 0),
-                             (1, -1, 0), (-1, 1, 0),
-                             (1, 0, 1), (-1, 0, -1),
-                             (1, 0, -1), (-1, 0, 1),
-                             (0, 1, 1), (0, -1, -1),
-                             (0, 1, -1), (0, -1, 1)),
+                             (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                             (1, 1, 0), (-1, -1, 0), (1, -1, 0), (-1, 1, 0),
+                             (1, 0, 1), (-1, 0, -1), (1, 0, -1), (-1, 0, 1),
+                             (0, 1, 1), (0, -1, -1), (0, 1, -1), (0, -1, 1)),
             'premnath': ((0, 0, 0),
                          (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
                          (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
                          (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
                          (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
+            'lehmann': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
+                        (0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
+                        (1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1)),
         },
         'D3Q27': {
             'walberla': ((0, 0, 0),
@@ -226,6 +253,13 @@ def _predefined_stencils(stencil: str, ordering: str):
                         (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
                         (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1),
                         (0, -1, -1)),
+            'lehmann': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
+                        (0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
+                        (1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1),
+                        (1, 1, 1), (-1, -1, -1), (1, 1, -1), (-1, -1, 1), (1, -1, 1), (-1, 1, -1), (-1, 1, 1),
+                        (1, -1, -1)),
         }
     }
 

lbmpy_tests/test_version_string.py

@@ -5,7 +5,7 @@ def test_version_string():
     version = lbmpy.__version__
     print(version)
 
-    numeric_version = [int(x, 10) for x in version.split('.')[0:2]]
+    numeric_version = [int(x, 10) for x in version.split('.')[0:1]]
     test_version = sum(x * (100 ** i) for i, x in enumerate(numeric_version))
 
-    assert test_version >= 200
+    assert test_version >= 1

setup.py

@@ -84,7 +84,7 @@ setup(name='lbmpy',
       description='Code Generation for Lattice Boltzmann Methods',
       long_description=readme(),
       long_description_content_type="text/markdown",
-      author='Martin Bauer, Markus Holzer',
+      author='Martin Bauer, Markus Holzer, Frederik Hennig',
       license='AGPLv3',
       author_email='cs10-codegen@fau.de',
       url='https://i10git.cs.fau.de/pycodegen/lbmpy/',