Improved simplification for moment-based method with force model

forcemodels.py

@@ -7,9 +7,9 @@ import sympy as sp
 
 
 class Simple:
-    """
+    r"""
     A simple force model which introduces the following additional force term in the
-    collision process: ::math:`3 * w_i * e_i * f_i` (often: force = rho * acceleration)
+    collision process :math:`3  w_i  e_i  f_i` (often: force = rho * acceleration)
     Should only be used with constant forces!
     Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
     """
@@ -23,12 +23,12 @@ class Simple:
 
 
 class Luo:
-    """
+    r"""
     Force model by Luo with the following forcing term
 
     .. math ::
 
-            F_i = w_i * \left( \frac{c_i - u}{c_s^2} + \frac{c_i * (c_i * u)}{c_s^4} \right) * F
+            F_i = w_i  \left( \frac{c_i - u}{c_s^2} + \frac{c_i  \left<c_i, u\right>}{c_s^4} \right)  F
 
     Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
     """
@@ -50,12 +50,12 @@ class Luo:
 
 
 class Guo:
-    """
+    r"""
      Force model by Guo with the following term:
 
     .. math ::
 
-        F_i = w_i * ( 1 - \frac{1}{2 * tau} ) * \left( \frac{c_i - u}{c_s^2} + \frac{c_i * (c_i * u)}{c_s^4} \right) * F
+        F_i = w_i  ( 1 - \frac{\omega}{2} )  \left( \frac{c_i - u}{c_s^2} + \frac{c_i \left<c_i, u\right>}{c_s^4} \right)  F
 
     Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity (both shifted by F/2)!
     """

lbmmethod.py

-import sympy as sp
-
-
-
-class MomentBasedLbmMethod(AbstractLbmMethod):
-    def __init__(self, stencil, moments, relaxationRates):
-        pass
-
-
-class FullStencilCumulantMethod(AbstractLbmMethod):
-    def __init__(self):
-        pass
-
-
-class FullStencilMomentMethod():
-
-    def __init__(self, ):
-        pass
-
-
-# -------------------------------------------------------------------------------------------------------------------
-
-
-class RelaxationScheme:
-
-    def doRelaxation(self, equationCollection):
-        pass
-
-    def equilibriumInCollisionSpace(self):
-        """Returns a map from collision space component to its equilibrium value"""
-        pass
-
-    @property
-    def conservedQuantities(self, name):
-        """Returns """
-        pass
-
-    @property
-    def symbolicRelaxationRates(self):
-        pass
-
-
-#class LbmMethod:
-#    """
-#    - splitting Relaxation and lattice model is bad - since relaxation may depend on set of methods
-#
-#    Class that holds all information about a lattice Boltzmann collision scheme:
-#        - discretization of velocity space i.e. stencil
-#        - equations to transform into collision space and back
-#        - relaxation scheme and equilibrium
-#        - conserved quantities and equations to compute them
-#
-#    Interface to force models and boundary conditions:
-#        - getMacroscopicQuantity('velocity')
-#        - getMacroscopicQuantity('density')
-#        - getMacroscopicQuantitySymbol('velocity')
-#
-#    General Idea:
-#        - separate collision space transformation (e.g. raw moments, central moments and cumulant)
-#          from relaxation schemes (e.g. various orthogonalization methods)
-#        - question: who defines macroscopic/conserved values
-#
-#    Open Questions: TODO
-#        - forcemodel as member? or can it be supplied in getCollisionRule()
-#          and getMacroscopicQuantity() -> probably a member is better
-#
-#    ForceModel Interface:
-#        - pass velocity & compressibility instead of complete lbmMethod? not sure...
-#        -
-#    Boundary Interface:
-#        - pass full lbm method to boundary function
-#
-#    """
-#
-#    def __init__(self, stencil, collisionSpaceTransformation, relaxationScheme,
-#                 inverseCollisionSpaceTransformation=None):
-#        """
-#        Create a LbmMethod
-#        :param stencil:
-#                   sequence of directions, each direction is a tuple of integers of length equal to dimension
-#        :param collisionSpaceTransformation:
-#                    EquationCollection object, defining a transformation of distribution space (pdfs)
-#                    into collision space. Collision space can be for example a moment or cumulant representation.
-#                    The equation collection must have as many free symbols as there are stencil entries.
-#                    The free symbols are termed collisionSpaceSymbols and are passed to the relaxation scheme.
-#        :param relaxationScheme:
-#                    a relaxation scheme object, providing information about the relaxation process and the
-#                    equilibrium in collision space
-#        :param inverseCollisionSpaceTransformation:
-#                    if passed None, the inverse transformation is determined by
-#                    inverting the collisionSpaceTransformation using sympy. There are
-#                    cases where sympy is slow or unable to do that, in this case the
-#                    inverse can be specified here.
-#        """
-#
-#    @property
-#    def stencil(self):
-#        pass # TODO
-#
-#    #@property
-#    #def compressible(self):
-#    #    pass # TODO
-#
-#    @property
-#    def zeroCentered(self):
-#        pass
-#
-#    @property
-#    def dim(self):
-#        return len(self.stencil[0])
-#
-#    def preCollisionSymbols(self):
-#        return sp.symbols("f_:%d" % (len(self.stencil),))
-#
-#    def collisionSpaceSymbols(self):
-#        return [sp.Symbol("c_%d" % (i,)) for i in range(len(self.stencil))]
-#
-#    def postCollisionSymbols(self):
-#        return [sp.Symbol("d_%d" % (i,)) for i in range(len(self.stencil))]
-#
-#    def equilibrium(self):
-#        """Returns equation collection, defining the equilibrium in collision space"""
-#
-#    def conservedQuantities(self):
-#        """Returns equation collection defining conserved quantities"""
-#        pass
-#
-#    def transformToCollisionSpace(self, equationCollection):
-#        pass
-#
-#    def transformToPdfSpace(self, equationCollection):
-#        pass

methods/momentbased.py

@@ -104,18 +104,16 @@ class MomentBasedLbMethod(AbstractLbMethod):
                 raise NotImplementedError("Shear moments seem to be not relaxed separately - "
                                           "Can not determine their relaxation rate automatically")
 
-    def getEquilibrium(self, conservedQuantityEquations=None):
+    def getEquilibrium(self, conservedQuantityEquations=None, includeForceTerms=True):
         D = sp.eye(len(self.relaxationRates))
-        return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations=conservedQuantityEquations)
+        return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations=conservedQuantityEquations,
+                                                          includeForceTerms=includeForceTerms)
 
     def getCollisionRule(self, conservedQuantityEquations=None):
         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)
+        eqColl = self._getCollisionRuleWithRelaxationMatrix(D, relaxationRateSubExpressions,
+                                                            True, conservedQuantityEquations)
         return eqColl
 
     def setFirstMomentRelaxationRate(self, relaxationRate):
@@ -154,7 +152,9 @@ class MomentBasedLbMethod(AbstractLbMethod):
 
     def _computeWeights(self):
         replacements = self._conservedQuantityComputation.defaultValues
-        eqColl = self.getEquilibrium().copyWithSubstitutionsApplied(replacements).insertSubexpressions()
+        eqColl = self.getEquilibrium(includeForceTerms=False).copyWithSubstitutionsApplied(replacements)
+        eqColl = eqColl.insertSubexpressions()
+
         newMainEqs = [sp.Eq(e.lhs,
                             replaceAdditive(e.rhs, 1, sum(self.preCollisionPdfSymbols), requiredMatchReplacement=1.0))
                       for e in eqColl.mainEquations]
@@ -167,7 +167,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
             weights.append(value)
         return weights
 
-    def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=(), conservedQuantityEquations=None):
+    def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=(), includeForceTerms=True,
+                                              conservedQuantityEquations=None):
         f = sp.Matrix(self.preCollisionPdfSymbols)
         M = momentMatrix(self.moments, self.stencil)
         m_eq = sp.Matrix(self.momentEquilibriumValues)
@@ -183,6 +184,17 @@ class MomentBasedLbMethod(AbstractLbMethod):
         simplificationHints['relaxationRates'] = D.atoms(sp.Symbol)
 
         allSubexpressions = list(additionalSubexpressions) + conservedQuantityEquations.allEquations
+
+        if self._forceModel is not None and includeForceTerms:
+            forceModelTerms = self._forceModel(self)
+            forceTermSymbols = sp.symbols("forceTerm_:%d" % (len(forceModelTerms,)))
+            forceSubexpressions = [sp.Eq(sym, forceModelTerm)
+                                   for sym, forceModelTerm in zip(forceTermSymbols, forceModelTerms)]
+            allSubexpressions += forceSubexpressions
+            collisionEqs = [sp.Eq(eq.lhs, eq.rhs + forceTermSymbol)
+                            for eq, forceTermSymbol in zip(collisionEqs, forceTermSymbols)]
+            simplificationHints['forceTerms'] = forceTermSymbols
+
         return LbmCollisionRule(self, collisionEqs, allSubexpressions,
                                 simplificationHints)
 

methods/momentbasedsimplifications.py

@@ -222,6 +222,9 @@ def __getCommonQuadraticAndConstantTerms(lbmCollisionEqs):
     for fa in pdfSymbols:
         t = t.subs(fa, 0)
 
+    if 'forceTerms' in sh:
+        t = t.subs({ft: 0 for ft in sh['forceTerms']})
+
     weight = t
 
     for u in sh['velocity']: