From ed79ce1b57c0b55a6d3b14d69d044d559b89f258 Mon Sep 17 00:00:00 2001
From: Martin Bauer <bauer_martin@gmx.de>
Date: Tue, 31 Jan 2017 14:47:25 +0100
Subject: [PATCH] new lbm: moment based methods - equivalence to waLBerla
checked
---
boundaries/boundaryconditions.py | 26 ++-
creationfunctions.py | 5 -
methods/abstractlbmethod.py | 4 +
methods/conservedquantitycomputation.py | 9 +-
methods/cumulantbased.py | 95 +++++++++++
methods/momentbased.py | 210 ++++++++++++------------
6 files changed, 226 insertions(+), 123 deletions(-)
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 9650abdc..d65ad424 100644
--- a/boundaries/boundaryconditions.py
+++ b/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):
diff --git a/creationfunctions.py b/creationfunctions.py
index 89e59848..245d9f28 100644
--- a/creationfunctions.py
+++ b/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
diff --git a/methods/abstractlbmethod.py b/methods/abstractlbmethod.py
index 201ab054..de64dac5 100644
--- a/methods/abstractlbmethod.py
+++ b/methods/abstractlbmethod.py
@@ -1,8 +1,12 @@
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)
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 2cd989ba..40b00fca 100644
--- a/methods/conservedquantitycomputation.py
+++ b/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
diff --git a/methods/cumulantbased.py b/methods/cumulantbased.py
index e69de29b..86a09e4a 100644
--- a/methods/cumulantbased.py
+++ b/methods/cumulantbased.py
@@ -0,0 +1,95 @@
+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
+
diff --git a/methods/momentbased.py b/methods/momentbased.py
index 34b0988e..abd3903f 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -1,63 +1,16 @@
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
--
GitLab