From 8bd55fe316b6c69d54830cd593d1bcd7378fffdf Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Mon, 30 Jan 2017 09:53:24 +0100
Subject: [PATCH] fixed compressible bug
---
boundaries/boundaryconditions.py | 5 ++++
creationfunctions.py | 8 ++++++-
methods/conservedquantitycomputation.py | 31 ++++++++++++++++++++-----
serialscenario.py | 11 +++++----
4 files changed, 43 insertions(+), 12 deletions(-)
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 35b3709e..9650abdc 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -20,6 +20,11 @@ def ubb(pdfField, direction, lbMethod, velocity):
neighbor = offsetFromDir(direction, lbMethod.dim)
inverseDir = invDir(direction)
+ # TODO adapt velocity to force
+ # TODO compute density
+
+ densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
+
velTerm = 6 * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
return [sp.Eq(pdfField[neighbor](inverseDir),
pdfField(direction) - velTerm)]
diff --git a/creationfunctions.py b/creationfunctions.py
index 1409f05b..89e59848 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -26,7 +26,7 @@ def _getParams(params, optParams):
defaultOptimizationDescription = {
'doCseInOpposingDirections': False,
'doOverallCse': False,
- 'split': True,
+ 'split': False,
'fieldSize': None,
'fieldLayout': 'reverseNumpy', # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
@@ -80,6 +80,7 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
if params['target'] == 'cpu':
from pystencils.cpu import createKernel
if 'splitGroups' in updateRule.simplificationHints:
+ print("splitting!")
splitGroups = updateRule.simplificationHints['splitGroups']
else:
splitGroups = ()
@@ -91,6 +92,11 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
return ValueError("'target' has to be either 'cpu' or 'gpu'")
res.method = updateRule.method
+ #TODO debug begin
+ from pystencils.cpu import generateC
+ from pystencils.display_utils import highlightCpp
+ print(generateC(res))
+ #TODO debug end
return res
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 5e07c2aa..2cd989ba 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -1,4 +1,6 @@
import abc
+from collections import OrderedDict
+
import sympy as sp
from pystencils.equationcollection import EquationCollection
@@ -139,8 +141,6 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
def outputEquationsFromPdfs(self, pdfs, outputQuantityNamesToSymbols):
dim = len(self._stencil[0])
- symbolsToExtract = set()
-
eqColl = getEquationsForZerothAndFirstOrderMoment(self._stencil, pdfs, self._symbolOrder0, self._symbolsOrder1)
if self._compressible:
@@ -151,14 +151,22 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
eqColl = applyForceModelShift('macroscopicVelocityShift', dim, eqColl, self._forceModel, self._compressible)
mainEquations = []
+ eqs = OrderedDict([(eq.lhs, eq.rhs) for eq in eqColl.allEquations])
+
if 'density' in outputQuantityNamesToSymbols:
densityOutputSymbol = outputQuantityNamesToSymbols['density']
- mainEquations.append(sp.Eq(densityOutputSymbol, self._symbolOrder0))
- symbolsToExtract.add(densityOutputSymbol)
+ if densityOutputSymbol != self._symbolOrder0:
+ mainEquations.append(sp.Eq(densityOutputSymbol, self._symbolOrder0))
+ else:
+ mainEquations.append(sp.Eq(self._symbolOrder0, eqs[self._symbolOrder0]))
+ del eqs[self._symbolOrder0]
if 'velocity' in outputQuantityNamesToSymbols:
velOutputSymbols = outputQuantityNamesToSymbols['velocity']
- mainEquations += [sp.Eq(a, b) for a, b in zip(velOutputSymbols, self._symbolsOrder1)]
- symbolsToExtract.update(velOutputSymbols)
+ if tuple(velOutputSymbols) != tuple(self._symbolsOrder1):
+ mainEquations += [sp.Eq(a, b) for a, b in zip(velOutputSymbols, self._symbolsOrder1)]
+ else:
+ # TODO
+ pass
eqColl = eqColl.copy(mainEquations, eqColl.allEquations)
return eqColl.newWithoutUnusedSubexpressions()
@@ -271,3 +279,14 @@ def applyForceModelShift(shiftMemberName, dim, equationCollection, forceModel, c
+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
diff --git a/serialscenario.py b/serialscenario.py
index 1bb25fc8..96fd5fac 100644
--- a/serialscenario.py
+++ b/serialscenario.py
@@ -18,6 +18,7 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
# 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)
@@ -38,7 +39,7 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
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 = compileMacroscopicValuesSetter(method, {'density': 1.0, 'velocity': [0] * D}, pdfArr=pdfSrc)
setMacroscopic(pdfs=pdfSrc)
# Run simulation
@@ -50,6 +51,7 @@ def runScenario(domainSize, boundarySetupFunction, methodParameters, optimizatio
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
@@ -64,6 +66,7 @@ def runLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParameters={},
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)
@@ -74,11 +77,11 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
if radius is not None:
assert length is not None
if dim == 3:
- domainSize = (length, 2*radius+1, 2*radius+1)
+ domainSize = (length, 2 * radius + 1, 2 * radius + 1)
roundChannel = True
else:
if domainSize is None:
- domainSize = (length, 2*radius)
+ domainSize = (length, 2 * radius)
else:
roundChannel = False
@@ -108,5 +111,3 @@ def runForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None,
return runScenario(domainSize, boundarySetupFunction, kwargs, optimizationParameters,
preUpdateFunctions=[periodicity])
-
-
--
GitLab