SerialScenarios

- square channel benchmark for evaluating spurious currents in Poseuille flow
- force model support in cumulant methods
- option to set initial velocity field for simple serial scenarios
- bugfix in simple force model

cache.py

@@ -16,4 +16,4 @@ except ImportError:
 import sys
 calledBySphinx = 'sphinx' in sys.modules
 if calledBySphinx:
-    diskcache = memorycache(maxsize=64)
\ No newline at end of file
+    diskcache = memorycache(maxsize=64)

creationfunctions.py

@@ -17,12 +17,13 @@ def _getParams(params, optParams):
         'method': 'srt',  # can be srt, trt or mrt
         'relaxationRates': sp.symbols("omega_:10"),
         'compressible': False,
-        'order': 2,
+        'equilibriumAccuracyOrder': 2,
 
         'useContinuousMaxwellianEquilibrium': False,
         'cumulant': False,
         'forceModel': 'none',  # can be 'simple', 'luo' or 'guo'
         'force': (0, 0, 0),
+        'initialVelocity': None,
     }
 
     defaultOptimizationDescription = {
@@ -39,7 +40,7 @@ def _getParams(params, optParams):
     unknownParams = [k for k in params.keys() if k not in defaultMethodDescription]
     unknownOptParams = [k for k in optParams.keys() if k not in defaultOptimizationDescription]
     if unknownParams:
-        raise ValueError("Unknown parameter(s): " + ",".join(unknownParams))
+        raise ValueError("Unknown parameter(s): " + ", ".join(unknownParams))
     if unknownOptParams:
         raise ValueError("Unknown optimization parameter(s): " + ",".join(unknownOptParams))
 
@@ -155,7 +156,7 @@ def createLatticeBoltzmannMethod(**params):
 
     commonParams = {
         'compressible': params['compressible'],
-        'equilibriumAccuracyOrder': params['order'],
+        'equilibriumAccuracyOrder': params['equilibriumAccuracyOrder'],
         'forceModel': forceModel,
         'useContinuousMaxwellianEquilibrium': params['useContinuousMaxwellianEquilibrium'],
         'cumulant': params['cumulant'],

forcemodels.py

@@ -25,6 +25,9 @@ class Simple(object):
         return [3 * w_i * scalarProduct(self._force, direction)
                 for direction, w_i in zip(lbMethod.stencil, lbMethod.weights)]
 
+    def macroscopicVelocityShift(self, density):
+        return defaultVelocityShift(density, self._force)
+
 
 class Luo(object):
     r"""

methods/cumulantbased.py

@@ -26,6 +26,14 @@ class CumulantBasedLbMethod(AbstractLbMethod):
     def relaxationInfoDict(self):
         return self._cumulantToRelaxationInfoDict
 
+    @property
+    def zerothOrderEquilibriumMomentSymbol(self, ):
+        return self._conservedQuantityComputation.zerothOrderMomentSymbol
+
+    @property
+    def firstOrderEquilibriumMomentSymbols(self, ):
+        return self._conservedQuantityComputation.firstOrderMomentSymbols
+
     def setFirstMomentRelaxationRate(self, relaxationRate):
         for e in MOMENT_SYMBOLS[:self.dim]:
             assert e in self._cumulantToRelaxationInfoDict, "First cumulants are not relaxed separately by this method"
@@ -84,7 +92,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
 
     def getEquilibrium(self, conservedQuantityEquations=None):
         D = sp.eye(len(self.relaxationRates))
-        return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations, False, False, False)
+        return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations, False, False, False, False)
 
     def getCollisionRule(self, conservedQuantityEquations=None, momentSubexpressions=False,
                          preCollisionSubexpressions=True, postCollisionSubexpressions=False):
@@ -109,7 +117,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
 
     def _getCollisionRuleWithRelaxationMatrix(self, relaxationMatrix, conservedQuantityEquations=None,
                                               momentSubexpressions=False, preCollisionSubexpressions=True,
-                                              postCollisionSubexpressions=False):
+                                              postCollisionSubexpressions=False, includeForceTerms=True):
         def tupleToSymbol(exp, prefix):
             dim = len(exp)
             formatString = prefix + "_" + "_".join(["%d"]*dim)
@@ -177,6 +185,16 @@ class CumulantBasedLbMethod(AbstractLbMethod):
         result = momentTransformationMatrix.inv() * sp.Matrix(collidedMoments)
         mainEquations = [sp.Eq(sym, val) for sym, val in zip(self.postCollisionPdfSymbols, result)]
 
+        # 6) Add forcing terms
+        if self._forceModel is not None and includeForceTerms:
+            forceModelTerms = self._forceModel(self)
+            forceTermSymbols = sp.symbols("forceTerm_:%d" % (len(forceModelTerms,)))
+            forceSubexpressions = [sp.Eq(sym, forceModelTerm)
+                                   for sym, forceModelTerm in zip(forceTermSymbols, forceModelTerms)]
+            subexpressions += forceSubexpressions
+            mainEquations = [sp.Eq(eq.lhs, eq.rhs + forceTermSymbol)
+                             for eq, forceTermSymbol in zip(mainEquations, forceTermSymbols)]
+
         return LbmCollisionRule(self, mainEquations, subexpressions, simplificationHints={})
 
 

serialscenario.py

@@ -9,7 +9,7 @@ from lbmpy.stencils import getStencil
 
 
 def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, lbmKernel=None,
-                preUpdateFunctions=[]):
+                initialVelocity=None, preUpdateFunctions=[]):
     ghostLayers = 1
     domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
     D = len(domainSize)
@@ -42,7 +42,11 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
     densityArr = [np.zeros(domainSizeWithGhostLayer)]
     velocityArr = [np.zeros(domainSizeWithGhostLayer + (D,))]
     getMacroscopic = compileMacroscopicValuesGetter(method, ['density', 'velocity'], pdfArr=pdfArrays[0])
-    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0] * D}, pdfArr=pdfArrays[0])
+
+    if initialVelocity is None:
+        initialVelocity = [0] * D
+    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': initialVelocity},
+                                                    pdfArr=pdfArrays[0])
     setMacroscopic(pdfs=pdfArrays[0])
 
     # Run simulation
@@ -121,7 +125,7 @@ def runPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, r
 
 
 def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, lbmKernel=None,
-                          optimizationParameters={}, **kwargs):
+                          optimizationParameters={}, initialVelocity=None, **kwargs):
     assert dim in (2, 3)
     kwargs['force'] = tuple([force, 0, 0][:dim])
 
@@ -161,5 +165,5 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
         kwargs['forceModel'] = 'guo'
 
     return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel,
-                       preUpdateFunctions=[periodicity])
+                       initialVelocity=initialVelocity, preUpdateFunctions=[periodicity])