diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index 882a2a77a45db4124562ff582c5b218c60881091..327e1bb2ce0cb7690f8e51b1a8ed811a1f3c3e39 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -1,2 +1,2 @@
-from lbmpy.boundaries.boundaryconditions import noSlip, ubb, fixedDensity
+from lbmpy.boundaries.boundaryconditions import NoSlip, NoSlipFullWay, UBB, FixedDensity
from lbmpy.boundaries.boundaryhandling import BoundaryHandling
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 3076b41671c28d67e8c6ea402c9ccab2a743e541..1fd31309a40af2581153e71cc1f15023d89ccfbc 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -4,81 +4,160 @@ from lbmpy.simplificationfactory import createSimplificationStrategy
from pystencils.sympyextensions import getSymmetricPart
from pystencils import Field
from lbmpy.boundaries.boundaryhandling import offsetFromDir, weightOfDirection, invDir
+from pystencils.types import createTypeFromString
-def noSlip(pdfField, direction, lbMethod):
- """No-Slip, simple bounce back boundary condition, enforcing zero velocity at obstacle"""
- neighbor = offsetFromDir(direction, lbMethod.dim)
- inverseDir = invDir(direction)
- return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(direction))]
+class Boundary(object):
+ """Base class all boundaries should derive from"""
+ def __call__(self, pdfField, directionSymbol, lbMethod, indexField):
+ """
+ This function defines the boundary behavior and must therefore be implemented by all boundaries.
+ Here the boundary is defined as a list of sympy equations, from which a boundary kernel is generated.
+ :param pdfField: pystencils field describing the pdf. The current cell is cell next to the boundary,
+ which is influenced by the boundary cell i.e. has a link from the boundary cell to
+ itself.
+ :param directionSymbol: a sympy symbol that can be used as index to the pdfField. It describes
+ the direction pointing from the fluid to the boundary cell
+ :param lbMethod: an instance of the LB method used. Use this to adapt the boundary to the method
+ (e.g. compressiblity)
+ :param indexField: the boundary index field that can be used to retrieve and update boundary data
+ :return: list of sympy equations
+ """
+ raise NotImplementedError("Boundary class has to overwrite __call__")
+
+ @property
+ def additionalData(self):
+ return []
+
+ @property
+ def name(self):
+ return type(self).__name__
+
+ def additionalDataInit(self, idxArray):
+ return []
+
+
+class NoSlip(Boundary):
+ """No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
+ def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+ neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
+ inverseDir = invDir(directionSymbol)
+ return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(directionSymbol))]
+
+ def __hash__(self):
+ # All boundaries of these class behave equal -> should also be equal
+ return hash("NoSlip")
+
+ def __eq__(self, other):
+ return type(other) == NoSlip
+
+
+class NoSlipFullWay(Boundary):
+ """Full-way bounce back"""
+
+ @property
+ def additionalData(self):
+ return [('lastValue', createTypeFromString("double"))]
+
+ def additionalDataInit(self, pdfField, directionSymbol, indexField, **kwargs):
+ return [sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
+
+ def __call__(self, pdfField, directionSymbol, lbMethod, indexField, **kwargs):
+ neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
+ inverseDir = invDir(directionSymbol)
+ return [sp.Eq(pdfField[neighbor](inverseDir), indexField('lastValue')),
+ sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
+
+ def __hash__(self):
+ # All boundaries of these class behave equal -> should also be equal
+ return hash("NoSlipFullWay")
+
+ def __eq__(self, other):
+ return type(other) == NoSlipFullWay
+
+
+class UBB(Boundary):
-def ubb(pdfField, direction, lbMethod, velocity, adaptVelocityToForce=False):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
- 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)
+ def __init__(self, velocity, adaptVelocityToForce=False):
+ self._velocity = velocity
+ self._adaptVelocityToForce = adaptVelocityToForce
- velocity = tuple(v_i.getShifted(*neighbor) if isinstance(v_i, Field.Access) else v_i for v_i in velocity)
+ def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+ vel = self._velocity
+ direction = directionSymbol
- if adaptVelocityToForce:
- cqc = lbMethod.conservedQuantityComputation
- shiftedVelEqs = cqc.equilibriumInputEquationsFromInitValues(velocity=velocity)
- velocity = [eq.rhs for eq in shiftedVelEqs.extract(cqc.firstOrderMomentSymbols).mainEquations]
+ assert len(vel) == lbMethod.dim, \
+ "Dimension of velocity (%d) does not match dimension of LB method (%d)" % (len(vel), lbMethod.dim)
+ neighbor = offsetFromDir(direction, lbMethod.dim)
+ inverseDir = invDir(direction)
+
+ velocity = tuple(v_i.getShifted(*neighbor) if isinstance(v_i, Field.Access) else v_i for v_i in vel)
+
+ if self._adaptVelocityToForce:
+ cqc = lbMethod.conservedQuantityComputation
+ shiftedVelEqs = cqc.equilibriumInputEquationsFromInitValues(velocity=velocity)
+ velocity = [eq.rhs for eq in shiftedVelEqs.extract(cqc.firstOrderMomentSymbols).mainEquations]
+
+ c_s_sq = sp.Rational(1, 3)
+ velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
+
+ # Better alternative: in conserved value computation
+ # rename what is currently called density to "virtualDensity"
+ # provide a new quantity density, which is constant in case of incompressible models
+ if not lbMethod.conservedQuantityComputation.zeroCenteredPdfs:
+ cqc = lbMethod.conservedQuantityComputation
+ densitySymbol = sp.Symbol("rho")
+ pdfFieldAccesses = [pdfField(i) for i in range(len(lbMethod.stencil))]
+ densityEquations = cqc.outputEquationsFromPdfs(pdfFieldAccesses, {'density': densitySymbol})
+ densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
+ result = densityEquations.allEquations
+ result += [sp.Eq(pdfField[neighbor](inverseDir),
+ pdfField(direction) - velTerm * densitySymbol)]
+ return result
+ else:
+ return [sp.Eq(pdfField[neighbor](inverseDir),
+ pdfField(direction) - velTerm)]
+
+
+class FixedDensity(Boundary):
+
+ def __init__(self, density):
+ self._density = density
- c_s_sq = sp.Rational(1, 3)
- velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
+ def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+ """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(directionSymbol, lbMethod.dim)
+ inverseDir = invDir(directionSymbol)
- # Better alternative: in conserved value computation
- # rename what is currently called density to "virtualDensity"
- # provide a new quantity density, which is constant in case of incompressible models
- if not lbMethod.conservedQuantityComputation.zeroCenteredPdfs:
cqc = lbMethod.conservedQuantityComputation
- densitySymbol = sp.Symbol("rho")
- pdfFieldAccesses = [pdfField(i) for i in range(len(lbMethod.stencil))]
- densityEquations = cqc.outputEquationsFromPdfs(pdfFieldAccesses, {'density': densitySymbol})
- densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
- result = densityEquations.allEquations
- result += [sp.Eq(pdfField[neighbor](inverseDir),
- pdfField(direction) - velTerm * densitySymbol)]
- return result
- else:
- return [sp.Eq(pdfField[neighbor](inverseDir),
- pdfField(direction) - velTerm)]
-
-
-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, lbMethod.dim)
- inverseDir = invDir(direction)
-
- cqc = lbMethod.conservedQuantityComputation
- velocity = cqc.definedSymbols()['velocity']
- symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
- substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
- symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
-
- simplification = createSimplificationStrategy(lbMethod)
- symmetricEq = simplification(symmetricEq)
-
- densitySymbol = cqc.definedSymbols()['density']
-
- densityEq = cqc.equilibriumInputEquationsFromInitValues(density=density).insertSubexpressions().mainEquations[0]
- assert densityEq.lhs == densitySymbol
- transformedDensity = densityEq.rhs
-
- conditions = [(eq_i.rhs, sp.Equality(direction, i))
- for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
- eq_component = sp.Piecewise(*conditions)
-
- subExprs = [sp.Eq(eq.lhs, transformedDensity if eq.lhs == densitySymbol else eq.rhs)
- for eq in symmetricEq.subexpressions]
- return subExprs + [sp.Eq(pdfField[neighbor](inverseDir), 2 * eq_component - pdfField(direction))]
+ velocity = cqc.definedSymbols()['velocity']
+ symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
+ substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
+ symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
+
+ simplification = createSimplificationStrategy(lbMethod)
+ symmetricEq = simplification(symmetricEq)
+
+ densitySymbol = cqc.definedSymbols()['density']
+
+ density = self._density
+ densityEq = cqc.equilibriumInputEquationsFromInitValues(density=density).insertSubexpressions().mainEquations[0]
+ assert densityEq.lhs == densitySymbol
+ transformedDensity = densityEq.rhs
+
+ conditions = [(eq_i.rhs, sp.Equality(directionSymbol, i))
+ for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
+ eq_component = sp.Piecewise(*conditions)
+
+ subExprs = [sp.Eq(eq.lhs, transformedDensity if eq.lhs == densitySymbol else eq.rhs)
+ for eq in symmetricEq.subexpressions]
+ return subExprs + [sp.Eq(pdfField[neighbor](inverseDir), 2 * eq_component - pdfField(directionSymbol))]
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 148120668ab4d0699a906efcb884da25950bd591..5863b4b39aa1451608c7bb08629dc1807d5812cd 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -13,10 +13,9 @@ WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
class BoundaryHandling(object):
class BoundaryInfo(object):
- def __init__(self, name, flag, function, kernel, ast):
- self.name = name
+ def __init__(self, flag, object, kernel, ast):
self.flag = flag
- self.function = function
+ self.object = object
self.kernel = kernel
self.ast = ast
@@ -24,29 +23,24 @@ class BoundaryHandling(object):
"""
Class for managing boundary kernels
- :param pdfField: either pdf numpy array (including ghost layers), or pystencils.Field
+ :param pdfField: pdf numpy array including ghost layers
:param domainShape: domain size without ghost layers
:param lbMethod: lattice Boltzmann method
:param ghostLayers: number of ghost layers
:param target: either 'cpu' or 'gpu'
"""
- if isinstance(pdfField, np.ndarray):
- symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
- assert pdfField.shape[:-1] == tuple(d + 2*ghostLayers for d in domainShape)
- elif isinstance(pdfField, Field):
- symbolicPdfField = pdfField
- else:
- raise ValueError("pdfField has to be either a numpy array or a pystencils.Field")
+ symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
+ assert pdfField.shape[:-1] == tuple(d + 2*ghostLayers for d in domainShape)
self._symbolicPdfField = symbolicPdfField
+ self._pdfField = pdfField
self._shapeWithGhostLayers = [d + 2 * ghostLayers for d in domainShape]
self._fluidFlag = 2 ** 0
self.flagField = np.full(self._shapeWithGhostLayers, self._fluidFlag, dtype=np.int32)
self._ghostLayers = ghostLayers
self._lbMethod = lbMethod
- self._boundaryInfos = []
- self._nameToBoundary = {}
+ self._boundaryInfos = {}
self._periodicityKernels = []
self._dirty = False
@@ -66,13 +60,12 @@ class BoundaryHandling(object):
"""Flag that is set where the lattice Boltzmann update should happen"""
return self._fluidFlag
- def getFlag(self, name):
- """Flag that represents the boundary with given name. Raises KeyError if no such boundary exists."""
- return self._nameToBoundary[name].flag
+ def getFlag(self, boundaryObject):
+ """Flag that represents the given boundary."""
+ return self._boundaryInfos[boundaryObject].flag
- def getBoundaryNames(self):
- """List of names of all registered boundary conditions"""
- return [b.name for b in self._boundaryInfos]
+ def getBoundaries(self):
+ return [b.object for b in self._boundaryInfos.values()]
def setPeriodicity(self, x=False, y=False, z=False):
"""Enable periodic boundary conditions at the border of the domain"""
@@ -82,37 +75,24 @@ class BoundaryHandling(object):
self._periodicity = [x, y, z]
self._compilePeriodicityKernels()
- def hasBoundary(self, name):
+ def hasBoundary(self, boundaryObject):
"""Returns boolean indicating if a boundary with that name exists"""
- return name in self._nameToBoundary
+ return boundaryObject in self._boundaryInfos
- def setBoundary(self, function, indexExpr=None, maskArr=None, name=None):
+ def setBoundary(self, boundaryObject, indexExpr=None, maskArr=None):
"""
Sets boundary in a rectangular region (slice)
- :param function: boundary function or the string 'fluid' to remove boundary conditions
+ :param boundaryObject: boundary condition object or the string 'fluid' to remove boundary conditions
:param indexExpr: slice expression, where boundary should be set, see :mod:`pystencils.slicing`
:param maskArr: optional boolean (masked) array specifying where the boundary should be set
- :param name: name of the boundary
"""
if indexExpr is None:
indexExpr = [slice(None, None, None)] * len(self.flagField.shape)
- if function == 'fluid':
+ if boundaryObject == 'fluid':
flag = self._fluidFlag
else:
- if name is None:
- 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 name is not specified")
-
- if not self.hasBoundary(name):
- self.addBoundary(function, name)
-
- flag = self.getFlag(name)
+ flag = self.addBoundary(boundaryObject)
assert flag != self._fluidFlag
indexExpr = normalizeSlice(indexExpr, self._shapeWithGhostLayers)
@@ -123,26 +103,20 @@ class BoundaryHandling(object):
flagFieldView[maskArr] = flag
self._dirty = True
- def addBoundary(self, boundaryFunction, name=None):
+ def addBoundary(self, boundaryObject):
"""
Adds a boundary condition, i.e. reserves a flog in the flag field and returns that flag
- If a boundary with that name already exists, the existing flag is returned.
+ If the boundary already exists, the existing flag is returned.
This flag can be logicalled or'ed to the boundaryHandling.flagField
- :param boundaryFunction: boundary condition function, see :mod:`lbmpy.boundaries.boundaryconditions`
- :param name: boundaries with different name are considered different. If not given
- ```boundaryFunction.__name__`` is used
+ :param boundaryObject: boundary condition object, see :mod:`lbmpy.boundaries.boundaryconditions`
"""
- if name is None:
- name = boundaryFunction.__name__
-
- if self.hasBoundary(name):
- return self._boundaryInfos[name].flag
+ if boundaryObject in self._boundaryInfos:
+ return self._boundaryInfos[boundaryObject].flag
newIdx = len(self._boundaryInfos) + 1 # +1 because 2**0 is reserved for fluid flag
- boundaryInfo = self.BoundaryInfo(name, 2 ** newIdx, boundaryFunction, None, None)
- self._boundaryInfos.append(boundaryInfo)
- self._nameToBoundary[name] = boundaryInfo
+ boundaryInfo = self.BoundaryInfo(2 ** newIdx, boundaryObject, None, None)
+ self._boundaryInfos[boundaryObject] = boundaryInfo
self._dirty = True
return boundaryInfo.flag
@@ -157,7 +131,7 @@ class BoundaryHandling(object):
"""Run the boundary handling, all keyword args are passed through to the boundary sweeps"""
if self._dirty:
self.prepare()
- for boundary in self._boundaryInfos:
+ for boundary in self._boundaryInfos.values():
boundary.kernel(**kwargs)
for k in self._periodicityKernels:
k(**kwargs)
@@ -166,27 +140,47 @@ class BoundaryHandling(object):
"""Removes all boundaries and fills the domain with fluid"""
self.flagField.fill(self._fluidFlag)
self._dirty = False
- self._boundaryInfos = []
- self._nameToBoundary = {}
+ self._boundaryInfos = {}
def prepare(self):
"""Fills data structures to speed up the boundary handling and compiles all boundary kernels.
This is automatically done when first called. With this function this can be triggered before running."""
- for boundary in self._boundaryInfos:
+ for boundary in self._boundaryInfos.values():
assert boundary.flag != self._fluidFlag
- idxField = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
+ idxArray = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
boundary.flag, self._fluidFlag, self._ghostLayers)
- ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._lbMethod, boundary.function,
- target=self._target)
+
+ dim = self._lbMethod.dim
+
+ if boundary.object.additionalData:
+ coordinateNames = ["x", "y", "z"][:dim]
+ indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] +
+ [('dir', np.int32)] +
+ [(i[0], i[1].numpyDtype) for i in boundary.object.additionalData])
+ extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype)
+ for prop in coordinateNames + ['dir']:
+ extendedIdxField[prop] = idxArray[prop]
+
+ idxArray = extendedIdxField
+ if boundary.object.additionalDataInit:
+ initKernelAst = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod,
+ boundary.object, target='cpu',
+ createInitializationKernel=True)
+ from pystencils.cpu import makePythonFunction as makePythonCpuFunction
+ initKernel = makePythonCpuFunction(initKernelAst, {'indexField': idxArray, 'pdfs': self._pdfField})
+ initKernel()
+
+ ast = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod, boundary.object,
+ target=self._target)
boundary.ast = ast
if self._target == 'cpu':
from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
addOpenMP(ast, numThreads=self.openMP)
- boundary.kernel = makePythonCpuFunction(ast, {'indexField': idxField})
+ boundary.kernel = makePythonCpuFunction(ast, {'indexField': idxArray})
elif self._target == 'gpu':
from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
import pycuda.gpuarray as gpuarray
- idxGpuField = gpuarray.to_gpu(idxField)
+ idxGpuField = gpuarray.to_gpu(idxArray)
boundary.kernel = makePythonGpuFunction(ast, {'indexField': idxGpuField})
else:
assert False
@@ -276,19 +270,20 @@ class LbmMethodInfo(CustomCppCode):
super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
-def generateBoundaryHandling(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu'):
+def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu',
+ createInitializationKernel=False):
indexField = Field.createFromNumpyArray("indexField", indexArr)
elements = [LbmMethodInfo(lbMethod)]
dirSymbol = TypedSymbol("dir", indexArr.dtype.fields['dir'][0])
boundaryEqList = [sp.Eq(dirSymbol, indexField[0]('dir'))]
- boundaryEqList += boundaryFunctor(pdfField, dirSymbol, lbMethod)
- if type(boundaryEqList) is tuple:
- boundaryEqList, additionalNodes = boundaryEqList
- elements += boundaryEqList
- elements += additionalNodes
+ if createInitializationKernel:
+ boundaryEqList += boundaryFunctor.additionalDataInit(pdfField=pdfField, directionSymbol=dirSymbol,
+ lbMethod=lbMethod, indexField=indexField)
else:
- elements += boundaryEqList
+ boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
+ indexField=indexField)
+ elements += boundaryEqList
if target == 'cpu':
from pystencils.cpu import createIndexedKernel
diff --git a/plot2d.py b/plot2d.py
index f49f7fe46c6ebf5cc369486fc884a11cb1016db8..320c891343042ac21096f6db582b483de371fec0 100644
--- a/plot2d.py
+++ b/plot2d.py
@@ -65,9 +65,9 @@ def plotBoundaryHandling(boundaryHandling, boundaryNameToColor=None):
'noSlip': '#000000'
}
- boundaryNames = ['fluid'] + boundaryHandling.getBoundaryNames()
- flagValues = [boundaryHandling.fluidFlag] + [boundaryHandling.getFlag(n)
- for n in boundaryHandling.getBoundaryNames()]
+ boundaryNames = ['fluid'] + [b.name for b in boundaryHandling.getBoundaries()]
+ flagValues = [boundaryHandling.fluidFlag] + [boundaryHandling.getFlag(b)
+ for b in boundaryHandling.getBoundaries()]
defaultCycler = matplotlib.rcParams['axes.prop_cycle']
colorValues = [fixedColors[name] if name in fixedColors else cycle['color']
diff --git a/scenarios.py b/scenarios.py
index d13483753f945ec6c2e5f0321537ebc9ccecba58..9f6ed78545be27d4e92556d37cebb402e2300e0d 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -11,9 +11,9 @@ It is a good starting point if you are new to lbmpy.
All scenarios can be modified, for example you can create a simple channel first, then place an object in it:
->>> from lbmpy.boundaries import noSlip
+>>> from lbmpy.boundaries import NoSlip
>>> from pystencils.slicing import makeSlice
->>> scenario.boundaryHandling.setBoundary(noSlip, makeSlice[0.3:0.4, 0.0:0.3])
+>>> scenario.boundaryHandling.setBoundary(NoSlip(), makeSlice[0.3:0.4, 0.0:0.3])
Functions for scenario setup:
-----------------------------
@@ -25,12 +25,11 @@ that defines the viscosity of the fluid (valid values being between 0 and 2).
"""
import numpy as np
import sympy as sp
-from functools import partial
from pystencils.field import getLayoutOfArray, createNumpyArrayWithLayout
from pystencils.slicing import sliceFromDirection, addGhostLayers, removeGhostLayers, normalizeSlice, makeSlice
from lbmpy.creationfunctions import createLatticeBoltzmannFunction, updateWithDefaultParameters
from lbmpy.macroscopic_value_kernels import compileMacroscopicValuesGetter, compileMacroscopicValuesSetter
-from lbmpy.boundaries import BoundaryHandling, noSlip, ubb, fixedDensity
+from lbmpy.boundaries import BoundaryHandling, NoSlip, NoSlipFullWay, UBB, FixedDensity
from lbmpy.stencils import getStencil
from lbmpy.updatekernels import createPdfArray
@@ -78,18 +77,17 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={},
"""
scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel, kernelParams=kernelParams)
- myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:scenario.method.dim])
- myUbb.name = 'ubb'
+ myUbb = UBB(velocity=[lidVelocity, 0, 0][:scenario.method.dim])
dim = scenario.method.dim
scenario.boundaryHandling.setBoundary(myUbb, sliceFromDirection('N', dim))
for direction in ('W', 'E', 'S') if scenario.method.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
- scenario.boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+ scenario.boundaryHandling.setBoundary(NoSlip(), sliceFromDirection(direction, dim))
return scenario
def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, length=None, initialVelocity=None,
- optimizationParams={}, lbmKernel=None, kernelParams={}, **kwargs):
+ optimizationParams={}, lbmKernel=None, kernelParams={}, halfWayBounceBack=True, **kwargs):
"""
Creates a channel flow in x direction, which is driven by a constant force along the x axis
@@ -104,9 +102,12 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
:param optimizationParams: see :mod:`lbmpy.creationfunctions`
:param lbmKernel: a LBM function, which would otherwise automatically created
:param kernelParams: additional parameters passed to the sweep
+ :param halfWayBounceBack: determines wall boundary condition: if true half-way bounce back, if false full-way
:param kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
:return: instance of :class:`Scenario`
"""
+ wallBoundary = NoSlip() if halfWayBounceBack else NoSlipFullWay()
+
if domainSize is not None:
dim = len(domainSize)
else:
@@ -137,18 +138,18 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
boundaryHandling.setPeriodicity(True, False, False)
if dim == 2:
for direction in ('N', 'S'):
- boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+ boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
elif dim == 3:
if roundChannel:
- noSlipIdx = boundaryHandling.addBoundary(noSlip)
+ wallIdx = boundaryHandling.addBoundary(wallBoundary)
ff = boundaryHandling.flagField
yMid = ff.shape[1] // 2
zMid = ff.shape[2] // 2
y, z = np.meshgrid(range(ff.shape[1]), range(ff.shape[2]))
- ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
+ ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = wallIdx
else:
for direction in ('N', 'S', 'T', 'B'):
- boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+ boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
assert domainSize is not None
if 'forceModel' not in kwargs:
@@ -159,7 +160,7 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim=2, radius=None, length=None,
initialVelocity=None, optimizationParams={},
- lbmKernel=None, kernelParams={}, **kwargs):
+ lbmKernel=None, kernelParams={}, halfWayBounceBack=True, **kwargs):
"""
Creates a channel flow in x direction, which is driven by two pressure boundaries.
Consider using :func:`createForceDrivenChannel` which does not have artifacts an inflow and outflow.
@@ -175,9 +176,12 @@ def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim
:param optimizationParams: see :mod:`lbmpy.creationfunctions`
:param lbmKernel: a LBM function, which would otherwise automatically created
:param kernelParams: additional parameters passed to the sweep
+ :param halfWayBounceBack: determines wall boundary condition: if true half-way bounce back, if false full-way
:param kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
:return: instance of :class:`Scenario`
"""
+ wallBoundary = NoSlip() if halfWayBounceBack else NoSlipFullWay()
+
if domainSize is not None:
dim = len(domainSize)
else:
@@ -202,20 +206,18 @@ def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim
scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel,
initialVelocity=initialVelocity, kernelParams=kernelParams)
boundaryHandling = scenario.boundaryHandling
- pressureBoundaryInflow = partial(fixedDensity, density=1.0 + pressureDifference)
- pressureBoundaryInflow.__name__ = "Inflow"
+ pressureBoundaryInflow = FixedDensity(density=1.0 + pressureDifference)
- pressureBoundaryOutflow = partial(fixedDensity, density=1.0)
- pressureBoundaryOutflow.__name__ = "Outflow"
+ pressureBoundaryOutflow = FixedDensity(density=1.0)
boundaryHandling.setBoundary(pressureBoundaryInflow, sliceFromDirection('W', dim))
boundaryHandling.setBoundary(pressureBoundaryOutflow, sliceFromDirection('E', dim))
if dim == 2:
for direction in ('N', 'S'):
- boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+ boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
elif dim == 3:
if roundChannel:
- noSlipIdx = boundaryHandling.addBoundary(noSlip)
+ noSlipIdx = boundaryHandling.addBoundary(wallBoundary)
ff = boundaryHandling.flagField
yMid = ff.shape[1] // 2
zMid = ff.shape[2] // 2
@@ -223,7 +225,7 @@ def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim
ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
else:
for direction in ('N', 'S', 'T', 'B'):
- boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+ boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
return scenario