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

more CUDA bugfixes & periodicity kernels for CUDA

parent 6a335105
No related branches found
No related tags found
No related merge requests found
import sympy as sp import sympy as sp
import numpy as np import numpy as np
from lbmpy.stencils import getStencil
from pystencils import TypedSymbol, Field from pystencils import TypedSymbol, Field
from pystencils.backends.cbackend import CustomCppCode from pystencils.backends.cbackend import CustomCppCode
from lbmpy.boundaries.createindexlist import createBoundaryIndexList from lbmpy.boundaries.createindexlist import createBoundaryIndexList
...@@ -37,10 +39,19 @@ class BoundaryHandling(object): ...@@ -37,10 +39,19 @@ class BoundaryHandling(object):
self._boundaryFunctions = [] self._boundaryFunctions = []
self._nameToIndex = {} self._nameToIndex = {}
self._boundarySweeps = [] self._boundarySweeps = []
self._periodicity = [False, False, False]
self._target = target self._target = target
if target not in ('cpu', 'gpu'): if target not in ('cpu', 'gpu'):
raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,)) 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): def setBoundary(self, function, indexExpr, maskArr=None, name=None):
""" """
Sets boundary in a rectangular region (slice) Sets boundary in a rectangular region (slice)
...@@ -122,10 +133,42 @@ class BoundaryHandling(object): ...@@ -122,10 +133,42 @@ class BoundaryHandling(object):
self._boundarySweeps.append(makePythonGpuFunction(ast, {'indexField': idxGpuField})) self._boundarySweeps.append(makePythonGpuFunction(ast, {'indexField': idxGpuField}))
else: else:
assert False 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): def __call__(self, **kwargs):
if len(self._boundaryFunctions) == 0:
return
if len(self._boundarySweeps) == 0: if len(self._boundarySweeps) == 0:
self.prepare() self.prepare()
for boundarySweep in self._boundarySweeps: for boundarySweep in self._boundarySweeps:
......
...@@ -67,7 +67,7 @@ class Scenario(object): ...@@ -67,7 +67,7 @@ class Scenario(object):
target=optimizationParams['target']) target=optimizationParams['target'])
self._preUpdateFunctions = preUpdateFunctions self._preUpdateFunctions = preUpdateFunctions
self._kernelParams = kernelParams self.kernelParams = kernelParams
self._pdfGpuArrays = [] self._pdfGpuArrays = []
if initialVelocity is None: if initialVelocity is None:
...@@ -179,7 +179,7 @@ class Scenario(object): ...@@ -179,7 +179,7 @@ class Scenario(object):
f(pdfArrays[0]) f(pdfArrays[0])
if self._boundaryHandling is not None: if self._boundaryHandling is not None:
self._boundaryHandling(pdfs=pdfArrays[0]) 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 pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0] # swap
...@@ -208,12 +208,10 @@ def createFullyPeriodicFlow(initialVelocity, optimizationParams={}, lbmKernel=No ...@@ -208,12 +208,10 @@ def createFullyPeriodicFlow(initialVelocity, optimizationParams={}, lbmKernel=No
""" """
Creates a fully periodic setup with prescribed velocity field Creates a fully periodic setup with prescribed velocity field
""" """
domainSize = initialVelocity.shape domainSize = initialVelocity.shape[:-1]
scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel, initialVelocity, kernelParams)
stencil = getStencil(kwargs['stencil']) scenario.boundaryHandling.setPeriodicity(True, True, True)
periodicity = getPeriodicBoundaryFunctor(stencil) return scenario
return Scenario(domainSize, kwargs, optimizationParams, lbmKernel, initialVelocity, kernelParams, [periodicity])
def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None, def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None,
...@@ -307,16 +305,14 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No ...@@ -307,16 +305,14 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
else: else:
roundChannel = False roundChannel = False
def periodicity(pdfArr): if 'forceModel' not in kwargs:
pdfArr[0, :, :] = pdfArr[-2, :, :] kwargs['forceModel'] = 'guo'
pdfArr[-1, :, :] = pdfArr[1, :, :]
scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel, scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel,
initialVelocity=initialVelocity, preUpdateFunctions=[periodicity], initialVelocity=initialVelocity, kernelParams=kernelParams)
kernelParams=kernelParams)
boundaryHandling = scenario.boundaryHandling boundaryHandling = scenario.boundaryHandling
boundaryHandling.setPeriodicity(True, False, False)
if dim == 2: if dim == 2:
for direction in ('N', 'S'): for direction in ('N', 'S'):
boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim)) boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
...@@ -339,3 +335,21 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No ...@@ -339,3 +335,21 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
kwargs['forceModel'] = 'guo' kwargs['forceModel'] = 'guo'
return scenario 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)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment