Pressure driven channel

boundaries/boundaryconditions.py

@@ -52,7 +52,11 @@ def fixedDensity(pdfField, direction, lbMethod, density):
     neighbor = offsetFromDir(direction, lbMethod.dim)
     inverseDir = invDir(direction)
 
-    symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium())
+    velocity = lbMethod.conservedQuantityComputation.definedSymbols()['velocity']
+    symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
+    substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
+    symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
+
     simplification = createSimplificationStrategy(lbMethod)
     symmetricEq = simplification(symmetricEq)
 

creationfunctions.py

@@ -34,6 +34,7 @@ def _getParams(params, optParams):
         'fieldLayout': 'reverseNumpy',  # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
 
         'openMP': True,
+        'pdfArr': None,
     }
     unknownParams = [k for k in params.keys() if k not in defaultMethodDescription]
     unknownOptParams = [k for k in optParams.keys() if k not in defaultOptimizationDescription]
@@ -70,6 +71,7 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
         return ValueError("'target' has to be either 'cpu' or 'gpu'")
 
     res.method = ast.method
+    res.ast = ast
     return res
 
 
@@ -116,11 +118,14 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
         npField = createPdfArray(optParams['fieldSize'], len(stencil), layout=optParams['fieldLayout'])
         updateRule = createStreamPullKernel(collisionRule, numpyField=npField)
     else:
-        layoutName = optParams['fieldLayout']
-        if layoutName == 'fzyx' or 'zyxf':
-            dim = len(stencil[0])
-            layoutName = tuple(reversed(range(dim)))
-        updateRule = createStreamPullKernel(collisionRule, genericLayout=layoutName)
+        if 'pdfArr' in optParams:
+            updateRule = createStreamPullKernel(collisionRule, numpyField=optParams['pdfArr'])
+        else:
+            layoutName = optParams['fieldLayout']
+            if layoutName == 'fzyx' or 'zyxf':
+                dim = len(stencil[0])
+                layoutName = tuple(reversed(range(dim)))
+            updateRule = createStreamPullKernel(collisionRule, genericLayout=layoutName)
 
     return updateRule
 

macroscopicValueKernels.py

 import sympy as sp
 from copy import deepcopy
-from pystencils.field import Field
+from pystencils.field import Field, getLayoutFromNumpyArray
 from lbmpy.simplificationfactory import createSimplificationStrategy
 
 
@@ -30,6 +30,7 @@ def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fiel
         pdfField = Field.createGeneric('pdfs', lbMethod.dim, indexDimensions=1, layout=fieldLayout)
     else:
         pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
+        fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape)-1])
 
     outputMapping = {}
     for outputQuantity in outputQuantities:

methods/creationfunctions.py

@@ -232,7 +232,7 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, useContinuousMaxwell
 # ----------------------------------------- Comparison view for notebooks ----------------------------------------------
 
 
-def compareMomentBasedLbMethods(reference, other):
+def compareMomentBasedLbMethods(reference, other, showDeviationsOnly=False):
     import ipy_table
     table = []
     captionRows = [len(table)]
@@ -244,10 +244,13 @@ def compareMomentBasedLbMethods(reference, other):
         referenceValue = reference.momentToRelaxationInfoDict[moment].equilibriumValue
         otherValue = other.momentToRelaxationInfoDict[moment].equilibriumValue
         diff = sp.simplify(referenceValue - otherValue)
-        table.append(["$%s$" % (sp.latex(moment), ),
-                      "$%s$" % (sp.latex(referenceValue), ),
-                      "$%s$" % (sp.latex(otherValue), ),
-                      "$%s$" % (sp.latex(diff),)])
+        if showDeviationsOnly and diff == 0:
+            pass
+        else:
+            table.append(["$%s$" % (sp.latex(moment), ),
+                          "$%s$" % (sp.latex(referenceValue), ),
+                          "$%s$" % (sp.latex(otherValue), ),
+                          "$%s$" % (sp.latex(diff),)])
 
     onlyInRef = referenceMoments - otherMoments
     if onlyInRef:

methods/entropic.py

@@ -192,7 +192,6 @@ def determineRelaxationRateByEntropyConditionIterative(updateRule, omega_s, omeg
 
 
 def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False, velocityRelaxation=None,
-                                   shearRelaxationRate=sp.Symbol("omega"),
                                    higherOrderRelaxationRate=sp.Symbol("omega_h"),
                                    fixedOmega=None, **kwargs):
     def product(iterable):
@@ -230,7 +229,7 @@ def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False
     else:
         raise ValueError("Unknown model. Supported models KBC-Nx")
 
-    omega_s, omega_h = shearRelaxationRate, higherOrderRelaxationRate
+    omega_s, omega_h = sp.Symbol("omega"), higherOrderRelaxationRate
     shearPart, restPart = decomposition
 
     velRelaxation = omega_s if velocityRelaxation is None else velocityRelaxation
@@ -253,7 +252,7 @@ def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False
         collisionRule = determineRelaxationRateByEntropyCondition(collisionRule, omega_s, omega_h)
 
     if fixedOmega:
-        collisionRule = collisionRule.newWithSubstitutions({omega_s: fixedOmega})
+        collisionRule = collisionRule.copyWithSubstitutionsApplied({omega_s: fixedOmega})
 
     return collisionRule
 

methods/momentbased.py

@@ -8,7 +8,7 @@ from pystencils.sympyextensions import replaceAdditive
 
 
 class MomentBasedLbMethod(AbstractLbMethod):
-    def __init__(self, stencil, momentToRelaxationInfoDict, conservedQuantityComputation, forceModel=None):
+    def __init__(self, stencil, momentToRelaxationInfoDict, conservedQuantityComputation=None, forceModel=None):
         """
         Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
         These methods work by transforming the pdfs into moment space using a linear transformation. In the moment

serialscenario.py

@@ -4,10 +4,11 @@ from pystencils import Field
 from pystencils.slicing import sliceFromDirection
 from lbmpy.creationfunctions import createLatticeBoltzmannFunction
 from lbmpy.macroscopicValueKernels import compileMacroscopicValuesGetter, compileMacroscopicValuesSetter
-from lbmpy.boundaries import BoundaryHandling, noSlip, ubb
+from lbmpy.boundaries import BoundaryHandling, noSlip, ubb, fixedDensity
 
 
-def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, preUpdateFunctions=[]):
+def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, lbmKernel=None,
+                preUpdateFunctions=[]):
     ghostLayers = 1
     domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
     D = len(domainSize)
@@ -16,7 +17,8 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
         methodParameters['stencil'] = 'D2Q9' if D == 2 else 'D3Q27'
 
     # Create kernel
-    lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
+    if lbmKernel is None:
+        lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
     method = lbmKernel.method
 
     assert D == method.dim, "Domain size and stencil do not match"
@@ -56,10 +58,12 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
         getMacroscopic(pdfs=pdfSrc, density=densityArr, velocity=velocityArr)
         return pdfSrc, densityArr, velocityArr
 
+    timeLoop.kernel = lbmKernel
+
     return timeLoop
 
 
-def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, **kwargs):
+def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, lbmKernel=None, **kwargs):
     def boundarySetupFunction(boundaryHandling, method):
         myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:method.dim])
         myUbb.name = 'ubb'
@@ -67,10 +71,57 @@ def runLidDrivenCavity(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 runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters)
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel)
 
 
-def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, optimizationParameters={}, **kwargs):
+def runPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None, lbmKernel=None,
+                                     optimizationParameters={}, **kwargs):
+    assert dim in (2, 3)
+
+    if radius is not None:
+        assert length is not None
+        if dim == 3:
+            domainSize = (length, 2 * radius + 1, 2 * radius + 1)
+            roundChannel = True
+        else:
+            if domainSize is None:
+                domainSize = (length, 2 * radius)
+    else:
+        roundChannel = False
+
+    def boundarySetupFunction(boundaryHandling, method):
+        pressureBoundaryInflow = partial(fixedDensity, density=1.0 + pressureDifference)
+        pressureBoundaryInflow.__name__ = "Inflow"
+
+        pressureBoundaryOutflow = partial(fixedDensity, density=1.0)
+        pressureBoundaryOutflow.__name__ = "Outflow"
+        boundaryHandling.setBoundary(pressureBoundaryInflow, sliceFromDirection('W', method.dim))
+        boundaryHandling.setBoundary(pressureBoundaryOutflow, sliceFromDirection('E', method.dim))
+
+        if method.dim == 2:
+            for direction in ('N', 'S'):
+                boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+        elif method.dim == 3:
+            if roundChannel:
+                noSlipIdx = boundaryHandling.addBoundary(noSlip)
+                ff = boundaryHandling.flagField
+                yMid = ff.shape[1] // 2
+                zMid = ff.shape[2] // 2
+                y, z = np.meshgrid(range(ff.shape[1]), range(ff.shape[2]))
+                ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
+            else:
+                for direction in ('N', 'S', 'T', 'B'):
+                    boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+
+    assert domainSize is not None
+    if 'forceModel' not in kwargs:
+        kwargs['forceModel'] = 'guo'
+
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel)
+
+
+def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, lbmKernel=None,
+                          optimizationParameters={}, **kwargs):
     assert dim in (2, 3)
     kwargs['force'] = tuple([force, 0, 0][:dim])
 
@@ -109,5 +160,16 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
     if 'forceModel' not in kwargs:
         kwargs['forceModel'] = 'guo'
 
-    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters,
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel,
                        preUpdateFunctions=[periodicity])
+
+if __name__ == '__main__':
+    import sympy as sp
+    from pystencils.display_utils import highlightCpp
+    from pystencils.cpu.cpujit import generateC
+    from lbmpy.serialscenario import runPressureGradientDrivenChannel
+    import lbmpy.plot2d as plt
+    timeloop = runPressureGradientDrivenChannel(radius=10, length=30, pressureDifference=0.001,
+                                                relaxationRates=[1.9],
+                                                dim=2)
+    pdfs, rho, vel = timeloop(20)

simplificationfactory.py

@@ -30,11 +30,12 @@ def createSimplificationStrategy(lbmMethod, doCseInOpposingDirections=False, doO
         if splitInnerLoop:
             s.add(createLbmSplitGroups)
 
+    s.add(addSubexpressionsForDivisions)
+
     if doCseInOpposingDirections:
         s.add(cseInOpposingDirections)
     if doOverallCse:
         s.add(sympyCSE)
 
-    s.add(addSubexpressionsForDivisions)
 
     return s

stencils.py

@@ -114,6 +114,21 @@ def stencilsHaveSameEntries(s1, s2):
 
 # -------------------------------------- Visualization -----------------------------------------------------------------
 
+def visualizeStencil(stencil, **kwargs):
+    dim = len(stencil[0])
+    if dim == 2:
+        visualizeStencil2D(stencil, **kwargs)
+    else:
+        slicing = False
+        if 'slice' in kwargs:
+            slicing = kwargs['slice']
+            del kwargs['slice']
+
+        if slicing:
+            visualizeStencil3DBySlicing(stencil, **kwargs)
+        else:
+            visualizeStencil3D(stencil, **kwargs)
+
 
 def visualizeStencil2D(stencil, axes=None, data=None):
     """