kida vortex flow for SRT

creationfunctions.py

@@ -68,7 +68,7 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
         if 'openMP' in optParams:
             if isinstance(optParams['openMP'], bool) and optParams['openMP']:
                 addOpenMP(ast)
-            elif isinstance(optParams['openMP'], int):
+            elif not isinstance(optParams['openMP'], bool) and isinstance(optParams['openMP'], int):
                 addOpenMP(ast, numThreads=optParams['openMP'])
         res = makePythonCpuFunction(ast)
     elif optParams['target'] == 'gpu':

macroscopic_value_kernels.py

@@ -37,7 +37,6 @@ def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fiel
 
         indDims = 0 if numberOfElements <= 1 else 1
         if pdfArr is None:
-            fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape) - 1])
             outputField = Field.createGeneric(outputQuantity, lbMethod.dim, layout=fieldLayout, indexDimensions=indDims)
         else:
             outputFieldShape = pdfArr.shape[:-1]
@@ -142,7 +141,9 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field
     return setter
 
 
-def createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityRelaxationRate=1.3):
+def createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velocityRelaxationRate=0.8):
+
+    lbMethod = collisionRule.method
     velocityField = Field.createFromNumpyArray('velInput', velocityArray, indexDimensions=1)
 
     cqc = lbMethod.conservedQuantityComputation
@@ -162,30 +163,44 @@ def createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityR
     # set first order relaxation rate
     lbMethod = deepcopy(lbMethod)
     lbMethod.setFirstMomentRelaxationRate(velocityRelaxationRate)
+    lbMethod.setZerothMomentRelaxationRate(0)
+
+    simplificationStrategy = createSimplificationStrategy(lbMethod)
+    newCollisionRule = simplificationStrategy(lbMethod.getCollisionRule(eqInput))
+
+    # if the original collision rule used an entropy condition, the initialization should use one as well
+    # the simplification hints contain the entropy condition parameters
+    sh = collisionRule.simplificationHints
+    if 'entropic' in sh and sh['entropic']:
+        iterations = sh['entropicNewtonIterations']
+        if iterations:
+            from lbmpy.methods.entropic import addIterativeEntropyCondition
+            newCollisionRule = addIterativeEntropyCondition(newCollisionRule, newtonIterations=iterations)
+        else:
+            from lbmpy.methods.entropic import addEntropyCondition
+            newCollisionRule = addEntropyCondition(newCollisionRule)
 
-    return lbMethod.getCollisionRule(eqInput)
+    return newCollisionRule
 
 
-def compileAdvancedVelocitySetter(lbMethod, velocityArray, velocityRelaxationRate=1.3, pdfArr=None,
-                                  fieldLayout='numpy', optimizationParameters={}):
+def compileAdvancedVelocitySetter(collisionRule, velocityArray, velocityRelaxationRate=0.8, pdfArr=None,
+                                  fieldLayout='numpy', optimizationParams={}):
     """
     Advanced initialization of velocity field through iteration procedure according to
     Mei, Luo, Lallemand and Humieres: Consistent initial conditions for LBM simulations, 2005
 
-    :param lbMethod:
+    :param collisionRule:
     :param velocityArray: array with velocity field
     :param velocityRelaxationRate: relaxation rate for the velocity moments - determines convergence behaviour
                                    of the initialization scheme
     :param pdfArr: optional numpy array for pdf field - used to get optimal loop structure for kernel
     :param fieldLayout: layout of the pdf field if pdfArr was not given
-    :param optimizationParameters: dictionary with optimization hints
+    :param optimizationParams: dictionary with optimization hints
     :return: stream-collide update function
     """
     from lbmpy.updatekernels import createStreamPullKernel
-    from lbmpy.simplificationfactory import createSimplificationStrategy
     from lbmpy.creationfunctions import createLatticeBoltzmannAst, createLatticeBoltzmannFunction
-    collisionRule = createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityRelaxationRate)
-    simp = createSimplificationStrategy(collisionRule.method)
-    updateRule = createStreamPullKernel(simp(collisionRule), pdfArr, genericLayout=fieldLayout)
-    ast = createLatticeBoltzmannAst(updateRule, optimizationParameters)
-    return createLatticeBoltzmannFunction(ast, optimizationParameters)
+    newCollisionRule = createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velocityRelaxationRate)
+    updateRule = createStreamPullKernel(newCollisionRule, pdfArr, genericLayout=fieldLayout)
+    ast = createLatticeBoltzmannAst(updateRule, optimizationParams)
+    return createLatticeBoltzmannFunction(ast, optimizationParams)

methods/entropic.py

@@ -50,7 +50,10 @@ def addEntropyCondition(collisionRule):
         index = collisionRule.method.postCollisionPdfSymbols.index(updateEq.lhs)
         newEq = sp.Eq(updateEq.lhs, fSymbols[index] + omega_s * dsSymbols[index] + omega_h * dhSymbols[index])
         newUpdateEquations.append(newEq)
-    return collisionRule.copy(newUpdateEquations, collisionRule.subexpressions + subexprs)
+    newCollisionRule = collisionRule.copy(newUpdateEquations, collisionRule.subexpressions + subexprs)
+    newCollisionRule.simplificationHints['entropic'] = True
+    newCollisionRule.simplificationHints['entropicNewtonIterations'] = None
+    return newCollisionRule
 
 
 def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue=1):
@@ -109,7 +112,10 @@ def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue
 
     # final update equations
     newSubexpressions = collisionRule.subexpressions + rrFactorDefinitions + fEqEqs + newtonIterationEquations
-    return collisionRule.copy(newUpdateEquations, newSubexpressions)
+    newCollisionRule = collisionRule.copy(newUpdateEquations, newSubexpressions)
+    newCollisionRule.simplificationHints['entropic'] = True
+    newCollisionRule.simplificationHints['entropicNewtonIterations'] = newtonIterations
+    return newCollisionRule
 
 
 # --------------------------------- Helper Functions and Classes -------------------------------------------------------

methods/momentbased.py

@@ -116,6 +116,12 @@ class MomentBasedLbMethod(AbstractLbMethod):
                                                             True, conservedQuantityEquations)
         return eqColl
 
+    def setZerothMomentRelaxationRate(self, relaxationRate):
+        e = sp.Rational(1, 1)
+        prevEntry = self._momentToRelaxationInfoDict[e]
+        newEntry = RelaxationInfo(prevEntry[0], relaxationRate)
+        self._momentToRelaxationInfoDict[e] = newEntry
+
     def setFirstMomentRelaxationRate(self, relaxationRate):
         for e in MOMENT_SYMBOLS[:self.dim]:
             assert e in self._momentToRelaxationInfoDict, "First moments are not relaxed separately by this method"

serialscenario.py

@@ -8,10 +8,10 @@ from lbmpy.boundaries import BoundaryHandling, noSlip, ubb, fixedDensity
 from lbmpy.stencils import getStencil
 
 
-def createScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, lbmKernel=None,
+def createScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParams, lbmKernel=None,
                    initialVelocity=None, preUpdateFunctions=[]):
-    if 'target' not in optimizationParameters:
-        optimizationParameters['target'] = 'cpu'
+    if 'target' not in optimizationParams:
+        optimizationParams['target'] = 'cpu'
 
     ghostLayers = 1
     domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
@@ -26,7 +26,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
 
     # Create kernel
     if lbmKernel is None:
-        methodParameters['optimizationParams'] = optimizationParameters
+        methodParameters['optimizationParams'] = optimizationParams
         lbmKernel = createLatticeBoltzmannFunction(**methodParameters)
     method = lbmKernel.method
 
@@ -36,7 +36,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
     if boundarySetupFunction is not None:
         symPdfField = Field.createFromNumpyArray('pdfs', pdfArrays[0], indexDimensions=1)
         boundaryHandling = BoundaryHandling(symPdfField, domainSize, lbmKernel.method,
-                                            target=optimizationParameters['target'])
+                                            target=optimizationParams['target'])
         boundarySetupFunction(boundaryHandling=boundaryHandling, method=method)
         boundaryHandling.prepare()
     else:
@@ -53,7 +53,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
                                                     pdfArr=pdfArrays[0], target='cpu')
     setMacroscopic(pdfs=pdfArrays[0])
 
-    if optimizationParameters['target'] == 'gpu':
+    if optimizationParams['target'] == 'gpu':
         import pycuda.gpuarray as gpuarray
         pdfGpuArrays = [gpuarray.to_gpu(a) for a in pdfArrays]
     else:
@@ -95,10 +95,10 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
     cpuTimeLoop.kernel = lbmKernel
     gpuTimeLoop.kernel = lbmKernel
 
-    return gpuTimeLoop if optimizationParameters['target'] == 'gpu' else cpuTimeLoop
+    return gpuTimeLoop if optimizationParams['target'] == 'gpu' else cpuTimeLoop
 
 
-def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, lbmKernel=None, **kwargs):
+def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None, **kwargs):
     def boundarySetupFunction(boundaryHandling, method):
         myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:method.dim])
         myUbb.name = 'ubb'
@@ -106,11 +106,11 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters=
         for direction in ('W', 'E', 'S') if method.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
             boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
 
-    return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel)
+    return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel)
 
 
 def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None,
-                                        lbmKernel=None, optimizationParameters={}, **kwargs):
+                                        lbmKernel=None, optimizationParams={}, **kwargs):
     assert dim in (2, 3)
 
     if radius is not None:
@@ -152,11 +152,11 @@ def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None
     if 'forceModel' not in kwargs:
         kwargs['forceModel'] = 'guo'
 
-    return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel)
+    return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel)
 
 
 def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, lbmKernel=None,
-                             optimizationParameters={}, initialVelocity=None, boundarySetupFunctions=[], **kwargs):
+                             optimizationParams={}, initialVelocity=None, boundarySetupFunctions=[], **kwargs):
     assert dim in (2, 3)
     kwargs['force'] = tuple([force, 0, 0][:dim])
 
@@ -197,6 +197,6 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
     if 'forceModel' not in kwargs:
         kwargs['forceModel'] = 'guo'
 
-    return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel,
+    return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel,
                           initialVelocity=initialVelocity, preUpdateFunctions=[periodicity])