lbmpy: geometry for pipe setup

- bugfixes in parallel walberla scenario

geometry.py

+import numpy as np
+from lbmpy.boundaries import NoSlip
+
+
+def addPipe(boundaryHandling, diameter, boundary=NoSlip()):
+    """
+    Sets boundary for the wall of a pipe with flow in x direction.
+
+    :param boundaryHandling: boundary handling object, works for serial and parallel versions 
+    :param diameter: pipe diameter, can be either a constant value or a callback function.
+                     the callback function has the signature (xCoordArray, domainShapeInCells) and has to return
+                     a array of same shape as the received xCoordArray, with the diameter for each x position
+    :param boundary: boundary object that is set at the wall, defaults to NoSlip (bounce back)
+    """
+    domainShape = boundaryHandling.domainShape
+    dim = len(domainShape)
+    assert dim in (2, 3)
+
+    def callback(*coordinates):
+        flowCoord = coordinates[0]
+
+        if callable(diameter):
+            flowCoordLine = flowCoord[:, 0, 0] if dim == 3 else flowCoord[:, 0]
+            diameterValue = diameter(flowCoordLine, domainShape)
+            diameterValue = diameterValue[:, np.newaxis, np.newaxis] if dim == 3 else diameterValue[:, np.newaxis]
+        else:
+            diameterValue = diameter
+
+        radiusSq = (diameterValue / 2) ** 2
+
+        mid = [domainShape[i] // 2 for i in range(1, dim)]
+        distance = sum((c_i - mid_i) ** 2 for c_i, mid_i in zip(coordinates[1:], mid))
+        return distance > radiusSq
+
+    boundaryHandling.setBoundary(boundary, maskCallback=callback)

parallel/boundaryhandling.py

-import numpy as np
 import waLBerla as wlb
 from lbmpy.boundaries.handlinginterface import GenericBoundaryHandling, FlagFieldInterface
 from lbmpy.parallel.blockiteration import slicedBlockIteration
-from pystencils.slicing import normalizeSlice
 
 
 class BoundaryHandling(object):
@@ -19,6 +17,7 @@ class BoundaryHandling(object):
         self._pdfFieldId = pdfFieldId
         self._blocks = blocks
         self.dim = lbMethod.dim
+        self.domainShape = [i + 1 for i in self._blocks.getDomainCellBB().max]
         blocks.addBlockData(boundaryId, addBoundaryHandling)
 
     def __call__(self, **kwargs):
@@ -36,7 +35,7 @@ class BoundaryHandling(object):
 
     def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, sliceNormalizationGhostLayers=1):
         if indexExpr is None:
-            indexExpr = [slice(None, None, None)] * self.dim
+            indexExpr = [slice(None, None, None)] * 3
 
         if len(self._blocks) == 0:
             return
@@ -50,29 +49,6 @@ class BoundaryHandling(object):
             mask = maskCallback(*indexArrs) if maskCallback else None
             block[self._boundaryId].setBoundary(boundaryObject, indexExpr=localSlice, maskArr=mask)
 
-        #domainCellBB = self._blocks.getDomainCellBB()
-        #domainExtent = [s + 2 * sliceNormalizationGhostLayers for s in domainCellBB.size]
-        #indexExpr = normalizeSlice(indexExpr, domainExtent)
-        #targetCellBB = wlb.CellInterval.fromSlice(indexExpr)
-        #targetCellBB.shift(*[a - sliceNormalizationGhostLayers for a in domainCellBB.min])
-        #
-        #for block in self._blocks:
-        #    boundaryHandling = block[self._boundaryId]
-        #    ghostLayers = boundaryHandling.ghostLayers
-        #    intersection = self._blocks.getBlockCellBB(block).getExpanded(ghostLayers)
-        #    intersection.intersect(targetCellBB)
-        #    if not intersection.empty():
-        #        if maskCallback is not None:
-        #            meshGridParams = [offset + 0.5 + np.arange(width)
-        #                              for offset, width in zip(intersection.min, intersection.size)]
-        #            indexArr = np.meshgrid(*meshGridParams, indexing='ij')
-        #            mask = maskCallback(*indexArr)
-        #        else:
-        #            mask = None
-        #        localTargetBB = self._blocks.transformGlobalToLocal(block, intersection)
-        #        localTargetBB.shift(ghostLayers, ghostLayers, ghostLayers)
-        #        block[self._boundaryId].setBoundary(boundaryObject, indexExpr=localTargetBB.toSlice(), maskArr=mask)
-
 
 # ----------------------------------------------------------------------------------------------------------------------
 
@@ -115,7 +91,7 @@ class WalberlaFlagFieldInterface(FlagFieldInterface):
         raise NotImplementedError()
 
 
-if __name__ == '__main__':
+def test():
     from lbmpy.creationfunctions import createLatticeBoltzmannMethod
     from lbmpy.boundaries.boundaryconditions import NoSlip
     from pystencils.slicing import makeSlice
@@ -129,7 +105,6 @@ if __name__ == '__main__':
 
     bh = BoundaryHandling(blocks, lbMethod, 'pdfField', 'flagField')
 
-
     def maskCallback(x, y, z):
         midPoint = (4, 4, 4)
         radius = 2.0
@@ -139,7 +114,6 @@ if __name__ == '__main__':
                      (z - midPoint[2]) ** 2 <= radius ** 2
         return resultMask
 
-
     bh.setBoundary(NoSlip(), makeSlice[:, :, :], maskCallback, sliceNormalizationGhostLayers=0)
     vtkOutput = wlb.vtk.makeOutput(blocks, "vtkOutput")
     flagFieldWriter = wlb.field.createVTKWriter(blocks, "flagField")

parallel/scenarios.py

@@ -9,6 +9,7 @@ from lbmpy.parallel.blockiteration import slicedBlockIteration
 
 
 class Scenario(object):
+    vtkCounter = 0
 
     def __init__(self, blocks, methodParameters, optimizationParams, lbmKernel=None,
                  preUpdateFunctions=[], kernelParams={}, blockDataPrefix='', directCommunication=False):
@@ -17,6 +18,10 @@ class Scenario(object):
         self._blockDataPrefix = blockDataPrefix
         self.kernelParams = kernelParams
 
+        if 'stencil' not in methodParameters:
+            domainShape = [i + 1 for i in blocks.getDomainCellBB().max]
+            methodParameters['stencil'] = 'D2Q9' if domainShape[2] == 1 else 'D3Q27'
+
         methodParameters, optimizationParams = updateWithDefaultParameters(methodParameters, optimizationParams)
         self._target = optimizationParams['target']
         if optimizationParams['fieldLayout'] in ('f', 'reverseNumpy'):
@@ -36,7 +41,7 @@ class Scenario(object):
             switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams)
             methodParameters['optimizationParams'] = optimizationParams
             self._lbmKernel = createLatticeBoltzmannFunction(**methodParameters)
-        else:for
+        else:
             self._lbmKernel = lbmKernel
 
         # ----- Add fields
@@ -78,7 +83,8 @@ class Scenario(object):
         self._macroscopicValueSetter = None
         self._macroscopicValueGetter = None
 
-        self._vtkOutput = wlb.vtk.makeOutput(blocks, "vtk")
+        self._vtkOutput = wlb.vtk.makeOutput(blocks, "vtk_%02d" % (Scenario.vtkCounter,))
+        Scenario.vtkCounter += 1
         self._vtkOutput.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, self.flagFieldId))
         self._vtkOutput.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, self.densityFieldId))
         self._vtkOutput.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, self.velocityFieldId))