diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index 6d1d37c07b7eef741c663b488ac61b2b6ea8e799..4315452897b7a53ba34818a409e2f78bb8e70274 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -1,2 +1,2 @@
from lbmpy.boundaries.boundaryconditions import NoSlip, UBB, FixedDensity, NeumannByCopy
-from lbmpy.boundaries.boundary_handling import BoundaryHandling
+from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
diff --git a/boundaries/boundary_handling.py b/boundaries/boundary_handling.py
deleted file mode 100644
index 083f863b842573e0018c83809c0900ae5d9aa8e0..0000000000000000000000000000000000000000
--- a/boundaries/boundary_handling.py
+++ /dev/null
@@ -1,337 +0,0 @@
-import numpy as np
-
-from lbmpy.stencils import inverseDirection
-from pystencils.cache import memorycache
-from pystencils import Field
-from lbmpy.boundaries import NoSlip
-from lbmpy.boundaries.boundary_kernel import generateIndexBoundaryKernelGeneric
-from lbmpy.boundaries.createindexlist import createBoundaryIndexArray, numpyDataTypeForBoundaryObject
-
-
-class BoundaryHandling:
-
- def __init__(self, lbMethod, dataHandling, pdfFieldName, name="boundaryHandling", target='cpu', openMP=True):
- assert dataHandling.hasData(pdfFieldName)
-
- self._lbMethod = lbMethod
- self._dataHandling = dataHandling
- self._pdfArrayName = pdfFieldName
- self._flagFieldName = name + "Flags"
- self._indexArrayName = name + "IndexArrays"
- self._target = target
- self._openMP = openMP
- self._boundaryObjectToBoundaryInfo = {}
- self._fluidFlag = 1 << 0
- self._nextFreeFlag = 1
-
- 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 shape(self):
- return self._dataHandling.shape
-
- @property
- def dim(self):
- return self._lbMethod.dim
-
- @property
- def lbMethod(self):
- return self._lbMethod
-
- @property
- def dataHandling(self):
- return self._dataHandling
-
- @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 forceOnBoundary(self, boundaryObj):
- if isinstance(boundaryObj, NoSlip):
- return self._forceOnNoSlip(boundaryObj)
- else:
- self.__call__()
- return self._forceOnBoundary(boundaryObj)
-
- 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._pdfArrayName] = b[self._pdfArrayName]
- 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 _forceOnNoSlip(self, boundaryObj):
- dh = self._dataHandling
- ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
- method = self._lbMethod
- stencil = np.array(method.stencil)
-
- result = np.zeros(self.dim)
-
- for b in dh.iterate(ghostLayers=ffGhostLayers):
- objToIndList = b[self._indexArrayName].boundaryObjectToIndexList
- pdfArray = b[self._pdfArrayName]
- if boundaryObj in objToIndList:
- indArr = objToIndList[boundaryObj]
- indices = [indArr[name] for name in ('x', 'y', 'z')[:method.dim]]
- indices.append(indArr['dir'])
- values = 2 * pdfArray[tuple(indices)]
- forces = stencil[indArr['dir']] * values[:, np.newaxis]
- result += forces.sum(axis=0)
- return dh.reduceFloatSequence(list(result), 'sum')
-
- def _forceOnBoundary(self, boundaryObj):
- dh = self._dataHandling
- ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
- method = self._lbMethod
- stencil = np.array(method.stencil)
- invDirection = np.array([method.stencil.index(inverseDirection(d))
- for d in method.stencil])
-
- result = np.zeros(self.dim)
-
- for b in dh.iterate(ghostLayers=ffGhostLayers):
- objToIndList = b[self._indexArrayName].boundaryObjectToIndexList
- pdfArray = b[self._pdfArrayName]
- if boundaryObj in objToIndList:
- indArr = objToIndList[boundaryObj]
- indices = [indArr[name] for name in ('x', 'y', 'z')[:method.dim]]
- offsets = stencil[indArr['dir']]
- invDir = invDirection[indArr['dir']]
- fluidValues = pdfArray[tuple(indices) + (indArr['dir'],)]
- boundaryValues = pdfArray[tuple(ind + offsets[:, i] for i, ind in enumerate(indices)) + (invDir,)]
- values = fluidValues + boundaryValues
- forces = stencil[indArr['dir']] * values[:, np.newaxis]
- result += forces.sum(axis=0)
-
- return dh.reduceFloatSequence(list(result), 'sum')
-
- def _getFlagForBoundary(self, boundaryObject):
- if boundaryObject not in self._boundaryObjectToBoundaryInfo:
- symbolicIndexField = Field.createGeneric('indexField', spatialDimensions=1,
- dtype=numpyDataTypeForBoundaryObject(boundaryObject, self.dim))
-
- spPdfField = self._dataHandling.fields[self._pdfArrayName]
- ast = generateIndexBoundaryKernelGeneric(spPdfField, symbolicIndexField, self._lbMethod,
- boundaryObject, target=self._target, openMP=self._openMP)
-
- boundaryInfo = self.BoundaryInfo(boundaryObject, flag=1 << self._nextFreeFlag, kernel=ast.compile())
-
- self._nextFreeFlag += 1
- self._boundaryObjectToBoundaryInfo[boundaryObject] = boundaryInfo
- return self._boundaryObjectToBoundaryInfo[boundaryObject].flag
-
- 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._pdfArrayName]
- indexArrayBD = b[self._indexArrayName]
- indexArrayBD.clear()
- for bInfo in self._boundaryObjectToBoundaryInfo.values():
- idxArr = createBoundaryIndexArray(flagArr, self.lbMethod.stencil, bInfo.flag, self._fluidFlag,
- bInfo.boundaryObject, dh.ghostLayersOfField(self._flagFieldName))
- if idxArr.size == 0:
- continue
-
- # TODO test boundary data setter (especially offset)
- boundaryDataSetter = BoundaryDataSetter(idxArr, b.offset, self.lbMethod.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]
\ No newline at end of file
diff --git a/boundaries/boundary_kernel.py b/boundaries/boundary_kernel.py
deleted file mode 100644
index d9e0f89c623856ef6dc03b73a39c93e2f470949b..0000000000000000000000000000000000000000
--- a/boundaries/boundary_kernel.py
+++ /dev/null
@@ -1,76 +0,0 @@
-import sympy as sp
-import numpy as np
-from pystencils.backends.cbackend import CustomCppCode
-from pystencils import TypedSymbol, Field
-from pystencils.cpu import addOpenMP
-from pystencils.data_types import castFunc, createType
-
-INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
-WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
-
-
-# -------------------------------------- Helper Functions --------------------------------------------------------------
-
-
-def offsetSymbols(dim):
- return [TypedSymbol("c_%d" % (d,), createType(np.int64)) for d in range(dim)]
-
-
-def offsetFromDir(dirIdx, dim):
- return tuple([sp.IndexedBase(symbol, shape=(1,))[dirIdx] for symbol in offsetSymbols(dim)])
-
-
-def invDir(dirIdx):
- return sp.IndexedBase(INV_DIR_SYMBOL, shape=(1,))[dirIdx]
-
-
-def weightOfDirection(dirIdx):
- return sp.IndexedBase(WEIGHTS_SYMBOL, shape=(1,))[dirIdx]
-
-
-# ------------------------------------- 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 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" % (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', openMP=True):
- indexField = Field.createFromNumpyArray("indexField", indexArr)
- return generateIndexBoundaryKernelGeneric(pdfField, indexField, lbMethod, boundaryFunctor, target, openMP)
-
-
-def generateIndexBoundaryKernelGeneric(pdfField, indexField, lbMethod, boundaryFunctor, target='cpu', openMP=True):
- elements = [LbmMethodInfo(lbMethod)]
- indexArrDtype = indexField.dtype.numpyDtype
- dirSymbol = TypedSymbol("dir", indexArrDtype.fields['dir'][0])
- boundaryEqList = [sp.Eq(dirSymbol, indexField[0]('dir'))]
- boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
- indexField=indexField)
- elements += boundaryEqList
-
- if target == 'cpu':
- from pystencils.cpu import createIndexedKernel
- ast = createIndexedKernel(elements, [indexField])
- addOpenMP(ast, numThreads=openMP)
- return ast
- elif target == 'gpu':
- from pystencils.gpucuda import createdIndexedCUDAKernel
- return createdIndexedCUDAKernel(elements, [indexField])
\ No newline at end of file
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index ab8caf25270836d276b871a462d5556b202a0c4f..3b84aeede7dcea506818f9f13e3c2ec5c2adaf27 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -1,10 +1,9 @@
import sympy as sp
-import numpy as np
from lbmpy.simplificationfactory import createSimplificationStrategy
from pystencils.astnodes import SympyAssignment
from pystencils.sympyextensions import getSymmetricPart
from pystencils import Field
-from lbmpy.boundaries.boundary_kernel import offsetFromDir, weightOfDirection, invDir
+from lbmpy.boundaries.boundaryhandling import BoundaryOffsetInfo, LbmWeightInfo
from pystencils.data_types import createType
@@ -63,8 +62,8 @@ class NoSlip(Boundary):
"""No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
- neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
- inverseDir = invDir(directionSymbol)
+ neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, lbMethod.dim)
+ inverseDir = BoundaryOffsetInfo.invDir(directionSymbol)
return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(directionSymbol))]
def __hash__(self):
@@ -118,8 +117,8 @@ class UBB(Boundary):
assert self.dim == lbMethod.dim, "Dimension of UBB (%d) does not match dimension of method (%d)" \
% (self.dim, lbMethod.dim)
- neighbor = offsetFromDir(direction, lbMethod.dim)
- inverseDir = invDir(direction)
+ neighbor = BoundaryOffsetInfo.offsetFromDir(direction, lbMethod.dim)
+ inverseDir = BoundaryOffsetInfo.invDir(direction)
velocity = tuple(v_i.getShifted(*neighbor) if isinstance(v_i, Field.Access) and not velFromIdxField else v_i
for v_i in vel)
@@ -130,6 +129,7 @@ class UBB(Boundary):
velocity = [eq.rhs for eq in shiftedVelEqs.extract(cqc.firstOrderMomentSymbols).mainEquations]
c_s_sq = sp.Rational(1, 3)
+ weightOfDirection = LbmWeightInfo.weightOfDirection
velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
# Better alternative: in conserved value computation
@@ -163,8 +163,8 @@ class FixedDensity(Boundary):
newMainEquations = [sp.Eq(e.lhs, getSymmetricPart(e.rhs, dofs)) for e in eqColl.mainEquations]
return eqColl.copy(newMainEquations)
- neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
- inverseDir = invDir(directionSymbol)
+ neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, lbMethod.dim)
+ inverseDir = BoundaryOffsetInfo.invDir(directionSymbol)
cqc = lbMethod.conservedQuantityComputation
velocity = cqc.definedSymbols()['velocity']
@@ -195,7 +195,7 @@ class FixedDensity(Boundary):
class NeumannByCopy(Boundary):
def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
- neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
+ neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, lbMethod.dim)
return [sp.Eq(pdfField[neighbor](directionSymbol), pdfField(directionSymbol))]
def __hash__(self):
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ab372a541d866b763316268ea92b48043d2b26a
--- /dev/null
+++ b/boundaries/boundaryhandling.py
@@ -0,0 +1,103 @@
+import numpy as np
+import sympy as sp
+from lbmpy.stencils import inverseDirection
+from pystencils import TypedSymbol, createIndexedKernel
+from pystencils.backends.cbackend import CustomCppCode
+from pystencils.boundaries import BoundaryHandling
+from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
+
+
+class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
+
+ def __init__(self, lbMethod, dataHandling, pdfFieldName, name="boundaryHandling", target='cpu', openMP=True):
+ self.lbMethod = lbMethod
+ super(LatticeBoltzmannBoundaryHandling, self).__init__(dataHandling, pdfFieldName, lbMethod.stencil,
+ name, target, openMP)
+
+ def forceOnBoundary(self, boundaryObj):
+ from lbmpy.boundaries import NoSlip
+ if isinstance(boundaryObj, NoSlip):
+ return self._forceOnNoSlip(boundaryObj)
+ else:
+ self.__call__()
+ return self._forceOnBoundary(boundaryObj)
+
+ # ------------------------------ Implementation Details ------------------------------------------------------------
+
+ def _forceOnNoSlip(self, boundaryObj):
+ dh = self._dataHandling
+ ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
+ method = self.lbMethod
+ stencil = np.array(method.stencil)
+
+ result = np.zeros(self.dim)
+
+ for b in dh.iterate(ghostLayers=ffGhostLayers):
+ objToIndList = b[self._indexArrayName].boundaryObjectToIndexList
+ pdfArray = b[self._fieldName]
+ if boundaryObj in objToIndList:
+ indArr = objToIndList[boundaryObj]
+ indices = [indArr[name] for name in ('x', 'y', 'z')[:method.dim]]
+ indices.append(indArr['dir'])
+ values = 2 * pdfArray[tuple(indices)]
+ forces = stencil[indArr['dir']] * values[:, np.newaxis]
+ result += forces.sum(axis=0)
+ return dh.reduceFloatSequence(list(result), 'sum')
+
+ def _forceOnBoundary(self, boundaryObj):
+ dh = self._dataHandling
+ ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
+ method = self.lbMethod
+ stencil = np.array(method.stencil)
+ invDirection = np.array([method.stencil.index(inverseDirection(d))
+ for d in method.stencil])
+
+ result = np.zeros(self.dim)
+
+ for b in dh.iterate(ghostLayers=ffGhostLayers):
+ objToIndList = b[self._indexArrayName].boundaryObjectToIndexList
+ pdfArray = b[self._fieldName]
+ if boundaryObj in objToIndList:
+ indArr = objToIndList[boundaryObj]
+ indices = [indArr[name] for name in ('x', 'y', 'z')[:method.dim]]
+ offsets = stencil[indArr['dir']]
+ invDir = invDirection[indArr['dir']]
+ fluidValues = pdfArray[tuple(indices) + (indArr['dir'],)]
+ boundaryValues = pdfArray[tuple(ind + offsets[:, i] for i, ind in enumerate(indices)) + (invDir,)]
+ values = fluidValues + boundaryValues
+ forces = stencil[indArr['dir']] * values[:, np.newaxis]
+ result += forces.sum(axis=0)
+
+ return dh.reduceFloatSequence(list(result), 'sum')
+
+ def _createBoundaryKernel(self, symbolicField, symbolicIndexField, boundaryObject):
+ return createLatticeBoltzmannBoundaryKernel(symbolicField, symbolicIndexField, self.lbMethod,
+ boundaryObject, target=self._target, openMP=self._openMP)
+
+
+class LbmWeightInfo(CustomCppCode):
+
+ # --------------------------- Functions to be used by boundaries --------------------------
+
+ @staticmethod
+ def weightOfDirection(dirIdx):
+ return sp.IndexedBase(LbmWeightInfo.WEIGHTS_SYMBOL, shape=(1,))[dirIdx]
+
+ # ---------------------------------- Internal ---------------------------------------------
+
+ WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
+
+ def __init__(self, lbMethod):
+ weights = [str(w.evalf()) for w in lbMethod.weights]
+ code = "const double %s [] = { %s };\n" % (LbmWeightInfo.WEIGHTS_SYMBOL.name, ",".join(weights))
+ super(LbmWeightInfo, self).__init__(code, symbolsRead=set(),
+ symbolsDefined=set([LbmWeightInfo.WEIGHTS_SYMBOL]))
+
+
+def createLatticeBoltzmannBoundaryKernel(pdfField, indexField, lbMethod, boundaryFunctor, target='cpu', openMP=True):
+ elements = [BoundaryOffsetInfo(lbMethod.stencil), LbmWeightInfo(lbMethod)]
+ indexArrDtype = indexField.dtype.numpyDtype
+ dirSymbol = TypedSymbol("dir", indexArrDtype.fields['dir'][0])
+ elements += [sp.Eq(dirSymbol, indexField[0]('dir'))]
+ elements += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod, indexField=indexField)
+ return createIndexedKernel(elements, [indexField], target=target, cpuOpenMP=openMP)
diff --git a/boundaries/createindexlist.py b/boundaries/createindexlist.py
deleted file mode 100644
index e3a971e261aeea4e106bc00bf59172556cc09ed4..0000000000000000000000000000000000000000
--- a/boundaries/createindexlist.py
+++ /dev/null
@@ -1,79 +0,0 @@
-import numpy as np
-import itertools
-import warnings
-
-try:
- import pyximport;
-
- pyximport.install()
- from lbmpy.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/lbstep.py b/lbstep.py
index 1cf28aae7d1d054247821002b2edd39cbb350d80..21da0682985ada3e394788ff08dc58285fab1c00 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -1,5 +1,5 @@
import numpy as np
-from lbmpy.boundaries.boundary_handling import BoundaryHandling
+from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
from lbmpy.creationfunctions import switchToSymbolicRelaxationRatesForEntropicMethods, createLatticeBoltzmannFunction, \
updateWithDefaultParameters
from lbmpy.simplificationfactory import createSimplificationStrategy
@@ -101,9 +101,9 @@ class LatticeBoltzmannStep:
# -- Boundary Handling & Synchronization ---
self._sync = dataHandling.synchronizationFunction([self._pdfArrName], methodParameters['stencil'], target)
- self._boundaryHandling = BoundaryHandling(self.method, self._dataHandling, self._pdfArrName,
- name=name + "_boundaryHandling",
- target=target, openMP=optimizationParams['openMP'])
+ self._boundaryHandling = LatticeBoltzmannBoundaryHandling(self.method, self._dataHandling, self._pdfArrName,
+ name=name + "_boundaryHandling",
+ target=target, openMP=optimizationParams['openMP'])
# -- Macroscopic Value Kernels
self._getterKernel, self._setterKernel = self._compilerMacroscopicSetterAndGetter()
diff --git a/scenarios.py b/scenarios.py
index 0ead321806614a33ae337e7d0a15d5b3096f0a57..4fe2ad872a1bdd7905e0d9fb96e0857f161d88b7 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -58,6 +58,7 @@ def createFullyPeriodicFlow(initialVelocity, periodicityInKernel=False, lbmKerne
step = LatticeBoltzmannStep(dataHandling=dataHandling, name="periodicScenario", lbmKernel=lbmKernel, **kwargs)
for b in step.dataHandling.iterate(ghostLayers=False):
np.copyto(b[step.velocityDataName], initialVelocity[b.globalSlice])
+ step.setPdfFieldsFromMacroscopicValues()
return step