Bugfixes in LBM phase field code

- 3 phase model works now with step sytem

boundaries/boundary_handling.py

@@ -32,7 +32,7 @@ class BoundaryHandling:
                              "If you want to add multiple handlings, choose a different name.")
 
         gpu = self._target == 'gpu'
-        dataHandling.addArray(self._flagFieldName, dtype=np.uint32, cpu=True, gpu=gpu)
+        dataHandling.addArray(self._flagFieldName, dtype=np.uint32, cpu=True, gpu=False)
         dataHandling.addCustomClass(self._indexArrayName, self.IndexFieldBlockData, cpu=True, gpu=gpu)
 
         ffGhostLayers = self._dataHandling.ghostLayersOfField(self._flagFieldName)

lbstep.py

@@ -42,6 +42,7 @@ class LatticeBoltzmannStep:
             Q = len(getStencil(methodParameters['stencil']))
         target = optimizationParams['target']
 
+        self.name = name
         self._dataHandling = dataHandling
         self._pdfArrName = name + "_pdfSrc"
         self._tmpArrName = name + "_pdfTmp"
@@ -67,11 +68,11 @@ class LatticeBoltzmannStep:
         if velocityInputArrayName is not None:
             methodParameters['velocityInput'] = self._dataHandling.fields[velocityInputArrayName]
 
-        self._kernelParams = kernelParams
+        self.kernelParams = kernelParams
 
         # --- Kernel creation ---
         if lbmKernel is None:
-            switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, self._kernelParams)
+            switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, self.kernelParams)
             optimizationParams['symbolicField'] = dataHandling.fields[self._pdfArrName]
             methodParameters['fieldName'] = self._pdfArrName
             methodParameters['secondFieldName'] = self._tmpArrName
@@ -82,10 +83,7 @@ class LatticeBoltzmannStep:
             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._sync = dataHandling.synchronizationFunction([self._pdfArrName], methodParameters['stencil'], target)
         self._boundaryHandling = BoundaryHandling(self._lbmKernel.method, self._dataHandling, self._pdfArrName,
                                                   name=name + "_boundaryHandling",
                                                   target=target, openMP=optimizationParams['openMP'])
@@ -98,6 +96,7 @@ class LatticeBoltzmannStep:
         for b in self._dataHandling.iterate():
             b[self.densityDataName].fill(1.0)
             b[self.velocityDataName].fill(0.0)
+        self.setPdfFieldsFromMacroscopicValues()
 
         # -- VTK output
         self.vtkWriter = self.dataHandling.vtkWriter(name + str(LatticeBoltzmannStep.vtkScenarioNrCounter),
@@ -180,7 +179,6 @@ class LatticeBoltzmannStep:
         return SlicedGetter(self.densitySlice)
 
     def preRun(self):
-        self._dataHandling.runKernel(self._setterKernel, **self._kernelParams)
         if self._gpu:
             self._dataHandling.toGpu(self._pdfArrName)
             if self._dataHandling.isOnGpu(self.velocityDataName):
@@ -188,17 +186,20 @@ class LatticeBoltzmannStep:
             if self._dataHandling.isOnGpu(self.densityDataName):
                 self._dataHandling.toGpu(self.densityDataName)
 
+    def setPdfFieldsFromMacroscopicValues(self):
+        self._dataHandling.runKernel(self._setterKernel, **self.kernelParams)
+
     def timeStep(self):
         self._sync()
-        self._boundaryHandling()
-        self._dataHandling.runKernel(self._lbmKernel, **self._kernelParams)
+        self._boundaryHandling(**self.kernelParams)
+        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, **self._kernelParams)
+        self._dataHandling.runKernel(self._getterKernel, **self.kernelParams)
 
     def run(self, timeSteps):
         self.preRun()

phasefield/analytical.py

@@ -73,17 +73,8 @@ def symmetricSymbolicSurfaceTension(i, j):
     return sp.Symbol("%s_%d_%d" % ((surfaceTensionSymbolName, ) + index))
 
 
-def symbolicOrderParameters(numPhases, symbolicLastOrderParameter=False):
-    """
-    Returns a tuple with numPhases entries, where the all but the last are numbered symbols and the last entry
-    is 1 - others
-    """
-    if symbolicLastOrderParameter:
-        return sp.symbols("%s_:%i" % (orderParameterSymbolName, numPhases))
-    else:
-        phi = sp.symbols("%s_:%i" % (orderParameterSymbolName, numPhases - 1))
-        phi = phi + (1 - sum(phi),)  # choose last order parameter as 1 - sum(others)
-        return phi
+def symbolicOrderParameters(numSymbols):
+    return sp.symbols("%s_:%i" % (orderParameterSymbolName, numSymbols))
 
 
 def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidthSymbol, transformed=True,
@@ -155,10 +146,12 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli
     """
     assert not (numPhases is None and orderParameters is None)
     if orderParameters is None:
-        phi = symbolicOrderParameters(numPhases)
+        phi = symbolicOrderParameters(numPhases-1)
     else:
         phi = orderParameters
-        numPhases = len(phi)
+        numPhases = len(phi) + 1
+
+    phi = tuple(phi) + (1 - sum(phi),)
 
     # 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
@@ -195,10 +188,10 @@ def analyticInterfaceProfile(x, interfaceWidth=interfaceWidthSymbol):
     get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order
     parameter, i.e. at a transition between two phases.
     >>> numPhases = 4
-    >>> x, phi = sp.Symbol("x"), symbolicOrderParameters(numPhases)
-    >>> F = freeEnergyFunctionalNPhases(numPhases)
+    >>> x, phi = sp.Symbol("x"), symbolicOrderParameters(numPhases-1)
+    >>> F = freeEnergyFunctionalNPhases(orderParameters=phi)
     >>> mu = chemicalPotentialsFromFreeEnergy(F)
-    >>> mu0 = mu[0].subs({p: 0 for p in phi[1:-1]})  # mu[0] as function of one order parameter only
+    >>> mu0 = mu[0].subs({p: 0 for p in phi[1:]})  # mu[0] as function of one order parameter only
     >>> solution = analyticInterfaceProfile(x)
     >>> solutionSubstitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) }
     >>> sp.expand(mu0.subs(solutionSubstitution))  # inserting solution should solve the mu_0=0 equation
@@ -220,7 +213,7 @@ def chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters=None):
     return sp.Matrix([functionalDerivative(freeEnergy, op, constants) for op in orderParameters])
 
 
-def createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiField, muField, dx=1):
+def createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiField, muField, dx=1, cse=True):
     """Reads from order parameter (phi) field and updates chemical potentials"""
     chemicalPotential = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters)
     laplaceDiscretization = {Diff(Diff(op)): discreteLaplace(phiField, i, dx)
@@ -229,7 +222,7 @@ def createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiFi
     chemicalPotential = chemicalPotential.subs({op: phiField(i) for i, op in enumerate(orderParameters)})
 
     muSweepEqs = [sp.Eq(muField(i), cp) for i, cp in enumerate(chemicalPotential)]
-    return muSweepEqs #sympyCseOnEquationList(muSweepEqs)
+    return muSweepEqs if not cse else sympyCseOnEquationList(muSweepEqs)
 
 
 def createForceUpdateEquations(forceField, phiField, muField, dx=1):

phasefield/scenario.py

-import numpy as np
-
-from lbmpy.boundaries.periodicityhandling import PeriodicityHandling
-from lbmpy.creationfunctions import updateWithDefaultParameters, createLatticeBoltzmannFunction
-from lbmpy.macroscopic_value_kernels import compileMacroscopicValuesSetter
-from lbmpy.phasefield.cahn_hilliard_lbm import createCahnHilliardLbFunction
-from lbmpy.stencils import getStencil
-from lbmpy.updatekernels import createPdfArray
-from pystencils.field import Field, getLayoutOfArray
-from pystencils.slicing import addGhostLayers, removeGhostLayers
-from lbmpy.phasefield.analytical import symbolicOrderParameters, freeEnergyFunctionalNPhases, \
-    createChemicalPotentialEvolutionEquations, createForceUpdateEquations
-
-
-class PhasefieldScenario(object):
-    def __init__(self, domainSize, numPhases, mobilityRelaxationRates=1.1,
-                 surfaceTensionCallback=lambda i, j: 1e-3 if i !=j else 0, interfaceWidth=3, dx=1, gamma=1,
-                 optimizationParams={}, initialVelocity=None, kernelParams={}, **kwargs):
-
-        self.numPhases = numPhases
-        self.timeStepsRun = 0
-        self.domainSize = domainSize
-
-        # ---- Parameter normalization
-        if not hasattr(mobilityRelaxationRates, '__len__'):
-            mobilityRelaxationRates = [mobilityRelaxationRates] * numPhases
-
-        D = len(domainSize)
-
-        ghostLayers = 1
-        domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
-
-        if 'stencil' not in kwargs:
-            kwargs['stencil'] = 'D2Q9' if D == 2 else 'D3Q27'
-
-        methodParameters, optimizationParams = updateWithDefaultParameters(kwargs, optimizationParams)
-
-        stencil = getStencil(methodParameters['stencil'])
-        fieldLayout = optimizationParams['fieldLayout']
-        Q = len(stencil)
-
-        if isinstance(initialVelocity, np.ndarray):
-            assert initialVelocity.shape[-1] == D
-            initialVelocity = addGhostLayers(initialVelocity, indexDimensions=1, ghostLayers=1,
-                                             layout=getLayoutOfArray(self._pdfArrays[0]))
-        elif initialVelocity is None:
-            initialVelocity = [0] * D
-
-        self.kernelParams = kernelParams
-
-        # ---- Arrays
-        self.velArr = np.zeros(domainSizeWithGhostLayer + (D,), order=fieldLayout)
-        self.muArr = np.zeros(domainSizeWithGhostLayer + (numPhases - 1,), order=fieldLayout)
-        self.phiArr = np.zeros(domainSizeWithGhostLayer + (numPhases - 1,), order=fieldLayout)
-        self.forceArr = np.zeros(domainSizeWithGhostLayer + (D,), order=fieldLayout)
-
-        self._pdfArrays = [[createPdfArray(domainSize, Q, layout=optimizationParams['fieldLayout'])
-                            for i in range(numPhases)],
-                           [createPdfArray(domainSize, Q, layout=optimizationParams['fieldLayout'])
-                            for i in range(numPhases)]]
-
-        # ---- Fields
-        velField = Field.createFromNumpyArray('vel', self.velArr, indexDimensions=1)
-        muField = Field.createFromNumpyArray('mu', self.muArr, indexDimensions=1)
-        phiField = Field.createFromNumpyArray('phi', self.phiArr, indexDimensions=1)
-        forceField = Field.createFromNumpyArray('F', self.forceArr, indexDimensions=1)
-
-        orderParameters = symbolicOrderParameters(numPhases)
-        freeEnergy = freeEnergyFunctionalNPhases(numPhases, surfaceTensionCallback, interfaceWidth, orderParameters)
-
-        # ---- Sweeps
-        muSweepEquations = createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiField, muField, dx)
-        forceSweepEquations = createForceUpdateEquations(numPhases, forceField, phiField, muField, dx)
-        if optimizationParams['target'] == 'cpu':
-            from pystencils.cpu import createKernel, makePythonFunction
-            self.muSweep = makePythonFunction(createKernel(muSweepEquations))
-            self.forceSweep = makePythonFunction(createKernel(forceSweepEquations))
-        else:
-            from pystencils.gpucuda import createCUDAKernel, makePythonFunction
-            self.muSweep = makePythonFunction(createCUDAKernel(muSweepEquations))
-            self.forceSweep = makePythonFunction(createCUDAKernel(forceSweepEquations))
-
-        optimizationParams['pdfArr'] = self._pdfArrays[0][0]
-
-        self.lbSweepHydro = createLatticeBoltzmannFunction(force=[forceField(i) for i in range(D)],
-                                                           output={'velocity': velField},
-                                                           optimizationParams=optimizationParams, **kwargs)
-
-        useFdForCahnHilliard = False
-        if useFdForCahnHilliard:
-            dt = 0.01
-            mobility = 1
-            from lbmpy.phasefield.analytical import cahnHilliardFdKernel
-            self.sweepsCH = [cahnHilliardFdKernel(i, phiField, muField, velField, mobility,
-                                                  dx, dt, optimizationParams['target'])
-                             for i in range(numPhases-1)]
-        else:
-            self.sweepsCH = [createCahnHilliardLbFunction(stencil, mobilityRelaxationRates[i],
-                                                          velField, muField(i), phiField(i), optimizationParams, gamma)
-                             for i in range(numPhases-1)]
-
-        self.lbSweeps = [self.lbSweepHydro] + self.sweepsCH
-
-        self._pdfPeriodicityHandler = PeriodicityHandling(self._pdfArrays[0][0].shape, (True, True, True),
-                                                          optimizationParams['target'])
-
-        assert self.muArr.shape == self.phiArr.shape
-        self._muPhiPeriodicityHandler = PeriodicityHandling(self.muArr.shape, (True, True, True),
-                                                            optimizationParams['target'])
-
-        # Pdf array initialization
-        hydroLbmInit = compileMacroscopicValuesSetter(self.lbSweepHydro.method,
-                                                      {'density': 1.0, 'velocity': initialVelocity},
-                                                      pdfArr=self._pdfArrays[0][0], target='cpu')
-        hydroLbmInit(pdfs=self._pdfArrays[0][0], F=self.forceArr, **self.kernelParams)
-        self.initializeCahnHilliardPdfsAccordingToPhi()
-
-        self._nonPdfArrays = {
-            'phiArr': self.phiArr,
-            'muArr': self.muArr,
-            'velArr': self.velArr,
-            'forceArr': self.forceArr,
-        }
-        self._nonPdfGpuArrays = None
-        self._pdfGpuArrays = None
-        self.target = optimizationParams['target']
-
-        self.hydroVelocitySetter = None
-
-    def updateHydroPdfsAccordingToVelocity(self):
-        if self.hydroVelocitySetter is None:
-            self.hydroVelocitySetter = compileMacroscopicValuesSetter(self.lbSweepHydro.method,
-                                                                      {'density': 1.0, 'velocity': self.velArr},
-                                                                      pdfArr=self._pdfArrays[0][0], target='cpu')
-        self.hydroVelocitySetter(pdfs=self._pdfArrays[0][0], F=self.forceArr, **self.kernelParams)
-
-    def _arraysFromCpuToGpu(self):
-        import pycuda.gpuarray as gpuarray
-        if self._nonPdfGpuArrays is None:
-            self._nonPdfGpuArrays = {name: gpuarray.to_gpu(arr) for name, arr in self._nonPdfArrays.items()}
-            self._pdfGpuArrays = [[gpuarray.to_gpu(arr) for arr in self._pdfArrays[0]],
-                                  [gpuarray.to_gpu(arr) for arr in self._pdfArrays[1]]]
-        else:
-            for name, arr in self._nonPdfArrays.items():
-                self._nonPdfGpuArrays[name].set(arr)
-            for i in range(2):
-                for cpuArr, gpuArr in zip(self._pdfArrays[i], self._pdfGpuArrays[i]):
-                    gpuArr.set(cpuArr)
-
-    def _arraysFromGpuToCpu(self):
-        for name, arr in self._nonPdfArrays.items():
-            self._nonPdfGpuArrays[name].get(arr)
-        for cpuArr, gpuArr in zip(self._pdfArrays[0], self._pdfGpuArrays[0]):
-            gpuArr.get(cpuArr)
-
-    def initializeCahnHilliardPdfsAccordingToPhi(self):
-        for i in range(1, self.numPhases):
-            self._pdfArrays[0][i].fill(0)
-            self._pdfArrays[0][i][..., 0] = self.phiArr[..., i-1]
-
-    def gaussianSmoothPhiFields(self, sigma):
-        from scipy.ndimage.filters import gaussian_filter
-        for i in range(self.phiArr.shape[-1]):
-            gaussian_filter(self.phi[..., i], sigma, output=self.phi[..., i], mode='wrap')
-
-    @property
-    def phi(self):
-        return removeGhostLayers(self.phiArr, indexDimensions=1)
-
-    @property
-    def mu(self):
-        return removeGhostLayers(self.muArr, indexDimensions=1)
-
-    @property
-    def velocity(self):
-        return removeGhostLayers(self.velArr, indexDimensions=1)
-
-    def run(self, timeSteps=1):
-        """Run the scenario for the given amount of time steps"""
-        if self.target == 'gpu':
-            self._arraysFromCpuToGpu()
-            self._timeLoop(self._pdfGpuArrays, timeSteps=timeSteps, **self._nonPdfGpuArrays)
-            self._arraysFromGpuToCpu()
-        else:
-            self._timeLoop(self._pdfArrays, timeSteps=timeSteps, **self._nonPdfArrays)
-        self.timeStepsRun += timeSteps
-
-    def _timeLoop(self, pdfArrays, phiArr, muArr, velArr, forceArr, timeSteps):
-        for t in range(timeSteps):
-            self._muPhiPeriodicityHandler(pdfs=phiArr)
-            self.muSweep(phi=phiArr, mu=muArr)
-
-            self._muPhiPeriodicityHandler(pdfs=muArr)
-            self.forceSweep(mu=muArr, phi=phiArr, F=forceArr)
-
-            for src in pdfArrays[0]:
-                self._pdfPeriodicityHandler(pdfs=src)
-
-            for sweep, src, dst in zip(self.lbSweeps, *pdfArrays):
-                sweep(src=src, dst=dst, F=forceArr, phi=phiArr, vel=velArr, mu=muArr)
-
-            pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0]