Added entropic SRT method

-> equilibrium is constructed according to maximum entropy principle
-> many sqrts in equilibrium formulation

chapman_enskog/chapman_enskog.py

@@ -623,7 +623,7 @@ class SteadyStateChapmanEnskogAnalysis(object):
                                          spatialDerivativeOrders=None,  # do not expand the differential operator itself
                                          pdfs=(['f', 0, order + 1],))  # expand only the 'f' terms
 
-        self.scaleHierarchy = [-B * epsDict[i] for i in range(0, 5)]
+        self.scaleHierarchy = [-B * epsDict[i] for i in range(0, order+1)]
         self.scaleHierarchyRaw = self.scaleHierarchy.copy()
 
         expandedPdfs = [feq, self.scaleHierarchy[1]]

creationfunctions.py

@@ -145,6 +145,7 @@ from functools import partial
 from lbmpy.methods import createSRT, createTRT, createOrthogonalMRT, createKBCTypeTRT, \
     createRawMRT, createThreeRelaxationRateMRT
 from lbmpy.methods.entropic import addIterativeEntropyCondition, addEntropyCondition
+from lbmpy.methods.entropic_eq_srt import createEntropicSRT
 from lbmpy.methods.relaxationrates import relaxationRateFromMagicNumber
 from lbmpy.stencils import getStencil
 import lbmpy.forcemodels as forceModels
@@ -406,6 +407,8 @@ def createLatticeBoltzmannMethod(**params):
             raise NotImplementedError("KBC type TRT methods can only be constructed for D2Q9 and D3Q27 stencils")
         methodNr = methodName[-1]
         method = createKBCTypeTRT(dim, relaxationRates[0], relaxationRates[1], 'KBC-N' + methodNr, **commonParams)
+    elif methodName.lower() == 'entropic-srt':
+        method = createEntropicSRT(stencil, relaxationRates[0], forceModel, params['compressible'])
     else:
         raise ValueError("Unknown method %s" % (methodName,))
 

methods/entropic_eq_srt.py

+import sympy as sp
+from lbmpy.maxwellian_equilibrium import getWeights
+from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule
+from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
+
+
+class EntropicEquilibriumSRT(AbstractLbMethod):
+    def __init__(self, stencil, relaxationRate, forceModel, conservedQuantityCalculation):
+        super(EntropicEquilibriumSRT, self).__init__(stencil)
+
+        self._cqc = conservedQuantityCalculation
+        self._weights = getWeights(stencil, c_s_sq=sp.Rational(1, 3))
+        self._relaxationRate = relaxationRate
+        self._forceModel = forceModel
+        self.shearRelaxationRate = relaxationRate
+
+    @property
+    def conservedQuantityComputation(self):
+        return self._cqc
+
+    @property
+    def weights(self):
+        return self._weights
+
+    @property
+    def zerothOrderEquilibriumMomentSymbol(self, ):
+        return self._cqc.zerothOrderMomentSymbol
+
+    @property
+    def firstOrderEquilibriumMomentSymbols(self, ):
+        return self._cqc.firstOrderMomentSymbols
+
+    def getEquilibrium(self, conservedQuantityEquations=None, includeForceTerms=False):
+        return self._getCollisionRuleWithRelaxationRate(1, conservedQuantityEquations=conservedQuantityEquations,
+                                                        includeForceTerms=includeForceTerms)
+
+    def _getCollisionRuleWithRelaxationRate(self, relaxationRate, includeForceTerms=True,
+                                            conservedQuantityEquations=None):
+        f = sp.Matrix(self.preCollisionPdfSymbols)
+        rho = self._cqc.zerothOrderMomentSymbol
+        u = self._cqc.firstOrderMomentSymbols
+
+        if conservedQuantityEquations is None:
+            conservedQuantityEquations = self._cqc.equilibriumInputEquationsFromPdfs(f)
+        allSubexpressions = conservedQuantityEquations.allEquations
+
+        eq = []
+        for w_i, direction in zip(self.weights, self.stencil):
+            f_i = rho * w_i
+            for u_a, e_ia in zip(u, direction):
+                B = sp.sqrt(1 + 3 * u_a ** 2)
+                f_i *= (2 - B) * ((2 * u_a + B) / (1 - u_a)) ** e_ia
+            eq.append(f_i)
+
+        collisionEqs = [sp.Eq(lhs, (1 - relaxationRate) * f_i + relaxationRate * eq_i)
+                        for lhs, f_i, eq_i in zip(self.postCollisionPdfSymbols, self.preCollisionPdfSymbols, eq)]
+
+        if (self._forceModel is not None) and includeForceTerms:
+            forceModelTerms = self._forceModel(self)
+            forceTermSymbols = sp.symbols("forceTerm_:%d" % (len(forceModelTerms, )))
+            forceSubexpressions = [sp.Eq(sym, forceModelTerm)
+                                   for sym, forceModelTerm in zip(forceTermSymbols, forceModelTerms)]
+            allSubexpressions += forceSubexpressions
+            collisionEqs = [sp.Eq(eq.lhs, eq.rhs + forceTermSymbol)
+                            for eq, forceTermSymbol in zip(collisionEqs, forceTermSymbols)]
+
+        return LbmCollisionRule(self, collisionEqs, allSubexpressions)
+
+    def getCollisionRule(self):
+        return self._getCollisionRuleWithRelaxationRate(self._relaxationRate)
+
+
+def createEntropicSRT(stencil, relaxationRate, forceModel, compressible):
+    if not compressible:
+        raise NotImplementedError("entropic-srt only implemented for compressible models")
+    densityVelocityComputation = DensityVelocityComputation(stencil, compressible, forceModel)
+    return EntropicEquilibriumSRT(stencil, relaxationRate, forceModel, densityVelocityComputation)

methods/relaxationrates.py

@@ -27,6 +27,9 @@ def getShearRelaxationRate(method):
     Shear moments in 3D are: x*y, x*z and y*z - in 2D its only x*y
     The shear relaxation rate determines the viscosity in hydrodynamic LBM schemes
     """
+    if hasattr(method, 'shearRelaxationRate'):
+        return method.shearRelaxationRate
+
     relaxationRates = set()
     for moment, relaxInfo in method.relaxationInfoDict.items():
         if isShearMoment(moment):