Commit 4d7380a8 authored by Martin Bauer's avatar Martin Bauer
Browse files

LB creation functions: more flexibility

- now an lbMethod or updateRule can be passed in directly
- made method/update rule/ast function independent i.e. no function
  uses the same parameters as on of the others
parent 70f4d7ad
......@@ -153,125 +153,16 @@ from lbmpy.methods.entropic import addIterativeEntropyCondition, addEntropyCondi
from lbmpy.methods.entropic_eq_srt import createEntropicSRT
from lbmpy.methods.relaxationrates import relaxationRateFromMagicNumber
from lbmpy.stencils import getStencil, stencilsHaveSameEntries
import lbmpy.forcemodels as forceModels
import lbmpy.forcemodels as forcemodels
from lbmpy.simplificationfactory import createSimplificationStrategy
from lbmpy.updatekernels import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAccessor, CollideOnlyInplaceAccessor, \
createLBMKernel
from pystencils.data_types import collateTypes
from pystencils.equationcollection.equationcollection import EquationCollection
from pystencils.field import getLayoutOfArray, Field
from pystencils import createKernel
def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
defaultMethodDescription = {
'stencil': 'D2Q9',
'method': 'srt', # can be srt, trt or mrt
'relaxationRates': None,
'compressible': False,
'equilibriumAccuracyOrder': 2,
'c_s_sq': sp.Rational(1, 3),
'forceModel': 'none',
'force': (0, 0, 0),
'useContinuousMaxwellianEquilibrium': True,
'cumulant': False,
'initialVelocity': None,
'entropic': False,
'entropicNewtonIterations': None,
'omegaOutputField': None,
'output': {},
'velocityInput': None,
'kernelType': 'streamPullCollide',
'fieldName': 'src',
'secondFieldName': 'dst',
}
defaultOptimizationDescription = {
'doCseInOpposingDirections': False,
'doOverallCse': False,
'split': False,
'fieldSize': None,
'fieldLayout': 'fzyx', # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
'symbolicField': None,
'target': 'cpu',
'openMP': False,
'doublePrecision': True,
'gpuIndexing': 'block',
'gpuIndexingParams': {},
'vectorization': None,
'builtinPeriodicity': (False, False, False),
}
if 'relaxationRate' in params:
if 'relaxationRates' not in params:
if 'entropic' in params and params['entropic']:
params['relaxationRates'] = [params['relaxationRate']]
else:
params['relaxationRates'] = [params['relaxationRate'],
relaxationRateFromMagicNumber(params['relaxationRate'])]
del params['relaxationRate']
if 'pdfArr' in optParams:
arr = optParams['pdfArr']
optParams['fieldSize'] = tuple(e - 2 for e in arr.shape[:-1])
optParams['fieldLayout'] = getLayoutOfArray(arr)
del optParams['pdfArr']
if failOnUnknownParameter:
unknownParams = [k for k in params.keys() if k not in defaultMethodDescription]
unknownOptParams = [k for k in optParams.keys() if k not in defaultOptimizationDescription]
if unknownParams:
raise ValueError("Unknown parameter(s): " + ", ".join(unknownParams))
if unknownOptParams:
raise ValueError("Unknown optimization parameter(s): " + ",".join(unknownOptParams))
paramsResult = copy(defaultMethodDescription)
paramsResult.update(params)
optParamsResult = copy(defaultOptimizationDescription)
optParamsResult.update(optParams)
if paramsResult['relaxationRates'] is None:
stencilParam = paramsResult['stencil']
if isinstance(stencilParam, tuple) or isinstance(stencilParam, list):
stencilEntries = stencilParam
else:
stencilEntries = getStencil(paramsResult['stencil'])
paramsResult['relaxationRates'] = sp.symbols("omega_:%d" % len(stencilEntries))
return paramsResult, optParamsResult
def switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams):
"""
For entropic kernels the relaxation rate has to be a variable. If a constant was passed a
new dummy variable is inserted and the value of this variable is later on passed to the kernel
"""
if methodParameters['entropic']:
valueToSymbolMap = {}
newRelaxationRates = []
for rr in methodParameters['relaxationRates']:
if not isinstance(rr, sp.Symbol):
if rr not in valueToSymbolMap:
valueToSymbolMap[rr] = sp.Dummy()
dummyVar = valueToSymbolMap[rr]
newRelaxationRates.append(dummyVar)
kernelParams[dummyVar.name] = rr
else:
newRelaxationRates.append(rr)
if len(newRelaxationRates) < 2:
newRelaxationRates.append(sp.Dummy())
methodParameters['relaxationRates'] = newRelaxationRates
def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
params, optParams = updateWithDefaultParameters(kwargs, optimizationParams)
......@@ -294,8 +185,8 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
params['optimizationParams'] = optimizationParams
updateRule = createLatticeBoltzmannUpdateRule(**params)
dtype = 'double' if optParams['doublePrecision'] else 'float32'
res = createKernel(updateRule, target=optParams['target'], dataType=dtype,
fieldTypes = set(fa.field.dtype for fa in updateRule.freeSymbols if isinstance(fa, Field.Access))
res = createKernel(updateRule, target=optParams['target'], dataType=collateTypes(fieldTypes),
cpuOpenMP=optParams['openMP'], cpuVectorizeInfo=optParams['vectorization'],
gpuIndexing=optParams['gpuIndexing'], gpuIndexingParams=optParams['gpuIndexingParams'],
ghostLayers=1)
......@@ -306,11 +197,9 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
@diskcacheNoFallback
def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwargs):
def createLatticeBoltzmannCollisionRule(lbMethod=None, optimizationParams={}, **kwargs):
params, optParams = updateWithDefaultParameters(kwargs, optimizationParams)
stencil = getStencil(params['stencil'])
if lbMethod is None:
lbMethod = createLatticeBoltzmannMethod(**params)
......@@ -335,11 +224,19 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
outputEqs = cqc.outputEquationsFromPdfs(lbMethod.preCollisionPdfSymbols, params['output'])
collisionRule = collisionRule.merge(outputEqs)
collisionRule = simplification(collisionRule)
return simplification(collisionRule)
@diskcacheNoFallback
def createLatticeBoltzmannUpdateRule(collisionRule=None, optimizationParams={}, **kwargs):
params, optParams = updateWithDefaultParameters(kwargs, optimizationParams)
if collisionRule is None:
collisionRule = createLatticeBoltzmannCollisionRule(**params, optimizationParams=optParams)
if params['entropic']:
if params['entropicNewtonIterations']:
if isinstance(params['entropicNewtonIterations'], bool) or params['cumulant']:
if isinstance(params['entropicNewtonIterations'], bool):
iterations = 3
else:
iterations = params['entropicNewtonIterations']
......@@ -348,20 +245,17 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
else:
collisionRule = addEntropyCondition(collisionRule, omegaOutputField=params['omegaOutputField'])
if 'doublePrecision' in optimizationParams:
fieldDtype = 'float64' if optimizationParams['doublePrecision'] else 'float32'
else:
fieldDtype = 'float64'
fieldDtype = 'float64' if optParams['doublePrecision'] else 'float32'
if optParams['symbolicField'] is not None:
srcField = optParams['symbolicField']
elif optParams['fieldSize']:
fieldSize = [s + 2 for s in optParams['fieldSize']] + [len(stencil)]
fieldSize = [s + 2 for s in optParams['fieldSize']] + [len(collisionRule.stencil)]
srcField = Field.createFixedSize(params['fieldName'], fieldSize, indexDimensions=1,
layout=optParams['fieldLayout'], dtype=fieldDtype)
else:
srcField = Field.createGeneric(params['fieldName'], spatialDimensions=lbMethod.dim, indexDimensions=1,
layout=optParams['fieldLayout'], dtype=fieldDtype)
srcField = Field.createGeneric(params['fieldName'], spatialDimensions=collisionRule.method.dim,
indexDimensions=1, layout=optParams['fieldLayout'], dtype=fieldDtype)
dstField = srcField.newFieldWithDifferentName(params['secondFieldName'])
......@@ -399,21 +293,7 @@ def createLatticeBoltzmannMethod(**params):
params['forceModel'] = 'guo'
if 'forceModel' in params:
forceModelName = params['forceModel']
if type(forceModelName) is not str:
forceModel = forceModelName
elif forceModelName.lower() == 'none':
forceModel = None
elif forceModelName.lower() == 'simple':
forceModel = forceModels.Simple(params['force'][:dim])
elif forceModelName.lower() == 'luo':
forceModel = forceModels.Luo(params['force'][:dim])
elif forceModelName.lower() == 'guo':
forceModel = forceModels.Guo(params['force'][:dim])
elif forceModelName.lower() == 'silva' or forceModelName.lower() == 'buick':
forceModel = forceModels.Buick(params['force'][:dim])
else:
raise ValueError("Unknown force model %s" % (forceModelName,))
forceModel = forceModelFromString(params['forceModel'], params['force'][:dim])
else:
forceModel = None
......@@ -462,3 +342,138 @@ def createLatticeBoltzmannMethod(**params):
return method
# ----------------------------------------------------------------------------------------------------------------------
def forceModelFromString(forceModelName, forceValues):
if type(forceModelName) is not str:
forceModel = forceModelName
elif forceModelName.lower() == 'none':
forceModel = None
elif forceModelName.lower() == 'simple':
forceModel = forcemodels.Simple(forceValues)
elif forceModelName.lower() == 'luo':
forceModel = forcemodels.Luo(forceValues)
elif forceModelName.lower() == 'guo':
forceModel = forcemodels.Guo(forceValues)
elif forceModelName.lower() == 'silva' or forceModelName.lower() == 'buick':
forceModel = forcemodels.Buick(forceValues)
else:
raise ValueError("Unknown force model %s" % (forceModelName,))
return forceModel
def switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams):
"""
For entropic kernels the relaxation rate has to be a variable. If a constant was passed a
new dummy variable is inserted and the value of this variable is later on passed to the kernel
"""
if methodParameters['entropic']:
valueToSymbolMap = {}
newRelaxationRates = []
for rr in methodParameters['relaxationRates']:
if not isinstance(rr, sp.Symbol):
if rr not in valueToSymbolMap:
valueToSymbolMap[rr] = sp.Dummy()
dummyVar = valueToSymbolMap[rr]
newRelaxationRates.append(dummyVar)
kernelParams[dummyVar.name] = rr
else:
newRelaxationRates.append(rr)
if len(newRelaxationRates) < 2:
newRelaxationRates.append(sp.Dummy())
methodParameters['relaxationRates'] = newRelaxationRates
def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
defaultMethodDescription = {
'stencil': 'D2Q9',
'method': 'srt', # can be srt, trt or mrt
'relaxationRates': None,
'compressible': False,
'equilibriumAccuracyOrder': 2,
'c_s_sq': sp.Rational(1, 3),
'forceModel': 'none',
'force': (0, 0, 0),
'useContinuousMaxwellianEquilibrium': True,
'cumulant': False,
'initialVelocity': None,
'entropic': False,
'entropicNewtonIterations': None,
'omegaOutputField': None,
'output': {},
'velocityInput': None,
'kernelType': 'streamPullCollide',
'fieldName': 'src',
'secondFieldName': 'dst',
'lbMethod': None,
'collisionRule': None,
'updateRule': None,
'ast': None,
}
defaultOptimizationDescription = {
'doCseInOpposingDirections': False,
'doOverallCse': False,
'split': False,
'fieldSize': None,
'fieldLayout': 'fzyx', # can be 'numpy' (='c'), 'reverseNumpy' (='f'), 'fzyx', 'zyxf'
'symbolicField': None,
'target': 'cpu',
'openMP': False,
'doublePrecision': True,
'gpuIndexing': 'block',
'gpuIndexingParams': {},
'vectorization': None,
'builtinPeriodicity': (False, False, False),
}
if 'relaxationRate' in params:
if 'relaxationRates' not in params:
if 'entropic' in params and params['entropic']:
params['relaxationRates'] = [params['relaxationRate']]
else:
params['relaxationRates'] = [params['relaxationRate'],
relaxationRateFromMagicNumber(params['relaxationRate'])]
del params['relaxationRate']
if 'pdfArr' in optParams:
arr = optParams['pdfArr']
optParams['fieldSize'] = tuple(e - 2 for e in arr.shape[:-1])
optParams['fieldLayout'] = getLayoutOfArray(arr)
del optParams['pdfArr']
if failOnUnknownParameter:
unknownParams = [k for k in params.keys() if k not in defaultMethodDescription]
unknownOptParams = [k for k in optParams.keys() if k not in defaultOptimizationDescription]
if unknownParams:
raise ValueError("Unknown parameter(s): " + ", ".join(unknownParams))
if unknownOptParams:
raise ValueError("Unknown optimization parameter(s): " + ",".join(unknownOptParams))
paramsResult = copy(defaultMethodDescription)
paramsResult.update(params)
optParamsResult = copy(defaultOptimizationDescription)
optParamsResult.update(optParams)
if paramsResult['relaxationRates'] is None:
stencilParam = paramsResult['stencil']
if isinstance(stencilParam, tuple) or isinstance(stencilParam, list):
stencilEntries = stencilParam
else:
stencilEntries = getStencil(paramsResult['stencil'])
paramsResult['relaxationRates'] = sp.symbols("omega_:%d" % len(stencilEntries))
return paramsResult, optParamsResult
\ No newline at end of file
......@@ -40,6 +40,8 @@ def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, c
:param c_s_sq: Speed of sound squared
:return: :class:`lbmpy.methods.MomentBasedLbMethod` instance
"""
if isinstance(stencil, str):
stencil = getStencil(stencil)
momToRrDict = OrderedDict(momentToRelaxationRateDict)
assert len(momToRrDict) == len(stencil), \
"The number of moments has to be the same as the number of stencil entries"
......@@ -72,6 +74,8 @@ def createWithContinuousMaxwellianEqMoments(stencil, momentToRelaxationRateDict,
For parameter description see :func:`lbmpy.methods.createWithDiscreteMaxwellianEqMoments`.
By using the continuous Maxwellian we automatically get a compressible model.
"""
if isinstance(stencil, str):
stencil = getStencil(stencil)
momToRrDict = OrderedDict(momentToRelaxationRateDict)
assert len(momToRrDict) == len(stencil), "The number of moments has to be the same as the number of stencil entries"
dim = len(stencil[0])
......@@ -125,10 +129,16 @@ def createFromEquilibrium(stencil, equilibrium, momentToRelaxationRateDict, comp
Creates a moment-based LB method using a given equilibrium distribution function
:param stencil: see createWithDiscreteMaxwellianEqMoments
:param equilibrium: list of equilibrium terms, dependent on rho and u, one for each stencil direction
:param momentToRelaxationRateDict: relaxation rate for each moment
:param momentToRelaxationRateDict: relaxation rate for each moment, or a symbol/float if all should relaxed with the
same rate
:param compressible: see createWithDiscreteMaxwellianEqMoments
:param forceModel: see createWithDiscreteMaxwellianEqMoments
"""
if isinstance(stencil, str):
stencil = getStencil(stencil)
if isinstance(momentToRelaxationRateDict, sp.Symbol) or isinstance(momentToRelaxationRateDict, float):
momentToRelaxationRateDict = {m: momentToRelaxationRateDict for m in getDefaultMomentSetForStencil(stencil)}
momToRrDict = OrderedDict(momentToRelaxationRateDict)
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)
......@@ -152,6 +162,8 @@ def createSRT(stencil, relaxationRate, useContinuousMaxwellianEquilibrium=False,
used to compute the equilibrium moments
:return: :class:`lbmpy.methods.MomentBasedLbMethod` instance
"""
if isinstance(stencil, str):
stencil = getStencil(stencil)
moments = getDefaultMomentSetForStencil(stencil)
rrDict = OrderedDict([(m, relaxationRate) for m in moments])
if useContinuousMaxwellianEquilibrium:
......@@ -171,6 +183,8 @@ def createTRT(stencil, relaxationRateEvenMoments, relaxationRateOddMoments,
two relaxation rates: one for even moments (determines viscosity) and one for odd moments.
If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.createTRTWithMagicNumber`.
"""
if isinstance(stencil, str):
stencil = getStencil(stencil)
moments = getDefaultMomentSetForStencil(stencil)
rrDict = OrderedDict([(m, relaxationRateEvenMoments if isEven(m) else relaxationRateOddMoments) for m in moments])
if useContinuousMaxwellianEquilibrium:
......@@ -193,6 +207,8 @@ def createRawMRT(stencil, relaxationRates, useContinuousMaxwellianEquilibrium=Fa
"""
Creates a MRT method using non-orthogonalized moments
"""
if isinstance(stencil, str):
stencil = getStencil(stencil)
moments = getDefaultMomentSetForStencil(stencil)
rrDict = OrderedDict(zip(moments, relaxationRates))
if useContinuousMaxwellianEquilibrium:
......@@ -208,6 +224,8 @@ def createThreeRelaxationRateMRT(stencil, relaxationRates, useContinuousMaxwelli
"""
def product(iterable):
return reduce(operator.mul, iterable, 1)
if isinstance(stencil, str):
stencil = getStencil(stencil)
dim = len(stencil[0])
theMoment = MOMENT_SYMBOLS[:dim]
......@@ -337,6 +355,8 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, useContinuousMaxwell
"""
if relaxationRateGetter is None:
relaxationRateGetter = defaultRelaxationRateNames()
if isinstance(stencil, str):
stencil = getStencil(stencil)
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
......@@ -411,6 +431,19 @@ def createOrthogonalMRT(stencil, relaxationRateGetter=None, useContinuousMaxwell
return createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, **kwargs)
def modifyRelaxationTable(method, modificationFunction):
"""
Creates a new method, based on an existing method with modified relaxation table
:param method: old method
:param modificationFunction: function receiving (moment, equilibriumValue, relaxationRate) tuple, i.e. on row
of the relaxation table, returning a modified tuple
"""
relaxationTable = (modificationFunction(m, eq, rr)
for m, eq, rr in zip(method.moments, method.momentEquilibriumValues, method.relaxationRates))
compressible = method.conservedQuantityComputation.compressible
cumulant = isinstance(method, CumulantBasedLbMethod)
return createGenericMRT(method.stencil, relaxationTable, compressible, method.forceModel, cumulant)
# ----------------------------------------- Comparison view for notebooks ----------------------------------------------
......
......@@ -64,6 +64,10 @@ class CumulantBasedLbMethod(AbstractLbMethod):
self._weights = self._computeWeights()
return self._weights
def overrideWeights(self, weights):
assert len(weights) == len(self.stencil)
self._weights = weights
def _repr_html_(self):
table = """
<table style="border:none; width: 100%">
......
......@@ -70,6 +70,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._weights = self._computeWeights()
return self._weights
def overrideWeights(self, weights):
assert len(weights) == len(self.stencil)
self._weights = weights
def getEquilibrium(self, conservedQuantityEquations=None, includeForceTerms=False):
D = sp.eye(len(self.relaxationRates))
return self._getCollisionRuleWithRelaxationMatrix(D, conservedQuantityEquations=conservedQuantityEquations,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment