From 1d1a192ae08bc5de03e900ceee09295acf86b1c0 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 26 Jan 2018 12:16:49 +0100
Subject: [PATCH] New boundary handling
- plotting
- calculation of force on boundary
- vtk outpu
---
datahandling/__init__.py | 2 +-
datahandling/datahandling_interface.py | 41 +++++++++++-
datahandling/parallel_datahandling.py | 63 +++++++++++++++++--
datahandling/serial_datahandling.py | 86 +++++++++++++++++++++++---
parallel/blockiteration.py | 6 +-
vtk.py | 10 +--
6 files changed, 185 insertions(+), 23 deletions(-)
diff --git a/datahandling/__init__.py b/datahandling/__init__.py
index e30d0a5a3..719573055 100644
--- a/datahandling/__init__.py
+++ b/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,)
diff --git a/datahandling/datahandling_interface.py b/datahandling/datahandling_interface.py
index ced806dd8..af69c67c8 100644
--- a/datahandling/datahandling_interface.py
+++ b/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"""
diff --git a/datahandling/parallel_datahandling.py b/datahandling/parallel_datahandling.py
index b92e7963b..20fd0a99f 100644
--- a/datahandling/parallel_datahandling.py
+++ b/datahandling/parallel_datahandling.py
@@ -1,5 +1,5 @@
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
+
diff --git a/datahandling/serial_datahandling.py b/datahandling/serial_datahandling.py
index 3fd69756a..ed149b9b3 100644
--- a/datahandling/serial_datahandling.py
+++ b/datahandling/serial_datahandling.py
@@ -1,8 +1,8 @@
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)
+
diff --git a/parallel/blockiteration.py b/parallel/blockiteration.py
index 6e12791a8..e6280d216 100644
--- a/parallel/blockiteration.py
+++ b/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):
diff --git a/vtk.py b/vtk.py
index e47f22eda..2b5408180 100644
--- a/vtk.py
+++ b/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]]
--
GitLab