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

@@ -6,21 +6,23 @@ import sympy as sp
 
 from lbmpy.moments import polynomial_to_exponent_representation
 from pystencils.cache import disk_cache, memorycache
-from pystencils.sympyextensions import complete_the_squares_in_exp
+from pystencils.sympyextensions import complete_the_squares_in_exp, scalar_product
 
 
 @memorycache()
-def moment_generating_function(generating_function, symbols, symbols_in_result):
+def moment_generating_function(generating_function, symbols, symbols_in_result, velocity=None):
     r"""
     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.
 
     Returns:
         transformation result F: an expression that depends now on symbols_in_result
@@ -55,9 +57,27 @@ def moment_generating_function(generating_function, symbols, symbols_in_result):
     return sp.simplify(result)
 
 
-def cumulant_generating_function(func, symbols, symbols_in_result):
+def central_moment_generating_function(func, symbols, symbols_in_result, velocity=sp.symbols("u_:3")):
+    r"""
+    Computes central moment generating func, which is defined as:
+
+    .. math ::
+        K( \mathbf{\Xi} ) = \exp ( - \mathbf{\Xi} \cdot \mathbf{u} ) M( \mathbf{\Xi} ).
+
+    For parameter description see :func:`moment_generating_function`.
     """
-    Computes cumulant generating func, which is the logarithm of the moment generating func.
+    argument = - scalar_product(symbols_in_result, velocity)
+
+    return sp.exp(argument) * moment_generating_function(func, symbols, symbols_in_result)
+
+
+def cumulant_generating_function(func, symbols, symbols_in_result, velocity=None):
+    r"""
+    Computes cumulant generating func, which is the logarithm of the moment generating func:
+
+    .. math ::
+        C(\mathbf{\Xi}) = \log M(\mathbf{\Xi})
+
     For parameter description see :func:`moment_generating_function`.
     """
     return sp.ln(moment_generating_function(func, symbols, symbols_in_result))
@@ -93,16 +113,16 @@ def multi_differentiation(generating_function, index, symbols):
 
 
 @memorycache(maxsize=512)
-def __continuous_moment_or_cumulant(func, moment, symbols, generating_function):
+def __continuous_moment_or_cumulant(func, moment, symbols, generating_function, velocity=sp.symbols("u_:3")):
     if type(moment) is tuple and not symbols:
         symbols = sp.symbols("xvar yvar zvar")
 
     dim = len(moment) if type(moment) is tuple else len(symbols)
 
     # not using sp.Dummy here - since it prohibits caching
-    t = tuple([sp.Symbol("tmpvar_%d" % i, ) for i in range(dim)])
+    t = sp.symbols(f"tmpvar_:{dim}")
     symbols = symbols[:dim]
-    generating_function = generating_function(func, symbols, t)
+    generating_function = generating_function(func, symbols, t, velocity=velocity)
 
     if type(moment) is tuple:
         return multi_differentiation(generating_function, moment, t)
@@ -128,6 +148,18 @@ def continuous_moment(func, moment, symbols=None):
     return __continuous_moment_or_cumulant(func, moment, symbols, moment_generating_function)
 
 
+def continuous_central_moment(func, moment, symbols=None, velocity=sp.symbols("u_:3")):
+    """Computes central moment of given function.
+
+    Args:
+        func: function to compute moments of
+        moment: tuple or polynomial describing the moment
+        symbols: if moment is given as polynomial, pass the moment symbols, i.e. the dof of the polynomial
+    """
+    return __continuous_moment_or_cumulant(func, moment, symbols, central_moment_generating_function,
+                                           velocity=velocity)
+
+
 def continuous_cumulant(func, moment, symbols=None):
     """Computes cumulant of continuous 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,10 +121,10 @@ 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 = tuple([sp.Symbol("Xi_%d" % (i,)) for i in range(dim)])
+    wave_numbers = sp.symbols(f"Xi_:{dim}")
 
     generating_function = __get_discrete_cumulant_generating_function(func, stencil, wave_numbers)
     if type(cumulant) is tuple:
@@ -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/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/methods/__init__.py

 from lbmpy.methods.creationfunctions import (
-    create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc,
+    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_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,13 +2,17 @@ import abc
 from collections import namedtuple
 
 import sympy as sp
+from sympy.core.numbers import Zero
 
-from pystencils import AssignmentCollection
+from pystencils import Assignment, AssignmentCollection
 
 RelaxationInfo = namedtuple('RelaxationInfo', ['equilibrium_value', 'relaxation_rate'])
 
 
 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
@@ -32,12 +36,41 @@ class AbstractLbMethod(abc.ABC):
     @property
     def pre_collision_pdf_symbols(self):
         """Tuple of symbols representing the pdf values before collision"""
-        return sp.symbols("f_:%d" % (len(self.stencil),))
+        return sp.symbols(f"f_:{len(self.stencil)}")
 
     @property
     def post_collision_pdf_symbols(self):
         """Tuple of symbols representing the pdf values after collision"""
-        return sp.symbols("d_:%d" % (len(self.stencil),))
+        return sp.symbols(f"d_:{len(self.stencil)}")
+
+    @property
+    @abc.abstractmethod
+    def relaxation_rates(self):
+        """Tuple containing the relaxation rates of each moment"""
+
+    @property
+    def relaxation_matrix(self):
+        """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
+
+    @property
+    def symbolic_relaxation_matrix(self):
+        """Returns a qxq diagonal matrix which contains the relaxation rate for each moment on the diagonal.
+           In contrast to the normal relaxation matrix all numeric values are replaced by sympy symbols"""
+        _, d = self._generate_symbolic_relaxation_matrix()
+        return d
+
+    @property
+    def subs_dict_relxation_rate(self):
+        """returns a dictonary which maps the replaced numerical relaxation rates to its original numerical value"""
+        result = dict()
+        for i in range(len(self.stencil)):
+            result[self.symbolic_relaxation_matrix[i, i]] = self.relaxation_matrix[i, i]
+        return result
 
     # ------------------------- Abstract Methods & Properties ----------------------------------------------------------
 
@@ -59,3 +92,29 @@ 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."""
+
+    # -------------------------------- Helper Functions ----------------------------------------------------------------
+
+    def _generate_symbolic_relaxation_matrix(self):
+        """
+        This function replaces the numbers in the relaxation matrix with symbols in this case, and returns also
+        the subexpressions, that assign the number to the newly introduced symbol
+        """
+        rr = [self.relaxation_matrix[i, i] for i in range(self.relaxation_matrix.rows)]
+        unique_relaxation_rates = set()
+        subexpressions = {}
+        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, sp.Symbol):
+                    rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
+                    subexpressions[relaxation_rate] = rt_symbol
+            unique_relaxation_rates.add(relaxation_rate)
+
+        new_rr = [subexpressions[sp.sympify(e)] if sp.sympify(e) in subexpressions else e for e in rr]
+        substitutions = [Assignment(e[1], e[0]) for e in subexpressions.items()]
+
+        return substitutions, sp.diag(*new_rr)

lbmpy/methods/centeredcumulant/__init__.py

 from .force_model import CenteredCumulantForceModel
 from .centeredcumulantmethod import CenteredCumulantBasedLbMethod
-from .centered_cumulants import get_default_polynomial_cumulants_for_stencil
 
-__all__ = ['CenteredCumulantForceModel', 'CenteredCumulantBasedLbMethod',
-           'get_default_polynomial_cumulants_for_stencil']
+__all__ = ['CenteredCumulantForceModel', 'CenteredCumulantBasedLbMethod']

lbmpy/methods/centeredcumulant/force_model.py

-from lbmpy.forcemodels import default_velocity_shift
+from lbmpy.forcemodels import AbstractForceModel, default_velocity_shift
 
 
 #   =========================== Centered Cumulant Force Model ==========================================================
 
 
-class CenteredCumulantForceModel:
+class CenteredCumulantForceModel(AbstractForceModel):
     """
     A force model to be used with the centered cumulant-based LB Method.
     It only shifts the macroscopic and equilibrium velocities and does not introduce a forcing term to the
@@ -16,14 +16,12 @@ class CenteredCumulantForceModel:
     """
 
     def __init__(self, force):
-        self._force = force
         self.override_momentum_relaxation_rate = 2
 
-    def __call__(self, lb_method, **kwargs):
-        raise Exception('This force model does not provide a forcing term.')
+        super(CenteredCumulantForceModel, self).__init__(force)
 
-    def macroscopic_velocity_shift(self, density):
-        return default_velocity_shift(density, self._force)
+    def __call__(self, lb_method):
+        raise Exception('This force model does not provide a forcing term.')
 
     def equilibrium_velocity_shift(self, density):
-        return default_velocity_shift(density, self._force)
+        return default_velocity_shift(density, self.symbolic_force_vector)

lbmpy/methods/centeredcumulant/galilean_correction.py

@@ -2,8 +2,7 @@ from pystencils.simp.assignment_collection import AssignmentCollection
 import sympy as sp
 from pystencils import Assignment
 
-from lbmpy.moments import MOMENT_SYMBOLS
-from lbmpy.methods.centeredcumulant.centered_cumulants import statistical_quantity_symbol
+from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol
 from lbmpy.methods.centeredcumulant.cumulant_transform import PRE_COLLISION_CUMULANT
 
 x, y, z = MOMENT_SYMBOLS