Merge branch 'master' into fix_macroscopic_value_kernels

lbmpy/chapman_enskog/chapman_enskog_steady_state.py

 import functools
 
+import numpy as np
 import sympy as sp
 
 from lbmpy.chapman_enskog.chapman_enskog import (
@@ -150,7 +151,8 @@ class SteadyStateChapmanEnskogAnalysis:
 
                 have_shape = hasattr(arg, 'shape') and hasattr(new_prod, 'shape')
                 if have_shape and arg.shape == new_prod.shape and arg.shape[1] == 1:
-                    new_prod = sp.matrix_multiply_elementwise(new_prod, arg)
+                    # since sympy 1.9 sp.matrix_multiply_elementwise does not work anymore in this case
+                    new_prod = sp.Matrix(np.multiply(new_prod, arg))
                 else:
                     new_prod = arg * new_prod
                 if new_prod == 0:

lbmpy/continuous_distribution_measures.py

@@ -15,11 +15,11 @@ def moment_generating_function(generating_function, symbols, symbols_in_result,
     Computes the moment generating function of a probability distribution. It is defined as:
 
     .. math ::
-        F[f(\mathbf{x})](\mathbf{t}) = \int e^{<\mathbf{x}, \mathbf{t}>} f(x)\; dx
+        F[f(\mathbf{x})](t) = \int e^{<\mathbf{x}, t>} f(\mathbf{x})\; dx
 
     Args:
         generating_function: sympy expression
-        symbols: a sequence of symbols forming the vector x
+        symbols: a sequence of symbols forming the vector :math:`\mathbf{x}`
         symbols_in_result: a sequence forming the vector t
         velocity: if the generating function generates central moments, the velocity needs to be substracted. Thus the
                   velocity symbols need to be passed. All generating functions need to have the same parameters.
@@ -62,7 +62,7 @@ def central_moment_generating_function(func, symbols, symbols_in_result, velocit
     Computes central moment generating func, which is defined as:
 
     .. math ::
-        K( \vec{\Xi} ) = \exp ( - \vec{\Xi} \cdot \vec{u} ) M( \vec{\Xi}.
+        K( \mathbf{\Xi} ) = \exp ( - \mathbf{\Xi} \cdot \mathbf{u} ) M( \mathbf{\Xi} ).
 
     For parameter description see :func:`moment_generating_function`.
     """
@@ -76,7 +76,7 @@ def cumulant_generating_function(func, symbols, symbols_in_result, velocity=None
     Computes cumulant generating func, which is the logarithm of the moment generating func:
 
     .. math ::
-        C(\vec{\Xi}) = \log M(\vec{\Xi})
+        C(\mathbf{\Xi}) = \log M(\mathbf{\Xi})
 
     For parameter description see :func:`moment_generating_function`.
     """

lbmpy/cumulants.py

@@ -102,7 +102,7 @@ def __cumulant_raw_moment_transform(index, dependent_var_dict, outer_function, d
 
 @memorycache(maxsize=16)
 def __get_discrete_cumulant_generating_function(func, stencil, wave_numbers):
-    assert len(stencil) == len(func)
+    assert stencil.Q == len(func)
 
     laplace_transformation = sum([factor * sp.exp(scalar_product(wave_numbers, e)) for factor, e in zip(func, stencil)])
     return sp.ln(laplace_transformation)
@@ -121,7 +121,7 @@ def discrete_cumulant(func, cumulant, stencil):
                   (similar to moment description)
         stencil: sequence of directions
     """
-    assert len(stencil) == len(func)
+    assert stencil.Q == len(func)
 
     dim = len(stencil[0])
     wave_numbers = sp.symbols(f"Xi_:{dim}")
@@ -157,9 +157,9 @@ def cumulants_from_pdfs(stencil, cumulant_indices=None, pdf_symbols=None):
     dim = len(stencil[0])
     if cumulant_indices is None:
         cumulant_indices = moments_up_to_component_order(2, dim=dim)
-    assert len(stencil) == len(cumulant_indices), "Stencil has to have same length as cumulant_indices sequence"
+    assert stencil.Q == len(cumulant_indices), "Stencil has to have same length as cumulant_indices sequence"
     if pdf_symbols is None:
-        pdf_symbols = __get_indexed_symbols(pdf_symbols, "f", range(len(stencil)))
+        pdf_symbols = __get_indexed_symbols(pdf_symbols, "f", range(stencil.Q))
     return {idx: discrete_cumulant(tuple(pdf_symbols), idx, stencil) for idx in cumulant_indices}
 
 

lbmpy/enums.py

@@ -2,37 +2,186 @@ from enum import Enum, auto
 
 
 class Stencil(Enum):
+    """
+    The Stencil enumeration represents all possible lattice Boltzmann stencils that are available in lbmpy.
+    It should be passed to :class:`lbmpy.stencils.LBStenil`. This class then creates a stencils representation
+    containing the concrete neighbour directions as a tuple of tuples.
+
+    The number of spatial dimensions *d* and the number of discrete velocities *q* are stated in the DdQq notation
+    """
     D2Q9 = auto()
+    """
+    A two dimensional stencil using 9 discrete velocities.
+    """
     D2V17 = auto()
+    """
+    A two dimensional stencil using 17 discrete velocities. (long range stencil).
+    """
     D2V37 = auto()
+    """
+    A two dimensional stencil using 37 discrete velocities. (long range stencil).
+    """
+    D3Q7 = auto()
+    """
+    A three dimensional stencil using 7 discrete velocities.
+    """
     D3Q15 = auto()
+    """
+    A three dimensional stencil using 15 discrete velocities.
+    """
     D3Q19 = auto()
+    """
+    A three dimensional stencil using 19 discrete velocities.
+    """
     D3Q27 = auto()
+    """
+    A three dimensional stencil using 27 discrete velocities.
+    """
 
 
 class Method(Enum):
+    """
+    The Method enumeration represents all possible lattice Boltzmann collision operators that are available in lbmpy.
+    It should be passed to :class:`lbmpy.creationfunctions.LBMConfig`. The LBM configuration *dataclass* then derives
+    the respective collision equations when passed to the creations functions in the `lbmpy.creationfunctions`
+    module of lbmpy.
+
+    Note here, when using a specific enumeration to derive a particular LBM collision operator,
+    different parameters of the :class:`lbmpy.creationfunctions.LBMConfig` might become necessary.
+    For example, it does not make sense to define *relaxation_rates* for a single relaxation rate method, which
+    is essential for multiple relaxation rate methods. Important specific parameters are listed below to the enum value.
+    A specific creation function is stated for each case which explains these parameters in detail.
+    """
     SRT = auto()
+    """
+    See :func:`lbmpy.methods.create_srt`, 
+    Single relaxation time method
+    """
     TRT = auto()
+    """
+    See :func:`lbmpy.methods.create_trt`, 
+    Two relaxation time, the first relaxation rate is for even moments and determines the
+    viscosity (as in SRT). The second relaxation rate is used for relaxing odd moments and controls the
+    bulk viscosity. For details in the TRT collision operator see :cite:`TRT`
+    """
+    MRT_RAW = auto()
+    """
+    See :func:`lbmpy.methods.create_mrt_raw`, 
+    Non-orthogonal MRT where all relaxation rates can be specified independently, i.e. there are as many relaxation 
+    rates as stencil entries. Look at the generated method in Jupyter to see which moment<->relaxation rate mapping.
+    Originally defined in :cite:`raw_moments`
+    """
     MRT = auto()
+    """
+    See :func:`lbmpy.methods.create_mrt_orthogonal`
+    Orthogonal multi relaxation time model, relaxation rates are used in this order for *shear modes*, *bulk modes*,
+    *third-order modes*, *fourth-order modes*, etc. Requires also a parameter *weighted* that should be `True` if the
+    moments should be orthogonal w.r.t. weighted scalar product using the lattice weights. If `False`, the normal
+    scalar product is used. For custom definition of the method, a *nested_moments* can be passed.
+    For example: [ [1, x, y], [x*y, x**2, y**2], ... ] that groups all moments together that should be relaxed
+    at the same rate. Literature values of this list can be obtained through 
+    :func:`lbmpy.methods.creationfunctions.mrt_orthogonal_modes_literature`.
+    WMRT collision operators are reported to be numerically more stable and more accurate, 
+    whilst also having a lower computational cos :cite:`FAKHARI201722`
+    """
     CENTRAL_MOMENT = auto()
-    MRT_RAW = auto()
+    """
+    See :func:`lbmpy.methods.create_central_moment`
+    Creates moment based LB method where the collision takes place in the central moment space. By default, 
+    a raw-moment set is used where the bulk and the shear viscosity are separated. An original derivation can be 
+    found in :cite:`Geier2006`
+    """
     TRT_KBC_N1 = auto()
+    """
+    See :func:`lbmpy.methods.create_trt_kbc`
+    Particular two-relaxation rate method. This is not the entropic method yet, only the relaxation pattern. 
+    To get the entropic method also *entropic* needs to be set to `True`. 
+    There are four KBC methods available in lbmpy. The naming is according to :cite:`karlin2015entropic`
+    """
     TRT_KBC_N2 = auto()
+    """
+    See :func:`lbmpy.methods.create_trt_kbc`
+    Particular two-relaxation rate method. This is not the entropic method yet, only the relaxation pattern. 
+    To get the entropic method also *entropic* needs to be set to `True`. 
+    There are four KBC methods available in lbmpy. The naming is according to :cite:`karlin2015entropic`
+    """
     TRT_KBC_N3 = auto()
+    """
+    See :func:`lbmpy.methods.create_trt_kbc`
+    Particular two-relaxation rate method. This is not the entropic method yet, only the relaxation pattern. 
+    To get the entropic method also *entropic* needs to be set to `True`. 
+    There are four KBC methods available in lbmpy. The naming is according to :cite:`karlin2015entropic`
+    """
     TRT_KBC_N4 = auto()
+    """
+    See :func:`lbmpy.methods.create_trt_kbc`
+    Particular two-relaxation rate method. This is not the entropic method yet, only the relaxation pattern. 
+    To get the entropic method also *entropic* needs to be set to `True`. 
+    There are four KBC methods available in lbmpy. The naming is according to :cite:`karlin2015entropic`
+    """
     ENTROPIC_SRT = auto()
+    """
+    See :func:`lbmpy.methods.create_srt_entropic`,
+    An entropic version of the isothermal lattice Boltzmann method with the simplicity and 
+    computational efficiency of the standard lattice Boltzmann model. For details see :cite:`Ansumali2003`
+    """
     CUMULANT = auto()
+    """
+    See :func:`lbmpy.methods.create_with_default_polynomial_cumulants`
+    Cumulant-based LB method which relaxes groups of polynomial cumulants chosen to optimize rotational invariance.
+    For details on the method see :cite:`geier2015`
+    """
     MONOMIAL_CUMULANT = auto()
+    """
+    See :func:`lbmpy.methods.create_with_monomial_cumulants`
+    Cumulant-based LB method which relaxes monomial cumulants.
+    For details on the method see :cite:`geier2015` and :cite:`Coreixas2019`
+    """
 
 
 class ForceModel(Enum):
+    """
+    The ForceModel enumeration defines which force model is used to introduce forcing terms in the collision operator
+    of the lattice Boltzmann method. A short summary of the theory behind is shown in `lbmpy.forcemodels`.
+    More precise definitions are given in Chapter 6 and 10 of :cite:`lbm_book`
+    """
     SIMPLE = auto()
+    """
+    See :class:`lbmpy.forcemodels.Simple`
+    """
     LUO = auto()
+    """
+    See :class:`lbmpy.forcemodels.Luo`
+    """
     GUO = auto()
+    """
+    See :class:`lbmpy.forcemodels.Guo`
+    """
     BUICK = auto()
+    """
+    See :class:`lbmpy.forcemodels.Buick`
+    """
     SILVA = auto()
+    """
+    See :class:`lbmpy.forcemodels.Buick`
+    """
     EDM = auto()
+    """
+    See :class:`lbmpy.forcemodels.EDM`
+    """
     KUPERSHTOKH = auto()
+    """
+    See :class:`lbmpy.forcemodels.EDM`
+    """
     CUMULANT = auto()
+    """
+    See :class:`lbmpy.methods.centeredcumulant.CenteredCumulantForceModel`
+    """
     HE = auto()
+    """
+    See :class:`lbmpy.forcemodels.He`
+    """
     SHANCHEN = auto()
+    """
+    See :class:`lbmpy.forcemodels.ShanChen`
+    """

lbmpy/fieldaccess.py

@@ -7,6 +7,9 @@ from pystencils import Field
 from pystencils.astnodes import LoopOverCoordinate
 from pystencils.stencil import inverse_direction
 
+from lbmpy.enums import Stencil
+from lbmpy.stencils import LBStencil
+
 __all__ = ['PdfFieldAccessor', 'CollideOnlyInplaceAccessor', 'StreamPullTwoFieldsAccessor',
            'AAEvenTimeStepAccessor', 'AAOddTimeStepAccessor',
            'PeriodicTwoFieldsAccessor', 'StreamPushTwoFieldsAccessor',
@@ -51,11 +54,11 @@ class CollideOnlyInplaceAccessor(PdfFieldAccessor):
 
     @staticmethod
     def read(field, stencil):
-        return [field(i) for i in range(len(stencil))]
+        return [field(i) for i in range(stencil.Q)]
 
     @staticmethod
     def write(field, stencil):
-        return [field(i) for i in range(len(stencil))]
+        return [field(i) for i in range(stencil.Q)]
 
 
 class StreamPullTwoFieldsAccessor(PdfFieldAccessor):
@@ -67,7 +70,7 @@ class StreamPullTwoFieldsAccessor(PdfFieldAccessor):
 
     @staticmethod
     def write(field, stencil):
-        return [field(i) for i in range(len(stencil))]
+        return [field(i) for i in range(stencil.Q)]
 
 
 class StreamPushTwoFieldsAccessor(PdfFieldAccessor):
@@ -75,7 +78,7 @@ class StreamPushTwoFieldsAccessor(PdfFieldAccessor):
 
     @staticmethod
     def read(field, stencil):
-        return [field(i) for i in range(len(stencil))]
+        return [field(i) for i in range(stencil.Q)]
 
     @staticmethod
     def write(field, stencil):
@@ -125,7 +128,7 @@ class PeriodicTwoFieldsAccessor(PdfFieldAccessor):
 
     @staticmethod
     def write(field, stencil):
-        return [field(i) for i in range(len(stencil))]
+        return [field(i) for i in range(stencil.Q)]
 
 
 class AAEvenTimeStepAccessor(PdfFieldAccessor):
@@ -133,7 +136,7 @@ class AAEvenTimeStepAccessor(PdfFieldAccessor):
 
     @staticmethod
     def read(field, stencil):
-        return [field(i) for i in range(len(stencil))]
+        return [field(i) for i in range(stencil.Q)]
 
     @staticmethod
     def write(field, stencil):
@@ -214,15 +217,18 @@ def visualize_field_mapping(axes, stencil, field_mapping, inverted=False, color=
     grid.draw(axes)
 
 
-def visualize_pdf_field_accessor(pdf_field_accessor, title=True, read_plot_params={}, write_plot_params={},
+def visualize_pdf_field_accessor(pdf_field_accessor, title=True, read_plot_params=None, write_plot_params=None,
                                  figure=None):
-    from lbmpy.stencils import get_stencil
 
+    if write_plot_params is None:
+        write_plot_params = {}
+    if read_plot_params is None:
+        read_plot_params = {}
     if figure is None:
         import matplotlib.pyplot as plt
         figure = plt.gcf()
 
-    stencil = get_stencil('D2Q9')
+    stencil = LBStencil(Stencil.D2Q9)
 
     figure.patch.set_facecolor('white')
 

lbmpy/forcemodels.py

@@ -17,16 +17,16 @@ Force models add a term :math:`C_F` to the collision equation:
 
 .. math ::
 
-    f(\pmb{x} + c_q \Delta t, t + \Delta t) - f(\pmb{x},t) = \Omega(f, f^{(\mathrm{eq})})
+    f(\mathbf{x} + c_q \Delta t, t + \Delta t) - f(\mathbf{x},t) = \Omega(f, f^{(\mathrm{eq})})
                                                             + \underbrace{F_q}_{\mbox{forcing term}}
 
 The form of this term depends on the concrete force model: the first moment of this forcing term is equal
-to the acceleration :math:`\pmb{a}` for all force models. Here :math:`\mathbf{F}` is the D dimensional force vector,
+to the acceleration :math:`\mathbf{a}` for all force models. Here :math:`\mathbf{F}` is the D dimensional force vector,
 which defines the force for each spatial dircetion.
 
 .. math ::
 
-    \sum_q \pmb{c}_q \mathbf{F} = \pmb{a}
+    \sum_q \mathbf{c}_q \mathbf{F} = \mathbf{a}
 
 
 The second order moment is different for the forcing models - if it is zero the model is suited for
@@ -57,7 +57,7 @@ For all force models the computation of the macroscopic velocity has to be adapt
     
     .. math ::
     
-        \pmb{u} &= \sum_q \pmb{c}_q f_q + S_{\mathrm{macro}}
+        \mathbf{u} &= \sum_q \mathbf{c}_q f_q + S_{\mathrm{macro}}
         
         S_{\mathrm{macro}} &= \frac{\Delta t}{2 \cdot \rho} \sum_q F_q
 
@@ -296,7 +296,7 @@ class He(AbstractForceModel):
                 F_x m^{\mathrm{eq}}_{\alpha+1,\beta,\gamma} 
                 + F_y m^{\mathrm{eq}}_{\alpha,\beta+1,\gamma} 
                 + F_z m^{\mathrm{eq}}_{\alpha,\beta,\gamma+1} 
-                - m^{eq}_{\alpha\beta\gamma} ( \mathbf{F} \cdot \vec{u} ) 
+                - m^{eq}_{\alpha\beta\gamma} ( \mathbf{F} \cdot \mathbf{u} )
             \right)
     """
 

lbmpy/lbstep.py

 from types import MappingProxyType
+from dataclasses import replace
 
 import numpy as np
 
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
-from lbmpy.creationfunctions import (
-    create_lb_function, update_with_default_parameters)
+from lbmpy.creationfunctions import (create_lb_function, update_with_default_parameters)
 from lbmpy.enums import Stencil
 from lbmpy.macroscopic_value_kernels import (
     create_advanced_velocity_setter_collision_rule, pdf_initialization_assignments)
@@ -18,14 +18,16 @@ from pystencils.timeloop import TimeLoop
 class LatticeBoltzmannStep:
 
     def __init__(self, domain_size=None, lbm_kernel=None, periodicity=False,
-                 kernel_params=MappingProxyType({}), data_handling=None, name="lbm", optimization={},
+                 kernel_params=MappingProxyType({}), data_handling=None, name="lbm", optimization=None,
                  velocity_data_name=None, density_data_name=None, density_data_index=None,
                  compute_velocity_in_every_step=False, compute_density_in_every_step=False,
                  velocity_input_array_name=None, time_step_order='stream_collide', flag_interface=None,
-                 alignment_if_vectorized=64, fixed_loop_sizes=True, fixed_relaxation_rates=True,
+                 alignment_if_vectorized=64, fixed_loop_sizes=True,
                  timeloop_creation_function=TimeLoop,
                  lbm_config=None, lbm_optimisation=None, config=None, **method_parameters):
 
+        if optimization is None:
+            optimization = {}
         self._timeloop_creation_function = timeloop_creation_function
 
         # --- Parameter normalization  ---
@@ -55,7 +57,9 @@ class LatticeBoltzmannStep:
                                                                               lbm_config, lbm_optimisation, config)
 
         # the parallel datahandling understands only numpy datatypes. Strings lead to an error.
-        field_dtype = np.float64 if config.data_type == 'double' else np.float32
+        field_dtype = np.float64
+        if config.data_type == 'float' or config.data_type == 'float32':
+            field_dtype = np.float32
 
         if lbm_kernel:
             q = lbm_kernel.method.stencil.Q
@@ -99,10 +103,12 @@ class LatticeBoltzmannStep:
                 density_field = density_field(density_data_index)
             lbm_config.output['density'] = density_field
         if velocity_input_array_name is not None:
-            lbm_config.velocity_input = self._data_handling.fields[velocity_input_array_name]
+            lbm_config = replace(lbm_config, velocity_input=self._data_handling.fields[velocity_input_array_name])
         if isinstance(lbm_config.omega_output_field, str):
-            lbm_config.omega_output_field = data_handling.add_array(lbm_config.omega_output_field,
-                                                                    dtype=field_dtype, alignment=alignment)
+            lbm_config = replace(lbm_config, omega_output_field=data_handling.add_array(lbm_config.omega_output_field,
+                                                                                        dtype=field_dtype,
+                                                                                        alignment=alignment,
+                                                                                        values_per_cell=1))
 
         self.kernel_params = kernel_params.copy()
 
@@ -110,9 +116,10 @@ class LatticeBoltzmannStep:
         if lbm_kernel is None:
 
             if fixed_loop_sizes:
-                lbm_optimisation.symbolic_field = data_handling.fields[self._pdf_arr_name]
-            lbm_config.field_name = self._pdf_arr_name
-            lbm_config.temporary_field_name = self._tmp_arr_name
+                lbm_optimisation = replace(lbm_optimisation, symbolic_field=data_handling.fields[self._pdf_arr_name])
+            lbm_config = replace(lbm_config, field_name=self._pdf_arr_name)
+            lbm_config = replace(lbm_config, temporary_field_name=self._tmp_arr_name)
+
             if time_step_order == 'stream_collide':
                 self._lbmKernels = [create_lb_function(lbm_config=lbm_config,
                                                        lbm_optimisation=lbm_optimisation,

lbmpy/maxwellian_equilibrium.py

@@ -7,7 +7,6 @@ Additionally functions are provided to compute moments and cumulants of these di
 import warnings
 
 import sympy as sp
-from sympy import Rational as R
 
 from pystencils.cache import disk_cache
 from pystencils.sympyextensions import remove_higher_order_terms
@@ -17,38 +16,40 @@ from lbmpy.continuous_distribution_measures import continuous_moment, continuous
 
 
 def get_weights(stencil, c_s_sq=sp.Rational(1, 3)):
-    q = len(stencil)
-
     if c_s_sq != sp.Rational(1, 3) and c_s_sq != sp.Symbol("c_s") ** 2:
         warnings.warn("Weights of discrete equilibrium are only valid if c_s^2 = 1/3")
 
     def weight_for_direction(direction):
         abs_sum = sum([abs(d) for d in direction])
-        return get_weights.weights[q][abs_sum]
+        return get_weights.weights[stencil.Q][abs_sum]
     return [weight_for_direction(d) for d in stencil]
 
 
 get_weights.weights = {
     9: {
-        0: R(4, 9),
-        1: R(1, 9),
-        2: R(1, 36),
+        0: sp.Rational(4, 9),
+        1: sp.Rational(1, 9),
+        2: sp.Rational(1, 36),
+    },
+    7: {
+        0: sp.simplify(0.0),
+        1: sp.Rational(1, 6),
     },
     15: {
-        0: R(2, 9),
-        1: R(1, 9),
-        3: R(1, 72),
+        0: sp.Rational(2, 9),
+        1: sp.Rational(1, 9),
+        3: sp.Rational(1, 72),
     },
     19: {
-        0: R(1, 3),
-        1: R(1, 18),
-        2: R(1, 36),
+        0: sp.Rational(1, 3),
+        1: sp.Rational(1, 18),
+        2: sp.Rational(1, 36),
     },
     27: {
-        0: R(8, 27),
-        1: R(2, 27),
-        2: R(1, 54),
-        3: R(1, 216),
+        0: sp.Rational(8, 27),
+        1: sp.Rational(2, 27),
+        2: sp.Rational(1, 54),
+        3: sp.Rational(1, 216),
     }
 }
 
@@ -68,10 +69,9 @@ def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=sp.symbols(
         compressible: compressibility
     """
     weights = get_weights(stencil, c_s_sq)
-    assert len(stencil) == len(weights)
+    assert stencil.Q == len(weights)
 
-    dim = len(stencil[0])
-    u = u[:dim]
+    u = u[:stencil.D]
 
     rho_outside = rho if compressible else sp.Rational(1, 1)
     rho_inside = rho if not compressible else sp.Rational(1, 1)
@@ -113,14 +113,12 @@ def generate_equilibrium_by_matching_moments(stencil, moments, rho=sp.Symbol("rh
     see :func:`get_equilibrium_values_of_maxwell_boltzmann_function`
     """
     from lbmpy.moments import moment_matrix
-    dim = len(stencil[0])
-    Q = len(stencil)
-    assert len(moments) == Q, "Moment count(%d) does not match stencil size(%d)" % (len(moments), Q)
-    continuous_moments_vector = get_equilibrium_values_of_maxwell_boltzmann_function(moments, dim, rho, u, c_s_sq,
+    assert len(moments) == stencil.Q, f"Moment count({len(moments)}) does not match stencil size({stencil.Q})"
+    continuous_moments_vector = get_equilibrium_values_of_maxwell_boltzmann_function(moments, stencil.D, rho, u, c_s_sq,
                                                                                      order, space="moment")
     continuous_moments_vector = sp.Matrix(continuous_moments_vector)
     M = moment_matrix(moments, stencil)
-    assert M.rank() == Q, "Rank of moment matrix (%d) does not match stencil size (%d)" % (M.rank(), Q)
+    assert M.rank() == stencil.Q, f"Rank of moment matrix ({M.rank()}) does not match stencil size ({stencil.Q})"
     return M.inv() * continuous_moments_vector
 
 

lbmpy/methods/__init__.py

 from lbmpy.methods.creationfunctions import (
     create_mrt_orthogonal, create_mrt_raw, create_central_moment, create_srt, create_trt, create_trt_kbc,
     create_trt_with_magic_number, create_with_continuous_maxwellian_eq_moments,
-    create_with_discrete_maxwellian_eq_moments, mrt_orthogonal_modes_literature,
+    create_with_discrete_maxwellian_eq_moments,
     create_centered_cumulant_model, create_with_default_polynomial_cumulants,
     create_with_polynomial_cumulants, create_with_monomial_cumulants)
 
-from lbmpy.methods.abstractlbmethod import AbstractLbMethod, RelaxationInfo
+from lbmpy.methods.default_moment_sets import mrt_orthogonal_modes_literature, cascaded_moment_sets_literature
+
+from lbmpy.methods.abstractlbmethod import LbmCollisionRule, AbstractLbMethod, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
 
 from .conservedquantitycomputation import DensityVelocityComputation
 
-__all__ = ['RelaxationInfo', 'AbstractLbMethod',
+__all__ = ['RelaxationInfo', 'AbstractLbMethod', 'LbmCollisionRule',
            'AbstractConservedQuantityComputation', 'DensityVelocityComputation',
            'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
            'create_mrt_orthogonal', 'create_mrt_raw', 'create_central_moment',
            'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments',
-           'mrt_orthogonal_modes_literature', 'create_centered_cumulant_model',
-           'create_with_default_polynomial_cumulants', 'create_with_polynomial_cumulants',
-           'create_with_monomial_cumulants']
+           'mrt_orthogonal_modes_literature', 'cascaded_moment_sets_literature',
+           'create_centered_cumulant_model', 'create_with_default_polynomial_cumulants',
+           'create_with_polynomial_cumulants', 'create_with_monomial_cumulants']

lbmpy/methods/abstractlbmethod.py

@@ -2,6 +2,7 @@ import abc
 from collections import namedtuple
 
 import sympy as sp
+from sympy.core.numbers import Zero
 
 from pystencils import Assignment, AssignmentCollection
 
@@ -9,6 +10,9 @@ RelaxationInfo = namedtuple('RelaxationInfo', ['equilibrium_value', 'relaxation_
 
 
 class LbmCollisionRule(AssignmentCollection):
+    """
+    A pystencils AssignmentCollection that additionally holds an `AbstractLbMethod`
+    """
     def __init__(self, lb_method, *args, **kwargs):
         super(LbmCollisionRule, self).__init__(*args, **kwargs)
         self.method = lb_method
@@ -49,6 +53,7 @@ class AbstractLbMethod(abc.ABC):
         """Returns a qxq diagonal matrix which contains the relaxation rate for each moment on the diagonal"""
         d = sp.zeros(len(self.relaxation_rates))
         for i in range(0, len(self.relaxation_rates)):
+            # note that 0.0 is converted to sp.Zero here. It is not possible to prevent this.
             d[i, i] = self.relaxation_rates[i]
         return d
 
@@ -101,6 +106,9 @@ class AbstractLbMethod(abc.ABC):
         for relaxation_rate in rr:
             if relaxation_rate not in unique_relaxation_rates:
                 relaxation_rate = sp.sympify(relaxation_rate)
+                # special treatment for zero, sp.Zero would be an integer ..
+                if isinstance(relaxation_rate, Zero):
+                    relaxation_rate = 0.0
                 if not isinstance(relaxation_rate