From ffcf6991cbbb1c7862229e9000990cf3b7be1bc9 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Thu, 15 Feb 2018 18:55:43 +0100
Subject: [PATCH] Boundary handling for Finite Differences
- splitted existing LBM boundary handling into two parts:
-> generic part , that is used for FD as well and moved it to pystencils
-> LBM specific part - remained in lbmpy
- bugfixes
---
boundaries/__init__.py | 1 +
boundaries/boundaryconditions.py | 68 ++++++
boundaries/boundaryhandling.py | 328 +++++++++++++++++++++++++++
boundaries/createindexlist.py | 79 +++++++
boundaries/createindexlistcython.pyx | 63 +++++
kernelcreation.py | 2 +-
sympyextensions.py | 47 ++--
7 files changed, 573 insertions(+), 15 deletions(-)
create mode 100644 boundaries/__init__.py
create mode 100644 boundaries/boundaryconditions.py
create mode 100644 boundaries/boundaryhandling.py
create mode 100644 boundaries/createindexlist.py
create mode 100644 boundaries/createindexlistcython.pyx
diff --git a/boundaries/__init__.py b/boundaries/__init__.py
new file mode 100644
index 000000000..a7ab44b75
--- /dev/null
+++ b/boundaries/__init__.py
@@ -0,0 +1 @@
+from pystencils.boundaries.boundaryhandling import BoundaryHandling
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
new file mode 100644
index 000000000..80906feb7
--- /dev/null
+++ b/boundaries/boundaryconditions.py
@@ -0,0 +1,68 @@
+import sympy as sp
+from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
+
+
+class Boundary(object):
+ """Base class all boundaries should derive from"""
+
+ def __init__(self, name=None):
+ self._name = name
+
+ def __call__(self, field, directionSymbol, indexField):
+ """
+ This function defines the boundary behavior and must therefore be implemented by all boundaries.
+ Here the boundary is defined as a list of sympy equations, from which a boundary kernel is generated.
+ :param field: pystencils field where boundary condition should be applied.
+ The current cell is cell next to the boundary, which is influenced by the boundary
+ cell i.e. has a link from the boundary cell to itself.
+ :param directionSymbol: a sympy symbol that can be used as index to the pdfField. It describes
+ the direction pointing from the fluid to the boundary cell
+ :param indexField: the boundary index field that can be used to retrieve and update boundary data
+ :return: list of sympy equations
+ """
+ raise NotImplementedError("Boundary class has to overwrite __call__")
+
+ @property
+ def additionalData(self):
+ """Return a list of (name, type) tuples for additional data items required in this boundary
+ These data items can either be initialized in separate kernel see additionalDataKernelInit or by
+ Python callbacks - see additionalDataCallback """
+ return []
+
+ @property
+ def additionalDataInitCallback(self):
+ """Return a callback function called with a boundary data setter object and returning a dict of
+ data-name to data for each element that should be initialized"""
+ return None
+
+ @property
+ def name(self):
+ if self._name:
+ return self._name
+ else:
+ return type(self).__name__
+
+ @name.setter
+ def name(self, newValue):
+ self._name = newValue
+
+
+class Neumann(Boundary):
+ def __call__(self, field, directionSymbol, **kwargs):
+
+ neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, field.spatialDimensions)
+ if field.indexDimensions == 0:
+ return [sp.Eq(field[neighbor], field.center)]
+ else:
+ from itertools import product
+ if not field.hasFixedIndexShape:
+ raise NotImplementedError("Neumann boundary works only for fields with fixed index shape")
+ indexIter = product(*(range(i) for i in field.indexShape))
+ return [sp.Eq(field[neighbor](idx), field(idx)) for idx in indexIter]
+
+ def __hash__(self):
+ # All boundaries of these class behave equal -> should also be equal
+ return hash("Neumann")
+
+ def __eq__(self, other):
+ return type(other) == Neumann
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
new file mode 100644
index 000000000..9591acad7
--- /dev/null
+++ b/boundaries/boundaryhandling.py
@@ -0,0 +1,328 @@
+import numpy as np
+import sympy as sp
+from pystencils import Field, TypedSymbol, createIndexedKernel
+from pystencils.backends.cbackend import CustomCppCode
+from pystencils.boundaries.createindexlist import numpyDataTypeForBoundaryObject, createBoundaryIndexArray
+from pystencils.cache import memorycache
+from pystencils.data_types import createType
+
+
+class BoundaryHandling:
+
+ def __init__(self, dataHandling, fieldName, stencil, name="boundaryHandling", target='cpu', openMP=True):
+ assert dataHandling.hasData(fieldName)
+
+ self._dataHandling = dataHandling
+ self._fieldName = fieldName
+ self._flagFieldName = name + "Flags"
+ self._indexArrayName = name + "IndexArrays"
+ self._target = target
+ self._openMP = openMP
+ self._boundaryObjectToBoundaryInfo = {}
+ self._fluidFlag = 1 << 0
+ self._nextFreeFlag = 1
+ self.stencil = stencil
+ self._dirty = True
+
+ # Add flag field to data handling if it does not yet exist
+ if dataHandling.hasData(self._flagFieldName) or dataHandling.hasData(self._indexArrayName):
+ raise ValueError("There is already a boundary handling registered at the data handling."
+ "If you want to add multiple handlings, choose a different name.")
+
+ gpu = self._target == 'gpu'
+ dataHandling.addArray(self._flagFieldName, dtype=np.uint32, cpu=True, gpu=False)
+ dataHandling.addCustomClass(self._indexArrayName, self.IndexFieldBlockData, cpu=True, gpu=gpu)
+
+ ffGhostLayers = self._dataHandling.ghostLayersOfField(self._flagFieldName)
+ for b in self._dataHandling.iterate(ghostLayers=ffGhostLayers):
+ b[self._flagFieldName].fill(self._fluidFlag)
+
+ @property
+ def dataHandling(self):
+ return self._dataHandling
+
+ @property
+ def shape(self):
+ return self._dataHandling.shape
+
+ @property
+ def dim(self):
+ return self._dataHandling.dim
+
+ @property
+ def boundaryObjects(self):
+ return tuple(self._boundaryObjectToName.keys())
+
+ @property
+ def flagArrayName(self):
+ return self._flagFieldName
+
+ def getBoundaryNameToFlagDict(self):
+ result = {bObj.name: bInfo.flag for bObj, bInfo in self._boundaryObjectToBoundaryInfo.items()}
+ result['fluid'] = self._fluidFlag
+ return result
+
+ def getMask(self, sliceObj, boundaryObj, inverse=False):
+ if isinstance(boundaryObj, str) and boundaryObj.lower() == 'fluid':
+ flag = self._fluidFlag
+ else:
+ flag = self._boundaryObjectToBoundaryInfo[boundaryObj].flag
+
+ arr = self.dataHandling.gatherArray(self.flagArrayName, sliceObj)
+ if arr is None:
+ return None
+ else:
+ result = np.bitwise_and(arr, flag)
+ if inverse:
+ result = np.logical_not(result)
+ return result
+
+ def setBoundary(self, boundaryObject, sliceObj=None, maskCallback=None, ghostLayers=True, innerGhostLayers=True):
+ """
+ Sets boundary using either a rectangular slice, a boolean mask or a combination of both
+
+ :param boundaryObject: instance of a boundary object that should be set
+ :param sliceObj: a slice object (can be created with makeSlice[]) that selects a part of the domain where
+ the boundary should be set. If none, the complete domain is selected which makes only sense
+ if a maskCallback is passed. The slice can have ':' placeholders, which are interpreted
+ depending on the 'includeGhostLayers' parameter i.e. if it is True, the slice extends
+ into the ghost layers
+ :param maskCallback: callback function getting x,y (z) parameters of the cell midpoints and returning a
+ boolean mask with True entries where boundary cells should be set.
+ The x, y, z arrays have 2D/3D shape such that they can be used directly
+ to create the boolean return array. i.e return x < 10 sets boundaries in cells with
+ midpoint x coordinate smaller than 10.
+ :param ghostLayers see DataHandling.iterate()
+ """
+ if isinstance(boundaryObject, str) and boundaryObject.lower() == 'fluid':
+ flag = self._fluidFlag
+ else:
+ flag = self._getFlagForBoundary(boundaryObject)
+
+ for b in self._dataHandling.iterate(sliceObj, ghostLayers=ghostLayers, innerGhostLayers=innerGhostLayers):
+ flagArr = b[self._flagFieldName]
+ if maskCallback is not None:
+ mask = maskCallback(*b.midpointArrays)
+ flagArr[mask] = flag
+ else:
+ flagArr.fill(flag)
+
+ self._dirty = True
+
+ def prepare(self):
+ if not self._dirty:
+ return
+ self._createIndexFields()
+ self._dirty = False
+
+ def triggerReinitializationOfBoundaryData(self, **kwargs):
+ if self._dirty:
+ self.prepare()
+ else:
+ ffGhostLayers = self._dataHandling.ghostLayersOfField(self._flagFieldName)
+ for b in self._dataHandling.iterate(ghostLayers=ffGhostLayers):
+ for bObj, setter in b[self._indexArrayName].boundaryObjectToDataSetter.items():
+ self._boundaryDataInitialization(bObj, setter, **kwargs)
+
+ def __call__(self, **kwargs):
+ if self._dirty:
+ self.prepare()
+
+ for b in self._dataHandling.iterate(gpu=self._target == 'gpu'):
+ for bObj, idxArr in b[self._indexArrayName].boundaryObjectToIndexList.items():
+ kwargs[self._fieldName] = b[self._fieldName]
+ self._boundaryObjectToBoundaryInfo[bObj].kernel(indexField=idxArr, **kwargs)
+
+ def geometryToVTK(self, fileName='geometry', boundaries='all', ghostLayers=False):
+ """
+ Writes a VTK field where each cell with the given boundary is marked with 1, other cells are 0
+ This can be used to display the simulation geometry in Paraview
+ :param fileName: vtk filename
+ :param boundaries: boundary object, or special string 'fluid' for fluid cells or special string 'all' for all
+ boundary conditions.
+ can also be a sequence, to write multiple boundaries to VTK file
+ :param ghostLayers: number of ghost layers to write, or True for all, False for none
+ """
+ if boundaries == 'all':
+ boundaries = list(self._boundaryObjectToBoundaryInfo.keys()) + ['fluid']
+ elif not hasattr(boundaries, "__len__"):
+ boundaries = [boundaries]
+
+ masksToName = {}
+ for b in boundaries:
+ if b == 'fluid':
+ masksToName[self._fluidFlag] = 'fluid'
+ else:
+ masksToName[self._boundaryObjectToBoundaryInfo[b].flag] = b.name
+
+ writer = self.dataHandling.vtkWriterFlags(fileName, self._flagFieldName, masksToName, ghostLayers=ghostLayers)
+ writer(1)
+
+ # ------------------------------ Implementation Details ------------------------------------------------------------
+
+ def _getFlagForBoundary(self, boundaryObject):
+ if boundaryObject not in self._boundaryObjectToBoundaryInfo:
+ symbolicIndexField = Field.createGeneric('indexField', spatialDimensions=1,
+ dtype=numpyDataTypeForBoundaryObject(boundaryObject, self.dim))
+ ast = self._createBoundaryKernel(self._dataHandling.fields[self._fieldName],
+ symbolicIndexField, boundaryObject)
+ boundaryInfo = self.BoundaryInfo(boundaryObject, flag=1 << self._nextFreeFlag, kernel=ast.compile())
+
+ self._nextFreeFlag += 1
+ self._boundaryObjectToBoundaryInfo[boundaryObject] = boundaryInfo
+ return self._boundaryObjectToBoundaryInfo[boundaryObject].flag
+
+ def _createBoundaryKernel(self, symbolicField, symbolicIndexField, boundaryObject):
+ return createBoundaryKernel(symbolicField, symbolicIndexField, self.stencil, boundaryObject,
+ target=self._target, openMP=self._openMP)
+
+ def _createIndexFields(self):
+ dh = self._dataHandling
+ ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
+ for b in dh.iterate(ghostLayers=ffGhostLayers):
+ flagArr = b[self._flagFieldName]
+ pdfArr = b[self._fieldName]
+ indexArrayBD = b[self._indexArrayName]
+ indexArrayBD.clear()
+ for bInfo in self._boundaryObjectToBoundaryInfo.values():
+ idxArr = createBoundaryIndexArray(flagArr, self.stencil, bInfo.flag, self._fluidFlag,
+ bInfo.boundaryObject, dh.ghostLayersOfField(self._flagFieldName))
+ if idxArr.size == 0:
+ continue
+
+ boundaryDataSetter = BoundaryDataSetter(idxArr, b.offset, self.stencil, ffGhostLayers, pdfArr)
+ indexArrayBD.boundaryObjectToIndexList[bInfo.boundaryObject] = idxArr
+ indexArrayBD.boundaryObjectToDataSetter[bInfo.boundaryObject] = boundaryDataSetter
+ self._boundaryDataInitialization(bInfo.boundaryObject, boundaryDataSetter)
+
+ def _boundaryDataInitialization(self, boundaryObject, boundaryDataSetter, **kwargs):
+ if boundaryObject.additionalDataInitCallback:
+ boundaryObject.additionalDataInitCallback(boundaryDataSetter, **kwargs)
+ if self._target == 'gpu':
+ self._dataHandling.toGpu(self._indexArrayName)
+
+ class BoundaryInfo(object):
+ def __init__(self, boundaryObject, flag, kernel):
+ self.boundaryObject = boundaryObject
+ self.flag = flag
+ self.kernel = kernel
+
+ class IndexFieldBlockData:
+ def __init__(self, *args, **kwargs):
+ self.boundaryObjectToIndexList = {}
+ self.boundaryObjectToDataSetter = {}
+
+ def clear(self):
+ self.boundaryObjectToIndexList.clear()
+ self.boundaryObjectToDataSetter.clear()
+
+ @staticmethod
+ def toCpu(gpuVersion, cpuVersion):
+ gpuVersion = gpuVersion.boundaryObjectToIndexList
+ cpuVersion = cpuVersion.boundaryObjectToIndexList
+ for obj, cpuArr in cpuVersion.values():
+ gpuVersion[obj].get(cpuArr)
+
+ @staticmethod
+ def toGpu(gpuVersion, cpuVersion):
+ from pycuda import gpuarray
+ gpuVersion = gpuVersion.boundaryObjectToIndexList
+ cpuVersion = cpuVersion.boundaryObjectToIndexList
+ for obj, cpuArr in cpuVersion.items():
+ if obj not in gpuVersion:
+ gpuVersion[obj] = gpuarray.to_gpu(cpuArr)
+ else:
+ gpuVersion[obj].set(cpuArr)
+
+
+class BoundaryDataSetter:
+
+ def __init__(self, indexArray, offset, stencil, ghostLayers, pdfArray):
+ self.indexArray = indexArray
+ self.offset = offset
+ self.stencil = np.array(stencil)
+ self.pdfArray = pdfArray.view()
+ self.pdfArray.flags.writeable = False
+
+ arrFieldNames = indexArray.dtype.names
+ self.dim = 3 if 'z' in arrFieldNames else 2
+ assert 'x' in arrFieldNames and 'y' in arrFieldNames and 'dir' in arrFieldNames, str(arrFieldNames)
+ self.boundaryDataNames = set(self.indexArray.dtype.names) - set(['x', 'y', 'z', 'dir'])
+ self.coordMap = {0: 'x', 1: 'y', 2: 'z'}
+ self.ghostLayers = ghostLayers
+
+ def fluidCellPositions(self, coord):
+ assert coord < self.dim
+ return self.indexArray[self.coordMap[coord]] + self.offset[coord] - self.ghostLayers + 0.5
+
+ @memorycache()
+ def linkOffsets(self):
+ return self.stencil[self.indexArray['dir']]
+
+ @memorycache()
+ def linkPositions(self, coord):
+ return self.fluidCellPositions(coord) + 0.5 * self.linkOffsets()[:, coord]
+
+ @memorycache()
+ def boundaryCellPositions(self, coord):
+ return self.fluidCellPositions(coord) + self.linkOffsets()[:, coord]
+
+ def __setitem__(self, key, value):
+ if key not in self.boundaryDataNames:
+ raise KeyError("Invalid boundary data name %s. Allowed are %s" % (key, self.boundaryDataNames))
+ self.indexArray[key] = value
+
+ def __getitem__(self, item):
+ if item not in self.boundaryDataNames:
+ raise KeyError("Invalid boundary data name %s. Allowed are %s" % (item, self.boundaryDataNames))
+ return self.indexArray[item]
+
+
+class BoundaryOffsetInfo(CustomCppCode):
+
+ # --------------------------- Functions to be used by boundaries --------------------------
+
+ @staticmethod
+ def offsetFromDir(dirIdx, dim):
+ return tuple([sp.IndexedBase(symbol, shape=(1,))[dirIdx]
+ for symbol in BoundaryOffsetInfo._offsetSymbols(dim)])
+
+ @staticmethod
+ def invDir(dirIdx):
+ return sp.IndexedBase(BoundaryOffsetInfo.INV_DIR_SYMBOL, shape=(1,))[dirIdx]
+
+ # ---------------------------------- Internal ---------------------------------------------
+
+ def __init__(self, stencil):
+ dim = len(stencil[0])
+
+ offsetSym = BoundaryOffsetInfo._offsetSymbols(dim)
+ code = "\n"
+ for i in range(dim):
+ offsetStr = ", ".join([str(d[i]) for d in stencil])
+ code += "const int64_t %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" % (self.INV_DIR_SYMBOL.name, ", ".join(invDirs))
+ offsetSymbols = BoundaryOffsetInfo._offsetSymbols(dim)
+ super(BoundaryOffsetInfo, self).__init__(code, symbolsRead=set(),
+ symbolsDefined=set(offsetSymbols + [self.INV_DIR_SYMBOL]))
+
+ @staticmethod
+ def _offsetSymbols(dim):
+ return [TypedSymbol("c_%d" % (d,), createType(np.int64)) for d in range(dim)]
+
+ INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
+
+
+def createBoundaryKernel(field, indexField, stencil, boundaryFunctor, target='cpu', openMP=True):
+ elements = [BoundaryOffsetInfo(stencil)]
+ indexArrDtype = indexField.dtype.numpyDtype
+ dirSymbol = TypedSymbol("dir", indexArrDtype.fields['dir'][0])
+ elements += [sp.Eq(dirSymbol, indexField[0]('dir'))]
+ elements += boundaryFunctor(field, directionSymbol=dirSymbol, indexField=indexField)
+ return createIndexedKernel(elements, [indexField], target=target, cpuOpenMP=openMP)
diff --git a/boundaries/createindexlist.py b/boundaries/createindexlist.py
new file mode 100644
index 000000000..ee46b0b42
--- /dev/null
+++ b/boundaries/createindexlist.py
@@ -0,0 +1,79 @@
+import numpy as np
+import itertools
+import warnings
+
+try:
+ import pyximport;
+
+ pyximport.install()
+ from pystencils.boundaries.createindexlistcython import createBoundaryIndexList2D, createBoundaryIndexList3D
+
+ cythonFuncsAvailable = True
+except Exception:
+ cythonFuncsAvailable = False
+ createBoundaryIndexList2D = None
+ 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 = boundaryIndexArrayCoordinateNames[:len(flagFieldArr.shape)]
+ indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [(directionMemberName, np.int32)])
+
+ result = []
+ gl = nrOfGhostLayers
+ for cell in itertools.product(*reversed([range(gl, i-gl) for i in flagFieldArr.shape])):
+ cell = cell[::-1]
+ if not flagFieldArr[cell] & fluidMask:
+ continue
+ for dirIdx, direction in enumerate(stencil):
+ neighborCell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
+ if flagFieldArr[neighborCell] & boundaryMask:
+ result.append(cell + (dirIdx,))
+
+ return np.array(result, dtype=indexArrDtype)
+
+
+def createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGhostLayers=1):
+ dim = len(flagField.shape)
+ 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)
+ if dim == 2:
+ idxList = createBoundaryIndexList2D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
+ elif dim == 3:
+ idxList = createBoundaryIndexList3D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
+ else:
+ raise ValueError("Flag field has to be a 2 or 3 dimensional numpy array")
+ return np.array(idxList, dtype=indexArrDtype)
+ else:
+ if flagField.size > 1e6:
+ warnings.warn("Boundary setup may take very long! Consider installing cython to speed it up")
+ return _createBoundaryIndexListPython(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
+
+
+def createBoundaryIndexArray(flagField, stencil, boundaryMask, fluidMask, boundaryObject, nrOfGhostLayers=1):
+ idxArray = createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGhostLayers)
+ dim = len(flagField.shape)
+
+ if boundaryObject.additionalData:
+ coordinateNames = boundaryIndexArrayCoordinateNames[:dim]
+ indexArrDtype = numpyDataTypeForBoundaryObject(boundaryObject, dim)
+ extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype)
+ for prop in coordinateNames + ['dir']:
+ extendedIdxField[prop] = idxArray[prop]
+
+ idxArray = extendedIdxField
+
+ return idxArray
diff --git a/boundaries/createindexlistcython.pyx b/boundaries/createindexlistcython.pyx
new file mode 100644
index 000000000..dfa84069d
--- /dev/null
+++ b/boundaries/createindexlistcython.pyx
@@ -0,0 +1,63 @@
+# Workaround for cython bug
+# see https://stackoverflow.com/questions/8024805/cython-compiled-c-extension-importerror-dynamic-module-does-not-define-init-fu
+WORKAROUND = "Something"
+
+import cython
+
+ctypedef fused IntegerType:
+ short
+ int
+ long
+ long long
+ unsigned short
+ unsigned int
+ unsigned long
+
+@cython.boundscheck(False) # turn off bounds-checking for entire function
+@cython.wraparound(False) # turn off negative index wrapping for entire function
+def createBoundaryIndexList2D(object[IntegerType, ndim=2] flagField,
+ int nrOfGhostLayers, IntegerType boundaryMask, IntegerType fluidMask,
+ object[int, ndim=2] stencil):
+ cdef int xs, ys, x, y
+ cdef int dirIdx, numDirections, dx, dy
+
+ xs, ys = flagField.shape
+ boundaryIndexList = []
+ numDirections = stencil.shape[0]
+
+ for y in range(nrOfGhostLayers,ys-nrOfGhostLayers):
+ for x in range(nrOfGhostLayers,xs-nrOfGhostLayers):
+ if flagField[x,y] & fluidMask:
+ for dirIdx in range(1, numDirections):
+ dx = stencil[dirIdx,0]
+ dy = stencil[dirIdx,1]
+ if flagField[x+dx, y+dy] & boundaryMask:
+ boundaryIndexList.append((x,y, dirIdx))
+ return boundaryIndexList
+
+
+@cython.boundscheck(False) # turn off bounds-checking for entire function
+@cython.wraparound(False) # turn off negative index wrapping for entire function
+def createBoundaryIndexList3D(object[IntegerType, ndim=3] flagField,
+ int nrOfGhostLayers, IntegerType boundaryMask, IntegerType fluidMask,
+ object[int, ndim=2] stencil):
+ cdef int xs, ys, zs, x, y, z
+ cdef int dirIdx, numDirections, dx, dy, dz
+
+ xs, ys, zs = flagField.shape
+ boundaryIndexList = []
+ numDirections = stencil.shape[0]
+
+ for z in range(nrOfGhostLayers, zs-nrOfGhostLayers):
+ for y in range(nrOfGhostLayers,ys-nrOfGhostLayers):
+ for x in range(nrOfGhostLayers,xs-nrOfGhostLayers):
+ if flagField[x, y, z] & fluidMask:
+ for dirIdx in range(1, numDirections):
+ dx = stencil[dirIdx,0]
+ dy = stencil[dirIdx,1]
+ dz = stencil[dirIdx,2]
+ if flagField[x + dx, y + dy, z + dz] & boundaryMask:
+ boundaryIndexList.append((x,y,z, dirIdx))
+ return boundaryIndexList
+
+
diff --git a/kernelcreation.py b/kernelcreation.py
index 65760903e..7d27ed6c9 100644
--- a/kernelcreation.py
+++ b/kernelcreation.py
@@ -88,7 +88,7 @@ def createIndexedKernel(equations, indexFields, target='cpu', dataType="double",
if target == 'cpu':
from pystencils.cpu import createIndexedKernel
from pystencils.cpu import addOpenMP
- ast = createIndexedKernel(equations, indexField=indexFields, typeForSymbol=dataType,
+ ast = createIndexedKernel(equations, indexFields=indexFields, typeForSymbol=dataType,
coordinateNames=coordinateNames)
if cpuOpenMP:
addOpenMP(ast, numThreads=cpuOpenMP)
diff --git a/sympyextensions.py b/sympyextensions.py
index f56781e67..566cc7441 100644
--- a/sympyextensions.py
+++ b/sympyextensions.py
@@ -5,6 +5,8 @@ import itertools
import warnings
import sympy as sp
+from pystencils.data_types import getTypeOfExpression, getBaseType
+
def prod(seq):
"""Takes a sequence and returns the product of all elements"""
@@ -414,7 +416,7 @@ def mostCommonTermFactorization(term):
return sp.Mul(commonFactor, term, evaluate=False)
-def countNumberOfOperations(term):
+def countNumberOfOperations(term, onlyType='real'):
"""
Counts the number of additions, multiplications and division
:param term: a sympy term, equation or sequence of terms/equations
@@ -433,15 +435,31 @@ def countNumberOfOperations(term):
term = term.evalf()
+ def checkType(e):
+ if onlyType is None:
+ return True
+ try:
+ type = getBaseType(getTypeOfExpression(e))
+ except ValueError:
+ return False
+ if onlyType == 'int' and (type.is_int() or type.is_uint()):
+ return True
+ if onlyType == 'real' and (type.is_float()):
+ return True
+ else:
+ return type == onlyType
+
def visit(t):
visitChildren = True
if t.func is sp.Add:
- result['adds'] += len(t.args) - 1
+ if checkType(t):
+ result['adds'] += len(t.args) - 1
elif t.func is sp.Mul:
- result['muls'] += len(t.args) - 1
- for a in t.args:
- if a == 1 or a == -1:
- result['muls'] -= 1
+ if checkType(t):
+ result['muls'] += len(t.args) - 1
+ for a in t.args:
+ if a == 1 or a == -1:
+ result['muls'] -= 1
elif t.func is sp.Float:
pass
elif isinstance(t, sp.Symbol):
@@ -451,14 +469,15 @@ def countNumberOfOperations(term):
elif t.is_integer:
pass
elif t.func is sp.Pow:
- visitChildren = False
- if t.exp.is_integer and t.exp.is_number:
- if t.exp >= 0:
- result['muls'] += int(t.exp) - 1
- else:
- result['muls'] -= 1
- result['divs'] += 1
- result['muls'] += (-int(t.exp)) - 1
+ if checkType(t):
+ visitChildren = False
+ if t.exp.is_integer and t.exp.is_number:
+ if t.exp >= 0:
+ result['muls'] += int(t.exp) - 1
+ else:
+ result['muls'] -= 1
+ result['divs'] += 1
+ result['muls'] += (-int(t.exp)) - 1
else:
warnings.warn("Counting operations: only integer exponents are supported in Pow, "
"counting will be inaccurate")
--
GitLab