diff --git a/methods/momentbased.py b/methods/momentbased.py index 0495980d3f89f0108c4473a1d37b7921fa50738f..34b0988e88a1427bf7bf7be5a3d116547eeb8c00 100644 --- a/methods/momentbased.py +++ b/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) diff --git a/serialscenario.py b/serialscenario.py new file mode 100644 index 0000000000000000000000000000000000000000..1bb25fc885ea235d20bcbc0a4eeef4040a35e9a7 --- /dev/null +++ b/serialscenario.py @@ -0,0 +1,112 @@ +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]) + + diff --git a/viscosityComputation.py b/viscosityComputation.py new file mode 100644 index 0000000000000000000000000000000000000000..28f0f52f0136897506951d7c6c2d465c16a54e99 --- /dev/null +++ b/viscosityComputation.py @@ -0,0 +1,8 @@ + + +def relaxationRateFromLatticeViscosity(nu): + return 1.0 / (3 * nu + 0.5) + + +def latticeViscosityFromRelaxationRate(omega): + return (1/omega - 1/2) / 3