Commit 1913fa42 authored by Martin Bauer's avatar Martin Bauer
Browse files

boundary generatlization

parent fde2c8ce
...@@ -9,6 +9,8 @@ def noSlip(pdfField, direction, lbmMethod): ...@@ -9,6 +9,8 @@ def noSlip(pdfField, direction, lbmMethod):
def ubb(pdfField, direction, lbmMethod, velocity): 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) neighbor = offsetFromDir(direction, lbmMethod.dim)
inverseDir = invDir(direction) inverseDir = invDir(direction)
......
...@@ -6,7 +6,7 @@ from pystencils.ast import Block, SympyAssignment, LoopOverCoordinate, KernelFun ...@@ -6,7 +6,7 @@ from pystencils.ast import Block, SympyAssignment, LoopOverCoordinate, KernelFun
from pystencils.transformations import moveConstantsBeforeLoop, resolveFieldAccesses, typingFromSympyInspection, \ from pystencils.transformations import moveConstantsBeforeLoop, resolveFieldAccesses, typingFromSympyInspection, \
typeAllEquations typeAllEquations
from pystencils.cpu import makePythonFunction from pystencils.cpu import makePythonFunction
from lbmpy.boundaries.createindexlist import createBoundaryIndexList
INV_DIR_SYMBOL = TypedSymbol("invDir", "int") INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
WEIGHTS_SYMBOL = TypedSymbol("weights", "double") WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
...@@ -41,16 +41,16 @@ class BoundaryHandling: ...@@ -41,16 +41,16 @@ class BoundaryHandling:
def getFlag(self, name): def getFlag(self, name):
return 2 ** self._nameToIndex[name] return 2 ** self._nameToIndex[name]
def setBoundary(self, name, indexExpr, clearOtherBoundaries=True): def setBoundary(self, function, indexExpr, clearOtherBoundaries=True):
if not isinstance(name, str): if hasattr(function, '__name__'):
function = name name = function.__name__
if hasattr(function, '__name__'): elif hasattr(function, 'name'):
name = function.__name__ name = function.name
else: else:
name = function.name raise ValueError("Boundary function has to have a '__name__' or 'name' attribute")
if function not in self._boundaryFunctions: if function not in self._boundaryFunctions:
self.addBoundary(function, name) self.addBoundary(function, name)
flag = self.getFlag(name) flag = self.getFlag(name)
if clearOtherBoundaries: if clearOtherBoundaries:
...@@ -66,8 +66,8 @@ class BoundaryHandling: ...@@ -66,8 +66,8 @@ class BoundaryHandling:
def prepare(self): def prepare(self):
self.invalidateIndexCache() self.invalidateIndexCache()
for boundaryIdx, boundaryFunc in enumerate(self._boundaryFunctions): for boundaryIdx, boundaryFunc in enumerate(self._boundaryFunctions):
idxField = createBoundaryIndexList(self.flagField, self._ghostLayers, self._latticeModel.stencil, idxField = createBoundaryIndexList(self.flagField, self._latticeModel.stencil,
2 ** boundaryIdx, self._fluidFlag) 2 ** boundaryIdx, self._fluidFlag, self._ghostLayers)
ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._latticeModel, boundaryFunc) ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._latticeModel, boundaryFunc)
self._boundarySweeps.append(makePythonFunction(ast, {'indexField': idxField})) self._boundarySweeps.append(makePythonFunction(ast, {'indexField': idxField}))
...@@ -99,14 +99,14 @@ def weightOfDirection(dirIdx): ...@@ -99,14 +99,14 @@ def weightOfDirection(dirIdx):
# ------------------------------------- Kernel Generation -------------------------------------------------------------- # ------------------------------------- Kernel Generation --------------------------------------------------------------
class LatticeModelInfo(CustomCppCode): class LbmMethodInfo(CustomCppCode):
def __init__(self, latticeModel): def __init__(self, lbMethod):
stencil = latticeModel.stencil stencil = lbMethod.stencil
symbolsDefined = set(offsetSymbols(latticeModel.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL]) symbolsDefined = set(offsetSymbols(lbMethod.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL])
offsetSym = offsetSymbols(latticeModel.dim) offsetSym = offsetSymbols(lbMethod.dim)
code = "\n" code = "\n"
for i in range(latticeModel.dim): for i in range(lbMethod.dim):
offsetStr = ", ".join([str(d[i]) for d in stencil]) offsetStr = ", ".join([str(d[i]) for d in stencil])
code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr) code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr)
...@@ -116,9 +116,9 @@ class LatticeModelInfo(CustomCppCode): ...@@ -116,9 +116,9 @@ class LatticeModelInfo(CustomCppCode):
invDirs.append(str(stencil.index(inverseDir))) invDirs.append(str(stencil.index(inverseDir)))
code += "static const int %s [] = { %s };\n" % (INV_DIR_SYMBOL.name, ", ".join(invDirs)) 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)) 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): def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
...@@ -158,7 +158,7 @@ def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor): ...@@ -158,7 +158,7 @@ def generateBoundaryHandling(pdfField, indexArr, latticeModel, boundaryFunctor):
for node in additionalNodes: for node in additionalNodes:
loop.body.append(node) loop.body.append(node)
functionBody.insertFront(LatticeModelInfo(latticeModel)) functionBody.insertFront(LbmMethodInfo(latticeModel))
fixedCoordinateMapping = {f.name: coordinateSymbols[:dim] for f in fieldsAccessed} fixedCoordinateMapping = {f.name: coordinateSymbols[:dim] for f in fieldsAccessed}
resolveFieldAccesses(ast, set(['indexField']), fieldToFixedCoordinates=fixedCoordinateMapping) resolveFieldAccesses(ast, set(['indexField']), fieldToFixedCoordinates=fixedCoordinateMapping)
......
import numpy as np import numpy as np
import itertools import itertools
#try: try:
if True:
import pyximport; import pyximport;
pyximport.install() pyximport.install()
from lbmpy.boundaries.createindexlistcython import createBoundaryIndexList2D, createBoundaryIndexList3D from lbmpy.boundaries.createindexlistcython import createBoundaryIndexList2D, createBoundaryIndexList3D
cythonFuncsAvailable = True cythonFuncsAvailable = True
#except Exception: except Exception:
# cythonFuncsAvailable = False cythonFuncsAvailable = False
# createBoundaryIndexList2D = None createBoundaryIndexList2D = None
# createBoundaryIndexList3D = None createBoundaryIndexList3D = None
def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, fluidMask, stencil): def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, fluidMask, stencil):
......
...@@ -144,12 +144,14 @@ class MomentBasedLbmMethod(AbstractLbmMethod): ...@@ -144,12 +144,14 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
"Can not determine their relaxation rate automatically") "Can not determine their relaxation rate automatically")
def getEquilibrium(self): def getEquilibrium(self):
sp.factor
D = sp.eye(len(self._relaxationRates)) D = sp.eye(len(self._relaxationRates))
return self._getCollisionRuleWithRelaxationMatrix(D) return self._getCollisionRuleWithRelaxationMatrix(D)
def getCollisionRule(self): def getCollisionRule(self):
D = sp.diag(*self._relaxationRates) 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: if self._forceModel is not None:
forceModelTerms = self._forceModel(self) forceModelTerms = self._forceModel(self)
newEqs = [sp.Eq(eq.lhs, eq.rhs + fmt) for eq, fmt in zip(eqColl.mainEquations, forceModelTerms)] newEqs = [sp.Eq(eq.lhs, eq.rhs + fmt) for eq, fmt in zip(eqColl.mainEquations, forceModelTerms)]
...@@ -160,13 +162,11 @@ class MomentBasedLbmMethod(AbstractLbmMethod): ...@@ -160,13 +162,11 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
def conservedQuantityComputation(self): def conservedQuantityComputation(self):
return self._conservedQuantityComputation return self._conservedQuantityComputation
def _getCollisionRuleWithRelaxationMatrix(self, D): def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=[]):
f = sp.Matrix(self.preCollisionPdfSymbols) f = sp.Matrix(self.preCollisionPdfSymbols)
M = self._momentMatrix M = self._momentMatrix
m_eq = self._equilibriumMoments m_eq = self._equilibriumMoments
relaxationRateSubExpressions, D = self._generateRelaxationMatrix(D)
collisionRule = f + M.inv() * D * (m_eq - M * f) collisionRule = f + M.inv() * D * (m_eq - M * f)
collisionEqs = [sp.Eq(lhs, rhs) for lhs, rhs in zip(self.postCollisionPdfSymbols, collisionRule)] collisionEqs = [sp.Eq(lhs, rhs) for lhs, rhs in zip(self.postCollisionPdfSymbols, collisionRule)]
...@@ -175,7 +175,7 @@ class MomentBasedLbmMethod(AbstractLbmMethod): ...@@ -175,7 +175,7 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
simplificationHints.update(self._conservedQuantityComputation.definedSymbols()) simplificationHints.update(self._conservedQuantityComputation.definedSymbols())
simplificationHints['relaxationRates'] = D.atoms(sp.Symbol) simplificationHints['relaxationRates'] = D.atoms(sp.Symbol)
allSubexpressions = relaxationRateSubExpressions + eqValueEqs.subexpressions + eqValueEqs.mainEquations allSubexpressions = additionalSubexpressions + eqValueEqs.subexpressions + eqValueEqs.mainEquations
return LbmCollisionRule(self, collisionEqs, allSubexpressions, return LbmCollisionRule(self, collisionEqs, allSubexpressions,
simplificationHints) simplificationHints)
......
...@@ -36,6 +36,7 @@ def createSimplificationStrategy(lbmMethod, doCseInOpposingDirections=False, doO ...@@ -36,6 +36,7 @@ def createSimplificationStrategy(lbmMethod, doCseInOpposingDirections=False, doO
return s return s
if __name__ == '__main__': if __name__ == '__main__':
from lbmpy.stencils import getStencil from lbmpy.stencils import getStencil
from lbmpy.methods.momentbased import createOrthogonalMRT from lbmpy.methods.momentbased import createOrthogonalMRT
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment