diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index d43c00c944d84aaab02c1522011f600d7f69f12f..41e77228233b7ed056fbfff6d6b8b8a7fbc029ed 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -1,5 +1,7 @@
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
@@ -37,10 +39,19 @@ class BoundaryHandling(object):
self._boundaryFunctions = []
self._nameToIndex = {}
self._boundarySweeps = []
+ self._periodicity = [False, False, False]
self._target = target
if target not in ('cpu', 'gpu'):
raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
+ @property
+ def periodicity(self):
+ return self._periodicity
+
+ def setPeriodicity(self, x=False, y=False, z=False):
+ self._periodicity = [x, y, z]
+ self.invalidateIndexCache()
+
def setBoundary(self, function, indexExpr, maskArr=None, name=None):
"""
Sets boundary in a rectangular region (slice)
@@ -122,10 +133,42 @@ class BoundaryHandling(object):
self._boundarySweeps.append(makePythonGpuFunction(ast, {'indexField': idxGpuField}))
else:
assert False
+ self._addPeriodicityHandlers()
+
+ def _addPeriodicityHandlers(self):
+ 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._boundarySweeps.append(getPeriodicBoundaryFunctor(filteredStencil, ghostLayers=1))
+ elif self._target == 'gpu':
+ from pystencils.gpucuda.periodicity import getPeriodicBoundaryFunctor
+ self._boundarySweeps.append(getPeriodicBoundaryFunctor(filteredStencil, self.flagField.shape,
+ indexDimensions=1,
+ indexDimShape=len(self._lbMethod.stencil),
+ ghostLayers=1))
+ else:
+ assert False
def __call__(self, **kwargs):
- if len(self._boundaryFunctions) == 0:
- return
if len(self._boundarySweeps) == 0:
self.prepare()
for boundarySweep in self._boundarySweeps:
diff --git a/scenarios.py b/scenarios.py
index 2cbd3cf325753ffa6710241041d9a45224d91541..1f07e4706ae9d71492756c5579190a118a4e1b83 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -67,7 +67,7 @@ class Scenario(object):
target=optimizationParams['target'])
self._preUpdateFunctions = preUpdateFunctions
- self._kernelParams = kernelParams
+ self.kernelParams = kernelParams
self._pdfGpuArrays = []
if initialVelocity is None:
@@ -179,7 +179,7 @@ class Scenario(object):
f(pdfArrays[0])
if self._boundaryHandling is not None:
self._boundaryHandling(pdfs=pdfArrays[0])
- self._lbmKernel(src=pdfArrays[0], dst=pdfArrays[1], **self._kernelParams)
+ self._lbmKernel(src=pdfArrays[0], dst=pdfArrays[1], **self.kernelParams)
pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0] # swap
@@ -208,12 +208,10 @@ def createFullyPeriodicFlow(initialVelocity, optimizationParams={}, lbmKernel=No
"""
Creates a fully periodic setup with prescribed velocity field
"""
- domainSize = initialVelocity.shape
-
- stencil = getStencil(kwargs['stencil'])
- periodicity = getPeriodicBoundaryFunctor(stencil)
-
- return Scenario(domainSize, kwargs, optimizationParams, lbmKernel, initialVelocity, kernelParams, [periodicity])
+ domainSize = initialVelocity.shape[:-1]
+ scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel, initialVelocity, kernelParams)
+ scenario.boundaryHandling.setPeriodicity(True, True, True)
+ return scenario
def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None,
@@ -307,16 +305,14 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
else:
roundChannel = False
- def periodicity(pdfArr):
- pdfArr[0, :, :] = pdfArr[-2, :, :]
- pdfArr[-1, :, :] = pdfArr[1, :, :]
+ if 'forceModel' not in kwargs:
+ kwargs['forceModel'] = 'guo'
scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel,
- initialVelocity=initialVelocity, preUpdateFunctions=[periodicity],
- kernelParams=kernelParams)
+ initialVelocity=initialVelocity, kernelParams=kernelParams)
boundaryHandling = scenario.boundaryHandling
-
+ boundaryHandling.setPeriodicity(True, False, False)
if dim == 2:
for direction in ('N', 'S'):
boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
@@ -339,3 +335,21 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
kwargs['forceModel'] = 'guo'
return scenario
+
+if __name__ == '__main__':
+ from lbmpy.scenarios import createForceDrivenChannel
+ from lbmpy.boundaries.geometry import BlackAndWhiteImageBoundary
+ from pystencils.slicing import makeSlice
+ from lbmpy.boundaries import noSlip
+ import numpy as np
+ domainSize = (10, 10)
+ scenario = createForceDrivenChannel(dim=2, method='srt', force=0.000002,
+ domainSize=domainSize,
+ relaxationRates=[1.92], forceModel='guo',
+ compressible=True,
+ optimizationParams={'target': 'gpu'})
+ imageSetup = BlackAndWhiteImageBoundary("/home/staff/bauer/opensource.png",
+ noSlip, targetSlice=makeSlice[2:-2, 1:-2])
+ imageSetup(scenario.boundaryHandling, scenario.method, domainSize)
+ scenario.boundaryHandling.prepare()
+ scenario.run(1)