From fd57b76555b3ffe1490eb60cda258de23180182f Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Wed, 25 Jan 2017 12:09:23 +0100
Subject: [PATCH] new lbm module: macroscopic values & bugfixes
---
boundaries/boundaryconditions.py | 22 +++++-----
boundaries/boundaryhandling.py | 22 +++++-----
boundaries/createindexlist.py | 5 ++-
boundaries/createindexlistcython.pyx | 1 +
creationfunctions.py | 37 +++++++++++------
macroscopicValueKernels.py | 62 +++++++++++++++++-----------
methods/momentbased.py | 55 +++++++++++++-----------
updatekernels.py | 3 +-
8 files changed, 122 insertions(+), 85 deletions(-)
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 293f3ad9..35b3709e 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -5,19 +5,19 @@ from pystencils.sympyextensions import getSymmetricPart
from lbmpy.boundaries.boundaryhandling import offsetFromDir, weightOfDirection, invDir
-def noSlip(pdfField, direction, lbmMethod):
+def noSlip(pdfField, direction, lbMethod):
"""No-Slip, simple bounce back boundary condition, enforcing zero velocity at obstacle"""
- neighbor = offsetFromDir(direction, lbmMethod.dim)
+ neighbor = offsetFromDir(direction, lbMethod.dim)
inverseDir = invDir(direction)
return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(direction))]
-def ubb(pdfField, direction, lbmMethod, velocity):
+def ubb(pdfField, direction, lbMethod, velocity):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
- assert len(velocity) == lbmMethod.dim, \
- "Dimension of velocity (%d) does not match dimension of LB method (%d)" % (len(velocity, lbmMethod.dim))
- neighbor = offsetFromDir(direction, lbmMethod.dim)
+ assert len(velocity) == lbMethod.dim, \
+ "Dimension of velocity (%d) does not match dimension of LB method (%d)" % (len(velocity), lbMethod.dim)
+ neighbor = offsetFromDir(direction, lbMethod.dim)
inverseDir = invDir(direction)
velTerm = 6 * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
@@ -25,21 +25,21 @@ def ubb(pdfField, direction, lbmMethod, velocity):
pdfField(direction) - velTerm)]
-def fixedDensity(pdfField, direction, lbmMethod, density):
+def fixedDensity(pdfField, direction, lbMethod, density):
"""Boundary condition that fixes the density/pressure at the obstacle"""
def removeAsymmetricPartOfMainEquations(eqColl, dofs):
newMainEquations = [sp.Eq(e.lhs, getSymmetricPart(e.rhs, dofs)) for e in eqColl.mainEquations]
return eqColl.copy(newMainEquations)
- neighbor = offsetFromDir(direction, lbmMethod.dim)
+ neighbor = offsetFromDir(direction, lbMethod.dim)
inverseDir = invDir(direction)
- symmetricEq = removeAsymmetricPartOfMainEquations(lbmMethod.getEquilibrium())
- simplification = createSimplificationStrategy(lbmMethod)
+ symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium())
+ simplification = createSimplificationStrategy(lbMethod)
symmetricEq = simplification(symmetricEq)
- densitySymbol = lbmMethod.conservedQuantityComputation.definedSymbols()['density']
+ densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
conditions = [(eq_i.rhs, sp.Equality(direction, i))
for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index cc2dbb9b..2f6ecd27 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -12,13 +12,13 @@ WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
class BoundaryHandling:
- def __init__(self, symbolicPdfField, domainShape, latticeModel, ghostLayers=1, target='cpu'):
+ def __init__(self, symbolicPdfField, domainShape, lbMethod, ghostLayers=1, target='cpu'):
self._symbolicPdfField = symbolicPdfField
self._shapeWithGhostLayers = [d + 2 * ghostLayers for d in domainShape]
- self._fluidFlag = 2 ** 31
+ self._fluidFlag = 2 ** 30
self.flagField = np.full(self._shapeWithGhostLayers, self._fluidFlag, dtype=np.int32)
self._ghostLayers = ghostLayers
- self._latticeModel = latticeModel
+ self._lbMethod = lbMethod
self._boundaryFunctions = []
self._nameToIndex = {}
self._boundarySweeps = []
@@ -30,8 +30,10 @@ class BoundaryHandling:
if name is None:
name = boundaryFunction.__name__
- self._nameToIndex[name] = len(self._boundaryFunctions)
+ newIdx = len(self._boundaryFunctions)
+ self._nameToIndex[name] = newIdx
self._boundaryFunctions.append(boundaryFunction)
+ return 2 ** newIdx
def invalidateIndexCache(self):
self._boundarySweeps = []
@@ -68,9 +70,9 @@ class BoundaryHandling:
def prepare(self):
self.invalidateIndexCache()
for boundaryIdx, boundaryFunc in enumerate(self._boundaryFunctions):
- idxField = createBoundaryIndexList(self.flagField, self._latticeModel.stencil,
+ idxField = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
2 ** boundaryIdx, self._fluidFlag, self._ghostLayers)
- ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._latticeModel, boundaryFunc)
+ ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._lbMethod, boundaryFunc)
if self._target == 'cpu':
from pystencils.cpu import makePythonFunction as makePythonCpuFunction
@@ -131,8 +133,8 @@ class LbmMethodInfo(CustomCppCode):
super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
-def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
- dim = latticeModel.dim
+def generateBoundaryHandling(pdfField, indexArr, lbMethod, boundaryFunctor):
+ dim = lbMethod.dim
cellLoopBody = Block([])
cellLoop = LoopOverCoordinate(cellLoopBody, coordinateToLoopOver=0, start=0, stop=indexArr.shape[0])
@@ -145,7 +147,7 @@ def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
dirSymbol = TypedSymbol("dir", "int")
cellLoopBody.append(SympyAssignment(dirSymbol, indexField[0](dim)))
- boundaryEqList = boundaryFunctor(pdfField, dirSymbol, latticeModel)
+ boundaryEqList = boundaryFunctor(pdfField, dirSymbol, lbMethod)
if type(boundaryEqList) is tuple:
boundaryEqList, additionalNodes = boundaryEqList
else:
@@ -168,7 +170,7 @@ def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
for node in additionalNodes:
loop.body.append(node)
- functionBody.insertFront(LbmMethodInfo(latticeModel))
+ functionBody.insertFront(LbmMethodInfo(lbMethod))
fixedCoordinateMapping = {f.name: coordinateSymbols[:dim] for f in fieldsAccessed}
resolveFieldAccesses(ast, set(['indexField']), fieldToFixedCoordinates=fixedCoordinateMapping)
diff --git a/boundaries/createindexlist.py b/boundaries/createindexlist.py
index 52b52662..83ce8575 100644
--- a/boundaries/createindexlist.py
+++ b/boundaries/createindexlist.py
@@ -30,10 +30,11 @@ def createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGho
if cythonFuncsAvailable:
stencil = np.array(stencil, dtype=np.int32)
if len(flagField.shape) == 2:
- return np.array(createBoundaryIndexList2D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil))
+ idxList = createBoundaryIndexList2D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
elif len(flagField.shape) == 3:
- return np.array(createBoundaryIndexList3D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil))
+ idxList = createBoundaryIndexList3D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
else:
raise ValueError("Flag field has to be a 2 or 3 dimensional numpy array")
+ return np.array(idxList, dtype=np.int32)
else:
return _createBoundaryIndexListPython(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
diff --git a/boundaries/createindexlistcython.pyx b/boundaries/createindexlistcython.pyx
index 02034844..94acf552 100644
--- a/boundaries/createindexlistcython.pyx
+++ b/boundaries/createindexlistcython.pyx
@@ -4,6 +4,7 @@ ctypedef fused IntegerType:
short
int
long
+ long long
unsigned short
unsigned int
unsigned long
diff --git a/creationfunctions.py b/creationfunctions.py
index 2ed6ec43..1409f05b 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -47,10 +47,12 @@ def _getParams(params, optParams):
return paramsResult, optParamsResult
-def createLatticeBoltzmannFunction(optimizationParams={}, **kwargs):
+def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
params, optParams = _getParams(kwargs, optimizationParams)
- ast = createLatticeBoltzmannKernel(**params, optimizationParams=optParams)
+ if ast is None:
+ ast = createLatticeBoltzmannAst(**params, optimizationParams=optParams)
+
if params['target'] == 'cpu':
from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
if 'openMP' in optParams:
@@ -58,18 +60,22 @@ def createLatticeBoltzmannFunction(optimizationParams={}, **kwargs):
addOpenMP(ast)
elif isinstance(optParams['openMP'], int):
addOpenMP(ast, numThreads=optParams['openMP'])
- return makePythonCpuFunction(ast)
+ res = makePythonCpuFunction(ast)
elif params['target'] == 'gpu':
from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
- return makePythonGpuFunction(ast)
+ res = makePythonGpuFunction(ast)
else:
return ValueError("'target' has to be either 'cpu' or 'gpu'")
+ res.method = ast.method
+ return res
+
-def createLatticeBoltzmannKernel(optimizationParams={}, **kwargs):
+def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
params, optParams = _getParams(kwargs, optimizationParams)
- updateRule = createLatticeBoltzmannUpdateRule(**params, optimizationParams=optimizationParams)
+ if updateRule is None:
+ updateRule = createLatticeBoltzmannUpdateRule(**params, optimizationParams=optimizationParams)
if params['target'] == 'cpu':
from pystencils.cpu import createKernel
@@ -77,26 +83,31 @@ def createLatticeBoltzmannKernel(optimizationParams={}, **kwargs):
splitGroups = updateRule.simplificationHints['splitGroups']
else:
splitGroups = ()
- return createKernel(updateRule.allEquations, splitGroups=splitGroups)
+ res = createKernel(updateRule.allEquations, splitGroups=splitGroups)
elif params['target'] == 'gpu':
from pystencils.gpucuda import createCUDAKernel
- return createCUDAKernel(updateRule.allEquations)
+ res = createCUDAKernel(updateRule.allEquations)
else:
return ValueError("'target' has to be either 'cpu' or 'gpu'")
+ res.method = updateRule.method
+ return res
-def createLatticeBoltzmannUpdateRule(optimizationParams={}, **kwargs):
+
+def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwargs):
params, optParams = _getParams(kwargs, optimizationParams)
stencil = getStencil(params['stencil'])
- method = createLatticeBoltzmannCollisionRule(**params)
+
+ if lbMethod is None:
+ lbMethod = createLatticeBoltzmannMethod(**params)
splitInnerLoop = 'split' in optParams and optParams['split']
dirCSE = 'doCseInOpposingDirections'
doCseInOpposingDirections = False if dirCSE not in optParams else optParams[dirCSE]
doOverallCse = False if 'doOverallCse' not in optParams else optParams['doOverallCse']
- simplification = createSimplificationStrategy(method, doCseInOpposingDirections, doOverallCse, splitInnerLoop)
- collisionRule = simplification(method.getCollisionRule())
+ simplification = createSimplificationStrategy(lbMethod, doCseInOpposingDirections, doOverallCse, splitInnerLoop)
+ collisionRule = simplification(lbMethod.getCollisionRule())
if 'fieldSize' in optParams and optParams['fieldSize']:
npField = createPdfArray(optParams['fieldSize'], len(stencil), layout=optParams['fieldLayout'])
@@ -111,7 +122,7 @@ def createLatticeBoltzmannUpdateRule(optimizationParams={}, **kwargs):
return updateRule
-def createLatticeBoltzmannCollisionRule(**params):
+def createLatticeBoltzmannMethod(**params):
params, _ = _getParams(params, {})
stencil = getStencil(params['stencil'])
diff --git a/macroscopicValueKernels.py b/macroscopicValueKernels.py
index 2e911a4e..ea03843b 100644
--- a/macroscopicValueKernels.py
+++ b/macroscopicValueKernels.py
@@ -1,7 +1,7 @@
import sympy as sp
-
-from lbmpy.simplificationfactory import createSimplificationStrategy
+from copy import deepcopy
from pystencils.field import Field
+from lbmpy.simplificationfactory import createSimplificationStrategy
def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fieldLayout='numpy', target='cpu'):
@@ -75,6 +75,9 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field
:param lbMethod: instance of :class:`lbmpy.methods.AbstractLbMethod`
:param quantitiesToSet: map from conserved quantity name to fixed value or numpy array
: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 target: 'cpu' or 'gpu'
+ :return: function taking pdf array as single argument and which sets the field to the given values
"""
if pdfArr is not None:
pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
@@ -82,22 +85,28 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field
pdfField = Field.createGeneric('pdfs', lbMethod.dim, indexDimensions=1, layout=fieldLayout)
fixedKernelParameters = {}
+ cqc = lbMethod.conservedQuantityComputation
valueMap = {}
atLeastOneFieldInput = False
for quantityName, value in quantitiesToSet.items():
if hasattr(value, 'shape'):
fixedKernelParameters[quantityName] = value
- value = Field.createFromNumpyArray(quantityName, value)
atLeastOneFieldInput = True
+ numComponents = cqc.conservedQuantities[quantityName]
+ field = Field.createFromNumpyArray(quantityName, value, indexDimensions=0 if numComponents <= 1 else 1)
+ if numComponents == 1:
+ value = field(0)
+ else:
+ value = [field(i) for i in range(numComponents)]
+
valueMap[quantityName] = value
- cqc = lbMethod.conservedQuantityComputation
cqEquations = cqc.equilibriumInputEquationsFromInitValues(**valueMap)
eq = lbMethod.getEquilibrium(conservedQuantityEquations=cqEquations)
if atLeastOneFieldInput:
- simplification = createSimplificationStrategy(eq)
+ simplification = createSimplificationStrategy(lbMethod)
eq = simplification(eq)
else:
eq = eq.insertSubexpressions()
@@ -123,32 +132,37 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field
return setter
-def compileAdvancedVelocitySetter(collisionRule, velocityArray, pdfArr=None):
+def compileAdvancedVelocitySetter(lbMethod, velocityArray, velocityRelaxationRate=1.3, pdfArr=None):
"""
Advanced initialization of velocity field through iteration procedure according to
Mei, Luo, Lallemand and Humieres: Consistent initial conditions for LBM simulations, 2005
- Important: this procedure only works if a non-zero relaxation rate was used for the velocity moments!
-
- :param collisionRule: unsimplified collision rule
+ :param lbMethod:
:param velocityArray: array with velocity field
- :param pdfArr: optional array, to compile kernel with fixed layout and shape
- :return: function, that has to be called multiple times, with a pdf field (src/dst) until convergence
- similar to the normal streamCollide step, also with boundary handling
+ :param velocityRelaxationRate: relaxation rate for the velocity moments - determines convergence behaviour
+ of the initialization scheme
+ :return: collision rule
"""
- velocityField = Field.createFromNumpyArray('vel', velocityArray, indexDimensions=1)
+ velocityField = Field.createFromNumpyArray('velInput', velocityArray, indexDimensions=1)
+
+ cqc = lbMethod.conservedQuantityComputation
+ densitySymbol = cqc.definedSymbols(order=0)[1]
+ velocitySymbols = cqc.definedSymbols(order=1)[1]
+
+ # density is computed from pdfs
+ eqInputFromPdfs = cqc.equilibriumInputEquationsFromPdfs(lbMethod.preCollisionPdfSymbols)
+ eqInputFromPdfs = eqInputFromPdfs.extract([densitySymbol])
+ # velocity is read from input field
+ velSymbols = [velocityField(i) for i in range(lbMethod.dim)]
+ eqInputFromField = cqc.equilibriumInputEquationsFromInitValues(velocity=velSymbols)
+ eqInputFromField = eqInputFromField.extract(velocitySymbols)
+ # then both are merged together
+ eqInput = eqInputFromPdfs.merge(eqInputFromField)
- # create normal LBM kernel and replace velocity by expressions of velocity field
- from lbmpy_old.simplifications import sympyCSE
- latticeModel = collisionRule.latticeModel
- collisionRule = sympyCSE(collisionRule)
- collisionRule = streamPullWithSourceAndDestinationFields(collisionRule, pdfArr)
+ # set first order relaxation rate
+ lbMethod = deepcopy(lbMethod)
+ lbMethod.setFirstMomentRelaxationRate(velocityRelaxationRate)
- replacements = {u_i: sp.Eq(u_i, velocityField(i)) for i, u_i in enumerate(latticeModel.symbolicVelocity)}
+ return lbMethod.getCollisionRule(eqInput)
- newSubExpressions = [replacements[eq.lhs] if eq.lhs in replacements else eq for eq in collisionRule.subexpressions]
- newCollisionRule = LbmCollisionRule(collisionRule.updateEquations, newSubExpressions,
- latticeModel, collisionRule.updateEquationDirections)
- kernelAst = createKernel(newCollisionRule.equations)
- return makePythonFunction(kernelAst, {'vel': velocityArray})
diff --git a/methods/momentbased.py b/methods/momentbased.py
index 05075f76..0495980d 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -35,23 +35,14 @@ class MomentBasedLbMethod(AbstractLbMethod):
assert isinstance(conservedQuantityComputation, AbstractConservedQuantityComputation)
- moments = []
- relaxationRates = []
- equilibriumMoments = []
- for moment, relaxInfo in momentToRelaxationInfoDict.items():
- moments.append(moment)
- relaxationRates.append(relaxInfo.relaxationRate)
- equilibriumMoments.append(relaxInfo.equilibriumValue)
-
self._forceModel = forceModel
- self._moments = moments
- self._momentToRelaxationInfoDict = momentToRelaxationInfoDict
- self._momentMatrix = momentMatrix(moments, self.stencil)
- self._relaxationRates = sp.Matrix(relaxationRates)
- self._equilibriumMoments = sp.Matrix(equilibriumMoments)
+ self._momentToRelaxationInfoDict = OrderedDict(momentToRelaxationInfoDict.items())
self._conservedQuantityComputation = conservedQuantityComputation
- symbolsInEquilibriumMoments = self._equilibriumMoments.atoms(sp.Symbol)
+ equilibriumMoments = []
+ for moment, relaxInfo in momentToRelaxationInfoDict.items():
+ equilibriumMoments.append(relaxInfo.equilibriumValue)
+ symbolsInEquilibriumMoments = sp.Matrix(equilibriumMoments).atoms(sp.Symbol)
conservedQuantities = set()
for v in self._conservedQuantityComputation.definedSymbols().values():
if isinstance(v, collections.Sequence):
@@ -65,6 +56,14 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._weights = None
+ 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"
+ for e in MOMENT_SYMBOLS[:self.dim]:
+ prevEntry = self._momentToRelaxationInfoDict[e]
+ newEntry = RelaxationInfo(prevEntry[0], relaxationRate)
+ self._momentToRelaxationInfoDict[e] = newEntry
+
def _repr_html_(self):
table = """
<table style="border:none; width: 100%">
@@ -77,7 +76,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
</table>
"""
content = ""
- for rr, moment, eqValue in zip(self._relaxationRates, self._moments, self._equilibriumMoments):
+ for moment, (eqValue, rr) in self._momentToRelaxationInfoDict.items():
vals = {
'rr': sp.latex(rr),
'moment': sp.latex(moment),
@@ -93,7 +92,15 @@ class MomentBasedLbMethod(AbstractLbMethod):
@property
def moments(self):
- return self._moments
+ return tuple(self._momentToRelaxationInfoDict.keys())
+
+ @property
+ def momentEquilibriumValues(self):
+ return tuple([e.equilibriumValue for e in self._momentToRelaxationInfoDict.values()])
+
+ @property
+ def relaxationRates(self):
+ return tuple([e.relaxationRate for e in self._momentToRelaxationInfoDict.values()])
@property
def zerothOrderEquilibriumMomentSymbol(self, ):
@@ -144,13 +151,13 @@ class MomentBasedLbMethod(AbstractLbMethod):
"Can not determine their relaxation rate automatically")
def getEquilibrium(self, conservedQuantityEquations=None):
- D = sp.eye(len(self._relaxationRates))
+ D = sp.eye(len(self.relaxationRates))
return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations=conservedQuantityEquations)
- def getCollisionRule(self):
- D = sp.diag(*self._relaxationRates)
+ def getCollisionRule(self, conservedQuantityEquations=None):
+ D = sp.diag(*self.relaxationRates)
relaxationRateSubExpressions, D = self._generateRelaxationMatrix(D)
- eqColl = self._getCollisionRuleWithRelaxationMatrix(D, relaxationRateSubExpressions)
+ eqColl = self._getCollisionRuleWithRelaxationMatrix(D, relaxationRateSubExpressions, conservedQuantityEquations)
if self._forceModel is not None:
forceModelTerms = self._forceModel(self)
newEqs = [sp.Eq(eq.lhs, eq.rhs + fmt) for eq, fmt in zip(eqColl.mainEquations, forceModelTerms)]
@@ -163,8 +170,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=[], conservedQuantityEquations=None):
f = sp.Matrix(self.preCollisionPdfSymbols)
- M = self._momentMatrix
- m_eq = self._equilibriumMoments
+ M = momentMatrix(self.moments, self.stencil)
+ m_eq = sp.Matrix(self.momentEquilibriumValues)
collisionRule = f + M.inv() * D * (m_eq - M * f)
collisionEqs = [sp.Eq(lhs, rhs) for lhs, rhs in zip(self.postCollisionPdfSymbols, collisionRule)]
@@ -316,7 +323,7 @@ def createSRT(stencil, relaxationRate, compressible=False, forceModel=None, equi
:return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
"""
moments = getDefaultMomentSetForStencil(stencil)
- rrDict = {m: relaxationRate for m in moments}
+ rrDict = OrderedDict([(m, relaxationRate) for m in moments])
return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, compressible, forceModel, equilibriumAccuracyOrder)
@@ -332,7 +339,7 @@ def createTRT(stencil, relaxationRateEvenMoments, relaxationRateOddMoments, comp
If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.createTRTWithMagicNumber`.
"""
moments = getDefaultMomentSetForStencil(stencil)
- rrDict = {m: relaxationRateEvenMoments if isEven(m) else relaxationRateOddMoments for m in moments}
+ rrDict = OrderedDict([(m, relaxationRateEvenMoments if isEven(m) else relaxationRateOddMoments) for m in moments])
return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, compressible, forceModel, equilibriumAccuracyOrder)
diff --git a/updatekernels.py b/updatekernels.py
index b3607e00..cab9c760 100644
--- a/updatekernels.py
+++ b/updatekernels.py
@@ -57,7 +57,8 @@ def createStreamPullKernel(collisionRule, numpyField=None, srcFieldName="src", d
"""
dim = collisionRule.method.dim
if numpyField is not None:
- assert len(numpyField.shape) == dim + 1
+ assert len(numpyField.shape) == dim + 1, "Field dimension mismatch: dimension is %s, should be %d" % \
+ (len(numpyField.shape), dim + 1)
if numpyField is None:
src = Field.createGeneric(srcFieldName, dim, indexDimensions=1, layout=genericLayout, dtype=genericFieldType)
--
GitLab