diff --git a/lbstep.py b/lbstep.py
index bcafc392b5453fed1d19507493574045162d046b..1cf28aae7d1d054247821002b2edd39cbb350d80 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -12,8 +12,6 @@ from pystencils.timeloop import TimeLoop
 
 class LatticeBoltzmannStep:
 
-    vtkScenarioNrCounter = 0
-
     def __init__(self, domainSize=None, lbmKernel=None, periodicity=False,
                  kernelParams={}, dataHandling=None, name="lbm", optimizationParams={},
                  velocityDataName=None, densityDataName=None, densityDataIndex=None,
@@ -116,9 +114,7 @@ class LatticeBoltzmannStep:
         self.setPdfFieldsFromMacroscopicValues()
 
         # -- VTK output
-        self.vtkWriter = self.dataHandling.vtkWriter(name + str(LatticeBoltzmannStep.vtkScenarioNrCounter),
-                                                     [self.velocityDataName, self.densityDataName])
-        LatticeBoltzmannStep.vtkScenarioNrCounter += 1
+        self.vtkWriter = self.dataHandling.vtkWriter(name, [self.velocityDataName, self.densityDataName])
         self.timeStepsRun = 0
 
     @property
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 819d6f818be919426275e2d3adbbd90769130fec..cb981e83de4ea66fdff53d86b2e05f82855eb6bc 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -29,7 +29,7 @@ def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidt
     kappaPrime = tuple(interfaceWidth**2 * k for k in kappa)
     C = sp.symbols("C_:3")
 
-    bulkFreeEnergy = sum(k * C_i ** 2 * (1 - C_i) ** 2  / 2 for k, C_i in zip(kappa, C))
+    bulkFreeEnergy = sum(k * C_i ** 2 * (1 - C_i) ** 2 / 2 for k, C_i in zip(kappa, C))
     surfaceFreeEnergy = sum(k * Diff(C_i) ** 2 / 2 for k, C_i in zip(kappaPrime, C))
 
     F = 0
@@ -45,18 +45,18 @@ def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidt
         rho, phi, psi = orderParameters
     else:
         rho, phi, psi = sp.symbols("rho phi psi")
-    rhoDef = C[0] + C[1] + C[2]
-    phiDef = C[0] - C[1]
-    psiDef = C[2]
 
-    concentrationToOrderParamRelation = {rho: rhoDef, phi: phiDef, psi: psiDef}
+    transformationMatrix = sp.Matrix([[1,  1, 1],
+                                      [1, -1, 0],
+                                      [0,  0, 1]])
+    rhoDef, phiDef, psiDef = transformationMatrix * sp.Matrix(C)
     orderParamToConcentrationRelation = sp.solve([rhoDef - rho, phiDef - phi, psiDef - psi], C)
 
     F = F.subs(orderParamToConcentrationRelation)
     if expandDerivatives:
         F = expandUsingLinearity(F, functions=orderParameters)
 
-    return F
+    return F, transformationMatrix
 
 
 def freeEnergyFunctionalNPhasesPenaltyTerm(orderParameters, interfaceWidth=interfaceWidthSymbol, kappa=None,
@@ -237,7 +237,6 @@ def functionalDerivative(functional, v, constants=None):
     diffs = functional.atoms(Diff)
 
     diffV = Diff(v)
-    #assert diffV in diffs  # not necessary in general, but for this use case this should be true
 
     nonDiffPart = functional.subs({d: 0 for d in diffs})
 
diff --git a/phasefield/experiments1D.py b/phasefield/experiments1D.py
index 867e09024034113efec2097e6171131a40a5ea02..fa46de2e53b7a2e4149783ad47fe2cf291195248 100644
--- a/phasefield/experiments1D.py
+++ b/phasefield/experiments1D.py
@@ -36,11 +36,12 @@ def plotStatus(phaseFieldStep, fromX=None, toX=None):
     plt.legend()
 
 
-def plotFreeEnergyBulkContours(freeEnergy, orderParameters, phase0=0, phase1=1, **kwargs):
+def plotFreeEnergyBulkContours(freeEnergy, orderParameters, phase0=0, phase1=1,
+                               xRange=(-0.2, 1.2), yRange=(-0.2, 1.2), **kwargs):
     import lbmpy.plot2d as plt
 
-    x = np.linspace(-.2, 1.2, 100)
-    y = np.linspace(-.2, 1.2, 100)
+    x = np.linspace(xRange[0], xRange[1], 100)
+    y = np.linspace(yRange[0], yRange[1], 100)
     xg, yg = np.meshgrid(x, y)
     substitutions = {op: 0 for i, op in enumerate(orderParameters) if i not in (phase0, phase1)}
     substitutions.update({d: 0 for d in freeEnergy.atoms(Diff)})  # remove interface components of free energy
diff --git a/phasefield/experiments2D.py b/phasefield/experiments2D.py
index dc94cd426fa6057eaf4f529a0d6c73c2c7aeba0d..7b73c9f00cb43266e295a5459163c1f17e4b94e5 100644
--- a/phasefield/experiments2D.py
+++ b/phasefield/experiments2D.py
@@ -1,94 +1,8 @@
-import numpy as np
-from itertools import cycle
-
 from lbmpy.boundaries import NoSlip
+from lbmpy.plot2d import phasePlot
 from pystencils import makeSlice
 
 
-def setPhase(pfStep, phaseIndex, sliceObj, value=1.0, others=0.0):
-    for b in pfStep.dataHandling.iterate(sliceObj):
-        b[pfStep.phiFieldName][..., phaseIndex].fill(value)
-        if others is not None:
-            for i in range(pfStep.numOrderParameters):
-                if i == phaseIndex:
-                    continue
-                else:
-                    b[pfStep.phiFieldName][..., i].fill(others)
-
-
-def boxBetweenPhases(pfStep, phaseIndices=(0, 1, 2), values=((1, 0, 0), (0, 1, 0), (0, 0, 1)),
-                     boxX=(0.25, 0.75), boxY=(0.25, 0.75)):
-    assert all(idx < pfStep.numOrderParameters for idx in phaseIndices)
-    assert pfStep.dataHandling.dim == 2
-
-    domainSize = pfStep.dataHandling.shape
-    x1, x2 = (int(c * domainSize[0]) for c in boxX)
-    y1, y2 = (int(c * domainSize[1]) for c in boxY)
-    
-    pfStep.reset()
-    
-    for b in pfStep.dataHandling.iterate():
-        x, y = b.cellIndexArrays[0:2]
-        phi = b[pfStep.phiFieldName]
-
-        innerBox = np.logical_and(np.logical_and(x < x2, x > x1),
-                                  np.logical_and(y < y2, y > y1))
-
-        upperValues, lowerValues, midValues = values
-        for phaseIdx, val in zip(phaseIndices, lowerValues):
-            phi[y <= domainSize[1] //2, phaseIdx] = val
-        for phaseIdx, val in zip(phaseIndices, upperValues):
-            phi[y > domainSize[1] // 2, phaseIdx] = val
-        for phaseIdx, val in zip(phaseIndices, midValues):
-            phi[innerBox, phaseIdx] = val
-
-    pfStep.setPdfFieldsFromMacroscopicValues()
-
-
-def initThreeBoxes(pfStep, phaseIndices=(0, 1, 2), values=((1, 0, 0), (0, 1, 0), (0, 0, 1))):
-    assert all(idx < pfStep.numOrderParameters for idx in phaseIndices)
-    assert pfStep.dataHandling.dim == 2
-
-    domainSize = pfStep.dataHandling.shape
-    x1, x2 = (domainSize[0] // 3), (domainSize[0] // 3) * 2
-    y1, y2 = (domainSize[1] // 3), (domainSize[1] // 3) * 2
-
-    pfStep.reset()
-
-    for b in pfStep.dataHandling.iterate():
-        x, y = b.cellIndexArrays[0:2]
-        phi = b[pfStep.phiFieldName]
-        xMid = np.logical_and(x1 < x, x < x2)
-        yMid = np.logical_and(y1 < y, y < y2)
-
-        center = np.logical_and(xMid, yMid)
-        side = np.logical_and(np.logical_not(xMid),yMid)
-        topBottom = np.logical_not(yMid)
-
-        midValues, sideValues, topBottomValues = values
-        for phaseIdx, val in zip(phaseIndices, midValues):
-            phi[center, phaseIdx] = val
-        for phaseIdx, val in zip(phaseIndices, sideValues):
-            phi[side, phaseIdx] = val
-        for phaseIdx, val in zip(phaseIndices, topBottomValues):
-            phi[topBottom, phaseIdx] = val
-
-    pfStep.setPdfFieldsFromMacroscopicValues()
-
-
-def phasePlot(pfStep, sliceObj=makeSlice[:, :], linewidth=1.0):
-    import lbmpy.plot2d as plt
-    colors = ['#fe0002', '#00fe00', '#0000ff', '#ffa800', '#f600ff']
-    colorCycle = cycle(colors)
-
-    for i in range(pfStep.numOrderParameters):
-        s = sliceObj + (i,)
-        plt.scalarFieldAlphaValue(pfStep.phi[s], next(colorCycle), clip=True, interpolation='bilinear')
-    for i in range(pfStep.numOrderParameters):
-        s = sliceObj + (i,)
-        plt.scalarFieldContour(pfStep.phi[s], levels=[0.5], colors='k', linewidth=linewidth)
-
-
 def dropsBetweenTwoPhase():
     import lbmpy.plot2d as plt
     import sympy as sp
@@ -97,15 +11,12 @@ def dropsBetweenTwoPhase():
     c = sp.symbols("c_:4")
     F = freeEnergyFunctionalNPhasesPenaltyTerm(c, 1, [0.01, 0.02, 0.001, 0.02], 0.0001)
     sc = PhaseFieldStep(F, c, domainSize=(2*200, 2*70), openMP=4, hydroDynamicRelaxationRate=1.9)
-    setPhase(sc, 0, makeSlice[:, 0.5:])
-    setPhase(sc, 1, makeSlice[:, :0.5])
-    setPhase(sc, 2, makeSlice[0.2:0.4, 0.3:0.7])
-    setPhase(sc, 3, makeSlice[0.7:0.8, 0.3:0.7])
+    sc.setConcentration(makeSlice[:, 0.5:], [1, 0, 0, 0])
+    sc.setConcentration(makeSlice[:, :0.5], [0, 1, 0, 0])
+    sc.setConcentration(makeSlice[0.2:0.4, 0.3:0.7], [0, 0, 1, 0])
+    sc.setConcentration(makeSlice[0.7:0.8, 0.3:0.7], [0, 0, 0, 1])
     sc.setPdfFieldsFromMacroscopicValues()
-    #sc.run(1)
-    #for i in range(3):
-    #    plt.subplot(3, 1, i+1)
-    #    plt.scalarField(sc.phi[:, :, i])
+
     for i in range(500000):
         print("Step", i)
         plt.figure(figsize=(14, 7))
@@ -124,22 +35,18 @@ def fallingDrop():
     from lbmpy.phasefield.phasefieldstep import PhaseFieldStep
     from lbmpy.phasefield.analytical import freeEnergyFunctionalNPhasesPenaltyTerm
     c = sp.symbols("c_:3")
-    F = freeEnergyFunctionalNPhasesPenaltyTerm(c, 1, [0.0008, 0.0008, 0.00002])
-    gravity = -1e-5
+    F = freeEnergyFunctionalNPhasesPenaltyTerm(c, 1, [0.001, 0.001, 0.0005])
+    gravity = -0.5e-5
     sc = PhaseFieldStep(F, c, domainSize=(160, 200), openMP=4,
                         hydroDynamicRelaxationRate=1.9,
-                        orderParameterForce={2: (0, gravity), 1: (0, gravity)})
-    setPhase(sc, 0, makeSlice[:, 0.4:])
-    setPhase(sc, 1, makeSlice[:, :0.4])
-    setPhase(sc, 2, makeSlice[0.45:0.55, 0.8:0.9])
+                        orderParameterForce={2: (0, gravity), 1: (0, 0)})
+    sc.setConcentration(makeSlice[:, 0.4:], [1, 0, 0])
+    sc.setConcentration(makeSlice[:, :0.4], [0, 1, 0])
+    sc.setConcentration(makeSlice[0.45:0.55, 0.8:0.9], [0, 0, 1])
     sc.hydroLbmStep.boundaryHandling.setBoundary(NoSlip(), makeSlice[:, 0])
-    sc.hydroLbmStep.boundaryHandling.setBoundary(NoSlip(), makeSlice[:, -1])
+    #sc.hydroLbmStep.boundaryHandling.setBoundary(NoSlip(), makeSlice[:, -1])
 
     sc.setPdfFieldsFromMacroscopicValues()
-    #sc.run(1)
-    #for i in range(3):
-    #    plt.subplot(3, 1, i+1)
-    #    plt.scalarField(sc.phi[:, :, i])
     for i in range(650):
         print("Step", i)
         plt.figure(figsize=(14, 10))
@@ -157,4 +64,4 @@ def fallingDrop():
 
 
 if __name__ == '__main__':
-    fallingDrop()
\ No newline at end of file
+    fallingDrop()
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index fdee27900d92b98bd896d6019762f4e5f9079427..ed3043be0b628b500c219f1048eb3462b2554934 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -1,4 +1,5 @@
 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, \
@@ -15,13 +16,20 @@ 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,
-                 orderParameterForce=None):
+                 optimizationParams=None, kernelParams={}, dx=1, dt=1, solveCahnHilliardWithFiniteDifferences=False,
+                 orderParameterForce=None, concentrationToOrderParameters=None, orderParametersToConcentrations=None):
+
+        if optimizationParams is None:
+            optimizationParams = {'openMP': False, 'target': 'cpu'}
+        openMP, target = optimizationParams['openMP'], optimizationParams['target']
 
         if dataHandling is None:
             dataHandling = SerialDataHandling(domainSize, periodicity=True)
 
         self.freeEnergy = freeEnergy
+        self.concentrationToOrderParameter = concentrationToOrderParameters
+        self.orderParametersToConcentrations = orderParametersToConcentrations
+
         self.chemicalPotentials = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters)
 
         # ------------------ Adding arrays ---------------------
@@ -67,7 +75,16 @@ class PhaseFieldStep:
 
         # Hydrodynamic LBM
         if densityOrderParameter is not None:
-            hydroLbmParameters['output'] = {'density': self.phiField(orderParameters.index(densityOrderParameter))}
+            densityIdx = orderParameters.index(densityOrderParameter)
+            hydroLbmParameters['computeDensityInEveryStep'] = True
+            hydroLbmParameters['densityDataName'] = self.phiFieldName
+            hydroLbmParameters['densityDataIndex'] = densityIdx
+
+        if 'optimizationParams' not in hydroLbmParameters:
+            hydroLbmParameters['optimizationParams'] = optimizationParams
+        else:
+            hydroLbmParameters['optimizationParams'].update(optimizationParams)
+
         self.hydroLbmStep = LatticeBoltzmannStep(dataHandling=dataHandling, name=name + '_hydroLBM',
                                                  relaxationRate=hydroDynamicRelaxationRate,
                                                  computeVelocityInEveryStep=True, force=self.forceField,
@@ -83,7 +100,8 @@ class PhaseFieldStep:
 
         if solveCahnHilliardWithFiniteDifferences:
             if densityOrderParameter is not None:
-                raise NotImplementedError("densityOrderParameter not supported when CH is solved with finite differences")
+                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)
@@ -97,16 +115,24 @@ class PhaseFieldStep:
                 chStep = LatticeBoltzmannStep(dataHandling=dataHandling, relaxationRate=1, lbMethod=chMethod,
                                               velocityInputArrayName=self.velField.name,
                                               densityDataName=self.phiField.name,
+                                              stencil='D3Q19' if self.dataHandling.dim == 3 else 'D2Q9',
                                               computeDensityInEveryStep=True,
                                               densityDataIndex=i,
-                                              name=name + "_chLbm_%d" % (i,), )
+                                              name=name + "_chLbm_%d" % (i,),
+                                              optimizationParams=optimizationParams)
                 self.cahnHilliardSteps.append(chStep)
 
+        self.vtkWriter = self.dataHandling.vtkWriter(name, [self.phiFieldName, self.muFieldName, self.velFieldName,
+                                                            self.forceFieldName])
+
         self.runHydroLbm = True
         self.densityOrderParameter = densityOrderParameter
         self.timeStepsRun = 0
         self.reset()
 
+    def writeVTK(self):
+        self.vtkWriter(self.timeStepsRun)
+
     def reset(self):
         # Init φ and μ
         self.dataHandling.fill(self.phiFieldName, 0.0)
@@ -125,6 +151,7 @@ class PhaseFieldStep:
 
     def preRun(self):
         if self.gpu:
+            self.dataHandling.toGpu(self.phiFieldName)
             self.dataHandling.toGpu(self.muFieldName)
             self.dataHandling.toGpu(self.forceFieldName)
         self.hydroLbmStep.preRun()
@@ -133,6 +160,7 @@ class PhaseFieldStep:
 
     def postRun(self):
         if self.gpu:
+            self.dataHandling.toCpu(self.phiFieldName)
             self.dataHandling.toCpu(self.muFieldName)
             self.dataHandling.toCpu(self.forceFieldName)
         if self.runHydroLbm:
@@ -154,6 +182,21 @@ class PhaseFieldStep:
 
         self.timeStepsRun += 1
 
+    def setConcentration(self, sliceObj, concentration):
+        if self.concentrationToOrderParameter is not None:
+            phi = self.concentrationToOrderParameter(concentration)
+        else:
+            phi = np.array(concentration)
+
+        for b in self.dataHandling.iterate(sliceObj):
+            for i in range(phi.shape[-1]):
+                b[self.phiFieldName][..., i] = phi[i]
+
+    def setDensity(self, sliceObj, value):
+        for b in self.dataHandling.iterate(sliceObj):
+            for i in range(self.numOrderParameters):
+                b[self.hydroLbmStep.densityDataName].fill(value)
+
     def run(self, timeSteps):
         self.preRun()
         for i in range(timeSteps):
@@ -168,6 +211,10 @@ class PhaseFieldStep:
     def phiSlice(self, sliceObj=None):
         return self._getSlice(self.phiFieldName, sliceObj)
 
+    def concentrationSlice(self, sliceObj=None):
+        phi = self.phiSlice(sliceObj)
+        return phi if self.orderParametersToConcentrations is None else self.orderParametersToConcentrations(phi)
+
     def muSlice(self, sliceObj=None):
         return self._getSlice(self.muFieldName, sliceObj)
 
@@ -181,6 +228,10 @@ class PhaseFieldStep:
     def phi(self):
         return SlicedGetter(self.phiSlice)
 
+    @property
+    def concentration(self):
+        return SlicedGetter(self.concentrationSlice)
+
     @property
     def mu(self):
         return SlicedGetter(self.muSlice)
diff --git a/phasefield/scenarios.py b/phasefield/scenarios.py
index 3e31f43093a85e67f764cc4213420344dcdfaf6c..1f7f4920b92b700d5320cac06cb6270d8e8b0846 100644
--- a/phasefield/scenarios.py
+++ b/phasefield/scenarios.py
@@ -1,6 +1,8 @@
 import sympy as sp
+import numpy as np
 from lbmpy.phasefield.phasefieldstep import PhaseFieldStep
-from lbmpy.phasefield.analytical import freeEnergyFunction3Phases, freeEnergyFunctionalNPhases, symbolicOrderParameters
+from lbmpy.phasefield.analytical import freeEnergyFunction3Phases, freeEnergyFunctionalNPhases, symbolicOrderParameters, \
+    freeEnergyFunctionalNPhasesPenaltyTerm
 
 
 def createThreePhaseModel(alpha=1, kappa=(0.015, 0.015, 0.015), includeRho=True, **kwargs):
@@ -10,15 +12,51 @@ def createThreePhaseModel(alpha=1, kappa=(0.015, 0.015, 0.015), includeRho=True,
                   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)
+        freeEnergy, transformationMatrix = freeEnergyFunction3Phases(orderParameters)
+        freeEnergy = freeEnergy.subs(parameters)
+        M = np.array(transformationMatrix).astype(float)
+        Minv = np.array(transformationMatrix.inv()).astype(float)
+
+        def concentrationToOrderParameters(c):
+            phi = np.tensordot(c, M, axes=([-1], [1]))
+            return phi
+
+        return PhaseFieldStep(freeEnergy, orderParameters, densityOrderParameter=orderParameters[0],
+                              concentrationToOrderParameters=concentrationToOrderParameters,
+                              orderParametersToConcentrations=lambda phi: np.tensordot(phi, Minv, axes=([-1], [1])),
+                              **kwargs)
     else:
         orderParameters = sp.symbols("phi psi")
-        freeEnergy = freeEnergyFunction3Phases((1,) + orderParameters).subs(parameters)
-        return PhaseFieldStep(freeEnergy, orderParameters, densityOrderParameter=None, **kwargs)
+        freeEnergy, transformationMatrix = freeEnergyFunction3Phases((1,) + orderParameters)
+        freeEnergy = freeEnergy.subs(parameters)
+        M = transformationMatrix.copy()
+        M.row_del(0)  # rho is assumed to be 1 - is not required
+        M = np.array(M).astype(float)
+        reverse = transformationMatrix.inv() * sp.Matrix(sp.symbols("rho phi psi"))
+        MinvTrafo = sp.lambdify(sp.symbols("phi psi"), reverse.subs(sp.Symbol("rho"), 1))
+
+        def orderParametersToConcentrations(phi):
+            phi = np.array(phi)
+            transformed = MinvTrafo(phi[..., 0], phi[..., 1])
+            return np.moveaxis(transformed[:, 0], 0, -1)
+
+        def concentrationToOrderParameters(c):
+            phi = np.tensordot(c, M, axes=([-1], [1]))
+            return phi
+
+        return PhaseFieldStep(freeEnergy, orderParameters, densityOrderParameter=None,
+                              concentrationToOrderParameters=concentrationToOrderParameters,
+                              orderParametersToConcentrations=orderParametersToConcentrations,
+                              **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)
+
+
+def createNPhaseModelPenaltyTerm(alpha=1, numPhases=4, kappa=0.015, penaltyTermFactor=0.01, **kwargs):
+    orderParameters = symbolicOrderParameters(numPhases)
+    freeEnergy = freeEnergyFunctionalNPhasesPenaltyTerm(orderParameters, alpha, kappa, penaltyTermFactor)
+    return PhaseFieldStep(freeEnergy, orderParameters, **kwargs)
diff --git a/plot2d.py b/plot2d.py
index 71a927eca1c028f23b7e8c132fcea38c7741c293..05845f98af1800b5fd5cd6a0df2f469695cff528 100644
--- a/plot2d.py
+++ b/plot2d.py
@@ -1,3 +1,6 @@
+from itertools import cycle
+
+from pystencils import makeSlice
 from pystencils.plot2d import *
 
 
@@ -55,3 +58,16 @@ def boundaryHandling(boundaryHandlingObj, indexExpr=None, boundaryNameToColor=No
     plt.axis('equal')
     if showLegend:
         plt.legend(handles=patches, bbox_to_anchor=(1.02, 0.5), loc=2, borderaxespad=0.)
+
+
+def phasePlot(pfStep, sliceObj=makeSlice[:, :], linewidth=1.0):
+    import lbmpy.plot2d as plt
+    colors = ['#fe0002', '#00fe00', '#0000ff', '#ffa800', '#f600ff']
+    colorCycle = cycle(colors)
+
+    concentrations = pfStep.concentration[sliceObj]
+
+    for i in range(concentrations.shape[-1]):
+        plt.scalarFieldAlphaValue(concentrations[..., i], next(colorCycle), clip=True, interpolation='bilinear')
+    for i in range(concentrations.shape[-1]):
+        plt.scalarFieldContour(concentrations[..., i], levels=[0.5], colors='k', linewidths=[linewidth])
\ No newline at end of file