diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 3b84aeede7dcea506818f9f13e3c2ec5c2adaf27..f551f3aeede52ee6929d7af118902eb3055fe30e 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -196,7 +196,9 @@ class FixedDensity(Boundary):
class NeumannByCopy(Boundary):
def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, lbMethod.dim)
- return [sp.Eq(pdfField[neighbor](directionSymbol), pdfField(directionSymbol))]
+ inverseDir = BoundaryOffsetInfo.invDir(directionSymbol)
+ return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(inverseDir)),
+ sp.Eq(pdfField[neighbor](directionSymbol), pdfField(directionSymbol))]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal
@@ -205,3 +207,38 @@ class NeumannByCopy(Boundary):
def __eq__(self, other):
return type(other) == NeumannByCopy
+
+class StreamInZero(Boundary):
+ def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+ neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, lbMethod.dim)
+ inverseDir = BoundaryOffsetInfo.invDir(directionSymbol)
+ return [sp.Eq(pdfField[neighbor](inverseDir), 0),
+ sp.Eq(pdfField[neighbor](directionSymbol), 0)]
+
+ def __hash__(self):
+ # All boundaries of these class behave equal -> should also be equal
+ return hash("StreamInZero")
+
+ def __eq__(self, other):
+ return type(other) == StreamInZero
+
+
+class AntiBounceBack(Boundary):
+
+ """No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
+ def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+ neighbor = BoundaryOffsetInfo.offsetFromDir(directionSymbol, lbMethod.dim)
+ inverseDir = BoundaryOffsetInfo.invDir(directionSymbol)
+ rhs = pdfField(directionSymbol)
+ t = -rhs
+ return [SympyAssignment(pdfField[neighbor](inverseDir), t)]
+
+ def __hash__(self):
+ # All boundaries of these class behave equal -> should also be equal (as long as name is equal)
+ return hash(self.name)
+
+ def __eq__(self, other):
+ if not isinstance(other, AntiBounceBack):
+ return False
+ return self.name == other.name
+
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 1ab372a541d866b763316268ea92b48043d2b26a..990b17e2c1bd30d09e7129f2ad2504b6b4821664 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -9,10 +9,11 @@ from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
- def __init__(self, lbMethod, dataHandling, pdfFieldName, name="boundaryHandling", target='cpu', openMP=True):
+ def __init__(self, lbMethod, dataHandling, pdfFieldName, name="boundaryHandling", flagInterface=None,
+ target='cpu', openMP=True):
self.lbMethod = lbMethod
super(LatticeBoltzmannBoundaryHandling, self).__init__(dataHandling, pdfFieldName, lbMethod.stencil,
- name, target, openMP)
+ name, flagInterface, target, openMP)
def forceOnBoundary(self, boundaryObj):
from lbmpy.boundaries import NoSlip
@@ -26,7 +27,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
def _forceOnNoSlip(self, boundaryObj):
dh = self._dataHandling
- ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
+ ffGhostLayers = dh.ghostLayersOfField(self.flagInterface.flagFieldName)
method = self.lbMethod
stencil = np.array(method.stencil)
@@ -46,7 +47,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
def _forceOnBoundary(self, boundaryObj):
dh = self._dataHandling
- ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
+ ffGhostLayers = dh.ghostLayersOfField(self.flagInterface.flagFieldName)
method = self.lbMethod
stencil = np.array(method.stencil)
invDirection = np.array([method.stencil.index(inverseDirection(d))
diff --git a/geometry.py b/geometry.py
index 5bbd3aaa15f90c2103c376c2e97a389a4cc05018..9dfe7bd07778382bf1a1dc5d54711dfc231a7c37 100644
--- a/geometry.py
+++ b/geometry.py
@@ -25,12 +25,13 @@ def getParabolicInitialVelocity(domainSize, u_max, velCoord=0, diameter=None):
return u
-def addBox(boundaryHandling, boundary=NoSlip()):
+def addBox(boundaryHandling, boundary=NoSlip(), replace=True):
borders = ['N', 'S', 'E', 'W']
if boundaryHandling.dim == 3:
borders += ['T', 'B']
for d in borders:
- boundaryHandling.setBoundary(boundary, sliceFromDirection(d, boundaryHandling.dim))
+ flag = boundaryHandling.setBoundary(boundary, sliceFromDirection(d, boundaryHandling.dim), replace=replace)
+ return flag
def addParabolicVelocityInflow(boundaryHandling, u_max, indexExpr, velCoord=0, diameter=None):
diff --git a/lbstep.py b/lbstep.py
index 21da0682985ada3e394788ff08dc58285fab1c00..46baa5b05f734c085ef94201f3dffd9c5ae68fbf 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -16,7 +16,7 @@ class LatticeBoltzmannStep:
kernelParams={}, dataHandling=None, name="lbm", optimizationParams={},
velocityDataName=None, densityDataName=None, densityDataIndex=None,
computeVelocityInEveryStep=False, computeDensityInEveryStep=False,
- velocityInputArrayName=None, timeStepOrder='streamCollide',
+ velocityInputArrayName=None, timeStepOrder='streamCollide', flagInterface=None,
**methodParameters):
# --- Parameter normalization ---
@@ -103,6 +103,7 @@ class LatticeBoltzmannStep:
self._sync = dataHandling.synchronizationFunction([self._pdfArrName], methodParameters['stencil'], target)
self._boundaryHandling = LatticeBoltzmannBoundaryHandling(self.method, self._dataHandling, self._pdfArrName,
name=name + "_boundaryHandling",
+ flagInterface=flagInterface,
target=target, openMP=optimizationParams['openMP'])
# -- Macroscopic Value Kernels
@@ -154,7 +155,7 @@ class LatticeBoltzmannStep:
return
if masked:
- mask = self.boundaryHandling.getMask(sliceObj[:self.dim], 'fluid', True)
+ mask = self.boundaryHandling.getMask(sliceObj[:self.dim], 'domain', True)
if len(mask.shape) < len(result.shape):
assert len(mask.shape) + 1 == len(result.shape)
mask = np.repeat(mask[..., np.newaxis], result.shape[-1], axis=2)
@@ -195,7 +196,7 @@ class LatticeBoltzmannStep:
self._sync()
self._boundaryHandling(**self.kernelParams)
self._dataHandling.runKernel(self._lbmKernels[1], **self.kernelParams)
- else: # stream collide
+ else: # stream collide
self._sync()
self._boundaryHandling(**self.kernelParams)
self._dataHandling.runKernel(self._lbmKernels[0], **self.kernelParams)
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index 10400fde2456ce4031a6917ed434ca3e98b31dea..56f84746ba0d5a9441d808b309f79fad2e34bb57 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -62,7 +62,7 @@ def cahnHilliardFdEq(phaseIdx, phi, mu, velocity, mobility, dx, dt):
class CahnHilliardFDStep:
def __init__(self, dataHandling, phiFieldName, muFieldName, velocityFieldName, name='ch_fd', target='cpu',
- dx=1, dt=1, mobilities=1):
+ dx=1, dt=1, mobilities=1, equationModifier=lambda eqs: eqs):
from pystencils import createKernel
self.dataHandling = dataHandling
@@ -80,13 +80,14 @@ class CahnHilliardFDStep:
rhs = cahnHilliardFdEq(i, self.phiField, muField, velField, mobilities[i], dx, dt)
updateEqs.append(sp.Eq(self.tmpField(i), rhs))
self.updateEqs = updateEqs
- self.kernel = createKernel(updateEqs, target=target).compile()
+ self.updateEqs = equationModifier(updateEqs)
+ self.kernel = createKernel(self.updateEqs, target=target).compile()
self.sync = self.dataHandling.synchronizationFunction([phiFieldName, velocityFieldName, muFieldName],
target=target)
- def timeStep(self):
+ def timeStep(self, **kwargs):
self.sync()
- self.dataHandling.runKernel(self.kernel)
+ self.dataHandling.runKernel(self.kernel, **kwargs)
self.dataHandling.swap(self.phiField.name, self.tmpField.name)
def setPdfFieldsFromMacroscopicValues(self):
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index ed3043be0b628b500c219f1048eb3462b2554934..dc52d57cce75b9d6c94ad96df8a55c4ebbc3dc46 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -1,11 +1,14 @@
import sympy as sp
import numpy as np
+
from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.phasefield.cahn_hilliard_lbm import cahnHilliardLbmMethod
-from lbmpy.phasefield.kerneleqs import muKernel, forceKernelUsingMu, CahnHilliardFDStep, pressureTensorKernel, \
+from lbmpy.phasefield.kerneleqs import muKernel, CahnHilliardFDStep, pressureTensorKernel, \
forceKernelUsingPressureTensor
from pystencils import createKernel
from lbmpy.phasefield.analytical import chemicalPotentialsFromFreeEnergy, symmetricTensorLinearization
+from pystencils.boundaries.boundaryhandling import FlagInterface
+from pystencils.boundaries.inkernel import addNeumannBoundary
from pystencils.datahandling import SerialDataHandling
from pystencils.equationcollection.simplifications import sympyCseOnEquationList
from pystencils.slicing import makeSlice, SlicedGetter
@@ -17,7 +20,8 @@ class PhaseFieldStep:
name='pf', hydroLbmParameters={},
hydroDynamicRelaxationRate=1.0, cahnHilliardRelaxationRates=1.0, densityOrderParameter=None,
optimizationParams=None, kernelParams={}, dx=1, dt=1, solveCahnHilliardWithFiniteDifferences=False,
- orderParameterForce=None, concentrationToOrderParameters=None, orderParametersToConcentrations=None):
+ orderParameterForce=None, concentrationToOrderParameters=None, orderParametersToConcentrations=None,
+ activateHomogenousNeumannBoundaries=False):
if optimizationParams is None:
optimizationParams = {'openMP': False, 'target': 'cpu'}
@@ -50,17 +54,31 @@ class PhaseFieldStep:
self.forceField = dataHandling.addArray(self.forceFieldName, fSize=dataHandling.dim, gpu=gpu, latexName="F")
self.pressureTensorField = dataHandling.addArray(self.pressureTensorFieldName,
fSize=pressureTensorSize, latexName='P')
+ self.flagInterface = FlagInterface(dataHandling, 'flags')
# ------------------ Creating kernels ------------------
phi = tuple(self.phiField(i) for i in range(len(orderParameters)))
F = self.freeEnergy.subs({old: new for old, new in zip(orderParameters, phi)})
+ if activateHomogenousNeumannBoundaries:
+ def applyNeumannBoundaries(eqs):
+ fields = [dataHandling.fields[self.phiFieldName],
+ dataHandling.fields[self.pressureTensorFieldName],
+ ]
+ flagField = dataHandling.fields[self.flagInterface.flagFieldName]
+ return addNeumannBoundary(eqs, fields, flagField, "neumannFlag", inverseFlag=False)
+ else:
+ def applyNeumannBoundaries(eqs):
+ return eqs
+
# μ and pressure tensor update
self.phiSync = dataHandling.synchronizationFunction([self.phiFieldName], target=target)
self.muEqs = muKernel(F, phi, self.phiField, self.muField, dx)
self.pressureTensorEqs = pressureTensorKernel(self.freeEnergy, orderParameters,
self.phiField, self.pressureTensorField, dx)
- self.muAndPressureTensorKernel = createKernel(sympyCseOnEquationList(self.muEqs + self.pressureTensorEqs),
+ muAndPressureTensorEqs = self.muEqs + self.pressureTensorEqs
+ muAndPressureTensorEqs = applyNeumannBoundaries(muAndPressureTensorEqs)
+ self.muAndPressureTensorKernel = createKernel(sympyCseOnEquationList(muAndPressureTensorEqs),
target=target, cpuOpenMP=openMP).compile()
# F Kernel
@@ -70,7 +88,8 @@ class PhaseFieldStep:
extraForce += self.phiField(orderParameterIdx) * sp.Matrix(force)
self.forceEqs = forceKernelUsingPressureTensor(self.forceField, self.pressureTensorField, dx=dx,
extraForce=extraForce)
- self.forceFromPressureTensorKernel = createKernel(self.forceEqs, target=target, cpuOpenMP=openMP).compile()
+ self.forceFromPressureTensorKernel = createKernel(applyNeumannBoundaries(self.forceEqs),
+ target=target, cpuOpenMP=openMP).compile()
self.pressureTensorSync = dataHandling.synchronizationFunction([self.pressureTensorFieldName], target=target)
# Hydrodynamic LBM
@@ -89,6 +108,7 @@ class PhaseFieldStep:
relaxationRate=hydroDynamicRelaxationRate,
computeVelocityInEveryStep=True, force=self.forceField,
velocityDataName=self.velFieldName, kernelParams=kernelParams,
+ flagInterface=self.flagInterface,
timeStepOrder='collideStream',
**hydroLbmParameters)
@@ -103,7 +123,8 @@ class PhaseFieldStep:
raise NotImplementedError("densityOrderParameter not supported when "
"CH is solved with finite differences")
chStep = CahnHilliardFDStep(self.dataHandling, self.phiFieldName, self.muFieldName, self.velFieldName,
- target=target, dx=dx, dt=dt, mobilities=1)
+ target=target, dx=dx, dt=dt, mobilities=1,
+ equationModifier=applyNeumannBoundaries)
self.cahnHilliardSteps.append(chStep)
else:
for i, op in enumerate(orderParameters):
@@ -118,6 +139,7 @@ class PhaseFieldStep:
stencil='D3Q19' if self.dataHandling.dim == 3 else 'D2Q9',
computeDensityInEveryStep=True,
densityDataIndex=i,
+ flagInterface=self.hydroLbmStep.boundaryHandling.flagInterface,
name=name + "_chLbm_%d" % (i,),
optimizationParams=optimizationParams)
self.cahnHilliardSteps.append(chStep)
@@ -130,6 +152,8 @@ class PhaseFieldStep:
self.timeStepsRun = 0
self.reset()
+ self.neumannFlag = 0
+
def writeVTK(self):
self.vtkWriter(self.timeStepsRun)
@@ -169,19 +193,32 @@ class PhaseFieldStep:
chStep.postRun()
def timeStep(self):
+ neumannFlag = self.neumannFlag
+ #for b in self.dataHandling.iterate(sliceObj=makeSlice[:, 0]):
+ # b[self.phiFieldName][..., 0] = 0.0
+ # b[self.phiFieldName][..., 1] = 1.0
+ #for b in self.dataHandling.iterate(sliceObj=makeSlice[0, :]):
+ # b[self.phiFieldName][..., 0] = 1.0
+ # b[self.phiFieldName][..., 1] = 0.0
+
self.phiSync()
- self.dataHandling.runKernel(self.muAndPressureTensorKernel)
+ self.dataHandling.runKernel(self.muAndPressureTensorKernel, neumannFlag=neumannFlag)
self.pressureTensorSync()
- self.dataHandling.runKernel(self.forceFromPressureTensorKernel)
+ self.dataHandling.runKernel(self.forceFromPressureTensorKernel, neumannFlag=neumannFlag)
if self.runHydroLbm:
self.hydroLbmStep.timeStep()
for chLbm in self.cahnHilliardSteps:
+ #chLbm.timeStep(neumannFlag=neumannFlag)
chLbm.timeStep()
self.timeStepsRun += 1
+ @property
+ def boundaryHandling(self):
+ return self.hydroLbmStep.boundaryHandling
+
def setConcentration(self, sliceObj, concentration):
if self.concentrationToOrderParameter is not None:
phi = self.concentrationToOrderParameter(concentration)
@@ -242,4 +279,4 @@ class PhaseFieldStep:
@property
def force(self):
- return SlicedGetter(self.forceSlice)
\ No newline at end of file
+ return SlicedGetter(self.forceSlice)
diff --git a/scenarios.py b/scenarios.py
index 4fe2ad872a1bdd7905e0d9fb96e0857f161d88b7..b47aa9ee3f0b0049d4bfc62b83cd16ce499df89b 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -13,7 +13,9 @@ All scenarios can be modified, for example you can create a simple channel first
>>> from lbmpy.boundaries import NoSlip
>>> from pystencils.slicing import makeSlice
->>> scenario.boundaryHandling.setBoundary(NoSlip(), makeSlice[0.3:0.4, 0.0:0.3])
+>>> flag = scenario.boundaryHandling.setBoundary(NoSlip(), makeSlice[0.3:0.4, 0.0:0.3])
+
+
Functions for scenario setup:
---- -------------------------