Central moments

This MR implements the central moment collision operator. Furthermore, some clean-up in the derivation pipeline is provided.

• derivation of the central moments from the continuous equilibrium
• derivation of the central moments via shift matrix for the discrete equilibrium
• implementation of a central moment class
• finish implementation of the collision
• provide test cases for the operator
• add documentation

lbmpy/methods/momentbased/centralmomentbasedmethod.py

 import sympy as sp
-from warnings import warn
 from collections import OrderedDict
 
-from pystencils import Assignment, AssignmentCollection
+from pystencils import Assignment
 
-from pystencils.simp.simplifications import sympy_cse
+from pystencils.sympyextensions import subs_additive
+
+from lbmpy.maxwellian_equilibrium import get_weights
 
-from lbmpy.stencils import get_stencil
 from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
+from lbmpy.methods.momentbased.moment_transforms import (FastCentralMomentTransform)
 
-from lbmpy.moments import (
-    moments_up_to_order, get_order,
-    monomial_to_polynomial_transformation_matrix,
-    moment_sort_key, polynomial_to_exponent_representation, extract_monomials, MOMENT_SYMBOLS, set_up_shift_matrix)
+from lbmpy.moments import (polynomial_to_exponent_representation, MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix)
 
-from lbmpy.methods.momentbased.moment_transforms import (
-    PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT,
-    FastCentralMomentTransform)
 #   =============================== LB Method Implementation ===========================================================
 
+
 class CentralMomentBasedLbMethod(AbstractLbMethod):
     """
     Central Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT)
@@ -277,17 +273,16 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
             _, exponent_tuple = polynomial_to_exponent_representation(moment)[0]
             moments_as_exponents.add(exponent_tuple[:dim])
 
-
         #   1) Get Forward Transformation from PDFs to central moments
         pdfs_to_k_transform = self._central_moment_transform_class(
             stencil, moments_as_exponents, density, velocity, conserved_quantity_equations=cqe)
         pdfs_to_k_eqs = pdfs_to_k_transform.forward_transform(f, simplification=pre_simplification)
 
-        #   6) Get backward transformation from central moments to PDFs
+        #   2) Get backward transformation from central moments to PDFs
         d = self.post_collision_pdf_symbols
         k_post_to_pdfs_eqs = pdfs_to_k_transform.backward_transform(d, simplification=pre_simplification)
 
-        #   7) That's all. Now, put it all together.
+        #   3) Now, put it all together.
         all_acs = [] if pdfs_to_k_transform.absorbs_conserved_quantity_equations else [cqe]
         all_acs += [pdfs_to_k_eqs]
         subexpressions = [ac.all_assignments for ac in all_acs]
@@ -295,15 +290,15 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
         main_assignments = k_post_to_pdfs_eqs.main_assignments
 
         #   8) Maybe add forcing terms if CenteredCumulantForceModel was not used
-        if self._force_model is not None and \
-                not isinstance(self._force_model, CenteredCumulantForceModel) and include_force_terms:
-            force_model_terms = self._force_model(self)
-            force_term_symbols = sp.symbols(f"forceTerm_:{len(force_model_terms)}")
-            force_subexpressions = [Assignment(sym, force_model_term)
-                                    for sym, force_model_term in zip(force_term_symbols, force_model_terms)]
-            subexpressions += force_subexpressions
-            main_assignments = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
-                                for eq, force_term_symbol in zip(main_assignments, force_term_symbols)]
+        # if self._force_model is not None and \
+        #         not isinstance(self._force_model, CenteredCumulantForceModel) and include_force_terms:
+        #     force_model_terms = self._force_model(self)
+        #     force_term_symbols = sp.symbols(f"forceTerm_:{len(force_model_terms)}")
+        #     force_subexpressions = [Assignment(sym, force_model_term)
+        #                             for sym, force_model_term in zip(force_term_symbols, force_model_terms)]
+        #     subexpressions += force_subexpressions
+        #     main_assignments = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
+        #                         for eq, force_term_symbol in zip(main_assignments, force_term_symbols)]
 
         #   Aaaaaand we're done.
         return LbmCollisionRule(self, main_assignments, subexpressions)