Bugfixes in boundary data handling for GPU & new rod scenario

boundaries/handlinginterface.py

@@ -53,7 +53,7 @@ class BoundaryDataSetter:
 
     def fluidCellPositions(self, coord):
         assert coord < self.dim
-        return self.indexArray[self.coordMap[coord]] + self.offset[coord] - self.ghostLayers
+        return self.indexArray[self.coordMap[coord]] + self.offset[coord] - self.ghostLayers + 0.5
 
     @memorycache()
     def linkOffsets(self):
@@ -150,9 +150,6 @@ class GenericBoundaryHandling(object):
         else:
             for boundaryObject, boundaryInfo in self._boundaryInfos.items():
                 self.__boundaryDataInitialization(boundaryInfo, boundaryObject, **kwargs)
-                if self._target == 'gpu':
-                    import pycuda.gpuarray as gpuarray
-                    boundaryInfo.idxArrayForExecution = gpuarray.to_gpu(boundaryInfo.indexArray)
 
     def prepare(self):
         """Compiles boundary kernels according to flag field. When setting boundaries the cache is automatically
@@ -297,3 +294,9 @@ class GenericBoundaryHandling(object):
             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

geometry.py

@@ -7,17 +7,20 @@ def addParabolicVelocityInflow(boundaryHandling, u_max, indexExpr, velCoord=0, d
     def velocityInfoCallback(boundaryData):
         for i, name in enumerate(['vel_0', 'vel_1', 'vel_2']):
             if i != velCoord:
-                boundaryData[name] = 0
+                boundaryData[name] = 0.0
         if diameter is None:
             radius = min(sh for i, sh in enumerate(boundaryHandling.domainShape) if i != velCoord) // 2
         else:
             radius = diameter // 2
-        print("radius", radius)
-        y, z = boundaryData.linkPositions(1), boundaryData.linkPositions(2)
-        centeredY = y - radius
-        centeredZ = z - radius
-        distToCenter = np.sqrt(centeredY ** 2 + centeredZ ** 2)
-        boundaryData['vel_%d' % (velCoord,)] = u_max * (1 - distToCenter / radius)
+
+        normalCoord1 = (velCoord + 1) % 3
+        normalCoord2 = (velCoord + 2) % 3
+        y, z = boundaryData.linkPositions(normalCoord1), boundaryData.linkPositions(normalCoord2)
+        centeredNormal1 = y - radius
+        centeredNormal2 = z - radius
+        distToCenter = np.sqrt(centeredNormal1 ** 2 + centeredNormal2 ** 2)
+        velProfile = u_max * (1 - distToCenter / radius)
+        boundaryData['vel_%d' % (velCoord,)] = velProfile
 
     inflow = UBB(velocityInfoCallback, dim=boundaryHandling.dim)
     boundaryHandling.setBoundary(inflow, indexExpr=indexExpr, includeGhostLayers=False)

scenarios.py

@@ -253,7 +253,7 @@ class Scenario(object):
     :param kernelParams: additional parameters passed to the sweep
     """
 
-    def __init__(self, domainSize, methodParameters, optimizationParams, lbmKernel=None,
+    def __init__(self, domainSize, methodParameters, optimizationParams={}, lbmKernel=None,
                  initialVelocity=None, preUpdateFunctions=[], kernelParams={}):
         ghostLayers = 1
         domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
@@ -329,6 +329,19 @@ class Scenario(object):
         mlups = self.numberOfCells / durationOfTimeStep * 1e-6
         return mlups
 
+    def writeVTK(self):
+        from pyevtk.hl import gridToVTK
+        x = np.arange(self.domainSize[0])
+        y = np.arange(self.domainSize[1])
+        z = np.arange(self.domainSize[2])
+        gridToVTK("vtk_%06d" % (self.timeStepsRun,), x, y, z,
+                  cellData={
+                      'velocity_0': self.velocity[..., 0].filled(0.0),
+                      'velocity_1': self.velocity[..., 1].filled(0.0),
+                      'velocity_2': self.velocity[..., 2].filled(0.0),
+                      'density': self.density.filled(0.0),
+                  })
+
     @property
     def numberOfCells(self):
         result = 1