Still buggy periodic pipe test case

lbmpy_tests/advanced_streaming/test_fully_periodic.ipynb

@@ -63,7 +63,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -71,7 +71,7 @@
     "q = len(stencil)\n",
     "dim = len(stencil[0])\n",
     "\n",
-    "streaming_pattern = 'esotwist'\n",
+    "streaming_pattern = 'push'\n",
     "inplace = is_inplace(streaming_pattern)\n",
     "timesteps = ['even', 'odd'] if inplace else ['both']\n",
     "zeroth_timestep = timesteps[0]"
@@ -86,7 +86,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -109,7 +109,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -151,6 +151,69 @@
     "                                         **method_params))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from lbmpy.creationfunctions import create_lb_update_rule\n",
+    "push_update = create_lb_update_rule(\n",
+    "                collision_rule=lb_collision, optimization=optimization, kernel_type='collide_stream_push', **method_params)\n",
+    "esotwist_update = create_lb_update_rule(\n",
+    "                    collision_rule=lb_collision, optimization=optimization, kernel_type='esotwist_even', **method_params)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": [
+       "Assignment(pdfs_tmp_N^1, pdfs_C^1 + rr_0*(-pdfs_C^1 + f_eq_common/9 + u_1**2/2 + u_1/3))"
+      ],
+      "text/latex": "$\\displaystyle {{pdfs_tmp}_{(0,1)}^{1}} \\leftarrow {{pdfs}_{(0,0)}^{1}} + rr_{0} \\left(- {{pdfs}_{(0,0)}^{1}} + \\frac{f_{eq common}}{9} + \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{3}\\right)$"
+     },
+     "metadata": {},
+     "execution_count": 24
+    }
+   ],
+   "source": [
+    "push_update.main_assignments[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": [
+       "Assignment(pdfs_N^2, rr_0*(f_eq_common/9 + u_1**2/2 + u_1/3 - xi_8) + xi_8)"
+      ],
+      "text/latex": "$\\displaystyle {{pdfs}_{(0,1)}^{2}} \\leftarrow rr_{0} \\left(\\frac{f_{eq common}}{9} + \\frac{u_{1}^{2}}{2} + \\frac{u_{1}}{3} - \\xi_{8}\\right) + \\xi_{8}$"
+     },
+     "metadata": {},
+     "execution_count": 25
+    }
+   ],
+   "source": [
+    "esotwist_update.main_assignments[1]"
+   ]
+  },
+  {
+   "source": [
+    "...damnit!\n",
+    "`pystencils` does different algebraic modifications to the update equations depending on the streaming pattern. Floating-point addition is not associative! This causes minor numerical differences between the streaming patterns, as visible in the test case."
+   ],
+   "cell_type": "markdown",
+   "metadata": {}
+  },
   {
    "source": [
     "## Macroscopic Values"

lbmpy_tests/advanced_streaming/test_fully_periodic_flow.py

@@ -146,8 +146,11 @@ def test_fully_periodic_flow(stencil, streaming_pattern):
     #   Equal to the steady-state velocity field up to numerical errors
     assert_allclose(u, u_ref)
 
-    #   Exactly equal to other streaming patterns!
+    #   Flow must be equal up to numerical error for all streaming patterns
     global all_results
-    for prev_pattern, prev_u in all_results.items():
-        assert_array_equal(
-            u, prev_u, err_msg=f'Velocity field for {streaming_pattern} differed from {prev_pattern}!')
+    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
+
+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 = ['even', 'odd'] if self.inplace else ['both']
+        self.between_timesteps = ['even_to_odd', 'odd_to_even'] if self.inplace else ['both']
+        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(
+            stencil, self.dh, self.pdfs.name,
+            streaming_pattern=self.streaming_pattern, zeroth_timestep=self.zeroth_timestep)
+
+        #   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.t_modulus = 0
+
+    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.t_modulus = 0
+        self.dh.run_kernel(self.init_kernel)
+
+    def step(self):
+        # Boundaries
+        self.bh(between_timesteps=self.between_timesteps[self.t_modulus])
+
+        # Periodicty
+        self.periodicity_handler(self.timesteps[self.t_modulus])
+
+        # Macroscopic Values
+        self.dh.run_kernel(self.getter_kernels[self.t_modulus])
+
+        # Here, the next time step begins
+        self.t_modulus = (self.t_modulus + 1) % len(self.timesteps)
+
+        # LBM Step
+        self.dh.run_kernel(self.lb_kernels[self.t_modulus])
+
+        # Field Swaps
+        if not self.inplace:
+            self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
+
+    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,
+                err_msg=f'Velocity field for {streaming_pattern} differed from {prev_pattern}!')
+    all_results[(stencil, streaming_pattern)] = u