From 7dbbf0d97fe61950b1a641e9ef791ecc78349064 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Tue, 23 Jan 2018 13:37:20 +0100
Subject: [PATCH] Implemented LBStep - works for basic scenarios

---
 boundary_handling.py | 22 ++++++++++++++--------
 lbstep.py            |  8 +++++---
 scenarios.py         | 36 +++++++++++++++++++++++++-----------
 3 files changed, 44 insertions(+), 22 deletions(-)

diff --git a/boundary_handling.py b/boundary_handling.py
index 81b41d64..07e0c29a 100644
--- a/boundary_handling.py
+++ b/boundary_handling.py
@@ -145,11 +145,10 @@ class BoundaryHandling:
         if self._dirty:
             self.prepare()
 
-        for b in self._dataHandling.iterate():
-            for bInfo in self._boundaryObjectToBoundaryInfo.values():
-                idxArr = b[self._indexArrayName].boundaryObjectToIndexList[bInfo.boundaryObject]
+        for b in self._dataHandling.iterate(gpu=self._target == 'gpu'):
+            for bObj, idxArr in b[self._indexArrayName].boundaryObjectToIndexList.items():
                 kwargs[self._pdfFieldName] = b[self._pdfFieldName]
-                bInfo.kernel(indexField=idxArr, **kwargs)
+                self._boundaryObjectToBoundaryInfo[bObj].kernel(indexField=idxArr, **kwargs)
 
     # ------------------------------ Implementation Details ------------------------------------------------------------
 
@@ -173,14 +172,17 @@ class BoundaryHandling:
         ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
         for b in dh.iterate(ghostLayers=ffGhostLayers):
             flagArr = b[self._flagFieldName]
+            pdfArr = b[self._pdfFieldName]
+            indexArrayBD = b[self._indexArrayName]
+            indexArrayBD.clear()
             for bInfo in self._boundaryObjectToBoundaryInfo.values():
                 idxArr = createBoundaryIndexArray(flagArr, self.lbMethod.stencil, bInfo.flag, self._fluidFlag,
                                                   bInfo.boundaryObject, dh.ghostLayersOfField(self._flagFieldName))
+                if idxArr.size == 0:
+                    continue
 
-                pdfArr = b[self._pdfFieldName]
-                # TODO test that offset is used correctly here
+                # TODO test boundary data setter (especially offset)
                 boundaryDataSetter = BoundaryDataSetter(idxArr, b.offset, self.lbMethod.stencil, ffGhostLayers, pdfArr)
-                indexArrayBD = b[self._indexArrayName]
                 indexArrayBD.boundaryObjectToIndexList[bInfo.boundaryObject] = idxArr
                 indexArrayBD.boundaryObjectToDataSetter[bInfo.boundaryObject] = boundaryDataSetter
                 self._boundaryDataInitialization(bInfo.boundaryObject, boundaryDataSetter)
@@ -198,10 +200,14 @@ class BoundaryHandling:
             self.kernel = kernel
 
     class IndexFieldBlockData:
-        def __init__(self):
+        def __init__(self, *args, **kwargs):
             self.boundaryObjectToIndexList = {}
             self.boundaryObjectToDataSetter = {}
 
+        def clear(self):
+            self.boundaryObjectToIndexList.clear()
+            self.boundaryObjectToDataSetter.clear()
+
         @staticmethod
         def toCpu(gpuVersion, cpuVersion):
             gpuVersion = gpuVersion.boundaryObjectToIndexList
diff --git a/lbstep.py b/lbstep.py
index 01b45ac0..b3a9a55d 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -1,9 +1,10 @@
+import numpy as np
 from lbmpy.boundary_handling import BoundaryHandling
 from lbmpy.creationfunctions import switchToSymbolicRelaxationRatesForEntropicMethods, createLatticeBoltzmannFunction, \
     updateWithDefaultParameters
 from lbmpy.simplificationfactory import createSimplificationStrategy
 from lbmpy.stencils import getStencil
-from pystencils.datahandling import SerialDataHandling
+from pystencils.datahandling.serial_datahandling import SerialDataHandling
 from pystencils import createKernel, makeSlice
 
 
@@ -99,7 +100,7 @@ class LatticeBoltzmannStep:
         if sliceObj is None:
             sliceObj = makeSlice[:, :] if self.dim == 2 else makeSlice[:, :, 0.5]
         for arr in self._dataHandling.gatherArray(dataName, sliceObj):
-            return arr
+            return np.squeeze(arr)
         return None
 
     def velocitySlice(self, sliceObj=None):
@@ -171,4 +172,5 @@ if __name__ == '__main__':
     bh.setBoundary(movingWall, makeSlice[:, -1])
     bh.prepare()
 
-    step.run(5000)
+    step.run(4)
+    step.run(100)
diff --git a/scenarios.py b/scenarios.py
index 2cb5b234..a3f393ec 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -16,7 +16,7 @@ All scenarios can be modified, for example you can create a simple channel first
 >>> scenario.boundaryHandling.setBoundary(NoSlip(), makeSlice[0.3:0.4, 0.0:0.3])
 
 Functions for scenario setup:
------------------------------
+----    -------------------------
 
 All of the following scenario creation functions take keyword arguments specifying which LBM method should be used
 and a ``optimizationParams`` dictionary, defining performance related options. These parameters are documented
@@ -31,7 +31,8 @@ from pystencils.slicing import sliceFromDirection
 from lbmpy.boundaries import NoSlip, UBB, FixedDensity
 
 
-def createFullyPeriodicFlow(initialVelocity, periodicityInKernel=False, lbmKernel=None, parallel=False, **kwargs):
+def createFullyPeriodicFlow(initialVelocity, periodicityInKernel=False, lbmKernel=None,
+                            dataHandling=None, parallel=False, **kwargs):
     """
     Creates a fully periodic setup with prescribed velocity field
 
@@ -51,14 +52,17 @@ def createFullyPeriodicFlow(initialVelocity, periodicityInKernel=False, lbmKerne
     if periodicityInKernel:
         kwargs['optimizationParams']['builtinPeriodicity'] = (True, True, True)
 
-    dataHandling = createDataHandling(parallel, domainSize, periodicity=not periodicityInKernel, defaultGhostLayers=1)
+    if dataHandling is None:
+        dataHandling = createDataHandling(parallel, domainSize, periodicity=not periodicityInKernel,
+                                          defaultGhostLayers=1)
     step = LatticeBoltzmannStep(dataHandling=dataHandling, lbmKernel=lbmKernel, **kwargs)
     for b in step.dataHandling.iterate(ghostLayers=False):
         np.copyto(b[step.velocityDataName], initialVelocity[b.globalSlice])
     return step
 
 
-def createLidDrivenCavity(domainSize, lidVelocity=0.005, lbmKernel=None, parallel=False, **kwargs):
+def createLidDrivenCavity(domainSize=None, lidVelocity=0.005, lbmKernel=None, parallel=False,
+                          dataHandling=None, **kwargs):
     """
     Creates a lid driven cavity scenario
 
@@ -69,7 +73,9 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, lbmKernel=None, paralle
     :param parallel: True for distributed memory parallelization with waLBerla
     :return: instance of :class:`Scenario`
     """
-    dataHandling = createDataHandling(parallel, domainSize, periodicity=False, defaultGhostLayers=1)
+    assert domainSize is not None or dataHandling is not None
+    if dataHandling is None:
+        dataHandling = createDataHandling(parallel, domainSize, periodicity=False, defaultGhostLayers=1)
     step = LatticeBoltzmannStep(dataHandling=dataHandling, lbmKernel=lbmKernel, **kwargs)
 
     myUbb = UBB(velocity=[lidVelocity, 0, 0][:step.method.dim])
@@ -80,8 +86,8 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, lbmKernel=None, paralle
     return step
 
 
-def createChannel(domainSize, force=None, pressureDifference=None, u_max=None,
-                  diameterCallback=None, duct=False, wallBoundary=NoSlip(), parallel=False, **kwargs):
+def createChannel(domainSize, force=None, pressureDifference=None, u_max=None, diameterCallback=None,
+                  duct=False, wallBoundary=NoSlip(), parallel=False, dataHandling=None, **kwargs):
     """
     Create a channel scenario (2D or 3D)
     :param domainSize: size of the simulation domain. First coordinate is the flow direction.
@@ -107,7 +113,8 @@ def createChannel(domainSize, force=None, pressureDifference=None, u_max=None,
         raise ValueError("Please specify exactly one of the parameters 'force', 'pressureDifference' or 'u_max'")
 
     periodicity = (True, False, False) if force else (False, False, False)
-    dataHandling = createDataHandling(parallel, domainSize, periodicity=periodicity[:dim], defaultGhostLayers=1)
+    if dataHandling is None:
+        dataHandling = createDataHandling(parallel, domainSize, periodicity=periodicity[:dim], defaultGhostLayers=1)
     step = LatticeBoltzmannStep(dataHandling=dataHandling, **kwargs)
     boundaryHandling = step.boundaryHandling
     if force:
@@ -136,7 +143,14 @@ def createChannel(domainSize, force=None, pressureDifference=None, u_max=None,
 
 if __name__ == '__main__':
     import pycuda.autoinit
-    ldc = createLidDrivenCavity((30, 30, 30), relaxationRate=1.5,
-                                optimizationParams={'openMP': 2, 'target': 'gpu', 'fieldLayout': 'c'})
+    import waLBerla as wlb
+    from pystencils.datahandling import ParallelDataHandling
+    blocks = wlb.createUniformBlockGrid(blocks=(2, 2, 1), cellsPerBlock=(20, 20, 1), oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks, dim=2)
+    ldc = createLidDrivenCavity(relaxationRate=1.5, dataHandling=dh,
+                                optimizationParams={'openMP': 2, 'target': 'gpu', 'fieldLayout': 'f',
+                                                    'gpuIndexingParams': {'blockSize': (8, 8, 2)}})
     ldc.boundaryHandling.prepare()
-    ldc.run(3)
\ No newline at end of file
+
+    ldc.run(4)
+    ldc.run(100)
-- 
GitLab