fixed compressible bug

boundaries/boundaryconditions.py

@@ -20,6 +20,11 @@ def ubb(pdfField, direction, lbMethod, velocity):
     neighbor = offsetFromDir(direction, lbMethod.dim)
     inverseDir = invDir(direction)
 
+    # TODO adapt velocity to force
+    # TODO compute density
+
+    densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
+
     velTerm = 6 * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
     return [sp.Eq(pdfField[neighbor](inverseDir),
                   pdfField(direction) - velTerm)]

creationfunctions.py

@@ -26,7 +26,7 @@ def _getParams(params, optParams):
     defaultOptimizationDescription = {
         'doCseInOpposingDirections': False,
         'doOverallCse': False,
-        'split': True,
+        'split': False,
 
         'fieldSize': None,
         'fieldLayout': 'reverseNumpy',  # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
@@ -80,6 +80,7 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
     if params['target'] == 'cpu':
         from pystencils.cpu import createKernel
         if 'splitGroups' in updateRule.simplificationHints:
+            print("splitting!")
             splitGroups = updateRule.simplificationHints['splitGroups']
         else:
             splitGroups = ()
@@ -91,6 +92,11 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
         return ValueError("'target' has to be either 'cpu' or 'gpu'")
 
     res.method = updateRule.method
+    #TODO debug begin
+    from pystencils.cpu import generateC
+    from pystencils.display_utils import highlightCpp
+    print(generateC(res))
+    #TODO debug end
     return res
 
 

methods/conservedquantitycomputation.py

 import abc
+from collections import OrderedDict
+
 import sympy as sp
 from pystencils.equationcollection import EquationCollection
 
@@ -139,8 +141,6 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
     def outputEquationsFromPdfs(self, pdfs, outputQuantityNamesToSymbols):
         dim = len(self._stencil[0])
 
-        symbolsToExtract = set()
-
         eqColl = getEquationsForZerothAndFirstOrderMoment(self._stencil, pdfs, self._symbolOrder0, self._symbolsOrder1)
 
         if self._compressible:
@@ -151,14 +151,22 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
         eqColl = applyForceModelShift('macroscopicVelocityShift', dim, eqColl, self._forceModel, self._compressible)
 
         mainEquations = []
+        eqs = OrderedDict([(eq.lhs, eq.rhs) for eq in eqColl.allEquations])
+
         if 'density' in outputQuantityNamesToSymbols:
             densityOutputSymbol = outputQuantityNamesToSymbols['density']
-            mainEquations.append(sp.Eq(densityOutputSymbol, self._symbolOrder0))
-            symbolsToExtract.add(densityOutputSymbol)
+            if densityOutputSymbol != self._symbolOrder0:
+                mainEquations.append(sp.Eq(densityOutputSymbol, self._symbolOrder0))
+            else:
+                mainEquations.append(sp.Eq(self._symbolOrder0, eqs[self._symbolOrder0]))
+                del eqs[self._symbolOrder0]
         if 'velocity' in outputQuantityNamesToSymbols:
             velOutputSymbols = outputQuantityNamesToSymbols['velocity']
-            mainEquations += [sp.Eq(a, b) for a, b in zip(velOutputSymbols, self._symbolsOrder1)]
-            symbolsToExtract.update(velOutputSymbols)
+            if tuple(velOutputSymbols) != tuple(self._symbolsOrder1):
+                mainEquations += [sp.Eq(a, b) for a, b in zip(velOutputSymbols, self._symbolsOrder1)]
+            else:
+                # TODO
+                pass
 
         eqColl = eqColl.copy(mainEquations, eqColl.allEquations)
         return eqColl.newWithoutUnusedSubexpressions()
@@ -271,3 +279,14 @@ def applyForceModelShift(shiftMemberName, dim, equationCollection, forceModel, c
 
 
 
+if __name__ == '__main__':
+    from lbmpy.creationfunctions import createLatticeBoltzmannMethod
+    from lbmpy.simplificationfactory import createSimplificationStrategy
+    from lbmpy.stencils import getStencil
+    from lbmpy_old.lbmgenerator import createStreamCollideUpdateRule
+    from lbmpy_old.latticemodel import makeSRT
+    import sympy as sp
+    methodNew = createLatticeBoltzmannMethod(compressible=True)
+    newSimp = createSimplificationStrategy(methodNew)
+    cqc = methodNew.conservedQuantityComputation
+    cqc.outputEquationsFromPdfs(sp.symbols("f_:9"), {'density': sp.Symbol("rho_out")})
\ No newline at end of file

serialscenario.py

@@ -18,6 +18,7 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
     # Create kernel
     lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
     method = lbmKernel.method
+
     assert D == method.dim, "Domain size and stencil do not match"
     Q = len(method.stencil)
 
@@ -38,7 +39,7 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
     densityArr = np.zeros(domainSizeWithGhostLayer)
     velocityArr = np.zeros(domainSizeWithGhostLayer + (D,))
     getMacroscopic = compileMacroscopicValuesGetter(method, ['density', 'velocity'], pdfArr=pdfSrc)
-    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0]*D}, pdfArr=pdfSrc)
+    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0] * D}, pdfArr=pdfSrc)
     setMacroscopic(pdfs=pdfSrc)
 
     # Run simulation
@@ -50,6 +51,7 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
             if boundaryHandling is not None:
                 boundaryHandling(pdfs=pdfSrc)
             lbmKernel(src=pdfSrc, dst=pdfDst)
+
             pdfSrc, pdfDst = pdfDst, pdfSrc
         getMacroscopic(pdfs=pdfSrc, density=densityArr, velocity=velocityArr)
         return pdfSrc, densityArr, velocityArr
@@ -64,6 +66,7 @@ def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={},
         boundaryHandling.setBoundary(myUbb, sliceFromDirection('N', method.dim))
         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)
 
 
@@ -74,11 +77,11 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
     if radius is not None:
         assert length is not None
         if dim == 3:
-            domainSize = (length, 2*radius+1, 2*radius+1)
+            domainSize = (length, 2 * radius + 1, 2 * radius + 1)
             roundChannel = True
         else:
             if domainSize is None:
-                domainSize = (length, 2*radius)
+                domainSize = (length, 2 * radius)
     else:
         roundChannel = False
 
@@ -108,5 +111,3 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
 
     return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters,
                        preUpdateFunctions=[periodicity])
-
-