diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index b80a281b17f44375197e84c060931d572dae7675..a9688fd097ece5e55f4dc9fb93937400d529eb8d 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -9,6 +9,8 @@ def noSlip(pdfField, direction, lbmMethod):
def ubb(pdfField, direction, lbmMethod, velocity):
+ 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)
inverseDir = invDir(direction)
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index dba1bd2db4eed5af19674fd04ac7b7bb93d1818e..a71a28e3a65d5664883482ecd1d816a6d5c7bdf3 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -6,7 +6,7 @@ from pystencils.ast import Block, SympyAssignment, LoopOverCoordinate, KernelFun
from pystencils.transformations import moveConstantsBeforeLoop, resolveFieldAccesses, typingFromSympyInspection, \
typeAllEquations
from pystencils.cpu import makePythonFunction
-
+from lbmpy.boundaries.createindexlist import createBoundaryIndexList
INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
@@ -41,16 +41,16 @@ class BoundaryHandling:
def getFlag(self, name):
return 2 ** self._nameToIndex[name]
- def setBoundary(self, name, indexExpr, clearOtherBoundaries=True):
- if not isinstance(name, str):
- function = name
- if hasattr(function, '__name__'):
- name = function.__name__
- else:
- name = function.name
+ def setBoundary(self, function, indexExpr, clearOtherBoundaries=True):
+ if hasattr(function, '__name__'):
+ name = function.__name__
+ elif hasattr(function, 'name'):
+ name = function.name
+ else:
+ raise ValueError("Boundary function has to have a '__name__' or 'name' attribute")
- if function not in self._boundaryFunctions:
- self.addBoundary(function, name)
+ if function not in self._boundaryFunctions:
+ self.addBoundary(function, name)
flag = self.getFlag(name)
if clearOtherBoundaries:
@@ -66,8 +66,8 @@ class BoundaryHandling:
def prepare(self):
self.invalidateIndexCache()
for boundaryIdx, boundaryFunc in enumerate(self._boundaryFunctions):
- idxField = createBoundaryIndexList(self.flagField, self._ghostLayers, self._latticeModel.stencil,
- 2 ** boundaryIdx, self._fluidFlag)
+ idxField = createBoundaryIndexList(self.flagField, self._latticeModel.stencil,
+ 2 ** boundaryIdx, self._fluidFlag, self._ghostLayers)
ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._latticeModel, boundaryFunc)
self._boundarySweeps.append(makePythonFunction(ast, {'indexField': idxField}))
@@ -99,14 +99,14 @@ def weightOfDirection(dirIdx):
# ------------------------------------- Kernel Generation --------------------------------------------------------------
-class LatticeModelInfo(CustomCppCode):
- def __init__(self, latticeModel):
- stencil = latticeModel.stencil
- symbolsDefined = set(offsetSymbols(latticeModel.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL])
+class LbmMethodInfo(CustomCppCode):
+ def __init__(self, lbMethod):
+ stencil = lbMethod.stencil
+ symbolsDefined = set(offsetSymbols(lbMethod.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL])
- offsetSym = offsetSymbols(latticeModel.dim)
+ offsetSym = offsetSymbols(lbMethod.dim)
code = "\n"
- for i in range(latticeModel.dim):
+ for i in range(lbMethod.dim):
offsetStr = ", ".join([str(d[i]) for d in stencil])
code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr)
@@ -116,9 +116,9 @@ class LatticeModelInfo(CustomCppCode):
invDirs.append(str(stencil.index(inverseDir)))
code += "static const int %s [] = { %s };\n" % (INV_DIR_SYMBOL.name, ", ".join(invDirs))
- weights = [str(w.evalf()) for w in latticeModel.weights]
+ weights = [str(w.evalf()) for w in lbMethod.weights]
code += "static const double %s [] = { %s };\n" % (WEIGHTS_SYMBOL.name, ",".join(weights))
- super(LatticeModelInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
+ super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
@@ -158,7 +158,7 @@ def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
for node in additionalNodes:
loop.body.append(node)
- functionBody.insertFront(LatticeModelInfo(latticeModel))
+ functionBody.insertFront(LbmMethodInfo(latticeModel))
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 8f3fb6e772ba10276054830590094c9f7f854e51..52b526625bc2721b8d293ba2ec9e79252a801851 100644
--- a/boundaries/createindexlist.py
+++ b/boundaries/createindexlist.py
@@ -1,16 +1,15 @@
import numpy as np
import itertools
-#try:
-if True:
+try:
import pyximport;
pyximport.install()
from lbmpy.boundaries.createindexlistcython import createBoundaryIndexList2D, createBoundaryIndexList3D
cythonFuncsAvailable = True
-#except Exception:
-# cythonFuncsAvailable = False
-# createBoundaryIndexList2D = None
-# createBoundaryIndexList3D = None
+except Exception:
+ cythonFuncsAvailable = False
+ createBoundaryIndexList2D = None
+ createBoundaryIndexList3D = None
def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, fluidMask, stencil):
diff --git a/methods/momentbased.py b/methods/momentbased.py
index d4f166c42f0d6ae630c3ff25f0f40fcd54bb6f16..282dd8820a5cb65363a750f105c0b950ce07e198 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -144,12 +144,14 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
"Can not determine their relaxation rate automatically")
def getEquilibrium(self):
+ sp.factor
D = sp.eye(len(self._relaxationRates))
return self._getCollisionRuleWithRelaxationMatrix(D)
def getCollisionRule(self):
D = sp.diag(*self._relaxationRates)
- eqColl = self._getCollisionRuleWithRelaxationMatrix(D)
+ relaxationRateSubExpressions, D = self._generateRelaxationMatrix(D)
+ eqColl = self._getCollisionRuleWithRelaxationMatrix(D, relaxationRateSubExpressions)
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)]
@@ -160,13 +162,11 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
def conservedQuantityComputation(self):
return self._conservedQuantityComputation
- def _getCollisionRuleWithRelaxationMatrix(self, D):
+ def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=[]):
f = sp.Matrix(self.preCollisionPdfSymbols)
M = self._momentMatrix
m_eq = self._equilibriumMoments
- relaxationRateSubExpressions, D = self._generateRelaxationMatrix(D)
-
collisionRule = f + M.inv() * D * (m_eq - M * f)
collisionEqs = [sp.Eq(lhs, rhs) for lhs, rhs in zip(self.postCollisionPdfSymbols, collisionRule)]
@@ -175,7 +175,7 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
simplificationHints.update(self._conservedQuantityComputation.definedSymbols())
simplificationHints['relaxationRates'] = D.atoms(sp.Symbol)
- allSubexpressions = relaxationRateSubExpressions + eqValueEqs.subexpressions + eqValueEqs.mainEquations
+ allSubexpressions = additionalSubexpressions + eqValueEqs.subexpressions + eqValueEqs.mainEquations
return LbmCollisionRule(self, collisionEqs, allSubexpressions,
simplificationHints)
diff --git a/simplificationfactory.py b/simplificationfactory.py
index 2a3c01d0f0e9d97ef969784ea2a2c1f6ceacafa3..152d927e9b4a8ea15817b7d090a072c65fc3c66e 100644
--- a/simplificationfactory.py
+++ b/simplificationfactory.py
@@ -36,6 +36,7 @@ def createSimplificationStrategy(lbmMethod, doCseInOpposingDirections=False, doO
return s
+
if __name__ == '__main__':
from lbmpy.stencils import getStencil
from lbmpy.methods.momentbased import createOrthogonalMRT