new lbm: moment based methods - equivalence to waLBerla checked

boundaries/boundaryconditions.py

@@ -21,13 +21,25 @@ def ubb(pdfField, direction, lbMethod, velocity):
     inverseDir = invDir(direction)
 
     # TODO adapt velocity to force
-    # TODO compute density
-
-    densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
-
-    velTerm = 6 * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
-    return [sp.Eq(pdfField[neighbor](inverseDir),
-                  pdfField(direction) - velTerm)]
+    c_s_sq = sp.Rational(1, 3)
+    velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
+
+    # in conserved value computation
+    # rename what is currently called density to "virtualDensity"
+    # provide a new quantity density, which is constant in case of incompressible models
+    if lbMethod.conservedQuantityComputation._compressible:  # TODO
+        cqc = lbMethod.conservedQuantityComputation
+        densitySymbol = sp.Symbol("rho")
+        pdfFieldAccesses = [pdfField(i) for i in range(len(lbMethod.stencil))]
+        densityEquations = cqc.outputEquationsFromPdfs(pdfFieldAccesses, {'density': densitySymbol})
+        densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
+        result = densityEquations.allEquations
+        result += [sp.Eq(pdfField[neighbor](inverseDir),
+                         pdfField(direction) - velTerm * densitySymbol)]
+        return result
+    else:
+        return [sp.Eq(pdfField[neighbor](inverseDir),
+                      pdfField(direction) - velTerm)]
 
 
 def fixedDensity(pdfField, direction, lbMethod, density):

creationfunctions.py

@@ -92,11 +92,6 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
         return ValueError("'target' has to be either 'cpu' or 'gpu'")
 
     res.method = updateRule.method
-    #TODO debug begin
-    from pystencils.cpu import generateC
-    from pystencils.display_utils import highlightCpp
-    print(generateC(res))
-    #TODO debug end
     return res
 
 

methods/abstractlbmethod.py

 import abc
 import sympy as sp
+from collections import namedtuple
 from pystencils.equationcollection import EquationCollection
 
 
+RelaxationInfo = namedtuple('Relaxationinfo', ['equilibriumValue', 'relaxationRate'])
+
+
 class LbmCollisionRule(EquationCollection):
     def __init__(self, lbmMethod, *args, **kwargs):
         super(LbmCollisionRule, self).__init__(*args, **kwargs)

methods/conservedquantitycomputation.py

@@ -165,10 +165,11 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
             if tuple(velOutputSymbols) != tuple(self._symbolsOrder1):
                 mainEquations += [sp.Eq(a, b) for a, b in zip(velOutputSymbols, self._symbolsOrder1)]
             else:
-                # TODO
-                pass
+                for u_i in self._symbolsOrder1:
+                    mainEquations.append(sp.Eq(u_i, eqs[u_i]))
+                    del eqs[u_i]
 
-        eqColl = eqColl.copy(mainEquations, eqColl.allEquations)
+        eqColl = eqColl.copy(mainEquations, [sp.Eq(a, b) for a, b in eqs.items()])
         return eqColl.newWithoutUnusedSubexpressions()
 
     def __repr__(self):
@@ -277,8 +278,6 @@ def applyForceModelShift(shiftMemberName, dim, equationCollection, forceModel, c
         return equationCollection
 
 
-
-
 if __name__ == '__main__':
     from lbmpy.creationfunctions import createLatticeBoltzmannMethod
     from lbmpy.simplificationfactory import createSimplificationStrategy

methods/cumulantbased.py

+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
+
+
+class CumulantBasedLbMethod(AbstractLbMethod):
+
+    def __init__(self, stencil, cumulantToRelaxationInfoDict, conservedQuantityComputation, forceModel=None):
+        assert isinstance(conservedQuantityComputation, AbstractConservedQuantityComputation)
+        super(CumulantBasedLbMethod, self).__init__(stencil)
+
+        self._forceModel = forceModel
+        self._cumulantToRelaxationInfoDict = OrderedDict(cumulantToRelaxationInfoDict.items())
+        self._conservedQuantityComputation = conservedQuantityComputation
+
+    @property
+    def cumulantToRelaxationInfoDict(self):
+        return self._cumulantToRelaxationInfoDict
+
+    def setFirstMomentRelaxationRate(self, relaxationRate):
+        for e in MOMENT_SYMBOLS[:self.dim]:
+            assert e in self._cumulantToRelaxationInfoDict, "First cumulants are not relaxed separately by this method"
+        for e in MOMENT_SYMBOLS[:self.dim]:
+            prevEntry = self._cumulantToRelaxationInfoDict[e]
+            newEntry = RelaxationInfo(prevEntry[0], relaxationRate)
+            self._cumulantToRelaxationInfoDict[e] = newEntry
+
+    @property
+    def moments(self):
+        return tuple(self._cumulantToRelaxationInfoDict.keys())
+
+    @property
+    def momentEquilibriumValues(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()])
+
+    def _repr_html_(self):
+        table = """
+        <table style="border:none; width: 100%">
+            <tr {nb}>
+                <th {nb} >Cumulant</th>
+                <th {nb} >Eq. Value </th>
+                <th {nb} >Relaxation Time</th>
+            </tr>
+            {content}
+        </table>
+        """
+        content = ""
+        for cumulant, (eqValue, rr) in self._cumulantToRelaxationInfoDict.items():
+            vals = {
+                'rr': sp.latex(rr),
+                'cumulant': sp.latex(cumulant),
+                'eqValue': sp.latex(eqValue),
+                'nb': 'style="border:none"',
+            }
+            content += """<tr {nb}>
+                            <td {nb}>${cumulant}$</td>
+                            <td {nb}>${eqValue}$</td>
+                            <td {nb}>${rr}$</td>
+                         </tr>\n""".format(**vals)
+        return table.format(content=content, nb='style="border:none"')
+
+    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
+
+        # 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

 import sympy as sp
 import collections
-from collections import namedtuple, OrderedDict, defaultdict
+from collections import OrderedDict, defaultdict
 
 from lbmpy.stencils import stencilsHaveSameEntries, getStencil
 from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
     getMomentsOfContinuousMaxwellianEquilibrium
-from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule
+from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
 from lbmpy.moments import MOMENT_SYMBOLS, momentMatrix, isShearMoment, \
     isEven, gramSchmidt, getOrder, getDefaultMomentSetForStencil
 from pystencils.sympyextensions import commonDenominator, replaceAdditive
 
-RelaxationInfo = namedtuple('Relaxationinfo', ['equilibriumValue', 'relaxationRate'])
-
-
-def compareMomentBasedLbMethods(reference, other):
-    import ipy_table
-    table = []
-    captionRows = [len(table)]
-    table.append(['Shared Moment', 'ref', 'other', 'difference'])
-
-    referenceMoments = set(reference.moments)
-    otherMoments = set(other.moments)
-    for moment in referenceMoments.intersection(otherMoments):
-        referenceValue = reference.momentToRelaxationInfoDict[moment].equilibriumValue
-        otherValue = other.momentToRelaxationInfoDict[moment].equilibriumValue
-        diff = sp.simplify(referenceValue - otherValue)
-        table.append(["$%s$" % (sp.latex(moment), ),
-                      "$%s$" % (sp.latex(referenceValue), ),
-                      "$%s$" % (sp.latex(otherValue), ),
-                      "$%s$" % (sp.latex(diff),)])
-
-    onlyInRef = referenceMoments - otherMoments
-    if onlyInRef:
-        captionRows.append(len(table))
-        table.append(['Only in Ref', 'value', '', ''])
-        for moment in onlyInRef:
-            val = reference.momentToRelaxationInfoDict[moment].equilibriumValue
-            table.append(["$%s$" % (sp.latex(moment),),
-                          "$%s$" % (sp.latex(val),),
-                          " ", " "])
-
-    onlyInOther = otherMoments - referenceMoments
-    if onlyInOther:
-        captionRows.append(len(table))
-        table.append(['Only in Other', '', 'value', ''])
-        for moment in onlyInOther:
-            val = other.momentToRelaxationInfoDict[moment].equilibriumValue
-            table.append(["$%s$" % (sp.latex(moment),),
-                          " ",
-                          "$%s$" % (sp.latex(val),),
-                          " "])
-
-    tableDisplay = ipy_table.make_table(table)
-    for rowIdx in captionRows:
-        for col in range(4):
-            ipy_table.set_cell_style(rowIdx, col, color='#bbbbbb')
-    return tableDisplay
-
 
 class MomentBasedLbMethod(AbstractLbMethod):
     def __init__(self, stencil, momentToRelaxationInfoDict, conservedQuantityComputation, forceModel=None):
@@ -76,13 +29,13 @@ class MomentBasedLbMethod(AbstractLbMethod):
                                              the symbols used in the equilibrium moments like e.g. density and velocity
         :param forceModel: force model instance, or None if no forcing terms are required
         """
-        super(MomentBasedLbMethod, self).__init__(stencil)
-
         assert isinstance(conservedQuantityComputation, AbstractConservedQuantityComputation)
+        super(MomentBasedLbMethod, self).__init__(stencil)
 
         self._forceModel = forceModel
         self._momentToRelaxationInfoDict = OrderedDict(momentToRelaxationInfoDict.items())
         self._conservedQuantityComputation = conservedQuantityComputation
+        self._weights = None
 
         equilibriumMoments = []
         for moment, relaxInfo in momentToRelaxationInfoDict.items():
@@ -99,45 +52,13 @@ class MomentBasedLbMethod(AbstractLbMethod):
         assert len(undefinedEquilibriumSymbols) == 0, "Undefined symbol(s) in equilibrium moment: %s" % \
                                                       (undefinedEquilibriumSymbols,)
 
-        self._weights = None
-
     @property
     def momentToRelaxationInfoDict(self):
         return self._momentToRelaxationInfoDict
 
-    def setFirstMomentRelaxationRate(self, relaxationRate):
-        for e in MOMENT_SYMBOLS[:self.dim]:
-            assert e in self._momentToRelaxationInfoDict, "First moments are not relaxed separately by this method"
-        for e in MOMENT_SYMBOLS[:self.dim]:
-            prevEntry = self._momentToRelaxationInfoDict[e]
-            newEntry = RelaxationInfo(prevEntry[0], relaxationRate)
-            self._momentToRelaxationInfoDict[e] = newEntry
-
-    def _repr_html_(self):
-        table = """
-        <table style="border:none; width: 100%">
-            <tr {nb}>
-                <th {nb} >Moment</th>
-                <th {nb} >Eq. Value </th>
-                <th {nb} >Relaxation Time</th>
-            </tr>
-            {content}
-        </table>
-        """
-        content = ""
-        for moment, (eqValue, rr) in self._momentToRelaxationInfoDict.items():
-            vals = {
-                'rr': sp.latex(rr),
-                'moment': sp.latex(moment),
-                'eqValue': sp.latex(eqValue),
-                'nb': 'style="border:none"',
-            }
-            content += """<tr {nb}>
-                            <td {nb}>${moment}$</td>
-                            <td {nb}>${eqValue}$</td>
-                            <td {nb}>${rr}$</td>
-                         </tr>\n""".format(**vals)
-        return table.format(content=content, nb='style="border:none"')
+    @property
+    def conservedQuantityComputation(self):
+        return self._conservedQuantityComputation
 
     @property
     def moments(self):
@@ -165,21 +86,6 @@ class MomentBasedLbMethod(AbstractLbMethod):
             self._weights = self._computeWeights()
         return self._weights
 
-    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 getShearRelaxationRate(self):
         """
         Assumes that all shear moments are relaxed with same rate - returns this rate
@@ -213,11 +119,56 @@ class MomentBasedLbMethod(AbstractLbMethod):
             eqColl = eqColl.copy(newEqs)
         return eqColl
 
-    @property
-    def conservedQuantityComputation(self):
-        return self._conservedQuantityComputation
+    def setFirstMomentRelaxationRate(self, relaxationRate):
+        for e in MOMENT_SYMBOLS[:self.dim]:
+            assert e in self._momentToRelaxationInfoDict, "First moments are not relaxed separately by this method"
+        for e in MOMENT_SYMBOLS[:self.dim]:
+            prevEntry = self._momentToRelaxationInfoDict[e]
+            newEntry = RelaxationInfo(prevEntry[0], relaxationRate)
+            self._momentToRelaxationInfoDict[e] = newEntry
+
+    def _repr_html_(self):
+        table = """
+        <table style="border:none; width: 100%">
+            <tr {nb}>
+                <th {nb} >Moment</th>
+                <th {nb} >Eq. Value </th>
+                <th {nb} >Relaxation Time</th>
+            </tr>
+            {content}
+        </table>
+        """
+        content = ""
+        for moment, (eqValue, rr) in self._momentToRelaxationInfoDict.items():
+            vals = {
+                'rr': sp.latex(rr),
+                'moment': sp.latex(moment),
+                'eqValue': sp.latex(eqValue),
+                'nb': 'style="border:none"',
+            }
+            content += """<tr {nb}>
+                            <td {nb}>${moment}$</td>
+                            <td {nb}>${eqValue}$</td>
+                            <td {nb}>${rr}$</td>
+                         </tr>\n""".format(**vals)
+        return table.format(content=content, nb='style="border:none"')
+
+    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, D, additionalSubexpressions=[], conservedQuantityEquations=None):
+    def _getCollisionRuleWithRelaxationMatrix(self, D, additionalSubexpressions=(), conservedQuantityEquations=None):
         f = sp.Matrix(self.preCollisionPdfSymbols)
         M = momentMatrix(self.moments, self.stencil)
         m_eq = sp.Matrix(self.momentEquilibriumValues)
@@ -232,7 +183,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
         simplificationHints.update(self._conservedQuantityComputation.definedSymbols())
         simplificationHints['relaxationRates'] = D.atoms(sp.Symbol)
 
-        allSubexpressions = additionalSubexpressions + conservedQuantityEquations.allEquations
+        allSubexpressions = list(additionalSubexpressions) + conservedQuantityEquations.allEquations
         return LbmCollisionRule(self, collisionEqs, allSubexpressions,
                                 simplificationHints)
 
@@ -517,3 +468,50 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, compressible=False,
     return createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible, forceModel,
                                                  equilibriumAccuracyOrder)
 
+
+# ----------------------------------------- Comparison view for notebooks ----------------------------------------------
+
+
+def compareMomentBasedLbMethods(reference, other):
+    import ipy_table
+    table = []
+    captionRows = [len(table)]
+    table.append(['Shared Moment', 'ref', 'other', 'difference'])
+
+    referenceMoments = set(reference.moments)
+    otherMoments = set(other.moments)
+    for moment in referenceMoments.intersection(otherMoments):
+        referenceValue = reference.momentToRelaxationInfoDict[moment].equilibriumValue
+        otherValue = other.momentToRelaxationInfoDict[moment].equilibriumValue
+        diff = sp.simplify(referenceValue - otherValue)
+        table.append(["$%s$" % (sp.latex(moment), ),
+                      "$%s$" % (sp.latex(referenceValue), ),
+                      "$%s$" % (sp.latex(otherValue), ),
+                      "$%s$" % (sp.latex(diff),)])
+
+    onlyInRef = referenceMoments - otherMoments
+    if onlyInRef:
+        captionRows.append(len(table))
+        table.append(['Only in Ref', 'value', '', ''])
+        for moment in onlyInRef:
+            val = reference.momentToRelaxationInfoDict[moment].equilibriumValue
+            table.append(["$%s$" % (sp.latex(moment),),
+                          "$%s$" % (sp.latex(val),),
+                          " ", " "])
+
+    onlyInOther = otherMoments - referenceMoments
+    if onlyInOther:
+        captionRows.append(len(table))
+        table.append(['Only in Other', '', 'value', ''])
+        for moment in onlyInOther:
+            val = other.momentToRelaxationInfoDict[moment].equilibriumValue
+            table.append(["$%s$" % (sp.latex(moment),),
+                          " ",
+                          "$%s$" % (sp.latex(val),),
+                          " "])
+
+    tableDisplay = ipy_table.make_table(table)
+    for rowIdx in captionRows:
+        for col in range(4):
+            ipy_table.set_cell_style(rowIdx, col, color='#bbbbbb')
+    return tableDisplay