Wolf Gladrow equilibirum construction

chapman_enskog/__init__.py

 from lbmpy.chapman_enskog.derivative import DiffOperator, Diff, expandUsingLinearity, expandUsingProductRule, \
     normalizeDiffOrder, chapmanEnskogDerivativeExpansion, chapmanEnskogDerivativeRecombination
 from lbmpy.chapman_enskog.chapman_enskog import getExpandedName, LbMethodEqMoments, insertMoments, takeMoments, \
-    CeMoment, chainSolveAndSubstitute, timeDiffSelector, momentSelector
+    CeMoment, chainSolveAndSubstitute, timeDiffSelector, momentSelector, ChapmanEnskogAnalysis

cumulants.py

@@ -4,6 +4,7 @@ Additionally functions are provided to compute cumulants from moments and vice v
 """
 
 import sympy as sp
+from pystencils.sympyextensions import scalarProduct
 from lbmpy.moments import momentsUpToComponentOrder
 from lbmpy.continuous_distribution_measures import multiDifferentiation
 from lbmpy.cache import memorycache
@@ -104,9 +105,6 @@ def __cumulantRawMomentTransform(index, dependentVarDict, outerFunction, default
 def __getDiscreteCumulantGeneratingFunction(function, stencil, waveNumbers):
     assert len(stencil) == len(function)
 
-    def scalarProduct(a, b):
-        return sum(a_i * b_i for a_i, b_i in zip(a, b))
-
     laplaceTrafo = sum([factor * sp.exp(scalarProduct(waveNumbers, e)) for factor, e in zip(function, stencil)])
     return sp.ln(laplaceTrafo)
 

maxwellian_equilibrium.py

@@ -192,6 +192,24 @@ def getMomentsOfDiscreteMaxwellianEquilibrium(stencil, moments,
     return tuple([discreteMoment(mb, moment, stencil).expand() for moment in moments])
 
 
+def compressibleToIncompressibleMomentValue(term, rho, u):
+    term = sp.sympify(term)
+    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
+
+
 # -------------------------------- Equilibrium moments -----------------------------------------------------------------
 
 

methods/conservedquantitycomputation.py

@@ -288,16 +288,3 @@ def applyForceModelShift(shiftMemberName, dim, equationCollection, forceModel, c
         return equationCollection.copy(newEqs)
     else:
         return equationCollection
-
-
-if __name__ == '__main__':
-    from lbmpy.creationfunctions import createLatticeBoltzmannMethod
-    from lbmpy.simplificationfactory import createSimplificationStrategy
-    from lbmpy.stencils import getStencil
-    from lbmpy_old.lbmgenerator import createStreamCollideUpdateRule
-    from lbmpy_old.latticemodel import makeSRT
-    import sympy as sp
-    methodNew = createLatticeBoltzmannMethod(compressible=True)
-    newSimp = createSimplificationStrategy(methodNew)
-    cqc = methodNew.conservedQuantityComputation
-    cqc.outputEquationsFromPdfs(sp.symbols("f_:9"), {'density': sp.Symbol("rho_out")})
\ No newline at end of file

methods/creationfunctions.py

@@ -16,7 +16,7 @@ from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputatio
 from lbmpy.methods.abstractlbmethod import RelaxationInfo
 from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
     getMomentsOfContinuousMaxwellianEquilibrium, getCumulantsOfDiscreteMaxwellianEquilibrium, \
-    getCumulantsOfContinuousMaxwellianEquilibrium
+    getCumulantsOfContinuousMaxwellianEquilibrium, compressibleToIncompressibleMomentValue
 from lbmpy.methods.relaxationrates import relaxationRateFromMagicNumber, defaultRelaxationRateNames
 
 
@@ -423,23 +423,3 @@ def compareMomentBasedLbMethods(reference, other, showDeviationsOnly=False):
             ipy_table.set_cell_style(rowIdx, col, color='#bbbbbb')
     return tableDisplay
 
-
-# ------------------------------------ Helper Functions ----------------------------------------------------------------
-
-
-def compressibleToIncompressibleMomentValue(term, rho, u):
-    term = sp.sympify(term)
-    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

moments.py

@@ -245,6 +245,21 @@ def isEven(moment):
         return sp.expand(moment - opposite) == 0
 
 
+def getMomentIndices(momentExponentTuple):
+    """Returns indices for a given exponent tuple:
+    
+    Example:
+        >>> getMomentIndices((2,1,0))
+        [0, 0, 1]
+        >>> getMomentIndices((0,0,3))
+        [2, 2, 2]
+    """
+    result = []
+    for i, element in enumerate(momentExponentTuple):
+        result += [i] * element
+    return result
+
+
 def getOrder(moment):
     """
     Computes polynomial order of given moment

quadratic_equilibrium_construction.py

+"""
+Method to construct a quadratic equilibrium using a generic quadratic ansatz according to the book of 
+Wolf-Gladrow, section 5.4
+"""
+
+import sympy as sp
+import numpy as np
+from pystencils.sympyextensions import scalarProduct
+from lbmpy.moments import discreteMoment
+from lbmpy.maxwellian_equilibrium import compressibleToIncompressibleMomentValue
+
+
+def genericEquilibriumAnsatz(stencil, u=sp.symbols("u_:3")):
+    """Returns a generic quadratic equilibrium with coefficients A, B, C, D according to 
+    Wolf Gladrow Book equation (5.4.1) """
+    dim = len(stencil[0])
+    u = u[:dim]
+
+    equilibrium = []
+
+    for direction in stencil:
+        speed = np.abs(direction).sum()
+        weight, linear, mixQuadratic, quadratic = getParameterSymbols(speed)
+        uTimesD = scalarProduct(u, direction)
+        eq = weight + linear * uTimesD + mixQuadratic * uTimesD ** 2 + quadratic * scalarProduct(u, u)
+        equilibrium.append(eq)
+    return tuple(equilibrium)
+
+
+def matchGenericEquilibriumAnsatz(stencil, equilibrium, u=sp.symbols("u_:3")):
+    """Given a quadratic equilibrium, the generic coefficients A,B,C,D are determined. 
+    Returns a dict that maps these coefficients to their values. If the equilibrium does not have a
+    generic quadratic form, a ValueError is raised"""
+    dim = len(stencil[0])
+    u = u[:dim]
+
+    result = dict()
+    for direction, actualEquilibrium in zip(stencil, equilibrium):
+        speed = np.abs(direction).sum()
+        A, B, C, D = getParameterSymbols(speed)
+        uTimesD = scalarProduct(u, direction)
+        genericEquation = A + B * uTimesD + C * uTimesD ** 2 + D * scalarProduct(u, u)
+
+        equations = sp.poly(actualEquilibrium - genericEquation, *u).coeffs()
+        solveRes = sp.solve(equations, [A, B, C, D])
+        if not solveRes:
+            raise ValueError("This equilibrium does not match the generic quadratic standard form")
+        for dof, value in solveRes.items():
+            if dof in result and result[dof] != value:
+                raise ValueError("This equilibrium does not match the generic quadratic standard form")
+            result[dof] = value
+
+    return result
+
+
+def momentConstraintEquations(stencil, equilibrium, momentToValueDict, u=sp.symbols("u_:3")):
+    """Returns a set of equations that have to be fulfilled for a generic equilibrium match moment conditions 
+    passed in momentToValueDict. This dict is expected to map moment tuples to values."""
+    dim = len(stencil[0])
+    u = u[:dim]
+    constraintEquations = set()
+    for moment, desiredValue in momentToValueDict.items():
+        genericMoment = discreteMoment(equilibrium, moment, stencil)
+        equations = sp.poly(genericMoment - desiredValue, *u).coeffs()
+        constraintEquations.update(equations)
+    return constraintEquations
+
+
+def hydrodynamicMomentValues(upToOrder=3, dim=3, compressible=True):
+    """Returns the values of moments that are required to approximate Navier Stokes (if upToOrder=3)"""
+    from lbmpy.maxwellian_equilibrium import getMomentsOfContinuousMaxwellianEquilibrium
+    from lbmpy.moments import momentsUpToOrder
+
+    moms = momentsUpToOrder(upToOrder, dim)
+    c_s_sq = sp.Symbol("p") / sp.Symbol("rho")
+    momValues = getMomentsOfContinuousMaxwellianEquilibrium(moms, dim=dim, c_s_sq=c_s_sq, order=2)
+    if not compressible:
+        momValues = [compressibleToIncompressibleMomentValue(m, sp.Symbol("rho"), sp.symbols("u_:3")[:dim])
+                     for m in momValues]
+
+    return {a: b.expand() for a, b in zip(moms, momValues)}
+
+
+def getParameterSymbols(i):
+    return sp.symbols("A_%d B_%d C_%d D_%d" % (i, i, i, i))