diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index f2902c8481185baaae12cf46ac31db9d9bbc712c..f67d31cd96895253e02643b635a3c2bd3e21ab8d 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -5,7 +5,7 @@ from pystencils.astnodes import SympyAssignment
 from pystencils.sympyextensions import getSymmetricPart
 from pystencils import Field
 from lbmpy.boundaries.boundary_kernel import offsetFromDir, weightOfDirection, invDir
-from pystencils.data_types import createTypeFromString
+from pystencils.data_types import createType
 
 
 class Boundary(object):
@@ -87,7 +87,7 @@ class NoSlipFullWay(Boundary):
 
     @property
     def additionalData(self):
-        return [('lastValue', createTypeFromString("double"))]
+        return [('lastValue', createType("double"))]
 
     @property
     def additionalDataInitKernelEquations(self):
@@ -133,7 +133,7 @@ class UBB(Boundary):
     @property
     def additionalData(self):
         if callable(self._velocity):
-            return [('vel_%d' % (i,), createTypeFromString("double")) for i in range(self.dim)]
+            return [('vel_%d' % (i,), createType("double")) for i in range(self.dim)]
         else:
             return []
 
diff --git a/creationfunctions.py b/creationfunctions.py
index aacbe746dc765c86738bc10173d00720c3103973..9839eaf58fcfa9fe2ae2ccc64510a5c96fca646f 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -146,9 +146,7 @@ For example, to modify the AST one can run::
 """
 import sympy as sp
 from copy import copy
-from functools import partial
-import json
-from pystencils.cache import diskcache, SympyJSONDecoder, SympyJSONEncoder
+from pystencils.cache import diskcache
 from lbmpy.methods import createSRT, createTRT, createOrthogonalMRT, createKBCTypeTRT, \
     createRawMRT, createThreeRelaxationRateMRT
 from lbmpy.methods.entropic import addIterativeEntropyCondition, addEntropyCondition
@@ -161,6 +159,7 @@ from lbmpy.updatekernels import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAc
     createLBMKernel
 from pystencils.equationcollection.equationcollection import EquationCollection
 from pystencils.field import getLayoutOfArray, Field
+from pystencils import createKernel
 
 
 def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
@@ -276,22 +275,7 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
         params['optimizationParams'] = optParams
         ast = createLatticeBoltzmannAst(**params)
 
-    if optParams['target'] == 'cpu':
-        if optParams['vectorization']:
-            import pystencils.backends.simd_instruction_sets as vec
-            from pystencils.vectorization import vectorize
-            vecParams = optParams['vectorization']
-            vec.selectedInstructionSet = vec.x86VectorInstructionSet(instructionSet=vecParams[0], dataType=vecParams[1])
-            vectorize(ast)
-
-        from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
-        addOpenMP(ast, numThreads=optParams['openMP'])
-        res = makePythonCpuFunction(ast)
-    elif optParams['target'] == 'gpu':
-        from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
-        res = makePythonGpuFunction(ast)
-    else:
-        return ValueError("'target' has to be either 'cpu' or 'gpu'")
+    res = ast.compile()
 
     res.method = ast.method
     res.updateRule = ast.updateRule
@@ -306,36 +290,10 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
         params['optimizationParams'] = optimizationParams
         updateRule = createLatticeBoltzmannUpdateRule(**params)
 
-    if optParams['target'] == 'cpu':
-        from pystencils.cpu import createKernel
-        if 'splitGroups' in updateRule.simplificationHints:
-            splitGroups = updateRule.simplificationHints['splitGroups']
-        else:
-            splitGroups = ()
-        res = createKernel(updateRule.allEquations, splitGroups=splitGroups,
-                           typeForSymbol='double' if optParams['doublePrecision'] else 'float32',
-                           ghostLayers=1)
-    elif optParams['target'] == 'gpu':
-        from pystencils.gpucuda import createCUDAKernel
-        from pystencils.gpucuda.indexing import LineIndexing, BlockIndexing
-        assert optParams['gpuIndexing'] in ('line', 'block')
-        indexingCreator = LineIndexing if optParams['gpuIndexing'] == 'line' else BlockIndexing
-        if optParams['gpuIndexingParams']:
-            indexingCreator = partial(indexingCreator, **optParams['gpuIndexingParams'])
-        res = createCUDAKernel(updateRule.allEquations,
-                               typeForSymbol='double' if optParams['doublePrecision'] else 'float',
-                               indexingCreator=indexingCreator, ghostLayers=1)
-    elif optParams['target'] == 'llvm':
-        from pystencils.llvm import createKernel
-        if 'splitGroups' in updateRule.simplificationHints:
-            splitGroups = updateRule.simplificationHints['splitGroups']
-        else:
-            splitGroups = ()
-        res = createKernel(updateRule.allEquations, splitGroups=splitGroups,
-                           typeForSymbol='double' if optParams['doublePrecision'] else 'float',
-                           ghostLayers=1)
-    else:
-        return ValueError("'target' has to be either 'cpu' or 'gpu'")
+    dtype = 'double' if optParams['doublePrecision'] else 'float32'
+    res = createKernel(updateRule, target=optParams['target'], dataType=dtype,
+                       cpuOpenMP=optParams['openMP'], cpuVectorizeInfo=optParams['vectorization'],
+                       gpuIndexing=optParams['gpuIndexing'], gpuIndexingParams=optParams['gpuIndexingParams'])
 
     res.method = updateRule.method
     res.updateRule = updateRule
diff --git a/forcemodels.py b/forcemodels.py
index dbc93edf4762d60f3be2ddefd2eb936e8c2b9ebe..2f903e4ad8f5ff0628288b4dba5a0844b28cfaf8 100644
--- a/forcemodels.py
+++ b/forcemodels.py
@@ -89,6 +89,14 @@ import sympy as sp
 from lbmpy.methods.relaxationrates import getShearRelaxationRate
 
 
+class ScalarSource(object):
+    def __init__(self, source):
+        self._source = source
+
+    def __call__(self, lbMethod, **kwargs):
+        return [w_i * self._source for w_i in lbMethod.weights]
+
+
 class Simple(object):
     r"""
     A simple force model which introduces the following additional force term in the