diff --git a/boundaries/boundary_kernel.py b/boundaries/boundary_kernel.py index a7e86cc0768d713a64e2e930ce2f86671ba47afd..ed6f40075b25fa094353123b0801c51f93e3d7ef 100644 --- a/boundaries/boundary_kernel.py +++ b/boundaries/boundary_kernel.py @@ -1,6 +1,8 @@ import sympy as sp +import numpy as np from pystencils.backends.cbackend import CustomCppCode from pystencils import TypedSymbol, Field +from pystencils.data_types import castFunc, createType INV_DIR_SYMBOL = TypedSymbol("invDir", "int") WEIGHTS_SYMBOL = TypedSymbol("weights", "double") @@ -10,7 +12,7 @@ WEIGHTS_SYMBOL = TypedSymbol("weights", "double") def offsetSymbols(dim): - return [TypedSymbol("c_%d" % (d,), "int") for d in range(dim)] + return [TypedSymbol("c_%d" % (d,), createType(np.int64)) for d in range(dim)] def offsetFromDir(dirIdx, dim): @@ -36,7 +38,7 @@ class LbmMethodInfo(CustomCppCode): code = "\n" for i in range(lbMethod.dim): offsetStr = ", ".join([str(d[i]) for d in stencil]) - code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr) + code += "const int64_t %s [] = { %s };\n" % (offsetSym[i].name, offsetStr) invDirs = [] for direction in stencil: @@ -52,9 +54,15 @@ class LbmMethodInfo(CustomCppCode): def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu', createInitializationKernel=False): indexField = Field.createFromNumpyArray("indexField", indexArr) + return generateIndexBoundaryKernelGeneric(pdfField, indexField, indexArr.dtype, lbMethod, boundaryFunctor, target, + createInitializationKernel) + + +def generateIndexBoundaryKernelGeneric(pdfField, indexField, indexArrDtype, lbMethod, boundaryFunctor, target='cpu', + createInitializationKernel=False): elements = [LbmMethodInfo(lbMethod)] - dirSymbol = TypedSymbol("dir", indexArr.dtype.fields['dir'][0]) + dirSymbol = TypedSymbol("dir", indexArrDtype.fields['dir'][0]) boundaryEqList = [sp.Eq(dirSymbol, indexField[0]('dir'))] if createInitializationKernel: boundaryEqList += boundaryFunctor.additionalDataInitKernelEquations(pdfField=pdfField, directionSymbol=dirSymbol, @@ -69,5 +77,4 @@ def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, t return createIndexedKernel(elements, [indexField]) elif target == 'gpu': from pystencils.gpucuda import createdIndexedCUDAKernel - return createdIndexedCUDAKernel(elements, [indexField]) - + return createdIndexedCUDAKernel(elements, [indexField]) \ No newline at end of file diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py index c49e613927d13204189246619d012790e80aa6c3..f2902c8481185baaae12cf46ac31db9d9bbc712c 100644 --- a/boundaries/boundaryconditions.py +++ b/boundaries/boundaryconditions.py @@ -1,5 +1,5 @@ import sympy as sp - +import numpy as np from lbmpy.simplificationfactory import createSimplificationStrategy from pystencils.astnodes import SympyAssignment from pystencils.sympyextensions import getSymmetricPart @@ -11,6 +11,9 @@ from pystencils.data_types import createTypeFromString class Boundary(object): """Base class all boundaries should derive from""" + def __init__(self, name=None): + self._name = name + def __call__(self, pdfField, directionSymbol, lbMethod, indexField): """ This function defines the boundary behavior and must therefore be implemented by all boundaries. @@ -47,14 +50,21 @@ class Boundary(object): @property def name(self): - return type(self).__name__ + if self._name: + return self._name + else: + return type(self).__name__ + + @name.setter + def name(self, newValue): + self._name = newValue class NoSlip(Boundary): def __init__(self, name=None): """Set an optional name here, to mark boundaries, for example for force evaluations""" - self._name = name + super(NoSlip, self).__init__(name) """No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle""" def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs): @@ -71,13 +81,6 @@ class NoSlip(Boundary): 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): """Full-way bounce back""" @@ -110,7 +113,7 @@ class UBB(Boundary): """Velocity bounce back boundary condition, enforcing specified velocity at obstacle""" - def __init__(self, velocity, adaptVelocityToForce=False, dim=None): + def __init__(self, velocity, adaptVelocityToForce=False, dim=None, name=None): """ :param velocity: can either be a constant, an access into a field, or a callback function. @@ -118,6 +121,7 @@ class UBB(Boundary): and 'velocity' which has to be set to the desired velocity of the corresponding link :param adaptVelocityToForce: """ + super(UBB, self).__init__(name) self._velocity = velocity self._adaptVelocityToForce = adaptVelocityToForce if callable(self._velocity) and not dim: @@ -180,7 +184,8 @@ class UBB(Boundary): class FixedDensity(Boundary): - def __init__(self, density): + def __init__(self, density, name=None): + super(FixedDensity, self).__init__(name) self._density = density def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs): diff --git a/boundaries/boundaryhandling_old.py b/boundaries/boundaryhandling_old.py new file mode 100644 index 0000000000000000000000000000000000000000..3be5b7089d4a717c6780d905929a6d824669e576 --- /dev/null +++ b/boundaries/boundaryhandling_old.py @@ -0,0 +1,280 @@ +import sympy as sp +import numpy as np + +from lbmpy.stencils import getStencil +from pystencils import TypedSymbol, Field +from pystencils.backends.cbackend import CustomCppCode +from lbmpy.boundaries.createindexlist import createBoundaryIndexList +from pystencils.slicing import normalizeSlice, makeSlice + +INV_DIR_SYMBOL = TypedSymbol("invDir", "int") +WEIGHTS_SYMBOL = TypedSymbol("weights", "double") + + +class BoundaryHandling(object): + class BoundaryInfo(object): + def __init__(self, flag, object, kernel, ast): + self.flag = flag + self.object = object + self.kernel = kernel + self.ast = ast + + def __init__(self, pdfField, domainShape, lbMethod, ghostLayers=1, target='cpu', openMP=True): + """ + Class for managing boundary kernels + + :param pdfField: pdf numpy array including ghost layers + :param domainShape: domain size without ghost layers + :param lbMethod: lattice Boltzmann method + :param ghostLayers: number of ghost layers + :param target: either 'cpu' or 'gpu' + """ + symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1) + assert pdfField.shape[:-1] == tuple(d + 2*ghostLayers for d in domainShape) + + self._symbolicPdfField = symbolicPdfField + self._pdfField = pdfField + self._shapeWithGhostLayers = [d + 2 * ghostLayers for d in domainShape] + self._fluidFlag = 2 ** 0 + self.flagField = np.full(self._shapeWithGhostLayers, self._fluidFlag, dtype=np.int32) + self._ghostLayers = ghostLayers + self._lbMethod = lbMethod + + self._boundaryInfos = {} + self._periodicityKernels = [] + + self._dirty = False + self._periodicity = [False, False, False] + self._target = target + self.openMP = openMP + if target not in ('cpu', 'gpu'): + raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,)) + + @property + def periodicity(self): + """List that indicates for x,y (z) coordinate if domain is periodic in that direction""" + return self._periodicity + + @property + def fluidFlag(self): + """Flag that is set where the lattice Boltzmann update should happen""" + return self._fluidFlag + + def getFlag(self, boundaryObject): + """Flag that represents the given boundary.""" + return self._boundaryInfos[boundaryObject].flag + + def getBoundaries(self): + return [b.object for b in self._boundaryInfos.values()] + + def setPeriodicity(self, x=False, y=False, z=False): + """Enable periodic boundary conditions at the border of the domain""" + for d in (x, y, z): + assert isinstance(d, bool) + + self._periodicity = [x, y, z] + self._compilePeriodicityKernels() + + def hasBoundary(self, boundaryObject): + """Returns boolean indicating if a boundary with that name exists""" + return boundaryObject in self._boundaryInfos + + def setBoundary(self, boundaryObject, indexExpr=None, maskArr=None): + """ + Sets boundary in a rectangular region (slice) + + :param boundaryObject: boundary condition object or the string 'fluid' to remove boundary conditions + :param indexExpr: slice expression, where boundary should be set, see :mod:`pystencils.slicing` + :param maskArr: optional boolean (masked) array specifying where the boundary should be set + """ + if indexExpr is None: + indexExpr = [slice(None, None, None)] * len(self.flagField.shape) + if boundaryObject == 'fluid': + flag = self._fluidFlag + else: + flag = self.addBoundary(boundaryObject) + assert flag != self._fluidFlag + + indexExpr = normalizeSlice(indexExpr, self._shapeWithGhostLayers) + if maskArr is None: + self.flagField[indexExpr] = flag + else: + flagFieldView = self.flagField[indexExpr] + flagFieldView[maskArr] = flag + self._dirty = True + + def addBoundary(self, boundaryObject): + """ + Adds a boundary condition, i.e. reserves a flog in the flag field and returns that flag + If the boundary already exists, the existing flag is returned. + This flag can be logicalled or'ed to the boundaryHandling.flagField + + :param boundaryObject: boundary condition object, see :mod:`lbmpy.boundaries.boundaryconditions` + """ + if boundaryObject in self._boundaryInfos: + return self._boundaryInfos[boundaryObject].flag + + newIdx = len(self._boundaryInfos) + 1 # +1 because 2**0 is reserved for fluid flag + boundaryInfo = self.BoundaryInfo(2 ** newIdx, boundaryObject, None, None) + self._boundaryInfos[boundaryObject] = boundaryInfo + self._dirty = True + return boundaryInfo.flag + + def indices(self, dx=1.0, includeGhostLayers=False): + if not includeGhostLayers: + params = [np.arange(start=-1, stop=s-1) * dx for s in self.flagField.shape] + else: + params = [np.arange(s) * dx for s in self.flagField.shape] + return np.meshgrid(*params, indexing='ij') + + def __call__(self, **kwargs): + """Run the boundary handling, all keyword args are passed through to the boundary sweeps""" + if self._dirty: + self.prepare() + for boundary in self._boundaryInfos.values(): + boundary.kernel(**kwargs) + for k in self._periodicityKernels: + k(**kwargs) + + def clear(self): + """Removes all boundaries and fills the domain with fluid""" + self.flagField.fill(self._fluidFlag) + self._dirty = False + self._boundaryInfos = {} + + def prepare(self): + """Fills data structures to speed up the boundary handling and compiles all boundary kernels. + This is automatically done when first called. With this function this can be triggered before running.""" + for boundary in self._boundaryInfos.values(): + assert boundary.flag != self._fluidFlag + idxArray = createBoundaryIndexList(self.flagField, self._lbMethod.stencil, + boundary.flag, self._fluidFlag, self._ghostLayers) + + dim = self._lbMethod.dim + + if boundary.object.additionalData: + coordinateNames = ["x", "y", "z"][:dim] + indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + + [('dir', np.int32)] + + [(i[0], i[1].numpyDtype) for i in boundary.object.additionalData]) + extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype) + for prop in coordinateNames + ['dir']: + extendedIdxField[prop] = idxArray[prop] + + idxArray = extendedIdxField + if boundary.object.additionalDataInitKernelEquations: + initKernelAst = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod, + boundary.object, target='cpu', + createInitializationKernel=True) + from pystencils.cpu import makePythonFunction as makePythonCpuFunction + initKernel = makePythonCpuFunction(initKernelAst, {'indexField': idxArray, 'pdfs': self._pdfField}) + initKernel() + + ast = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod, boundary.object, + target=self._target) + boundary.ast = ast + if self._target == 'cpu': + from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP + addOpenMP(ast, numThreads=self.openMP) + boundary.kernel = makePythonCpuFunction(ast, {'indexField': idxArray}) + elif self._target == 'gpu': + from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction + import pycuda.gpuarray as gpuarray + idxGpuField = gpuarray.to_gpu(idxArray) + boundary.kernel = makePythonGpuFunction(ast, {'indexField': idxGpuField}) + else: + assert False + self._dirty = False + + def invalidateIndexCache(self): + """Invalidates the cache for optimization data structures. When setting boundaries the cache is automatically + invalidated, so there is no need to call this function manually, as long as the flag field is not manually + modified.""" + self._dirty = True + + def _compilePeriodicityKernels(self): + self._periodicityKernels = [] + dim = len(self.flagField.shape) + if dim == 2: + stencil = getStencil("D2Q9") + elif dim == 3: + stencil = getStencil("D3Q27") + else: + assert False + + filteredStencil = [] + for direction in stencil: + useDirection = True + if direction == (0, 0) or direction == (0, 0, 0): + useDirection = False + for component, periodicity in zip(direction, self._periodicity): + if not periodicity and component != 0: + useDirection = False + if useDirection: + filteredStencil.append(direction) + + if len(filteredStencil) > 0: + if self._target == 'cpu': + from pystencils.slicing import getPeriodicBoundaryFunctor + self._periodicityKernels.append(getPeriodicBoundaryFunctor(filteredStencil, ghostLayers=1)) + elif self._target == 'gpu': + from pystencils.gpucuda.periodicity import getPeriodicBoundaryFunctor + self._periodicityKernels.append(getPeriodicBoundaryFunctor(filteredStencil, self.flagField.shape, + indexDimensions=1, + indexDimShape=len(self._lbMethod.stencil), + ghostLayers=1)) + else: + assert False + + +# -------------------------------------- Helper Functions -------------------------------------------------------------- + + + + +# ------------------------------------- Kernel Generation -------------------------------------------------------------- + +class LbmMethodInfo(CustomCppCode): + def __init__(self, lbMethod): + stencil = lbMethod.stencil + symbolsDefined = set(offsetSymbols(lbMethod.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL]) + + offsetSym = offsetSymbols(lbMethod.dim) + code = "\n" + for i in range(lbMethod.dim): + offsetStr = ", ".join([str(d[i]) for d in stencil]) + code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr) + + invDirs = [] + for direction in stencil: + inverseDir = tuple([-i for i in direction]) + invDirs.append(str(stencil.index(inverseDir))) + + code += "const int %s [] = { %s };\n" % (INV_DIR_SYMBOL.name, ", ".join(invDirs)) + weights = [str(w.evalf()) for w in lbMethod.weights] + code += "const double %s [] = { %s };\n" % (WEIGHTS_SYMBOL.name, ",".join(weights)) + super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined) + + +def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu', + createInitializationKernel=False): + indexField = Field.createFromNumpyArray("indexField", indexArr) + + elements = [LbmMethodInfo(lbMethod)] + dirSymbol = TypedSymbol("dir", indexArr.dtype.fields['dir'][0]) + boundaryEqList = [sp.Eq(dirSymbol, indexField[0]('dir'))] + if createInitializationKernel: + boundaryEqList += boundaryFunctor.additionalDataInitKernelEquations(pdfField=pdfField, directionSymbol=dirSymbol, + lbMethod=lbMethod, indexField=indexField) + else: + boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod, + indexField=indexField) + elements += boundaryEqList + + if target == 'cpu': + from pystencils.cpu import createIndexedKernel + return createIndexedKernel(elements, [indexField]) + elif target == 'gpu': + from pystencils.gpucuda import createdIndexedCUDAKernel + return createdIndexedCUDAKernel(elements, [indexField]) + diff --git a/boundaries/createindexlist.py b/boundaries/createindexlist.py index a98baf9a8acfa5a17bcab82f06301b2093b43f27..69c5cb4dd9895c90660d764175a97c4f0492b15f 100644 --- a/boundaries/createindexlist.py +++ b/boundaries/createindexlist.py @@ -13,9 +13,20 @@ except Exception: createBoundaryIndexList3D = None +boundaryIndexArrayCoordinateNames = ["x", "y", "z"] +directionMemberName = "dir" + + +def numpyDataTypeForBoundaryObject(boundaryObject, dim): + coordinateNames = boundaryIndexArrayCoordinateNames[:dim] + return np.dtype([(name, np.int32) for name in coordinateNames] + + [(directionMemberName, np.int32)] + + [(i[0], i[1].numpyDtype) for i in boundaryObject.additionalData], align=True) + + def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, fluidMask, stencil): - coordinateNames = ("x", "y", "z")[:len(flagFieldArr.shape)] - indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [('dir', np.int32)]) + coordinateNames = boundaryIndexArrayCoordinateNames[:len(flagFieldArr.shape)] + indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [(directionMemberName, np.int32)]) result = [] gl = nrOfGhostLayers @@ -33,8 +44,8 @@ def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, def createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGhostLayers=1): dim = len(flagField.shape) - coordinateNames = ("x", "y", "z")[:dim] - indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [('dir', np.int32)]) + coordinateNames = boundaryIndexArrayCoordinateNames[:dim] + indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [(directionMemberName, np.int32)]) if cythonFuncsAvailable: stencil = np.array(stencil, dtype=np.int32) diff --git a/boundaries/handlinginterface.py b/boundaries/handlinginterface.py index 213824646ef93280d63ea1709bfa0e3938cba44d..0100cd98d4854175acacaa73812d803356384089 100644 --- a/boundaries/handlinginterface.py +++ b/boundaries/handlinginterface.py @@ -3,7 +3,8 @@ import numpy as np from pystencils.field import Field from pystencils.slicing import normalizeSlice, shiftSlice from lbmpy.boundaries.boundary_kernel import generateIndexBoundaryKernel -from lbmpy.boundaries.createindexlist import createBoundaryIndexList +from lbmpy.boundaries.createindexlist import createBoundaryIndexList, boundaryIndexArrayCoordinateNames,\ + numpyDataTypeForBoundaryObject from pystencils.cache import memorycache @@ -177,10 +178,8 @@ class GenericBoundaryHandling(object): continue if boundaryObject.additionalData: - coordinateNames = ["x", "y", "z"][:dim] - indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + - [('dir', np.int32)] + - [(i[0], i[1].numpyDtype) for i in boundaryObject.additionalData], align=True) + coordinateNames = boundaryIndexArrayCoordinateNames[:dim] + indexArrDtype = numpyDataTypeForBoundaryObject(boundaryObject, dim) extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype) for prop in coordinateNames + ['dir']: extendedIdxField[prop] = idxArray[prop]