From 54c2055c0945f33db76ec5f97f29b32bc5b08fb1 Mon Sep 17 00:00:00 2001
From: Martin Bauer <bauer_martin@gmx.de>
Date: Mon, 23 Jan 2017 13:57:17 +0100
Subject: [PATCH] new lbmpy: boundaries
---
boundaries/__init__.py | 2 ++
boundaries/boundaryconditions.py | 44 ++++++++++++++++----------------
boundaries/boundaryhandling.py | 16 +++++++++---
methods/momentbased.py | 1 -
4 files changed, 37 insertions(+), 26 deletions(-)
diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index e69de29b..882a2a77 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -0,0 +1,2 @@
+from lbmpy.boundaries.boundaryconditions import noSlip, ubb, fixedDensity
+from lbmpy.boundaries.boundaryhandling import BoundaryHandling
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index a9688fd0..293f3ad9 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -1,16 +1,22 @@
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))]
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index a71a28e3..407c6863 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -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:
diff --git a/methods/momentbased.py b/methods/momentbased.py
index 282dd882..baf71515 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -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)
--
GitLab