Tutorial on LB method specification

methods/creationfunctions.py

@@ -99,13 +99,20 @@ def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict,
         return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
 
 
-def createWithGivenEqMoments(stencil, momentToEqValueDict, compressible=False, forceModel=None,
-                             momentToRelaxationRateDict=defaultdict(lambda: sp.Symbol("omega"))):
+def createGenericMRT(stencil, momentEqValueRelaxationRateTuples, compressible=False, forceModel=None):
+    r"""
+    Creates a generic moment-based LB method
+    :param stencil: sequence of lattice velocities
+    :param momentEqValueRelaxationRateTuples: sequence of tuples containing (moment, equilibrium value, relax. rate)
+    :param compressible: compressibility, determines calculation of velocity for force models
+    :param forceModel: see createWithDiscreteMaxwellianEqMoments
+    """
     densityVelocityComputation = DensityVelocityComputation(stencil, compressible, forceModel)
 
     rrDict = OrderedDict()
-    for moment in momentToEqValueDict.keys():
-        rrDict[moment] = RelaxationInfo(momentToEqValueDict[moment], momentToRelaxationRateDict[moment])
+    for moment, eqValue, rr in momentEqValueRelaxationRateTuples:
+        moment = sp.sympify(moment)
+        rrDict[moment] = RelaxationInfo(eqValue, rr)
     return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
 
 

methods/cumulantbased.py

@@ -136,7 +136,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
 
         subexpressions = conservedQuantityEquations.allEquations
 
-        # 1) Determine monomial indices, and arange them such that the zeroth and first order indices come first
+        # 1) Determine monomial indices, and arrange 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)]

methods/momentbased.py

@@ -35,7 +35,6 @@ class MomentBasedLbMethod(AbstractLbMethod):
         equilibriumMoments = []
         for moment, relaxInfo in momentToRelaxationInfoDict.items():
             equilibriumMoments.append(relaxInfo.equilibriumValue)
-        symbolsInEquilibriumMoments = sp.Matrix(equilibriumMoments).atoms(sp.Symbol)
         conservedQuantities = set()
         for v in self._conservedQuantityComputation.definedSymbols().values():
             if isinstance(v, Sequence):
@@ -43,11 +42,6 @@ class MomentBasedLbMethod(AbstractLbMethod):
             else:
                 conservedQuantities.add(v)
 
-        undefinedEquilibriumSymbols = symbolsInEquilibriumMoments - conservedQuantities
-        #TODO
-        #assert len(undefinedEquilibriumSymbols) == 0, "Undefined symbol(s) in equilibrium moment: %s" % \
-        #                                             (undefinedEquilibriumSymbols,)
-
     @property
     def forceModel(self):
         return self._forceModel
@@ -116,6 +110,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
             newEntry = RelaxationInfo(prevEntry[0], relaxationRate)
             self._momentToRelaxationInfoDict[e] = newEntry
 
+    @property
+    def momentMatrix(self):
+        return momentMatrix(self.moments, self.stencil)
+
     def _repr_html_(self):
         table = """
         <table style="border:none; width: 100%">
@@ -162,7 +160,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
     def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=(), includeForceTerms=True,
                                               conservedQuantityEquations=None):
         f = sp.Matrix(self.preCollisionPdfSymbols)
-        M = momentMatrix(self.moments, self.stencil)
+        M = self.momentMatrix
         m_eq = sp.Matrix(self.momentEquilibriumValues)
 
         collisionRule = f + M.inv() * D * (m_eq - M * f)