Finished central moments

lbmpy/creationfunctions.py

@@ -455,12 +455,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 +522,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)
 
+
 # ----------------------------------------------------------------------------------------------------------------------
 
 
@@ -643,7 +647,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/forcemodels.py

@@ -86,6 +86,7 @@ second moment nonzero   :class:`Luo`         :class:`Guo`
 =====================  ====================  =================
   
 """
+from warnings import warn
 import abc
 import sympy as sp
 
@@ -244,7 +245,9 @@ class Guo(AbstractForceModel):
         correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
         moments = (lb_method.moment_matrix * sp.Matrix(luo(lb_method)))
         central_moments = correction_factor * (lb_method.shift_matrix * moments)
-        return central_moments.expand()
+
+        # remove almost zero terms due to correction factor
+        return sp.nsimplify(central_moments.expand(), tolerance=1e-16)
 
     def equilibrium_velocity_shift(self, density):
         return default_velocity_shift(density, self._force)
@@ -300,7 +303,9 @@ class He(AbstractForceModel):
         central_moments = correction_factor * (lb_method.shift_matrix * moments)
         reduced_central_moments = [remove_higher_order_terms(central_moment, order=2, symbols=u)
                                    for central_moment in central_moments]
-        return sp.Matrix(reduced_central_moments)
+
+        # remove almost zero terms due to correction factor
+        return sp.nsimplify(sp.Matrix(reduced_central_moments), tolerance=1e-16)
 
     def equilibrium_velocity_shift(self, density):
         return default_velocity_shift(density, self._force)
@@ -313,6 +318,8 @@ class Schiller(AbstractForceModel):
     """
 
     def __init__(self, force):
+        warn("The Schiller force model is depricated and probably wrong in some cases. Please use GUO instead!",
+             DeprecationWarning)
         super(Schiller, self).__init__(force)
 
     def __call__(self, lb_method):
@@ -360,7 +367,9 @@ class Buick(AbstractForceModel):
         correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
         moments = (lb_method.moment_matrix * sp.Matrix(simple(lb_method)))
         central_moments = correction_factor * (lb_method.shift_matrix * moments)
-        return central_moments.expand()
+
+        # remove almost zero terms due to correction factor
+        return sp.nsimplify(central_moments.expand(), tolerance=1e-16)
 
     def equilibrium_velocity_shift(self, density):
         return default_velocity_shift(density, self._force)

lbmpy/methods/centeredcumulant/centeredcumulantmethod.py

@@ -11,12 +11,12 @@ from lbmpy.stencils import get_stencil
 from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
 
-from lbmpy.moments import (
-    moments_up_to_order, get_order,
-    monomial_to_polynomial_transformation_matrix,
-    moment_sort_key, exponent_tuple_sort_key,
-    exponent_to_polynomial_representation, extract_monomials, MOMENT_SYMBOLS,
-    statistical_quantity_symbol)
+from lbmpy.moments import (exponents_to_polynomial_representations,
+                           moments_up_to_order, get_order,
+                           monomial_to_polynomial_transformation_matrix,
+                           moment_sort_key, exponent_tuple_sort_key,
+                           exponent_to_polynomial_representation, extract_monomials, MOMENT_SYMBOLS,
+                           statistical_quantity_symbol)
 
 #   Local Imports
 
@@ -54,7 +54,6 @@ def cached_backward_transform(transform_obj, *args, **kwargs):
 def relax_lower_order_central_moments(moment_indices, pre_collision_values,
                                       relaxation_rates, equilibrium_values,
                                       post_collision_base=POST_COLLISION_CENTRAL_MOMENT):
-
     post_collision_symbols = [statistical_quantity_symbol(post_collision_base, i) for i in moment_indices]
     equilibrium_vec = sp.Matrix(equilibrium_values)
     moment_vec = sp.Matrix(pre_collision_values)
@@ -74,7 +73,6 @@ def relax_polynomial_cumulants(monomial_exponents, polynomials, relaxation_rates
                                pre_collision_base=PRE_COLLISION_CUMULANT,
                                post_collision_base=POST_COLLISION_CUMULANT,
                                subexpression_base='sub_col'):
-
     mon_to_poly_matrix = monomial_to_polynomial_transformation_matrix(monomial_exponents, polynomials)
     mon_vec = sp.Matrix([statistical_quantity_symbol(pre_collision_base, exp) for exp in monomial_exponents])
     equilibrium_vec = sp.Matrix(equilibrium_values)
@@ -373,7 +371,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
         k_to_c_eqs = cached_forward_transform(k_to_c_transform, simplification=pre_simplification)
         c_post_to_k_post_eqs = cached_backward_transform(
             k_to_c_transform, simplification=pre_simplification, omit_conserved_moments=True)
-        central_moments = k_to_c_transform.required_central_moments
+        central_moments = exponents_to_polynomial_representations(k_to_c_transform.required_central_moments)
         assert len(central_moments) == len(stencil), 'Number of required central moments must match stencil size.'
 
         #   3) Get Forward Transformation from PDFs to central moments

lbmpy/methods/creationfunctions.py

@@ -296,7 +296,8 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
         return create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
 
 
-def create_central_moment(stencil, relaxation_rates, maxwellian_moments=False, **kwargs):
+def create_central_moment(stencil, relaxation_rates, nested_moments=None,
+                          maxwellian_moments=False, **kwargs):
     r"""
     Creates moment based LB method where the collision takes place in the central moment space.
 
@@ -311,32 +312,42 @@ def create_central_moment(stencil, relaxation_rates, maxwellian_moments=False, *
     if isinstance(stencil, str):
         stencil = get_stencil(stencil)
 
-    dim = len(stencil[0])
+    dim = len(stencil)
 
-    moments = get_default_moment_set_for_stencil(stencil)
-    x, y, z = MOMENT_SYMBOLS
-    if dim == 2:
-        diagonal_viscous_moments = [x ** 2 + y ** 2, x ** 2 - y ** 2]
-    for i, d in enumerate(MOMENT_SYMBOLS[:dim]):
-        moments[moments.index(d ** 2)] = diagonal_viscous_moments[i]
-
-    sorted_moments = sort_moments_into_groups_of_same_order(moments)
-    if len(relaxation_rates) == len(sorted_moments) - 2:
-        relaxation_rates = [0, 0, *relaxation_rates]
-
-    if len(relaxation_rates) == len(moments):
-        rr_dict = OrderedDict(zip(moments, relaxation_rates))
-    elif len(relaxation_rates) == len(sorted_moments):
-        full_relaxation_rates_list = list()
-        for i in sorted_moments:
-            full_relaxation_rates_list.extend([relaxation_rates[i]] * len(sorted_moments[i]))
-        rr_dict = OrderedDict(zip(moments, full_relaxation_rates_list))
+    if nested_moments and not isinstance(nested_moments[0], list):
+        nested_moments = list(sort_moments_into_groups_of_same_order(nested_moments).values())
+        second_order_moments = nested_moments[2]
+        bulk_moment = [m for m in second_order_moments if is_bulk_moment(m, dim)]
+        shear_moments = [m for m in second_order_moments if is_shear_moment(m, dim)]
+        assert len(shear_moments) + len(bulk_moment) == len(second_order_moments)
+        nested_moments[2] = shear_moments
+        nested_moments.insert(3, bulk_moment)
+
+    if not nested_moments:
+        nested_moments = get_default_polynomial_cumulants_for_stencil(stencil)
+
+    if len(relaxation_rates) == 1:
+        r_rates = [relaxation_rates[0]]     # For correct viscosity
+        r_rates += [1] * (len(nested_moments) - 3)
     else:
-        raise ValueError(f"The number of relaxation rates does not fit to the method. "
-                         f"You can either choose {len(moments)} relaxation rates to relax every central moment with "
-                         f"a specific value or {len(sorted_moments)} relaxation rates to relax each order of "
-                         f"central moments or {len(sorted_moments) - 2} relaxation rates to relax the conserved "
-                         f"moments with zero and the higher order moments with the defined values.")
+        assert len(relaxation_rates) >= len(nested_moments) - 2, \
+            f"Number of relaxation rates must at least match the number of non-conserved moment groups. " \
+            f"For this stencil we have {len(nested_moments) - 2} such moment groups"
+        r_rates = relaxation_rates
+
+    rr_dict = OrderedDict()
+    rr_iter = iter(r_rates)
+    for group in nested_moments:
+        if all(get_order(c) <= 1 for c in group):
+            for moment in group:
+                rr_dict[moment] = 0
+        else:
+            try:
+                rr = next(rr_iter)
+                for moment in group:
+                    rr_dict[moment] = rr
+            except StopIteration:
+                raise ValueError('Not enough relaxation rates specified.')
 
     if maxwellian_moments:
         return create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, central_moment_space=True, **kwargs)

lbmpy/methods/momentbased/centralmomentbasedmethod.py

@@ -10,17 +10,18 @@ from lbmpy.moment_transforms import (FastCentralMomentTransform,
                                      PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT)
 
 from lbmpy.moments import (MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix,
-                           statistical_quantity_symbol, extract_monomials, exponent_tuple_sort_key)
+                           statistical_quantity_symbol)
 
 
 def relax_central_moments(moment_indices, pre_collision_values,
                           relaxation_rates, equilibrium_values,
+                          force_terms,
                           post_collision_base=POST_COLLISION_CENTRAL_MOMENT):
     post_collision_symbols = [sp.Symbol(f'{post_collision_base}_{"".join(str(i) for i in m)}') for m in moment_indices]
     equilibrium_vec = sp.Matrix(equilibrium_values)
     moment_vec = sp.Matrix(pre_collision_values)
     relaxation_matrix = sp.diag(*relaxation_rates)
-    moment_vec = moment_vec + relaxation_matrix * (equilibrium_vec - moment_vec)
+    moment_vec = moment_vec + relaxation_matrix * (equilibrium_vec - moment_vec) + force_terms
     main_assignments = [Assignment(s, eq) for s, eq in zip(post_collision_symbols, moment_vec)]
 
     return AssignmentCollection(main_assignments)
@@ -254,15 +255,20 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
         if cqe is None:
             cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
 
-        moments_as_exponents = list(extract_monomials(list(self.moments), dim=2))
-        moments_as_exponents = sorted(moments_as_exponents, key=exponent_tuple_sort_key)
+        if self._force_model is None:
+            include_force_terms = False
+
+        moment_space_forcing = False
+
+        if include_force_terms:
+            moment_space_forcing = self.force_model.has_central_moment_space_forcing
 
         #   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)
+        pdfs_to_c_transform = self.central_moment_transform_class(
+            stencil, self.moments, density, velocity, conserved_quantity_equations=cqe)
+        pdfs_to_c_eqs = pdfs_to_c_transform.forward_transform(f, simplification=pre_simplification)
 
-        moments_as_exponents = pdfs_to_k_transform.moment_exponents
+        moments_as_exponents = pdfs_to_c_transform.moment_exponents
 
         #   2) Collision
         moment_symbols = [statistical_quantity_symbol(PRE_COLLISION_CENTRAL_MOMENT, exp)
@@ -272,23 +278,29 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
         relaxation_rates = [info.relaxation_rate for info in relaxation_infos]
         equilibrium_value = [info.equilibrium_value for info in relaxation_infos]
 
-        collision_eqs = relax_central_moments(
-            moments_as_exponents, tuple(moment_symbols),
-            tuple(relaxation_rates), tuple(equilibrium_value))
+        if moment_space_forcing:
+            force_model_terms = self._force_model.central_moment_space_forcing(self)
+        else:
+            force_model_terms = sp.Matrix([0] * len(stencil))
+
+        collision_eqs = relax_central_moments(moments_as_exponents, tuple(moment_symbols),
+                                              tuple(relaxation_rates), tuple(equilibrium_value),
+                                              force_terms=force_model_terms)
 
         #   3) 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)
+        post_collision_values = self.post_collision_pdf_symbols
+        c_post_to_pdfs_eqs = pdfs_to_c_transform.backward_transform(post_collision_values,
+                                                                    simplification=pre_simplification)
 
         #   4) Now, put it all together.
-        all_acs = [] if pdfs_to_k_transform.absorbs_conserved_quantity_equations else [cqe]
-        all_acs += [pdfs_to_k_eqs, collision_eqs]
+        all_acs = [] if pdfs_to_c_transform.absorbs_conserved_quantity_equations else [cqe]
+        all_acs += [pdfs_to_c_eqs, collision_eqs]
         subexpressions = [ac.all_assignments for ac in all_acs]
-        subexpressions += k_post_to_pdfs_eqs.subexpressions
-        main_assignments = k_post_to_pdfs_eqs.main_assignments
+        subexpressions += c_post_to_pdfs_eqs.subexpressions
+        main_assignments = c_post_to_pdfs_eqs.main_assignments
 
         #   5) Maybe add forcing terms.
-        if self._force_model is not None and include_force_terms:
+        if include_force_terms and not moment_space_forcing:
             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)

lbmpy/moment_transforms/centralmomenttransforms.py

@@ -6,7 +6,8 @@ from pystencils.simp import (
 from pystencils.simp.assignment_collection import SymbolGen
 from pystencils.sympyextensions import subs_additive, fast_subs
 
-from lbmpy.moments import moment_matrix, set_up_shift_matrix, contained_moments, moments_up_to_order
+from lbmpy.moments import moment_matrix, monomial_to_polynomial_transformation_matrix,\
+    set_up_shift_matrix, contained_moments, moments_up_to_order
 from lbmpy.moments import statistical_quantity_symbol as sq_sym
 
 from .abstractmomenttransform import (
@@ -19,17 +20,17 @@ from .momenttransforms import PdfsToMomentsByChimeraTransform
 
 class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
 
-    def __init__(self, stencil, moment_exponents,
+    def __init__(self, stencil, moment_polynomials,
                  equilibrium_density,
                  equilibrium_velocity,
                  pre_collision_central_moment_base=PRE_COLLISION_CENTRAL_MOMENT,
                  post_collision_central_moment_base=POST_COLLISION_CENTRAL_MOMENT,
                  **kwargs):
-        assert len(moment_exponents) == len(stencil), 'Number of moments must match stencil'
+        assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
 
         super(PdfsToCentralMomentsByMatrix, self).__init__(
             stencil, equilibrium_density, equilibrium_velocity,
-            moment_exponents=moment_exponents, **kwargs)
+            moment_polynomials=moment_polynomials, **kwargs)
 
         moment_matrix_without_shift = moment_matrix(self.moment_exponents, self.stencil)
         shift_matrix = set_up_shift_matrix(self.moment_exponents, self.stencil, equilibrium_velocity)
@@ -37,6 +38,10 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
         self.forward_matrix = moment_matrix(self.moment_exponents, self.stencil, equilibrium_velocity)
         self.backward_matrix = moment_matrix_without_shift.inv() * shift_matrix.inv()
 
+        self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.moment_exponents,
+                                                                                self.moment_polynomials)
+        self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
+
         self.kappa_pre = pre_collision_central_moment_base
         self.kappa_post = post_collision_central_moment_base
 
@@ -52,7 +57,7 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
         simplification = self._get_simp_strategy(simplification)
 
         f_vec = sp.Matrix(pdf_symbols)
-        central_moments = self.forward_matrix * f_vec
+        central_moments = self.mono_to_poly_matrix * self.forward_matrix * f_vec
         main_assignments = [Assignment(sq_sym(self.kappa_pre, e), eq)
                             for e, eq in zip(self.moment_exponents, central_moments)]
         symbol_gen = SymbolGen(subexpression_base)
@@ -67,7 +72,7 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
 
         moments = [sq_sym(self.kappa_post, exp) for exp in self.moment_exponents]
         moment_vec = sp.Matrix(moments)
-        pdfs_from_moments = self.backward_matrix * moment_vec
+        pdfs_from_moments = self.backward_matrix * self.poly_to_mono_matrix * moment_vec
         main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, pdfs_from_moments)]
         symbol_gen = SymbolGen(subexpression_base)
 
@@ -87,25 +92,28 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
 class FastCentralMomentTransform(AbstractMomentTransform):
 
     def __init__(self, stencil,
-                 moment_exponents,
+                 moment_polynomials,
                  equilibrium_density,
                  equilibrium_velocity,
                  conserved_quantity_equations=None,
                  pre_collision_central_moment_base=PRE_COLLISION_CENTRAL_MOMENT,
                  post_collision_central_moment_base=POST_COLLISION_CENTRAL_MOMENT,
                  **kwargs):
-        assert len(moment_exponents) == len(stencil), 'Number of moments must match stencil'
+        assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
 
         super(FastCentralMomentTransform, self).__init__(
             stencil, equilibrium_density, equilibrium_velocity,
             conserved_quantity_equations=conserved_quantity_equations,
-            moment_exponents=moment_exponents, **kwargs)
+            moment_polynomials=moment_polynomials, **kwargs)
+
+        self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.moment_exponents,
+                                                                                self.moment_polynomials)
 
         self.kappa_pre = pre_collision_central_moment_base
         self.kappa_post = post_collision_central_moment_base
 
         self.mat_transform = PdfsToCentralMomentsByMatrix(
-            stencil, moment_exponents, equilibrium_density, equilibrium_velocity,
+            stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
             conserved_quantity_equations=conserved_quantity_equations,
             pre_collision_central_moment_base=pre_collision_central_moment_base,
             post_collision_central_moment_base=post_collision_central_moment_base,
@@ -130,6 +138,8 @@ class FastCentralMomentTransform(AbstractMomentTransform):
 
         subexpressions_dict = dict()
         main_assignments = []
+        central_moment_eqs = []
+        central_moment_symbols = []
 
         def collect_partial_sums(exponents, dimension=0, fixed_directions=tuple()):
             if dimension == self.dim:
@@ -148,7 +158,8 @@ class FastCentralMomentTransform(AbstractMomentTransform):
 
                 if dimension == 0:
                     lhs_symbol = sq_sym(moment_symbol_base, exponents)
-                    main_assignments.append(Assignment(lhs_symbol, summation))
+                    central_moment_eqs.append(summation)
+                    central_moment_symbols.append(lhs_symbol)
                 else:
                     lhs_symbol = _partial_kappa_symbol(fixed_directions, exponents[dimension:])
                     subexpressions_dict[lhs_symbol] = summation
@@ -157,6 +168,9 @@ class FastCentralMomentTransform(AbstractMomentTransform):
         for e in self.moment_exponents:
             collect_partial_sums(e)
 
+        central_moment_eqs = self.mono_to_poly_matrix * sp.Matrix(central_moment_eqs)
+        main_assignments = [Assignment(lhs, rhs) for lhs, rhs in zip(central_moment_symbols, central_moment_eqs)]
+
         subexpressions = [Assignment(lhs, rhs) for lhs, rhs in subexpressions_dict.items()]
         symbol_gen = SymbolGen(subexpression_base)
         ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
@@ -295,24 +309,23 @@ class FastCentralMomentTransform(AbstractMomentTransform):
 
 
 class PdfsToCentralMomentsByShiftMatrix(AbstractMomentTransform):
-    def __init__(self, stencil, moment_exponents,
+    def __init__(self, stencil, moment_polynomials,
                  equilibrium_density,
                  equilibrium_velocity,
                  conserved_quantity_equations=None,
                  pre_collision_central_moment_base=PRE_COLLISION_CENTRAL_MOMENT,
                  post_collision_central_moment_base=POST_COLLISION_CENTRAL_MOMENT,
                  **kwargs):
-        assert len(moment_exponents) == len(stencil), 'Number of moments must match stencil'
+        assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
 
         super(PdfsToCentralMomentsByShiftMatrix, self).__init__(
             stencil, equilibrium_density, equilibrium_velocity,
             conserved_quantity_equations=conserved_quantity_equations,
-            moment_exponents=moment_exponents, **kwargs)
+            moment_polynomials=moment_polynomials, **kwargs)
 
         self.raw_moment_transform = PdfsToMomentsByChimeraTransform(
-            stencil, None, equilibrium_density, equilibrium_velocity,
+            stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
             conserved_quantity_equations=conserved_quantity_equations,
-            moment_exponents=moment_exponents,
             **kwargs)
 
         self.kappa_pre = pre_collision_central_moment_base
@@ -321,6 +334,10 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractMomentTransform):
         self.shift_matrix = set_up_shift_matrix(self.moment_exponents, self.stencil, self.equilibrium_velocity)
         self.inv_shift_matrix = self.shift_matrix.inv()
 
+        self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.moment_exponents,
+                                                                                self.moment_polynomials)
+        self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
+
     @property
     def absorbs_conserved_quantity_equations(self):
         return True
@@ -361,9 +378,15 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractMomentTransform):
                 rm_to_cm_dict, self.moment_exponents, raw_moment_base, central_moment_base)
             rm_to_cm_dict = self._undo_remaining_cq_subexpressions(rm_to_cm_dict, cq_subs)
 
+        central_moment_eqs = list(rm_to_cm_dict.values())
+        central_moment_symbols = list(rm_to_cm_dict.keys())
+
+        central_moment_eqs = self.mono_to_poly_matrix * sp.Matrix(central_moment_eqs)
+        main_assignments = [Assignment(lhs, rhs) for lhs, rhs in zip(central_moment_symbols, central_moment_eqs)]
+
         subexpressions = rm_ac.all_assignments
         symbol_gen = SymbolGen(subexpression_base, dtype=float)
-        ac = AssignmentCollection(rm_to_cm_dict, subexpressions=subexpressions,
+        ac = AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions,
                                   subexpression_symbol_generator=symbol_gen)
 
         if simplification:
@@ -379,7 +402,7 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractMomentTransform):
         symbolic_rms = [sq_sym(raw_moment_base, e) for e in self.moment_exponents]
         symbolic_cms = [sq_sym(central_moment_base, e) for e in self.moment_exponents]
 
-        cm_to_rm_vec = self.inv_shift_matrix * sp.Matrix(symbolic_cms)
+        cm_to_rm_vec = self.inv_shift_matrix * self.poly_to_mono_matrix * sp.Matrix(symbolic_cms)
         cm_to_rm_dict = {rm: eq for rm, eq in zip(symbolic_rms, cm_to_rm_vec)}
 
         if simplification: