Skip to content
GitLab
Commit fd57b765 authored by Martin Bauer's avatar Martin Bauer
Browse files

new lbm module: macroscopic values & bugfixes

parent f5ad4934
Hide whitespace changes
Inline Side-by-side
boundaries/boundaryconditions.py
View file @ fd57b765
......@@ -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)]
......
boundaries/boundaryhandling.py
View file @ fd57b765
......@@ -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)
......
boundaries/createindexlist.py
View file @ fd57b765
......@@ -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)
boundaries/createindexlistcython.pyx
View file @ fd57b765
......@@ -4,6 +4,7 @@ ctypedef fused IntegerType:
short
int
long
long long
unsigned short
unsigned int
unsigned long
......
creationfunctions.py
View file @ fd57b765
......@@ -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'])
......
macroscopicValueKernels.py
View file @ fd57b765
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})
methods/momentbased.py
View file @ fd57b765
......@@ -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)
......
updatekernels.py
View file @ fd57b765
......@@ -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)
......
Supports Markdown
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