diff --git a/lbmpy/chapman_enskog/chapman_enskog_steady_state.py b/lbmpy/chapman_enskog/chapman_enskog_steady_state.py
index fcbc8aeb3363e443ae23b12ef19150c90329b61e..65496f2fd6161391212755a07e9d8970988a9f42 100644
--- a/lbmpy/chapman_enskog/chapman_enskog_steady_state.py
+++ b/lbmpy/chapman_enskog/chapman_enskog_steady_state.py
@@ -272,8 +272,9 @@ class SteadyStateChapmanEnskogAnalysisSRT:
 
         The bulk viscosity is predicted differently than by the normal Navier Stokes analysis...why??
 
-        :param coordinate: which momentum equation to use i.e. x,y or z, to approximate Navier Stokes
-                           all have to return the same result
+        Args:
+            coordinate: which momentum equation to use i.e. x,y or z, to approximate Navier Stokes
+                        all have to return the same result
         """
         dim = self.method.dim
 
diff --git a/lbmpy/maxwellian_equilibrium.py b/lbmpy/maxwellian_equilibrium.py
index 98b73d88f2d61517720d406a886ac798069ffb16..4cf856977cd424bd2c2ae47ecde1441f9f74f133 100644
--- a/lbmpy/maxwellian_equilibrium.py
+++ b/lbmpy/maxwellian_equilibrium.py
@@ -127,11 +127,12 @@ def continuous_maxwellian_equilibrium(dim=3, rho=sp.Symbol("rho"),
     """
     Returns sympy expression of Maxwell Boltzmann distribution
 
-    :param dim: number of space dimensions
-    :param rho: sympy symbol for the density
-    :param u: symbols for macroscopic velocity (expected value for velocity)
-    :param v: symbols for particle velocity
-    :param c_s_sq: symbol for speed of sound squared, defaults to symbol c_s**2
+    Args:
+        dim: number of space dimensions
+        rho: sympy symbol for the density
+        u: symbols for macroscopic velocity (expected value for velocity)
+        v: symbols for particle velocity
+        c_s_sq: symbol for speed of sound squared, defaults to symbol c_s**2
     """
     u = u[:dim]
     v = v[:dim]
@@ -150,13 +151,14 @@ def get_moments_of_continuous_maxwellian_equilibrium(moments, dim, rho=sp.Symbol
     """
     Computes moments of the continuous Maxwell Boltzmann equilibrium distribution
 
-    :param moments: moments to compute, either in polynomial or exponent-tuple form
-    :param dim: dimension (2 or 3)
-    :param rho: symbol or value for the density
-    :param u: symbols or values for the macroscopic velocity
-    :param c_s_sq: symbol for speed of sound squared, defaults to symbol c_s**2
-    :param order: if this parameter is not None, terms that have a higher polynomial order in the macroscopic velocity
-                  are removed
+    Args:
+        moments: moments to compute, either in polynomial or exponent-tuple form
+        dim: dimension (2 or 3)
+        rho: symbol or value for the density
+        u: symbols or values for the macroscopic velocity
+        c_s_sq: symbol for speed of sound squared, defaults to symbol c_s**2
+        order: if this parameter is not None, terms that have a higher polynomial order in the macroscopic velocity
+               are removed
 
     >>> get_moments_of_continuous_maxwellian_equilibrium( ( (0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0) ), dim=3 )
     [rho, rho*u_0, rho*u_1, rho*u_2, rho*(c_s**2 + u_0**2)]
diff --git a/lbmpy/methods/conservedquantitycomputation.py b/lbmpy/methods/conservedquantitycomputation.py
index d34939b7e89477ebe2b2f4f531978a9bc2acb184..4d9c8afadbfb99cab8bc02cf701cd7b7acfa2fa9 100644
--- a/lbmpy/methods/conservedquantitycomputation.py
+++ b/lbmpy/methods/conservedquantitycomputation.py
@@ -54,7 +54,8 @@ class AbstractConservedQuantityComputation(abc.ABC):
         of the pdfs.
         For hydrodynamic LBM schemes this is usually the density and velocity.
 
-        :param pdfs: values or symbols for the pdf values
+        Args:
+            pdfs: values or symbols for the pdf values
         """
 
     @abc.abstractmethod
@@ -63,9 +64,10 @@ class AbstractConservedQuantityComputation(abc.ABC):
         Returns an equation collection that defines conserved quantities for output. These conserved quantities might
         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
+        Args:
+            pdfs: values for the pdf entries
+            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
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index e5a37bffae5afc51d4907d53271de977b3ab3be1..93174eb7ed3bbfee0e4b92b048b680ed1da65e46 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -294,13 +294,14 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
     one for higher order moments. In entropic models this second relaxation rate is chosen subject to an entropy
     condition. Which moments are relaxed by which rate is determined by the method_name
 
-    :param dim: 2 or 3, leads to stencil D2Q9 or D3Q27
-    :param shear_relaxation_rate: relaxation rate that determines viscosity
-    :param higher_order_relaxation_rate: relaxation rate for higher order moments
-    :param method_name: string 'KBC-Nx' where x can be an number from 1 to 4, for details see
-                       "Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows"
-    :param maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
-                           used to compute the equilibrium moments
+    Args:
+        dim: 2 or 3, leads to stencil D2Q9 or D3Q27
+        shear_relaxation_rate: relaxation rate that determines viscosity
+        higher_order_relaxation_rate: relaxation rate for higher order moments
+        method_name: string 'KBC-Nx' where x can be an number from 1 to 4, for details see
+                    "Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows"
+        maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
+                        used to compute the equilibrium moments
     """
     def product(iterable):
         return reduce(operator.mul, iterable, 1)
diff --git a/lbmpy/methods/entropic.py b/lbmpy/methods/entropic.py
index 2e400a093bed6ed1297f72245874093ccdc5aa23..79a99296c7075f3a611cd04c4019f9ed5bc909d6 100644
--- a/lbmpy/methods/entropic.py
+++ b/lbmpy/methods/entropic.py
@@ -15,9 +15,12 @@ def add_entropy_condition(collision_rule, omega_output_field=None):
     The entropy is approximated such that the optimality condition can be written explicitly, no Newton iterations
     have to be done.
 
-    :param collision_rule: collision rule with two relaxation times
-    :param omega_output_field: pystencils field where computed omegas are stored
-    :return: new collision rule which only one relaxation rate
+    Args:
+        collision_rule: collision rule with two relaxation times
+        omega_output_field: pystencils field where computed omegas are stored
+
+    Returns:
+        new collision rule which only one relaxation rate
     """
     if collision_rule.method.conserved_quantity_computation.zero_centered_pdfs:
         raise NotImplementedError("Entropic Methods only implemented for models where pdfs are centered around 1. "
@@ -81,13 +84,16 @@ def add_iterative_entropy_condition(collision_rule, free_omega=None, newton_iter
 
     A fixed number of Newton iterations is used to determine the maximum entropy relaxation rate.
 
-    :param collision_rule: collision rule with two relaxation times
-    :param free_omega: relaxation rate which should be determined by entropy condition. If left to None, the
-                      relaxation rate is automatically detected, which works only if there are 2 relaxation times
-    :param newton_iterations: (integer) number of newton iterations
-    :param initial_value: initial value of the relaxation rate
-    :param omega_output_field: pystencils field where computed omegas are stored
-    :return: new collision rule which only one relaxation rate
+    Args:
+        collision_rule: collision rule with two relaxation times
+        free_omega: relaxation rate which should be determined by entropy condition. If left to None, the
+                   relaxation rate is automatically detected, which works only if there are 2 relaxation times
+        newton_iterations: (integer) number of newton iterations
+        initial_value: initial value of the relaxation rate
+        omega_output_field: pystencils field where computed omegas are stored
+
+    Returns:
+        new collision rule which only one relaxation rate
     """
 
     if collision_rule.method.conserved_quantity_computation.zero_centered_pdfs:
diff --git a/lbmpy/methods/momentbased.py b/lbmpy/methods/momentbased.py
index 2c197c4fcc555a4eafa03c8fb9b68602cbd45942..67d0ec2d6fdbe3a1cb873bc11fa726a23ca182c1 100644
--- a/lbmpy/methods/momentbased.py
+++ b/lbmpy/methods/momentbased.py
@@ -17,14 +17,15 @@ class MomentBasedLbMethod(AbstractLbMethod):
         space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
         equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
 
-        :param stencil: see :func:`lbmpy.stencils.get_stencil`
-        :param moment_to_relaxation_info_dict: a dictionary mapping moments in either tuple or polynomial formulation
-                                           to a RelaxationInfo, which consists of the corresponding equilibrium moment
-                                           and a relaxation rate
-        :param conserved_quantity_computation: instance of :class:`lbmpy.methods.AbstractConservedQuantityComputation`.
-                                             This determines how conserved quantities are computed, and defines
-                                             the symbols used in the equilibrium moments like e.g. density and velocity
-        :param force_model: force model instance, or None if no forcing terms are required
+        Args:
+            stencil: see :func:`lbmpy.stencils.get_stencil`
+            moment_to_relaxation_info_dict: a dictionary mapping moments in either tuple or polynomial formulation
+                                        to a RelaxationInfo, which consists of the corresponding equilibrium moment
+                                        and a relaxation rate
+            conserved_quantity_computation: instance of :class:`lbmpy.methods.AbstractConservedQuantityComputation`.
+                                          This determines how conserved quantities are computed, and defines
+                                          the symbols used in the equilibrium moments like e.g. density and velocity
+            force_model: force model instance, or None if no forcing terms are required
         """
         assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
         super(MomentBasedLbMethod, self).__init__(stencil)
diff --git a/lbmpy/moments.py b/lbmpy/moments.py
index 0a5da08c2e724287cfc58e5cb86b47ab200e871c..c5fe1f76213d7f5f5a61d8fd2f18d776648f99ff 100644
--- a/lbmpy/moments.py
+++ b/lbmpy/moments.py
@@ -418,8 +418,10 @@ def get_default_moment_set_for_stencil(stencil):
 def extract_monomials(sequence_of_polynomials, dim=3):
     """
     Returns a set of exponent tuples of all monomials contained in the given set of polynomials
-    :param sequence_of_polynomials: sequence of polynomials in the MOMENT_SYMBOLS
-    :param dim: length of returned exponent tuples
+
+    Args:
+        sequence_of_polynomials: sequence of polynomials in the MOMENT_SYMBOLS
+        dim: length of returned exponent tuples
 
     >>> x, y, z = MOMENT_SYMBOLS
     >>> extract_monomials([x**2 + y**2 + y, y + y**2])
@@ -437,8 +439,10 @@ def extract_monomials(sequence_of_polynomials, dim=3):
 def monomial_to_polynomial_transformation_matrix(monomials, polynomials):
     """
     Returns a transformation matrix from a monomial to a polynomial representation
-    :param monomials: sequence of exponent tuples
-    :param polynomials: sequence of polynomials in the MOMENT_SYMBOLS
+
+    Args:
+        monomials: sequence of exponent tuples
+        polynomials: sequence of polynomials in the MOMENT_SYMBOLS
 
     >>> x, y, z = MOMENT_SYMBOLS
     >>> polys = [7 * x**2 + 3 * x + 2 * y **2, \
diff --git a/lbmpy/phasefield/analytical.py b/lbmpy/phasefield/analytical.py
index dfd23718196b30658330c982ac030ebb3c6af569..db8a88a95877e1610dd869252f38bdbbd34fd5dc 100644
--- a/lbmpy/phasefield/analytical.py
+++ b/lbmpy/phasefield/analytical.py
@@ -250,8 +250,10 @@ def force_from_phi_and_mu(order_parameters, dim, mu=None):
 def substitute_laplacian_by_sum(eq, dim):
     """Substitutes abstract Laplacian represented by ∂∂ by a sum over all dimensions
     i.e. in case of 3D: ∂∂ is replaced by ∂0∂0 + ∂1∂1 + ∂2∂2
-    :param eq: the term where the substitutions should be made
-    :param dim: spatial dimension, in example above, 3
+
+    Args:
+        eq: the term where the substitutions should be made
+        dim: spatial dimension, in example above, 3
     """
     functions = [d.args[0] for d in eq.atoms(Diff)]
     substitutions = {Diff(Diff(op)): sum(Diff(Diff(op, i), i) for i in range(dim))
diff --git a/lbmpy/stencils.py b/lbmpy/stencils.py
index 98bd3b1e68610568c09bba2aa36e7dc8597c11d8..a0d9001c782bb5e7043257f40885731ff828ec9e 100644
--- a/lbmpy/stencils.py
+++ b/lbmpy/stencils.py
@@ -3,11 +3,12 @@ def get_stencil(name, ordering='walberla'):
     Stencils are tuples of discrete velocities. They are commonly labeled in the 'DxQy' notation, where d is the
     dimension (length of the velocity tuples) and y is number of discrete velocities.
 
-    :param name: DxQy notation
-    :param ordering: the LBM literature does not use a common order of the discrete velocities, therefore here
-                     different common orderings are available. All orderings lead to the same method, it just has
-                     to be used consistently. Here more orderings are available to compare intermediate results with
-                     the literature.
+    Args:
+        name: DxQy notation
+        ordering: the LBM literature does not use a common order of the discrete velocities, therefore here
+                  different common orderings are available. All orderings lead to the same method, it just has
+                  to be used consistently. Here more orderings are available to compare intermediate results with
+                  the literature.
     """
     try:
         return get_stencil.data[name.upper()][ordering.lower()]