Rework force models and central moments transform

README.md

@@ -56,17 +56,26 @@ Documentation
 Read the docs [here](http://pycodegen.pages.i10git.cs.fau.de/lbmpy) and
 check out the Jupyter notebooks in `doc/notebooks`. 
 
-Authors
+Contributing
 -------
+To see how to open issues, submit bug reports, create feature requests or submit your additions to lbmpy please refer to
+[contribution documentation](https://i10git.cs.fau.de/pycodegen/pystencils/-/blob/master/CONTRIBUTING.md) of pystencils since lbmpy is heavily build on pystencils.
+
+
+Many thanks go to the [contributors](https://i10git.cs.fau.de/pycodegen/lbmpy/-/blob/master/AUTHORS.txt) of lbmpy.
 
-Many thanks go to the [contributors](AUTHORS.txt) of lbmpy.
 
 ### Please cite us
 
 If you use lbmpy in a publication, please cite the following articles:
 
 Overview:
-  - M. Bauer et al, lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods. Journal of Computational Science, 2021. https://doi.org/10.1016/j.jocs.2020.101269
+  - M. Bauer et al, lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods. Journal of Computational Science, 2021. https://doi.org/10.1016/j.jocs.2020.101269 ([Preprint](https://arxiv.org/abs/2001.11806))
 
 Multiphase:
    - M. Holzer et al, Highly efficient lattice Boltzmann multiphase simulations of immiscible fluids at high-density ratios on CPUs and GPUs through code generation. The International Journal of High Performance Computing Applications, 2021. https://doi.org/10.1177/10943420211016525
+
+
+### Further Reading
+
+- F. Hennig et al, Automatic Code Generation for the Cumulant Lattice Boltzmann Method. ICMMES, 2021. [Poster Link](https://www.researchgate.net/publication/353224406_Automatic_Code_Generation_for_the_Cumulant_Lattice_Boltzmann_Method)

doc/conf.py

@@ -21,6 +21,9 @@ extensions = [
     'sphinx_autodoc_typehints',
 ]
 
+# set mathjax v3 path according to https://www.sphinx-doc.org/en/master/usage/extensions/math.html
+mathjax_path="https://cdn.jsdelivr.net/npm/mathjax@2/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
+
 add_module_names = False
 templates_path = ['_templates']
 source_suffix = '.rst'

doc/sphinx/lbmpy.bib

@@ -6,7 +6,8 @@
   number={4},
   pages={043309},
   year={2015},
-  publisher={APS}
+  publisher={APS},
+  doi={10.1103/PhysRevE.92.043309},
 }
 
 @PHDTHESIS{luo1993lattice,
@@ -16,7 +17,7 @@
    school = {GEORGIA INSTITUTE OF TECHNOLOGY.},
      year = 1993,
    adsurl = {http://adsabs.harvard.edu/abs/1993PhDT.......233L},
-  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System},
 }
 
 @Article{guo2002discrete,
@@ -28,8 +29,24 @@
   number    = {4},
   pages     = {046308},
   publisher = {APS},
+  doi = {10.1103/PhysRevE.65.046308},
 }
 
+@article{HeForce,
+  title = {Discrete Boltzmann equation model for nonideal gases},
+  author = {He, Xiaoyi and Shan, Xiaowen and Doolen, Gary D.},
+  journal = {Physical Review E},
+  volume = {57},
+  issue = {1},
+  pages = {R13--R16},
+  numpages = {0},
+  year = {1998},
+  month = {1},
+  publisher = {APS},
+  doi = {10.1103/PhysRevE.57.R13}
+}
+
+
 @article{buick2000gravity,
   title={Gravity in a lattice Boltzmann model},
   author={Buick, JM and Greated, CA},
@@ -38,7 +55,8 @@
   number={5},
   pages={5307},
   year={2000},
-  publisher={APS}
+  publisher={APS},
+  doi = {10.1103/PhysRevE.61.5307},
 }
 
 @phdthesis{schiller2008thermal,
@@ -120,3 +138,35 @@ issn = {0309-1708},
 doi = {10.1016/j.advwatres.2018.02.005},
 author = {Fakhari, A. and Li, Y. and Bolster, D. and Christensen, K. T.},
 }
+
+@article{silva2010,
+title = {A study on the inclusion of body forces in the lattice Boltzmann BGK equation to recover steady-state hydrodynamics},
+journal = {Physica A: Statistical Mechanics and its Applications},
+volume = {390},
+number = {6},
+pages = {1085-1095},
+year = {2011},
+doi = {10.1016/j.physa.2010.11.037},
+author = {Silva, G. and Semiao, V.},
+keywords = {{Lattice Boltzmann} method, {BGK} collision operator, Steady-state flows, Body force driven flows}
+}
+
+@article{silva2020,
+  title = {Force methods for the two-relaxation-times lattice {Boltzmann}},
+  author = {Postma, B. and Silva, G.},
+  journal = {Phys. Rev. E},
+  volume = {102},
+  issue = {6},
+  year = {2020},
+  month = {12},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.102.063307},
+}
+
+@book{lbm_book,
+  doi = {10.1007/978-3-319-44649-3},
+  year = {2017},
+  publisher = {Springer International Publishing},
+  author = {Kr\"{u}ger, T. and Kusumaatmaja, H. and Kuzmin, A. and Shardt, O. and Silva, G. and Viggen, E. M.},
+  title = {The Lattice {Boltzmann} Method}
+}
\ No newline at end of file

doc/sphinx/methods.rst

@@ -68,8 +68,6 @@ Creation Functions
 Utility
 -------
 
-.. autofunction:: lbmpy.methods.centeredcumulant.get_default_polynomial_cumulants_for_stencil
-
 .. autoclass:: lbmpy.methods.centeredcumulant.CenteredCumulantForceModel
     :members:
 

doc/sphinx/moment_transforms.rst

@@ -2,31 +2,241 @@
 Moment Transforms (lbmpy.moment_transforms)
 *******************************************
 
+The ``moment_transforms`` module provides an ecosystem for transformation of quantities between
+lattice Boltzmann population space and other spaces of observable quantities. Currently, transforms
+to raw and central moment space are available.
+
+The common base class `lbmpy.moment_transforms.AbstractMomentTransform` defines the interface all transform classes share.
+This interface, together with the fundamental principles all transforms adhere to, is explained 
+in the following.
+
+Principles of Collision Space Transforms
+========================================
+
+Each class of this module implements a bijective map :math:`\mathcal{M}` from population space 
+to raw moment or central moment space, capable of transforming the particle distribution 
+function :math:`\mathbf{f}` to (almost) arbitrary user-defined moment sets.
+
+Monomial and Polynomial Moments
+-------------------------------
+
+We discriminate *monomial* and *polynomial* moments. The monomial moments are the canonical basis of their
+corresponding space. Polynomial moments are defined as linear combinations of this basis. Monomial moments are
+typically expressed by exponent tuples :math:`(\alpha, \beta, \gamma)`. The monomial raw moments are defined as
+
+.. math::
+    
+    m_{\alpha \beta \gamma} 
+        = \sum_{i = 0}^{q - 1} f_i c_{i,x}^{\alpha} c_{i,y}^{\beta} c_{i,z}^{\gamma}
+
+and the monomial central moments are given by
+
+.. math::
+    
+    \kappa_{\alpha \beta \gamma} 
+        = \sum_{i = 0}^{q - 1} 
+            f_i 
+            (c_{i,x} - u_x)^{\alpha} 
+            (c_{i,y} - u_y)^{\beta} 
+            (c_{i,z} - u_z)^{\gamma}.
+
+Polynomial moments are, in turn, expressed by
+polynomial expressions :math:`p` in the coordinate symbols :math:`x`, :math:`y` and :math:`z`.
+An exponent tuple :math:`(\alpha, \beta, \gamma)` corresponds directly 
+to the mixed monomial expression :math:`x^{\alpha} y^{\beta} z^{\gamma}`. Polynomial moments are then
+constructed as linear combinations of these monomials. For example, the polynomial
+
+.. math::
+
+    p(x,y,z) = x^2 + y^2 + z^2 + 1.
+
+defines both the polynomial raw moment
+
+.. math::
+
+    M = m_{200} + m_{020} + m_{002}
+
+and central moment
+
+.. math::
+
+    K = \kappa_{200} + \kappa_{020} + \kappa_{002}.
+
+
+Transformation Matrices
+-----------------------
+
+The collision space basis for an MRT LB method on a :math:`DdQq` lattice is spanned by a set of :math:`q`
+polynomial quantities, given by polynomials :math:`\left( p_i \right)_{i=0, \dots, q-1}`.
+Through the polynomials, a polynomialization matrix :math:`P` is defined such that
+
+.. math::
+
+    \mathbf{M} &= P \mathbf{m} \\
+    \mathbf{K} &= P \mathbf{\kappa}
+
+Both rules can also be written using matrix multiplication, such that
+
+.. math::
+    \mathbf{m} = M \mathbf{f} 
+    \qquad 
+    \mathbf{\kappa} = K \mathbf{f}.
+
+Further, there exists a mapping from raw to central moment space using (monomial or polynomial)
+shift matrices (see `set_up_shift_matrix`) such that
+
+.. math::
+    \mathbf{\kappa} = N \mathbf{m}
+    \qquad
+    \mathbf{K} = N^P \mathbf{M}.
+
+Working with the transformation matrices, those mappings can easily be inverted.
+This module provides classes that derive efficient implementations of these transformations.
+
+Moment Aliasing
+---------------
+
+For a collision space transform to be invertible, exactly :math:`q` independent polynomial quantities of
+the collision space must be chosen. In that case, since all transforms discussed here are linear, the space of
+populations is isomorphic to the chosen moment space. The important word here is 'independent'. Even if :math:`q`
+distinct moment polynomials are chosen, due to discrete lattice artifacts, they might not span the entire collision
+space if chosen wrongly. The discrete lattice structure gives rise to *moment aliasing* effects. The most simple such
+effect occurs in the monomial raw moments, where are all nonzero even and all odd exponents are essentially the same.
+For example, we have :math:`m_{400} = m_{200}` or :math:`m_{305} = m_{101}`. This rule holds on every direct-neighborhood
+stencil. Sparse stencils, like :math:`D3Q15`, may introduce additional aliases. On the :math:`D3Q15` stencil, due to its
+structure, we have also :math:`m_{112} = m_{110}` and even :math:`m_{202} = m_{220} = m_{022} = m_{222}`.
+
+Including aliases in a set of monomial moments will lead to a non-invertible transform and is hence not possible.
+In polynomials, however, aliases may occur without breaking the transform. Several established sets of polynomial
+moments, like in orthogonal raw moment space MRT methods, will comprise :math:`q` polynomials :math:`\mathbf{M}` that in turn are built
+out of :math:`r > q` monomial moments :math:`\mathbf{m}`. In that set of monomials, non-independent moments have to be included by definition.
+Since the full transformation matrix :math:`M^P = PM` is still invertible, this is not a problem in general. If, however, the
+two steps of the transform are computed separately, it becomes problematic, as the matrices :math:`M \in \mathbb{R}^{r \times q}`
+and :math:`P \in \mathbb{R}^{q \times r}` are not invertible on their own. 
+
+But there is a remedy. By eliminating aliases from the moment polynomials, they can be reduced to a new set of polynomials based
+on a non-aliased reduced vector of monomial moments :math:`\tilde{\mathbf{m}} \in \mathbb{R}^{q}`, expressed through the reduced
+matrix :math:`\tilde{P}`.
+
+
+Interfaces and Usage
+====================
+
+Construction
+------------
+
+All moment transform classes expect either a sequence of exponent tuples or a sequence of polynomials that define
+the set of quantities spanning the destination space. If polynomials are given, the exponent tuples of the underlying
+set of monomials are extracted automatically. Depending on the implementation, moment aliases may be eliminated in the
+process, altering both the polynomials and the set of monomials.
+
+
+Forward and Backward Transform
+------------------------------
+
+The core functionality of the transform classes is given by the methods ``forward_transform`` and ``backward_transform``.
+They return assignment collections containing the equations to compute moments from populations, and vice versa.
+
+Symbols Used
+^^^^^^^^^^^^
+
+Unless otherwise specified by the user, monomial and polynomial quantities occur in the transformation equations according
+to the naming conventions listed in the following table:
+
++--------------------------------+---------------------------------------------+----------------------+
+|                                |              Monomial                       |    Polynomial        |
++--------------------------------+---------------------------------------------+----------------------+
+| Pre-Collision Raw Moment       | :math:`m_{\alpha \beta \gamma}`             | :math:`M_{i}`        |
++--------------------------------+---------------------------------------------+----------------------+
+| Post-Collision Raw Moment      | :math:`m_{post, \alpha \beta \gamma}`       | :math:`M_{post, i}`  |
++--------------------------------+---------------------------------------------+----------------------+
+| Pre-Collision Central Moment   | :math:`\kappa_{\alpha \beta \gamma}`        | :math:`K_{i}`        |
++--------------------------------+---------------------------------------------+----------------------+
+| Post-Collision Central Moment  | :math:`\kappa_{post, \alpha \beta \gamma}`  | :math:`K_{post, i}`  |
++--------------------------------+---------------------------------------------+----------------------+
+
+These symbols are also exposed by the member properties `pre_collision_symbols`, `post_collision_symbols`, 
+`pre_collision_monomial_symbols` and `post_collision_monomial_symbols`.
+
+Forward Transform
+^^^^^^^^^^^^^^^^^
+
+Implementations of the `lbmpy.moment_transforms.AbstractMomentTransform.forward_transform` method 
+derive equations for the transform from population to moment space, to compute pre-collision moments
+from pre-collision populations. The returned ``AssignmentCollection`` has the `pre_collision_symbols` 
+as left-hand sides of its main assignments, computed from the given ``pdf_symbols`` and, potentially,
+the macroscopic density and velocity. Depending on the implementation, the `pre_collision_monomial_symbols`
+may be part of the assignment collection in the form of subexpressions. 
+
+Backward Transform
+^^^^^^^^^^^^^^^^^^
+
+Implementations of `lbmpy.moment_transforms.AbstractMomentTransform.backward_transform` contain the post-collision
+polynomial quantities as free symbols of their equation right-hand sides, and compute the post-collision populations
+from those. The PDF symbols are the left-hand sides of the main assignments.
+
+
+Absorption of Conserved Quantity Equations
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Transformations from the population space to any space of observable quantities may *absorb* the equations 
+defining the macroscopic quantities entering the equilibrium (typically the density :math:`\rho` and the
+velocity :math:`\mathbf{u}`). This means that the ``forward_transform`` will possibly rewrite the
+assignments given in the constructor argument ``conserved_quantity_equations`` to reduce
+the total operation count. For example, in the transformation step from populations to
+raw moments (see `lbmpy.moment_transforms.PdfsToMomentsByChimeraTransform`), :math:`\rho` can be aliased as the zeroth-order moment
+:math:`m_{000}`. Assignments to the conserved quantities will then be part of the AssignmentCollection
+returned by ``forward_transform`` and need not be added to the collision rule separately. 
+
+Simplification
+^^^^^^^^^^^^^^
+
+Both ``forward_transform`` and ``backward_transform`` expect a keyword argument ``simplification``
+which can be used to direct simplification steps applied during the derivation of the transformation
+equations. Possible values are:
+
+- `False` or ``'none'``: No simplification is to be applied
+- `True` or ``'default'``: A default simplification strategy specific to the implementation is applied.
+    The actual simplification steps depend strongly on the nature of the equations. They are defined by
+    the implementation. It is the responsibility of the implementation to select the most effective
+    simplification strategy.
+- ``'default_with_cse'``: Same as ``'default'``, but with an additional pass of common subexpression elimination.
+
+
+Working With Monomials
+^^^^^^^^^^^^^^^^^^^^^^
+
+In certain situations, we want the ``forward_transform`` to yield equations for the monomial symbols :math:`m_{\alpha \beta \gamma}`
+and :math:`\kappa_{\alpha \beta \gamma}` only, and the ``backward_transform`` to return equations that also expect these symbols as input. 
+In this case, it is not sufficient to pass a set of monomials or exponent tuples to the constructor, as those are still treated as 
+polynomials internally. Instead, both transform methods expose keyword arguments ``return_monomials`` and ``start_from_monomials``, respectively.
+If set to true, equations in the monomial moments are returned. They are best used *only* together with the ``exponent_tuples`` constructor argument
+to have full control over the monomials. If polynomials are passed to the constructor, the behaviour of these flags is generally not well-defined,
+especially in the presence of aliases.
+
+
+The Transform Classes
+=====================
+
 Abstract Base Class
-===================
+-------------------
 
 .. autoclass:: lbmpy.moment_transforms.AbstractMomentTransform
     :members:
 
 
 Moment Space Transforms
-===========================
-
-By Matrix-Vector Multiplication
--------------------------------
+-----------------------
 
 .. autoclass:: lbmpy.moment_transforms.PdfsToMomentsByMatrixTransform
     :members:
 
-By Chimera-Transform
---------------------
-
 .. autoclass:: lbmpy.moment_transforms.PdfsToMomentsByChimeraTransform
     :members:
 
 
 Central Moment Space Transforms
-===============================
+-------------------------------
 
 .. autoclass:: lbmpy.moment_transforms.PdfsToCentralMomentsByMatrix
     :members:
@@ -38,7 +248,7 @@ Central Moment Space Transforms
     :members:
 
 Cumulant Space Transforms
-=========================
+-------------------------
 
 .. autoclass:: lbmpy.methods.centeredcumulant.cumulant_transform.CentralMomentsToCumulantsByGeneratingFunc
     :members:

lbmpy/boundaries/boundaryconditions.py

+import abc
 from warnings import warn
 
 from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
@@ -13,7 +14,7 @@ from pystencils.stencil import offset_to_direction_string, direction_string_to_o
 import sympy as sp
 
 
-class LbBoundary:
+class LbBoundary(abc.ABC):
     """Base class that all boundaries should derive from.
 
     Args:
@@ -307,6 +308,7 @@ class UBB(LbBoundary):
         if self._adaptVelocityToForce:
             cqc = lb_method.conserved_quantity_computation
             shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
+            shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
             velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.first_order_moment_symbols).main_assignments]
 
         c_s_sq = sp.Rational(1, 3)

lbmpy/creationfunctions.py

@@ -46,8 +46,9 @@ General:
   compressible.
 - ``equilibrium_order=2``: order in velocity, at which the equilibrium moment/cumulant approximation is
   truncated. Order 2 is sufficient to approximate Navier-Stokes
-- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``, ``'buick'``, ``'schiller'``, 
-  ``'cumulant'`` or an instance of a class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. 
+- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``/``'schiller'``,
+  ``'buick'``/``'silva'``, ``'edm'``/``'kupershtokh'``, ``'he'``, ``'shanchen'``, ``'cumulant'``,
+  or an instance of a class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. 
   For details, see :mod:`lbmpy.forcemodels`
 - ``force=(0,0,0)``: either constant force or a symbolic expression depending on field value
 - ``maxwellian_moments=True``: way to compute equilibrium moments/cumulants, if False the standard
@@ -408,7 +409,7 @@ def create_lb_method(**params):
 
     no_force_model = 'force_model' not in params or params['force_model'] == 'none' or params['force_model'] is None
     if not force_is_zero and no_force_model:
-        params['force_model'] = 'cumulant' if method_name.lower().endswith('cumulant') else 'schiller'
+        params['force_model'] = 'cumulant' if method_name.lower().endswith('cumulant') else 'guo'
 
     if 'force_model' in params:
         force_model = force_model_from_string(params['force_model'], params['force'][:dim])
@@ -455,12 +456,15 @@ def create_lb_method(**params):
             except IndexError:
                 raise ValueError("Too few relaxation rates specified")
             return res
+
         weighted = params['weighted'] if 'weighted' in params else True
         nested_moments = params['nested_moments'] if 'nested_moments' in params else None
         method = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, weighted=weighted,
                                        nested_moments=nested_moments, **common_params)
     elif method_name.lower() == 'central_moment':
-        method = create_central_moment(stencil_entries, relaxation_rates, **common_params)
+        nested_moments = params['nested_moments'] if 'nested_moments' in params else None
+        method = create_central_moment(stencil_entries, relaxation_rates,
+                                       nested_moments=nested_moments, **common_params)
     elif method_name.lower() == 'mrt_raw':
         method = create_mrt_raw(stencil_entries, relaxation_rates, **common_params)
     elif method_name.lower().startswith('trt-kbc-n'):
@@ -519,6 +523,7 @@ def create_lb_method_from_existing(method, modification_function):
                             zip(method.moments, method.moment_equilibrium_values, method.relaxation_rates))
         return create_generic_mrt(method.stencil, relaxation_table, compressible, method.force_model)
 
+
 # ----------------------------------------------------------------------------------------------------------------------
 
 
@@ -532,11 +537,14 @@ def force_model_from_string(force_model_name, force_values):
         'simple': forcemodels.Simple,
         'luo': forcemodels.Luo,
         'guo': forcemodels.Guo,
+        'schiller': forcemodels.Guo,
         'buick': forcemodels.Buick,
         'silva': forcemodels.Buick,
         'edm': forcemodels.EDM,
-        'schiller': forcemodels.Schiller,
+        'kupershtokh': forcemodels.EDM,
         'cumulant': cumulant_force_model.CenteredCumulantForceModel,
+        'he': forcemodels.He,
+        'shanchen': forcemodels.ShanChen
     }
     if force_model_name.lower() not in force_model_dict:
         raise ValueError("Unknown force model %s" % (force_model_name,))
@@ -642,7 +650,8 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
         if 'relaxation_rates' not in params:
             if 'entropic' in params and params['entropic']:
                 params['relaxation_rates'] = [params['relaxation_rate']]
-            elif 'method' in params and params['method'].endswith('cumulant'):
+            elif 'method' in params and (params['method'].endswith('cumulant')
+                                         or params['method'].endswith('central_moment')):
                 params['relaxation_rates'] = [params['relaxation_rate']]
             else:
                 params['relaxation_rates'] = [params['relaxation_rate'],

lbmpy/macroscopic_value_kernels.py

@@ -19,7 +19,7 @@ def pdf_initialization_assignments(lb_method, density, velocity, pdfs,
                          + f"initialization assignments for streaming pattern {streaming_pattern}.")
 
     cqc = lb_method.conserved_quantity_computation
-    inp_eqs = cqc.equilibrium_input_equations_from_init_values(density, velocity)
+    inp_eqs = cqc.equilibrium_input_equations_from_init_values(density, velocity, force_substitution=False)
     setter_eqs = lb_method.get_equilibrium(conserved_quantity_equations=inp_eqs)
     setter_eqs = setter_eqs.new_with_substitutions({sym: field_accesses[i]
                                                     for i, sym in enumerate(lb_method.post_collision_pdf_symbols)})
@@ -186,7 +186,7 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
 
         value_map[quantity_name] = value
 
-    cq_equations = cqc.equilibrium_input_equations_from_init_values(**value_map)
+    cq_equations = cqc.equilibrium_input_equations_from_init_values(**value_map, force_substitution=False)
 
     eq = lb_method.get_equilibrium(conserved_quantity_equations=cq_equations)
     if at_least_one_field_input:
@@ -240,11 +240,13 @@ def create_advanced_velocity_setter_collision_rule(method, velocity_field: Field
     velocity_symbols = cqc.defined_symbols(order=1)[1]
 
     # density is computed from pdfs
-    eq_input_from_pdfs = cqc.equilibrium_input_equations_from_pdfs(method.pre_collision_pdf_symbols)
+    eq_input_from_pdfs = cqc.equilibrium_input_equations_from_pdfs(
+        method.pre_collision_pdf_symbols, force_substitution=False)
     eq_input_from_pdfs = eq_input_from_pdfs.new_filtered([density_symbol])
     # velocity is read from input field
     vel_symbols = [velocity_field(i) for i in range(method.dim)]
-    eq_input_from_field = cqc.equilibrium_input_equations_from_init_values(velocity=vel_symbols)
+    eq_input_from_field = cqc.equilibrium_input_equations_from_init_values(
+        velocity=vel_symbols, force_substitution=False)
     eq_input_from_field = eq_input_from_field.new_filtered(velocity_symbols)
     # then both are merged together
     eq_input = eq_input_from_pdfs.new_merged(eq_input_from_field)

lbmpy/methods/abstractlbmethod.py

@@ -3,7 +3,7 @@ from collections import namedtuple
 
 import sympy as sp
 
-from pystencils import AssignmentCollection
+from pystencils import Assignment, AssignmentCollection
 
 RelaxationInfo = namedtuple('RelaxationInfo', ['equilibrium_value', 'relaxation_rate'])
 
@@ -39,6 +39,34 @@ class AbstractLbMethod(abc.ABC):
         """Tuple of symbols representing the pdf values after collision"""
         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)):
+            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 ----------------------------------------------------------
 
     @abc.abstractmethod
@@ -59,3 +87,26 @@ 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)
+                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/centered_cumulants.py

-from lbmpy.stencils import get_stencil
-import sympy as sp
-
-from pystencils.stencil import have_same_entries
-
-from lbmpy.moments import MOMENT_SYMBOLS
-
-
-def get_default_polynomial_cumulants_for_stencil(stencil):
-    """
-    Returns default groups of cumulants to be relaxed with common relaxation rates as stated in literature.