diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 213fc1f61b2c3747d59cf929763ee100e72e3c1d..4cfa82d51fee33a3b4dc8182ec07c23c06c82a53 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -11,7 +11,7 @@ INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
 WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
 
 
-class BoundaryHandling:
+class BoundaryHandling(object):
     def __init__(self, symbolicPdfField, domainShape, lbMethod, ghostLayers=1, target='cpu'):
         self._symbolicPdfField = symbolicPdfField
         self._shapeWithGhostLayers = [d + 2 * ghostLayers for d in domainShape]
diff --git a/diskcache.py b/cache.py
similarity index 65%
rename from diskcache.py
rename to cache.py
index 75e7f8f3b7ca20990a7fb84eb25a771650c43baa..96e46531a872bfdd6bcf074593a318adece85c11 100644
--- a/diskcache.py
+++ b/cache.py
@@ -1,10 +1,14 @@
+try:
+    from functools import lru_cache as memorycache
+except ImportError:
+    from backports.functools_lru_cache import lru_cache as memorycache
+
 try:
     from joblib import Memory
     diskcache = Memory(cachedir="/tmp/lbmpy", verbose=False).cache
 except ImportError:
     # fallback to in-memory caching if joblib is not available
-    import functools
-    diskcache = functools.lru_cache(maxsize=64)
+    diskcache = memorycache(maxsize=64)
 
 
 # joblibs Memory decorator does not play nicely with sphinx autodoc (decorated functions do not occur in autodoc)
@@ -12,5 +16,4 @@ except ImportError:
 import sys
 calledBySphinx = 'sphinx' in sys.modules
 if calledBySphinx:
-    import functools
-    diskcache = functools.lru_cache(maxsize=64)
\ No newline at end of file
+    diskcache = memorycache(maxsize=64)
\ No newline at end of file
diff --git a/continuous_distribution_measures.py b/continuous_distribution_measures.py
index df38428282dbb68b8ee4c746b7ec305c5b40e8ee..acf6cf05094db5a1914dd9f60e1cfd8604e91231 100644
--- a/continuous_distribution_measures.py
+++ b/continuous_distribution_measures.py
@@ -3,12 +3,11 @@
 """
 
 import sympy as sp
-import functools
 from pystencils.sympyextensions import makeExponentialFuncArgumentSquares
-from lbmpy.diskcache import diskcache
+from lbmpy.cache import diskcache, memorycache
 
 
-@functools.lru_cache()
+@memorycache()
 def momentGeneratingFunction(function, symbols, symbolsInResult):
     """
     Computes the moment generating function of a probability distribution. It is defined as:
@@ -87,7 +86,7 @@ def multiDifferentiation(generatingFunction, index, symbols):
     return r
 
 
-@functools.lru_cache(maxsize=512)
+@memorycache(maxsize=512)
 def __continuousMomentOrCumulant(function, moment, symbols, generatingFunction):
 
     if type(moment) is tuple:
diff --git a/creationfunctions.py b/creationfunctions.py
index 6c0626b6643caf243f8b53040a791badba8b7b45..e4a2a7a5bf4b68fd9344e4c9b98dcdc8dcd06459 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -54,7 +54,8 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
     params, optParams = _getParams(kwargs, optimizationParams)
 
     if ast is None:
-        ast = createLatticeBoltzmannAst(**params, optimizationParams=optParams)
+        params['optimizationParams'] = optParams
+        ast = createLatticeBoltzmannAst(**params)
 
     if params['target'] == 'cpu':
         from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
@@ -79,7 +80,8 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
     params, optParams = _getParams(kwargs, optimizationParams)
 
     if updateRule is None:
-        updateRule = createLatticeBoltzmannUpdateRule(**params, optimizationParams=optimizationParams)
+        params['optimizationParams'] = optimizationParams
+        updateRule = createLatticeBoltzmannUpdateRule(**params)
 
     if params['target'] == 'cpu':
         from pystencils.cpu import createKernel
@@ -168,12 +170,11 @@ def createLatticeBoltzmannMethod(**params):
         assert len(relaxationRates) >= 2, "Not enough relaxation rates"
         method = createTRT(stencil, relaxationRates[0], relaxationRates[1], **commonParams)
     elif methodName.lower() == 'mrt':
-        nextRelaxationRate = 0
+        nextRelaxationRate = [0]
 
         def relaxationRateGetter(momentGroup):
-            nonlocal nextRelaxationRate
-            res = relaxationRates[nextRelaxationRate]
-            nextRelaxationRate += 1
+            res = relaxationRates[nextRelaxationRate[0]]
+            nextRelaxationRate[0] += 1
             return res
         method = createOrthogonalMRT(stencil, relaxationRateGetter, **commonParams)
     else:
diff --git a/cumulants.py b/cumulants.py
index 1cec741ed8dcaace6bc2f0e7d191b408b9dcd9ba..fcafd342e5482893df01e99de20c6e19450d1826 100644
--- a/cumulants.py
+++ b/cumulants.py
@@ -3,12 +3,10 @@ This module provides functions to compute cumulants of discrete functions.
 Additionally functions are provided to compute cumulants from moments and vice versa
 """
 
-import functools
 import sympy as sp
-
 from lbmpy.moments import momentsUpToComponentOrder
 from lbmpy.continuous_distribution_measures import multiDifferentiation
-
+from lbmpy.cache import memorycache
 
 # ------------------------------------------- Internal Functions -------------------------------------------------------
 from pystencils.sympyextensions import fastSubs
@@ -60,7 +58,6 @@ def __cumulantRawMomentTransform(index, dependentVarDict, outerFunction, default
     subsDict = {}
 
     def createMomentSymbol(idx):
-        nonlocal subsDict
         idx = tuple(idx)
         resultSymbol = sp.Symbol(defaultPrefix + "_" + "_".join(["%d"] * len(idx)) % idx)
         if dependentVarDict is not None and idx in dependentVarDict:
@@ -103,7 +100,7 @@ def __cumulantRawMomentTransform(index, dependentVarDict, outerFunction, default
     return fastSubs(result, subsDict)
 
 
-@functools.lru_cache(maxsize=16)
+@memorycache(maxsize=16)
 def __getDiscreteCumulantGeneratingFunction(function, stencil, waveNumbers):
     assert len(stencil) == len(function)
 
@@ -117,7 +114,7 @@ def __getDiscreteCumulantGeneratingFunction(function, stencil, waveNumbers):
 # ------------------------------------------- Public Functions ---------------------------------------------------------
 
 
-@functools.lru_cache(maxsize=64)
+@memorycache(maxsize=64)
 def discreteCumulant(function, cumulant, stencil):
     """
     Computes cumulant of discrete function
@@ -146,7 +143,7 @@ def discreteCumulant(function, cumulant, stencil):
         return result
 
 
-@functools.lru_cache(maxsize=8)
+@memorycache(maxsize=8)
 def cumulantsFromPdfs(stencil, cumulantIndices=None, pdfSymbols=None):
     """
     Transformation of pdfs (or other discrete function on a stencil) to cumulant space
diff --git a/fieldaccess.py b/fieldaccess.py
index a444d577af77a935fdc94b5b94c693a5999b6ef4..6b04c389d316a9ac077b9404fcbc81217a16945d 100644
--- a/fieldaccess.py
+++ b/fieldaccess.py
@@ -6,7 +6,7 @@ import abc
 # ------------------------------------------------ Interface -----------------------------------------------------------
 
 
-class PdfFieldAccessor:
+class PdfFieldAccessor(object):
     """
     Defines how data is read and written in an LBM time step.
 
diff --git a/forcemodels.py b/forcemodels.py
index ca0c42038b7a9b50755e956e9769a6a8aacf3e09..8d0cf3d717d510a93462a8904b50431dcb4747df 100644
--- a/forcemodels.py
+++ b/forcemodels.py
@@ -6,7 +6,7 @@
 import sympy as sp
 
 
-class Simple:
+class Simple(object):
     r"""
     A simple force model which introduces the following additional force term in the
     collision process :math:`3  w_i  e_i  f_i` (often: force = rho * acceleration)
@@ -26,7 +26,7 @@ class Simple:
                 for direction, w_i in zip(lbMethod.stencil, lbMethod.weights)]
 
 
-class Luo:
+class Luo(object):
     r"""
     Force model by Luo with the following forcing term
 
@@ -53,7 +53,7 @@ class Luo:
         return defaultVelocityShift(density, self._force)
 
 
-class Guo:
+class Guo(object):
     r"""
      Force model by Guo with the following term:
 
diff --git a/gridvisualization.py b/gridvisualization.py
index ffca51cdedd0b4d3259917b7dd1095d381c555c9..5cf887db35983754c31e59c6c045d0f89962b584 100644
--- a/gridvisualization.py
+++ b/gridvisualization.py
@@ -1,7 +1,7 @@
 import matplotlib.patches as patches
 
 
-class Grid:
+class Grid(object):
     """Visualizes a 2D LBM grid with matplotlib by drawing cells and pdf arrows"""
 
     def __init__(self, xCells, yCells):
diff --git a/maxwellian_equilibrium.py b/maxwellian_equilibrium.py
index 3fc185cfe99eb35aa7308582e2ab482873b216d8..4b03b82015a7ca779024526a49e97eb12f44f900 100644
--- a/maxwellian_equilibrium.py
+++ b/maxwellian_equilibrium.py
@@ -6,7 +6,7 @@ Additionally functions are provided to compute moments and cumulants of these di
 
 import sympy as sp
 from sympy import Rational as R
-from lbmpy.diskcache import diskcache
+from lbmpy.cache import diskcache
 import warnings
 
 
diff --git a/methods/abstractlbmethod.py b/methods/abstractlbmethod.py
index de64dac564ccbaffc8a99e82d13e1954c9df3122..ed32c2c36cd2b9cbf6f8f39d99da194fff2f544c 100644
--- a/methods/abstractlbmethod.py
+++ b/methods/abstractlbmethod.py
@@ -12,8 +12,7 @@ class LbmCollisionRule(EquationCollection):
         super(LbmCollisionRule, self).__init__(*args, **kwargs)
         self.method = lbmMethod
 
-
-class AbstractLbMethod(metaclass=abc.ABCMeta):
+class AbstractLbMethod(abc.ABCMeta('ABC', (object,), {})):
     """
     Abstract base class for all LBM methods
     """
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 41f31460a80b8caeaa77524de567aaf229c23e4a..de59c5dcb0e522b3f2ac6789b6366fa816348dba 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -5,7 +5,7 @@ import sympy as sp
 from pystencils.equationcollection import EquationCollection
 
 
-class AbstractConservedQuantityComputation(metaclass=abc.ABCMeta):
+class AbstractConservedQuantityComputation(abc.ABCMeta('ABC', (object,), {})):
     """
 
     This class defines how conserved quantities are computed as functions of the pdfs.
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 420d6f324dc7816de2af6e3c0b430c78160f92e0..55eebd626f5f1f992879fb2ced2620e854927a20 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -33,7 +33,7 @@ def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, c
     :param forceModel: force model instance, or None if no external forces
     :param equilibriumAccuracyOrder: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
     :param cumulant: if True relax cumulants instead of moments
-    :return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
+    :return: :class:`lbmpy.methods.MomentBasedLbMethod` instance
     """
     momToRrDict = OrderedDict(momentToRelaxationRateDict)
     assert len(momToRrDict) == len(stencil), \
@@ -107,7 +107,7 @@ def createSRT(stencil, relaxationRate, useContinuousMaxwellianEquilibrium=False,
                            usually called :math:`\omega` in LBM literature
     :param useContinuousMaxwellianEquilibrium: determines if the discrete or continuous maxwellian equilibrium is
                            used to compute the equilibrium moments
-    :return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
+    :return: :class:`lbmpy.methods.MomentBasedLbMethod` instance
     """
     moments = getDefaultMomentSetForStencil(stencil)
     rrDict = OrderedDict([(m, relaxationRate) for m in moments])
@@ -371,10 +371,9 @@ def sortMomentsIntoGroupsOfSameOrder(moments):
 
 
 def defaultRelaxationRateNames():
-    nextIndex = 0
+    nextIndex = [0]
 
     def result(momentList):
-        nonlocal nextIndex
         shearMomentInside = False
         allConservedMoments = True
         for m in momentList:
@@ -388,8 +387,8 @@ def defaultRelaxationRateNames():
         elif allConservedMoments:
             return 0
         else:
-            nextIndex += 1
-            return sp.Symbol("omega_%d" % (nextIndex,))
+            nextIndex[0] += 1
+            return sp.Symbol("omega_%d" % (nextIndex[0],))
 
     return result
 
diff --git a/methods/entropic.py b/methods/entropic.py
index 490d0194bcdeee228a75276a614b2e058df8e17c..4adc69cdb6737986eb7f62c4e766387a252e8be1 100644
--- a/methods/entropic.py
+++ b/methods/entropic.py
@@ -86,7 +86,7 @@ def decompositionByRelaxationRate(updateRule, relaxationRate):
     return affineTerms, linearTerms, quadraticTerms
 
 
-class RelaxationRatePolynomialDecomposition:
+class RelaxationRatePolynomialDecomposition(object):
 
     def __init__(self, collisionRule, freeRelaxationRates, fixedRelaxationRates):
         self._collisionRule = collisionRule
diff --git a/methods/momentbasedsimplifications.py b/methods/momentbasedsimplifications.py
index 6d653d3134651be5c9c02a9a7892144b7f63b826..4f46fafe10c320f4c05609d1fb0ecc2260824ac3 100644
--- a/methods/momentbasedsimplifications.py
+++ b/methods/momentbasedsimplifications.py
@@ -1,7 +1,7 @@
 """
 This module holds special transformations for simplifying the collision equations of moment-based methods.
 All of these transformations operate on :class:`pystencils.equationcollection.EquationCollection` and need special
-simplification hints, which are set by the MomentBasedLbmMethod.
+simplification hints, which are set by the MomentBasedLbMethod.
 """
 import sympy as sp
 from pystencils.sympyextensions import replaceAdditive, replaceSecondOrderProducts, extractMostCommonFactor
diff --git a/moments.py b/moments.py
index 3a7c5579ad0c2705ce61352a622ef63ca7455731..6c7ce3c9c91cac28c1ec1f105fbba6cdbf0f4df3 100644
--- a/moments.py
+++ b/moments.py
@@ -37,12 +37,12 @@ Functions
 import sympy as sp
 import math
 import itertools
-import functools
 from copy import copy
 from collections import Counter
 
 from lbmpy.continuous_distribution_measures import continuousMoment
-from lbmpy_old.transformations import removeHigherOrderTerms
+from lbmpy.cache import memorycache
+from pystencils.sympyextensions import removeHigherOrderTerms
 
 MOMENT_SYMBOLS = sp.symbols("x y z")
 
@@ -84,7 +84,8 @@ def __generateFixedSumTuples(tupleLength, tupleSum, allowedValues=None, ordered=
                 newSum = totalSum - i
                 if newSum < 0:
                     continue
-                yield from recursive_helper(currentList, newPosition, newSum)
+                for item in recursive_helper(currentList, newPosition, newSum):
+                    yield item
         else:
             if totalSum in allowedValues:
                 currentList[-1] = totalSum
@@ -132,7 +133,10 @@ def momentPermutations(exponentTuple):
 
 def momentsOfOrder(order, dim=3, includePermutations=True):
     """All tuples of length 'dim' which sum equals 'order'"""
-    yield from __generateFixedSumTuples(dim, order, ordered=not includePermutations)
+    for item in __generateFixedSumTuples(dim, order, ordered=not includePermutations):
+        assert(len(item) == dim)
+        assert(sum(item) == order)
+        yield item
 
 
 def momentsUpToOrder(order, dim=3, includePermutations=True):
@@ -264,7 +268,7 @@ def isShearMoment(moment):
 isShearMoment.shearMoments = set([c[0] * c[1] for c in itertools.combinations(MOMENT_SYMBOLS, 2)])
 
 
-@functools.lru_cache(maxsize=512)
+@memorycache(maxsize=512)
 def discreteMoment(function, moment, stencil):
     """
     Computes discrete moment of given distribution function
diff --git a/serialscenario.py b/serialscenario.py
index 4224507ca14b56804bf65ae42c048210979e1b88..a897c366205ec8bcf43ceaa40f036de35f332a4a 100644
--- a/serialscenario.py
+++ b/serialscenario.py
@@ -18,19 +18,20 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
         methodParameters['stencil'] = 'D2Q9' if D == 2 else 'D3Q27'
 
     Q = len(getStencil(methodParameters['stencil']))
-    pdfSrc = np.zeros(domainSizeWithGhostLayer + (Q,))
-    pdfDst = np.zeros_like(pdfSrc)
+    pdfArrays = [np.zeros(domainSizeWithGhostLayer + (Q,)),
+                 np.zeros(domainSizeWithGhostLayer + (Q,))]
 
     # Create kernel
     if lbmKernel is None:
-        lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
+        methodParameters['optimizationParams'] = optimizationParameters
+        lbmKernel = createLatticeBoltzmannFunction(**methodParameters)
     method = lbmKernel.method
 
     assert D == method.dim, "Domain size and stencil do not match"
 
     # Boundary setup
     if boundarySetupFunction is not None:
-        symPdfField = Field.createFromNumpyArray('pdfs', pdfSrc, indexDimensions=1)
+        symPdfField = Field.createFromNumpyArray('pdfs', pdfArrays[0], indexDimensions=1)
         boundaryHandling = BoundaryHandling(symPdfField, domainSize, lbmKernel.method)
         boundarySetupFunction(boundaryHandling=boundaryHandling, method=method)
         boundaryHandling.prepare()
@@ -38,25 +39,24 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
         boundaryHandling = None
 
     # Macroscopic value input/output
-    densityArr = np.zeros(domainSizeWithGhostLayer)
-    velocityArr = np.zeros(domainSizeWithGhostLayer + (D,))
-    getMacroscopic = compileMacroscopicValuesGetter(method, ['density', 'velocity'], pdfArr=pdfSrc)
-    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0] * D}, pdfArr=pdfSrc)
-    setMacroscopic(pdfs=pdfSrc)
+    densityArr = [np.zeros(domainSizeWithGhostLayer)]
+    velocityArr = [np.zeros(domainSizeWithGhostLayer + (D,))]
+    getMacroscopic = compileMacroscopicValuesGetter(method, ['density', 'velocity'], pdfArr=pdfArrays[0])
+    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0] * D}, pdfArr=pdfArrays[0])
+    setMacroscopic(pdfs=pdfArrays[0])
 
     # Run simulation
     def timeLoop(timeSteps):
-        nonlocal pdfSrc, pdfDst, densityArr, velocityArr
         for t in range(timeSteps):
             for f in preUpdateFunctions:
-                f(pdfSrc)
+                f(pdfArrays[0])
             if boundaryHandling is not None:
-                boundaryHandling(pdfs=pdfSrc)
-            lbmKernel(src=pdfSrc, dst=pdfDst)
+                boundaryHandling(pdfs=pdfArrays[0])
+            lbmKernel(src=pdfArrays[0], dst=pdfArrays[1])
 
-            pdfSrc, pdfDst = pdfDst, pdfSrc
-        getMacroscopic(pdfs=pdfSrc, density=densityArr, velocity=velocityArr)
-        return pdfSrc, densityArr, velocityArr
+            pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0]
+        getMacroscopic(pdfs=pdfArrays[0], density=densityArr[0], velocity=velocityArr[0])
+        return pdfArrays[0], densityArr[0], velocityArr[0]
 
     timeLoop.kernel = lbmKernel
 
@@ -163,13 +163,3 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
     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)