diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 2f6ecd27dbeacfbc9f1e540a4d081e1cb300091d..14e40146d7509f3ee2b28dda456a57151508ac5c 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -2,7 +2,7 @@ import sympy as sp
import numpy as np
from pystencils import TypedSymbol, Field
from pystencils.backends.cbackend import CustomCppCode
-from pystencils.ast import Block, SympyAssignment, LoopOverCoordinate, KernelFunction
+from pystencils.astnodes import Block, SympyAssignment, LoopOverCoordinate, KernelFunction
from pystencils.transformations import moveConstantsBeforeLoop, resolveFieldAccesses, typingFromSympyInspection, \
typeAllEquations
from lbmpy.boundaries.createindexlist import createBoundaryIndexList
diff --git a/creationfunctions.py b/creationfunctions.py
index 245d9f2804124da8d3be9a6227f25559a8f9bebb..6d2e8121ff52a452c25510d7400e3aa0451beae5 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -4,7 +4,7 @@ Factory functions for standard LBM methods
import sympy as sp
from copy import copy
from lbmpy.stencils import getStencil
-from lbmpy.methods.momentbased import createSRT, createTRT, createOrthogonalMRT
+from lbmpy.methods import createSRT, createTRT, createOrthogonalMRT
import lbmpy.forcemodels as forceModels
from lbmpy.simplificationfactory import createSimplificationStrategy
from lbmpy.updatekernels import createStreamPullKernel, createPdfArray
@@ -19,6 +19,8 @@ def _getParams(params, optParams):
'compressible': False,
'order': 2,
+ 'useContinuousMaxwellianEquilibrium': False,
+ 'cumulant': False,
'forceModel': 'none', # can be 'simple', 'luo' or 'guo'
'force': (0, 0, 0),
}
@@ -148,6 +150,8 @@ def createLatticeBoltzmannMethod(**params):
'compressible': params['compressible'],
'equilibriumAccuracyOrder': params['order'],
'forceModel': forceModel,
+ 'useContinuousMaxwellianEquilibrium': params['useContinuousMaxwellianEquilibrium'],
+ 'cumulant': params['cumulant'],
}
methodName = params['method']
relaxationRates = params['relaxationRates']
diff --git a/methods/__init__.py b/methods/__init__.py
index f5231bd7ed67bb704e7b5164b649c5f9243743f0..cf454ae62242ca87ce713d80d8089e760114dfbe 100644
--- a/methods/__init__.py
+++ b/methods/__init__.py
@@ -1,5 +1,5 @@
from lbmpy.methods.abstractlbmethod import AbstractLbMethod
from lbmpy.methods.momentbased import MomentBasedLbMethod, RelaxationInfo
-from lbmpy.methods.momentbased import createSRT, createTRT, createTRTWithMagicNumber, createOrthogonalMRT, \
+from lbmpy.methods.creationfunctions import createSRT, createTRT, createTRTWithMagicNumber, createOrthogonalMRT, \
createWithContinuousMaxwellianEqMoments, createWithDiscreteMaxwellianEqMoments
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index daf6bbf6d662611cbfb8697739e4466085705fb1..340605285f033d86c784a1aea95048a43f577a27 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -108,7 +108,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
return self._symbolOrder0
@property
- def firstOrderMomentSymbol(self):
+ def firstOrderMomentSymbols(self):
return self._symbolsOrder1
@property
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
new file mode 100644
index 0000000000000000000000000000000000000000..12a84f339ee1e1b1570fe488c18d60edfa9c6ff0
--- /dev/null
+++ b/methods/creationfunctions.py
@@ -0,0 +1,335 @@
+import sympy as sp
+from collections import OrderedDict, defaultdict
+from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
+from lbmpy.methods.momentbased import MomentBasedLbMethod
+from lbmpy.stencils import stencilsHaveSameEntries, getStencil
+from lbmpy.moments import isEven, gramSchmidt, getDefaultMomentSetForStencil, MOMENT_SYMBOLS, getOrder, isShearMoment
+from pystencils.sympyextensions import commonDenominator
+from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
+from lbmpy.methods.abstractlbmethod import RelaxationInfo
+from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
+ getMomentsOfContinuousMaxwellianEquilibrium, getCumulantsOfDiscreteMaxwellianEquilibrium, \
+ getCumulantsOfContinuousMaxwellianEquilibrium
+
+
+def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=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 discrete Maxwellian distribution.
+
+ :param stencil: nested tuple defining the discrete velocity space. See `func:lbmpy.stencils.getStencil`
+ :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: 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)
+ assert len(momToRrDict) == len(stencil), \
+ "The number of moments has to be the same as the number of stencil entries"
+
+ densityVelocityComputation = DensityVelocityComputation(stencil, compressible, forceModel)
+
+ 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(), 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=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`.
+ By using the continuous Maxwellian we automatically get a compressible model.
+ """
+ momToRrDict = OrderedDict(momentToRelaxationRateDict)
+ assert len(momToRrDict) == len(
+ 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)
+
+ 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]
+ 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(), eqValues)])
+ if cumulant:
+ return CumulantBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+ else:
+ return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+
+
+# ------------------------------------ SRT / TRT/ MRT Creators ---------------------------------------------------------
+
+
+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
+ :return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
+ """
+ moments = getDefaultMomentSetForStencil(stencil)
+ rrDict = OrderedDict([(m, relaxationRate) for m in moments])
+ if useContinuousMaxwellianEquilibrium:
+ return createWithContinuousMaxwellianEqMoments(stencil, rrDict, **kwargs)
+ else:
+ return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **kwargs)
+
+
+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
+ have this problem.
+
+ Parameters are similar to :func:`lbmpy.methods.createSRT`, but instead of one relaxation rate there are
+ two relaxation rates: one for even moments (determines viscosity) and one for odd moments.
+ If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.createTRTWithMagicNumber`.
+ """
+ moments = getDefaultMomentSetForStencil(stencil)
+ rrDict = OrderedDict([(m, relaxationRateEvenMoments if isEven(m) else relaxationRateOddMoments) for m in moments])
+ if useContinuousMaxwellianEquilibrium:
+ return createWithContinuousMaxwellianEqMoments(stencil, rrDict, **kwargs)
+ else:
+ return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **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, **kwargs)
+
+
+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
+ which differ by the linear combination of moments and the grouping into equal relaxation times.
+ To create a generic MRT method use :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`
+
+ :param stencil: nested tuple defining the discrete velocity space. See `func:lbmpy.stencils.getStencil`
+ :param relaxationRateGetter: function getting a list of moments as argument, returning the associated relaxation
+ time. The default returns:
+
+ - 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
+ """
+ if relaxationRateGetter is None:
+ relaxationRateGetter = defaultRelaxationRateNames()
+
+ x, y, z = MOMENT_SYMBOLS
+ one = sp.Rational(1, 1)
+
+ momentToRelaxationRateDict = OrderedDict()
+ if stencilsHaveSameEntries(stencil, getStencil("D2Q9")):
+ moments = getDefaultMomentSetForStencil(stencil)
+ orthogonalMoments = gramSchmidt(moments, stencil)
+ orthogonalMomentsScaled = [e * commonDenominator(e) for e in orthogonalMoments]
+ nestedMoments = list(sortMomentsIntoGroupsOfSameOrder(orthogonalMomentsScaled).values())
+ elif stencilsHaveSameEntries(stencil, getStencil("D3Q15")):
+ sq = x ** 2 + y ** 2 + z ** 2
+ nestedMoments = [
+ [one, x, y, z], # [0, 3, 5, 7]
+ [sq - 1], # [1]
+ [3 * sq ** 2 - 6 * sq + 1], # [2]
+ [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
+ [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 10, 11, 12, 13]
+ [x * y * z]
+ ]
+ elif stencilsHaveSameEntries(stencil, getStencil("D3Q19")):
+ sq = x ** 2 + y ** 2 + z ** 2
+ nestedMoments = [
+ [one, x, y, z], # [0, 3, 5, 7]
+ [sq - 1], # [1]
+ [3 * sq ** 2 - 6 * sq + 1], # [2]
+ [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
+ [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
+ [(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
+ [(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
+ ]
+ elif stencilsHaveSameEntries(stencil, getStencil("D3Q27")):
+ xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
+ allMoments = [
+ sp.Rational(1, 1), # 0
+ x, y, z, # 1, 2, 3
+ x * y, x * z, y * z, # 4, 5, 6
+ xsq - ysq, # 7
+ (xsq + ysq + zsq) - 3 * zsq, # 8
+ (xsq + ysq + zsq) - 2, # 9
+ 3 * (x * ysq + x * zsq) - 4 * x, # 10
+ 3 * (xsq * y + y * zsq) - 4 * y, # 11
+ 3 * (xsq * z + ysq * z) - 4 * z, # 12
+ x * ysq - x * zsq, # 13
+ xsq * y - y * zsq, # 14
+ xsq * z - ysq * z, # 15
+ x * y * z, # 16
+ 3 * (xsq * ysq + xsq * zsq + ysq * zsq) - 4 * (xsq + ysq + zsq) + 4, # 17
+ 3 * (xsq * ysq + xsq * zsq - 2 * ysq * zsq) - 2 * (2 * xsq - ysq - zsq), # 18
+ 3 * (xsq * ysq - xsq * zsq) - 2 * (ysq - zsq), # 19
+ 3 * (xsq * y * z) - 2 * (y * z), # 20
+ 3 * (x * ysq * z) - 2 * (x * z), # 21
+ 3 * (x * y * zsq) - 2 * (x * y), # 22
+ 9 * (x * ysq * zsq) - 6 * (x * ysq + x * zsq) + 4 * x, # 23
+ 9 * (xsq * y * zsq) - 6 * (xsq * y + y * zsq) + 4 * y, # 24
+ 9 * (xsq * ysq * z) - 6 * (xsq * z + ysq * z) + 4 * z, # 25
+ 27 * (xsq * ysq * zsq) - 18 * (xsq * ysq + xsq * zsq + ysq * zsq) + 12 * (xsq + ysq + zsq) - 8, # 26
+ ]
+ nestedMoments = list(sortMomentsIntoGroupsOfSameOrder(allMoments).values())
+ else:
+ raise NotImplementedError("No MRT model is available (yet) for this stencil. "
+ "Create a custom MRT using 'createWithDiscreteMaxwellianEqMoments'")
+
+ for momentList in nestedMoments:
+ rr = relaxationRateGetter(momentList)
+ for m in momentList:
+ momentToRelaxationRateDict[m] = rr
+
+ if useContinuousMaxwellianEquilibrium:
+ return createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
+ else:
+ return createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
+
+
+# ----------------------------------------- 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
+
+
+# ------------------------------------ Helper Functions ----------------------------------------------------------------
+
+
+def sortMomentsIntoGroupsOfSameOrder(moments):
+ """Returns a dictionary mapping the order (int) to a list of moments with that order."""
+ result = defaultdict(list)
+ for i, moment in enumerate(moments):
+ order = getOrder(moment)
+ result[order].append(moment)
+ return result
+
+
+def defaultRelaxationRateNames():
+ nextIndex = 0
+
+ def result(momentList):
+ nonlocal nextIndex
+ shearMomentInside = False
+ allConservedMoments = True
+ for m in momentList:
+ if isShearMoment(m):
+ shearMomentInside = True
+ if not (getOrder(m) == 0 or getOrder(m) == 1):
+ allConservedMoments = False
+
+ if shearMomentInside:
+ return sp.Symbol("omega")
+ elif allConservedMoments:
+ return 0
+ else:
+ nextIndex += 1
+ return sp.Symbol("omega_%d" % (nextIndex,))
+
+ return result
+
+
+def relaxationRateFromMagicNumber(hydrodynamicRelaxationRate, magicNumber):
+ omega = hydrodynamicRelaxationRate
+ return (4 - 2 * omega) / (4 * magicNumber * omega + 2 - omega)
+
+
+def compressibleToIncompressibleMomentValue(term, rho, u):
+ term = term.expand()
+ if term.func != sp.Add:
+ args = [term, ]
+ else:
+ args = term.args
+
+ res = 0
+ for t in args:
+ containedSymbols = t.atoms(sp.Symbol)
+ if rho in containedSymbols and len(containedSymbols.intersection(set(u))) > 0:
+ res += t / rho
+ else:
+ res += t
+ return res
diff --git a/methods/entropic.py b/methods/entropic.py
new file mode 100644
index 0000000000000000000000000000000000000000..40600ce88792dc1430a4c958fc80a3f22d4c24ee
--- /dev/null
+++ b/methods/entropic.py
@@ -0,0 +1,262 @@
+from collections import OrderedDict
+from functools import reduce
+import itertools
+import operator
+import sympy as sp
+from lbmpy.methods.creationfunctions import createWithDiscreteMaxwellianEqMoments
+from pystencils.transformations import fastSubs
+from lbmpy.stencils import getStencil
+from lbmpy.simplificationfactory import createSimplificationStrategy
+from lbmpy.moments import momentsUpToComponentOrder, MOMENT_SYMBOLS, momentsOfOrder, \
+ exponentsToPolynomialRepresentations
+
+
+def discreteEntropy(function, f_eq):
+ return -sum([f_i * sp.ln(f_i / w_i) for f_i, w_i in zip(function, f_eq)])
+
+
+def discreteApproxEntropy(function, f_eq):
+ return -sum([f_i * ((f_i / w_i)-1) for f_i, w_i in zip(function, f_eq)])
+
+
+def splitUpdateEquationsInDeltasPlusRest(updateEqsRhs, relaxationRate):
+ deltas = [ue.expand().collect(relaxationRate).coeff(relaxationRate)
+ for ue in updateEqsRhs]
+ rest = [ue.expand().collect(relaxationRate) - relaxationRate * delta for ue, delta in zip(updateEqsRhs, deltas)]
+
+ return rest, deltas
+
+
+def findEntropyMaximizingOmega(omega_s, f_eq, ds, dh):
+ dsdh = sum([ds_i * dh_i / f_eq_i for ds_i, dh_i, f_eq_i in zip(ds, dh, f_eq)])
+ dhdh = sum([dh_i * dh_i / f_eq_i for dh_i, f_eq_i in zip(dh, f_eq)])
+ return 1 - ((omega_s - 1) * dsdh / dhdh)
+
+
+def determineRelaxationRateByEntropyCondition(updateRule, omega_s, omega_h):
+ stencil = updateRule.method.stencil
+ Q = len(stencil)
+ fSymbols = updateRule.method.preCollisionPdfSymbols
+
+ updateEqsRhs = [e.rhs for e in updateRule.mainEquations]
+ _, ds = splitUpdateEquationsInDeltasPlusRest(updateEqsRhs, omega_s)
+ _, dh = splitUpdateEquationsInDeltasPlusRest(updateEqsRhs, omega_h)
+ dsSymbols = [sp.Symbol("entropicDs_%d" % (i,)) for i in range(Q)]
+ dhSymbols = [sp.Symbol("entropicDh_%d" % (i,)) for i in range(Q)]
+ feqSymbols = [sp.Symbol("entropicFeq_%d" % (i,)) for i in range(Q)]
+
+ subexprs = [sp.Eq(a, b) for a, b in zip(dsSymbols, ds)] + \
+ [sp.Eq(a, b) for a, b in zip(dhSymbols, dh)] + \
+ [sp.Eq(a, f_i + ds_i + dh_i) for a, f_i, ds_i, dh_i in zip(feqSymbols, fSymbols, dsSymbols, dhSymbols)]
+
+ optimalOmegaH = findEntropyMaximizingOmega(omega_s, feqSymbols, dsSymbols, dhSymbols)
+
+ subexprs += [sp.Eq(omega_h, optimalOmegaH)]
+
+ newUpdateEquations = []
+ for updateEq in updateRule.mainEquations:
+ index = updateRule.method.postCollisionPdfSymbols.index(updateEq.lhs)
+ newEq = sp.Eq(updateEq.lhs, fSymbols[index] + omega_s * dsSymbols[index] + omega_h * dhSymbols[index])
+ newUpdateEquations.append(newEq)
+ return updateRule.copy(newUpdateEquations, updateRule.subexpressions + subexprs)
+
+
+def decompositionByRelaxationRate(updateRule, relaxationRate):
+ lm = updateRule.latticeModel
+ stencil = lm.stencil
+
+ affineTerms = [0] * len(stencil)
+ linearTerms = [0] * len(stencil)
+ quadraticTerms = [0] * len(stencil)
+
+ for updateEquation in updateRule.updateEquations:
+ index = lm.pdfDestinationSymbols.index(updateEquation.lhs)
+ rhs = updateEquation.rhs
+ linearTerms[index] = rhs.coeff(relaxationRate)
+ quadraticTerms[index] = rhs.coeff(relaxationRate**2)
+ affineTerms[index] = rhs - relaxationRate * linearTerms[index] - relaxationRate**2 * quadraticTerms[index]
+
+ if relaxationRate in affineTerms[index].atoms(sp.Symbol):
+ raise ValueError("Relaxation Rate decomposition failed (affine part) - run simplification first")
+ if relaxationRate in linearTerms[index].atoms(sp.Symbol):
+ raise ValueError("Relaxation Rate decomposition failed (linear part) - run simplification first")
+ if relaxationRate in quadraticTerms[index].atoms(sp.Symbol):
+ raise ValueError("Relaxation Rate decomposition failed (quadratic part) - run simplification first")
+
+ return affineTerms, linearTerms, quadraticTerms
+
+
+class RelaxationRatePolynomialDecomposition:
+
+ def __init__(self, collisionRule, freeRelaxationRates, fixedRelaxationRates):
+ self._collisionRule = collisionRule
+ self._freeRelaxationRates = freeRelaxationRates
+ self._fixedRelaxationRates = fixedRelaxationRates
+ self._allRelaxationRates = fixedRelaxationRates + freeRelaxationRates
+ for se in collisionRule.subexpressions:
+ for rr in freeRelaxationRates:
+ assert rr not in se.rhs.atoms(sp.Symbol), \
+ "Decomposition not possible since free relaxation rates are already in subexpressions"
+
+ def symbolicRelaxationRateFactors(self, relaxationRate, power):
+ Q = len(self._collisionRule.method.stencil)
+ omegaIdx = self._allRelaxationRates.index(relaxationRate)
+ return [sp.Symbol("entFacOmega_%d_%d_%d" % (i, omegaIdx, power)) for i in range(Q)]
+
+ def relaxationRateFactors(self, relaxationRate):
+ updateEquations = self._collisionRule.mainEquations
+
+ result = []
+ for updateEquation in updateEquations:
+ factors = []
+ rhs = updateEquation.rhs
+ power = 0
+ while True:
+ power += 1
+ factor = rhs.coeff(relaxationRate ** power)
+ if factor != 0:
+ if relaxationRate in factor.atoms(sp.Symbol):
+ raise ValueError("Relaxation Rate decomposition failed - run simplification first")
+ factors.append(factor)
+ else:
+ break
+
+ result.append(factors)
+
+ return result
+
+ def constantExprs(self):
+ subsDict = {rr: 0 for rr in self._freeRelaxationRates}
+ subsDict.update({rr: 0 for rr in self._fixedRelaxationRates})
+ updateEquations = self._collisionRule.mainEquations
+ return [fastSubs(eq.rhs, subsDict) for eq in updateEquations]
+
+ def equilibriumExprs(self):
+ subsDict = {rr: 1 for rr in self._freeRelaxationRates}
+ subsDict.update({rr: 1 for rr in self._fixedRelaxationRates})
+ updateEquations = self._collisionRule.mainEquations
+ return [fastSubs(eq.rhs, subsDict) for eq in updateEquations]
+
+ def symbolicEquilibrium(self):
+ Q = len(self._collisionRule.method.stencil)
+ return [sp.Symbol("entFeq_%d" % (i,)) for i in range(Q)]
+
+
+def determineRelaxationRateByEntropyConditionIterative(updateRule, omega_s, omega_h,
+ newtonIterations=2, initialValue=1):
+ method = updateRule.method
+ decomp = RelaxationRatePolynomialDecomposition(updateRule, [omega_h], [omega_s])
+
+ # compute and assign relaxation rate factors
+ newUpdateEquations = []
+ fEqEqs = []
+ rrFactorDefinitions = []
+ relaxationRates = [omega_s, omega_h]
+
+ for i, constantExpr in enumerate(decomp.constantExprs()):
+ updateEqRhs = constantExpr
+ fEqRhs = constantExpr
+ for rr in relaxationRates:
+ factors = decomp.relaxationRateFactors(rr)
+ for idx, f in enumerate(factors[i]):
+ power = idx + 1
+ symbolicFactor = decomp.symbolicRelaxationRateFactors(rr, power)[i]
+ rrFactorDefinitions.append(sp.Eq(symbolicFactor, f))
+ updateEqRhs += rr ** power * symbolicFactor
+ fEqRhs += symbolicFactor
+ newUpdateEquations.append(sp.Eq(method.postCollisionPdfSymbols[i], updateEqRhs))
+ fEqEqs.append(sp.Eq(decomp.symbolicEquilibrium()[i], fEqRhs))
+
+ # newton iterations to determine free omega
+ intermediateOmegas = [sp.Symbol("omega_iter_%i" % (i,)) for i in range(newtonIterations+1)]
+ intermediateOmegas[0] = initialValue
+ intermediateOmegas[-1] = omega_h
+
+ newtonIterationEquations = []
+ for omega_idx in range(len(intermediateOmegas)-1):
+ rhsOmega = intermediateOmegas[omega_idx]
+ lhsOmega = intermediateOmegas[omega_idx+1]
+ updateEqsRhs = [e.rhs for e in newUpdateEquations]
+ entropy = discreteApproxEntropy(updateEqsRhs, [e.lhs for e in fEqEqs])
+ entropyDiff = sp.diff(entropy, omega_h)
+ entropySecondDiff = sp.diff(entropyDiff, omega_h)
+ entropyDiff = entropyDiff.subs(omega_h, rhsOmega)
+ entropySecondDiff = entropySecondDiff.subs(omega_h, rhsOmega)
+
+ newtonEq = sp.Eq(lhsOmega, rhsOmega - entropyDiff / entropySecondDiff)
+ newtonIterationEquations.append(newtonEq)
+
+ # final update equations
+ newSubexpressions = updateRule.subexpressions + rrFactorDefinitions + fEqEqs + newtonIterationEquations
+ return updateRule.copy(newUpdateEquations, newSubexpressions)
+
+
+def createKbcEntropicCollisionRule(dim, name='KBC-N4', useNewtonIterations=False, velocityRelaxation=None,
+ shearRelaxationRate=sp.Symbol("omega"),
+ higherOrderRelaxationRate=sp.Symbol("omega_h"),
+ fixedOmega=None, **kwargs):
+ def product(iterable):
+ return reduce(operator.mul, iterable, 1)
+
+ theMoment = MOMENT_SYMBOLS[:dim]
+
+ rho = [sp.Rational(1, 1)]
+ velocity = list(theMoment)
+
+ shearTensorOffDiagonal = [product(t) for t in itertools.combinations(theMoment, 2)]
+ shearTensorDiagonal = [m_i * m_i for m_i in theMoment]
+ shearTensorTrace = sum(shearTensorDiagonal)
+ shearTensorTracefreeDiagonal = [dim * d - shearTensorTrace for d in shearTensorDiagonal]
+
+ energyTransportTensor = list(exponentsToPolynomialRepresentations([a for a in momentsOfOrder(3, dim, True)
+ if 3 not in a]))
+
+ explicitlyDefined = set(rho + velocity + shearTensorOffDiagonal + shearTensorDiagonal + energyTransportTensor)
+ rest = list(set(exponentsToPolynomialRepresentations(momentsUpToComponentOrder(2, dim))) - explicitlyDefined)
+ assert len(rest) + len(explicitlyDefined) == 3**dim
+
+ # naming according to paper Karlin 2015: Entropic multirelaxation lattice Boltzmann models for turbulent flows
+ D = shearTensorOffDiagonal + shearTensorTracefreeDiagonal[:-1]
+ T = [shearTensorTrace]
+ Q = energyTransportTensor
+ if name == 'KBC-N1':
+ decomposition = [D, T+Q+rest]
+ elif name == 'KBC-N2':
+ decomposition = [D+T, Q+rest]
+ elif name == 'KBC-N3':
+ decomposition = [D+Q, T+rest]
+ elif name == 'KBC-N4':
+ decomposition = [D+T+Q, rest]
+ else:
+ raise ValueError("Unknown model. Supported models KBC-Nx")
+
+ omega_s, omega_h = shearRelaxationRate, higherOrderRelaxationRate
+ shearPart, restPart = decomposition
+
+ velRelaxation = omega_s if velocityRelaxation is None else velocityRelaxation
+ relaxationRates = [omega_s] + \
+ [velRelaxation] * len(velocity) + \
+ [omega_s] * len(shearPart) + \
+ [omega_h] * len(restPart)
+
+ stencil = getStencil("D2Q9") if dim == 2 else getStencil("D3Q27")
+ allMoments = rho + velocity + shearPart + restPart
+ momentToRr = OrderedDict((m, rr) for m, rr in zip(allMoments, relaxationRates))
+ method = createWithDiscreteMaxwellianEqMoments(stencil, momentToRr, cumulant=False, **kwargs)
+
+ simplify = createSimplificationStrategy(method)
+ collisionRule = simplify(method.getCollisionRule())
+
+ if useNewtonIterations:
+ collisionRule = determineRelaxationRateByEntropyConditionIterative(collisionRule, omega_s, omega_h, 4)
+ else:
+ collisionRule = determineRelaxationRateByEntropyCondition(collisionRule, omega_s, omega_h)
+
+ if fixedOmega:
+ collisionRule = collisionRule.newWithSubstitutions({omega_s: fixedOmega})
+
+ return collisionRule
+
+
+if __name__ == "__main__":
+ ur = createKbcEntropicCollisionRule(2, useNewtonIterations=False)
diff --git a/methods/momentbased.py b/methods/momentbased.py
index 35626f81bb93dacfda9b74078a30088fc478dc8a..1ad8be9d9993bfdc0804ceb3227e4e822b6bc084 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -1,16 +1,10 @@
import sympy as sp
-import collections
-from collections import OrderedDict, defaultdict
+from collections import OrderedDict, Sequence
-from lbmpy.stencils import stencilsHaveSameEntries, getStencil
-from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
- 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, \
- isEven, gramSchmidt, getOrder, getDefaultMomentSetForStencil
-from pystencils.sympyextensions import commonDenominator, replaceAdditive
+from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
+from lbmpy.moments import MOMENT_SYMBOLS, momentMatrix, isShearMoment
+from pystencils.sympyextensions import replaceAdditive
class MomentBasedLbMethod(AbstractLbMethod):
@@ -44,7 +38,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
symbolsInEquilibriumMoments = sp.Matrix(equilibriumMoments).atoms(sp.Symbol)
conservedQuantities = set()
for v in self._conservedQuantityComputation.definedSymbols().values():
- if isinstance(v, collections.Sequence):
+ if isinstance(v, Sequence):
conservedQuantities.update(v)
else:
conservedQuantities.add(v)
@@ -79,7 +73,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
@property
def firstOrderEquilibriumMomentSymbols(self, ):
- return self._conservedQuantityComputation.firstOrderMomentsSymbol
+ return self._conservedQuantityComputation.firstOrderMomentSymbols
@property
def weights(self):
@@ -214,329 +208,5 @@ class MomentBasedLbMethod(AbstractLbMethod):
return [], relaxationMatrix
-# ------------------------------------ Helper Functions ----------------------------------------------------------------
-def sortMomentsIntoGroupsOfSameOrder(moments):
- """Returns a dictionary mapping the order (int) to a list of moments with that order."""
- result = defaultdict(list)
- for i, moment in enumerate(moments):
- order = getOrder(moment)
- result[order].append(moment)
- return result
-
-
-def defaultRelaxationRateNames():
- nextIndex = 0
-
- def result(momentList):
- nonlocal nextIndex
- shearMomentInside = False
- allConservedMoments = True
- for m in momentList:
- if isShearMoment(m):
- shearMomentInside = True
- if not (getOrder(m) == 0 or getOrder(m) == 1):
- allConservedMoments = False
-
- if shearMomentInside:
- return sp.Symbol("omega")
- elif allConservedMoments:
- return 0
- else:
- nextIndex += 1
- return sp.Symbol("omega_%d" % (nextIndex,))
-
- return result
-
-
-def relaxationRateFromMagicNumber(hydrodynamicRelaxationRate, magicNumber):
- omega = hydrodynamicRelaxationRate
- return (4 - 2 * omega) / (4 * magicNumber * omega + 2 - omega)
-
-
-# -------------------- Generic Creators by matching equilibrium moments ------------------------------------------------
-
-def compressibleToIncompressibleMomentValue(term, rho, u):
- term = term.expand()
- if term.func != sp.Add:
- args = [term, ]
- else:
- args = term.args
-
- res = 0
- for t in args:
- containedSymbols = t.atoms(sp.Symbol)
- if rho in containedSymbols and len(containedSymbols.intersection(set(u))) > 0:
- res += t / rho
- else:
- res += t
- return res
-
-
-from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
-
-
-def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=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 discrete Maxwellian distribution.
-
- :param stencil: nested tuple defining the discrete velocity space. See `func:lbmpy.stencils.getStencil`
- :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: 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)
- assert len(momToRrDict) == len(stencil), \
- "The number of moments has to be the same as the number of stencil entries"
-
- densityVelocityComputation = DensityVelocityComputation(stencil, compressible, forceModel)
-
- 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(), 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=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`.
- By using the continuous Maxwellian we automatically get a compressible model.
- """
- momToRrDict = OrderedDict(momentToRelaxationRateDict)
- assert len(momToRrDict) == len(
- 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)
-
- 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]
- 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(), eqValues)])
- if cumulant:
- return CumulantBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
- else:
- return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
-
-
-# ------------------------------------ SRT / TRT/ MRT Creators ---------------------------------------------------------
-
-
-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
- :return: :class:`lbmpy.methods.MomentBasedLbmMethod` instance
- """
- moments = getDefaultMomentSetForStencil(stencil)
- rrDict = OrderedDict([(m, relaxationRate) for m in moments])
- if useContinuousMaxwellianEquilibrium:
- return createWithContinuousMaxwellianEqMoments(stencil, rrDict, **kwargs)
- else:
- return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **kwargs)
-
-
-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
- have this problem.
-
- Parameters are similar to :func:`lbmpy.methods.createSRT`, but instead of one relaxation rate there are
- two relaxation rates: one for even moments (determines viscosity) and one for odd moments.
- If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.createTRTWithMagicNumber`.
- """
- moments = getDefaultMomentSetForStencil(stencil)
- rrDict = OrderedDict([(m, relaxationRateEvenMoments if isEven(m) else relaxationRateOddMoments) for m in moments])
- if useContinuousMaxwellianEquilibrium:
- return createWithContinuousMaxwellianEqMoments(stencil, rrDict, **kwargs)
- else:
- return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **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, **kwargs)
-
-
-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
- which differ by the linear combination of moments and the grouping into equal relaxation times.
- To create a generic MRT method use :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`
-
- :param stencil: nested tuple defining the discrete velocity space. See `func:lbmpy.stencils.getStencil`
- :param relaxationRateGetter: function getting a list of moments as argument, returning the associated relaxation
- time. The default returns:
-
- - 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
- """
- if relaxationRateGetter is None:
- relaxationRateGetter = defaultRelaxationRateNames()
-
- x, y, z = MOMENT_SYMBOLS
- one = sp.Rational(1, 1)
-
- momentToRelaxationRateDict = OrderedDict()
- if stencilsHaveSameEntries(stencil, getStencil("D2Q9")):
- moments = getDefaultMomentSetForStencil(stencil)
- orthogonalMoments = gramSchmidt(moments, stencil)
- orthogonalMomentsScaled = [e * commonDenominator(e) for e in orthogonalMoments]
- nestedMoments = list(sortMomentsIntoGroupsOfSameOrder(orthogonalMomentsScaled).values())
- elif stencilsHaveSameEntries(stencil, getStencil("D3Q15")):
- sq = x ** 2 + y ** 2 + z ** 2
- nestedMoments = [
- [one, x, y, z], # [0, 3, 5, 7]
- [sq - 1], # [1]
- [3 * sq ** 2 - 6 * sq + 1], # [2]
- [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
- [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 10, 11, 12, 13]
- [x * y * z]
- ]
- elif stencilsHaveSameEntries(stencil, getStencil("D3Q19")):
- sq = x ** 2 + y ** 2 + z ** 2
- nestedMoments = [
- [one, x, y, z], # [0, 3, 5, 7]
- [sq - 1], # [1]
- [3 * sq ** 2 - 6 * sq + 1], # [2]
- [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
- [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
- [(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
- [(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
- ]
- elif stencilsHaveSameEntries(stencil, getStencil("D3Q27")):
- xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
- allMoments = [
- sp.Rational(1, 1), # 0
- x, y, z, # 1, 2, 3
- x * y, x * z, y * z, # 4, 5, 6
- xsq - ysq, # 7
- (xsq + ysq + zsq) - 3 * zsq, # 8
- (xsq + ysq + zsq) - 2, # 9
- 3 * (x * ysq + x * zsq) - 4 * x, # 10
- 3 * (xsq * y + y * zsq) - 4 * y, # 11
- 3 * (xsq * z + ysq * z) - 4 * z, # 12
- x * ysq - x * zsq, # 13
- xsq * y - y * zsq, # 14
- xsq * z - ysq * z, # 15
- x * y * z, # 16
- 3 * (xsq * ysq + xsq * zsq + ysq * zsq) - 4 * (xsq + ysq + zsq) + 4, # 17
- 3 * (xsq * ysq + xsq * zsq - 2 * ysq * zsq) - 2 * (2 * xsq - ysq - zsq), # 18
- 3 * (xsq * ysq - xsq * zsq) - 2 * (ysq - zsq), # 19
- 3 * (xsq * y * z) - 2 * (y * z), # 20
- 3 * (x * ysq * z) - 2 * (x * z), # 21
- 3 * (x * y * zsq) - 2 * (x * y), # 22
- 9 * (x * ysq * zsq) - 6 * (x * ysq + x * zsq) + 4 * x, # 23
- 9 * (xsq * y * zsq) - 6 * (xsq * y + y * zsq) + 4 * y, # 24
- 9 * (xsq * ysq * z) - 6 * (xsq * z + ysq * z) + 4 * z, # 25
- 27 * (xsq * ysq * zsq) - 18 * (xsq * ysq + xsq * zsq + ysq * zsq) + 12 * (xsq + ysq + zsq) - 8, # 26
- ]
- nestedMoments = list(sortMomentsIntoGroupsOfSameOrder(allMoments).values())
- else:
- raise NotImplementedError("No MRT model is available (yet) for this stencil. "
- "Create a custom MRT using 'createWithDiscreteMaxwellianEqMoments'")
-
- for momentList in nestedMoments:
- rr = relaxationRateGetter(momentList)
- for m in momentList:
- momentToRelaxationRateDict[m] = rr
-
- if useContinuousMaxwellianEquilibrium:
- return createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
- else:
- return createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
-
-
-# ----------------------------------------- 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
diff --git a/updatekernels.py b/updatekernels.py
index cab9c76027f38197e3d8e266231eae4d88219781..d27ed4f6ad0f4bad8632a90c0cff90ad5ee49910 100644
--- a/updatekernels.py
+++ b/updatekernels.py
@@ -1,4 +1,5 @@
import numpy as np
+import sympy as sp
from pystencils import Field
from pystencils.sympyextensions import fastSubs
from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
@@ -109,3 +110,9 @@ def addOutputFieldForConservedQuantities(collisionRule, conservedQuantitiesToOut
conservedQuantitiesToOutputFieldDict)
return collisionRule.merge(cqc)
+
+def writeQuantitiesToField(collisionRule, symbols, outputField):
+ if not hasattr(symbols, "__len__"):
+ symbols = [symbols]
+ eqs = [sp.Eq(outputField(i), s) for i, s in enumerate(symbols)]
+ return collisionRule.copy(collisionRule.mainEquations + eqs)