testing equivalence to waLBerla

methods/momentbased.py

@@ -14,6 +14,51 @@ 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):
         """
@@ -56,6 +101,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
 
         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"
@@ -256,6 +305,22 @@ def relaxationRateFromMagicNumber(hydrodynamicRelaxationRate, magicNumber):
 
 # -------------------- 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
+
 
 def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=None,
                                           equilibriumAccuracyOrder=2):
@@ -284,7 +349,7 @@ def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, c
     return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
 
 
-def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, forceModel=None,
+def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=None,
                                             equilibriumAccuracyOrder=None):
     """
     Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate. These moments are
@@ -299,6 +364,12 @@ def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict,
     densityVelocityComputation = DensityVelocityComputation(stencil, True, forceModel)
     eqMoments = 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]
+
     rrDict = OrderedDict([(mom, RelaxationInfo(eqMom, rr))
                           for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqMoments)])
     return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)

serialscenario.py

+from functools import partial
+import numpy as np
+from pystencils import Field
+from pystencils.slicing import sliceFromDirection
+from lbmpy.creationfunctions import createLatticeBoltzmannFunction
+from lbmpy.macroscopicValueKernels import compileMacroscopicValuesGetter, compileMacroscopicValuesSetter
+from lbmpy.boundaries import BoundaryHandling, noSlip, ubb
+
+
+def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, preUpdateFunctions=[]):
+    ghostLayers = 1
+    domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
+    D = len(domainSize)
+
+    if 'stencil' not in methodParameters:
+        methodParameters['stencil'] = 'D2Q9' if D == 2 else 'D3Q27'
+
+    # Create kernel
+    lbmKernel = createLatticeBoltzmannFunction(**methodParameters, optimizationParams=optimizationParameters)
+    method = lbmKernel.method
+    assert D == method.dim, "Domain size and stencil do not match"
+    Q = len(method.stencil)
+
+    # Create pdf fields
+    pdfSrc = np.zeros(domainSizeWithGhostLayer + (Q,))
+    pdfDst = np.zeros_like(pdfSrc)
+
+    # Boundary setup
+    if boundarySetupFunction is not None:
+        symPdfField = Field.createFromNumpyArray('pdfs', pdfSrc, indexDimensions=1)
+        boundaryHandling = BoundaryHandling(symPdfField, domainSize, lbmKernel.method)
+        boundarySetupFunction(boundaryHandling=boundaryHandling, method=method)
+        boundaryHandling.prepare()
+    else:
+        boundaryHandling = None
+
+    # Macroscopic value input/output
+    densityArr = np.zeros(domainSizeWithGhostLayer)
+    velocityArr = np.zeros(domainSizeWithGhostLayer + (D,))
+    getMacroscopic = compileMacroscopicValuesGetter(method, ['density', 'velocity'], pdfArr=pdfSrc)
+    setMacroscopic = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0]*D}, pdfArr=pdfSrc)
+    setMacroscopic(pdfs=pdfSrc)
+
+    # Run simulation
+    def timeLoop(timeSteps):
+        nonlocal pdfSrc, pdfDst, densityArr, velocityArr
+        for t in range(timeSteps):
+            for f in preUpdateFunctions:
+                f(pdfSrc)
+            if boundaryHandling is not None:
+                boundaryHandling(pdfs=pdfSrc)
+            lbmKernel(src=pdfSrc, dst=pdfDst)
+            pdfSrc, pdfDst = pdfDst, pdfSrc
+        getMacroscopic(pdfs=pdfSrc, density=densityArr, velocity=velocityArr)
+        return pdfSrc, densityArr, velocityArr
+
+    return timeLoop
+
+
+def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, **kwargs):
+    def boundarySetupFunction(boundaryHandling, method):
+        myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:method.dim])
+        myUbb.name = 'ubb'
+        boundaryHandling.setBoundary(myUbb, sliceFromDirection('N', method.dim))
+        for direction in ('W', 'E', 'S') if method.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
+            boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters)
+
+
+def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, optimizationParameters={}, **kwargs):
+    assert dim in (2, 3)
+    kwargs['force'] = tuple([force, 0, 0][:dim])
+
+    if radius is not None:
+        assert length is not None
+        if dim == 3:
+            domainSize = (length, 2*radius+1, 2*radius+1)
+            roundChannel = True
+        else:
+            if domainSize is None:
+                domainSize = (length, 2*radius)
+    else:
+        roundChannel = False
+
+    def boundarySetupFunction(boundaryHandling, method):
+        if method.dim == 2:
+            for direction in ('N', 'S'):
+                boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+        elif method.dim == 3:
+            if roundChannel:
+                noSlipIdx = boundaryHandling.addBoundary(noSlip)
+                ff = boundaryHandling.flagField
+                yMid = ff.shape[1] // 2
+                zMid = ff.shape[2] // 2
+                y, z = np.meshgrid(range(ff.shape[1]), range(ff.shape[2]))
+                ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
+            else:
+                for direction in ('N', 'S', 'T', 'B'):
+                    boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
+
+    def periodicity(pdfArr):
+        pdfArr[0, :, :] = pdfArr[-2, :, :]
+        pdfArr[-1, :, :] = pdfArr[1, :, :]
+
+    assert domainSize is not None
+    if 'forceModel' not in kwargs:
+        kwargs['forceModel'] = 'guo'
+
+    return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters,
+                       preUpdateFunctions=[periodicity])
+
+

viscosityComputation.py

+
+
+def relaxationRateFromLatticeViscosity(nu):
+    return 1.0 / (3 * nu + 0.5)
+
+
+def latticeViscosityFromRelaxationRate(omega):
+    return (1/omega - 1/2) / 3