diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6a369f6d7658ae2552ef454154daef1530cc209d..1df61b213441351c3bbd1b58149ef8f155300607 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -86,6 +86,8 @@ ubuntu:
       - $ENABLE_NIGHTLY_BUILDS
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/ubuntu
   script:
+    - pip3 install `grep -Eo 'sympy[>=]+[0-9\.]+' setup.py | sed 's/>/=/g'`
+    - pip3 install `grep -Eo 'numpy[>=]+[0-9\.]+' setup.py | sed 's/>/=/g'`
     - mkdir -p ~/.config/matplotlib
     - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
     - pip3 install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
diff --git a/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py b/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
index 72b0069e0a28e2a1bb229a500ef996ab903854be..0c7e59dbb9c75793edff202ef5b32c97afcd17cf 100644
--- a/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
+++ b/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
@@ -29,8 +29,6 @@ from lbmpy.methods.momentbased.moment_transforms import (
     PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT,
     FastCentralMomentTransform)
 
-from lbmpy.methods.centeredcumulant.simplification import insert_aliases, insert_zeros
-
 from lbmpy.methods.centeredcumulant.force_model import CenteredCumulantForceModel
 from lbmpy.methods.centeredcumulant.galilean_correction import (
     contains_corrected_polynomials,
@@ -413,12 +411,5 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
             main_assignments = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
                                 for eq, force_term_symbol in zip(main_assignments, force_term_symbols)]
 
-        #   9) Clean up the subexpression tree
-        ac = AssignmentCollection(main_assignments, subexpressions)
-
-        if pre_simplification and pre_simplification != 'none':
-            ac = insert_aliases(insert_zeros(ac))
-            ac = ac.new_without_unused_subexpressions()
-
         #   Aaaaaand we're done.
-        return LbmCollisionRule(self, ac.main_assignments, ac.subexpressions)
+        return LbmCollisionRule(self, main_assignments, subexpressions)
diff --git a/lbmpy/methods/centeredcumulant/galilean_correction.py b/lbmpy/methods/centeredcumulant/galilean_correction.py
index 413b08f3ed3595a302ab0f3aa0e992f099bfe439..788289ab148610a6d3336118e462a71cb6f38a27 100644
--- a/lbmpy/methods/centeredcumulant/galilean_correction.py
+++ b/lbmpy/methods/centeredcumulant/galilean_correction.py
@@ -11,7 +11,7 @@ corrected_polynomials = [x**2 - y**2, x**2 - z**2, x**2 + y**2 + z**2]
 
 
 def contains_corrected_polynomials(polynomials):
-    return all(polynomials.containts(p) for p in corrected_polynomials)
+    return all(cp in polynomials for cp in corrected_polynomials)
 
 
 def add_galilean_correction(poly_relaxation_eqs, polynomials, correction_terms):
diff --git a/lbmpy/simplificationfactory.py b/lbmpy/simplificationfactory.py
index 7c929fb2fedadfc073485a2a221b9ae5754dba68..cb62817e04b71e75059b47b857ba04437a2b6744 100644
--- a/lbmpy/simplificationfactory.py
+++ b/lbmpy/simplificationfactory.py
@@ -2,6 +2,8 @@ import sympy as sp
 
 from lbmpy.innerloopsplit import create_lbm_split_groups
 from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
+from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
+from lbmpy.methods.centeredcumulant.simplification import insert_aliases, insert_zeros
 from lbmpy.methods.momentbased.momentbasedsimplifications import (
     factor_density_after_factoring_relaxation_times, factor_relaxation_rates,
     replace_common_quadratic_and_constant_term, replace_density_and_velocity, replace_second_order_velocity_products)
@@ -31,4 +33,8 @@ def create_simplification_strategy(lb_method, split_inner_loop=False):
             s.add(subexpression_substitution_in_main_assignments)
             if split_inner_loop:
                 s.add(create_lbm_split_groups)
+    elif isinstance(lb_method, CenteredCumulantBasedLbMethod):
+        s.add(insert_zeros)
+        s.add(insert_aliases)
+        s.add(lambda ac: ac.new_without_unused_subexpressions())
     return s
diff --git a/lbmpy/turbulence_models.py b/lbmpy/turbulence_models.py
index 6b67fa3c9126d4506ffba9d91b5b5f54c08a1eb3..1a9bf911b38ff8d43e6507e4183a70e90e9988bc 100644
--- a/lbmpy/turbulence_models.py
+++ b/lbmpy/turbulence_models.py
@@ -20,6 +20,8 @@ def frobenius_norm(matrix, factor=1):
 def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None):
     method = collision_rule.method
     omega_s = get_shear_relaxation_rate(method)
+    if isinstance(omega_s, float) or isinstance(omega_s, int):
+        raise ValueError("For the smagorinsky model the shear relaxation rate has to be a symbol")
     f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
 
     tau_0 = sp.Symbol("tau_0_")