diff --git a/continuous_distribution_measures.py b/continuous_distribution_measures.py
index 1afb4859f6ecdf749d5c8cd9333928b880cba23d..866093c5e89a243a942b23f25b3df8f3ba518ad6 100644
--- a/continuous_distribution_measures.py
+++ b/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)
diff --git a/cumulants.py b/cumulants.py
index 4e8093788f3fd186e1ead084282469f066411963..1cec741ed8dcaace6bc2f0e7d191b408b9dcd9ba 100644
--- a/cumulants.py
+++ b/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}
diff --git a/forcemodels.py b/forcemodels.py
index ff767b4ff5df728e38bc5ac0899704043c75fa35..f12b470255fa9f58d3dedcb4d67f0edf038b85ef 100644
--- a/forcemodels.py
+++ b/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)]
diff --git a/maxwellian_equilibrium.py b/maxwellian_equilibrium.py
index 9af585f2a7113997988cf51fba9528c74f12e510..3fc185cfe99eb35aa7308582e2ab482873b216d8 100644
--- a/maxwellian_equilibrium.py
+++ b/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])
+
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 40b00fca6f926eaf7a783ca1c6831ec01a219a17..daf6bbf6d662611cbfb8697739e4466085705fb1 100644
--- a/methods/conservedquantitycomputation.py
+++ b/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}
diff --git a/methods/cumulantbased.py b/methods/cumulantbased.py
index 86a09e4a6fe4219fac76c9b4e3650c76491df2e2..4b6161df462b226d19a336813dec463d0e43f0dd 100644
--- a/methods/cumulantbased.py
+++ b/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
diff --git a/methods/momentbased.py b/methods/momentbased.py
index abd3903f83af20330dc623c3588a3160698e9f52..35626f81bb93dacfda9b74078a30088fc478dc8a 100644
--- a/methods/momentbased.py
+++ b/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 ----------------------------------------------
diff --git a/moments.py b/moments.py
index 600aa5c2db95be3b04d4cfdc9a2634e6e2663fd7..e7a0eff29075bf63933b7fd9234e285ee5fbb188 100644
--- a/moments.py
+++ b/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