diff --git a/creationfunctions.py b/creationfunctions.py index 390bf5642a526847478afe1acae37579469edc6d..fe5ca7a0415cec4a03d139c52399d051fd9b09e8 100644 --- a/creationfunctions.py +++ b/creationfunctions.py @@ -68,7 +68,7 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs): if 'openMP' in optParams: if isinstance(optParams['openMP'], bool) and optParams['openMP']: addOpenMP(ast) - elif isinstance(optParams['openMP'], int): + elif not isinstance(optParams['openMP'], bool) and isinstance(optParams['openMP'], int): addOpenMP(ast, numThreads=optParams['openMP']) res = makePythonCpuFunction(ast) elif optParams['target'] == 'gpu': diff --git a/macroscopic_value_kernels.py b/macroscopic_value_kernels.py index c05a457b56a06ee28e1569394d14db70db605c7a..92ec1b73e91b6c83b7aa4fe545691d273cb1aee1 100644 --- a/macroscopic_value_kernels.py +++ b/macroscopic_value_kernels.py @@ -37,7 +37,6 @@ def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fiel indDims = 0 if numberOfElements <= 1 else 1 if pdfArr is None: - fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape) - 1]) outputField = Field.createGeneric(outputQuantity, lbMethod.dim, layout=fieldLayout, indexDimensions=indDims) else: outputFieldShape = pdfArr.shape[:-1] @@ -142,7 +141,9 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field return setter -def createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityRelaxationRate=1.3): +def createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velocityRelaxationRate=0.8): + + lbMethod = collisionRule.method velocityField = Field.createFromNumpyArray('velInput', velocityArray, indexDimensions=1) cqc = lbMethod.conservedQuantityComputation @@ -162,30 +163,44 @@ def createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityR # set first order relaxation rate lbMethod = deepcopy(lbMethod) lbMethod.setFirstMomentRelaxationRate(velocityRelaxationRate) + lbMethod.setZerothMomentRelaxationRate(0) + + simplificationStrategy = createSimplificationStrategy(lbMethod) + newCollisionRule = simplificationStrategy(lbMethod.getCollisionRule(eqInput)) + + # if the original collision rule used an entropy condition, the initialization should use one as well + # the simplification hints contain the entropy condition parameters + sh = collisionRule.simplificationHints + if 'entropic' in sh and sh['entropic']: + iterations = sh['entropicNewtonIterations'] + if iterations: + from lbmpy.methods.entropic import addIterativeEntropyCondition + newCollisionRule = addIterativeEntropyCondition(newCollisionRule, newtonIterations=iterations) + else: + from lbmpy.methods.entropic import addEntropyCondition + newCollisionRule = addEntropyCondition(newCollisionRule) - return lbMethod.getCollisionRule(eqInput) + return newCollisionRule -def compileAdvancedVelocitySetter(lbMethod, velocityArray, velocityRelaxationRate=1.3, pdfArr=None, - fieldLayout='numpy', optimizationParameters={}): +def compileAdvancedVelocitySetter(collisionRule, velocityArray, velocityRelaxationRate=0.8, pdfArr=None, + fieldLayout='numpy', optimizationParams={}): """ Advanced initialization of velocity field through iteration procedure according to Mei, Luo, Lallemand and Humieres: Consistent initial conditions for LBM simulations, 2005 - :param lbMethod: + :param collisionRule: :param velocityArray: array with velocity field :param velocityRelaxationRate: relaxation rate for the velocity moments - determines convergence behaviour of the initialization scheme :param pdfArr: optional numpy array for pdf field - used to get optimal loop structure for kernel :param fieldLayout: layout of the pdf field if pdfArr was not given - :param optimizationParameters: dictionary with optimization hints + :param optimizationParams: dictionary with optimization hints :return: stream-collide update function """ from lbmpy.updatekernels import createStreamPullKernel - from lbmpy.simplificationfactory import createSimplificationStrategy from lbmpy.creationfunctions import createLatticeBoltzmannAst, createLatticeBoltzmannFunction - collisionRule = createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityRelaxationRate) - simp = createSimplificationStrategy(collisionRule.method) - updateRule = createStreamPullKernel(simp(collisionRule), pdfArr, genericLayout=fieldLayout) - ast = createLatticeBoltzmannAst(updateRule, optimizationParameters) - return createLatticeBoltzmannFunction(ast, optimizationParameters) + newCollisionRule = createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velocityRelaxationRate) + updateRule = createStreamPullKernel(newCollisionRule, pdfArr, genericLayout=fieldLayout) + ast = createLatticeBoltzmannAst(updateRule, optimizationParams) + return createLatticeBoltzmannFunction(ast, optimizationParams) diff --git a/methods/entropic.py b/methods/entropic.py index d1556bbca8068e02631099d332f2b1aaaf6c7406..751c6ef592e2ca61629f6b27db40534525d7d986 100644 --- a/methods/entropic.py +++ b/methods/entropic.py @@ -50,7 +50,10 @@ def addEntropyCondition(collisionRule): index = collisionRule.method.postCollisionPdfSymbols.index(updateEq.lhs) newEq = sp.Eq(updateEq.lhs, fSymbols[index] + omega_s * dsSymbols[index] + omega_h * dhSymbols[index]) newUpdateEquations.append(newEq) - return collisionRule.copy(newUpdateEquations, collisionRule.subexpressions + subexprs) + newCollisionRule = collisionRule.copy(newUpdateEquations, collisionRule.subexpressions + subexprs) + newCollisionRule.simplificationHints['entropic'] = True + newCollisionRule.simplificationHints['entropicNewtonIterations'] = None + return newCollisionRule def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue=1): @@ -109,7 +112,10 @@ def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue # final update equations newSubexpressions = collisionRule.subexpressions + rrFactorDefinitions + fEqEqs + newtonIterationEquations - return collisionRule.copy(newUpdateEquations, newSubexpressions) + newCollisionRule = collisionRule.copy(newUpdateEquations, newSubexpressions) + newCollisionRule.simplificationHints['entropic'] = True + newCollisionRule.simplificationHints['entropicNewtonIterations'] = newtonIterations + return newCollisionRule # --------------------------------- Helper Functions and Classes ------------------------------------------------------- diff --git a/methods/momentbased.py b/methods/momentbased.py index 0a72bed4893cc0712229cd5c2edcf52cbc89a93f..e45debf40c208f14c392fedba3f19a797534a844 100644 --- a/methods/momentbased.py +++ b/methods/momentbased.py @@ -116,6 +116,12 @@ class MomentBasedLbMethod(AbstractLbMethod): True, conservedQuantityEquations) return eqColl + def setZerothMomentRelaxationRate(self, relaxationRate): + e = sp.Rational(1, 1) + prevEntry = self._momentToRelaxationInfoDict[e] + newEntry = RelaxationInfo(prevEntry[0], relaxationRate) + self._momentToRelaxationInfoDict[e] = newEntry + 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" diff --git a/serialscenario.py b/serialscenario.py index f7eddfc489f5a757a24a1e5d3b3893734036eb49..73d342e29d799c8dd4d944727593e44e40c5ce13 100644 --- a/serialscenario.py +++ b/serialscenario.py @@ -8,10 +8,10 @@ from lbmpy.boundaries import BoundaryHandling, noSlip, ubb, fixedDensity from lbmpy.stencils import getStencil -def createScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParameters, lbmKernel=None, +def createScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParams, lbmKernel=None, initialVelocity=None, preUpdateFunctions=[]): - if 'target' not in optimizationParameters: - optimizationParameters['target'] = 'cpu' + if 'target' not in optimizationParams: + optimizationParams['target'] = 'cpu' ghostLayers = 1 domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize]) @@ -26,7 +26,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza # Create kernel if lbmKernel is None: - methodParameters['optimizationParams'] = optimizationParameters + methodParameters['optimizationParams'] = optimizationParams lbmKernel = createLatticeBoltzmannFunction(**methodParameters) method = lbmKernel.method @@ -36,7 +36,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza if boundarySetupFunction is not None: symPdfField = Field.createFromNumpyArray('pdfs', pdfArrays[0], indexDimensions=1) boundaryHandling = BoundaryHandling(symPdfField, domainSize, lbmKernel.method, - target=optimizationParameters['target']) + target=optimizationParams['target']) boundarySetupFunction(boundaryHandling=boundaryHandling, method=method) boundaryHandling.prepare() else: @@ -53,7 +53,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza pdfArr=pdfArrays[0], target='cpu') setMacroscopic(pdfs=pdfArrays[0]) - if optimizationParameters['target'] == 'gpu': + if optimizationParams['target'] == 'gpu': import pycuda.gpuarray as gpuarray pdfGpuArrays = [gpuarray.to_gpu(a) for a in pdfArrays] else: @@ -95,10 +95,10 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza cpuTimeLoop.kernel = lbmKernel gpuTimeLoop.kernel = lbmKernel - return gpuTimeLoop if optimizationParameters['target'] == 'gpu' else cpuTimeLoop + return gpuTimeLoop if optimizationParams['target'] == 'gpu' else cpuTimeLoop -def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={}, lbmKernel=None, **kwargs): +def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None, **kwargs): def boundarySetupFunction(boundaryHandling, method): myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:method.dim]) myUbb.name = 'ubb' @@ -106,11 +106,11 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters= for direction in ('W', 'E', 'S') if method.dim == 2 else ('W', 'E', 'S', 'T', 'B'): boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim)) - return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel) + return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel) def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None, - lbmKernel=None, optimizationParameters={}, **kwargs): + lbmKernel=None, optimizationParams={}, **kwargs): assert dim in (2, 3) if radius is not None: @@ -152,11 +152,11 @@ def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None if 'forceModel' not in kwargs: kwargs['forceModel'] = 'guo' - return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel) + return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel) def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, lbmKernel=None, - optimizationParameters={}, initialVelocity=None, boundarySetupFunctions=[], **kwargs): + optimizationParams={}, initialVelocity=None, boundarySetupFunctions=[], **kwargs): assert dim in (2, 3) kwargs['force'] = tuple([force, 0, 0][:dim]) @@ -197,6 +197,6 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No if 'forceModel' not in kwargs: kwargs['forceModel'] = 'guo' - return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters, lbmKernel=lbmKernel, + return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel, initialVelocity=initialVelocity, preUpdateFunctions=[periodicity])