diff --git a/creationfunctions.py b/creationfunctions.py
index eb85f66578d9d4ff463f3142a92835713e51da01..15cf5e5a14ea4e5a9d453d6b8970f50763acf4f6 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -2,12 +2,8 @@ r"""
 Creating LBM kernels
 ====================
 
-Parameters
-----------
-
 The following list describes common parameters for the creation functions. They have to be passed as keyword parameters.
 
-
 Method parameters
 ^^^^^^^^^^^^^^^^^
 
@@ -173,6 +169,7 @@ from lbmpy.updatekernels import create_lbm_kernel, create_stream_pull_with_outpu
 
 
 def create_lb_function(ast=None, optimization={}, **kwargs):
+    """Creates a Python function for the LB method"""
     params, opt_params = update_with_default_parameters(kwargs, optimization)
 
     if ast is None:
@@ -188,6 +185,7 @@ def create_lb_function(ast=None, optimization={}, **kwargs):
 
 
 def create_lb_ast(update_rule=None, optimization={}, **kwargs):
+    """Creates a pystencils AST for the LB method"""
     params, opt_params = update_with_default_parameters(kwargs, optimization)
 
     if update_rule is None:
@@ -207,6 +205,7 @@ def create_lb_ast(update_rule=None, optimization={}, **kwargs):
 
 @disk_cache_no_fallback
 def create_lb_update_rule(collision_rule=None, optimization={}, **kwargs):
+    """Creates an update rule (list of Assignments) for a LB method that describe a full sweep"""
     params, opt_params = update_with_default_parameters(kwargs, optimization)
 
     if collision_rule is None:
@@ -263,6 +262,7 @@ def create_lb_update_rule(collision_rule=None, optimization={}, **kwargs):
 
 @disk_cache_no_fallback
 def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
+    """Creates a collision rule (list of Assignments) for a LB method describing the collision operator (no stream)"""
     params, opt_params = update_with_default_parameters(kwargs, optimization)
 
     if lb_method is None:
@@ -291,6 +291,7 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
 
 
 def create_lb_method(**params):
+    """Creates a LB method, defined by moments/cumulants for collision space, equilibrium and relaxation rates."""
     params, _ = update_with_default_parameters(params, {}, fail_on_unknown_parameter=False)
 
     if isinstance(params['stencil'], tuple) or isinstance(params['stencil'], list):
diff --git a/methods/__init__.py b/methods/__init__.py
index fa6975c9444c04d91be267274df2539cde1a0b44..3e068b6f0eef5c41de27d172f4cb3cc3c5ee1135 100644
--- a/methods/__init__.py
+++ b/methods/__init__.py
@@ -1,9 +1,12 @@
-from lbmpy.methods.momentbased import RelaxationInfo
+from lbmpy.methods.momentbased import RelaxationInfo, AbstractLbMethod, AbstractConservedQuantityComputation, \
+    MomentBasedLbMethod
+from .conservedquantitycomputation import DensityVelocityComputation
 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',
+__all__ = ['RelaxationInfo', 'AbstractLbMethod',
+           'AbstractConservedQuantityComputation', 'DensityVelocityComputation', 'MomentBasedLbMethod',
            '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/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index fd3c87834fe92ccc5b4d0533117929d87f6ee7c6..84e2df2c42d9cbbea65cc4f41a2caf1c887a2e13 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -17,7 +17,7 @@ class AbstractConservedQuantityComputation(abc.ABC):
     For example in zero centered hydrodynamic schemes with force model, the density has
     to be decreased by one, and the given velocity has to be shifted dependent on the force.
 
-    .. image:: moment_shift.svg
+    .. image:: /img/moment_shift.svg
 
     """
 
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 328bd409cbe46293fe2e74edd09067be4960ef83..c12cd5d16ca01ef4ba8451d4ccb8c887c1a55730 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -23,26 +23,26 @@ from lbmpy.relaxationrates import relaxation_rate_from_magic_number, default_rel
 def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
                                                force_model=None, equilibrium_order=2,
                                                cumulant=False, c_s_sq=sp.Rational(1, 3)):
-    r"""
-    Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate. These moments are
-    relaxed against the moments of the discrete Maxwellian distribution.
+    r"""Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate.
+
+    These moments are relaxed against the moments of the discrete Maxwellian distribution.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See `func:lbmpy.stencils.get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
         moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
-                                    represented by an exponent tuple or in polynomial form
-                                    (see `lbmpy.moments`), is mapped to a relaxation rate.
+                                        represented by an exponent tuple or in polynomial form
+                                        (see `lbmpy.moments`), is mapped to a relaxation rate.
         compressible: incompressible LBM methods split the density into :math:`\rho = \rho_0 + \Delta \rho`
-        where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
-        This approximates the incompressible Navier-Stokes equations better than the standard
-        compressible model.
+                      where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
+                      This approximates the incompressible Navier-Stokes equations better than the standard
+                      compressible model.
         force_model: force model instance, or None if no external forces
         equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
         cumulant: if True relax cumulants instead of moments
         c_s_sq: Speed of sound squared
 
     Returns:
-        :class:`lbmpy.methods.MomentBasedLbMethod` instance
+        `lbmpy.methods.MomentBasedLbMethod` instance
     """
     if isinstance(stencil, str):
         stencil = get_stencil(stencil)
@@ -357,16 +357,17 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
     Returns a orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
     These MRT methods are just one specific version - there are many MRT methods possible for all these stencils
     which differ by the linear combination of moments and the grouping into equal relaxation times.
-    To create a generic MRT method use :func:`lbmpy.methods.create_with_discrete_maxwellian_eq_moments`
+    To create a generic MRT method use `create_with_discrete_maxwellian_eq_moments`
 
-    :param stencil: nested tuple defining the discrete velocity space. See `func:lbmpy.stencils.get_stencil`
-    :param relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
-                                 time. The default returns:
+    Args:
+        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+        relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
+                              time. The default returns:
 
-                                    - 0 for moments of order 0 and 1 (conserved)
-                                    - :math:`\omega`: from moments of order 2 (rate that determines viscosity)
-                                    - numbered :math:`\omega_i` for the rest
-    :param maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
+                                 - 0 for moments of order 0 and 1 (conserved)
+                                 - :math:`\omega`: from moments of order 2 (rate that determines viscosity)
+                                 - numbered :math:`\omega_i` for the rest
+        maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                                                used to compute the equilibrium moments
     """
     if relaxation_rate_getter is None:
diff --git a/parameterization.py b/parameterization.py
index 5cb3ad25a698f0e9a78e4164b26ac444f34542fa..b6392188da9e01003188ce945f74377a20352ed8 100644
--- a/parameterization.py
+++ b/parameterization.py
@@ -1,3 +1,4 @@
+from collections import namedtuple
 import ipywidgets.widgets as widgets
 from IPython.display import display
 from ipywidgets.widgets import FloatText, Label, VBox, HBox, Select, BoundedFloatText, Button
@@ -100,13 +101,13 @@ class ScalingWidget:
             self._update_dt_from_relaxation_rate_viscosity_and_dx()
             self._update_lattice_velocity_from_dx_dt_and_physical_velocity()
         elif self.scaling_type.value == r'acoustic (fixed dt)':
-            self._update_omega_from_viscosity_and_dt_dx()
-            self._update_dt_from_dx_and_lattice_velocity()
-            self.max_lattice_velocity.disabled = False
-        elif self.scaling_type.value == r'fixed lattice velocity':
             self._update_omega_from_viscosity_and_dt_dx()
             self.dt.disabled = False
             self._update_lattice_velocity_from_dx_dt_and_physical_velocity()
+        elif self.scaling_type.value == r'fixed lattice velocity':
+            self._update_omega_from_viscosity_and_dt_dx()
+            self._update_dt_from_dx_and_lattice_velocity()
+            self.max_lattice_velocity.disabled = False
         else:
             raise ValueError("Unknown Scaling Type")
 
@@ -118,10 +119,9 @@ class ScalingWidget:
     def _update_re(self, _=None):
         if self.kinematic_viscosity.value == 0:
             return
-        L = self.physical_length.value
-        u = self.max_physical_velocity.value
-        nu = self.kinematic_viscosity.value * 1e-6
-        self.re.value = round(L * u / nu, 7)
+        viscosity = self.kinematic_viscosity.value * 1e-6
+        re = reynolds_number(self.physical_length.value, self.max_physical_velocity.value, viscosity)
+        self.re.value = round(re, 7)
 
     def _on_dx_change(self, _):
         if self.processing_update:
@@ -159,3 +159,131 @@ class ScalingWidget:
         </style>
         """))
         return self.form
+
+
+def reynolds_number(length, velocity, kinematic_viscosity):
+    """Computes the Reynolds number.
+
+    All arguments have to be in the same set of units, for example all in SI units.
+    """
+    return length / kinematic_viscosity * velocity
+
+
+class Scaling:
+    """Class to convert physical parameters into lattice units.
+
+    The scaling is created with the central physical parameters length, velocity and viscosity
+    These parameters fix the Reynolds number.
+
+    Grid information is specified by fixing the number of cells to resolve the physical length.
+    The three lattice parameters relaxation rate, lattice velocity and time step length can be computed, if one of
+    them is given (see the '*_scaling' methods, that receive one parameter and compute the other two.
+
+    Args:
+        physical_length: typical physical length [m]
+        physical_velocity: (maximum) physical velocity [m/s]
+                            usually the maximum velocity occurring in the simulation domain should be passed here,
+                            such that the fixed lattice velocity scaling can ensure a maximum lattice velocity
+        kinematic_viscosity: kinematic viscosity in physical units [m*m/s]
+        cells_per_length: number of cells to resolve the physical length with
+    """
+    def __init__(self, physical_length, physical_velocity, kinematic_viscosity, cells_per_length):
+        self._physical_length = physical_length
+        self._physical_velocity = physical_velocity
+        self._kinematic_viscosity = kinematic_viscosity
+        self.cells_per_length = cells_per_length
+        self._reynolds_number = None
+        self._parameter_update()
+
+    def _parameter_update(self):
+        self._reynolds_number = reynolds_number(self._physical_length, self._physical_velocity,
+                                                self._kinematic_viscosity)
+
+    @property
+    def dx(self):
+        return self.physical_length / self.cells_per_length
+
+    @property
+    def reynolds_number(self):
+        return self._reynolds_number
+
+    @property
+    def physical_velocity(self):
+        return self._physical_velocity
+
+    @physical_velocity.setter
+    def physical_velocity(self, val):
+        self._physical_velocity = val
+        self._parameter_update()
+
+    @property
+    def kinematic_viscosity(self):
+        return self._kinematic_viscosity
+
+    @kinematic_viscosity.setter
+    def kinematic_viscosity(self, val):
+        self._kinematic_viscosity = val
+        self._parameter_update()
+
+    @property
+    def physical_length(self):
+        return self._physical_length
+
+    @physical_length.setter
+    def physical_length(self, val):
+        self._physical_length = val
+        self._parameter_update()
+
+    def fixed_lattice_velocity_scaling(self, lattice_velocity=0.1):
+        """Computes relaxation rate and time step length from lattice velocity.
+
+        Check the returned relaxation rate! If it is very close to 2, the simulation might get unstable.
+        In that case increase the resolution (more cells_per_length)
+        All physical quantities have to be passed in the same set of units.
+
+        Args:
+            lattice_velocity: Lattice velocity corresponding to physical_velocity.
+                                  Maximum velocity should not be larger than 0.1 i.e. the fluid should not move
+                                  faster than 0.1 cells per time step
+
+        Returns:
+            relaxation_rate, dt
+        """
+        dt = lattice_velocity / self.physical_velocity * self.dx
+        lattice_viscosity = self.kinematic_viscosity * dt / (self.dx ** 2)
+        relaxation_rate = relaxation_rate_from_lattice_viscosity(lattice_viscosity)
+        ResultType = namedtuple("FixedLatticeVelocityScalingResult", ['relaxation_rate', 'dt'])
+        return ResultType(relaxation_rate, dt)
+
+    def diffusive_scaling(self, relaxation_rate):
+        """Computes time step length and lattice velocity from relaxation rate
+
+        For stable simulations the lattice velocity should be smaller than 0.1
+
+        Args:
+            relaxation_rate: relaxation rate (between 0 and 2)
+
+        Returns:
+            dt, lattice velocity
+        """
+        lattice_viscosity = lattice_viscosity_from_relaxation_rate(relaxation_rate)
+        dx = self.physical_length / self.cells_per_length
+        dt = lattice_viscosity / self.kinematic_viscosity * dx ** 2
+        lattice_velocity = dt / dx * self.physical_velocity
+        ResultType = namedtuple("DiffusiveScalingResult", ['dt', 'lattice_velocity'])
+        return ResultType(dt, lattice_velocity)
+
+    def acoustic_scaling(self, dt):
+        """Computes relaxation rate and lattice velocity from time step length.
+
+        Args:
+            dt: time step length
+
+        Returns:
+            relaxation_rate, lattice_velocity
+        """
+        lattice_velocity = dt / self.dx * self.physical_velocity
+        lattice_viscosity = self.kinematic_viscosity * dt / (self.dx ** 2)
+        relaxation_rate = relaxation_rate_from_lattice_viscosity(lattice_viscosity)
+        ResultType = namedtuple("AcousticScalingResult", ['relaxation_rate', 'lattice_velocity'])
+        return ResultType(relaxation_rate, lattice_velocity)
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index f9061b9aa1b807d8979ffbd0ae17ec6ca5adac22..0ba3b8668d4404e68117b2ecc13c864e0a338a2f 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -58,7 +58,7 @@ def force_kernel_using_pressure_tensor(force_field, pressure_tensor_field, extra
 
 
 def cahn_hilliard_fd_eq(phase_idx, phi, mu, velocity, mobility, dx, dt):
-    from pystencils.fd.finitedifferences import transient, advection, diffusion
+    from pystencils.fd import transient, advection, diffusion
     cahn_hilliard = transient(phi, phase_idx) + advection(phi, velocity, phase_idx) - diffusion(mu, mobility, phase_idx)
     return Discretization2ndOrder(dx, dt)(cahn_hilliard)