Skip to content
Snippets Groups Projects
Commit 32197bd0 authored by Martin Bauer's avatar Martin Bauer
Browse files

lbmpy: now can build periodicity into kernel

parent d42bd4f0
Branches
Tags
No related merge requests found
...@@ -72,7 +72,10 @@ Simplifications / Transformations: ...@@ -72,7 +72,10 @@ Simplifications / Transformations:
- ``doOverallCse=False``: run common subexpression elimination after all other simplifications have been executed - ``doOverallCse=False``: run common subexpression elimination after all other simplifications have been executed
- ``split=False``: split innermost loop, to handle only 2 directions per loop. This reduces the number of parallel - ``split=False``: split innermost loop, to handle only 2 directions per loop. This reduces the number of parallel
load/store streams and thus speeds up the kernel on most architectures load/store streams and thus speeds up the kernel on most architectures
- ``builtinPeriodicity=(False,False,False)``: instead of handling periodicity by copying ghost layers, the periodicity
is built into the kernel. This parameters specifies if the domain is periodic in (x,y,z) direction. Even if the
periodicity is built into the kernel, the fields have one ghost layer to be consistent with other functions.
Field size information: Field size information:
...@@ -183,6 +186,8 @@ def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True): ...@@ -183,6 +186,8 @@ def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
'gpuIndexing': 'block', 'gpuIndexing': 'block',
'gpuIndexingParams': {}, 'gpuIndexingParams': {},
'builtinPeriodicity': (False, False, False),
} }
if 'relaxationRate' in params: if 'relaxationRate' in params:
if 'relaxationRates' not in params: if 'relaxationRates' not in params:
...@@ -242,7 +247,8 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs): ...@@ -242,7 +247,8 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
else: else:
splitGroups = () splitGroups = ()
res = createKernel(updateRule.allEquations, splitGroups=splitGroups, res = createKernel(updateRule.allEquations, splitGroups=splitGroups,
typeForSymbol='double' if optParams['doublePrecision'] else 'float') typeForSymbol='double' if optParams['doublePrecision'] else 'float',
ghostLayers=1)
elif optParams['target'] == 'gpu': elif optParams['target'] == 'gpu':
from pystencils.gpucuda import createCUDAKernel from pystencils.gpucuda import createCUDAKernel
from pystencils.gpucuda.indexing import LineIndexing, BlockIndexing from pystencils.gpucuda.indexing import LineIndexing, BlockIndexing
...@@ -252,7 +258,7 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs): ...@@ -252,7 +258,7 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
indexingCreator = partial(indexingCreator, **optParams['gpuIndexingParams']) indexingCreator = partial(indexingCreator, **optParams['gpuIndexingParams'])
res = createCUDAKernel(updateRule.allEquations, res = createCUDAKernel(updateRule.allEquations,
typeForSymbol='double' if optParams['doublePrecision'] else 'float', typeForSymbol='double' if optParams['doublePrecision'] else 'float',
indexingCreator=indexingCreator) indexingCreator=indexingCreator, ghostLayers=1)
else: else:
return ValueError("'target' has to be either 'cpu' or 'gpu'") return ValueError("'target' has to be either 'cpu' or 'gpu'")
...@@ -289,16 +295,19 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa ...@@ -289,16 +295,19 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
if 'fieldSize' in optParams and optParams['fieldSize']: if 'fieldSize' in optParams and optParams['fieldSize']:
npField = createPdfArray(optParams['fieldSize'], len(stencil), layout=optParams['fieldLayout']) npField = createPdfArray(optParams['fieldSize'], len(stencil), layout=optParams['fieldLayout'])
updateRule = createStreamPullKernel(collisionRule, numpyField=npField) updateRule = createStreamPullKernel(collisionRule, numpyField=npField,
builtinPeriodicity=optParams['builtinPeriodicity'])
else: else:
if 'pdfArr' in optParams: if 'pdfArr' in optParams:
updateRule = createStreamPullKernel(collisionRule, numpyField=optParams['pdfArr']) updateRule = createStreamPullKernel(collisionRule, numpyField=optParams['pdfArr'],
builtinPeriodicity=optParams['builtinPeriodicity'])
else: else:
layoutName = optParams['fieldLayout'] layoutName = optParams['fieldLayout']
if layoutName == 'fzyx' or 'zyxf': if layoutName == 'fzyx' or 'zyxf':
dim = len(stencil[0]) dim = len(stencil[0])
layoutName = tuple(reversed(range(dim))) layoutName = tuple(reversed(range(dim)))
updateRule = createStreamPullKernel(collisionRule, genericLayout=layoutName) updateRule = createStreamPullKernel(collisionRule, genericLayout=layoutName,
builtinPeriodicity=optParams['builtinPeriodicity'])
return updateRule return updateRule
......
import sympy as sp
import abc
from lbmpy.stencils import inverseDirection from lbmpy.stencils import inverseDirection
from pystencils import Field from pystencils import Field
import abc
# ------------------------------------------------ Interface ----------------------------------------------------------- # ------------------------------------------------ Interface -----------------------------------------------------------
from pystencils.astnodes import LoopOverCoordinate
class PdfFieldAccessor(object): class PdfFieldAccessor(object):
...@@ -50,6 +52,67 @@ class StreamPullTwoFieldsAccessor(PdfFieldAccessor): ...@@ -50,6 +52,67 @@ class StreamPullTwoFieldsAccessor(PdfFieldAccessor):
return [field(i) for i in range(len(stencil))] return [field(i) for i in range(len(stencil))]
class Pseudo2DTwoFieldsAccessor(PdfFieldAccessor):
"""Useful if a 3D simulation of a domain with size (x,y,1) is done and the dimensions with size 1
is periodic. In this case no periodicity exchange has to be done"""
def __init__(self, collapsedDim):
self._collapsedDim = collapsedDim
def read(self, field, stencil):
result = []
for i, d in enumerate(stencil):
direction = list(d)
direction[self._collapsedDim] = 0
result.append(field[inverseDirection(tuple(direction))](i))
return result
@staticmethod
def write(field, stencil):
return [field(i) for i in range(len(stencil))]
class PeriodicTwoFieldsAccessor(PdfFieldAccessor):
"""Access scheme that builds periodicity into the kernel, by introducing a condition on every load,
such that at the borders the periodic value is loaded. The periodicity is specified as a tuple of booleans, one for
each direction. The second parameter `ghostLayers` specifies the number of assumed ghost layers of the field.
For the periodic kernel itself no ghost layers are required, however other kernels might need them.
"""
def __init__(self, periodicity, ghostLayers=0):
self._periodicity = periodicity
self._ghostLayers = ghostLayers
def read(self, field, stencil):
result = []
for i, d in enumerate(stencil):
pullDirection = inverseDirection(d)
periodicPullDirection = []
for coordId, dirElement in enumerate(pullDirection):
if not self._periodicity[coordId]:
periodicPullDirection.append(dirElement)
continue
lowerLimit = self._ghostLayers
upperLimit = field.spatialShape[coordId] - 1 - self._ghostLayers
limitDiff = upperLimit - lowerLimit
loopCounter = LoopOverCoordinate.getLoopCounterSymbol(coordId)
if dirElement == 0:
periodicPullDirection.append(0)
elif dirElement == 1:
newDirElement = sp.Piecewise((dirElement, loopCounter < upperLimit), (-limitDiff, True))
periodicPullDirection.append(newDirElement)
elif dirElement == -1:
newDirElement = sp.Piecewise((dirElement, loopCounter > lowerLimit), (limitDiff, True))
periodicPullDirection.append(newDirElement)
else:
raise NotImplementedError("This accessor supports only nearest neighbor stencils")
result.append(field[tuple(periodicPullDirection)](i))
return result
@staticmethod
def write(field, stencil):
return [field(i) for i in range(len(stencil))]
class AABBEvenTimeStepAccessor(PdfFieldAccessor): class AABBEvenTimeStepAccessor(PdfFieldAccessor):
@staticmethod @staticmethod
def read(field, stencil): def read(field, stencil):
......
...@@ -37,21 +37,28 @@ from lbmpy.updatekernels import createPdfArray ...@@ -37,21 +37,28 @@ from lbmpy.updatekernels import createPdfArray
# ---------------------------------------- Example Scenarios ----------------------------------------------------------- # ---------------------------------------- Example Scenarios -----------------------------------------------------------
def createFullyPeriodicFlow(initialVelocity, optimizationParams={}, lbmKernel=None, kernelParams={}, **kwargs): def createFullyPeriodicFlow(initialVelocity, periodicityInKernel=False,
optimizationParams={}, lbmKernel=None, kernelParams={}, **kwargs):
""" """
Creates a fully periodic setup with prescribed velocity field Creates a fully periodic setup with prescribed velocity field
:param initialVelocity: numpy array that defines an initial velocity for each cell. The shape of this :param initialVelocity: numpy array that defines an initial velocity for each cell. The shape of this
array determines the domain size. array determines the domain size.
:param periodicityInKernel: don't use boundary handling for periodicity, but directly generate the kernel periodic
:param optimizationParams: see :mod:`lbmpy.creationfunctions` :param optimizationParams: see :mod:`lbmpy.creationfunctions`
:param lbmKernel: a LBM function, which would otherwise automatically created :param lbmKernel: a LBM function, which would otherwise automatically created
:param kernelParams: additional parameters passed to the sweep :param kernelParams: additional parameters passed to the sweep
:param kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions` :param kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
:return: instance of :class:`Scenario` :return: instance of :class:`Scenario`
""" """
optimizationParams = optimizationParams.copy()
domainSize = initialVelocity.shape[:-1] domainSize = initialVelocity.shape[:-1]
if periodicityInKernel:
optimizationParams['builtinPeriodicity'] = (True, True, True)
scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel, initialVelocity, kernelParams=kernelParams) scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel, initialVelocity, kernelParams=kernelParams)
scenario.boundaryHandling.setPeriodicity(True, True, True)
if not periodicityInKernel:
scenario.boundaryHandling.setPeriodicity(True, True, True)
return scenario return scenario
......
...@@ -3,7 +3,7 @@ import sympy as sp ...@@ -3,7 +3,7 @@ import sympy as sp
from pystencils import Field from pystencils import Field
from pystencils.field import createNumpyArrayWithLayout from pystencils.field import createNumpyArrayWithLayout
from pystencils.sympyextensions import fastSubs from pystencils.sympyextensions import fastSubs
from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor, Pseudo2DTwoFieldsAccessor, PeriodicTwoFieldsAccessor
# -------------------------------------------- LBM Kernel Creation ----------------------------------------------------- # -------------------------------------------- LBM Kernel Creation -----------------------------------------------------
...@@ -44,7 +44,8 @@ def createLBMKernel(collisionRule, inputField, outputField, accessor): ...@@ -44,7 +44,8 @@ def createLBMKernel(collisionRule, inputField, outputField, accessor):
def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", dstFieldName="dst", def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", dstFieldName="dst",
genericLayout='numpy', genericFieldType=np.float64): genericLayout='numpy', genericFieldType=np.float64,
builtinPeriodicity=(False, False, False)):
""" """
Implements a stream-pull scheme, where values are read from source and written to destination field Implements a stream-pull scheme, where values are read from source and written to destination field
:param collisionRule: a collision rule created by lbm method :param collisionRule: a collision rule created by lbm method
...@@ -55,6 +56,7 @@ def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", d ...@@ -55,6 +56,7 @@ def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", d
:param genericLayout: if no numpyField is given to determine the layout, a variable sized field with the given :param genericLayout: if no numpyField is given to determine the layout, a variable sized field with the given
genericLayout is used genericLayout is used
:param genericFieldType: if no numpyField is given, this data type is used for the fields :param genericFieldType: if no numpyField is given, this data type is used for the fields
:param builtinPeriodicity: periodicity that should be built into the kernel
:return: lbm update rule, where generic pdf references are replaced by field accesses :return: lbm update rule, where generic pdf references are replaced by field accesses
""" """
dim = collisionRule.method.dim dim = collisionRule.method.dim
...@@ -69,7 +71,11 @@ def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", d ...@@ -69,7 +71,11 @@ def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", d
src = Field.createFromNumpyArray(srcFieldName, numpyField, indexDimensions=1) src = Field.createFromNumpyArray(srcFieldName, numpyField, indexDimensions=1)
dst = Field.createFromNumpyArray(dstFieldName, numpyField, indexDimensions=1) dst = Field.createFromNumpyArray(dstFieldName, numpyField, indexDimensions=1)
return createLBMKernel(collisionRule, src, dst, StreamPullTwoFieldsAccessor) accessor = StreamPullTwoFieldsAccessor
if any(builtinPeriodicity):
accessor = PeriodicTwoFieldsAccessor(builtinPeriodicity, ghostLayers=1)
return createLBMKernel(collisionRule, src, dst, accessor)
# ---------------------------------- Pdf array creation for various layouts -------------------------------------------- # ---------------------------------- Pdf array creation for various layouts --------------------------------------------
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment