diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index 4315452897b7a53ba34818a409e2f78bb8e70274..528d5878ccebefc1b6559cbd6d171c7409d1d6e4 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -1,2 +1,4 @@
 from lbmpy.boundaries.boundaryconditions import NoSlip, UBB, FixedDensity, NeumannByCopy
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
+
+__all__ = ['NoSlip', 'UBB', 'FixedDensity', 'NeumannByCopy', 'LatticeBoltzmannBoundaryHandling']
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 218af24697d15e04e4ca28a2802af37190a5d733..7499a4fe8379953c0a5493292f25acd3dc9c3197 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -181,7 +181,8 @@ class FixedDensity(Boundary):
         density_symbol = cqc.defined_symbols()['density']
 
         density = self._density
-        equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density).new_without_subexpressions()
+        equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
+        equilibrium_input = equilibrium_input.new_without_subexpressions()
         density_eq = equilibrium_input.main_assignments[0]
         assert density_eq.lhs == density_symbol
         transformed_density = density_eq.rhs
diff --git a/chapman_enskog/chapman_enskog.py b/chapman_enskog/chapman_enskog.py
index 904121b1be05d32bcf37b989649e183df4d2f9ce..8c42f1d2e9481031c28002b7a79e764b739dd6df 100644
--- a/chapman_enskog/chapman_enskog.py
+++ b/chapman_enskog/chapman_enskog.py
@@ -544,7 +544,7 @@ def chapman_enskog_ansatz(equation, time_derivative_orders=(1, 3), spatial_deriv
         max_expansion_order = max(max_expansion_order, stop_order)
     equation = equation.subs(subs_dict)
     equation = expand_using_linearity(equation, functions=expanded_pdf_symbols).expand().collect(eps)
-    result = {eps_order: equation.coeff(eps ** eps_order) for eps_order in range(1, 2*max_expansion_order)}
+    result = {eps_order: equation.coeff(eps ** eps_order) for eps_order in range(1, 2 * max_expansion_order)}
     result[0] = equation.subs(eps, 0)
     return result
 
diff --git a/chapman_enskog/chapman_enskog_higher_order.py b/chapman_enskog/chapman_enskog_higher_order.py
index a4280eff63ea533d04c4858658f23f213e3e2f0f..2eb54b8ca85eb6fd73dc9f72b56aa59e6b157081 100644
--- a/chapman_enskog/chapman_enskog_higher_order.py
+++ b/chapman_enskog/chapman_enskog_higher_order.py
@@ -1,6 +1,6 @@
 import sympy as sp
+from pystencils.derivative import full_diff_expand
 from lbmpy.chapman_enskog.chapman_enskog import CeMoment, Diff, take_moments
-from lbmpy.chapman_enskog.derivative import full_diff_expand
 from lbmpy.chapman_enskog.chapman_enskog import expanded_symbol
 from lbmpy.moments import moments_up_to_order, moments_of_order
 from collections import OrderedDict, namedtuple
@@ -35,7 +35,7 @@ def get_solvability_conditions(dim, order):
     solvability_conditions = {}
     for name in ["\\Pi", "\\Upsilon"]:
         for moment_tuple in moments_up_to_order(1, dim=dim):
-            for k in range(order+1):
+            for k in range(order + 1):
                 solvability_conditions[CeMoment(name, superscript=k, moment_tuple=moment_tuple)] = 0
     return solvability_conditions
 
@@ -72,14 +72,15 @@ def determine_higher_order_moments(epsilon_hierarchy, relaxation_rates, moment_c
     for eps_order in range(1, order):
         eps_eq = epsilon_hierarchy[eps_order]
 
-        for order in range(order+1):
+        for order in range(order + 1):
             eqs = sp.Matrix([full_expand(tm(eps_eq * m)) for m in poly_moments(order, dim)])
             unknown_moments = [m for m in eqs.atoms(CeMoment)
                                if m.superscript == eps_order and sum(m.moment_tuple) == order]
             if len(unknown_moments) == 0:
                 for eq in eqs:
-                    time_diffs_in_expr = [d for d in eq.atoms(Diff) if
-                                          (d.target == 't' or d.target == sp.Symbol("t")) and d.superscript == eps_order]
+                    t = sp.Symbol("t")
+                    time_diffs_in_expr = [d for d in eq.atoms(Diff)
+                                          if (d.target == 't' or d.target == t) and d.superscript == eps_order]
                     if len(time_diffs_in_expr) == 0:
                         continue
                     assert len(time_diffs_in_expr) == 1, \
diff --git a/chapman_enskog/chapman_enskog_steady_state.py b/chapman_enskog/chapman_enskog_steady_state.py
index 8bd701eea25141ecb50ea896d2cf490a6870bc80..315bb14458e009462777341244530803ee908ba0 100644
--- a/chapman_enskog/chapman_enskog_steady_state.py
+++ b/chapman_enskog/chapman_enskog_steady_state.py
@@ -1,9 +1,9 @@
 import sympy as sp
 import functools
-from lbmpy.chapman_enskog.chapman_enskog import LbMethodEqMoments, CeMoment, take_moments, insert_moments
+from pystencils.derivative import Diff, DiffOperator, expand_using_linearity, normalize_diff_order
 from pystencils.derivative import collect_derivatives, create_nested_diff
 from pystencils.sympyextensions import normalize_product, multidimensional_sum, kronecker_delta
-from lbmpy.chapman_enskog.derivative import Diff, DiffOperator, expand_using_linearity, normalize_diff_order
+from lbmpy.chapman_enskog.chapman_enskog import LbMethodEqMoments, CeMoment, take_moments, insert_moments
 from lbmpy.chapman_enskog.chapman_enskog import expanded_symbol, chapman_enskog_ansatz, remove_higher_order_u
 
 
diff --git a/chapman_enskog/derivative.py b/chapman_enskog/derivative.py
index f37e2ca9305dc794fe0a1b824bb2caca091bb9d0..5ef5aea3d365784932d56ae2288ced3b696da049 100644
--- a/chapman_enskog/derivative.py
+++ b/chapman_enskog/derivative.py
@@ -1,4 +1,5 @@
-from pystencils.derivative import *
+import sympy as sp
+from pystencils.derivative import Diff
 
 
 def chapman_enskog_derivative_expansion(expr, label, eps=sp.Symbol("epsilon"), start_order=1, stop_order=4):
diff --git a/creationfunctions.py b/creationfunctions.py
index 9b6758bb8d67b6807b6e15993b5ce49164f3c77c..2d8ab4ad46191ddcdd95ff14fa2d8404fb4a919f 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -168,8 +168,8 @@ from lbmpy.relaxationrates import relaxation_rate_from_magic_number
 from lbmpy.stencils import get_stencil, stencils_have_same_entries
 import lbmpy.forcemodels as forcemodels
 from lbmpy.simplificationfactory import create_simplification_strategy
-from lbmpy.updatekernels import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAccessor, CollideOnlyInplaceAccessor, \
-    create_lbm_kernel, create_stream_pull_with_output_kernel
+from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAccessor, CollideOnlyInplaceAccessor
+from lbmpy.updatekernels import create_lbm_kernel, create_stream_pull_with_output_kernel
 
 
 def create_lb_function(ast=None, optimization={}, **kwargs):
diff --git a/cumulants.py b/cumulants.py
index b330f0eb1f134ce0ad03888509dcacc58d79b094..049d0244d44f673bc76e2942be74b563e75f2955 100644
--- a/cumulants.py
+++ b/cumulants.py
@@ -23,7 +23,7 @@ def __get_indexed_symbols(passed_symbols, prefix, indices):
     if passed_symbols is not None:
         return passed_symbols
     else:
-        format_string = "%s_" + "_".join(["%d"]*dim)
+        format_string = "%s_" + "_".join(["%d"] * dim)
         tuple_indices = []
         for i in indices:
             tuple_indices.append(i if type(i) is tuple else (i,))
@@ -65,7 +65,7 @@ def __cumulant_raw_moment_transform(index, dependent_var_dict, outer_function, d
             subs_dict[result_symbol] = dependent_var_dict[idx]
         return result_symbol
 
-    zeroth_moment = create_moment_symbol((0,)*dim)
+    zeroth_moment = create_moment_symbol((0,) * dim)
 
     def outer_function_derivative(n):
         x = zeroth_moment
diff --git a/fieldaccess.py b/fieldaccess.py
index 3f4309c84225d3e2ba54543425307efb064a2da8..95cbc1d5fe51953b93b2f2056b5ed524d49c6f92 100644
--- a/fieldaccess.py
+++ b/fieldaccess.py
@@ -55,9 +55,10 @@ class StreamPullTwoFieldsAccessor(PdfFieldAccessor):
 class PeriodicTwoFieldsAccessor(PdfFieldAccessor):
     """Access scheme that builds periodicity into the kernel.
 
-    Introduces a condition on every load, such that at the borders the periodic value is loaded. The periodicity is specified as a tuple of booleans, one for
-    each direction. The second parameter `ghost_layers` specifies the number of assumed ghost layers of the field.
-    For the periodic kernel itself no ghost layers are required, however other kernels might need them. 
+    Introduces a condition on every load, such that at the borders the periodic value is loaded. The periodicity is
+    specified as a tuple of booleans, one for each direction. The second parameter `ghost_layers` specifies the number
+    of assumed ghost layers of the field. For the periodic kernel itself no ghost layers are required,
+    however other kernels might need them.
     """
     def __init__(self, periodicity, ghost_layers=0):
         self._periodicity = periodicity
@@ -158,6 +159,3 @@ def visualize_pdf_field_accessor(pdf_field_accessor, figure=None):
 
     ax_left.set_title("Read")
     ax_right.set_title("Write")
-
-
-
diff --git a/forcemodels.py b/forcemodels.py
index 559f6dcd0e34f328ad0b1f6ca40959286b1d42c9..3cdbf5ad5b1858e48886cad0346e0637cf406600 100644
--- a/forcemodels.py
+++ b/forcemodels.py
@@ -6,8 +6,8 @@ r"""
 Get started:
 -----------
 
-This module offers different models to introduce a body force in the lattice Boltzmann scheme. 
-If you don't know which model to choose, use :class:`lbmpy.forcemodels.Guo`. 
+This module offers different models to introduce a body force in the lattice Boltzmann scheme.
+If you don't know which model to choose, use :class:`lbmpy.forcemodels.Guo`.
 For incompressible collision models the :class:`lbmpy.forcemodels.Buick` model can be better.
 
 
@@ -26,10 +26,10 @@ to the acceleration :math:`\pmb{a}` for all force models.
 
 .. math ::
 
-    \sum_q \pmb{c}_q F_q = \pmb{a} 
+    \sum_q \pmb{c}_q F_q = \pmb{a}
 
 
-The second order moment is different for the forcing models - if it is zero the model is suited for 
+The second order moment is different for the forcing models - if it is zero the model is suited for
 incompressible flows. For weakly compressible collision operators a force model with a corrected second order moment
 should be chosen.
 
@@ -52,7 +52,7 @@ Models with nonzero second moment have:
     F_q = \frac{w_q}{c_s^2} c_{qi} \; a_i + \frac{w_q}{c_s^4} (c_{qi} c_{qj} - c_s^2 \delta_{ij} ) u_j \, a_i
 
 
-For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term 
+For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term
 :math:`S_{macro}` that we call "macroscopic velocity shift"
     
     .. math ::
@@ -62,7 +62,7 @@ For all force models the computation of the macroscopic velocity has to be adapt
         S_{macro} = \frac{\Delta t}{2} \sum_q F_q
 
 
-Some models also shift the velocity entering the equilibrium distribution. 
+Some models also shift the velocity entering the equilibrium distribution.
 
 Comparison
 ----------
@@ -70,7 +70,8 @@ Comparison
 Force models can be distinguished by 2 options:
 
 **Option 1**:
-    :math:`C_F = 1` and equilibrium velocity is not shifted, or :math:`C_F=(1 - \frac{\omega}{2})` and equilibrum is shifted.
+    :math:`C_F = 1` and equilibrium velocity is not shifted, or :math:`C_F=(1 - \frac{\omega}{2})`
+    and equilibrum is shifted.
  
 **Option 2**: 
     second velocity moment is zero or :math:`F_i u_j + F_j u_i`
@@ -212,4 +213,3 @@ class EDM(object):
 
 def default_velocity_shift(density, force):
     return [f_i / (2 * density) for f_i in force]
-
diff --git a/geometry.py b/geometry.py
index 4f9508d79d94e4e097e7ac6199f1e7646f6b9c60..87f2fe05557fd47f7e434c5bb18e7ccd21a2b7fc 100644
--- a/geometry.py
+++ b/geometry.py
@@ -141,7 +141,7 @@ def get_pipe_velocity_field(domain_size, u_max, flow_direction=0, diameter=None)
         radius = int(round(diameter / 2))
 
     params = [np.arange(s) + 0.5 for s in domain_size]
-    grid = np.meshgrid(*params, indexing='ij')  # type: Tuple[np.array]
+    grid: Tuple[np.array] = np.meshgrid(*params, indexing='ij')
 
     dist = 0
     for i in range(len(domain_size)):
@@ -202,8 +202,8 @@ def add_black_and_white_image(boundary_handling, image_file, target_slice=None,
     # resize necessary if aspect ratio should be constant
     if zoomed_image.shape != target_size:
         resized_image = np.zeros(target_size, dtype=np.bool)
-        mid = [(ts - s)//2 for ts, s in zip(target_size, zoomed_image.shape)]
-        resized_image[mid[0]:zoomed_image.shape[0]+mid[0], mid[1]:zoomed_image.shape[1]+mid[1]] = zoomed_image
+        mid = [(ts - s) // 2 for ts, s in zip(target_size, zoomed_image.shape)]
+        resized_image[mid[0]:zoomed_image.shape[0] + mid[0], mid[1]:zoomed_image.shape[1] + mid[1]] = zoomed_image
         zoomed_image = resized_image
 
     def callback(*coordinates):
diff --git a/innerloopsplit.py b/innerloopsplit.py
index ea1524cf8d3c15421e82e70fbfd5cb18777e1e02..2a3849f19d9335c04de161b3bd97cac72287c495 100644
--- a/innerloopsplit.py
+++ b/innerloopsplit.py
@@ -46,7 +46,7 @@ def create_lbm_split_groups(cr: LbmCollisionRule):
     dim = len(stencil[0])
 
     for direction, eq in zip(stencil, cr.main_assignments):
-        if direction == tuple([0]*dim):
+        if direction == tuple([0] * dim):
             split_groups[0].append(eq.lhs)
             continue
 
diff --git a/methods/__init__.py b/methods/__init__.py
index b364fe560c98e323cd547f7831deb403d5e16ab5..fa6975c9444c04d91be267274df2539cde1a0b44 100644
--- a/methods/__init__.py
+++ b/methods/__init__.py
@@ -1,6 +1,9 @@
-from lbmpy.methods.abstractlbmethod import AbstractLbMethod
-from lbmpy.methods.momentbased import MomentBasedLbMethod, RelaxationInfo
-from lbmpy.methods.creationfunctions import create_srt, create_trt, create_trt_with_magic_number, create_mrt_orthogonal, \
-    create_with_continuous_maxwellian_eq_moments, create_with_discrete_maxwellian_eq_moments, create_trt_kbc, create_mrt_raw, \
-    create_mrt3
-from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
+from lbmpy.methods.momentbased import RelaxationInfo
+from lbmpy.methods.creationfunctions import create_srt, create_trt, create_trt_with_magic_number, create_trt_kbc, \
+    create_mrt_orthogonal, create_mrt_raw, create_mrt3, \
+    create_with_continuous_maxwellian_eq_moments, create_with_discrete_maxwellian_eq_moments
+
+__all__ = ['RelaxationInfo',
+           'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
+           'create_mrt_orthogonal', 'create_mrt_raw', 'create_mrt3',
+           'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments']
diff --git a/methods/abstractlbmethod.py b/methods/abstractlbmethod.py
index 12687314bdaba78d2f281e73400b8390ad64534a..20e7d177d54a78fa73ce1d368f3d6c003933f5ef 100644
--- a/methods/abstractlbmethod.py
+++ b/methods/abstractlbmethod.py
@@ -58,4 +58,3 @@ class AbstractLbMethod(abc.ABC):
     def get_collision_rule(self):
         """Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
          This collision rule defines the collision operator."""
-
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 020548abcc71bb574c50d7fb2e9e1ad51dcacd9c..f58b69eb9b566f1659ba66c97ccaddae1d3c1597 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -2,7 +2,7 @@ import abc
 import sympy as sp
 from collections import OrderedDict
 from pystencils.assignment_collection import AssignmentCollection
-from pystencils.field import Field, Assignment
+from pystencils import Field, Assignment
 
 
 class AbstractConservedQuantityComputation(abc.ABC):
@@ -63,8 +63,8 @@ class AbstractConservedQuantityComputation(abc.ABC):
         be slightly different that the ones used as input for the equilibrium e.g. due to a force model.
 
         :param pdfs: values for the pdf entries
-        :param output_quantity_names_to_symbols: dict mapping of conserved quantity names (See :func:`conserved_quantities`)
-                                            to symbols or field accesses where they should be written to
+        :param output_quantity_names_to_symbols: dict mapping of conserved quantity names
+                (See :func:`conserved_quantities`) to symbols or field accesses where they should be written to
         """
 
     @abc.abstractmethod
@@ -131,11 +131,12 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
     def equilibrium_input_equations_from_pdfs(self, pdfs):
         dim = len(self._stencil[0])
         eq_coll = get_equations_for_zeroth_and_first_order_moment(self._stencil, pdfs, self._symbolOrder0,
-                                                          self._symbolsOrder1[:dim])
+                                                                  self._symbolsOrder1[:dim])
         if self._compressible:
             eq_coll = divide_first_order_moments_by_rho(eq_coll, dim)
 
-        eq_coll = apply_force_model_shift('equilibrium_velocity_shift', dim, eq_coll, self._forceModel, self._compressible)
+        eq_coll = apply_force_model_shift('equilibrium_velocity_shift', dim, eq_coll,
+                                          self._forceModel, self._compressible)
         return eq_coll
 
     def equilibrium_input_equations_from_init_values(self, density=1, velocity=(0, 0, 0)):
@@ -284,8 +285,8 @@ def divide_first_order_moments_by_rho(assignment_collection, dim):
     """
     old_eqs = assignment_collection.main_assignments
     rho = old_eqs[0].lhs
-    new_first_order_moment_eq = [Assignment(eq.lhs, eq.rhs / rho) for eq in old_eqs[1:dim+1]]
-    new_eqs = [old_eqs[0]] + new_first_order_moment_eq + old_eqs[dim+1:]
+    new_first_order_moment_eq = [Assignment(eq.lhs, eq.rhs / rho) for eq in old_eqs[1:dim + 1]]
+    new_eqs = [old_eqs[0]] + new_first_order_moment_eq + old_eqs[dim + 1:]
     return assignment_collection.copy(new_eqs)
 
 
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 8493a290b92aceb1132398de951bd0a5a29b2eb0..328bd409cbe46293fe2e74edd09067be4960ef83 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -310,26 +310,27 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
     shear_tensor_trace = sum(shear_tensor_diagonal)
     shear_tensor_trace_free_diagonal = [dim * d - shear_tensor_trace for d in shear_tensor_diagonal]
 
-    energy_transport_tensor = list(exponents_to_polynomial_representations([a for a in moments_of_order(3, dim, True)
-                                                                            if 3 not in a]))
+    poly_repr = exponents_to_polynomial_representations
+    energy_transport_tensor = list(poly_repr([a for a in moments_of_order(3, dim, True)
+                                              if 3 not in a]))
 
     explicitly_defined = set(rho + velocity + shear_tensor_off_diagonal +
                              shear_tensor_diagonal + energy_transport_tensor)
-    rest = list(set(exponents_to_polynomial_representations(moments_up_to_component_order(2, dim))) - explicitly_defined)
+    rest = list(set(poly_repr(moments_up_to_component_order(2, dim))) - explicitly_defined)
     assert len(rest) + len(explicitly_defined) == 3**dim
 
-    # naming according to paper Karlin 2015: Entropic multirelaxation lattice Boltzmann models for turbulent flows
+    # naming according to paper Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows
     d = shear_tensor_off_diagonal + shear_tensor_trace_free_diagonal[:-1]
     t = [shear_tensor_trace]
     q = energy_transport_tensor
     if method_name == 'KBC-N1':
-        decomposition = [d, t+q+rest]
+        decomposition = [d, t + q + rest]
     elif method_name == 'KBC-N2':
-        decomposition = [d+t, q+rest]
+        decomposition = [d + t, q + rest]
     elif method_name == 'KBC-N3':
-        decomposition = [d+q, t+rest]
+        decomposition = [d + q, t + rest]
     elif method_name == 'KBC-N4':
-        decomposition = [d+t+q, rest]
+        decomposition = [d + t + q, rest]
     else:
         raise ValueError("Unknown model. Supported models KBC-Nx where x in (1,2,3,4)")
 
diff --git a/methods/cumulantbased.py b/methods/cumulantbased.py
index bbe3bc8449ad5b38508e498a9af369c5af600f1c..9756f34c5a77c0042c3ef79f2e702127c98ace66 100644
--- a/methods/cumulantbased.py
+++ b/methods/cumulantbased.py
@@ -134,7 +134,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
                                                    post_collision_subexpressions=False, include_force_terms=True):
         def tuple_to_symbol(exp, prefix):
             dim = len(exp)
-            format_string = prefix + "_" + "_".join(["%d"]*dim)
+            format_string = prefix + "_" + "_".join(["%d"] * dim)
             return sp.Symbol(format_string % exp)
 
         def substitute_conserved_quantities(expressions, cqe):
diff --git a/methods/entropic.py b/methods/entropic.py
index 7569848a6dbff26f87e40ac1e4d2c98161556deb..862403e6a0b214ba8cca36fd548c3900bbf626cf 100644
--- a/methods/entropic.py
+++ b/methods/entropic.py
@@ -60,7 +60,8 @@ def add_entropy_condition(collision_rule, omega_output_field=None):
     const_part = decomposition.constant_exprs()
     for update_eq in collision_rule.main_assignments:
         index = collision_rule.method.post_collision_pdf_symbols.index(update_eq.lhs)
-        new_eq = Assignment(update_eq.lhs, const_part[index] + omega_s * ds_symbols[index] + omega_h * dh_symbols[index])
+        new_eq = Assignment(update_eq.lhs,
+                            const_part[index] + omega_s * ds_symbols[index] + omega_h * dh_symbols[index])
         new_update_equations.append(new_eq)
     new_collision_rule = collision_rule.copy(new_update_equations, collision_rule.subexpressions + subexpressions)
     new_collision_rule.simplification_hints['entropic'] = True
@@ -118,26 +119,26 @@ def add_iterative_entropy_condition(collision_rule, free_omega=None, newton_iter
 
     # 2) get equilibrium from method and define subexpressions for it
     eq_terms = [eq.rhs for eq in collision_rule.method.get_equilibrium().main_assignments]
-    eq_symbols = sp.symbols("entropicFeq_:%d" % (len(eq_terms,)))
+    eq_symbols = sp.symbols("entropicFeq_:%d" % (len(eq_terms, )))
     eq_subexpressions = [Assignment(a, b) for a, b in zip(eq_symbols, eq_terms)]
 
     # 3) find coefficients of entropy derivatives
     entropy_diff = sp.diff(discrete_approx_entropy(rr_polynomials, eq_symbols), free_omega)
     coefficients_first_diff = [c.expand() for c in reversed(sp.poly(entropy_diff, free_omega).all_coeffs())]
-    sym_coeff_diff1 = sp.symbols("entropicDiffCoeff_:%d" % (len(coefficients_first_diff,)))
+    sym_coeff_diff1 = sp.symbols("entropicDiffCoeff_:%d" % (len(coefficients_first_diff, )))
     coefficient_eqs = [Assignment(a, b) for a, b in zip(sym_coeff_diff1, coefficients_first_diff)]
-    sym_coeff_diff2 = [(i+1) * coeff for i, coeff in enumerate(sym_coeff_diff1[1:])]
+    sym_coeff_diff2 = [(i + 1) * coeff for i, coeff in enumerate(sym_coeff_diff1[1:])]
 
     # 4) define Newtons method update iterations
     newton_iteration_equations = []
     intermediate_omegas = [sp.Symbol("omega_iter_%i" % (i,)) for i in range(newton_iterations + 1)]
     intermediate_omegas[0] = initial_value
     intermediate_omegas[-1] = free_omega
-    for omega_idx in range(len(intermediate_omegas)-1):
+    for omega_idx in range(len(intermediate_omegas) - 1):
         rhs_omega = intermediate_omegas[omega_idx]
-        lhs_omega = intermediate_omegas[omega_idx+1]
-        diff1_poly = sum([coeff * rhs_omega**i for i, coeff in enumerate(sym_coeff_diff1)])
-        diff2_poly = sum([coeff * rhs_omega**i for i, coeff in enumerate(sym_coeff_diff2)])
+        lhs_omega = intermediate_omegas[omega_idx + 1]
+        diff1_poly = sum([coeff * rhs_omega ** i for i, coeff in enumerate(sym_coeff_diff1)])
+        diff2_poly = sum([coeff * rhs_omega ** i for i, coeff in enumerate(sym_coeff_diff2)])
         newton_eq = Assignment(lhs_omega, rhs_omega - diff1_poly / diff2_poly)
         newton_iteration_equations.append(newton_eq)
 
@@ -179,7 +180,7 @@ def discrete_approx_entropy(func, reference):
     .. math ::
         S = - \sum_i f_i \left( \frac{f_i}{r_i} - 1 \right)
     """
-    return -sum([f_i * ((f_i / r_i)-1) for f_i, r_i in zip(func, reference)])
+    return -sum([f_i * ((f_i / r_i) - 1) for f_i, r_i in zip(func, reference)])
 
 
 def _get_entropy_maximizing_omega(omega_s, f_eq, ds, dh):
diff --git a/methods/momentbased.py b/methods/momentbased.py
index 42c545c354d814ef0733d03016ef60b0c153338e..55c8b0f769c1206849326648ad11146269c43e35 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -76,9 +76,9 @@ class MomentBasedLbMethod(AbstractLbMethod):
 
     def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
         relaxation_matrix = sp.eye(len(self.relaxation_rates))
-        return self._get_collision_rule_with_relaxation_matrix(relaxation_matrix,
-                                                               conserved_quantity_equations=conserved_quantity_equations,
-                                                               include_force_terms=include_force_terms)
+        return self._collision_rule_with_relaxation_matrix(relaxation_matrix,
+                                                           conserved_quantity_equations=conserved_quantity_equations,
+                                                           include_force_terms=include_force_terms)
 
     def get_equilibrium_terms(self):
         equilibrium = self.get_equilibrium()
@@ -87,8 +87,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
     def get_collision_rule(self, conserved_quantity_equations=None):
         d = sp.diag(*self.relaxation_rates)
         relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d)
-        ac = self._get_collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
-                                                             True, conserved_quantity_equations)
+        ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
+                                                         True, conserved_quantity_equations)
         return ac
 
     def set_zeroth_moment_relaxation_rate(self, relaxation_rate):
@@ -170,8 +170,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
             weights.append(value)
         return weights
 
-    def _get_collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
-                                                   conserved_quantity_equations=None):
+    def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
+                                               conserved_quantity_equations=None):
         f = sp.Matrix(self.pre_collision_pdf_symbols)
         pdf_to_moment = self.moment_matrix
         m_eq = sp.Matrix(self.moment_equilibrium_values)
diff --git a/phasefield/cahn_hilliard_lbm.py b/phasefield/cahn_hilliard_lbm.py
index 6039420136e5bbb2966cd0da1916cc1f053b2418..8481e4ab08580d0b0835422f587a55463b378c4f 100644
--- a/phasefield/cahn_hilliard_lbm.py
+++ b/phasefield/cahn_hilliard_lbm.py
@@ -42,4 +42,3 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
     rho = sp.Symbol("rho")
     equilibrium[0] = rho - sp.expand(sum(equilibrium[1:]))
     return create_from_equilibrium(stencil, tuple(equilibrium), relaxation_rate, compressible=True)
-
diff --git a/phasefield/experiments1D.py b/phasefield/experiments1D.py
index 8c44b9aa6053e341d3d25476eb0c87d012378e49..fa7e96b1fd5f149a0f3ae4d04ef27935746258d5 100644
--- a/phasefield/experiments1D.py
+++ b/phasefield/experiments1D.py
@@ -135,7 +135,7 @@ def galilean_invariance_test(pf_step, velocity=0.05, rounds=3, phase0=0, phase1=
     from lbmpy.phasefield.analytical import analytic_interface_profile
 
     domain_size = pf_step.data_handling.shape
-    round_time_steps = int((domain_size[0]+0.25) / velocity)
+    round_time_steps = int((domain_size[0] + 0.25) / velocity)
 
     print("Velocity:", velocity, " Time steps for round:", round_time_steps)
 
@@ -161,7 +161,7 @@ def galilean_invariance_test(pf_step, velocity=0.05, rounds=3, phase0=0, phase1=
 
     capture_profile()
     for rt in range(rounds):
-        print("Running round %d/%d" % (rt+1, rounds))
+        print("Running round %d/%d" % (rt + 1, rounds))
         pf_step.run(round_time_steps)
         capture_profile()
 
diff --git a/phasefield/plot.py b/phasefield/plot.py
deleted file mode 100644
index 6d7b00edfc576054fff6324701ef55b625dd9b05..0000000000000000000000000000000000000000
--- a/phasefield/plot.py
+++ /dev/null
@@ -1 +0,0 @@
-from lbmpy.plot2d import *
diff --git a/phasefield/scenarios.py b/phasefield/scenarios.py
index cc372f56cb770819d78a5acd25dccdd684a9307e..020045291c2f6044ace64c6d005443206fdb369f 100644
--- a/phasefield/scenarios.py
+++ b/phasefield/scenarios.py
@@ -1,8 +1,8 @@
 import sympy as sp
 import numpy as np
 from lbmpy.phasefield.phasefieldstep import PhaseFieldStep
-from lbmpy.phasefield.analytical import free_energy_functional_3_phases, free_energy_functional_n_phases, symbolic_order_parameters, \
-    free_energy_functional_n_phases_penalty_term
+from lbmpy.phasefield.analytical import free_energy_functional_3_phases, free_energy_functional_n_phases, \
+    symbolic_order_parameters, free_energy_functional_n_phases_penalty_term
 
 
 def create_three_phase_model(alpha=1, kappa=(0.015, 0.015, 0.015), include_rho=True, **kwargs):
diff --git a/plot2d.py b/plot2d.py
index 1bc774eadde3bf9a929d8954ef7de0fee4f2585c..871e7e4c48af1a0701a663aa3d513f6c22249943 100644
--- a/plot2d.py
+++ b/plot2d.py
@@ -6,14 +6,14 @@ from pystencils.plot2d import *
 
 
 def boundary_handling(boundary_handling_obj, slice_obj=None, boundary_name_to_color=None, show_legend=True):
-    """
-    Shows boundary cells
-
-    :param boundary_handling_obj: instance of :class:`lbmpy.boundaries.BoundaryHandling`
-    :param slice_obj: for 3D boundary handling a slice expression has to be passed here to define the plane that
-                      should be plotted
-    :param boundary_name_to_color: optional dictionary mapping boundary names to colors
-    :param show_legend: if True legend for color->boundary name is added
+    """Image plot that shows boundary markers of a 2D domain slice.
+
+    Args:
+        boundary_handling_obj: instance of :class:`lbmpy.boundaries.BoundaryHandling`
+        slice_obj: for 3D boundary handling a slice expression has to be passed here to define the plane that
+                   should be plotted
+        boundary_name_to_color: optional dictionary mapping boundary names to colors
+        show_legend: if True legend for color->boundary name is added
     """
     import matplotlib.pyplot as plt
 
diff --git a/postprocessing.py b/postprocessing.py
index 241234f1425a92f2958f206d3c2c874efe0d97c3..c8f84378cf31bf5849dda1ce7e5724bef72fe6a2 100644
--- a/postprocessing.py
+++ b/postprocessing.py
@@ -19,4 +19,3 @@ def vorticity_2d(velocity_field):
     grad_y_of_x = np.gradient(velocity_field[:, :, 0], axis=1)
     grad_x_of_y = np.gradient(velocity_field[:, :, 1], axis=0)
     return grad_x_of_y - grad_y_of_x
-
diff --git a/quadratic_equilibrium_construction.py b/quadratic_equilibrium_construction.py
index 1b1a045596b80d8cb06333953ae1d357d724e859..8cc1ed804cc31af0c26c6915871a669133260212 100644
--- a/quadratic_equilibrium_construction.py
+++ b/quadratic_equilibrium_construction.py
@@ -1,5 +1,5 @@
 """
-Method to construct a quadratic equilibrium using a generic quadratic ansatz according to the book of 
+Method to construct a quadratic equilibrium using a generic quadratic ansatz according to the book of
 Wolf-Gladrow, section 5.4
 """
 
@@ -11,7 +11,7 @@ from lbmpy.maxwellian_equilibrium import compressible_to_incompressible_moment_v
 
 
 def generic_equilibrium_ansatz(stencil, u=sp.symbols("u_:3")):
-    """Returns a generic quadratic equilibrium with coefficients A, B, C, D according to 
+    """Returns a generic quadratic equilibrium with coefficients A, B, C, D according to
     Wolf Gladrow Book equation (5.4.1) """
     dim = len(stencil[0])
     u = u[:dim]
@@ -38,7 +38,7 @@ def generic_equilibrium_ansatz_parameters(stencil):
 
 
 def match_generic_equilibrium_ansatz(stencil, equilibrium, u=sp.symbols("u_:3")):
-    """Given a quadratic equilibrium, the generic coefficients A,B,C,D are determined. 
+    """Given a quadratic equilibrium, the generic coefficients A,B,C,D are determined.
     Returns a dict that maps these coefficients to their values. If the equilibrium does not have a
     generic quadratic form, a ValueError is raised"""
     dim = len(stencil[0])
diff --git a/relaxationrates.py b/relaxationrates.py
index 434b19eac745c2ebfa369ef407c50b720d99fd21..d298a2b163d88f51bc0f5ac4d50b10ab0c5b7f55 100644
--- a/relaxationrates.py
+++ b/relaxationrates.py
@@ -8,7 +8,7 @@ def relaxation_rate_from_lattice_viscosity(nu):
 
 
 def lattice_viscosity_from_relaxation_rate(omega):
-    """Computes lattice viscosity from relaxation rate: 
+    """Computes lattice viscosity from relaxation rate:
     :math:`\nu_L=\frac{1}{3}\left(\frac{1}{\omega}-\frac{1}{2}\right)`"""
     return (2 - omega) / (6 * omega)
 
@@ -79,5 +79,3 @@ def default_relaxation_rate_names():
             return sp.Symbol("omega_%d" % (next_index[0],))
 
     return result
-
-
diff --git a/simplificationfactory.py b/simplificationfactory.py
index dbddacace1f96c399f57d5bdd70d516e80e6c327..306ab00cf28f2f53fb085d6b51d5ae1d514142ce 100644
--- a/simplificationfactory.py
+++ b/simplificationfactory.py
@@ -1,6 +1,4 @@
-from functools import partial
 import sympy as sp
-
 from lbmpy.innerloopsplit import create_lbm_split_groups
 from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
 from pystencils.assignment_collection.simplifications import apply_to_all_assignments, \
@@ -9,7 +7,7 @@ from pystencils.assignment_collection.simplifications import apply_to_all_assign
 
 def create_simplification_strategy(lb_method, cse_pdfs=False, cse_global=False, split_inner_loop=False):
     from pystencils.assignment_collection import SimplificationStrategy
-    from lbmpy.methods import MomentBasedLbMethod
+    from lbmpy.methods.momentbased import MomentBasedLbMethod
     from lbmpy.methods.momentbasedsimplifications import replace_second_order_velocity_products, \
         factor_density_after_factoring_relaxation_times, factor_relaxation_rates, cse_in_opposing_directions, \
         replace_common_quadratic_and_constant_term, replace_density_and_velocity
diff --git a/stencils.py b/stencils.py
index 0222d6b0f5cb6bca34ee728eaafd09840d819f02..0c2b062529a1ba35e8b3e885cdb6b5d16ebd2d06 100644
--- a/stencils.py
+++ b/stencils.py
@@ -51,16 +51,16 @@ get_stencil.data = {
                      (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
                      (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
                      (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1)),
-        'braunschweig':  ((0, 0, 0),
-                          (1, 0, 0), (-1, 0, 0),
-                          (0, 1, 0), (0, -1, 0),
-                          (0, 0, 1), (0, 0, -1),
-                          (1, 1, 0), (-1, -1, 0),
-                          (1, -1, 0), (-1, 1, 0),
-                          (1, 0, 1), (-1, 0, -1),
-                          (1, 0, -1), (-1, 0, 1),
-                          (0, 1, 1), (0, -1, -1),
-                          (0, 1, -1), (0, -1, 1)),
+        'braunschweig': ((0, 0, 0),
+                         (1, 0, 0), (-1, 0, 0),
+                         (0, 1, 0), (0, -1, 0),
+                         (0, 0, 1), (0, 0, -1),
+                         (1, 1, 0), (-1, -1, 0),
+                         (1, -1, 0), (-1, 1, 0),
+                         (1, 0, 1), (-1, 0, -1),
+                         (1, 0, -1), (-1, 0, 1),
+                         (0, 1, 1), (0, -1, -1),
+                         (0, 1, -1), (0, -1, 1)),
         'premnath': ((0, 0, 0),
                      (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
                      (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
@@ -174,24 +174,23 @@ def visualize_stencil_2d(stencil, axes=None, figure=None, data=None, textsize='1
             annotation = "$" + sp.latex(annotation) + "$"
         else:
             annotation = str(annotation)
-        axes.text(direction[0]*text_offset, direction[1]*text_offset, annotation, verticalalignment='center', zorder=30,
-                  horizontalalignment='center',
-                  size=textsize, bbox=dict(boxstyle=text_box_style, facecolor='#00b6eb', alpha=0.85, linewidth=0))
+        axes.text(direction[0] * text_offset, direction[1] * text_offset, annotation, verticalalignment='center',
+                  zorder=30, horizontalalignment='center', size=textsize,
+                  bbox=dict(boxstyle=text_box_style, facecolor='#00b6eb', alpha=0.85, linewidth=0))
 
     axes.set_axis_off()
     axes.set_aspect('equal')
-    axes.set_xlim([-text_offset*1.1, text_offset*1.1])
+    axes.set_xlim([-text_offset * 1.1, text_offset * 1.1])
     axes.set_ylim([-text_offset * 1.1, text_offset * 1.1])
 
 
 def visualize_stencil_3d_by_slicing(stencil, slice_axis=2, figure=None, data=None, **kwargs):
-    """
-    Visualizes a 3D, first-neighborhood stencil by plotting 3 slices along a given axis
+    """Visualizes a 3D, first-neighborhood stencil by plotting 3 slices along a given axis.
 
-    :param stencil: stencil as sequence of directions
-    :param slice_axis: 0, 1, or 2 indicating the axis to slice through
-    :param data: optional data to print as text besides the arrows
-    :return:
+    Args:
+        stencil: stencil as sequence of directions
+        slice_axis: 0, 1, or 2 indicating the axis to slice through
+        data: optional data to print as text besides the arrows
     """
     import matplotlib.pyplot as plt
 
@@ -202,7 +201,7 @@ def visualize_stencil_3d_by_slicing(stencil, slice_axis=2, figure=None, data=Non
     if figure is None:
         figure = plt.gcf()
 
-    axes = [figure.add_subplot(1, 3, i+1) for i in range(3)]
+    axes = [figure.add_subplot(1, 3, i + 1) for i in range(3)]
     splitted_directions = [[], [], []]
     splitted_data = [[], [], []]
     axes_names = ['x', 'y', 'z']
@@ -216,7 +215,7 @@ def visualize_stencil_3d_by_slicing(stencil, slice_axis=2, figure=None, data=Non
     for i in range(3):
         visualize_stencil_2d(splitted_directions[i], axes=axes[i], data=splitted_data[i], **kwargs)
     for i in [-1, 0, 1]:
-        axes[i+1].set_title("Cut at %s=%d" % (axes_names[slice_axis], i))
+        axes[i + 1].set_title("Cut at %s=%d" % (axes_names[slice_axis], i))
 
 
 def visualize_stencil_3d(stencil, figure=None, axes=None, data=None, textsize='8'):
@@ -262,7 +261,7 @@ def visualize_stencil_3d(stencil, figure=None, axes=None, data=None, textsize='8
 
     for d, annotation in zip(stencil, data):
         assert len(d) == 3, "Works only for 3D stencils"
-        if not(d[0] == 0 and d[1] == 0 and d[2] == 0):
+        if not (d[0] == 0 and d[1] == 0 and d[2] == 0):
             if d[0] == 0:
                 color = '#348abd'
             elif d[1] == 0:
@@ -281,11 +280,11 @@ def visualize_stencil_3d(stencil, figure=None, axes=None, data=None, textsize='8
             else:
                 annotation = str(annotation)
 
-            axes.text(d[0]*text_offset, d[1]*text_offset, d[2]*text_offset,
+            axes.text(d[0] * text_offset, d[1] * text_offset, d[2] * text_offset,
                       annotation, verticalalignment='center', zorder=30,
                       size=textsize, bbox=dict(boxstyle=text_box_style, facecolor='#777777', alpha=0.6, linewidth=0))
 
-    axes.set_xlim([-text_offset*1.1, text_offset*1.1])
+    axes.set_xlim([-text_offset * 1.1, text_offset * 1.1])
     axes.set_ylim([-text_offset * 1.1, text_offset * 1.1])
     axes.set_zlim([-text_offset * 1.1, text_offset * 1.1])
     axes.set_axis_off()
diff --git a/updatekernels.py b/updatekernels.py
index dc664ff9f598ff9233ec3c6e1707a58e851c6768..a63d0d115877f9890cc21ab7481a191a162ee1a7 100644
--- a/updatekernels.py
+++ b/updatekernels.py
@@ -6,7 +6,7 @@ from pystencils.assignment_collection.assignment_collection import AssignmentCol
 from pystencils.field import create_numpy_array_with_layout, layout_string_to_tuple
 from pystencils.sympyextensions import fast_subs
 from lbmpy.methods.abstractlbmethod import LbmCollisionRule
-from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAccessor, CollideOnlyInplaceAccessor
+from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
 
 
 # -------------------------------------------- LBM Kernel Creation -----------------------------------------------------