diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index d65ad42477ba257b57cadd0eea0c514f9462bc81..459c87992daa84ae6a48b3a11af357430e2b6965 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -52,7 +52,11 @@ def fixedDensity(pdfField, direction, lbMethod, density):
     neighbor = offsetFromDir(direction, lbMethod.dim)
     inverseDir = invDir(direction)
 
-    symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium())
+    velocity = lbMethod.conservedQuantityComputation.definedSymbols()['velocity']
+    symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
+    substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
+    symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
+
     simplification = createSimplificationStrategy(lbMethod)
     symmetricEq = simplification(symmetricEq)
 
diff --git a/creationfunctions.py b/creationfunctions.py
index 6d2e8121ff52a452c25510d7400e3aa0451beae5..6c0626b6643caf243f8b53040a791badba8b7b45 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -34,6 +34,7 @@ def _getParams(params, optParams):
         'fieldLayout': 'reverseNumpy',  # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
 
         'openMP': True,
+        'pdfArr': None,
     }
     unknownParams = [k for k in params.keys() if k not in defaultMethodDescription]
     unknownOptParams = [k for k in optParams.keys() if k not in defaultOptimizationDescription]
@@ -70,6 +71,7 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
         return ValueError("'target' has to be either 'cpu' or 'gpu'")
 
     res.method = ast.method
+    res.ast = ast
     return res
 
 
@@ -116,11 +118,14 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
         npField = createPdfArray(optParams['fieldSize'], len(stencil), layout=optParams['fieldLayout'])
         updateRule = createStreamPullKernel(collisionRule, numpyField=npField)
     else:
-        layoutName = optParams['fieldLayout']
-        if layoutName == 'fzyx' or 'zyxf':
-            dim = len(stencil[0])
-            layoutName = tuple(reversed(range(dim)))
-        updateRule = createStreamPullKernel(collisionRule, genericLayout=layoutName)
+        if 'pdfArr' in optParams:
+            updateRule = createStreamPullKernel(collisionRule, numpyField=optParams['pdfArr'])
+        else:
+            layoutName = optParams['fieldLayout']
+            if layoutName == 'fzyx' or 'zyxf':
+                dim = len(stencil[0])
+                layoutName = tuple(reversed(range(dim)))
+            updateRule = createStreamPullKernel(collisionRule, genericLayout=layoutName)
 
     return updateRule
 
diff --git a/macroscopicValueKernels.py b/macroscopicValueKernels.py
index ea03843b9cef4a8746f42509770442693526ed8a..084a6f6af90725598160c851030c1cfb310f97c8 100644
--- a/macroscopicValueKernels.py
+++ b/macroscopicValueKernels.py
@@ -1,6 +1,6 @@
 import sympy as sp
 from copy import deepcopy
-from pystencils.field import Field
+from pystencils.field import Field, getLayoutFromNumpyArray
 from lbmpy.simplificationfactory import createSimplificationStrategy
 
 
@@ -30,6 +30,7 @@ def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fiel
         pdfField = Field.createGeneric('pdfs', lbMethod.dim, indexDimensions=1, layout=fieldLayout)
     else:
         pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
+        fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape)-1])
 
     outputMapping = {}
     for outputQuantity in outputQuantities:
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 12a84f339ee1e1b1570fe488c18d60edfa9c6ff0..aca78f7b1126cc587730406b1da31555af79fdad 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -232,7 +232,7 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, useContinuousMaxwell
 # ----------------------------------------- Comparison view for notebooks ----------------------------------------------
 
 
-def compareMomentBasedLbMethods(reference, other):
+def compareMomentBasedLbMethods(reference, other, showDeviationsOnly=False):
     import ipy_table
     table = []
     captionRows = [len(table)]
@@ -244,10 +244,13 @@ def compareMomentBasedLbMethods(reference, other):
         referenceValue = reference.momentToRelaxationInfoDict[moment].equilibriumValue
         otherValue = other.momentToRelaxationInfoDict[moment].equilibriumValue
         diff = sp.simplify(referenceValue - otherValue)
-        table.append(["$%s$" % (sp.latex(moment), ),
-                      "$%s$" % (sp.latex(referenceValue), ),
-                      "$%s$" % (sp.latex(otherValue), ),
-                      "$%s$" % (sp.latex(diff),)])
+        if showDeviationsOnly and diff == 0:
+            pass
+        else:
+            table.append(["$%s$" % (sp.latex(moment), ),
+                          "$%s$" % (sp.latex(referenceValue), ),
+                          "$%s$" % (sp.latex(otherValue), ),
+                          "$%s$" % (sp.latex(diff),)])
 
     onlyInRef = referenceMoments - otherMoments
     if onlyInRef:
diff --git a/methods/entropic.py b/methods/entropic.py
index 40600ce88792dc1430a4c958fc80a3f22d4c24ee..490d0194bcdeee228a75276a614b2e058df8e17c 100644
--- a/methods/entropic.py
+++ b/methods/entropic.py
@@ -192,7 +192,6 @@ def determineRelaxationRateByEntropyConditionIterative(updateRule, omega_s, omeg
 
 
 def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False, velocityRelaxation=None,
-                                   shearRelaxationRate=sp.Symbol("omega"),
                                    higherOrderRelaxationRate=sp.Symbol("omega_h"),
                                    fixedOmega=None, **kwargs):
     def product(iterable):
@@ -230,7 +229,7 @@ def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False
     else:
         raise ValueError("Unknown model. Supported models KBC-Nx")
 
-    omega_s, omega_h = shearRelaxationRate, higherOrderRelaxationRate
+    omega_s, omega_h = sp.Symbol("omega"), higherOrderRelaxationRate
     shearPart, restPart = decomposition
 
     velRelaxation = omega_s if velocityRelaxation is None else velocityRelaxation
@@ -253,7 +252,7 @@ def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False
         collisionRule = determineRelaxationRateByEntropyCondition(collisionRule, omega_s, omega_h)
 
     if fixedOmega:
-        collisionRule = collisionRule.newWithSubstitutions({omega_s: fixedOmega})
+        collisionRule = collisionRule.copyWithSubstitutionsApplied({omega_s: fixedOmega})
 
     return collisionRule
 
diff --git a/methods/momentbased.py b/methods/momentbased.py
index 1ad8be9d9993bfdc0804ceb3227e4e822b6bc084..ed507c0f48c812792cbf3128dfb1c85670d65a40 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -8,7 +8,7 @@ from pystencils.sympyextensions import replaceAdditive
 
 
 class MomentBasedLbMethod(AbstractLbMethod):
-    def __init__(self, stencil, momentToRelaxationInfoDict, conservedQuantityComputation, forceModel=None):
+    def __init__(self, stencil, momentToRelaxationInfoDict, conservedQuantityComputation=None, forceModel=None):
         """
         Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
         These methods work by transforming the pdfs into moment space using a linear transformation. In the moment
diff --git a/serialscenario.py b/serialscenario.py
index 96fd5fac1b36e74b6868f82e936b98b48f01ae84..83677069a43530ba106c13e2a490c080fc76a3b0 100644
--- a/serialscenario.py
+++ b/serialscenario.py
@@ -4,10 +4,11 @@ from pystencils import Field
 from pystencils.slicing import sliceFromDirection
 from lbmpy.creationfunctions import createLatticeBoltzmannFunction
 from lbmpy.macroscopicValueKernels import compileMacroscopicValuesGetter, compileMacroscopicValuesSetter
-from lbmpy.boundaries import BoundaryHandling, noSlip, ubb
+from lbmpy.boundaries import BoundaryHandling, noSlip, ubb, fixedDensity
 
 
-def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, preUpdateFunctions=[]):
+def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, lbmKernel=None,
+                preUpdateFunctions=[]):
     ghostLayers = 1
     domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
     D = len(domainSize)
@@ -16,7 +17,8 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
         methodParameters['stencil'] = 'D2Q9' if D == 2 else 'D3Q27'
 
     # Create kernel
-    lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
+    if lbmKernel is None:
+        lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
     method = lbmKernel.method
 
     assert D == method.dim, "Domain size and stencil do not match"
@@ -56,10 +58,12 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
         getMacroscopic(pdfs=pdfSrc, density=densityArr, velocity=velocityArr)
         return pdfSrc, densityArr, velocityArr
 
+    timeLoop.kernel = lbmKernel
+
     return timeLoop
 
 
-def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, **kwargs):
+def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, lbmKernel=None, **kwargs):
     def boundarySetupFunction(boundaryHandling, method):
         myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:method.dim])
         myUbb.name = 'ubb'
@@ -67,10 +71,57 @@ def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={},
         for direction in ('W', 'E', 'S') if method.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
             boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
 
-    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters)
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel)
 
 
-def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, optimizationParameters={}, **kwargs):
+def runPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None, lbmKernel=None,
+                                     optimizationParameters={}, **kwargs):
+    assert dim in (2, 3)
+
+    if radius is not None:
+        assert length is not None
+        if dim == 3:
+            domainSize = (length, 2 * radius + 1, 2 * radius + 1)
+            roundChannel = True
+        else:
+            if domainSize is None:
+                domainSize = (length, 2 * radius)
+    else:
+        roundChannel = False
+
+    def boundarySetupFunction(boundaryHandling, method):
+        pressureBoundaryInflow = partial(fixedDensity, density=1.0 + pressureDifference)
+        pressureBoundaryInflow.__name__ = "Inflow"
+
+        pressureBoundaryOutflow = partial(fixedDensity, density=1.0)
+        pressureBoundaryOutflow.__name__ = "Outflow"
+        boundaryHandling.setBoundary(pressureBoundaryInflow, sliceFromDirection('W', method.dim))
+        boundaryHandling.setBoundary(pressureBoundaryOutflow, sliceFromDirection('E', method.dim))
+
+        if method.dim == 2:
+            for direction in ('N', 'S'):
+                boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+        elif method.dim == 3:
+            if roundChannel:
+                noSlipIdx = boundaryHandling.addBoundary(noSlip)
+                ff = boundaryHandling.flagField
+                yMid = ff.shape[1] // 2
+                zMid = ff.shape[2] // 2
+                y, z = np.meshgrid(range(ff.shape[1]), range(ff.shape[2]))
+                ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
+            else:
+                for direction in ('N', 'S', 'T', 'B'):
+                    boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+
+    assert domainSize is not None
+    if 'forceModel' not in kwargs:
+        kwargs['forceModel'] = 'guo'
+
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel)
+
+
+def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, lbmKernel=None,
+                          optimizationParameters={}, **kwargs):
     assert dim in (2, 3)
     kwargs['force'] = tuple([force, 0, 0][:dim])
 
@@ -109,5 +160,16 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
     if 'forceModel' not in kwargs:
         kwargs['forceModel'] = 'guo'
 
-    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters,
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel,
                        preUpdateFunctions=[periodicity])
+
+if __name__ == '__main__':
+    import sympy as sp
+    from pystencils.display_utils import highlightCpp
+    from pystencils.cpu.cpujit import generateC
+    from lbmpy.serialscenario import runPressureGradientDrivenChannel
+    import lbmpy.plot2d as plt
+    timeloop = runPressureGradientDrivenChannel(radius=10, length=30, pressureDifference=0.001,
+                                                relaxationRates=[1.9],
+                                                dim=2)
+    pdfs, rho, vel = timeloop(20)
diff --git a/simplificationfactory.py b/simplificationfactory.py
index ba8962f9ee1cb164d7a8c8c67a0bdab49b111cbd..eb3d7bb3be8d7af6cea552be627166cc1b1b7179 100644
--- a/simplificationfactory.py
+++ b/simplificationfactory.py
@@ -30,11 +30,12 @@ def createSimplificationStrategy(lbmMethod, doCseInOpposingDirections=False, doO
         if splitInnerLoop:
             s.add(createLbmSplitGroups)
 
+    s.add(addSubexpressionsForDivisions)
+
     if doCseInOpposingDirections:
         s.add(cseInOpposingDirections)
     if doOverallCse:
         s.add(sympyCSE)
 
-    s.add(addSubexpressionsForDivisions)
 
     return s
diff --git a/stencils.py b/stencils.py
index a0d3634a77b5f86072b37add894b0b678185d2bc..eca0b1ac5b051b58fa9826f36be233854d5062a5 100644
--- a/stencils.py
+++ b/stencils.py
@@ -114,6 +114,21 @@ def stencilsHaveSameEntries(s1, s2):
 
 # -------------------------------------- Visualization -----------------------------------------------------------------
 
+def visualizeStencil(stencil, **kwargs):
+    dim = len(stencil[0])
+    if dim == 2:
+        visualizeStencil2D(stencil, **kwargs)
+    else:
+        slicing = False
+        if 'slice' in kwargs:
+            slicing = kwargs['slice']
+            del kwargs['slice']
+
+        if slicing:
+            visualizeStencil3DBySlicing(stencil, **kwargs)
+        else:
+            visualizeStencil3D(stencil, **kwargs)
+
 
 def visualizeStencil2D(stencil, axes=None, data=None):
     """