diff --git a/createmethodfromdict.py b/creationfunctions.py
similarity index 100%
rename from createmethodfromdict.py
rename to creationfunctions.py
diff --git a/macroscopicValueKernels.py b/macroscopicValueKernels.py
new file mode 100644
index 0000000000000000000000000000000000000000..679f36843f62868070e291a5581349d21427d09c
--- /dev/null
+++ b/macroscopicValueKernels.py
@@ -0,0 +1,230 @@
+import sympy as sp
+
+from pystencils.field import Field
+import pystencils.cpu as cpu
+import pystencils.gpucuda as gpu
+
+
+def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fieldLayout='numpy', target='cpu'):
+ """
+
+ :param lbMethod:
+ :param outputQuantities:
+ :param pdfArr:
+ :param fieldLayout: layout for output field, also used for pdf field if pdfArr is not given
+ :param target:
+ :return: a function to compute macroscopic values:
+ - pdfArray, can be omitted if pdf array was already passed while compiling
+ - densityOut, if not None, density is written to that array
+ - velocityOut, if not None, velocity is written to that array
+ """
+ cqc = lbMethod.conservedQuantityComputation
+ unknownQuantities = [oq for oq in outputQuantities if oq not in cqc.conservedQuantities]
+ if unknownQuantities:
+ raise ValueError("No such conserved quantity: %s, conserved quantities are %s" %
+ (str(unknownQuantities), str(cqc.conservedQuantities.keys())))
+
+ if pdfArr is None:
+ pdfField = Field.createGeneric('pdfs', lbMethod.dim, indexDimensions=1, layout=fieldLayout)
+ else:
+ pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
+
+ outputMapping = {}
+ for outputQuantity in outputQuantities:
+ numberOfElements = cqc.conservedQuantityComputation[outputQuantity]
+ assert numberOfElements >= 1
+ outputField = Field.createGeneric(outputQuantity, lbMethod.dim, layout=fieldLayout,
+ indexDimensions=0 if numberOfElements <= 1 else 1)
+
+ outputMapping[outputQuantity] = [outputField(i) for i in range(numberOfElements)]
+ if len(outputMapping[outputQuantity]) == 1:
+ outputMapping[outputQuantity] = outputMapping[outputQuantity][0]
+
+ stencil = lbMethod.stencil
+ pdfSymbols = [pdfField(i) for i in range(len(stencil))]
+ eqs = cqc.outputEquationsFromPdfs(pdfSymbols, outputMapping)
+
+ if target == 'cpu':
+ kernel = cpu.makePythonFunction(cpu.createKernel(eqs))
+ elif target == 'gpu':
+ kernel = gpu.makePythonFunction(gpu.createCUDAKernel(eqs))
+ else:
+ raise ValueError("Unknown target '%s'. Possible targets are 'cpu' and 'gpu'" % (target,))
+
+ def getter(pdfs=None, **kwargs):
+ if pdfArr is not None:
+ assert pdfs.shape == pdfArr.shape and pdfs.strides == pdfArr.strides, \
+ "Pdf array not matching blueprint which was used to compile" + str(pdfs.shape) + str(pdfArr.shape)
+ if not set(kwargs.keys()).issubset(set(outputQuantities)):
+ raise ValueError("You have to specify the output field for each of the following quantities: %s"
+ % (str(outputQuantities),))
+ kernel(pdfs=pdfs, **kwargs)
+
+ return getter
+
+
+def compileMacroscopicValuesGetter(lbMethod, pdfArr=None, macroscopicFieldLayout='numpy'):
+ """
+ Creates a function that computes density and/or velocity and stores it into given arrays
+ :param lbMethod: lattice model (to get information about stencil, force velocity shift and compressibility)
+ :param pdfArr: array to set can already be specified here, or later when the returned function is called
+ :param macroscopicFieldLayout: layout specification for Field.createGeneric
+ :return: a function, which has three parameters:
+ - pdfArray, can be omitted if pdf array was already passed while compiling
+ - densityOut, if not None, density is written to that array
+ - velocityOut, if not None, velocity is written to that array
+ """
+ if pdfArr is None:
+ pdfField = Field.createGeneric('pdfs', lbMethod.dim, indexDimensions=1)
+ else:
+ pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
+
+ rhoField = Field.createGeneric('rho', lbMethod.dim, indexDimensions=0, layout=macroscopicFieldLayout)
+ velField = Field.createGeneric('vel', lbMethod.dim, indexDimensions=1, layout=macroscopicFieldLayout)
+
+ lm = lbMethod
+ Q = len(lm.stencil)
+ symPdfs = [pdfField(i) for i in range(Q)]
+ subexpressions, rho, velSubexpressions, u = getDensityVelocityExpressions(lm.stencil, symPdfs, lm.compressible)
+
+ uRhs = [u_i.rhs for u_i in u]
+ uLhs = [u_i.lhs for u_i in u]
+ if hasattr(lm.forceModel, "macroscopicVelocity"):
+ correctedVel = lm.forceModel.macroscopicVelocity(lm, uRhs, rho.lhs)
+ u = [sp.Eq(a, b) for a, b in zip(uLhs, correctedVel)]
+
+ eqsRhoKernel = subexpressions + [sp.Eq(rhoField(0), rho.rhs)]
+ eqsVelKernel = subexpressions + [rho] + velSubexpressions + [sp.Eq(velField(i), u[i].rhs) for i in range(lm.dim)]
+ eqsRhoAndVelKernel = eqsVelKernel + [sp.Eq(rhoField(0), rho.lhs)]
+
+ stencil = lbMethod.stencil
+ cqc = lbMethod.conservedQuantityComputation
+ pdfSymbols = [pdfField(i) for i in range(len(stencil))]
+ eqsDensity = cqc.outputEquationsFromPdfs(pdfSymbols, {'density': rhoField(0)})
+ eqsVelocity = cqc.outputEquationsFromPdfs(pdfSymbols, {'velocity': [velField(i) for i in range(lbMethod.dim)]})
+
+ kernelRho = makePythonFunction(createKernel(eqsRhoKernel))
+ kernelVel = makePythonFunction(createKernel(eqsVelKernel))
+ kernelRhoAndVel = makePythonFunction(createKernel(eqsRhoAndVelKernel))
+
+ def getter(pdfs=None, densityOut=None, velocityOut=None):
+ if pdfs is not None and pdfArr is not None:
+ assert pdfs.shape == pdfArr.shape and pdfs.strides == pdfArr.strides, \
+ "Pdf array not matching blueprint which was used to compile" + str(pdfs.shape) + str(pdfArr.shape)
+ if pdfs is None:
+ assert pdfArr is not None, "Pdf array has to be passed either when compiling or when calling."
+ pdfs = pdfArr
+
+ assert not (densityOut is None and velocityOut is None), \
+ "Specify either densityOut or velocityOut parameter or both"
+
+ if densityOut is not None and velocityOut is None:
+ kernelRho(pdfs=pdfs, rho=densityOut)
+ elif densityOut is None and velocityOut is not None:
+ kernelVel(pdfs=pdfs, vel=velocityOut)
+ else:
+ kernelRhoAndVel(pdfs=pdfs, rho=densityOut, vel=velocityOut)
+
+ return getter
+
+
+def compileAdvancedVelocitySetter(collisionRule, velocityArray, pdfArr=None):
+ """
+ Advanced initialization of velocity field through iteration procedure according to
+ Mei, Luo, Lallemand and Humieres: Consistent initial conditions for LBM simulations, 2005
+
+ Important: this procedure only works if a non-zero relaxation rate was used for the velocity moments!
+
+ :param collisionRule: unsimplified collision rule
+ :param velocityArray: array with velocity field
+ :param pdfArr: optional array, to compile kernel with fixed layout and shape
+ :return: function, that has to be called multiple times, with a pdf field (src/dst) until convergence
+ similar to the normal streamCollide step, also with boundary handling
+ """
+ velocityField = Field.createFromNumpyArray('vel', velocityArray, indexDimensions=1)
+
+ # create normal LBM kernel and replace velocity by expressions of velocity field
+ from lbmpy_old.simplifications import sympyCSE
+ latticeModel = collisionRule.latticeModel
+ collisionRule = sympyCSE(collisionRule)
+ collisionRule = streamPullWithSourceAndDestinationFields(collisionRule, pdfArr)
+
+ replacements = {u_i: sp.Eq(u_i, velocityField(i)) for i, u_i in enumerate(latticeModel.symbolicVelocity)}
+
+ newSubExpressions = [replacements[eq.lhs] if eq.lhs in replacements else eq for eq in collisionRule.subexpressions]
+ newCollisionRule = LbmCollisionRule(collisionRule.updateEquations, newSubExpressions,
+ latticeModel, collisionRule.updateEquationDirections)
+ kernelAst = createKernel(newCollisionRule.equations)
+ return makePythonFunction(kernelAst, {'vel': velocityArray})
+
+
+def compileMacroscopicValuesSetter(latticeModel, density=1, velocity=0, equilibriumOrder=2, pdfArr=None):
+ """
+ Creates a function that sets a pdf field to specified macroscopic quantities
+ The returned function can be called with the pdf field to set as single argument
+
+ :param latticeModel: lattice model (to get information about stencil, force velocity shift and compressibility)
+ :param density: density used for equilibrium. Can either be scalar (all cells are set to same density) or array
+ :param velocity: velocity for equilibrium. Can either be scalar (e.g. 0, sets all cells to (0,0,0) velocity)
+ or a tuple with D (2 or 3) entries to set same velocity in the complete domain, or an array
+ specifying the velocity for each cell
+ :param equilibriumOrder: approximation order of equilibrium
+ :param pdfArr: array to set can already be specified here, or later when the returned function is called
+ """
+ if pdfArr is not None:
+ pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
+ else:
+ pdfField = Field.createGeneric('pdfs', latticeModel.dim, indexDimensions=1)
+
+ noOffset = tuple([0] * latticeModel.dim)
+ kernelArguments = {}
+
+ if hasattr(density, 'shape'):
+ densityValue = Field.createFromNumpyArray('rho', density, indexDimensions=0)[noOffset]
+ kernelArguments['rho'] = density
+ else:
+ densityValue = density
+
+ if hasattr(velocity, 'shape'):
+ assert velocity.shape[-1] == latticeModel.dim, "Wrong shape of velocity array"
+ velocityValue = [Field.createFromNumpyArray('vel', velocity, indexDimensions=1)[noOffset](i)
+ for i in range(latticeModel.dim)]
+ kernelArguments['vel'] = velocity
+ else:
+ if not hasattr(velocity, "__len__"):
+ velocity = [velocity] * latticeModel.dim
+ velocityValue = tuple(velocity)
+
+ # force shift
+ if latticeModel.forceModel and hasattr(latticeModel.forceModel, "macroscopicVelocity"):
+ # force model knows only about one direction - use sympy to solve the shift equations to get backward
+ fm = latticeModel.forceModel
+ unshiftedVel = [sp.Symbol("v_unshifted_%d" % (i,)) for i in range(latticeModel.dim)]
+ shiftedVel = fm.macroscopicVelocity(latticeModel, unshiftedVel, densityValue)
+ velShiftEqs = [sp.Eq(a, b) for a, b in zip(velocityValue, shiftedVel)]
+ solveRes = sp.solve(velShiftEqs, unshiftedVel)
+ assert len(solveRes) == latticeModel.dim
+ velocityValue = [solveRes[unshiftedVel_i] for unshiftedVel_i in unshiftedVel]
+
+ eq = standardDiscreteEquilibrium(latticeModel.stencil, densityValue, velocityValue, equilibriumOrder,
+ c_s_sq=sp.Rational(1, 3), compressible=latticeModel.compressible)
+ updateEquations = [sp.Eq(pdfField(i), eq[i]) for i in range(len(latticeModel.stencil))]
+ f = makePythonFunction(createKernel(updateEquations), kernelArguments)
+
+ def setter(pdfs=None, **kwargs):
+ if pdfs is not None and pdfArr is not None:
+ assert pdfs.shape == pdfArr.shape and pdfs.strides == pdfArr.strides, \
+ "Pdf array not matching blueprint which was used to compile"
+
+ if pdfs is None:
+ assert pdfArr is not None, "Pdf array has to be passed either when compiling or when calling."
+ pdfs = pdfArr
+ assert pdfs.shape[-1] == len(latticeModel.stencil), "Wrong shape of pdf array"
+ assert len(pdfs.shape) == latticeModel.dim + 1, "Wrong shape of pdf array"
+ if hasattr(density, 'shape'):
+ assert pdfs.shape[:-1] == density.shape, "Density array shape does not match pdf array shape"
+ if hasattr(velocity, 'shape'):
+ assert pdfs.shape[:-1] == velocity.shape[:-1], "Velocity array shape does not match pdf array shape"
+ f(pdfs=pdfs, **kwargs) # kwargs passed here, since force model might need additional fields
+
+ return setter
diff --git a/methods/__init__.py b/methods/__init__.py
index 230ec443b76472993d249870a6d9332e43228ef0..e187d0ab43733ea92530fd0504937c1fe639a0ee 100644
--- a/methods/__init__.py
+++ b/methods/__init__.py
@@ -1,5 +1,5 @@
-from lbmpy.methods.abstractlbmmethod import AbstractLbmMethod
-from lbmpy.methods.momentbased import MomentBasedLbmMethod, RelaxationInfo
+from lbmpy.methods.abstractlbmmethod import AbstractLbMethod
+from lbmpy.methods.momentbased import MomentBasedLbMethod, RelaxationInfo
from lbmpy.methods.momentbased import createSRT, createTRT, createTRTWithMagicNumber, createOrthogonalMRT, \
createWithContinuousMaxwellianEqMoments, createWithDiscreteMaxwellianEqMoments
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
diff --git a/methods/abstractlbmmethod.py b/methods/abstractlbmethod.py
similarity index 97%
rename from methods/abstractlbmmethod.py
rename to methods/abstractlbmethod.py
index 8cd720ab08d3fb92e4575a71020b4b126a07b3e6..201ab0547a44c4b9b8f60c0410a5476c36dd56b7 100644
--- a/methods/abstractlbmmethod.py
+++ b/methods/abstractlbmethod.py
@@ -9,7 +9,7 @@ class LbmCollisionRule(EquationCollection):
self.method = lbmMethod
-class AbstractLbmMethod(metaclass=abc.ABCMeta):
+class AbstractLbMethod(metaclass=abc.ABCMeta):
"""
Abstract base class for all LBM methods
"""
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index a65d2a46804c73f6c5e444be135db7935f94172a..5e07c2aa91d74420bcc90a01ed9e6e9696213e0e 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -53,14 +53,14 @@ class AbstractConservedQuantityComputation(metaclass=abc.ABCMeta):
"""
@abc.abstractmethod
- def outputEquationsFromPdfs(self, pdfs, outputQuantityNames):
+ def outputEquationsFromPdfs(self, pdfs, outputQuantityNamesToSymbols):
"""
Returns an equation collection that defines conserved quantities for output. These conserved quantities might
be slightly different that the ones used as input for the equilibrium e.g. due to a force model.
:param pdfs: values for the pdf entries
- :param outputQuantityNames: list of conserved quantity names, defining which parameters should be written out.
- See :func:`conservedQuantities`
+ :param outputQuantityNamesToSymbols: dict mapping of conserved quantity names (See :func:`conservedQuantities`)
+ to symbols or field accesses where they should be written to
"""
@abc.abstractmethod
diff --git a/methods/momentbased.py b/methods/momentbased.py
index c5fe7d97acf6a56f54bd8412e8b6cb6859795d80..198ed4f3c31d1afbe3e6334eca763f8e74a3de61 100644
--- a/methods/momentbased.py
+++ b/methods/momentbased.py
@@ -5,7 +5,7 @@ from collections import namedtuple, OrderedDict, defaultdict
from lbmpy.stencils import stencilsHaveSameEntries, getStencil
from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
getMomentsOfContinuousMaxwellianEquilibrium
-from lbmpy.methods.abstractlbmmethod import AbstractLbmMethod, LbmCollisionRule
+from lbmpy.methods.abstractlbmmethod import AbstractLbMethod, LbmCollisionRule
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
from lbmpy.moments import MOMENT_SYMBOLS, momentMatrix, isShearMoment, \
isEven, gramSchmidt, getOrder, getDefaultMomentSetForStencil
@@ -14,7 +14,7 @@ from pystencils.sympyextensions import commonDenominator, replaceAdditive
RelaxationInfo = namedtuple('Relaxationinfo', ['equilibriumValue', 'relaxationRate'])
-class MomentBasedLbmMethod(AbstractLbmMethod):
+class MomentBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, momentToRelaxationInfoDict, conservedQuantityComputation, forceModel=None):
"""
Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
@@ -31,7 +31,7 @@ class MomentBasedLbmMethod(AbstractLbmMethod):
the symbols used in the equilibrium moments like e.g. density and velocity
:param forceModel: force model instance, or None if no forcing terms are required
"""
- super(MomentBasedLbmMethod, self).__init__(stencil)
+ super(MomentBasedLbMethod, self).__init__(stencil)
assert isinstance(conservedQuantityComputation, AbstractConservedQuantityComputation)
@@ -272,7 +272,7 @@ def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, c
compressible=compressible, order=equilibriumAccuracyOrder)
rrDict = OrderedDict([(mom, RelaxationInfo(eqMom, rr))
for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqMoments)])
- return MomentBasedLbmMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+ return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict, forceModel=None,
@@ -292,7 +292,7 @@ def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict,
order=equilibriumAccuracyOrder)
rrDict = OrderedDict([(mom, RelaxationInfo(eqMom, rr))
for mom, rr, eqMom in zip(momToRrDict.keys(), momToRrDict.values(), eqMoments)])
- return MomentBasedLbmMethod(stencil, rrDict, densityVelocityComputation, forceModel)
+ return MomentBasedLbMethod(stencil, rrDict, densityVelocityComputation, forceModel)
# ------------------------------------ SRT / TRT/ MRT Creators ---------------------------------------------------------
diff --git a/plot2d.py b/plot2d.py
new file mode 100644
index 0000000000000000000000000000000000000000..b1bb364721d143a77814256105d9cabdf847db2f
--- /dev/null
+++ b/plot2d.py
@@ -0,0 +1,47 @@
+from matplotlib.pyplot import *
+
+
+def removeGhostLayers(field):
+ return field[1:-1, 1:-1]
+
+
+def vectorField(field, step=2, **kwargs):
+ field = removeGhostLayers(field)
+ veln = field.swapaxes(0, 1)
+ quiver(veln[::step, ::step, 0], veln[::step, ::step, 1], **kwargs)
+
+
+def vectorFieldMagnitude(field, **kwargs):
+ from numpy.linalg import norm
+ field = norm(field, axis=2, ord=2)
+ return scalarField(field, **kwargs)
+
+
+def scalarField(field, **kwargs):
+ import numpy as np
+ field = removeGhostLayers(field)
+ field = np.swapaxes(field, 0, 1)
+ return imshow(field, origin='lower', **kwargs)
+
+
+def vectorFieldMagnitudeAnimation(runFunction, plotSetupFunction=lambda: None,
+ plotUpdateFunction=lambda: None, interval=30, frames=180, **kwargs):
+ import matplotlib.animation as animation
+ import numpy as np
+ from numpy.linalg import norm
+
+ fig = gcf()
+ im = None
+ field = runFunction()
+ im = vectorFieldMagnitude(field, **kwargs)
+ plotSetupFunction()
+
+ def updatefig(*args):
+ f = runFunction()
+ f = norm(f, axis=2, ord=2)
+ f = np.swapaxes(f, 0, 1)
+ im.set_array(f)
+ plotUpdateFunction()
+ return im,
+
+ return animation.FuncAnimation(fig, updatefig, interval=interval, frames=frames)
diff --git a/simplificationfactory.py b/simplificationfactory.py
index b6f42ac057ddadc4696799db0d57205938846bbd..ba8962f9ee1cb164d7a8c8c67a0bdab49b111cbd 100644
--- a/simplificationfactory.py
+++ b/simplificationfactory.py
@@ -8,14 +8,14 @@ from pystencils.equationcollection.simplifications import applyOnAllEquations, \
def createSimplificationStrategy(lbmMethod, doCseInOpposingDirections=False, doOverallCse=False, splitInnerLoop=False):
from pystencils.equationcollection import SimplificationStrategy
- from lbmpy.methods import MomentBasedLbmMethod
+ from lbmpy.methods import MomentBasedLbMethod
from lbmpy.methods.momentbasedsimplifications import replaceSecondOrderVelocityProducts, \
factorDensityAfterFactoringRelaxationTimes, factorRelaxationRates, cseInOpposingDirections, \
replaceCommonQuadraticAndConstantTerm, replaceDensityAndVelocity
s = SimplificationStrategy()
- if isinstance(lbmMethod, MomentBasedLbmMethod):
+ if isinstance(lbmMethod, MomentBasedLbMethod):
expand = partial(applyOnAllEquations, operation=sp.expand)
expand.__name__ = "expand"