Cumulant implemented & bugfixes

continuous_distribution_measures.py

@@ -31,7 +31,6 @@ def momentGeneratingFunction(function, symbols, symbolsInResult):
 
     """
     assert len(symbols) == len(symbolsInResult)
-
     for t_i, v_i in zip(symbolsInResult, symbols):
         function *= sp.exp(t_i * v_i)
 

cumulants.py

@@ -162,7 +162,8 @@ def cumulantsFromPdfs(stencil, cumulantIndices=None, pdfSymbols=None):
     if cumulantIndices is None:
         cumulantIndices = momentsUpToComponentOrder(2, dim=dim)
     assert len(stencil) == len(cumulantIndices), "Stencil has to have same length as cumulantIndices sequence"
-    pdfSymbols = __getIndexedSymbols(pdfSymbols, "f", range(len(stencil)))
+    if pdfSymbols is None:
+        pdfSymbols = __getIndexedSymbols(pdfSymbols, "f", range(len(stencil)))
     return {idx: discreteCumulant(tuple(pdfSymbols), idx, stencil) for idx in cumulantIndices}
 
 

forcemodels.py

@@ -64,7 +64,6 @@ class Guo:
 
     def __call__(self, lbmMethod):
         luo = Luo(self._force)
-        u = lbmMethod.firstOrderEquilibriumMomentSymbols
         shearRelaxationRate = lbmMethod.getShearRelaxationRate()
         correctionFactor = (1 - sp.Rational(1, 2) * shearRelaxationRate)
         return [correctionFactor * t for t in luo(lbmMethod)]

maxwellian_equilibrium.py

@@ -171,8 +171,8 @@ def getMomentsOfContinuousMaxwellianEquilibrium(moments, dim, rho=sp.Symbol("rho
 
 
 @diskcache
-def getMomentsOfDiscreteMaxwellianEquilibrium(stencil,
-                                              moments, rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
+def getMomentsOfDiscreteMaxwellianEquilibrium(stencil, moments,
+                                              rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
                                               c_s_sq=sp.Symbol("c_s") ** 2, order=None, compressible=True):
     """
     Compute moments of discrete maxwellian equilibrium
@@ -195,13 +195,30 @@ def getMomentsOfDiscreteMaxwellianEquilibrium(stencil,
 # -------------------------------- Equilibrium moments -----------------------------------------------------------------
 
 
-def getCumulantsOfContinuousMaxwellianEquilibrium(cumulants, rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
-                                                  c_s_sq=sp.Symbol("c_s") ** 2, dim=3):
+def getCumulantsOfContinuousMaxwellianEquilibrium(cumulants, dim, rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
+                                                  c_s_sq=sp.Symbol("c_s") ** 2, order=None):
     from lbmpy.moments import MOMENT_SYMBOLS
     from lbmpy.continuous_distribution_measures import continuousCumulant
+    from pystencils.sympyextensions import removeHigherOrderTerms
 
-    mb = continuousMaxwellianEquilibrium(dim, rho, u, MOMENT_SYMBOLS[:dim], c_s_sq)
-    result = [continuousCumulant(mb, cumulant, MOMENT_SYMBOLS[:dim]) for cumulant in cumulants]
+    # trick to speed up sympy integration (otherwise it takes multiple minutes, or aborts):
+    # use a positive, real symbol to represent c_s_sq -> then replace this symbol afterwards with the real c_s_sq
+    c_s_sq_helper = sp.Symbol("csqHelper", positive=True, real=True)
+    mb = continuousMaxwellianEquilibrium(dim, rho, u, MOMENT_SYMBOLS[:dim], c_s_sq_helper)
+    result = [continuousCumulant(mb, cumulant, MOMENT_SYMBOLS[:dim]).subs(c_s_sq_helper, c_s_sq) for cumulant in cumulants]
+    if order is not None:
+        result = [removeHigherOrderTerms(r, order, u) for r in result]
 
     return result
 
+
+@diskcache
+def getCumulantsOfDiscreteMaxwellianEquilibrium(stencil, cumulants,
+                                                rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
+                                                c_s_sq=sp.Symbol("c_s") ** 2, order=None, compressible=True):
+    from lbmpy.cumulants import discreteCumulant
+    if order is None:
+        order = 4
+    mb = discreteMaxwellianEquilibrium(stencil, rho, u, order, c_s_sq, compressible)
+    return tuple([discreteCumulant(mb, cumulant, stencil).expand() for cumulant in cumulants])
+

methods/conservedquantitycomputation.py

@@ -103,6 +103,14 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
         else:
             return None
 
+    @property
+    def zerothOrderMomentSymbol(self):
+        return self._symbolOrder0
+
+    @property
+    def firstOrderMomentSymbol(self):
+        return self._symbolsOrder1
+
     @property
     def defaultValues(self):
         result = {self._symbolOrder0: 1}

methods/cumulantbased.py

@@ -2,8 +2,9 @@ import sympy as sp
 from collections import OrderedDict
 from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
-from lbmpy.moments import MOMENT_SYMBOLS
-from lbmpy.cumulants import cumulantsFromPdfs
+from lbmpy.moments import MOMENT_SYMBOLS, extractMonomials, momentMatrix, monomialToPolynomialTransformationMatrix
+from lbmpy.cumulants import cumulantAsFunctionOfRawMoments, rawMomentAsFunctionOfCumulants
+from pystencils.sympyextensions import fastSubs, replaceAdditive
 
 
 class CumulantBasedLbMethod(AbstractLbMethod):
@@ -15,6 +16,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
         self._forceModel = forceModel
         self._cumulantToRelaxationInfoDict = OrderedDict(cumulantToRelaxationInfoDict.items())
         self._conservedQuantityComputation = conservedQuantityComputation
+        self._weights = None
 
     @property
     def cumulantToRelaxationInfoDict(self):
@@ -29,17 +31,27 @@ class CumulantBasedLbMethod(AbstractLbMethod):
             self._cumulantToRelaxationInfoDict[e] = newEntry
 
     @property
-    def moments(self):
+    def cumulants(self):
         return tuple(self._cumulantToRelaxationInfoDict.keys())
 
     @property
-    def momentEquilibriumValues(self):
+    def cumulantEquilibriumValues(self):
         return tuple([e.equilibriumValue for e in self._cumulantToRelaxationInfoDict.values()])
 
     @property
     def relaxationRates(self):
         return tuple([e.relaxationRate for e in self._cumulantToRelaxationInfoDict.values()])
 
+    @property
+    def conservedQuantityComputation(self):
+        return self._conservedQuantityComputation
+
+    @property
+    def weights(self):
+        if self._weights is None:
+            self._weights = self._computeWeights()
+        return self._weights
+
     def _repr_html_(self):
         table = """
         <table style="border:none; width: 100%">
@@ -68,28 +80,102 @@ class CumulantBasedLbMethod(AbstractLbMethod):
 
     def getEquilibrium(self, conservedQuantityEquations=None):
         D = sp.eye(len(self.relaxationRates))
-        return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations=conservedQuantityEquations)
-
-    def getCollisionRule(self, conservedQuantityEquations=None):
-        # Step 1: transform input into cumulant space
-
-        cumulantsFromPdfs(self.stencil, )
-        # Step 2: create linear transformation from basic cumulants to requested cumulants
-
-        # Step 3: relaxation
+        return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations, False, False, False)
+
+    def getCollisionRule(self, conservedQuantityEquations=None, momentSubexpressions=False,
+                         preCollisionSubexpressions=True, postCollisionSubexpressions=False):
+        return self._getCollisionRuleWithRelaxationMatrix(sp.diag(*self.relaxationRates), conservedQuantityEquations,
+                                                          momentSubexpressions, preCollisionSubexpressions,
+                                                          postCollisionSubexpressions)
+
+    def _computeWeights(self):
+        replacements = self._conservedQuantityComputation.defaultValues
+        eqColl = self.getEquilibrium().copyWithSubstitutionsApplied(replacements).insertSubexpressions()
+        newMainEqs = [sp.Eq(e.lhs,
+                            replaceAdditive(e.rhs, 1, sum(self.preCollisionPdfSymbols), requiredMatchReplacement=1.0))
+                      for e in eqColl.mainEquations]
+        eqColl = eqColl.copy(newMainEqs)
+
+        weights = []
+        for eq in eqColl.mainEquations:
+            value = eq.rhs.expand()
+            assert len(value.atoms(sp.Symbol)) == 0, "Failed to compute weights"
+            weights.append(value)
+        return weights
+
+    def _getCollisionRuleWithRelaxationMatrix(self, relaxationMatrix, conservedQuantityEquations=None,
+                                              momentSubexpressions=False, preCollisionSubexpressions=True,
+                                              postCollisionSubexpressions=False):
+        def tupleToSymbol(exp, prefix):
+            dim = len(exp)
+            formatString = prefix + "_" + "_".join(["%d"]*dim)
+            return sp.Symbol(formatString % exp)
+
+        def substituteConservedQuantities(expressions, cqe):
+            cqe = cqe.insertSubexpressions()
+            substitutionDict = {eq.rhs: eq.lhs for eq in cqe.mainEquations}
+            density = cqe.mainEquations[0].lhs
+            substitutionDict.update({density * eq.rhs: density * eq.lhs for eq in cqe.mainEquations[1:]})
+            return [fastSubs(e, substitutionDict) for e in expressions]
+
+        f = self.preCollisionPdfSymbols
+        if conservedQuantityEquations is None:
+            conservedQuantityEquations = self._conservedQuantityComputation.equilibriumInputEquationsFromPdfs(f)
+
+        subexpressions = conservedQuantityEquations.allEquations
+
+        # 1) Determine monomial indices, and arange them such that the zeroth and first order indices come first
+        indices = list(extractMonomials(self.cumulants, dim=len(self.stencil[0])))
+        zerothMomentExponent = (0,) * self.dim
+        firstMomentExponents = [tuple([1 if i == j else 0 for i in range(self.dim)]) for j in range(self.dim)]
+        lowerOrderIndices = [zerothMomentExponent] + firstMomentExponents
+        numLowerOrderIndices = len(lowerOrderIndices)
+        assert all(e in indices for e in lowerOrderIndices), \
+            "Cumulant system does not contain relaxation rules for zeroth and first order cumulants"
+        higherOrderIndices = [e for e in indices if e not in lowerOrderIndices]
+        indices = lowerOrderIndices + higherOrderIndices  # reorder
+
+        # 2) Transform pdfs to moments
+        momentTransformationMatrix = momentMatrix(indices, self.stencil)
+        moments = momentTransformationMatrix * sp.Matrix(f)
+        moments = substituteConservedQuantities(moments, conservedQuantityEquations)
+        if momentSubexpressions:
+            symbols = [tupleToSymbol(t, "m") for t in higherOrderIndices]
+            subexpressions += [sp.Eq(sym, moment) for sym, moment in zip(symbols, moments[numLowerOrderIndices:])]
+            moments = moments[:numLowerOrderIndices] + symbols
+
+        # 3) Transform moments to monomial cumulants
+        momentsDict = {idx: m for idx, m in zip(indices, moments)}
+        monomialCumulants = [cumulantAsFunctionOfRawMoments(idx, momentsDict) for idx in indices]
+
+        if preCollisionSubexpressions:
+            symbols = [tupleToSymbol(t, "preC") for t in higherOrderIndices]
+            subexpressions += [sp.Eq(sym, c)
+                               for sym, c in zip(symbols, monomialCumulants[numLowerOrderIndices:])]
+            monomialCumulants = monomialCumulants[:numLowerOrderIndices] + symbols
+
+        # 4) Transform monomial to polynomial cumulants which are then relaxed and transformed back
+        monToPoly = monomialToPolynomialTransformationMatrix(indices, self.cumulants)
+        polyValues = monToPoly * sp.Matrix(monomialCumulants)
+        eqValues = sp.Matrix(self.cumulantEquilibriumValues)
+        collidedPolyValues = polyValues + relaxationMatrix * (eqValues - polyValues)  # collision
+        relaxedMonomialCumulants = monToPoly.inv() * collidedPolyValues
+
+        if postCollisionSubexpressions:
+            symbols = [tupleToSymbol(t, "postC") for t in higherOrderIndices]
+            subexpressions += [sp.Eq(sym, c)
+                               for sym, c in zip(symbols, relaxedMonomialCumulants[numLowerOrderIndices:])]
+            relaxedMonomialCumulants = relaxedMonomialCumulants[:numLowerOrderIndices] + symbols
+
+        # 5) Transform post-collision cumulant back to moments and from there to pdfs
+        cumulantDict = {idx: value for idx, value in zip(indices, relaxedMonomialCumulants)}
+        collidedMoments = [rawMomentAsFunctionOfCumulants(idx, cumulantDict) for idx in indices]
+        result = momentTransformationMatrix.inv() * sp.Matrix(collidedMoments)
+        mainEquations = [sp.Eq(sym, val) for sym, val in zip(self.postCollisionPdfSymbols, result)]
+
+        return LbmCollisionRule(self, mainEquations, subexpressions, simplificationHints={})
 
-        # Step 4: transform back into standard cumulants
 
-        # Step 5: transform cumulants to moments
 
-        # Step 6: transform moments to pdfs
 
-        D = sp.diag(*self.relaxationRates)
-        relaxationRateSubExpressions, D = self._generateRelaxationMatrix(D)
-        eqColl = self._getCollisionRuleWithRelaxationMatrix(D, relaxationRateSubExpressions, conservedQuantityEquations)
-        if self._forceModel is not None:
-            forceModelTerms = self._forceModel(self)
-            newEqs = [sp.Eq(eq.lhs, eq.rhs + fmt) for eq, fmt in zip(eqColl.mainEquations, forceModelTerms)]
-            eqColl = eqColl.copy(newEqs)
-        return eqColl
 

methods/momentbased.py

@@ -4,7 +4,8 @@ from collections import OrderedDict, defaultdict
 
 from lbmpy.stencils import stencilsHaveSameEntries, getStencil
 from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
-    getMomentsOfContinuousMaxwellianEquilibrium
+    getMomentsOfContinuousMaxwellianEquilibrium, getCumulantsOfDiscreteMaxwellianEquilibrium, \
+    getCumulantsOfContinuousMaxwellianEquilibrium
 from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
 from lbmpy.moments import MOMENT_SYMBOLS, momentMatrix, isShearMoment, \
@@ -74,11 +75,11 @@ class MomentBasedLbMethod(AbstractLbMethod):
 
     @property
     def zerothOrderEquilibriumMomentSymbol(self, ):
-        return self._conservedQuantityComputation.definedSymbols(order=0)[1]
+        return self._conservedQuantityComputation.zerothOrderMomentSymbol
 
     @property
     def firstOrderEquilibriumMomentSymbols(self, ):
-        return self._conservedQuantityComputation.definedSymbols(order=1)[1]
+        return self._conservedQuantityComputation.firstOrderMomentsSymbol
 
     @property
     def weights(self):
@@ -273,9 +274,12 @@ def compressibleToIncompressibleMomentValue(term, rho, u):
     return res
 
 
+from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
+
+
 def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=None,
-                                          equilibriumAccuracyOrder=2):
-    """
+                                          equilibriumAccuracyOrder=2, cumulant=False):
+    r"""
     Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate. These moments are
     relaxed against the moments of the discrete Maxwellian distribution.
 
@@ -283,9 +287,13 @@ def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, c
     :param momentToRelaxationRateDict: dict that has as many entries as the stencil. Each moment, which can be
                                        represented by an exponent tuple or in polynomial form
                                        (see `lbmpy.moments`), is mapped to a relaxation rate.
-    :param compressible: using the compressible or incompressible discrete Maxwellian
+    :param compressible: incompressible LBM methods split the density into :math:`\rho = \rho_0 + \Delta \rho`
+         where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
+         This approximates the incompressible Navier-Stokes equations better than the standard
+         compressible model.
     :param forceModel: force model instance, or None if no external forces
     :param equilibriumAccuracyOrder: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
+    :param cumulant: if True relax cumulants instead of moments
     :return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
     """
     momToRrDict = OrderedDict(momentToRelaxationRateDict)
@@ -293,16 +301,27 @@ def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, c
         "The number of moments has to be the same as the number of stencil entries"
 
     densityVelocityComputation = DensityVelocityComputation(stencil, compressible, forceModel)
-    eqMoments = getMomentsOfDiscreteMaxwellianEquilibrium(stencil, list(momToRrDict.keys()), c_s_sq=sp.Rational(1, 3),
-                                                          compressible=compressible, order=equilibriumAccuracyOrder)
+
+    if cumulant:
+        eqValues = getCumulantsOfDiscreteMaxwellianEquilibrium(stencil, list(momToRrDict.keys()),
+                                                               c_s_sq=sp.Rational(1, 3), compressible=compressible,
+                                                               order=equilibriumAccuracyOrder)
+    else:
+        eqValues = getMomentsOfDiscreteMaxwellianEquilibrium(stencil, list(momToRrDict.keys()),
+                                                             c_s_sq=sp.Rational(1, 3), compressible=compressible,
+                                                             order=equilibriumAccuracyOrder)
+
     rrDict = OrderedDict([(mom, RelaxationInfo(eqMom, rr))
-                          for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqMoments)])
-    return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+                          for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqValues)])
+    if cumulant:
+        return CumulantBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+    else:
+        return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
 
 
 def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=None,
-                                            equilibriumAccuracyOrder=None):
-    """
+                                            equilibriumAccuracyOrder=2, cumulant=False):
+    r"""
     Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate. These moments are
     relaxed against the moments of the continuous Maxwellian distribution.
     For parameter description see :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`.
@@ -313,44 +332,50 @@ def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict,
         stencil), "The number of moments has to be the same as the number of stencil entries"
     dim = len(stencil[0])
     densityVelocityComputation = DensityVelocityComputation(stencil, True, forceModel)
-    eqMoments = getMomentsOfContinuousMaxwellianEquilibrium(list(momToRrDict.keys()), dim, c_s_sq=sp.Rational(1, 3),
-                                                            order=equilibriumAccuracyOrder)
+
+    if cumulant:
+        eqValues = getCumulantsOfContinuousMaxwellianEquilibrium(list(momToRrDict.keys()), dim,
+                                                                 c_s_sq=sp.Rational(1, 3),
+                                                                 order=equilibriumAccuracyOrder)
+    else:
+        eqValues = getMomentsOfContinuousMaxwellianEquilibrium(list(momToRrDict.keys()), dim, c_s_sq=sp.Rational(1, 3),
+                                                               order=equilibriumAccuracyOrder)
 
     if not compressible:
         rho = densityVelocityComputation.definedSymbols(order=0)[1]
         u = densityVelocityComputation.definedSymbols(order=1)[1]
-        eqMoments = [compressibleToIncompressibleMomentValue(em, rho, u) for em in eqMoments]
+        eqValues = [compressibleToIncompressibleMomentValue(em, rho, u) for em in eqValues]
 
     rrDict = OrderedDict([(mom, RelaxationInfo(eqMom, rr))
-                          for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqMoments)])
-    return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+                          for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqValues)])
+    if cumulant:
+        return CumulantBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+    else:
+        return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
 
 
 # ------------------------------------ SRT / TRT/ MRT Creators ---------------------------------------------------------
 
 
-def createSRT(stencil, relaxationRate, compressible=False, forceModel=None, equilibriumAccuracyOrder=2):
+def createSRT(stencil, relaxationRate, useContinuousMaxwellianEquilibrium=False, **kwargs):
     r"""
     Creates a single relaxation time (SRT) lattice Boltzmann model also known as BGK model.
 
     :param stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.getStencil`
     :param relaxationRate: relaxation rate (inverse of the relaxation time)
                            usually called :math:`\omega` in LBM literature
-    :param compressible: incompressible LBM methods split the density into :math:`\rho = \rho_0 + \Delta \rho`
-             where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
-             This approximates the incompressible Navier-Stokes equations better than the standard
-             compressible model.
-    :param forceModel: force model instance, or None if no external forces
-    :param equilibriumAccuracyOrder: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
     :return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
     """
     moments = getDefaultMomentSetForStencil(stencil)
     rrDict = OrderedDict([(m, relaxationRate) for m in moments])
-    return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, compressible, forceModel, equilibriumAccuracyOrder)
+    if useContinuousMaxwellianEquilibrium:
+        return createWithContinuousMaxwellianEqMoments(stencil, rrDict,  **kwargs)
+    else:
+        return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **kwargs)
 
 
-def createTRT(stencil, relaxationRateEvenMoments, relaxationRateOddMoments, compressible=False,
-              forceModel=None, equilibriumAccuracyOrder=2):
+def createTRT(stencil, relaxationRateEvenMoments, relaxationRateOddMoments,
+              useContinuousMaxwellianEquilibrium=False, **kwargs):
     """
     Creates a two relaxation time (TRT) lattice Boltzmann model, where even and odd moments are relaxed differently.
     In the SRT model the exact wall position of no-slip boundaries depends on the viscosity, the TRT method does not
@@ -362,21 +387,23 @@ def createTRT(stencil, relaxationRateEvenMoments, relaxationRateOddMoments, comp
     """
     moments = getDefaultMomentSetForStencil(stencil)
     rrDict = OrderedDict([(m, relaxationRateEvenMoments if isEven(m) else relaxationRateOddMoments) for m in moments])
-    return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, compressible, forceModel, equilibriumAccuracyOrder)
+    if useContinuousMaxwellianEquilibrium:
+        return createWithContinuousMaxwellianEqMoments(stencil, rrDict,  **kwargs)
+    else:
+        return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **kwargs)
 
 
-def createTRTWithMagicNumber(stencil, relaxationRate, magicNumber=sp.Rational(3, 16), *args, **kwargs):
+def createTRTWithMagicNumber(stencil, relaxationRate, magicNumber=sp.Rational(3, 16), **kwargs):
     """
     Creates a two relaxation time (TRT) lattice Boltzmann method, where the relaxation time for odd moments is
     determines from the even moment relaxation time and a "magic number".
     For possible parameters see :func:`lbmpy.methods.createTRT`
     """
     rrOdd = relaxationRateFromMagicNumber(relaxationRate, magicNumber)
-    return createTRT(stencil, relaxationRateEvenMoments=relaxationRate, relaxationRateOddMoments=rrOdd, *args, **kwargs)
+    return createTRT(stencil, relaxationRateEvenMoments=relaxationRate, relaxationRateOddMoments=rrOdd, **kwargs)
 
 
-def createOrthogonalMRT(stencil, relaxationRateGetter=None, compressible=False,
-                        forceModel=None, equilibriumAccuracyOrder=2):
+def createOrthogonalMRT(stencil, relaxationRateGetter=None, useContinuousMaxwellianEquilibrium=False, **kwargs):
     r"""
     Returns a orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
     These MRT methods are just one specific version - there are many MRT methods possible for all these stencils
@@ -390,10 +417,6 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, compressible=False,
                                     - 0 for moments of order 0 and 1 (conserved)
                                     - :math:`\omega`: from moments of order 2 (rate that determines viscosity)
                                     - numbered :math:`\omega_i` for the rest
-
-    :param compressible: see :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`
-    :param forceModel:  see :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`
-    :param equilibriumAccuracyOrder:  see :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`
     """
     if relaxationRateGetter is None:
         relaxationRateGetter = defaultRelaxationRateNames()
@@ -465,8 +488,10 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, compressible=False,
         for m in momentList:
             momentToRelaxationRateDict[m] = rr
 
-    return createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible, forceModel,
-                                                 equilibriumAccuracyOrder)
+    if useContinuousMaxwellianEquilibrium:
+        return createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
+    else:
+        return createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
 
 
 # ----------------------------------------- Comparison view for notebooks ----------------------------------------------

moments.py

@@ -514,3 +514,47 @@ def momentEqualityTableByStencil(nameToStencilDict, moments, truncateOrder=None)
         ipy_table.set_cell_style(cellIdx[0], cellIdx[1], color=color)
 
     return tableDisplay
+
+
+def extractMonomials(sequenceOfPolynomials, dim=3):
+    """
+    Returns a set of exponent tuples of all monomials contained in the given set of polynomials
+    :param sequenceOfPolynomials: sequence of polynomials in the MOMENT_SYMBOLS
+    :param dim: length of returned exponent tuples
+
+    >>> x, y, z = MOMENT_SYMBOLS
+    >>> extractMonomials([x**2 + y**2 + y, y + y**2])
+    {(0, 2, 0), (0, 1, 0), (2, 0, 0)}
+    >>> extractMonomials([x**2 + y**2 + y, y + y**2], dim=2)
+    {(0, 1), (2, 0), (0, 2)}
+    """
+    monomials = set()
+    for polynomial in sequenceOfPolynomials:
+        for factor, exponentTuple in polynomialToExponentRepresentation(polynomial):
+            monomials.add(exponentTuple[:dim])
+    return monomials
+
+
+def monomialToPolynomialTransformationMatrix(monomials, polynomials):
+    """
+    Returns a transformation matrix from a monomial to a polynomial representation
+    :param monomials: sequence of exponent tuples
+    :param polynomials: sequence of polynomials in the MOMENT_SYMBOLS
+
+    >>> x, y, z = MOMENT_SYMBOLS
+    >>> polys = [7 * x**2 + 3 * x + 2 * y **2, \
+                 9 * x**2 - 5 * x]
+    >>> mons = list(extractMonomials(polys, dim=2))
+    >>> monomialToPolynomialTransformationMatrix(mons, polys)
+    Matrix([
+    [7,  3, 2],
+    [9, -5, 0]])
+    """
+    dim = len(monomials[0])
+
+    result = sp.zeros(len(polynomials), len(monomials))
+    for polynomialIdx, polynomial in enumerate(polynomials):
+        for factor, exponentTuple in polynomialToExponentRepresentation(polynomial):
+            exponentTuple = exponentTuple[:dim]
+            result[polynomialIdx, monomials.index(exponentTuple)] = factor
+    return result