From cd11756dc7f25e770fbb58b968ff75ba13594359 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 26 Jan 2018 12:16:49 +0100
Subject: [PATCH] New boundary handling
- plotting
- calculation of force on boundary
- vtk outpu
---
boundaries/__init__.py | 2 +-
.../boundary_handling.py | 204 ++++++++----
boundaries/boundaryconditions.py | 29 --
boundaries/boundaryhandling.py | 61 ----
boundaries/boundaryhandling_old.py | 280 ----------------
boundaries/forceevaluation.py | 52 ---
boundaries/handlinginterface.py | 309 ------------------
boundaries/periodicityhandling.py | 69 ----
geometry.py | 14 +-
lbstep.py | 45 +--
plot2d.py | 86 +++--
scenarios.py | 41 +--
12 files changed, 230 insertions(+), 962 deletions(-)
rename boundary_handling.py => boundaries/boundary_handling.py (68%)
delete mode 100644 boundaries/boundaryhandling.py
delete mode 100644 boundaries/boundaryhandling_old.py
delete mode 100644 boundaries/forceevaluation.py
delete mode 100644 boundaries/handlinginterface.py
delete mode 100644 boundaries/periodicityhandling.py
diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index 7588d247..6d1d37c0 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.boundaryhandling import BoundaryHandling
\ No newline at end of file
+from lbmpy.boundaries.boundary_handling import BoundaryHandling
diff --git a/boundary_handling.py b/boundaries/boundary_handling.py
similarity index 68%
rename from boundary_handling.py
rename to boundaries/boundary_handling.py
index 07e0c29a..79eb83e3 100644
--- a/boundary_handling.py
+++ b/boundaries/boundary_handling.py
@@ -1,52 +1,11 @@
import numpy as np
-from lbmpy.boundaries.boundary_kernel import generateIndexBoundaryKernelGeneric
+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
-from pystencils.cache import memorycache
-
-
-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 BoundaryHandling:
@@ -56,7 +15,7 @@ class BoundaryHandling:
self._lbMethod = lbMethod
self._dataHandling = dataHandling
- self._pdfFieldName = pdfFieldName
+ self._pdfArrayName = pdfFieldName
self._flagFieldName = name + "Flags"
self._indexArrayName = name + "IndexArrays"
self._target = target
@@ -73,13 +32,17 @@ class BoundaryHandling:
"If you want to add multiple handlings, choose a different name.")
gpu = self._target == 'gpu'
- dataHandling.addArray(self._flagFieldName, dtype=np.int32, cpu=True, gpu=gpu)
+ dataHandling.addArray(self._flagFieldName, dtype=np.uint32, cpu=True, gpu=gpu)
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
@@ -88,11 +51,24 @@ class BoundaryHandling:
def lbMethod(self):
return self._lbMethod
+ @property
+ def dataHandling(self):
+ return self._dataHandling
+
@property
def boundaryObjects(self):
return tuple(self._boundaryObjectToName.keys())
- def setBoundary(self, boundaryObject, sliceObj=None, maskCallback=None, ghostLayers=True):
+ @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 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
@@ -111,7 +87,7 @@ class BoundaryHandling:
"""
flag = self._getFlagForBoundary(boundaryObject)
- for b in self._dataHandling.iterate(sliceObj=sliceObj, ghostLayers=ghostLayers):
+ for b in self._dataHandling.iterate(sliceObj, ghostLayers=ghostLayers, innerGhostLayers=innerGhostLayers):
flagArr = b[self._flagFieldName]
if maskCallback is not None:
mask = maskCallback(*b.midpointArrays)
@@ -121,10 +97,12 @@ class BoundaryHandling:
self._dirty = True
- def getBoundaryNameToFlagDict(self):
- result = {bObj.name: bInfo.flag for bObj, bInfo in self._boundaryObjectToBoundaryInfo.values()}
- result['fluid'] = self._fluidFlag
- return result
+ 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:
@@ -147,17 +125,88 @@ class BoundaryHandling:
for b in self._dataHandling.iterate(gpu=self._target == 'gpu'):
for bObj, idxArr in b[self._indexArrayName].boundaryObjectToIndexList.items():
- kwargs[self._pdfFieldName] = b[self._pdfFieldName]
+ 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._pdfFieldName]
+ spPdfField = self._dataHandling.fields[self._pdfArrayName]
ast = generateIndexBoundaryKernelGeneric(spPdfField, symbolicIndexField, self._lbMethod,
boundaryObject, target=self._target, openMP=self._openMP)
@@ -172,7 +221,7 @@ class BoundaryHandling:
ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
for b in dh.iterate(ghostLayers=ffGhostLayers):
flagArr = b[self._flagFieldName]
- pdfArr = b[self._pdfFieldName]
+ pdfArr = b[self._pdfArrayName]
indexArrayBD = b[self._indexArrayName]
indexArrayBD.clear()
for bInfo in self._boundaryObjectToBoundaryInfo.values():
@@ -225,3 +274,46 @@ class BoundaryHandling:
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/boundaryconditions.py b/boundaries/boundaryconditions.py
index 9557acc5..ab8caf25 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -77,35 +77,6 @@ class NoSlip(Boundary):
return self.name == other.name
-#class NoSlipFullWay(Boundary):
-# """Full-way bounce back"""
-#
-# @property
-# def additionalData(self):
-# return [('lastValue', createType("double"))]
-#
-# @property
-# def additionalDataInitKernelEquations(self):
-# # TODO this function is not longer availabel
-# # formulate as additionalDataInitCallback
-# def kernelEqGetter(pdfField, directionSymbol, indexField, **kwargs):
-# return [sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
-# return kernelEqGetter
-#
-# def __call__(self, pdfField, directionSymbol, lbMethod, indexField, **kwargs):
-# neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
-# inverseDir = invDir(directionSymbol)
-# return [sp.Eq(pdfField[neighbor](inverseDir), indexField('lastValue')),
-# sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
-#
-# def __hash__(self):
-# # All boundaries of these class behave equal -> should also be equal
-# return hash("NoSlipFullWay")
-#
-# def __eq__(self, other):
-# return type(other) == NoSlipFullWay
-
-
class UBB(Boundary):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
deleted file mode 100644
index 5897cd17..00000000
--- a/boundaries/boundaryhandling.py
+++ /dev/null
@@ -1,61 +0,0 @@
-import numpy as np
-from lbmpy.boundaries.handlinginterface import GenericBoundaryHandling, FlagFieldInterface
-from lbmpy.boundaries.periodicityhandling import PeriodicityHandling
-
-
-class BoundaryHandling(PeriodicityHandling, GenericBoundaryHandling): # important: first periodicity, then boundary
-
- def __init__(self, pdfField, domainShape, lbMethod, ghostLayers=1, target='cpu', openMP=True, flagDtype=np.uint32):
- shapeWithGl = [a + 2 * ghostLayers for a in domainShape]
- self.domainShape = domainShape
- self.domainShapeWithGhostLayers = shapeWithGl
- flagInterface = NumpyFlagFieldInterface(shapeWithGl, flagDtype)
-
- GenericBoundaryHandling.__init__(self, flagInterface, pdfField, lbMethod, None, ghostLayers, target, openMP)
- PeriodicityHandling.__init__(self, list(domainShape) + [len(lbMethod.stencil)], target=target)
-
- def __call__(self, *args, **kwargs):
- for cls in BoundaryHandling.__bases__:
- cls.__call__(self, *args, **kwargs)
-
- def prepare(self):
- for cls in BoundaryHandling.__bases__:
- cls.prepare(self)
-
-# ----------------------------------------------------------------------------------------------------------------------
-
-
-class NumpyFlagFieldInterface(FlagFieldInterface):
-
- def __init__(self, shape, dtype=np.uint32):
- self._array = np.ones(shape, dtype=dtype)
- self._boundaryObjectToName = {}
- self._boundaryObjectToFlag = {'fluid': dtype(2**0)}
- self._nextFreeExponent = 1
-
- @property
- def array(self):
- return self._array
-
- def getFlag(self, boundaryObject):
- if boundaryObject not in self._boundaryObjectToFlag:
- self._boundaryObjectToFlag[boundaryObject] = 2 ** self._nextFreeExponent
- self._nextFreeExponent += 1
- name = self._makeBoundaryName(boundaryObject, self._boundaryObjectToName.values())
- self._boundaryObjectToName[boundaryObject] = name
-
- return self._boundaryObjectToFlag[boundaryObject]
-
- def getName(self, boundaryObject):
- return self._boundaryObjectToName[boundaryObject]
-
- @property
- def boundaryObjects(self):
- return self._boundaryObjectToName.keys()
-
- def clear(self):
- self._array.fill(0)
- self._boundaryObjectToName = {}
- self._boundaryObjectToFlag = {'fluid': np.dtype(self._array.dtype)(2**0)}
- self._nextFreeExponent = 1
-
diff --git a/boundaries/boundaryhandling_old.py b/boundaries/boundaryhandling_old.py
deleted file mode 100644
index 3be5b708..00000000
--- a/boundaries/boundaryhandling_old.py
+++ /dev/null
@@ -1,280 +0,0 @@
-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/forceevaluation.py b/boundaries/forceevaluation.py
deleted file mode 100644
index 3e741d73..00000000
--- a/boundaries/forceevaluation.py
+++ /dev/null
@@ -1,52 +0,0 @@
-import numpy as np
-from lbmpy.stencils import inverseDirection
-
-
-def calculateForceOnNoSlipBoundary(boundaryObject, boundaryHandling, pdfArray):
- indArr = boundaryHandling.getBoundaryIndexArray(boundaryObject)
- method = boundaryHandling.lbMethod
-
- stencil = np.array(method.stencil)
-
- if method.dim == 2:
- x, y = indArr['x'], indArr['y']
- values = 2 * pdfArray[x, y, indArr['dir']]
- else:
- assert method.dim == 3
- x, y, z = indArr['x'], indArr['y'], indArr['z']
- values = 2 * pdfArray[x, y, z, indArr['dir']]
-
- forces = stencil[indArr['dir']] * values[:, np.newaxis]
- return forces.sum(axis=0)
-
-
-def calculateForceOnBoundary(boundaryObject, boundaryHandling, pdfArray):
- indArr = boundaryHandling.getBoundaryIndexArray(boundaryObject)
- method = boundaryHandling.lbMethod
-
- stencil = np.array(method.stencil)
- invDirection = np.array([method.stencil.index(inverseDirection(d))
- for d in method.stencil])
-
- if method.dim == 2:
- x, y = indArr['x'], indArr['y']
- offsets = stencil[indArr['dir']]
-
- fluidValues = pdfArray[x, y, indArr['dir']]
- boundaryValues = pdfArray[x + offsets[:, 0],
- y + offsets[:, 1],
- invDirection[indArr['dir']]]
- else:
- assert method.dim == 3
- x, y, z = indArr['x'], indArr['y'], indArr['z']
- offsets = stencil[indArr['dir']]
-
- fluidValues = pdfArray[x, y, z, indArr['dir']]
- boundaryValues = pdfArray[x + offsets[:, 0],
- y + offsets[:, 1],
- z + offsets[:, 2],
- invDirection[indArr['dir']]]
-
- values = fluidValues + boundaryValues
- forces = stencil[indArr['dir']] * values[:, np.newaxis]
- return forces.sum(axis=0)
diff --git a/boundaries/handlinginterface.py b/boundaries/handlinginterface.py
deleted file mode 100644
index 4653feda..00000000
--- a/boundaries/handlinginterface.py
+++ /dev/null
@@ -1,309 +0,0 @@
-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, boundaryIndexArrayCoordinateNames,\
- numpyDataTypeForBoundaryObject
-from pystencils.cache import memorycache
-
-
-class FlagFieldInterface(object):
-
- def getFlag(self, boundaryObject):
- raise NotImplementedError()
-
- def getName(self, boundaryObject):
- raise NotImplementedError()
-
- @property
- def array(self):
- raise NotImplementedError()
-
- @property
- def boundaryObjects(self):
- raise NotImplementedError()
-
- def clear(self):
- raise NotImplementedError()
-
- @staticmethod
- def _makeBoundaryName(boundaryObject, existingNames):
- baseName = boundaryObject.name
- name = baseName
- counter = 1
- while name in existingNames:
- name = "%s_%d" % (baseName, counter)
- counter += 1
-
- return name
-
-
-class BoundaryDataSetter:
-
- def __init__(self, indexArray, offset, stencil, ghostLayers):
- self.indexArray = indexArray
- self.offset = offset
- self.stencil = np.array(stencil)
- 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 GenericBoundaryHandling(object):
-
- class BoundaryInfo(object):
- def __init__(self, kernel, ast, indexArray, idxArrayForExecution, boundaryDataSetter):
- self.kernel = kernel
- self.ast = ast
- self.indexArray = indexArray
- self.idxArrayForExecution = idxArrayForExecution # is different for GPU kernels
- self.boundaryDataSetter = boundaryDataSetter
-
- def __init__(self, flagFieldInterface, pdfField, lbMethod, offset=None, ghostLayers=1, target='cpu', openMP=True):
- """
- :param flagFieldInterface: implementation of FlagFieldInterface
- :param pdfField: numpy array
- :param lbMethod:
- :param offset: offset that is added to all coordinates when calling callback functions to set up geometry
- or boundary data. This is used for waLBerla simulations to pass the block offset
- :param target: 'cpu' or 'gpu'
- :param openMP:
- """
- self._flagFieldInterface = flagFieldInterface
- self._pdfField = pdfField
- self._lbMethod = lbMethod
- self._target = target
- self._openMP = openMP
- self.ghostLayers = ghostLayers
- self._dirty = False
- self.offset = offset if offset else (0,) * lbMethod.dim
- self._boundaryInfos = {} # mapping of boundary object to boundary info
- self._symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
- self.dim = lbMethod.dim
-
- if target not in ('cpu', 'gpu'):
- raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
-
- @property
- def boundaryObjects(self):
- return self._flagFieldInterface.boundaryObjects
-
- @property
- def lbMethod(self):
- return self._lbMethod
-
- def getBoundaryNameToFlagDict(self):
- result = {self._flagFieldInterface.getName(o): self._flagFieldInterface.getFlag(o) for o in self._boundaryInfos}
- result['fluid'] = self._flagFieldInterface.getFlag('fluid')
- return result
-
- def hasBoundary(self, boundaryObject):
- """Returns boolean indicating if this boundary exists in that handling"""
- return boundaryObject in self._boundaryInfos
-
- def __call__(self, **kwargs):
- """
- Runs boundary handling
- :param kwargs: keyword arguments passed to boundary kernel
- :return:
- """
- if self._dirty:
- self.prepare()
- for boundary in self._boundaryInfos.values():
- if boundary.kernel:
- boundary.kernel(indexField=boundary.idxArrayForExecution, **kwargs)
-
- def getBoundaryIndexArray(self, boundaryObject):
- return self._boundaryInfos[boundaryObject].indexArray
-
- def clear(self):
- """Removes all boundaries and fills the domain with fluid"""
- self._flagFieldInterface.clear()
- self._dirty = False
- self._boundaryInfos = {}
-
- def triggerReinitializationOfBoundaryData(self, **kwargs):
- if self._dirty:
- self.prepare()
- return
- else:
- for boundaryObject, boundaryInfo in self._boundaryInfos.items():
- self.__boundaryDataInitialization(boundaryInfo, boundaryObject, **kwargs)
-
- def prepare(self):
- """Compiles boundary kernels according to flag field. 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."""
- if not self._dirty:
- return
-
- dim = self._lbMethod.dim
- fluidFlag = self._flagFieldInterface.getFlag('fluid')
-
- for boundaryObject in self._flagFieldInterface.boundaryObjects:
- boundaryFlag = self._flagFieldInterface.getFlag(boundaryObject)
- idxArray = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
- boundaryFlag, fluidFlag, self.ghostLayers)
- if len(idxArray) == 0:
- continue
-
- 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
-
- ast = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod, boundaryObject,
- target=self._target)
-
- if self._target == 'cpu':
- from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
- addOpenMP(ast, numThreads=self._openMP)
- idxArrayForExecution = idxArray
- kernel = makePythonCpuFunction(ast)
- elif self._target == 'gpu':
- from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
- import pycuda.gpuarray as gpuarray
- idxArrayForExecution = gpuarray.to_gpu(idxArray)
- kernel = makePythonGpuFunction(ast)
- else:
- assert False
-
- boundaryDataSetter = BoundaryDataSetter(idxArray, self.offset, self._lbMethod.stencil, self.ghostLayers)
- boundaryInfo = GenericBoundaryHandling.BoundaryInfo(kernel, ast, idxArray,
- idxArrayForExecution, boundaryDataSetter)
- self._boundaryInfos[boundaryObject] = boundaryInfo
-
- if boundaryObject.additionalData:
- self.__boundaryDataInitialization(boundaryInfo, boundaryObject)
-
- self._dirty = False
-
- def getMask(self, boundaryObject):
- return np.bitwise_and(self._flagFieldInterface.array, self._flagFieldInterface.getFlag(boundaryObject))
-
- @property
- def flagField(self):
- return self._flagFieldInterface.array
-
- def reserveFlag(self, boundaryObject):
- self._flagFieldInterface.getFlag(boundaryObject)
-
- def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, includeGhostLayers=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 indexExpr: 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 includeGhostLayers: if this parameter is False, boundaries can not be set in the ghost
- layer, because index 0 is the first inner layer and -1 is interpreted in the Python
- way as maximum. If this parameter is True, the lower ghost layers have index 0, and
- placeholders ':' in index expressions extend into the ghost layers.
- """
- if indexExpr is None:
- indexExpr = [slice(None, None, None)] * len(self.flagField.shape)
- if not includeGhostLayers:
- domainSize = [i - 2 * self.ghostLayers for i in self._flagFieldInterface.array.shape]
- indexExpr = normalizeSlice(indexExpr, domainSize)
- indexExpr = shiftSlice(indexExpr, self.ghostLayers)
- else:
- indexExpr = normalizeSlice(indexExpr, self._flagFieldInterface.array.shape)
-
- mask = None
- if maskCallback is not None:
- gridParams = []
- for s, offset in zip(indexExpr, self.offset):
- if isinstance(s, slice):
- gridParams.append(np.arange(s.start, s.stop) + offset + 0.5 - self.ghostLayers)
- else:
- gridParams.append(s + offset + 0.5 - self.ghostLayers)
- indices = np.meshgrid(*gridParams, indexing='ij')
- mask = maskCallback(*indices)
- return self._setBoundaryWithMaskArray(boundaryObject, indexExpr, mask)
-
- def _setBoundaryWithMaskArray(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. ghost layers are expected to have coord=0
- :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)
-
- flag = self._flagFieldInterface.getFlag(boundaryObject)
- flagField = self._flagFieldInterface.array
- indexExpr = normalizeSlice(indexExpr, flagField.shape)
-
- if maskArr is None:
- flagField[indexExpr] = flag
- else:
- flagFieldView = flagField[indexExpr].squeeze()
- maskArr = maskArr.squeeze()
- flagFieldView[maskArr] = flag
- self._dirty = True
-
- def _invalidateCache(self):
- self._dirty = True
-
- def __boundaryDataInitialization(self, boundaryInfo, boundaryObject, **kwargs):
- if boundaryObject.additionalDataInitCallback:
- boundaryObject.additionalDataInitCallback(boundaryInfo.boundaryDataSetter, **kwargs)
-
- if boundaryObject.additionalDataInitKernelEquations:
- initKernelAst = generateIndexBoundaryKernel(self._symbolicPdfField, boundaryInfo.indexArray, self._lbMethod,
- boundaryObject, target='cpu',
- createInitializationKernel=True)
- from pystencils.cpu import makePythonFunction as makePythonCpuFunction
- initKernel = makePythonCpuFunction(initKernelAst, {'indexField': boundaryInfo.indexArray,
- 'pdfs': self._pdfField})
- initKernel()
-
- if self._target == 'gpu':
- import pycuda.gpuarray as gpuarray
- boundaryInfo.idxArrayForExecution = gpuarray.to_gpu(boundaryInfo.indexArray)
- else:
- boundaryInfo.idxArrayForExecution = boundaryInfo.indexArray
diff --git a/boundaries/periodicityhandling.py b/boundaries/periodicityhandling.py
deleted file mode 100644
index 428afcae..00000000
--- a/boundaries/periodicityhandling.py
+++ /dev/null
@@ -1,69 +0,0 @@
-from lbmpy.stencils import getStencil
-
-
-class PeriodicityHandling(object):
- def __init__(self, fieldShape, periodicity=(False, False, False), target='cpu'):
- self._spatialShape = fieldShape[:-1]
- self._indexShape = fieldShape[-1]
- self._periodicity = list(periodicity)
- self._periodicityDirty = True
- self._periodicityKernels = []
- self._target = target
-
- @property
- def periodicity(self):
- """List that indicates for x,y (z) coordinate if domain is periodic in that direction"""
- return self._periodicity
-
- 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._periodicityDirty = True
-
- def __call__(self, **kwargs):
- if self._periodicityDirty:
- self.prepare()
- for k in self._periodicityKernels:
- k(**kwargs)
-
- def prepare(self):
- if not self._periodicityDirty:
- return
-
- self._periodicityKernels = []
- dim = len(self._spatialShape)
- 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._spatialShape,
- indexDimensions=1,
- indexDimShape=self._indexShape,
- ghostLayers=1))
- else:
- assert False
-
- self._periodicityDirty = False
diff --git a/geometry.py b/geometry.py
index 404aed0f..7533bf98 100644
--- a/geometry.py
+++ b/geometry.py
@@ -33,7 +33,7 @@ def addParabolicVelocityInflow(boundaryHandling, u_max, indexExpr, velCoord=0, d
if i != velCoord:
boundaryData[name] = 0.0
if diameter is None:
- radius = int(round(min(sh for i, sh in enumerate(boundaryHandling.domainShape) if i != velCoord) / 2))
+ radius = int(round(min(sh for i, sh in enumerate(boundaryHandling.shape) if i != velCoord) / 2))
else:
radius = int(round(diameter / 2))
@@ -41,8 +41,8 @@ def addParabolicVelocityInflow(boundaryHandling, u_max, indexExpr, velCoord=0, d
normalCoord1 = (velCoord + 1) % 3
normalCoord2 = (velCoord + 2) % 3
y, z = boundaryData.linkPositions(normalCoord1), boundaryData.linkPositions(normalCoord2)
- centeredNormal1 = y - int(round(boundaryHandling.domainShape[normalCoord1] / 2))
- centeredNormal2 = z - int(round(boundaryHandling.domainShape[normalCoord2] / 2))
+ centeredNormal1 = y - int(round(boundaryHandling.shape[normalCoord1] / 2))
+ centeredNormal2 = z - int(round(boundaryHandling.shape[normalCoord2] / 2))
distToCenter = np.sqrt(centeredNormal1 ** 2 + centeredNormal2 ** 2)
elif dim == 2:
normalCoord = (velCoord + 1) % 2
@@ -68,7 +68,7 @@ def setupChannelWalls(boundaryHandling, diameterCallback, duct=False, wallBounda
raise ValueError("For duct flows, passing a diameter callback does not make sense.")
if not duct:
- diameter = min(boundaryHandling.domainShape[1:])
+ diameter = min(boundaryHandling.shape[1:])
addPipe(boundaryHandling, diameterCallback if diameterCallback else diameter, wallBoundary)
@@ -82,7 +82,7 @@ def addPipe(boundaryHandling, diameter, boundary=NoSlip()):
a array of same shape as the received xCoordArray, with the diameter for each x position
:param boundary: boundary object that is set at the wall, defaults to NoSlip (bounce back)
"""
- domainShape = boundaryHandling.domainShape
+ domainShape = boundaryHandling.shape
dim = len(domainShape)
assert dim in (2, 3)
@@ -123,11 +123,11 @@ def addBlackAndWhiteImage(boundaryHandling, imageFile, targetSlice=None, plane=(
except ImportError:
raise ImportError("scipy image read could not be imported! Install 'scipy' and 'pillow'")
- domainSize = boundaryHandling.domainShape
+ domainSize = boundaryHandling.shape
if targetSlice is None:
targetSlice = [slice(None, None, None)] * len(domainSize)
- dim = len(boundaryHandling.domainShape)
+ dim = boundaryHandling.dim
imageSlice = normalizeSlice(targetSlice, domainSize)
targetSize = [imageSlice[i].stop - imageSlice[i].start for i in plane]
diff --git a/lbstep.py b/lbstep.py
index b3a9a55d..2008b434 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -1,5 +1,5 @@
import numpy as np
-from lbmpy.boundary_handling import BoundaryHandling
+from lbmpy.boundaries.boundary_handling import BoundaryHandling
from lbmpy.creationfunctions import switchToSymbolicRelaxationRatesForEntropicMethods, createLatticeBoltzmannFunction, \
updateWithDefaultParameters
from lbmpy.simplificationfactory import createSimplificationStrategy
@@ -11,7 +11,7 @@ from pystencils import createKernel, makeSlice
class LatticeBoltzmannStep:
def __init__(self, domainSize=None, lbmKernel=None, periodicity=False,
- kernelParams={}, dataHandling=None, dataPrefix="", optimizationParams={}, **methodParameters):
+ kernelParams={}, dataHandling=None, name="lbm", optimizationParams={}, **methodParameters):
# --- Parameter normalization ---
if dataHandling is not None:
@@ -36,10 +36,10 @@ class LatticeBoltzmannStep:
Q = len(getStencil(methodParameters['stencil']))
self._dataHandling = dataHandling
- self._pdfArrName = dataPrefix + "pdfSrc"
- self._tmpArrName = dataPrefix + "pdfTmp"
- self.velocityDataName = dataPrefix + "velocity"
- self.densityDataName = dataPrefix + "density"
+ self._pdfArrName = name + "_pdfSrc"
+ self._tmpArrName = name + "_pdfTmp"
+ self.velocityDataName = name + "_velocity"
+ self.densityDataName = name + "_density"
self._gpu = target == 'gpu'
layout = optimizationParams['fieldLayout']
@@ -67,7 +67,7 @@ class LatticeBoltzmannStep:
else:
self._sync = dataHandling.synchronizationFunctionCPU([self._pdfArrName], methodParameters['stencil'])
self._boundaryHandling = BoundaryHandling(self._lbmKernel.method, self._dataHandling, self._pdfArrName,
- name=dataPrefix + "_boundaryHandling",
+ name=name + "_boundaryHandling",
target=target, openMP=optimizationParams['openMP'])
# -- Macroscopic Value Kernels
@@ -79,6 +79,9 @@ class LatticeBoltzmannStep:
b[self.densityDataName].fill(1.0)
b[self.velocityDataName].fill(0.0)
+ # -- VTK output
+ self.vtkWriter = self.dataHandling.vtkWriter(name, [self.velocityDataName, self.densityDataName])
+
@property
def boundaryHandling(self):
"""Boundary handling instance of the scenario. Use this to change the boundary setup"""
@@ -96,6 +99,10 @@ class LatticeBoltzmannStep:
def method(self):
return self._lbmKernel.method
+ @property
+ def pdfArrayName(self):
+ return self._pdfArrName
+
def _getSlice(self, dataName, sliceObj):
if sliceObj is None:
sliceObj = makeSlice[:, :] if self.dim == 2 else makeSlice[:, :, 0.5]
@@ -132,6 +139,9 @@ class LatticeBoltzmannStep:
self.timeStep()
self.postRun()
+ def writeVTK(self):
+ self.vtkWriter(self.timeStepsRun)
+
def _compilerMacroscopicSetterAndGetter(self):
lbMethod = self._lbmKernel.method
D = lbMethod.dim
@@ -153,24 +163,3 @@ class LatticeBoltzmannStep:
setterEqs = createSimplificationStrategy(lbMethod)(setterEqs)
setterKernel = createKernel(setterEqs, target='cpu').compile()
return getterKernel, setterKernel
-
-
-if __name__ == '__main__':
- from pycuda import autoinit
- from lbmpy.boundaries import NoSlip, UBB
- from pystencils import makeSlice
- step = LatticeBoltzmannStep((30, 30), relaxationRate=1.8, periodicity=True,
- optimizationParams={'target': 'cpu', 'openMP': False})
-
- wall = NoSlip()
- movingWall = UBB((0.001, 0))
-
- bh = step.boundaryHandling
- bh.setBoundary(wall, makeSlice[0, :])
- bh.setBoundary(wall, makeSlice[-1, :])
- bh.setBoundary(wall, makeSlice[:, 0])
- bh.setBoundary(movingWall, makeSlice[:, -1])
- bh.prepare()
-
- step.run(4)
- step.run(100)
diff --git a/plot2d.py b/plot2d.py
index 8b0a3634..a08bbaa3 100644
--- a/plot2d.py
+++ b/plot2d.py
@@ -1,56 +1,54 @@
from pystencils.plot2d import *
-from pystencils.slicing import normalizeSlice
-def plotBoundaryHandling(boundaryHandling, indexExpr=None, boundaryNameToColor=None, legend=True):
+def boundaryHandling(boundaryHandlingObj, indexExpr=None, boundaryNameToColor=None, showLegend=True):
"""
Shows boundary cells
- :param boundaryHandling: instance of :class:`lbmpy.boundaries.BoundaryHandling`
+ :param boundaryHandlingObj: instance of :class:`lbmpy.boundaries.BoundaryHandling`
:param indexExpr: for 3D boundary handling a slice expression has to be passed here to define the plane that
should be plotted
:param boundaryNameToColor: optional dictionary mapping boundary names to colors
- :param legend: if True legend for color->boundary name is added
+ :param showLegend: if True legend for color->boundary name is added
"""
import matplotlib.pyplot as plt
- boundaryHandling.prepare()
-
- if len(boundaryHandling.flagField.shape) != 2 and indexExpr is None:
- raise ValueError("To plot a 3D boundary handling a slice has to be passed")
-
- if boundaryNameToColor:
- fixedColors = boundaryNameToColor
- else:
- fixedColors = {
- 'fluid': '#56b4e9',
- 'NoSlip': '#999999',
- 'UBB': '#d55e00',
- 'FixedDensity': '#009e73',
- }
-
- boundaryNames = []
- flagValues = []
- for name, flagName in sorted(boundaryHandling.getBoundaryNameToFlagDict().items(), key=lambda l: l[1]):
- boundaryNames.append(name)
- flagValues.append(flagName)
- defaultCycle = matplotlib.rcParams['axes.prop_cycle']
- colorValues = [fixedColors[name] if name in fixedColors else cycle['color']
- for cycle, name in zip(defaultCycle, boundaryNames)]
-
- cmap = matplotlib.colors.ListedColormap(colorValues)
- bounds = np.array(flagValues, dtype=float) - 0.5
- bounds = list(bounds) + [bounds[-1] + 1]
- norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
-
- flagField = boundaryHandling.flagField
- if indexExpr:
- flagField = flagField[normalizeSlice(indexExpr, flagField.shape)]
- flagField = flagField.swapaxes(0, 1)
- plt.imshow(flagField, interpolation='none', origin='lower',
- cmap=cmap, norm=norm)
-
- patches = [matplotlib.patches.Patch(color=color, label=name) for color, name in zip(colorValues, boundaryNames)]
- plt.axis('equal')
- if legend:
- plt.legend(handles=patches, bbox_to_anchor=(1.02, 0.5), loc=2, borderaxespad=0.)
+ boundaryHandlingObj.prepare()
+
+ dh = boundaryHandlingObj.dataHandling
+ for flagArr in dh.gatherArray(boundaryHandlingObj.flagArrayName, indexExpr):
+ if len(flagArr.shape) != 2 and indexExpr is None:
+ raise ValueError("To plot a 3D boundary handling a slice has to be passed")
+
+ if boundaryNameToColor:
+ fixedColors = boundaryNameToColor
+ else:
+ fixedColors = {
+ 'fluid': '#56b4e9',
+ 'NoSlip': '#999999',
+ 'UBB': '#d55e00',
+ 'FixedDensity': '#009e73',
+ }
+
+ boundaryNames = []
+ flagValues = []
+ for name, flagName in sorted(boundaryHandlingObj.getBoundaryNameToFlagDict().items(), key=lambda l: l[1]):
+ boundaryNames.append(name)
+ flagValues.append(flagName)
+ defaultCycle = matplotlib.rcParams['axes.prop_cycle']
+ colorValues = [fixedColors[name] if name in fixedColors else cycle['color']
+ for cycle, name in zip(defaultCycle, boundaryNames)]
+
+ cmap = matplotlib.colors.ListedColormap(colorValues)
+ bounds = np.array(flagValues, dtype=float) - 0.5
+ bounds = list(bounds) + [bounds[-1] + 1]
+ norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
+
+ flagArr = flagArr.swapaxes(0, 1)
+ plt.imshow(flagArr, interpolation='none', origin='lower',
+ cmap=cmap, norm=norm)
+
+ patches = [matplotlib.patches.Patch(color=color, label=name) for color, name in zip(colorValues, boundaryNames)]
+ plt.axis('equal')
+ if showLegend:
+ plt.legend(handles=patches, bbox_to_anchor=(1.02, 0.5), loc=2, borderaxespad=0.)
diff --git a/scenarios.py b/scenarios.py
index a3f393ec..39b3b75b 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -55,7 +55,7 @@ def createFullyPeriodicFlow(initialVelocity, periodicityInKernel=False, lbmKerne
if dataHandling is None:
dataHandling = createDataHandling(parallel, domainSize, periodicity=not periodicityInKernel,
defaultGhostLayers=1)
- step = LatticeBoltzmannStep(dataHandling=dataHandling, lbmKernel=lbmKernel, **kwargs)
+ 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])
return step
@@ -76,7 +76,7 @@ def createLidDrivenCavity(domainSize=None, lidVelocity=0.005, lbmKernel=None, pa
assert domainSize is not None or dataHandling is not None
if dataHandling is None:
dataHandling = createDataHandling(parallel, domainSize, periodicity=False, defaultGhostLayers=1)
- step = LatticeBoltzmannStep(dataHandling=dataHandling, lbmKernel=lbmKernel, **kwargs)
+ step = LatticeBoltzmannStep(dataHandling=dataHandling, lbmKernel=lbmKernel, name="lidDrivenCavity" **kwargs)
myUbb = UBB(velocity=[lidVelocity, 0, 0][:step.method.dim])
step.boundaryHandling.setBoundary(myUbb, sliceFromDirection('N', step.dim))
@@ -86,7 +86,7 @@ def createLidDrivenCavity(domainSize=None, lidVelocity=0.005, lbmKernel=None, pa
return step
-def createChannel(domainSize, force=None, pressureDifference=None, u_max=None, diameterCallback=None,
+def createChannel(domainSize=None, force=None, pressureDifference=None, u_max=None, diameterCallback=None,
duct=False, wallBoundary=NoSlip(), parallel=False, dataHandling=None, **kwargs):
"""
Create a channel scenario (2D or 3D)
@@ -106,6 +106,8 @@ def createChannel(domainSize, force=None, pressureDifference=None, u_max=None, d
:param parallel: True for distributed memory parallelization with waLBerla
:param kwargs: all other keyword parameters are passed directly to scenario class.
"""
+ assert domainSize is not None or dataHandling is not None
+
dim = len(domainSize)
assert dim in (2, 3)
@@ -115,42 +117,29 @@ def createChannel(domainSize, force=None, pressureDifference=None, u_max=None, d
periodicity = (True, False, False) if force else (False, False, False)
if dataHandling is None:
dataHandling = createDataHandling(parallel, domainSize, periodicity=periodicity[:dim], defaultGhostLayers=1)
- step = LatticeBoltzmannStep(dataHandling=dataHandling, **kwargs)
- boundaryHandling = step.boundaryHandling
+
if force:
kwargs['force'] = tuple([force, 0, 0][:dim])
- boundaryHandling.setPeriodicity(True, False, False)
+ assert dataHandling.periodicity[0]
+ step = LatticeBoltzmannStep(dataHandling=dataHandling, name="forceDrivenChannel", **kwargs)
elif pressureDifference:
inflow = FixedDensity(1.0 + pressureDifference)
outflow = FixedDensity(1.0)
- boundaryHandling.setBoundary(inflow, sliceFromDirection('W', dim))
- boundaryHandling.setBoundary(outflow, sliceFromDirection('E', dim))
+ step = LatticeBoltzmannStep(dataHandling=dataHandling, name="pressureDrivenChannel", **kwargs)
+ step.boundaryHandling.setBoundary(inflow, sliceFromDirection('W', dim))
+ step.boundaryHandling.setBoundary(outflow, sliceFromDirection('E', dim))
elif u_max:
if duct:
raise NotImplementedError("Velocity inflow for duct flows not yet implemented")
-
+ step = LatticeBoltzmannStep(dataHandling=dataHandling, name="velocityDrivenChannel", **kwargs)
diameter = diameterCallback(np.array([0]), domainSize)[0] if diameterCallback else min(domainSize[1:])
- addParabolicVelocityInflow(boundaryHandling, u_max, sliceFromDirection('W', dim),
+ addParabolicVelocityInflow(step.boundaryHandling, u_max, sliceFromDirection('W', dim),
velCoord=0, diameter=diameter)
outflow = FixedDensity(1.0)
- boundaryHandling.setBoundary(outflow, sliceFromDirection('E', dim))
+ step.boundaryHandling.setBoundary(outflow, sliceFromDirection('E', dim))
else:
assert False
- setupChannelWalls(boundaryHandling, diameterCallback, duct, wallBoundary)
+ setupChannelWalls(step.boundaryHandling, diameterCallback, duct, wallBoundary)
return step
-
-if __name__ == '__main__':
- import pycuda.autoinit
- import waLBerla as wlb
- from pystencils.datahandling import ParallelDataHandling
- blocks = wlb.createUniformBlockGrid(blocks=(2, 2, 1), cellsPerBlock=(20, 20, 1), oneBlockPerProcess=False)
- dh = ParallelDataHandling(blocks, dim=2)
- ldc = createLidDrivenCavity(relaxationRate=1.5, dataHandling=dh,
- optimizationParams={'openMP': 2, 'target': 'gpu', 'fieldLayout': 'f',
- 'gpuIndexingParams': {'blockSize': (8, 8, 2)}})
- ldc.boundaryHandling.prepare()
-
- ldc.run(4)
- ldc.run(100)
--
GitLab