Commit 54c2055c authored by Martin Bauer's avatar Martin Bauer Committed by Martin Bauer
Browse files

new lbmpy: boundaries

parent 1913fa42
from lbmpy.boundaries.boundaryconditions import noSlip, ubb, fixedDensity
from lbmpy.boundaries.boundaryhandling import BoundaryHandling
import sympy as sp
from lbmpy.simplificationfactory import createSimplificationStrategy
from pystencils.sympyextensions import getSymmetricPart
from lbmpy.boundaries.boundaryhandling import offsetFromDir, weightOfDirection, invDir
def noSlip(pdfField, direction, lbmMethod):
"""No-Slip, simple bounce back boundary condition, enforcing zero velocity at obstacle"""
neighbor = offsetFromDir(direction, lbmMethod.dim)
inverseDir = invDir(direction)
return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(direction))]
def ubb(pdfField, direction, lbmMethod, 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))
"Dimension of velocity (%d) does not match dimension of LB method (%d)" % (len(velocity, lbmMethod.dim))
neighbor = offsetFromDir(direction, lbmMethod.dim)
inverseDir = invDir(direction)
......@@ -19,32 +25,26 @@ def ubb(pdfField, direction, lbmMethod, velocity):
pdfField(direction) - velTerm)]
def fixedDensity(pdfField, direction, latticeModel, density):
from lbmpy_old.equilibria import standardDiscreteEquilibrium
neighbor = offsetFromDir(direction, latticeModel.dim)
inverseDir = invDir(direction)
stencil = latticeModel.stencil
def fixedDensity(pdfField, direction, lbmMethod, density):
"""Boundary condition that fixes the density/pressure at the obstacle"""
if not latticeModel.compressible:
density -= 1
def removeAsymmetricPartOfMainEquations(eqColl, dofs):
newMainEquations = [sp.Eq(e.lhs, getSymmetricPart(e.rhs, dofs)) for e in eqColl.mainEquations]
return eqColl.copy(newMainEquations)
eqParams = {'stencil': stencil,
'order': 2,
'c_s_sq': sp.Rational(1, 3),
'compressible': latticeModel.compressible,
'rho': density}
neighbor = offsetFromDir(direction, lbmMethod.dim)
inverseDir = invDir(direction)
u = sp.Matrix(latticeModel.symbolicVelocity)
symmetricEq = (standardDiscreteEquilibrium(u=u, **eqParams) + standardDiscreteEquilibrium(u=-u, **eqParams)) / 2
symmetricEq = removeAsymmetricPartOfMainEquations(lbmMethod.getEquilibrium())
simplification = createSimplificationStrategy(lbmMethod)
symmetricEq = simplification(symmetricEq)
subExpr1, rhoExpr, subExpr2, uExprs = getDensityVelocityExpressions(stencil,
[pdfField(i) for i in range(len(stencil))],
latticeModel.compressible)
subExprs = subExpr1 + [rhoExpr] + subExpr2 + uExprs
densitySymbol = lbmMethod.conservedQuantityComputation.definedSymbols()['density']
conditions = [(eq_i, sp.Equality(direction, i)) for i, eq_i in enumerate(symmetricEq)] + [(0, True)]
conditions = [(eq_i.rhs, sp.Equality(direction, i))
for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
eq_component = sp.Piecewise(*conditions)
return subExprs + [sp.Eq(pdfField[neighbor](inverseDir),
2 * eq_component - pdfField(direction))]
subExprs = [sp.Eq(eq.lhs, density if eq.lhs == densitySymbol else eq.rhs) for eq in symmetricEq.subexpressions]
return subExprs + [sp.Eq(pdfField[neighbor](inverseDir), 2 * eq_component - pdfField(direction))]
......@@ -5,7 +5,8 @@ from pystencils.backends.cbackend import CustomCppCode
from pystencils.ast import Block, SympyAssignment, LoopOverCoordinate, KernelFunction
from pystencils.transformations import moveConstantsBeforeLoop, resolveFieldAccesses, typingFromSympyInspection, \
typeAllEquations
from pystencils.cpu import makePythonFunction
from pystencils.cpu import makePythonFunction as makePythonCpuFunction
from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
from lbmpy.boundaries.createindexlist import createBoundaryIndexList
INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
......@@ -13,7 +14,7 @@ WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
class BoundaryHandling:
def __init__(self, symbolicPdfField, domainShape, latticeModel, ghostLayers=1):
def __init__(self, symbolicPdfField, domainShape, latticeModel, ghostLayers=1, target='cpu'):
self._symbolicPdfField = symbolicPdfField
self._shapeWithGhostLayers = [d + 2 * ghostLayers for d in domainShape]
self._fluidFlag = 2 ** 31
......@@ -23,6 +24,9 @@ class BoundaryHandling:
self._boundaryFunctions = []
self._nameToIndex = {}
self._boundarySweeps = []
self._target = target
if target not in ('cpu', 'gpu'):
raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
def addBoundary(self, boundaryFunction, name=None):
if name is None:
......@@ -69,7 +73,13 @@ class BoundaryHandling:
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}))
if self._target == 'cpu':
self._boundarySweeps.append(makePythonCpuFunction(ast, {'indexField': idxField}))
elif self._target == 'gpu':
self._boundarySweeps.append(makePythonGpuFunction(ast, {'indexField': idxField}))
else:
assert False
def __call__(self, **kwargs):
if len(self._boundarySweeps) == 0:
......
......@@ -144,7 +144,6 @@ 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)
......
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