From ad5a4d9bf84cf6ff0ce9e5e996a10e840df39fd0 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Wed, 26 Jul 2017 13:10:47 +0200
Subject: [PATCH] lbmpy: new parameters for input and output of macroscopic
quantities
- write out macroscopic quantities to field, directly from sweep
- use velocity from field ( for Advection-Diffusion schemes )
---
creationfunctions.py | 24 +++++++++++++++++++++++-
methods/conservedquantitycomputation.py | 17 ++++++++++++++---
methods/creationfunctions.py | 2 +-
3 files changed, 38 insertions(+), 5 deletions(-)
diff --git a/creationfunctions.py b/creationfunctions.py
index 35b5cea7..c6152bca 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -50,6 +50,9 @@ General:
- ``initialVelocity=None``: initial velocity in domain, can either be a tuple (x,y,z) velocity to set a constant
velocity everywhere, or a numpy array with the same size of the domain, with a last coordinate of shape dim to set
velocities on cell level
+- ``output={}``: a dictionary mapping macroscopic quantites e.g. the strings 'density' and 'velocity' to pystencils
+ fields. In each timestep the corresponding quantities are written to the given fields.
+- ``velocityInput``: symbolic field where the velocities are read from (for advection diffusion LBM)
Entropic methods:
@@ -151,6 +154,7 @@ from lbmpy.stencils import getStencil
import lbmpy.forcemodels as forceModels
from lbmpy.simplificationfactory import createSimplificationStrategy
from lbmpy.updatekernels import createStreamPullKernel, createPdfArray
+from pystencils.equationcollection.equationcollection import EquationCollection
def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
@@ -171,6 +175,9 @@ def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
'entropic': False,
'entropicNewtonIterations': None,
'omegaOutputField': None,
+
+ 'output': {},
+ 'velocityInput': None,
}
defaultOptimizationDescription = {
@@ -305,7 +312,22 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
doCseInOpposingDirections = False if dirCSE not in optParams else optParams[dirCSE]
doOverallCse = False if 'doOverallCse' not in optParams else optParams['doOverallCse']
simplification = createSimplificationStrategy(lbMethod, doCseInOpposingDirections, doOverallCse, splitInnerLoop)
- collisionRule = simplification(lbMethod.getCollisionRule())
+ cqc = lbMethod.conservedQuantityComputation
+
+ if params['velocityInput'] is not None:
+ eqs = [sp.Eq(cqc.zerothOrderMomentSymbol, sum(lbMethod.preCollisionPdfSymbols))]
+ velocityField = params['velocityInput']
+ eqs += [sp.Eq(uSym, velocityField(i)) for i, uSym in enumerate(cqc.firstOrderMomentSymbols)]
+ eqs = EquationCollection(eqs, [])
+ collisionRule = lbMethod.getCollisionRule(conservedQuantityEquations=eqs)
+ else:
+ collisionRule = lbMethod.getCollisionRule()
+
+ if params['output']:
+ outputEqs = cqc.outputEquationsFromPdfs(lbMethod.preCollisionPdfSymbols, params['output'])
+ collisionRule = collisionRule.merge(outputEqs)
+
+ collisionRule = simplification(collisionRule)
if params['entropic']:
if params['entropicNewtonIterations']:
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 60ff2c6d..e6c04922 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -3,6 +3,7 @@ from collections import OrderedDict
import sympy as sp
from pystencils.equationcollection import EquationCollection
+from pystencils.field import Field
class AbstractConservedQuantityComputation(abc.ABCMeta('ABC', (object,), {})):
@@ -132,11 +133,12 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
eqColl = applyForceModelShift('equilibriumVelocityShift', dim, eqColl, self._forceModel, self._compressible)
return eqColl
- def equilibriumInputEquationsFromInitValues(self, density=1, velocity=[0, 0, 0]):
+ def equilibriumInputEquationsFromInitValues(self, density=1, velocity=(0, 0, 0)):
dim = len(self._stencil[0])
zerothOrderMoment = density
firstOrderMoments = velocity[:dim]
velOffset = [0] * dim
+
if self._compressible:
if self._forceModel and hasattr(self._forceModel, 'macroscopicVelocityShift'):
velOffset = self._forceModel.macroscopicVelocityShift(zerothOrderMoment)
@@ -144,10 +146,14 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
if self._forceModel and hasattr(self._forceModel, 'macroscopicVelocityShift'):
velOffset = self._forceModel.macroscopicVelocityShift(sp.Rational(1, 1))
zerothOrderMoment -= sp.Rational(1, 1)
+ eqs = [sp.Eq(self._symbolOrder0, zerothOrderMoment)]
+
firstOrderMoments = [a - b for a, b in zip(firstOrderMoments, velOffset)]
+ if self._compressible:
+ eqs += [sp.expand(sp.Eq(l, r * density)) for l, r in zip(self._symbolsOrder1, firstOrderMoments)]
+ else:
+ eqs += [sp.Eq(l, r) for l, r in zip(self._symbolsOrder1, firstOrderMoments)]
- eqs = [sp.Eq(self._symbolOrder0, zerothOrderMoment)]
- eqs += [sp.Eq(l, r) for l, r in zip(self._symbolsOrder1, firstOrderMoments)]
return EquationCollection(eqs, [])
def outputEquationsFromPdfs(self, pdfs, outputQuantityNamesToSymbols):
@@ -167,6 +173,8 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
if 'density' in outputQuantityNamesToSymbols:
densityOutputSymbol = outputQuantityNamesToSymbols['density']
+ if isinstance(densityOutputSymbol, Field):
+ densityOutputSymbol = densityOutputSymbol()
if densityOutputSymbol != self._symbolOrder0:
mainEquations.append(sp.Eq(densityOutputSymbol, self._symbolOrder0))
else:
@@ -174,6 +182,9 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
del eqs[self._symbolOrder0]
if 'velocity' in outputQuantityNamesToSymbols:
velOutputSymbols = outputQuantityNamesToSymbols['velocity']
+ if isinstance(velOutputSymbols, Field):
+ field = velOutputSymbols
+ velOutputSymbols = [field(i) for i in range(len(self._symbolsOrder1))]
if tuple(velOutputSymbols) != tuple(self._symbolsOrder1):
mainEquations += [sp.Eq(a, b) for a, b in zip(velOutputSymbols, self._symbolsOrder1)]
else:
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 9b5e2bfb..0e8d0647 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -129,7 +129,7 @@ def createFromEquilibrium(stencil, equilibrium, momentToRelaxationRateDict, comp
assert len(momToRrDict) == len(stencil), "The number of moments has to be the same as the number of stencil entries"
densityVelocityComputation = DensityVelocityComputation(stencil, compressible, forceModel)
- rrDict = OrderedDict([(mom, RelaxationInfo(discreteMoment(equilibrium, mom, stencil), rr))
+ rrDict = OrderedDict([(mom, RelaxationInfo(discreteMoment(equilibrium, mom, stencil).expand(), rr))
for mom, rr in zip(momToRrDict.keys(), momToRrDict.values())])
return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
--
GitLab