Fix FreeSlip boundary condition

There is an issue with the proxy field substitution and the FreeSlip boundary condition. This is true in particular if the normal vector is not pre-determined as implemented in walberla/walberla!520 (merged).

Currently, with the pull pattern the following substitutions are performed:

f_in[0,0](invdir[dir]) ← f_in[{indexField}_{0}^{wnx},{indexField}_{0}^{wny}](indexField[0](ref_dir))
->
pdfs[{f_{in inv offsets x}}_{dir},{f_{in inv offsets y}}_{dir},{f_{in inv offsets z}}_{dir}](f_in_inv_dir_idx[dir]) ← pdfs[{indexVector}_{0}^{wnx} + {f_{in offsets x}}_{{indexVector}_{0}^{ref_dir}},{indexVector}_{0}^{wny} + {f_{in offsets y}}_{{indexVector}_{0}^{ref_dir}},{indexVector}_{0}^{wnz} + {f_{in offsets z}}_{{indexVector}_{0}^{ref_dir}}](indexVector[0](ref_dir))

Whereas the correct substitutions are:

f_in[0,0](invdir[dir]) ← f_in[{indexField}_{0}^{wnx},{indexField}_{0}^{wny}](indexField[0](ref_dir))
->
pdfs[{f_{in inv offsets x}}_{dir},{f_{in inv offsets y}}_{dir},{f_{in inv offsets z}}_{dir}](f_in_inv_dir_idx[dir]) ← pdfs[{indexVector}_{0}^{wnx} + {f_{in inv offsets x}}_{dir},{indexVector}_{0}^{wny} + {f_{in inv offsets y}}_{dir},{indexVector}_{0}^{wnz} + {f_{in inv offsets z}}_{dir}](indexVector[0](ref_dir))

The problem is that the population index (ref_dir) is used to determine a spatial offset. This MR in its current form fixes this issue for the particular case of the FreeSlip condition and the pull pattern. Other configurations break.

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