From 363112f3e0e27b7dcde62509d08601bbd24ba649 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Thu, 18 Jan 2018 09:11:41 +0100
Subject: [PATCH] New Boundary step system / fixes in lbmpy phasefield
- scaling interface width eta instead of surface tensions tau to correct
interface profile & surface tensions
---
boundaries/__init__.py | 2 +-
boundaries/boundary_kernel.py | 24 ++-
boundaries/boundaryconditions.py | 61 ++++----
boundaries/createindexlist.py | 16 ++
boundaries/createindexlistcython.pyx | 3 -
boundaries/handlinginterface.py | 2 +-
boundary_handling.py | 221 +++++++++++++++++++++++++++
creationfunctions.py | 5 +-
lbstep.py | 146 ++++++++++++++++++
phasefield/analytical.py | 9 +-
phasefield/postprocessing.py | 5 +-
11 files changed, 437 insertions(+), 57 deletions(-)
create mode 100644 boundary_handling.py
create mode 100644 lbstep.py
diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index a07c43c5..7588d247 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -1,2 +1,2 @@
-from lbmpy.boundaries.boundaryconditions import NoSlip, NoSlipFullWay, UBB, FixedDensity, NeumannByCopy
+from lbmpy.boundaries.boundaryconditions import NoSlip, UBB, FixedDensity, NeumannByCopy
from lbmpy.boundaries.boundaryhandling import BoundaryHandling
\ No newline at end of file
diff --git a/boundaries/boundary_kernel.py b/boundaries/boundary_kernel.py
index ed6f4007..d9e0f89c 100644
--- a/boundaries/boundary_kernel.py
+++ b/boundaries/boundary_kernel.py
@@ -2,6 +2,7 @@ import sympy as sp
import numpy as np
from pystencils.backends.cbackend import CustomCppCode
from pystencils import TypedSymbol, Field
+from pystencils.cpu import addOpenMP
from pystencils.data_types import castFunc, createType
INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
@@ -51,30 +52,25 @@ class LbmMethodInfo(CustomCppCode):
super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
-def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu',
- createInitializationKernel=False):
+def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu', openMP=True):
indexField = Field.createFromNumpyArray("indexField", indexArr)
- return generateIndexBoundaryKernelGeneric(pdfField, indexField, indexArr.dtype, lbMethod, boundaryFunctor, target,
- createInitializationKernel)
+ return generateIndexBoundaryKernelGeneric(pdfField, indexField, lbMethod, boundaryFunctor, target, openMP)
-def generateIndexBoundaryKernelGeneric(pdfField, indexField, indexArrDtype, lbMethod, boundaryFunctor, target='cpu',
- createInitializationKernel=False):
-
+def generateIndexBoundaryKernelGeneric(pdfField, indexField, lbMethod, boundaryFunctor, target='cpu', openMP=True):
elements = [LbmMethodInfo(lbMethod)]
+ indexArrDtype = indexField.dtype.numpyDtype
dirSymbol = TypedSymbol("dir", indexArrDtype.fields['dir'][0])
boundaryEqList = [sp.Eq(dirSymbol, indexField[0]('dir'))]
- if createInitializationKernel:
- boundaryEqList += boundaryFunctor.additionalDataInitKernelEquations(pdfField=pdfField, directionSymbol=dirSymbol,
- lbMethod=lbMethod, indexField=indexField)
- else:
- boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
- indexField=indexField)
+ boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
+ indexField=indexField)
elements += boundaryEqList
if target == 'cpu':
from pystencils.cpu import createIndexedKernel
- return createIndexedKernel(elements, [indexField])
+ ast = createIndexedKernel(elements, [indexField])
+ addOpenMP(ast, numThreads=openMP)
+ return ast
elif target == 'gpu':
from pystencils.gpucuda import createdIndexedCUDAKernel
return createdIndexedCUDAKernel(elements, [indexField])
\ No newline at end of file
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index f67d31cd..9557acc5 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -39,13 +39,8 @@ class Boundary(object):
@property
def additionalDataInitCallback(self):
- """Return a callback function called with (x, y, [z]), and returning a dict of data-name to data for each
- element that should be initialized"""
- return None
-
- @property
- def additionalDataInitKernelEquations(self):
- """Should return callback that returns kernel equations for the boundary data initialization kernel"""
+ """Return a callback function called with a boundary data setter object and returning a dict of
+ data-name to data for each element that should be initialized"""
return None
@property
@@ -82,31 +77,33 @@ class NoSlip(Boundary):
return self.name == other.name
-class NoSlipFullWay(Boundary):
- """Full-way bounce back"""
-
- @property
- def additionalData(self):
- return [('lastValue', createType("double"))]
-
- @property
- def additionalDataInitKernelEquations(self):
- def kernelEqGetter(pdfField, directionSymbol, indexField, **kwargs):
- return [sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
- return kernelEqGetter
-
- 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 NoSlipFullWay(Boundary):
+# """Full-way bounce back"""
+#
+# @property
+# def additionalData(self):
+# return [('lastValue', createType("double"))]
+#
+# @property
+# def additionalDataInitKernelEquations(self):
+# # TODO this function is not longer availabel
+# # formulate as additionalDataInitCallback
+# def kernelEqGetter(pdfField, directionSymbol, indexField, **kwargs):
+# return [sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
+# return kernelEqGetter
+#
+# 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):
diff --git a/boundaries/createindexlist.py b/boundaries/createindexlist.py
index 69c5cb4d..aa217305 100644
--- a/boundaries/createindexlist.py
+++ b/boundaries/createindexlist.py
@@ -60,3 +60,19 @@ def createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGho
if flagField.size > 1e6:
warnings.warn("Boundary setup may take very long! Consider installing cython to speed it up")
return _createBoundaryIndexListPython(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)
+
+
+def createBoundaryIndexArray(flagField, stencil, boundaryMask, fluidMask, boundaryObject, nrOfGhostLayers=1):
+ idxArray = createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGhostLayers)
+ dim = len(flagField.shape)
+
+ if boundaryObject.additionalData:
+ coordinateNames = boundaryIndexArrayCoordinateNames[:dim]
+ indexArrDtype = numpyDataTypeForBoundaryObject(boundaryObject, dim)
+ extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype)
+ for prop in coordinateNames + ['dir']:
+ extendedIdxField[prop] = idxArray[prop]
+
+ idxArray = extendedIdxField
+
+ return idxArray
diff --git a/boundaries/createindexlistcython.pyx b/boundaries/createindexlistcython.pyx
index b6126fd1..dfa84069 100644
--- a/boundaries/createindexlistcython.pyx
+++ b/boundaries/createindexlistcython.pyx
@@ -48,9 +48,6 @@ def createBoundaryIndexList3D(object[IntegerType, ndim=3] flagField,
boundaryIndexList = []
numDirections = stencil.shape[0]
- #for x in range(nrOfGhostLayers,xs-nrOfGhostLayers):
- # for y in range(nrOfGhostLayers,ys-nrOfGhostLayers):
- # for z in range(nrOfGhostLayers, zs-nrOfGhostLayers):
for z in range(nrOfGhostLayers, zs-nrOfGhostLayers):
for y in range(nrOfGhostLayers,ys-nrOfGhostLayers):
for x in range(nrOfGhostLayers,xs-nrOfGhostLayers):
diff --git a/boundaries/handlinginterface.py b/boundaries/handlinginterface.py
index 0100cd98..4653feda 100644
--- a/boundaries/handlinginterface.py
+++ b/boundaries/handlinginterface.py
@@ -223,7 +223,7 @@ class GenericBoundaryHandling(object):
self._flagFieldInterface.getFlag(boundaryObject)
def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, includeGhostLayers=True):
- """c
+ """
Sets boundary using either a rectangular slice, a boolean mask or a combination of both
:param boundaryObject: instance of a boundary object that should be set
diff --git a/boundary_handling.py b/boundary_handling.py
new file mode 100644
index 00000000..f981dbee
--- /dev/null
+++ b/boundary_handling.py
@@ -0,0 +1,221 @@
+import numpy as np
+
+from lbmpy.boundaries.boundary_kernel import generateIndexBoundaryKernelGeneric
+from pystencils import Field
+from lbmpy.boundaries.createindexlist import createBoundaryIndexArray, numpyDataTypeForBoundaryObject
+from pystencils.cache import memorycache
+
+
+class BoundaryDataSetter:
+
+ def __init__(self, indexArray, offset, stencil, ghostLayers, pdfArray):
+ self.indexArray = indexArray
+ self.offset = offset
+ self.stencil = np.array(stencil)
+ self.pdfArray = pdfArray.view()
+ self.pdfArray.flags.writeable = False
+
+ arrFieldNames = indexArray.dtype.names
+ self.dim = 3 if 'z' in arrFieldNames else 2
+ assert 'x' in arrFieldNames and 'y' in arrFieldNames and 'dir' in arrFieldNames, str(arrFieldNames)
+ self.boundaryDataNames = set(self.indexArray.dtype.names) - set(['x', 'y', 'z', 'dir'])
+ self.coordMap = {0: 'x', 1: 'y', 2: 'z'}
+ self.ghostLayers = ghostLayers
+
+ def fluidCellPositions(self, coord):
+ assert coord < self.dim
+ return self.indexArray[self.coordMap[coord]] + self.offset[coord] - self.ghostLayers + 0.5
+
+ @memorycache()
+ def linkOffsets(self):
+ return self.stencil[self.indexArray['dir']]
+
+ @memorycache()
+ def linkPositions(self, coord):
+ return self.fluidCellPositions(coord) + 0.5 * self.linkOffsets()[:, coord]
+
+ @memorycache()
+ def boundaryCellPositions(self, coord):
+ return self.fluidCellPositions(coord) + self.linkOffsets()[:, coord]
+
+ def __setitem__(self, key, value):
+ if key not in self.boundaryDataNames:
+ raise KeyError("Invalid boundary data name %s. Allowed are %s" % (key, self.boundaryDataNames))
+ self.indexArray[key] = value
+
+ def __getitem__(self, item):
+ if item not in self.boundaryDataNames:
+ raise KeyError("Invalid boundary data name %s. Allowed are %s" % (item, self.boundaryDataNames))
+ return self.indexArray[item]
+
+
+class BoundaryHandling:
+
+ def __init__(self, lbMethod, dataHandling, pdfFieldName, name="boundaryHandling", target='cpu', openMP=True):
+ assert dataHandling.hasData(pdfFieldName)
+
+ self._lbMethod = lbMethod
+ self._dataHandling = dataHandling
+ self._pdfFieldName = pdfFieldName
+ self._flagFieldName = name + "Flags"
+ self._indexArrayName = name + "IndexArrays"
+ self._target = target
+ self._openMP = openMP
+ self._boundaryObjectToBoundaryInfo = {}
+ self._fluidFlag = 1 << 0
+ self._nextFreeFlag = 1
+
+ self._dirty = True
+
+ # Add flag field to data handling if it does not yet exist
+ if dataHandling.hasData(self._flagFieldName) or dataHandling.hasData(self._indexArrayName):
+ raise ValueError("There is already a boundary handling registered at the data handling."
+ "If you want to add multiple handlings, choose a different name.")
+
+ gpu = self._target == 'gpu'
+ dataHandling.addArray(self._flagFieldName, dtype=np.int32, cpu=True, gpu=gpu)
+ dataHandling.addCustomClass(self._indexArrayName, self.IndexFieldBlockData, cpu=True, gpu=gpu)
+
+ ffGhostLayers = self._dataHandling.ghostLayersOfField(self._flagFieldName)
+ for b in self._dataHandling.iterate(ghostLayers=ffGhostLayers):
+ b[self._flagFieldName].fill(self._fluidFlag)
+
+ @property
+ def dim(self):
+ return self._lbMethod.dim
+
+ @property
+ def lbMethod(self):
+ return self._lbMethod
+
+ @property
+ def boundaryObjects(self):
+ return tuple(self._boundaryObjectToName.keys())
+
+ def setBoundary(self, boundaryObject, sliceObj=None, maskCallback=None, ghostLayers=True):
+ """
+ Sets boundary using either a rectangular slice, a boolean mask or a combination of both
+
+ :param boundaryObject: instance of a boundary object that should be set
+ :param sliceObj: a slice object (can be created with makeSlice[]) that selects a part of the domain where
+ the boundary should be set. If none, the complete domain is selected which makes only sense
+ if a maskCallback is passed. The slice can have ':' placeholders, which are interpreted
+ depending on the 'includeGhostLayers' parameter i.e. if it is True, the slice extends
+ into the ghost layers
+ :param maskCallback: callback function getting x,y (z) parameters of the cell midpoints and returning a
+ boolean mask with True entries where boundary cells should be set.
+ The x, y, z arrays have 2D/3D shape such that they can be used directly
+ to create the boolean return array. i.e return x < 10 sets boundaries in cells with
+ midpoint x coordinate smaller than 10.
+ :param ghostLayers see DataHandling.iterate()
+ """
+ flag = self._getFlagForBoundary(boundaryObject)
+
+ for b in self._dataHandling.iterate(sliceObj=sliceObj, ghostLayers=ghostLayers):
+ flagArr = b[self._flagFieldName]
+ if maskCallback is not None:
+ mask = maskCallback(*b.midpointArrays)
+ flagArr[mask] = flag
+ else:
+ flagArr.fill(flag)
+
+ self._dirty = True
+
+ def getBoundaryNameToFlagDict(self):
+ result = {bObj.name: bInfo.flag for bObj, bInfo in self._boundaryObjectToBoundaryInfo.values()}
+ result['fluid'] = self._fluidFlag
+ return result
+
+ def prepare(self):
+ if not self._dirty:
+ return
+ self._createIndexFields()
+ self._dirty = False
+
+ def triggerReinitializationOfBoundaryData(self, **kwargs):
+ if self._dirty:
+ self.prepare()
+ else:
+ ffGhostLayers = self._dataHandling.ghostLayersOfField(self._flagFieldName)
+ for b in self._dataHandling.iterate(ghostLayers=ffGhostLayers):
+ for bObj, setter in b[self._indexArrayName].boundaryObjectToDataSetter.items():
+ self._boundaryDataInitialization(bObj, setter, **kwargs)
+
+ def __call__(self, **kwargs):
+ if self._dirty:
+ self.prepare()
+
+ for b in self._dataHandling.iterate():
+ for bInfo in self._boundaryObjectToBoundaryInfo.values():
+ idxArr = b[self._indexArrayName].boundaryObjectToIndexList[bInfo.boundaryObject]
+ kwargs[self._pdfFieldName] = b[self._pdfFieldName]
+ bInfo.kernel(indexField=idxArr, **kwargs)
+
+ # ------------------------------ Implementation Details ------------------------------------------------------------
+
+ def _getFlagForBoundary(self, boundaryObject):
+ if boundaryObject not in self._boundaryObjectToBoundaryInfo:
+ symbolicIndexField = Field.createGeneric('indexField', spatialDimensions=1,
+ dtype=numpyDataTypeForBoundaryObject(boundaryObject, self.dim))
+
+ spPdfField = self._dataHandling.fields[self._pdfFieldName]
+ ast = generateIndexBoundaryKernelGeneric(spPdfField, symbolicIndexField, self._lbMethod,
+ boundaryObject, target=self._target, openMP=self._openMP)
+
+ boundaryInfo = self.BoundaryInfo(boundaryObject, flag=1 << self._nextFreeFlag, kernel=ast.compile())
+
+ self._nextFreeFlag += 1
+ self._boundaryObjectToBoundaryInfo[boundaryObject] = boundaryInfo
+ return self._boundaryObjectToBoundaryInfo[boundaryObject].flag
+
+ def _createIndexFields(self):
+ dh = self._dataHandling
+ ffGhostLayers = dh.ghostLayersOfField(self._flagFieldName)
+ for b in dh.iterate(ghostLayers=ffGhostLayers):
+ flagArr = b[self._flagFieldName]
+ for bInfo in self._boundaryObjectToBoundaryInfo.values():
+ idxArr = createBoundaryIndexArray(flagArr, self.lbMethod.stencil, bInfo.flag, self._fluidFlag,
+ bInfo.boundaryObject, dh.ghostLayersOfField(self._flagFieldName))
+
+ pdfArr = b[self._pdfFieldName]
+ # TODO test that offset is used correctly here
+ boundaryDataSetter = BoundaryDataSetter(idxArr, b.offset, self.lbMethod.stencil, ffGhostLayers, pdfArr)
+ indexArrayBD = b[self._indexArrayName]
+ indexArrayBD.boundaryObjectToIndexList[bInfo.boundaryObject] = idxArr
+ indexArrayBD.boundaryObjectToDataSetter[bInfo.boundaryObject] = boundaryDataSetter
+ self._boundaryDataInitialization(bInfo.boundaryObject, boundaryDataSetter)
+
+ def _boundaryDataInitialization(self, boundaryObject, boundaryDataSetter, **kwargs):
+ if boundaryObject.additionalDataInitCallback:
+ boundaryObject.additionalDataInitCallback(boundaryDataSetter, **kwargs)
+ if self._target == 'gpu':
+ self._dataHandling.toGpu(self._indexArrayName)
+
+ class BoundaryInfo(object):
+ def __init__(self, boundaryObject, flag, kernel):
+ self.boundaryObject = boundaryObject
+ self.flag = flag
+ self.kernel = kernel
+
+ class IndexFieldBlockData:
+ def __init__(self):
+ self.boundaryObjectToIndexList = {}
+ self.boundaryObjectToDataSetter = {}
+
+ @staticmethod
+ def toCpu(gpuVersion, cpuVersion):
+ gpuVersion = gpuVersion.boundaryObjectToIndexList
+ cpuVersion = cpuVersion.boundaryObjectToIndexList
+ for obj, cpuArr in cpuVersion.values():
+ gpuVersion[obj].get(cpuArr)
+
+ @staticmethod
+ def toGpu(gpuVersion, cpuVersion):
+ from pycuda import gpuarray
+ gpuVersion = gpuVersion.boundaryObjectToIndexList
+ cpuVersion = cpuVersion.boundaryObjectToIndexList
+ for obj, cpuArr in cpuVersion.values():
+ if obj not in gpuVersion:
+ gpuVersion[obj] = gpuarray.to_gpu(cpuArr)
+ else:
+ gpuVersion[obj].set(cpuArr)
\ No newline at end of file
diff --git a/creationfunctions.py b/creationfunctions.py
index 4e01a01d..d83f0f89 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -197,6 +197,7 @@ def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
'fieldSize': None,
'fieldLayout': 'fzyx', # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
+ 'symbolicField': None,
'target': 'cpu',
'openMP': True,
@@ -348,7 +349,9 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
else:
fieldDtype = 'float64'
- if optParams['fieldSize']:
+ if optParams['symbolicField'] is not None:
+ srcField = optParams['symbolicField']
+ elif optParams['fieldSize']:
fieldSize = [s + 2 for s in optParams['fieldSize']] + [len(stencil)]
srcField = Field.createFixedSize(params['fieldName'], fieldSize, indexDimensions=1,
layout=optParams['fieldLayout'], dtype=fieldDtype)
diff --git a/lbstep.py b/lbstep.py
new file mode 100644
index 00000000..496dc111
--- /dev/null
+++ b/lbstep.py
@@ -0,0 +1,146 @@
+from lbmpy.boundary_handling import BoundaryHandling
+from lbmpy.creationfunctions import switchToSymbolicRelaxationRatesForEntropicMethods, createLatticeBoltzmannFunction, \
+ updateWithDefaultParameters
+from lbmpy.simplificationfactory import createSimplificationStrategy
+from lbmpy.stencils import getStencil
+from pystencils.datahandling import SerialDataHandling
+from pystencils import createKernel
+
+
+class LatticeBoltzmannStep:
+
+ def __init__(self, domainSize=None, lbmKernel=None, periodicity=False,
+ kernelParams={}, dataHandling=None, dataPrefix="", optimizationParams={}, **methodParameters):
+
+ # --- Parameter normalization ---
+ if dataHandling is not None:
+ if domainSize is not None:
+ raise ValueError("When passing a dataHandling, the domainSize parameter can not be specified")
+
+ if dataHandling is None:
+ if domainSize is None:
+ raise ValueError("Specify either domainSize or dataHandling")
+ dataHandling = SerialDataHandling(domainSize, defaultGhostLayers=1, periodicity=periodicity)
+
+ methodParameters, optimizationParams = updateWithDefaultParameters(methodParameters, optimizationParams)
+ target = optimizationParams['target']
+
+ # --- Kernel creation ---
+ if lbmKernel:
+ Q = len(lbmKernel.method.stencil)
+ else:
+ Q = len(getStencil(methodParameters['stencil']))
+
+ self._dataHandling = dataHandling
+ self._pdfArrName = dataPrefix + "pdfSrc"
+ self._tmpArrName = dataPrefix + "pdfTmp"
+ self._velocityArrName = dataPrefix + "velocity"
+ self._densityArrName = dataPrefix + "density"
+
+ self._gpu = target == 'gpu'
+ layout = optimizationParams['fieldLayout']
+ self._dataHandling.addArray(self._pdfArrName, fSize=Q, gpu=self._gpu, layout=layout)
+ self._dataHandling.addArray(self._tmpArrName, fSize=Q, gpu=self._gpu, cpu=not self._gpu, layout=layout)
+ self._dataHandling.addArray(self._velocityArrName, fSize=self._dataHandling.dim, gpu=False, layout=layout)
+ self._dataHandling.addArray(self._densityArrName, fSize=1, gpu=False, layout=layout)
+
+ self._kernelParams = kernelParams
+
+ if lbmKernel is None:
+ switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, self._kernelParams)
+ optimizationParams['symbolicField'] = dataHandling.fields[self._pdfArrName]
+ methodParameters['fieldName'] = self._pdfArrName
+ methodParameters['secondFieldName'] = self._tmpArrName
+ self._lbmKernel = createLatticeBoltzmannFunction(optimizationParams=optimizationParams, **methodParameters)
+ else:
+ assert self._dataHandling.dim == self._lbmKernel.method.dim, \
+ "Error: %dD Kernel for %D domain" % (self._lbmKernel.method.dim, self._dataHandling.dim)
+ self._lbmKernel = lbmKernel
+
+ # -- Boundary Handling & Synchronization ---
+ if self._gpu:
+ self._sync = dataHandling.synchronizationFunctionGPU([self._pdfArrName], methodParameters['stencil'])
+ else:
+ self._sync = dataHandling.synchronizationFunctionCPU([self._pdfArrName], methodParameters['stencil'])
+ self._boundaryHandling = BoundaryHandling(self._lbmKernel.method, self._dataHandling, self._pdfArrName,
+ name=dataPrefix + "_boundaryHandling",
+ target=target, openMP=optimizationParams['openMP'])
+
+ # -- Macroscopic Value Kernels
+ self._getterKernel, self._setterKernel = self._compilerMacroscopicSetterAndGetter()
+
+ self.timeStepsRun = 0
+
+ for b in self._dataHandling.iterate():
+ b[self._densityArrName].fill(1.0)
+ b[self._velocityArrName].fill(0.0)
+
+ @property
+ def boundaryHandling(self):
+ """Boundary handling instance of the scenario. Use this to change the boundary setup"""
+ return self._boundaryHandling
+
+ def preRun(self):
+ self._dataHandling.runKernel(self._setterKernel)
+ if self._gpu:
+ self._dataHandling.toGpu(self._pdfArrName)
+
+ def timeStep(self):
+ self._sync()
+ self._boundaryHandling()
+ self._dataHandling.runKernel(self._lbmKernel, **self._kernelParams)
+ self._dataHandling.swap(self._pdfArrName, self._tmpArrName, self._gpu)
+ self.timeStepsRun += 1
+
+ def postRun(self):
+ if self._gpu:
+ self._dataHandling.toCpu(self._pdfArrName)
+ self._dataHandling.runKernel(self._getterKernel)
+
+ def run(self, timeSteps):
+ self.preRun()
+ for i in range(timeSteps):
+ self.timeStep()
+ self.postRun()
+
+ def _compilerMacroscopicSetterAndGetter(self):
+ lbMethod = self._lbmKernel.method
+ D = lbMethod.dim
+ Q = len(lbMethod.stencil)
+ cqc = lbMethod.conservedQuantityComputation
+ pdfField = self._dataHandling.fields[self._pdfArrName]
+ rhoField = self._dataHandling.fields[self._densityArrName]
+ velField = self._dataHandling.fields[self._velocityArrName]
+ pdfSymbols = [pdfField(i) for i in range(Q)]
+
+ getterEqs = cqc.outputEquationsFromPdfs(pdfSymbols, {'density': rhoField, 'velocity': velField})
+ getterKernel = createKernel(getterEqs, target='cpu').compile()
+
+ inpEqs = cqc.equilibriumInputEquationsFromInitValues(rhoField.center, [velField(i) for i in range(D)])
+ setterEqs = lbMethod.getEquilibrium(conservedQuantityEquations=inpEqs)
+ setterEqs = setterEqs.copyWithSubstitutionsApplied({sym: pdfField(i)
+ for i, sym in enumerate(lbMethod.postCollisionPdfSymbols)})
+
+ setterEqs = createSimplificationStrategy(lbMethod)(setterEqs)
+ setterKernel = createKernel(setterEqs, target='cpu').compile()
+ return getterKernel, setterKernel
+
+
+if __name__ == '__main__':
+ from pycuda import autoinit
+ from lbmpy.boundaries import NoSlip, UBB
+ from pystencils import makeSlice
+ step = LatticeBoltzmannStep((30, 30), relaxationRate=1.8, periodicity=True,
+ optimizationParams={'target': 'cpu', 'openMP': False})
+
+ wall = NoSlip()
+ movingWall = UBB((0.001, 0))
+
+ bh = step.boundaryHandling
+ bh.setBoundary(wall, makeSlice[0, :])
+ bh.setBoundary(wall, makeSlice[-1, :])
+ bh.setBoundary(wall, makeSlice[:, 0])
+ bh.setBoundary(movingWall, makeSlice[:, -1])
+ bh.prepare()
+
+ step.run(5000)
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 36c98083..6daad194 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -104,7 +104,10 @@ def freeEnergyFunctionalNPhases(numPhases, surfaceTensions=symmetricSymbolicSurf
else:
phi = orderParameters
- tauFactor = sp.Rational(1, 2) # originally this was 1 / sp.sqrt(2)
+ # Compared to handwritten notes we scale the interface width parameter here to obtain the correct
+ # equations for the interface profile and the surface tensions i.e. to pass tests
+ # test_analyticInterfaceSolution and test_surfaceTensionDerivation
+ interfaceWidth *= sp.sqrt(2)
def f(c):
return c ** 2 * (1 - c) ** 2
@@ -113,12 +116,12 @@ def freeEnergyFunctionalNPhases(numPhases, surfaceTensions=symmetricSymbolicSurf
N = numPhases - 1
if k == l:
assert surfaceTensions(l, l) == 0
- return 6 * tauFactor * interfaceWidth * (surfaceTensions(k, N) + surfaceTensions(l, N) - surfaceTensions(k, l))
+ return 3 / sp.sqrt(2) * interfaceWidth * (surfaceTensions(k, N) + surfaceTensions(l, N) - surfaceTensions(k, l))
def bulkTerm(i, j):
return surfaceTensions(i, j) / 2 * (f(phi[i]) + f(phi[j]) - f(phi[i] + phi[j]))
- F_bulk = 3 * tauFactor / interfaceWidth * sum(bulkTerm(i, j) for i, j in multiSum(2, numPhases) if i != j)
+ F_bulk = 3 / sp.sqrt(2) / interfaceWidth * sum(bulkTerm(i, j) for i, j in multiSum(2, numPhases) if i != j)
F_interface = sum(lambdaCoeff(i, j) / 2 * Diff(phi[i]) * Diff(phi[j]) for i, j in multiSum(2, numPhases))
return F_bulk + F_interface
diff --git a/phasefield/postprocessing.py b/phasefield/postprocessing.py
index a6c63bea..269fee3e 100644
--- a/phasefield/postprocessing.py
+++ b/phasefield/postprocessing.py
@@ -2,10 +2,11 @@ import numpy as np
from matplotlib.path import Path
import itertools
import scipy
+import scipy.spatial
import warnings
-def getIsolines(dataset, level=0.5, refinementFactor=1):
+def getIsolines(dataset, level=0.5, refinementFactor=4):
from matplotlib._contour import QuadContourGenerator
indexArrays = np.meshgrid(*[np.arange(s) for s in dataset.shape])
gen = QuadContourGenerator(*indexArrays, dataset, None, True, 0)
@@ -43,7 +44,7 @@ def findBranchingPoint(pathVertices1, pathVertices2, maxDistance=0.5, branchingP
def findAllBranchingPoints(phaseField1, phaseField2, maxDistance=0.1, branchingPointFilter=3):
result = []
- isoLines = [getIsolines(p, level=0.5, refinementFactor=4) for p in (phaseField1, phaseField2)]
+ isoLines = [getIsolines(p, level=0.5, refinementFactor=16) for p in (phaseField1, phaseField2)]
for path1, path2 in itertools.product(*isoLines):
bbs = findBranchingPoint(path1, path2, maxDistance, branchingPointFilter)
result += list(bbs)
--
GitLab