Fix FreeSlip boundary condition

lbmpy/boundaries/boundaryconditions.py

@@ -120,9 +120,10 @@ class FreeSlip(LbBoundary):
 
     Args:
         stencil: LBM stencil which is used for the simulation
-        normal_direction: optional normal direction. If the Free slip boundary is applied to a certain side in the
-                          domain it is not necessary to calculate the normal direction since it can be stated for all
-                          boundary cells. This reduces the memory space for the index array significantly.
+        normal_direction: optional normal direction pointing from wall to fluid.
+                          If the Free slip boundary is applied to a certain side in the domain it is not necessary
+                          to calculate the normal direction since it can be stated for all boundary cells.
+                          This reduces the memory space for the index array significantly.
         name: optional name of the boundary.
     """
 
@@ -182,7 +183,12 @@ class FreeSlip(LbBoundary):
                     normal_direction[i] = direction[i]
                     ref_direction = MirroredStencilDirections.mirror_stencil(ref_direction, i)
 
-            ref_direction = inverse_direction(ref_direction)
+            # convex corner special case:
+            if all(n == 0 for n in normal_direction):
+                normal_direction = direction
+            else:
+                ref_direction = inverse_direction(ref_direction)
+
             for i, cell_name in zip(range(dim), self.additional_data):
                 cell[cell_name[0]] = -normal_direction[i]
             cell['ref_dir'] = self.stencil.index(ref_direction)
@@ -208,13 +214,14 @@ class FreeSlip(LbBoundary):
 
     def get_additional_code_nodes(self, lb_method):
         if self.normal_direction:
-            return [MirroredStencilDirections(self.stencil, self.mirror_axis)]
+            return [MirroredStencilDirections(self.stencil, self.mirror_axis), NeighbourOffsetArrays(lb_method.stencil)]
         else:
-            return []
+            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)
         if self.normal_direction:
-            normal_direction = self.normal_direction
+            tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction))
             mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
             mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
         else:
@@ -222,10 +229,11 @@ class FreeSlip(LbBoundary):
             for i, cell_name in zip(range(self.dim), self.additional_data):
                 normal_direction.append(index_field[0](cell_name[0]))
             normal_direction = tuple(normal_direction)
+            tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, normal_direction))
 
             mirrored_direction = index_field[0]('ref_dir')
 
-        return Assignment(f_in(inv_dir[dir_symbol]), f_in[normal_direction](mirrored_direction))
+        return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](mirrored_direction))
 
 
 # end class FreeSlip

lbmpy_tests/test_boundary_handling.py

@@ -112,6 +112,122 @@ def test_simple(target):
     assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 8] == 5))
 
 
+@pytest.mark.parametrize("given_normal", [True, False])
+def test_free_slip(given_normal):
+    # check if Free slip BC is applied correctly
+
+    stencil = LBStencil(Stencil.D2Q9)
+    dh = create_data_handling(domain_size=(4, 4),)
+    src1 = dh.add_array('src1', values_per_cell=stencil.Q)
+    dh.fill('src1', 0.0, ghost_layers=True)
+
+    shape = dh.gather_array('src1', ghost_layers=True).shape
+
+    num = 0
+    for x in range(shape[0]):
+        for y in range(shape[1]):
+            for direction in range(shape[2]):
+                dh.cpu_arrays[src1.name][x, y, direction] = num
+                num += 1
+
+    method = create_lb_method(lbm_config=LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.8))
+
+    bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src1', name="bh1")
+    if given_normal:
+        free_slipN = FreeSlip(stencil=stencil, normal_direction=(0, -1))
+        free_slipS = FreeSlip(stencil=stencil, normal_direction=(0, 1))
+        free_slipE = FreeSlip(stencil=stencil, normal_direction=(-1, 0))
+        free_slipW = FreeSlip(stencil=stencil, normal_direction=(1, 0))
+
+        bh.set_boundary(free_slipN, slice_from_direction('N', dh.dim))
+        bh.set_boundary(free_slipS, slice_from_direction('S', dh.dim))
+        bh.set_boundary(free_slipE, slice_from_direction('E', dh.dim))
+        bh.set_boundary(free_slipW, slice_from_direction('W', dh.dim))
+    else:
+        free_slip = FreeSlip(stencil=stencil)
+
+        bh.set_boundary(free_slip, slice_from_direction('N', dh.dim))
+        bh.set_boundary(free_slip, slice_from_direction('S', dh.dim))
+        bh.set_boundary(free_slip, slice_from_direction('E', dh.dim))
+        bh.set_boundary(free_slip, slice_from_direction('W', dh.dim))
+
+    bh()
+
+    mirrored_dirN = {6: 8, 1: 2, 5: 7}
+    mirrored_dirS = {7: 5, 2: 1, 8: 6}
+    mirrored_dirE = {6: 5, 4: 3, 8: 7}
+    mirrored_dirW = {5: 6, 3: 4, 7: 8}
+
+    # check North
+    assert dh.cpu_arrays[src1.name][1, -1, mirrored_dirN[6]] == dh.cpu_arrays[src1.name][1, -2, 6]
+    assert dh.cpu_arrays[src1.name][1, -1, mirrored_dirN[1]] == dh.cpu_arrays[src1.name][1, -2, 1]
+
+    for i in range(2, 4):
+        assert dh.cpu_arrays[src1.name][i, -1, mirrored_dirN[6]] == dh.cpu_arrays[src1.name][i, -2, 6]
+        assert dh.cpu_arrays[src1.name][i, -1, mirrored_dirN[1]] == dh.cpu_arrays[src1.name][i, -2, 1]
+        assert dh.cpu_arrays[src1.name][i, -1, mirrored_dirN[5]] == dh.cpu_arrays[src1.name][i, -2, 5]
+
+    assert dh.cpu_arrays[src1.name][4, -1, mirrored_dirN[1]] == dh.cpu_arrays[src1.name][4, -2, 1]
+    assert dh.cpu_arrays[src1.name][4, -1, mirrored_dirN[5]] == dh.cpu_arrays[src1.name][4, -2, 5]
+
+    # check East
+    assert dh.cpu_arrays[src1.name][-1, 1, mirrored_dirE[6]] == dh.cpu_arrays[src1.name][-2, 1, 6]
+    assert dh.cpu_arrays[src1.name][-1, 1, mirrored_dirE[4]] == dh.cpu_arrays[src1.name][-2, 1, 4]
+
+    for i in range(2, 4):
+        assert dh.cpu_arrays[src1.name][-1, i, mirrored_dirE[6]] == dh.cpu_arrays[src1.name][-2, i, 6]
+        assert dh.cpu_arrays[src1.name][-1, i, mirrored_dirE[4]] == dh.cpu_arrays[src1.name][-2, i, 4]
+        assert dh.cpu_arrays[src1.name][-1, i, mirrored_dirE[8]] == dh.cpu_arrays[src1.name][-2, i, 8]
+
+    assert dh.cpu_arrays[src1.name][-1, 4, mirrored_dirE[4]] == dh.cpu_arrays[src1.name][-2, 4, 4]
+    assert dh.cpu_arrays[src1.name][-1, 4, mirrored_dirE[8]] == dh.cpu_arrays[src1.name][-2, 4, 8]
+
+    # check South
+    assert dh.cpu_arrays[src1.name][1, 0, mirrored_dirS[8]] == dh.cpu_arrays[src1.name][1, 1, 8]
+    assert dh.cpu_arrays[src1.name][1, 0, mirrored_dirS[2]] == dh.cpu_arrays[src1.name][1, 1, 2]
+
+    for i in range(2, 4):
+        assert dh.cpu_arrays[src1.name][i, 0, mirrored_dirS[7]] == dh.cpu_arrays[src1.name][i, 1, 7]
+        assert dh.cpu_arrays[src1.name][i, 0, mirrored_dirS[2]] == dh.cpu_arrays[src1.name][i, 1, 2]
+        assert dh.cpu_arrays[src1.name][i, 0, mirrored_dirS[8]] == dh.cpu_arrays[src1.name][i, 1, 8]
+
+    assert dh.cpu_arrays[src1.name][4, 0, mirrored_dirS[2]] == dh.cpu_arrays[src1.name][4, 1, 2]
+    assert dh.cpu_arrays[src1.name][4, 0, mirrored_dirS[7]] == dh.cpu_arrays[src1.name][4, 1, 7]
+
+    # check West
+    assert dh.cpu_arrays[src1.name][0, 1, mirrored_dirW[5]] == dh.cpu_arrays[src1.name][1, 1, 5]
+    assert dh.cpu_arrays[src1.name][0, 1, mirrored_dirW[3]] == dh.cpu_arrays[src1.name][1, 1, 3]
+
+    for i in range(2, 4):
+        assert dh.cpu_arrays[src1.name][0, i, mirrored_dirW[5]] == dh.cpu_arrays[src1.name][1, i, 5]
+        assert dh.cpu_arrays[src1.name][0, i, mirrored_dirW[3]] == dh.cpu_arrays[src1.name][1, i, 3]
+        assert dh.cpu_arrays[src1.name][0, i, mirrored_dirW[7]] == dh.cpu_arrays[src1.name][1, i, 7]
+
+    assert dh.cpu_arrays[src1.name][0, 4, mirrored_dirW[3]] == dh.cpu_arrays[src1.name][1, 4, 3]
+    assert dh.cpu_arrays[src1.name][0, 4, mirrored_dirW[7]] == dh.cpu_arrays[src1.name][1, 4, 7]
+
+    if given_normal:
+        # check corners --> determined by the last boundary applied there.
+        # SouthWest --> West
+        assert dh.cpu_arrays[src1.name][0, 0, mirrored_dirW[5]] == dh.cpu_arrays[src1.name][1, 0, 5]
+        # NorthWest --> West
+        assert dh.cpu_arrays[src1.name][0, -1, mirrored_dirW[7]] == dh.cpu_arrays[src1.name][1, -1, 7]
+        # NorthEast --> East
+        assert dh.cpu_arrays[src1.name][-1, -1, mirrored_dirE[8]] == dh.cpu_arrays[src1.name][-2, -1, 8]
+        # SouthEast --> East
+        assert dh.cpu_arrays[src1.name][-1, 0, mirrored_dirE[6]] == dh.cpu_arrays[src1.name][-2, 0, 6]
+    else:
+        # check corners --> this time the normals are calculated correctly in the corners
+        # SouthWest --> Normal = (1, 1); dir 7 --> 6
+        assert dh.cpu_arrays[src1.name][0, 0, 6] == dh.cpu_arrays[src1.name][1, 1, 7]
+        # NorthWest --> Normal = (1, -1); dir 8 --> 5
+        assert dh.cpu_arrays[src1.name][0, -1, 8] == dh.cpu_arrays[src1.name][1, -2, 5]
+        # NorthEast --> Normal = (-1, -1); dir 7 --> 6
+        assert dh.cpu_arrays[src1.name][-1, -1, 7] == dh.cpu_arrays[src1.name][-2, -2, 6]
+        # SouthEast --> Normal = (-1, 1); dir 5 --> 8
+        assert dh.cpu_arrays[src1.name][-1, 0, 5] == dh.cpu_arrays[src1.name][-2, 1, 8]
+
+
 def test_free_slip_index_list():
     stencil = LBStencil(Stencil.D2Q9)
     dh = create_data_handling(domain_size=(4, 4), periodicity=(False, False))
@@ -181,6 +297,50 @@ def test_free_slip_index_list():
             assert normal == normal_north_east
 
 
+def test_free_slip_index_list_convex_corner():
+    stencil = LBStencil(Stencil.D2Q9)
+    dh = create_data_handling(domain_size=(4, 4))
+    src = dh.add_array('src', values_per_cell=len(stencil))
+    dh.fill('src', 0.0, ghost_layers=True)
+
+    lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.8)
+    method = create_lb_method(lbm_config=lbm_config)
+
+    def bh_callback(x, y):
+        radius = 2
+        x_mid = 2
+        y_mid = 2
+        return (x - x_mid) ** 2 + (y - y_mid) ** 2 > radius ** 2
+
+    bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name="bh")
+
+    free_slip = FreeSlip(stencil=stencil)
+    bh.set_boundary(free_slip, mask_callback=bh_callback)
+
+    bh.prepare()
+    for b in dh.iterate():
+        for b_obj, idx_arr in b[bh._index_array_name].boundary_object_to_index_list.items():
+            index_array = idx_arr
+
+    # correct index array for this case with convex corners
+    test = [(2, 1, 2, 0, 1, 2), (2, 1, 3, 1, 0, 3), (2, 1, 7, 1, 1, 7),
+            (2, 1, 8, 0, 1, 7), (3, 1, 2, 0, 1, 2), (3, 1, 4, -1, 0, 4),
+            (3, 1, 7, 0, 1, 8), (3, 1, 8, -1, 1, 8), (1, 2, 2, 0, 1, 2),
+            (1, 2, 3, 1, 0, 3), (1, 2, 5, 1, 0, 7), (1, 2, 7, 1, 1, 7),
+            (2, 2, 7, 1, 1, 7), (3, 2, 8, -1, 1, 8), (4, 2, 2, 0, 1, 2),
+            (4, 2, 4, -1, 0, 4), (4, 2, 6, -1, 0, 8), (4, 2, 8, -1, 1, 8),
+            (1, 3, 1, 0, -1, 1), (1, 3, 3, 1, 0, 3), (1, 3, 5, 1, -1, 5),
+            (1, 3, 7, 1, 0, 5), (2, 3, 5, 1, -1, 5), (3, 3, 6, -1, -1, 6),
+            (4, 3, 1, 0, -1, 1), (4, 3, 4, -1, 0, 4), (4, 3, 6, -1, -1, 6),
+            (4, 3, 8, -1, 0, 6), (2, 4, 1, 0, -1, 1), (2, 4, 3, 1, 0, 3),
+            (2, 4, 5, 1, -1, 5), (2, 4, 6, 0, -1, 5), (3, 4, 1, 0, -1, 1),
+            (3, 4, 4, -1, 0, 4), (3, 4, 5, 0, -1, 6), (3, 4, 6, -1, -1, 6)]
+
+    for i, cell in enumerate(index_array):
+        for j in range(len(cell)):
+            assert cell[j] == test[i][j]
+
+
 def test_free_slip_equivalence():
     # check if Free slip BC does the same if the normal direction is specified or not
 
@@ -214,7 +374,7 @@ def test_free_slip_equivalence():
     bh1()
     bh2()
 
-    assert np.array_equal(dh.cpu_arrays['src1'], dh.cpu_arrays['src2'])
+    assert np.array_equal(dh.gather_array('src1'), dh.gather_array('src2'))
 
 
 def test_exotic_boundaries():