From 8d622127d2df020fd9d45a203e3bc53abaad5e3e Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Thu, 8 Feb 2018 09:03:41 +0100
Subject: [PATCH] lbmpy phasefield

- step class for LB phasefield generic enough to work with 3-phase and
  N-phase models
- cahn hilliard can either be solved by LBM or by finite differences
- 3 phase model can be solved with rho phase or without
---
 lbstep.py                       |  41 ++++----
 methods/creationfunctions.py    |   2 +-
 phasefield/analytical.py        |  50 +++++++--
 phasefield/cahn_hilliard_lbm.py |  37 ++-----
 phasefield/experiments1D.py     |  56 +++++++++++
 phasefield/phasefieldstep.py    | 173 ++++++++++++++++++++++++++++++++
 phasefield/scenarios.py         |  24 +++++
 7 files changed, 324 insertions(+), 59 deletions(-)
 create mode 100644 phasefield/experiments1D.py
 create mode 100644 phasefield/phasefieldstep.py
 create mode 100644 phasefield/scenarios.py

diff --git a/lbstep.py b/lbstep.py
index 66aac26e..048b67b2 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -16,7 +16,7 @@ class LatticeBoltzmannStep:
 
     def __init__(self, domainSize=None, lbmKernel=None, periodicity=False,
                  kernelParams={}, dataHandling=None, name="lbm", optimizationParams={},
-                 velocityDataName=None, densityDataName=None,
+                 velocityDataName=None, densityDataName=None, densityDataIndex=None,
                  computeVelocityInEveryStep=False, computeDensityInEveryStep=False,
                  velocityInputArrayName=None,
                  **methodParameters):
@@ -48,23 +48,28 @@ class LatticeBoltzmannStep:
         self._tmpArrName = name + "_pdfTmp"
         self.velocityDataName = name + "_velocity" if velocityDataName is None else velocityDataName
         self.densityDataName = name + "_density" if densityDataName is None else densityDataName
+        self.densityDataIndex = densityDataIndex
 
         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._pdfArrName, fSize=Q, gpu=self._gpu, layout=layout, latexName='src')
+        self._dataHandling.addArray(self._tmpArrName, fSize=Q, gpu=self._gpu, cpu=not self._gpu,
+                                    layout=layout, latexName='dst')
 
         if velocityDataName is None:
             self._dataHandling.addArray(self.velocityDataName, fSize=self._dataHandling.dim,
-                                        gpu=self._gpu and computeVelocityInEveryStep, layout=layout)
+                                        gpu=self._gpu and computeVelocityInEveryStep, layout=layout, latexName='u')
         if densityDataName is None:
             self._dataHandling.addArray(self.densityDataName, fSize=1,
-                                        gpu=self._gpu and computeDensityInEveryStep, layout=layout)
+                                        gpu=self._gpu and computeDensityInEveryStep, layout=layout, latexName='ρ')
 
         if computeVelocityInEveryStep:
             methodParameters['output']['velocity'] = self._dataHandling.fields[self.velocityDataName]
         if computeDensityInEveryStep:
-            methodParameters['output']['density'] = self._dataHandling.fields[self.densityDataName]
+            densityField = self._dataHandling.fields[self.densityDataName]
+            if self.densityDataIndex is not None:
+                densityField = densityField(densityDataIndex)
+            methodParameters['output']['density'] = densityField
         if velocityInputArrayName is not None:
             methodParameters['velocityInput'] = self._dataHandling.fields[velocityInputArrayName]
 
@@ -91,17 +96,16 @@ class LatticeBoltzmannStep:
         # -- Macroscopic Value Kernels
         self._getterKernel, self._setterKernel = self._compilerMacroscopicSetterAndGetter()
 
-        self.timeStepsRun = 0
-
-        for b in self._dataHandling.iterate():
-            b[self.densityDataName].fill(1.0)
-            b[self.velocityDataName].fill(0.0)
+        self._dataHandling.fill(self.densityDataName, 1.0, fValue=self.densityDataIndex,
+                                ghostLayers=True, innerGhostLayers=True)
+        self._dataHandling.fill(self.velocityDataName, 0.0, ghostLayers=True, innerGhostLayers=True)
         self.setPdfFieldsFromMacroscopicValues()
 
         # -- VTK output
         self.vtkWriter = self.dataHandling.vtkWriter(name + str(LatticeBoltzmannStep.vtkScenarioNrCounter),
                                                      [self.velocityDataName, self.densityDataName])
         LatticeBoltzmannStep.vtkScenarioNrCounter += 1
+        self.timeStepsRun = 0
 
     @property
     def boundaryHandling(self):
@@ -143,31 +147,25 @@ class LatticeBoltzmannStep:
         if sliceObj is None:
             sliceObj = makeSlice[:, :] if self.dim == 2 else makeSlice[:, :, 0.5]
 
-        indexSlice = None
-        if len(sliceObj) > self.dim:
-            indexSlice = sliceObj[self.dim:]
-            sliceObj = sliceObj[:self.dim]
-            assert len(indexSlice) == 1
-
         result = self._dataHandling.gatherArray(dataName, sliceObj)
         if result is None:
             return
 
         if masked:
-            mask = self.boundaryHandling.getMask(sliceObj, 'fluid', True)
+            mask = self.boundaryHandling.getMask(sliceObj[:self.dim], 'fluid', 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)
 
             result = np.ma.masked_array(result, mask)
-        if indexSlice:
-            result = result[..., indexSlice[-1]]
         return result.squeeze()
 
     def velocitySlice(self, sliceObj=None, masked=True):
         return self._getSlice(self.velocityDataName, sliceObj, masked)
 
     def densitySlice(self, sliceObj=None, masked=True):
+        if self.densityDataIndex is not None:
+            sliceObj += (self.densityDataIndex,)
         return self._getSlice(self.densityDataName, sliceObj, masked)
 
     @property
@@ -231,13 +229,14 @@ class LatticeBoltzmannStep:
         cqc = lbMethod.conservedQuantityComputation
         pdfField = self._dataHandling.fields[self._pdfArrName]
         rhoField = self._dataHandling.fields[self.densityDataName]
+        rhoField = rhoField.center if self.densityDataIndex is None else rhoField(self.densityDataIndex)
         velField = self._dataHandling.fields[self.velocityDataName]
         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)])
+        inpEqs = cqc.equilibriumInputEquationsFromInitValues(rhoField, [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)})
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 845b0a57..0bfe3d0e 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -136,7 +136,7 @@ def createFromEquilibrium(stencil, equilibrium, momentToRelaxationRateDict, comp
     """
     if isinstance(stencil, str):
         stencil = getStencil(stencil)
-    if isinstance(momentToRelaxationRateDict, sp.Symbol) or isinstance(momentToRelaxationRateDict, float):
+    if any(isinstance(momentToRelaxationRateDict, t) for t in (sp.Symbol, float, int)):
         momentToRelaxationRateDict = {m: momentToRelaxationRateDict for m in getDefaultMomentSetForStencil(stencil)}
 
     momToRrDict = OrderedDict(momentToRelaxationRateDict)
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 19faba63..c16028b6 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -201,9 +201,7 @@ def analyticInterfaceProfile(x, interfaceWidth=interfaceWidthSymbol):
 
 
 def chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters=None):
-    """
-    Computes chemical potentials as functional derivative of free energy
-    """
+    """Computes chemical potentials as functional derivative of free energy"""
     syms = freeEnergy.atoms(sp.Symbol)
     if orderParameters is None:
         orderParameters = [s for s in syms if s.name.startswith(orderParameterSymbolName)]
@@ -234,12 +232,52 @@ def createForceUpdateEquations(forceField, phiField, muField, dx=1):
         rhs = 0
         for i in range(muFSize):
             rhs -= phiField(i) * (muField.neighbor(d, 1)(i) - muField.neighbor(d, -1)(i)) / (2 * dx)
+            # In the C code this form is found: when commenting in make sure phi field is synced before!
+            #rhs += muField(i) * (phiField.neighbor(d, 1)(i) - phiField.neighbor(d, -1)(i)) / (2 * dx)
         forceSweepEqs.append(sp.Eq(forceField(d), rhs))
     return forceSweepEqs
 
 
-def cahnHilliardFdKernel(phaseIdx, phi, mu, velocity, mobility, dx, dt):
+def cahnHilliardFdEq(phaseIdx, phi, mu, velocity, mobility, dx, dt):
     from pystencils.finitedifferences import transient, advection, diffusion, Discretization2ndOrder
     cahnHilliard = transient(phi, phaseIdx) + advection(phi, velocity, phaseIdx) - diffusion(mu, mobility, phaseIdx)
-    discretizedEq = Discretization2ndOrder(dx, dt)(cahnHilliard)
-    return [sp.Eq(phi.newFieldWithDifferentName('dst')(phaseIdx), discretizedEq)]
+    return Discretization2ndOrder(dx, dt)(cahnHilliard)
+
+
+class CahnHilliardFDStep:
+    def __init__(self, dataHandling, phiFieldName, muFieldName, velocityFieldName, name='ch_fd', target='cpu',
+                 dx=1, dt=1, mobilities=1):
+        from pystencils import createKernel
+        self.dataHandling = dataHandling
+
+        muField = self.dataHandling.fields[muFieldName]
+        velField = self.dataHandling.fields[velocityFieldName]
+        self.phiField = self.dataHandling.fields[phiFieldName]
+        self.tmpField = self.dataHandling.addArrayLike(name + '_tmp', phiFieldName, latexName='tmp')
+
+        numPhases = self.dataHandling.fSize(phiFieldName)
+        if not hasattr(mobilities, '__len__'):
+            mobilities = [mobilities] * numPhases
+
+        updateEqs = []
+        for i in range(numPhases):
+            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.sync = self.dataHandling.synchronizationFunction([phiFieldName, velocityFieldName, muFieldName],
+                                                              target=target)
+
+    def timeStep(self):
+        self.sync()
+        self.dataHandling.runKernel(self.kernel)
+        self.dataHandling.swap(self.phiField.name, self.tmpField.name)
+
+    def setPdfFieldsFromMacroscopicValues(self):
+        pass
+
+    def preRun(self):
+        pass
+
+    def postRun(self):
+        pass
\ No newline at end of file
diff --git a/phasefield/cahn_hilliard_lbm.py b/phasefield/cahn_hilliard_lbm.py
index 44123ce3..16e05e3e 100644
--- a/phasefield/cahn_hilliard_lbm.py
+++ b/phasefield/cahn_hilliard_lbm.py
@@ -1,13 +1,11 @@
 import sympy as sp
-from lbmpy.creationfunctions import createLatticeBoltzmannFunction, createLatticeBoltzmannUpdateRule, \
-    createLatticeBoltzmannAst
 from lbmpy.methods.creationfunctions import createFromEquilibrium
+from lbmpy.stencils import getStencil
 from pystencils.sympyextensions import kroneckerDelta, multidimensionalSummation
-from lbmpy.moments import getDefaultMomentSetForStencil
 from lbmpy.maxwellian_equilibrium import getWeights
 
 
-def createCahnHilliardEquilibrium(stencil, mu, gamma=1):
+def cahnHilliardLbmMethod(stencil, mu, relaxationRate=sp.Symbol("omega"), gamma=1):
     """Returns LB equilibrium that solves the Cahn Hilliard equation
 
     ..math ::
@@ -16,8 +14,11 @@ def createCahnHilliardEquilibrium(stencil, mu, gamma=1):
     
     :param stencil: tuple of discrete directions
     :param mu: symbolic expression (field access) for the chemical potential
+    :param relaxationRate: relaxation rate of method
     :param gamma: tunable parameter affecting the second order equilibrium moment
     """
+    if isinstance(stencil, str):
+        stencil = getStencil(stencil)
     weights = getWeights(stencil, c_s_sq=sp.Rational(1, 3))
 
     kd = kroneckerDelta
@@ -39,31 +40,5 @@ def createCahnHilliardEquilibrium(stencil, mu, gamma=1):
 
     rho = sp.Symbol("rho")
     equilibrium[0] = rho - sp.expand(sum(equilibrium[1:]))
-    return tuple(equilibrium)
-
-
-def createCahnHilliardLbFunction(stencil, relaxationRate, velocityField, mu, orderParameterOut,
-                                 optimizationParams={}, gamma=1, srcFieldName='src', dstFieldName='dst'):
-    """
-    Update rule for a LB scheme that solves Cahn-Hilliard.
-
-    :param stencil:
-    :param relaxationRate: relaxation rate controls the mobility
-    :param velocityField: velocity field (output from N-S LBM)
-    :param mu: chemical potential field
-    :param orderParameterOut: field where order parameter :math:`\phi` is written to
-    :param optimizationParams: generic optimization parameters passed to creation functions
-    :param gamma: tunable equilibrium parameter
-    """
-    equilibrium = createCahnHilliardEquilibrium(stencil, mu, gamma)
-    rrRates = {m: relaxationRate for m in getDefaultMomentSetForStencil(stencil)}
-    method = createFromEquilibrium(stencil, tuple(equilibrium), rrRates, compressible=True)
-
-    updateRule = createLatticeBoltzmannUpdateRule(method, optimizationParams,
-                                                  output={'density': orderParameterOut},
-                                                  velocityInput=velocityField, fieldName=srcFieldName,
-                                                  secondFieldName=dstFieldName)
-
-    ast = createLatticeBoltzmannAst(updateRule=updateRule, optimizationParams=optimizationParams)
-    return createLatticeBoltzmannFunction(ast=ast, optimizationParams=optimizationParams)
+    return createFromEquilibrium(stencil, tuple(equilibrium), relaxationRate, compressible=True)
 
diff --git a/phasefield/experiments1D.py b/phasefield/experiments1D.py
new file mode 100644
index 00000000..ae304d6c
--- /dev/null
+++ b/phasefield/experiments1D.py
@@ -0,0 +1,56 @@
+import numpy as np
+from pystencils import makeSlice
+
+
+def plotStatus(phaseFieldStep, fromX=None, toX=None):
+    import lbmpy.plot2d as plt
+
+    domainSize = phaseFieldStep.dataHandling.shape
+    assert len(domainSize) == 2 and domainSize[1] == 1, "Not a 1D scenario"
+
+    dh = phaseFieldStep.dataHandling
+
+    numPhases = phaseFieldStep.numOrderParameters
+
+    plt.subplot(1, 3, 1)
+    plt.title('φ')
+    phiName = phaseFieldStep.phiFieldName
+    for i in range(numPhases):
+        plt.plot(dh.gatherArray(phiName, makeSlice[fromX:toX, 0, i]), marker='x', label='φ_%d' % (i,))
+    plt.legend()
+
+    plt.subplot(1, 3, 2)
+    plt.title("μ")
+    muName = phaseFieldStep.muFieldName
+    for i in range(numPhases):
+        plt.plot(dh.gatherArray(muName, makeSlice[fromX:toX, 0, i]), marker='x', label='μ_%d' % (i,));
+    plt.legend()
+
+    plt.subplot(1, 3, 3)
+    plt.title("Force and Velocity")
+    plt.plot(dh.gatherArray(phaseFieldStep.forceFieldName, makeSlice[fromX:toX, 0, 0]), label='F', marker='x')
+    plt.plot(dh.gatherArray(phaseFieldStep.velFieldName, makeSlice[fromX:toX, 0, 0]), label='u', marker='v')
+    plt.legend()
+
+
+def initSharpInterface(pfStep, phaseIdx1=1, phaseIdx2=2, x1=None, x2=None):
+    domainSize = pfStep.dataHandling.shape
+    if x1 is None:
+        x1 = domainSize[0] // 4
+    if x2 is None:
+        x2 = 3 * x1
+
+    for b in pfStep.dataHandling.iterate():
+        x = b.cellIndexArrays[0]
+        mid = np.logical_and(x1 < x, x < x2)
+
+        phi = b[pfStep.phiFieldName]
+
+        if phaseIdx1 is not None:
+            phi[..., phaseIdx1].fill(0)
+            phi[mid, phaseIdx1] = 1
+
+        if phaseIdx2 is not None:
+            phi[..., phaseIdx2].fill(1)
+            phi[mid, phaseIdx2] = 0
+    pfStep.setPdfFieldsFromMacroscopicValues()
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
new file mode 100644
index 00000000..3f04f33c
--- /dev/null
+++ b/phasefield/phasefieldstep.py
@@ -0,0 +1,173 @@
+from lbmpy.lbstep import LatticeBoltzmannStep
+from lbmpy.phasefield.cahn_hilliard_lbm import cahnHilliardLbmMethod
+from pystencils import createKernel
+from lbmpy.phasefield import createChemicalPotentialEvolutionEquations, createForceUpdateEquations
+from lbmpy.phasefield.analytical import chemicalPotentialsFromFreeEnergy, CahnHilliardFDStep
+from pystencils.datahandling import SerialDataHandling
+from pystencils.equationcollection.simplifications import sympyCseOnEquationList
+from pystencils.slicing import makeSlice, SlicedGetter
+
+
+class PhaseFieldStep:
+
+    def __init__(self, freeEnergy, orderParameters, domainSize=None, dataHandling=None,
+                 name='pf', hydroLbmParameters={},
+                 hydroDynamicRelaxationRate=1.0, cahnHilliardRelaxationRates=1.0, densityOrderParameter=None,
+                 target='cpu', openMP=False, kernelParams={}, dx=1, dt=1, solveCahnHilliardWithFiniteDifferences=False):
+
+        if dataHandling is None:
+            dataHandling = SerialDataHandling(domainSize)
+
+        self.freeEnergy = freeEnergy
+        self.chemicalPotentials = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters)
+
+        # ------------------ Adding arrays ---------------------
+        gpu = target == 'gpu'
+        self.gpu = gpu
+        self.numOrderParameters = len(orderParameters)
+        self.phiFieldName = name + "_phi"
+        self.muFieldName = name + "_mu"
+        self.velFieldName = name + "_u"
+        self.forceFieldName = name + "_F"
+        self.dataHandling = dataHandling
+        self.phiField = dataHandling.addArray(self.phiFieldName, fSize=len(orderParameters), gpu=gpu, latexName='φ')
+        self.muField = dataHandling.addArray(self.muFieldName, fSize=len(orderParameters), gpu=gpu, latexName="μ")
+        self.velField = dataHandling.addArray(self.velFieldName, fSize=dataHandling.dim, gpu=gpu, latexName="u")
+        self.forceField = dataHandling.addArray(self.forceFieldName, fSize=dataHandling.dim, gpu=gpu, latexName="F")
+
+        # ------------------ 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)})
+
+        # μ Kernel
+        self.muEqs = createChemicalPotentialEvolutionEquations(F, phi, self.phiField, self.muField, dx, cse=False)
+        self.muKernel = createKernel(sympyCseOnEquationList(self.muEqs), target=target, cpuOpenMP=openMP).compile()
+        self.phiSync = dataHandling.synchronizationFunction([self.phiFieldName], target=target)
+
+        # F Kernel
+        self.forceEqs = createForceUpdateEquations(self.forceField, self.phiField, self.muField, dx)
+        self.forceKernel = createKernel(self.forceEqs, target=target, cpuOpenMP=openMP).compile()
+        self.muSync = dataHandling.synchronizationFunction([self.muFieldName], target=target)
+
+        # Hydrodynamic LBM
+        if densityOrderParameter is not None:
+            hydroLbmParameters['output'] = {'density': self.phiField(orderParameters.index(densityOrderParameter))}
+        self.hydroLbmStep = LatticeBoltzmannStep(dataHandling=dataHandling, name=name + '_hydroLBM',
+                                                 relaxationRate=hydroDynamicRelaxationRate,
+                                                 computeVelocityInEveryStep=True, force=self.forceField,
+                                                 velocityDataName=self.velFieldName, kernelParams=kernelParams,
+                                                 **hydroLbmParameters)
+
+        # Cahn-Hilliard LBMs
+        if not hasattr(cahnHilliardRelaxationRates, '__len__'):
+            cahnHilliardRelaxationRates = [cahnHilliardRelaxationRates] * len(orderParameters)
+
+        self.cahnHilliardSteps = []
+
+        if solveCahnHilliardWithFiniteDifferences:
+            if densityOrderParameter is not None:
+                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)
+            self.cahnHilliardSteps.append(chStep)
+        else:
+            for i, op in enumerate(orderParameters):
+                if op == densityOrderParameter:
+                    continue
+
+                chMethod = cahnHilliardLbmMethod(self.hydroLbmStep.method.stencil, self.muField(i),
+                                                 relaxationRate=cahnHilliardRelaxationRates[i])
+                chStep = LatticeBoltzmannStep(dataHandling=dataHandling, relaxationRate=1, lbMethod=chMethod,
+                                              velocityInputArrayName=self.velField.name,
+                                              densityDataName=self.phiField.name,
+                                              computeDensityInEveryStep=True,
+                                              densityDataIndex=i,
+                                              name=name + "_chLbm_%d" % (i,), )
+                self.cahnHilliardSteps.append(chStep)
+
+        # Init φ and μ
+        self.dataHandling.fill(self.phiFieldName, 0.0)
+        self.dataHandling.fill(self.phiFieldName, 1.0 if densityOrderParameter is not None else 0.0, fValue=0)
+        self.dataHandling.fill(self.muFieldName, 0.0)
+        self.dataHandling.fill(self.forceFieldName, 0.0)
+        self.setPdfFieldsFromMacroscopicValues()
+
+        self.timeStepsRun = 0
+
+        self.runHydroLbm = True
+
+    def setPdfFieldsFromMacroscopicValues(self):
+        self.hydroLbmStep.setPdfFieldsFromMacroscopicValues()
+        for chStep in self.cahnHilliardSteps:
+            chStep.setPdfFieldsFromMacroscopicValues()
+
+    def preRun(self):
+        if self.gpu:
+            self.dataHandling.toGpu(self.muFieldName)
+            self.dataHandling.toGpu(self.forceFieldName)
+        self.hydroLbmStep.preRun()
+        for chStep in self.cahnHilliardSteps:
+            chStep.preRun()
+
+    def postRun(self):
+        if self.gpu:
+            self.dataHandling.toCpu(self.muFieldName)
+            self.dataHandling.toCpu(self.forceFieldName)
+        if self.runHydroLbm:
+            self.hydroLbmStep.postRun()
+        for chStep in self.cahnHilliardSteps:
+            chStep.postRun()
+
+    def timeStep(self):
+        self.phiSync()
+        self.dataHandling.runKernel(self.muKernel)
+
+        self.muSync()
+        self.dataHandling.runKernel(self.forceKernel)
+
+        if self.runHydroLbm:
+            self.hydroLbmStep.timeStep()
+
+        for chLbm in self.cahnHilliardSteps:
+            chLbm.timeStep()
+
+        self.timeStepsRun += 1
+
+    def run(self, timeSteps):
+        self.preRun()
+        for i in range(timeSteps):
+            self.timeStep()
+        self.postRun()
+
+    def _getSlice(self, dataName, sliceObj):
+        if sliceObj is None:
+            sliceObj = makeSlice[:, :] if self.dim == 2 else makeSlice[:, :, 0.5]
+        return self.dataHandling.gatherArray(dataName, sliceObj).squeeze()
+
+    def phiSlice(self, sliceObj=None):
+        return self._getSlice(self.phiFieldName, sliceObj)
+
+    def muSlice(self, sliceObj=None):
+        return self._getSlice(self.muFieldName, sliceObj)
+
+    def velocitySlice(self, sliceObj=None):
+        return self._getSlice(self.velFieldName, sliceObj)
+
+    def forceSlice(self, sliceObj=None):
+        return self._getSlice(self.forceFieldName, sliceObj)
+
+    @property
+    def phi(self):
+        return SlicedGetter(self.phiSlice)
+
+    @property
+    def mu(self):
+        return SlicedGetter(self.muSlice)
+
+    @property
+    def velocity(self):
+        return SlicedGetter(self.velocitySlice)
+
+    @property
+    def force(self):
+        return SlicedGetter(self.forceSlice)
\ No newline at end of file
diff --git a/phasefield/scenarios.py b/phasefield/scenarios.py
new file mode 100644
index 00000000..3e31f430
--- /dev/null
+++ b/phasefield/scenarios.py
@@ -0,0 +1,24 @@
+import sympy as sp
+from lbmpy.phasefield.phasefieldstep import PhaseFieldStep
+from lbmpy.phasefield.analytical import freeEnergyFunction3Phases, freeEnergyFunctionalNPhases, symbolicOrderParameters
+
+
+def createThreePhaseModel(alpha=1, kappa=(0.015, 0.015, 0.015), includeRho=True, **kwargs):
+    parameters = {sp.Symbol("alpha"): alpha,
+                  sp.Symbol("kappa_0"): kappa[0],
+                  sp.Symbol("kappa_1"): kappa[1],
+                  sp.Symbol("kappa_2"): kappa[2]}
+    if includeRho:
+        orderParameters = sp.symbols("rho phi psi")
+        freeEnergy = freeEnergyFunction3Phases(orderParameters).subs(parameters)
+        return PhaseFieldStep(freeEnergy, orderParameters, densityOrderParameter=orderParameters[0], **kwargs)
+    else:
+        orderParameters = sp.symbols("phi psi")
+        freeEnergy = freeEnergyFunction3Phases((1,) + orderParameters).subs(parameters)
+        return PhaseFieldStep(freeEnergy, orderParameters, densityOrderParameter=None, **kwargs)
+
+
+def createNPhaseModel(alpha=1, numPhases=4, surfaceTensions=lambda i, j: 0.005 if i != j else 0, **kwargs):
+    orderParameters = symbolicOrderParameters(numPhases-1)
+    freeEnergy = freeEnergyFunctionalNPhases(numPhases, surfaceTensions, alpha, orderParameters)
+    return PhaseFieldStep(freeEnergy, orderParameters, **kwargs)
-- 
GitLab