lbmpy: improved boundary setup + testcase for pressure driven channel

boundaries/boundaryhandling.py

@@ -30,6 +30,9 @@ class BoundaryHandling(object):
         if name is None:
             name = boundaryFunction.__name__
 
+        if name in self._nameToIndex:
+            return 2 ** self._nameToIndex[name]
+
         newIdx = len(self._boundaryFunctions)
         self._nameToIndex[name] = newIdx
         self._boundaryFunctions.append(boundaryFunction)
@@ -45,7 +48,7 @@ class BoundaryHandling(object):
     def getFlag(self, name):
         return 2 ** self._nameToIndex[name]
 
-    def setBoundary(self, function, indexExpr, clearOtherBoundaries=True):
+    def setBoundary(self, function, indexExpr, maskArr=None):
         if hasattr(function, '__name__'):
             name = function.__name__
         elif hasattr(function, 'name'):
@@ -57,13 +60,12 @@ class BoundaryHandling(object):
             self.addBoundary(function, name)
 
         flag = self.getFlag(name)
-        if clearOtherBoundaries:
+
+        if maskArr is None:
             self.flagField[indexExpr] = flag
         else:
-            # clear fluid flag
-            np.bitwise_and(self.flagField[indexExpr], np.invert(self._fluidFlag), self.flagField[indexExpr])
-            # add new boundary flag
-            np.bitwise_or(self.flagField[indexExpr], flag, self.flagField[indexExpr])
+            flagFieldView = self.flagField[indexExpr]
+            flagFieldView[maskArr] = flag
 
         self.invalidateIndexCache()
 

boundaries/geometry.py

+import numpy as np
+import scipy.misc
+from pystencils.slicing import makeSlice, normalizeSlice, shiftSlice
+
+
+class BlackAndWhiteImageBoundary:
+
+    def __init__(self, imagePath, boundaryFunction, targetSlice=makeSlice[:, :, :]):
+        self.imgArr = scipy.misc.imread(imagePath, flatten=True).astype(int)
+        self.imgArr = np.rot90(self.imgArr, 3)
+
+        self._boundaryFunction = boundaryFunction
+        self._targetSlice = targetSlice
+
+    def __call__(self, boundaryHandling, method, domainSize, **kwargs):
+        normalizedSlice = normalizeSlice(self._targetSlice, domainSize)
+        normalizedSlice = shiftSlice(normalizedSlice, 1)
+        targetSize = [s.stop - s.start for s in normalizedSlice]
+        img = scipy.misc.imresize(self.imgArr, size=targetSize)
+        img[img <= 254] = 0
+        img[img > 254] = 1
+        boundaryHandling.setBoundary(self._boundaryFunction, normalizedSlice, maskArr=(img == 0))
+

serialscenario.py

@@ -69,6 +69,9 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
 
             pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0]
         getMacroscopic(pdfs=pdfArrays[0], density=densityArr[0], velocity=velocityArr[0])
+        #for vComp in range(velocityArr[0].shape[-1]):
+        #    v = velocityArr[0][..., vComp]
+        #    v[boundaryHandling.flagField != boundaryHandling._fluidFlag] = 0
         return pdfArrays[0], densityArr[0], velocityArr[0]
 
     def gpuTimeLoop(timeSteps):
@@ -112,7 +115,8 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={},
 
 
 def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None,
-                                        lbmKernel=None, optimizationParams={}, kernelParams={}, **kwargs):
+                                        lbmKernel=None, optimizationParams={}, boundarySetupFunctions=[],
+                                        kernelParams={}, **kwargs):
     assert dim in (2, 3)
 
     if radius is not None:
@@ -150,6 +154,9 @@ def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None
                 for direction in ('N', 'S', 'T', 'B'):
                     boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
 
+        for userFunction in boundarySetupFunctions:
+            userFunction(boundaryHandling, method, domainSize)
+
     assert domainSize is not None
     return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel,
                           kernelParams=kernelParams)
@@ -188,7 +195,7 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
                 for direction in ('N', 'S', 'T', 'B'):
                     boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
         for userFunction in boundarySetupFunctions:
-            userFunction(boundaryHandling, method)
+            userFunction(boundaryHandling, method, domainSize)
 
     def periodicity(pdfArr):
         pdfArr[0, :, :] = pdfArr[-2, :, :]