diff --git a/geometry.py b/geometry.py new file mode 100644 index 0000000000000000000000000000000000000000..0b9f25fe9b97f71765cb4976b48af3f024168fe0 --- /dev/null +++ b/geometry.py @@ -0,0 +1,35 @@ +import numpy as np +from lbmpy.boundaries import NoSlip + + +def addPipe(boundaryHandling, diameter, boundary=NoSlip()): + """ + Sets boundary for the wall of a pipe with flow in x direction. + + :param boundaryHandling: boundary handling object, works for serial and parallel versions + :param diameter: pipe diameter, can be either a constant value or a callback function. + the callback function has the signature (xCoordArray, domainShapeInCells) and has to return + 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 + dim = len(domainShape) + assert dim in (2, 3) + + def callback(*coordinates): + flowCoord = coordinates[0] + + if callable(diameter): + flowCoordLine = flowCoord[:, 0, 0] if dim == 3 else flowCoord[:, 0] + diameterValue = diameter(flowCoordLine, domainShape) + diameterValue = diameterValue[:, np.newaxis, np.newaxis] if dim == 3 else diameterValue[:, np.newaxis] + else: + diameterValue = diameter + + radiusSq = (diameterValue / 2) ** 2 + + mid = [domainShape[i] // 2 for i in range(1, dim)] + distance = sum((c_i - mid_i) ** 2 for c_i, mid_i in zip(coordinates[1:], mid)) + return distance > radiusSq + + boundaryHandling.setBoundary(boundary, maskCallback=callback) diff --git a/parallel/boundaryhandling.py b/parallel/boundaryhandling.py index 3e59ec0abafb5f926bca4fffcb92f79be097348c..b33af206efdd0116af930553637be36b4cbcfd34 100644 --- a/parallel/boundaryhandling.py +++ b/parallel/boundaryhandling.py @@ -1,8 +1,6 @@ -import numpy as np import waLBerla as wlb from lbmpy.boundaries.handlinginterface import GenericBoundaryHandling, FlagFieldInterface from lbmpy.parallel.blockiteration import slicedBlockIteration -from pystencils.slicing import normalizeSlice class BoundaryHandling(object): @@ -19,6 +17,7 @@ class BoundaryHandling(object): self._pdfFieldId = pdfFieldId self._blocks = blocks self.dim = lbMethod.dim + self.domainShape = [i + 1 for i in self._blocks.getDomainCellBB().max] blocks.addBlockData(boundaryId, addBoundaryHandling) def __call__(self, **kwargs): @@ -36,7 +35,7 @@ class BoundaryHandling(object): def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, sliceNormalizationGhostLayers=1): if indexExpr is None: - indexExpr = [slice(None, None, None)] * self.dim + indexExpr = [slice(None, None, None)] * 3 if len(self._blocks) == 0: return @@ -50,29 +49,6 @@ class BoundaryHandling(object): mask = maskCallback(*indexArrs) if maskCallback else None block[self._boundaryId].setBoundary(boundaryObject, indexExpr=localSlice, maskArr=mask) - #domainCellBB = self._blocks.getDomainCellBB() - #domainExtent = [s + 2 * sliceNormalizationGhostLayers for s in domainCellBB.size] - #indexExpr = normalizeSlice(indexExpr, domainExtent) - #targetCellBB = wlb.CellInterval.fromSlice(indexExpr) - #targetCellBB.shift(*[a - sliceNormalizationGhostLayers for a in domainCellBB.min]) - # - #for block in self._blocks: - # boundaryHandling = block[self._boundaryId] - # ghostLayers = boundaryHandling.ghostLayers - # intersection = self._blocks.getBlockCellBB(block).getExpanded(ghostLayers) - # intersection.intersect(targetCellBB) - # if not intersection.empty(): - # if maskCallback is not None: - # meshGridParams = [offset + 0.5 + np.arange(width) - # for offset, width in zip(intersection.min, intersection.size)] - # indexArr = np.meshgrid(*meshGridParams, indexing='ij') - # mask = maskCallback(*indexArr) - # else: - # mask = None - # localTargetBB = self._blocks.transformGlobalToLocal(block, intersection) - # localTargetBB.shift(ghostLayers, ghostLayers, ghostLayers) - # block[self._boundaryId].setBoundary(boundaryObject, indexExpr=localTargetBB.toSlice(), maskArr=mask) - # ---------------------------------------------------------------------------------------------------------------------- @@ -115,7 +91,7 @@ class WalberlaFlagFieldInterface(FlagFieldInterface): raise NotImplementedError() -if __name__ == '__main__': +def test(): from lbmpy.creationfunctions import createLatticeBoltzmannMethod from lbmpy.boundaries.boundaryconditions import NoSlip from pystencils.slicing import makeSlice @@ -129,7 +105,6 @@ if __name__ == '__main__': bh = BoundaryHandling(blocks, lbMethod, 'pdfField', 'flagField') - def maskCallback(x, y, z): midPoint = (4, 4, 4) radius = 2.0 @@ -139,7 +114,6 @@ if __name__ == '__main__': (z - midPoint[2]) ** 2 <= radius ** 2 return resultMask - bh.setBoundary(NoSlip(), makeSlice[:, :, :], maskCallback, sliceNormalizationGhostLayers=0) vtkOutput = wlb.vtk.makeOutput(blocks, "vtkOutput") flagFieldWriter = wlb.field.createVTKWriter(blocks, "flagField") diff --git a/parallel/scenarios.py b/parallel/scenarios.py index 5c88699c578cff99f9debbe73b7d8c06e5afe7ba..c735c4842adea0f68018ad0c1acea417b678673f 100644 --- a/parallel/scenarios.py +++ b/parallel/scenarios.py @@ -9,6 +9,7 @@ from lbmpy.parallel.blockiteration import slicedBlockIteration class Scenario(object): + vtkCounter = 0 def __init__(self, blocks, methodParameters, optimizationParams, lbmKernel=None, preUpdateFunctions=[], kernelParams={}, blockDataPrefix='', directCommunication=False): @@ -17,6 +18,10 @@ class Scenario(object): self._blockDataPrefix = blockDataPrefix self.kernelParams = kernelParams + if 'stencil' not in methodParameters: + domainShape = [i + 1 for i in blocks.getDomainCellBB().max] + methodParameters['stencil'] = 'D2Q9' if domainShape[2] == 1 else 'D3Q27' + methodParameters, optimizationParams = updateWithDefaultParameters(methodParameters, optimizationParams) self._target = optimizationParams['target'] if optimizationParams['fieldLayout'] in ('f', 'reverseNumpy'): @@ -36,7 +41,7 @@ class Scenario(object): switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams) methodParameters['optimizationParams'] = optimizationParams self._lbmKernel = createLatticeBoltzmannFunction(**methodParameters) - else:for + else: self._lbmKernel = lbmKernel # ----- Add fields @@ -78,7 +83,8 @@ class Scenario(object): self._macroscopicValueSetter = None self._macroscopicValueGetter = None - self._vtkOutput = wlb.vtk.makeOutput(blocks, "vtk") + self._vtkOutput = wlb.vtk.makeOutput(blocks, "vtk_%02d" % (Scenario.vtkCounter,)) + Scenario.vtkCounter += 1 self._vtkOutput.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, self.flagFieldId)) self._vtkOutput.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, self.densityFieldId)) self._vtkOutput.addCellDataWriter(wlb.field.createVTKWriter(self.blocks, self.velocityFieldId))