From f9532aae7717e156e3198375f1518da45b2c19c6 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 2 Mar 2018 10:22:05 +0100
Subject: [PATCH] Boundary conditions

- in-kernel Neumann boundaries
- flag-interface for boundary handling makes one flag field multiple
  boundary handlings possible
- generator: support for bitwise logical operators
---
 boundaries/boundaryconditions.py | 39 ++++++++++++++++++++++-
 boundaries/boundaryhandling.py   |  9 +++---
 geometry.py                      |  5 +--
 lbstep.py                        |  7 +++--
 phasefield/kerneleqs.py          |  9 +++---
 phasefield/phasefieldstep.py     | 53 +++++++++++++++++++++++++++-----
 scenarios.py                     |  4 ++-
 7 files changed, 103 insertions(+), 23 deletions(-)

diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 3b84aeed..f551f3ae 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 1ab372a5..990b17e2 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 5bbd3aaa..9dfe7bd0 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 21da0682..46baa5b0 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 10400fde..56f84746 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 ed3043be..dc52d57c 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 4fe2ad87..b47aa9ee 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:
 ----    -------------------------
-- 
GitLab