Simple functions to compute force on boundary

- added functions to compute force on boundary for generic boundary and
  NoSlip
- simple calculations for serial scenarios
- to distinguish between different obstacles, NoSlip boundaries can now
  be named

boundaries/boundaryconditions.py

@@ -51,6 +51,11 @@ class Boundary(object):
 
 
 class NoSlip(Boundary):
+
+    def __init__(self, name=None):
+        """Set an optional name here, to mark boundaries, for example for force evaluations"""
+        self._name = name
+
     """No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
     def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
         neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
@@ -58,11 +63,20 @@ class NoSlip(Boundary):
         return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(directionSymbol))]
 
     def __hash__(self):
-        # All boundaries of these class behave equal -> should also be equal
-        return hash("NoSlip")
+        # All boundaries of these class behave equal -> should also be equal (as long as name is equal)
+        return hash(self.name)
 
     def __eq__(self, other):
-        return type(other) == NoSlip
+        if not isinstance(other, NoSlip):
+            return False
+        return self.name == other.name
+
+    @property
+    def name(self):
+        if self._name:
+            return self._name
+        else:
+            return type(self).__name__
 
 
 class NoSlipFullWay(Boundary):

boundaries/forceevaluation.py

+import numpy as np
+from lbmpy.stencils import inverseDirection
+
+
+def calculateForceOnNoSlipBoundary(boundaryObject, boundaryHandling, pdfArray):
+    indArr = boundaryHandling.getBoundaryIndexArray(boundaryObject)
+    method = boundaryHandling.lbMethod
+
+    stencil = np.array(method.stencil)
+
+    if method.dim == 2:
+        x, y = indArr['x'], indArr['y']
+        values = 2 * pdfArray[x, y,  indArr['dir']]
+    else:
+        assert method.dim == 3
+        x, y, z = indArr['x'], indArr['y'], indArr['z']
+        values = 2 * pdfArray[x, y, z, indArr['dir']]
+
+    forces = stencil[indArr['dir']] * values[:, np.newaxis]
+    return forces.sum(axis=0)
+
+
+def calculateForceOnBoundary(boundaryObject, boundaryHandling, pdfArray):
+    indArr = boundaryHandling.getBoundaryIndexArray(boundaryObject)
+    method = boundaryHandling.lbMethod
+
+    stencil = np.array(method.stencil)
+    invDirection = np.array([method.stencil.index(inverseDirection(d))
+                            for d in method.stencil])
+
+    if method.dim == 2:
+        x, y = indArr['x'], indArr['y']
+        offsets = stencil[indArr['dir']]
+
+        fluidValues = pdfArray[x, y,  indArr['dir']]
+        boundaryValues = pdfArray[x + offsets[:, 0],
+                                  y + offsets[:, 1],
+                                  invDirection[indArr['dir']]]
+    else:
+        assert method.dim == 3
+        x, y, z = indArr['x'], indArr['y'], indArr['z']
+        offsets = stencil[indArr['dir']]
+
+        fluidValues = pdfArray[x, y, z, indArr['dir']]
+        boundaryValues = pdfArray[x + offsets[:, 0],
+                                  y + offsets[:, 1],
+                                  z + offsets[:, 2],
+                                  invDirection[indArr['dir']]]
+
+    values = fluidValues + boundaryValues
+    forces = stencil[indArr['dir']] * values[:, np.newaxis]
+    return forces.sum(axis=0)

boundaries/handlinginterface.py

@@ -113,6 +113,10 @@ class GenericBoundaryHandling(object):
         if target not in ('cpu', 'gpu'):
             raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
 
+    @property
+    def lbMethod(self):
+        return self._lbMethod
+
     def getBoundaryNameToFlagDict(self):
         result = {self._flagFieldInterface.getName(o): self._flagFieldInterface.getFlag(o) for o in self._boundaryInfos}
         result['fluid'] = self._flagFieldInterface.getFlag('fluid')

scenarios.py

@@ -26,6 +26,7 @@ that defines the viscosity of the fluid (valid values being between 0 and 2).
 import numpy as np
 import sympy as sp
 
+from lbmpy.boundaries.forceevaluation import calculateForceOnBoundary, calculateForceOnNoSlipBoundary
 from lbmpy.geometry import setupChannelWalls, addParabolicVelocityInflow
 from pystencils.field import getLayoutOfArray, createNumpyArrayWithLayout
 from pystencils.slicing import sliceFromDirection, addGhostLayers, removeGhostLayers, normalizeSlice, makeSlice
@@ -346,6 +347,14 @@ class Scenario(object):
         """Equation collection defining the LBM update rule (already in simplified form)"""
         return self._lbmKernel.updateRule
 
+    def calculateForceOnBoundary(self, boundaryObject):
+        """Computes force on boundary using simple momentum exchange method"""
+        if isinstance(boundaryObject, NoSlip):
+            return calculateForceOnNoSlipBoundary(boundaryObject, self.boundaryHandling, self._pdfArrays[0])
+        else:
+            self.runBoundaryHandlingOnly()
+            return calculateForceOnBoundary(boundaryObject, self.boundaryHandling, self._pdfArrays[0])
+
     def animateVelocity(self, steps=10, **kwargs):
         import lbmpy.plot2d as plt
 
@@ -415,3 +424,11 @@ class Scenario(object):
         self._getMacroscopic(pdfs=self._pdfArrays[0], density=self._density, velocity=self._velocity,
                              **self.kernelParams)
 
+    def runBoundaryHandlingOnly(self):
+        isGpuSimulation = len(self._pdfGpuArrays) > 0
+        if isGpuSimulation:
+            self._pdfGpuArrays[0].set(self._pdfArrays[0])
+            self._boundaryHandling(pdfs=self._pdfGpuArrays[0], **self.kernelParams)
+            self._pdfGpuArrays[0].get(self._pdfArrays[0])
+        else:
+            self._boundaryHandling(pdfs=self._pdfArrays[0], **self.kernelParams)