New boundary handling

- plotting
- calculation of force on boundary
- vtk outpu

datahandling/__init__.py

@@ -21,7 +21,7 @@ def createDataHandling(parallel, domainSize, periodicity, defaultLayout='SoA', d
         elif periodicity is True:
             periodicity = (1, 1, 1)
         else:
-            periodicity = (int(bool(x)) for x in periodicity)
+            periodicity = tuple(int(bool(x)) for x in periodicity)
             if len(periodicity) == 2:
                 periodicity += (1,)
 

datahandling/datahandling_interface.py

@@ -19,6 +19,16 @@ class DataHandling(ABC):
     def dim(self):
         """Dimension of the domain, either 2 or 3"""
 
+    @property
+    @abstractmethod
+    def shape(self):
+        """Shape of outer bounding box"""
+
+    @property
+    @abstractmethod
+    def periodicity(self):
+        """Returns tuple of booleans for x,y,(z) directions with True if domain is periodic in that direction"""
+
     @abstractmethod
     def addArray(self, name, fSize=1, dtype=np.float64, latexName=None, ghostLayers=None, layout=None, cpu=True, gpu=False):
         """
@@ -86,7 +96,7 @@ class DataHandling(ABC):
         """Returns the number of ghost layers for a specific field/array"""
 
     @abstractmethod
-    def iterate(self, sliceObj=None, gpu=False, ghostLayers=None):
+    def iterate(self, sliceObj=None, gpu=False, ghostLayers=None, innerGhostLayers=True):
         """
         Iterate over local part of potentially distributed data structure.
         """
@@ -135,6 +145,26 @@ class DataHandling(ABC):
         """Copies data from CPU to GPU for all arrays that have a CPU and a GPU representation"""
         pass
 
+
+    @abstractmethod
+    def vtkWriter(self, fileName, dataNames, ghostLayers=False):
+        """VTK output for one or multiple arrays
+        :param fileName: base file name without extension for the VTK output
+        :param dataNames: list of array names that should be included in the vtk output
+        :param ghostLayers: true if ghost layer information should be written out as well
+        :return: a function that can be called with an integer time step to write the current state
+                i.e vtkWriter('someFile', ['velocity', 'density']) (1)
+        """
+    @abstractmethod
+    def vtkWriterFlags(self, fileName, dataName, masksToName, ghostLayers=False):
+        """VTK output for an unsigned integer field, where bits are intepreted as flags
+        :param fileName: see vtkWriter
+        :param dataName: name of an array with uint type
+        :param masksToName: dictionary mapping integer masks to a name in the output
+        :param ghostLayers: see vtkWriter
+        :returns: functor that can be called with time step
+         """
+
     # ------------------------------- Communication --------------------------------------------------------------------
 
     def synchronizationFunctionCPU(self, names, stencil=None, **kwargs):
@@ -152,3 +182,12 @@ class DataHandling(ABC):
         """
         Synchronization of GPU fields, for documentation see CPU version above
         """
+
+    def reduceFloatSequence(self, sequence, operation, allReduce=False):
+        """Takes a sequence of floating point values on each process and reduces it element wise to all
+        processes (allReduce=True) or only to the root process (allReduce=False).
+        Possible operations are 'sum', 'min', 'max'
+        """
+
+    def reduceIntSequence(self, sequence, operation, allReduce=False):
+        """See function reduceFloatSequence - this is the same for integers"""

datahandling/parallel_datahandling.py

 import numpy as np
-from pystencils import Field, makeSlice
+from pystencils import Field
 from pystencils.datahandling.datahandling_interface import DataHandling
 from pystencils.parallel.blockiteration import slicedBlockIteration, blockIteration
 from pystencils.utils import DotDict
@@ -33,6 +33,13 @@ class ParallelDataHandling(DataHandling):
         self._fieldInformation = {}
         self._cpuGpuPairs = []
         self._customDataTransferFunctions = {}
+
+        self._reduceMap = {
+            'sum': wlb.mpi.SUM,
+            'min': wlb.mpi.MIN,
+            'max': wlb.mpi.MAX,
+        }
+
         if self._dim == 2:
             assert self.blocks.getDomainCellBB().size[2] == 1
 
@@ -40,6 +47,14 @@ class ParallelDataHandling(DataHandling):
     def dim(self):
         return self._dim
 
+    @property
+    def shape(self):
+        return self.blocks.getDomainCellBB().size[:self.dim]
+
+    @property
+    def periodicity(self):
+        return self.blocks.periodic[:self._dim]
+
     @property
     def fields(self):
         return self._fields
@@ -120,22 +135,26 @@ class ParallelDataHandling(DataHandling):
         for block in self.blocks:
             block[name1].swapDataPointers(block[name2])
 
-    def iterate(self, sliceObj=None, gpu=False, ghostLayers=True):
+    def iterate(self, sliceObj=None, gpu=False, ghostLayers=True, innerGhostLayers=True):
         if ghostLayers is True:
             ghostLayers = self.defaultGhostLayers
         elif ghostLayers is False:
             ghostLayers = 0
+        if innerGhostLayers is True:
+            innerGhostLayers = self.defaultGhostLayers
+        elif innerGhostLayers is False:
+            innerGhostLayers = 0
 
         prefix = self.GPU_DATA_PREFIX if gpu else ""
         if sliceObj is not None:
-            yield from slicedBlockIteration(self.blocks, sliceObj, ghostLayers, ghostLayers,
+            yield from slicedBlockIteration(self.blocks, sliceObj, innerGhostLayers, ghostLayers,
                                             self.dim, prefix)
         else:
             yield from blockIteration(self.blocks, ghostLayers, self.dim, prefix)
 
     def gatherArray(self, name, sliceObj=None, allGather=False):
         if sliceObj is None:
-            sliceObj = makeSlice[:, :, :]
+            sliceObj = tuple([slice(None, None, None)] * self.dim)
         if self.dim == 2:
             sliceObj += (0.5,)
         for array in wlb.field.gatherGenerator(self.blocks, name, sliceObj, allGather):
@@ -223,3 +242,39 @@ class ParallelDataHandling(DataHandling):
             syncFunction.addDataToCommunicate(createPacking(self.blocks, name))
 
         return syncFunction
+
+    def reduceFloatSequence(self, sequence, operation, allReduce=False):
+        if allReduce:
+            return np.array(wlb.mpi.allreduceReal(sequence, self._reduceMap[operation.lower()]))
+        else:
+            return np.array(wlb.mpi.reduceReal(sequence, self._reduceMap[operation.lower()]))
+
+    def reduceIntSequence(self, sequence, operation, allReduce=False):
+        if allReduce:
+            return np.array(wlb.mpi.allreduceInt(sequence, self._reduceMap[operation.lower()]))
+        else:
+            return np.array(wlb.mpi.reduceInt(sequence, self._reduceMap[operation.lower()]))
+
+    def vtkWriter(self, fileName, dataNames, ghostLayers=False):
+        if ghostLayers is False:
+            ghostLayers = 0
+        if ghostLayers is True:
+            ghostLayers = min(self.ghostLayersOfField(n) for n in dataNames)
+
+        output = wlb.vtk.makeOutput(self.blocks, fileName, ghostLayers=ghostLayers)
+        for n in dataNames:
+            output.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, n))
+        return output
+
+    def vtkWriterFlags(self, fileName, dataName, masksToName, ghostLayers=False):
+        if ghostLayers is False:
+            ghostLayers = 0
+        if ghostLayers is True:
+            ghostLayers = self.ghostLayersOfField(dataName)
+
+        output = wlb.vtk.makeOutput(self.blocks, fileName, ghostLayers=ghostLayers)
+        for mask, name in masksToName.items():
+            w = wlb.field.createBinarizationVTKWriter(self.blocks, dataName, mask, name)
+            output.addCellDataWriter(w)
+        return output
+

datahandling/serial_datahandling.py

 import numpy as np
-from lbmpy.stencils import getStencil
+import itertools
 from pystencils import Field
 from pystencils.field import layoutStringToTuple, spatialLayoutStringToTuple, createNumpyArrayWithLayout
-from pystencils.parallel.blockiteration import Block, SerialBlock
+from pystencils.parallel.blockiteration import SerialBlock
 from pystencils.slicing import normalizeSlice, removeGhostLayers
 from pystencils.utils import DotDict
 from pystencils.datahandling.datahandling_interface import DataHandling
@@ -55,6 +55,14 @@ class SerialDataHandling(DataHandling):
     def dim(self):
         return len(self._domainSize)
 
+    @property
+    def shape(self):
+        return self._domainSize
+
+    @property
+    def periodicity(self):
+        return self._periodicity
+
     @property
     def fields(self):
         return self._fields
@@ -129,7 +137,7 @@ class SerialDataHandling(DataHandling):
     def addArrayLike(self, name, nameOfTemplateField, latexName=None, cpu=True, gpu=False):
         self.addArray(name, latexName=latexName, cpu=cpu, gpu=gpu, **self._fieldInformation[nameOfTemplateField])
 
-    def iterate(self, sliceObj=None, gpu=False, ghostLayers=True):
+    def iterate(self, sliceObj=None, gpu=False, ghostLayers=True, innerGhostLayers=True):
         if ghostLayers is True:
             ghostLayers = self.defaultGhostLayers
         elif ghostLayers is False:
@@ -208,19 +216,22 @@ class SerialDataHandling(DataHandling):
         return self._synchronizationFunctor(names, stencilName, 'gpu')
 
     def _synchronizationFunctor(self, names, stencil, target):
-        if stencil is None:
-            stencil = 'D3Q27' if self.dim == 3 else 'D2Q9'
-        if stencil == 'D3Q15' or stencil == 'D3Q19':
-            stencil = 'D3Q27'
-
-        assert stencil in ("D2Q9", 'D3Q27'), "Serial scenario support only D2Q9 or D3Q27 for periodicity sync"
 
         assert target in ('cpu', 'gpu')
         if not hasattr(names, '__len__') or type(names) is str:
             names = [names]
 
         filteredStencil = []
-        for direction in getStencil(stencil):
+        neighbors = [-1, 0, 1]
+
+        if stencil.startswith('D2'):
+            directions = itertools.product(*[neighbors] * 2)
+        elif stencil.startswith('D3'):
+            directions = itertools.product(*[neighbors] * 3)
+        else:
+            raise ValueError("Invalid stencil")
+
+        for direction in directions:
             useDirection = True
             if direction == (0, 0) or direction == (0, 0, 0):
                 useDirection = False
@@ -255,3 +266,58 @@ class SerialDataHandling(DataHandling):
                     func(pdfs=self.gpuArrays[name])
 
         return resultFunctor
+
+    @staticmethod
+    def reduceFloatSequence(sequence, operation, allReduce=False):
+        return np.array(sequence)
+
+    @staticmethod
+    def reduceIntSequence(sequence):
+        return np.array(sequence)
+
+    def vtkWriter(self, fileName, dataNames, ghostLayers=False):
+        from pystencils.vtk import imageToVTK
+
+        def writer(step):
+            fullFileName = "%s_%08d" % (fileName, step,)
+            cellData = {}
+            for name in dataNames:
+                field = self._getFieldWithGivenNumberOfGhostLayers(name, ghostLayers)
+                if self.dim == 2:
+                    field = field[:, :, np.newaxis]
+                if len(field.shape) == 3:
+                    field = np.ascontiguousarray(field)
+                elif len(field.shape) == 4:
+                    field = [np.ascontiguousarray(field[..., i]) for i in range(field.shape[-1])]
+                    if len(field) == 2:
+                        field.append(np.zeros_like(field[0]))
+                    field = tuple(field)
+                else:
+                    assert False
+                cellData[name] = field
+            imageToVTK(fullFileName, cellData=cellData)
+        return writer
+
+    def vtkWriterFlags(self, fileName, dataName, masksToName, ghostLayers=False):
+        from pystencils.vtk import imageToVTK
+
+        def writer(step):
+            fullFileName = "%s_%08d" % (fileName, step,)
+            field = self._getFieldWithGivenNumberOfGhostLayers(dataName, ghostLayers)
+            if self.dim == 2:
+                field = field[:, :, np.newaxis]
+            cellData = {name: np.ascontiguousarray(np.bitwise_and(field, mask) > 0, dtype=np.uint8)
+                        for mask, name in masksToName.items()}
+            imageToVTK(fullFileName, cellData=cellData)
+
+        return writer
+
+    def _getFieldWithGivenNumberOfGhostLayers(self, name, ghostLayers):
+        actualGhostLayers = self.ghostLayersOfField(name)
+        if ghostLayers is True:
+            ghostLayers = actualGhostLayers
+
+        glToRemove = actualGhostLayers - ghostLayers
+        indDims = 1 if self._fieldInformation[name]['fSize'] > 1 else 0
+        return removeGhostLayers(self.cpuArrays[name], indDims, glToRemove)
+

parallel/blockiteration.py

@@ -24,7 +24,7 @@ def blockIteration(blocks, ghostLayers, dim=3, accessPrefix=''):
         localSlice = [slice(0, w, None) for w in cellInterval.size]
         if dim == 2:
             localSlice[2] = ghostLayers
-        yield ParallelBlock(block, cellInterval.min, tuple(localSlice), ghostLayers, accessPrefix)
+        yield ParallelBlock(block, cellInterval.min[:dim], tuple(localSlice), ghostLayers, accessPrefix)
 
 
 def slicedBlockIteration(blocks, sliceObj=None, innerGhostLayers=1, outerGhostLayers=1, dim=3, accessPrefix=''):
@@ -69,7 +69,7 @@ def slicedBlockIteration(blocks, sliceObj=None, innerGhostLayers=1, outerGhostLa
         localSlice = localTargetBB.toSlice(False)
         if dim == 2:
             localSlice = (localSlice[0], localSlice[1], innerGhostLayers)
-        yield ParallelBlock(block, intersection.min, localSlice, innerGhostLayers, accessPrefix)
+        yield ParallelBlock(block, intersection.min[:dim], localSlice, innerGhostLayers, accessPrefix)
 
 
 class Block:
@@ -100,7 +100,7 @@ class Block:
     @property
     def shape(self):
         """Shape of the fields (potentially including ghost layers)"""
-        return tuple(s.stop - s.start for s in self._localSlice)
+        return tuple(s.stop - s.start for s in self._localSlice[:len(self._offset)])
 
     @property
     def globalSlice(self):

vtk.py

@@ -5,6 +5,7 @@ from pyevtk.hl import _addDataToFile, _appendDataToFile
 def imageToVTK(path, origin=(0.0, 0.0, 0.0), spacing=(1.0, 1.0, 1.0), cellData=None, pointData=None):
     """Patched version of same pyevtk function that also support vector data"""
     assert cellData != None or pointData != None
+
     # Extract dimensions
     start = (0, 0, 0)
     end = None
@@ -13,10 +14,11 @@ def imageToVTK(path, origin=(0.0, 0.0, 0.0), spacing=(1.0, 1.0, 1.0), cellData=N
         data = cellData[keys[0]]
         if hasattr(data, 'shape'):
             end = data.shape
-        elif data[0].ndim == 3 and data[1].ndim == 3 and data[0].ndim == 3:
-            keys = list(cellData.keys())
-            data = cellData[keys[0]]
-            end = data[0].shape
+        elif isinstance(data, tuple):
+            shapes = set(d.shape for d in data)
+            if len(shapes) > 1:
+                raise ValueError("All components have to have the same shape")
+            end = shapes.pop()
     elif pointData:
         keys = list(pointData.keys())
         data = pointData[keys[0]]