Commit ffcf6991 authored by Martin Bauer's avatar Martin Bauer
Browse files

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
parent a78c0879
from pystencils.boundaries.boundaryhandling import BoundaryHandling
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
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)
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
# 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
......@@ -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)
......
......@@ -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