diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index a5ead12f4a8c9bd11803166df1be2df8a841792f..0cb31f32c688201c8bd5c47ae181e3cd472d62f0 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -28,10 +28,10 @@ def ubb(pdfField, direction, lbMethod, velocity, adaptVelocityToForce=False):
c_s_sq = sp.Rational(1, 3)
velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
- # in conserved value computation
+ # Better alternative: in conserved value computation
# rename what is currently called density to "virtualDensity"
# provide a new quantity density, which is constant in case of incompressible models
- if lbMethod.conservedQuantityComputation._compressible: # TODO
+ if not lbMethod.conservedQuantityComputation.zeroCenteredPdfs:
cqc = lbMethod.conservedQuantityComputation
densitySymbol = sp.Symbol("rho")
pdfFieldAccesses = [pdfField(i) for i in range(len(lbMethod.stencil))]
diff --git a/creationfunctions.py b/creationfunctions.py
index fe5ca7a0415cec4a03d139c52399d051fd9b09e8..22b98e9c3521175f83d950a5adb546d3449c5294 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -4,7 +4,7 @@ Factory functions for standard LBM methods
import sympy as sp
from copy import copy
-from lbmpy.methods.creationfunctions import createKBCTypeTRT
+from lbmpy.methods.creationfunctions import createKBCTypeTRT, createRawMRT
from lbmpy.methods.entropic import addIterativeEntropyCondition, addEntropyCondition
from lbmpy.stencils import getStencil
from lbmpy.methods import createSRT, createTRT, createOrthogonalMRT
@@ -20,8 +20,10 @@ def _getParams(params, optParams):
'relaxationRates': sp.symbols("omega_:10"),
'compressible': False,
'equilibriumAccuracyOrder': 2,
+
'entropic': False,
'entropicNewtonIterations': None,
+ 'omegaOutputField': None,
'useContinuousMaxwellianEquilibrium': False,
'cumulant': False,
@@ -78,6 +80,7 @@ def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
return ValueError("'target' has to be either 'cpu' or 'gpu'")
res.method = ast.method
+ res.updateRule = ast.updateRule
res.ast = ast
return res
@@ -104,6 +107,7 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
return ValueError("'target' has to be either 'cpu' or 'gpu'")
res.method = updateRule.method
+ res.updateRule = updateRule
return res
@@ -128,9 +132,10 @@ def createLatticeBoltzmannUpdateRule(lbMethod=None, optimizationParams={}, **kwa
iterations = 3
else:
iterations = params['entropicNewtonIterations']
- collisionRule = addIterativeEntropyCondition(collisionRule, newtonIterations=iterations)
+ collisionRule = addIterativeEntropyCondition(collisionRule, newtonIterations=iterations,
+ omegaOutputField=params['omegaOutputField'])
else:
- collisionRule = addEntropyCondition(collisionRule)
+ collisionRule = addEntropyCondition(collisionRule, omegaOutputField=params['omegaOutputField'])
if 'fieldSize' in optParams and optParams['fieldSize']:
npField = createPdfArray(optParams['fieldSize'], len(stencil), layout=optParams['fieldLayout'])
@@ -193,6 +198,8 @@ def createLatticeBoltzmannMethod(**params):
nextRelaxationRate[0] += 1
return res
method = createOrthogonalMRT(stencil, relaxationRateGetter, **commonParams)
+ elif methodName.lower() == 'mrt_raw':
+ method = createRawMRT(stencil, relaxationRates, **commonParams)
elif methodName.lower().startswith('trt-kbc-n'):
if params['stencil'] == 'D2Q9':
dim = 2
diff --git a/forcemodels.py b/forcemodels.py
index 32d995e531ae6adb5c5dd69923298616ff4af8e7..8e8ebb9055d4781fa2b17a91292536476bc98156 100644
--- a/forcemodels.py
+++ b/forcemodels.py
@@ -4,6 +4,7 @@
"""
import sympy as sp
+from lbmpy.methods.relaxationrates import getShearRelaxationRate
class Simple(object):
@@ -71,7 +72,8 @@ class Guo(object):
def __call__(self, lbMethod):
luo = Luo(self._force)
- shearRelaxationRate = lbMethod.getShearRelaxationRate()
+
+ shearRelaxationRate = getShearRelaxationRate(lbMethod)
correctionFactor = (1 - sp.Rational(1, 2) * shearRelaxationRate)
return [correctionFactor * t for t in luo(lbMethod)]
diff --git a/macroscopic_value_kernels.py b/macroscopic_value_kernels.py
index 92ec1b73e91b6c83b7aa4fe545691d273cb1aee1..86ac145c5fe12f843812263f2a02b8e6614a3a2f 100644
--- a/macroscopic_value_kernels.py
+++ b/macroscopic_value_kernels.py
@@ -141,9 +141,8 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field
return setter
-def createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velocityRelaxationRate=0.8):
+def createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityRelaxationRate=0.8):
- lbMethod = collisionRule.method
velocityField = Field.createFromNumpyArray('velInput', velocityArray, indexDimensions=1)
cqc = lbMethod.conservedQuantityComputation
@@ -163,33 +162,20 @@ def createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velo
# 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 newCollisionRule
-def compileAdvancedVelocitySetter(collisionRule, velocityArray, velocityRelaxationRate=0.8, pdfArr=None,
+def compileAdvancedVelocitySetter(method, 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 collisionRule:
+ :param method:
:param velocityArray: array with velocity field
:param velocityRelaxationRate: relaxation rate for the velocity moments - determines convergence behaviour
of the initialization scheme
@@ -200,7 +186,7 @@ def compileAdvancedVelocitySetter(collisionRule, velocityArray, velocityRelaxati
"""
from lbmpy.updatekernels import createStreamPullKernel
from lbmpy.creationfunctions import createLatticeBoltzmannAst, createLatticeBoltzmannFunction
- newCollisionRule = createAdvancedVelocitySetterCollisionRule(collisionRule, velocityArray, velocityRelaxationRate)
+ newCollisionRule = createAdvancedVelocitySetterCollisionRule(method, velocityArray, velocityRelaxationRate)
updateRule = createStreamPullKernel(newCollisionRule, pdfArr, genericLayout=fieldLayout)
ast = createLatticeBoltzmannAst(updateRule, optimizationParams)
return createLatticeBoltzmannFunction(ast, optimizationParams)
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index de59c5dcb0e522b3f2ac6789b6366fa816348dba..152962de06e23f1cf09882212acd13e375bf385a 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -103,6 +103,10 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
else:
return None
+ @property
+ def zeroCenteredPdfs(self):
+ return not self._compressible
+
@property
def zerothOrderMomentSymbol(self):
return self._symbolOrder0
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 69b0b696d32ae9e3367eb16e26ef079a14295ff9..ea2183d3b4f9c2624d767563dd58ada8a596a0b0 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -1,19 +1,20 @@
import sympy as sp
-from collections import OrderedDict, defaultdict
+from collections import OrderedDict
from functools import reduce
import operator
import itertools
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
from lbmpy.methods.momentbased import MomentBasedLbMethod
from lbmpy.stencils import stencilsHaveSameEntries, getStencil
-from lbmpy.moments import isEven, gramSchmidt, getDefaultMomentSetForStencil, MOMENT_SYMBOLS, getOrder, isShearMoment, \
- exponentsToPolynomialRepresentations, momentsOfOrder, momentsUpToComponentOrder
+from lbmpy.moments import isEven, gramSchmidt, getDefaultMomentSetForStencil, MOMENT_SYMBOLS, \
+ exponentsToPolynomialRepresentations, momentsOfOrder, momentsUpToComponentOrder, sortMomentsIntoGroupsOfSameOrder
from pystencils.sympyextensions import commonDenominator
from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
from lbmpy.methods.abstractlbmethod import RelaxationInfo
from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium, \
getMomentsOfContinuousMaxwellianEquilibrium, getCumulantsOfDiscreteMaxwellianEquilibrium, \
getCumulantsOfContinuousMaxwellianEquilibrium
+from lbmpy.methods.relaxationrates import relaxationRateFromMagicNumber, defaultRelaxationRateNames
def createWithDiscreteMaxwellianEqMoments(stencil, momentToRelaxationRateDict, compressible=False, forceModel=None,
@@ -146,6 +147,18 @@ def createTRTWithMagicNumber(stencil, relaxationRate, magicNumber=sp.Rational(3,
return createTRT(stencil, relaxationRateEvenMoments=relaxationRate, relaxationRateOddMoments=rrOdd, **kwargs)
+def createRawMRT(stencil, relaxationRates, useContinuousMaxwellianEquilibrium=False, **kwargs):
+ """
+ Creates a MRT method using non-orthogonalized moments
+ """
+ moments = getDefaultMomentSetForStencil(stencil)
+ rrDict = OrderedDict(zip(moments, relaxationRates))
+ if useContinuousMaxwellianEquilibrium:
+ return createWithContinuousMaxwellianEqMoments(stencil, rrDict, **kwargs)
+ else:
+ return createWithDiscreteMaxwellianEqMoments(stencil, rrDict, **kwargs)
+
+
def createKBCTypeTRT(dim, shearRelaxationRate, higherOrderRelaxationRate, methodName='KBC-N4',
useContinuousMaxwellianEquilibrium=False, **kwargs):
"""
@@ -361,43 +374,6 @@ def compareMomentBasedLbMethods(reference, other, showDeviationsOnly=False):
# ------------------------------------ Helper Functions ----------------------------------------------------------------
-def sortMomentsIntoGroupsOfSameOrder(moments):
- """Returns a dictionary mapping the order (int) to a list of moments with that order."""
- result = defaultdict(list)
- for i, moment in enumerate(moments):
- order = getOrder(moment)
- result[order].append(moment)
- return result
-
-
-def defaultRelaxationRateNames():
- nextIndex = [0]
-
- def result(momentList):
- shearMomentInside = False
- allConservedMoments = True
- for m in momentList:
- if isShearMoment(m):
- shearMomentInside = True
- if not (getOrder(m) == 0 or getOrder(m) == 1):
- allConservedMoments = False
-
- if shearMomentInside:
- return sp.Symbol("omega")
- elif allConservedMoments:
- return 0
- else:
- nextIndex[0] += 1
- return sp.Symbol("omega_%d" % (nextIndex[0],))
-
- return result
-
-
-def relaxationRateFromMagicNumber(hydrodynamicRelaxationRate, magicNumber):
- omega = hydrodynamicRelaxationRate
- return (4 - 2 * omega) / (4 * magicNumber * omega + 2 - omega)
-
-
def compressibleToIncompressibleMomentValue(term, rho, u):
term = term.expand()
if term.func != sp.Add:
diff --git a/methods/entropic.py b/methods/entropic.py
index 751c6ef592e2ca61629f6b27db40534525d7d986..07940a58d4ee5cec8474f8910267445d98afb315 100644
--- a/methods/entropic.py
+++ b/methods/entropic.py
@@ -1,8 +1,9 @@
import sympy as sp
from pystencils.transformations import fastSubs
+from lbmpy.methods.relaxationrates import getShearRelaxationRate
-def addEntropyCondition(collisionRule):
+def addEntropyCondition(collisionRule, omegaOutputField=None):
"""
Transforms an update rule with two relaxation rate into a single relaxation rate rule, where the second
rate is locally chosen to maximize an entropy condition. This function works for update rules which are
@@ -13,8 +14,12 @@ def addEntropyCondition(collisionRule):
have to be done.
:param collisionRule: collision rule with two relaxation times
+ :param omegaOutputField: pystencils field where computed omegas are stored
:return: new collision rule which only one relaxation rate
"""
+ if collisionRule.method.conservedQuantityComputation.zeroCenteredPdfs:
+ raise NotImplementedError("Entropic Methods only implemented for models where pdfs are centered around 1")
+
omega_s, omega_h = _getRelaxationRates(collisionRule)
decomp = RelaxationRatePolynomialDecomposition(collisionRule, [omega_h], [omega_s])
@@ -25,8 +30,10 @@ def addEntropyCondition(collisionRule):
dh.append(entry[0])
ds = []
for entry in decomp.relaxationRateFactors(omega_s):
- assert len(entry) == 1, "The non-iterative entropic procedure works only for moment based methods, which have" \
+ assert len(entry) <= 1, "The non-iterative entropic procedure works only for moment based methods, which have" \
"an update rule linear in the relaxation rate."
+ if len(entry) == 0:
+ entry.append(0)
ds.append(entry[0])
stencil = collisionRule.method.stencil
@@ -53,10 +60,15 @@ def addEntropyCondition(collisionRule):
newCollisionRule = collisionRule.copy(newUpdateEquations, collisionRule.subexpressions + subexprs)
newCollisionRule.simplificationHints['entropic'] = True
newCollisionRule.simplificationHints['entropicNewtonIterations'] = None
+
+ if omegaOutputField:
+ from lbmpy.updatekernels import writeQuantitiesToField
+ newCollisionRule = writeQuantitiesToField(newCollisionRule, omega_h, omegaOutputField)
+
return newCollisionRule
-def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue=1):
+def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue=1, omegaOutputField=None):
"""
More generic, but slower version of :func:`addEntropyCondition`
@@ -65,8 +77,12 @@ def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue
:param collisionRule: collision rule with two relaxation times
:param newtonIterations: (integer) number of newton iterations
:param initialValue: initial value of the relaxation rate
+ :param omegaOutputField: pystencils field where computed omegas are stored
:return: new collision rule which only one relaxation rate
"""
+ if collisionRule.method.conservedQuantityComputation.zeroCenteredPdfs:
+ raise NotImplementedError("Entropic Methods only implemented for models where pdfs are centered around 1")
+
omega_s, omega_h = _getRelaxationRates(collisionRule)
decomp = RelaxationRatePolynomialDecomposition(collisionRule, [omega_h], [omega_s])
@@ -115,6 +131,11 @@ def addIterativeEntropyCondition(collisionRule, newtonIterations=3, initialValue
newCollisionRule = collisionRule.copy(newUpdateEquations, newSubexpressions)
newCollisionRule.simplificationHints['entropic'] = True
newCollisionRule.simplificationHints['entropicNewtonIterations'] = newtonIterations
+
+ if omegaOutputField:
+ from lbmpy.updatekernels import writeQuantitiesToField
+ newCollisionRule = writeQuantitiesToField(newCollisionRule, omega_h, omegaOutputField)
+
return newCollisionRule
@@ -219,7 +240,7 @@ def _getRelaxationRates(collisionRule):
"entropy condition")
method = collisionRule.method
- omega_s = method.getShearRelaxationRate()
+ omega_s = getShearRelaxationRate(method)
assert omega_s in relaxationRates
relaxationRatesWithoutOmegaS = relaxationRates - {omega_s}
diff --git a/methods/relaxationrates.py b/methods/relaxationrates.py
new file mode 100644
index 0000000000000000000000000000000000000000..a1033da84a99773a47d1b6f6c568491af0ffb9b7
--- /dev/null
+++ b/methods/relaxationrates.py
@@ -0,0 +1,53 @@
+import sympy as sp
+from lbmpy.moments import isShearMoment, getOrder
+
+
+def relaxationRateFromMagicNumber(hydrodynamicRelaxationRate, magicNumber=sp.Rational(3, 16)):
+ """
+ Computes second TRT relaxation rate from magic number
+ """
+ omega = hydrodynamicRelaxationRate
+ return (4 - 2 * omega) / (4 * magicNumber * omega + 2 - omega)
+
+
+def getShearRelaxationRate(method):
+ """
+ Assumes that all shear moments are relaxed with same rate - returns this rate
+ Shear moments in 3D are: x*y, x*z and y*z - in 2D its only x*y
+ The shear relaxation rate determines the viscosity in hydrodynamic LBM schemes
+ """
+ relaxationRates = set()
+ for moment, relaxInfo in method.relaxationInfoDict.items():
+ if isShearMoment(moment):
+ relaxationRates.add(relaxInfo.relaxationRate)
+ if len(relaxationRates) == 1:
+ return relaxationRates.pop()
+ else:
+ if len(relaxationRates) > 1:
+ raise ValueError("Shear moments are relaxed with different relaxation times: %s" % (relaxationRates,))
+ else:
+ raise NotImplementedError("Shear moments seem to be not relaxed separately - "
+ "Can not determine their relaxation rate automatically")
+
+
+def defaultRelaxationRateNames():
+ nextIndex = [0]
+
+ def result(momentList):
+ shearMomentInside = False
+ allConservedMoments = True
+ for m in momentList:
+ if isShearMoment(m):
+ shearMomentInside = True
+ if not (getOrder(m) == 0 or getOrder(m) == 1):
+ allConservedMoments = False
+
+ if shearMomentInside:
+ return sp.Symbol("omega")
+ elif allConservedMoments:
+ return 0
+ else:
+ nextIndex[0] += 1
+ return sp.Symbol("omega_%d" % (nextIndex[0],))
+
+ return result
\ No newline at end of file
diff --git a/moments.py b/moments.py
index b2209ec20854424a4bc08ce1aa42da385a0f5958..a7874d26b1c2b8c60af03a461a3b3baf946c1109 100644
--- a/moments.py
+++ b/moments.py
@@ -38,7 +38,7 @@ import sympy as sp
import math
import itertools
from copy import copy
-from collections import Counter
+from collections import Counter, defaultdict
from lbmpy.continuous_distribution_measures import continuousMoment
from lbmpy.cache import memorycache
@@ -210,6 +210,14 @@ def polynomialToExponentRepresentation(polynomial, dim=3):
return coeffExpTupleRepresentation
+def sortMomentsIntoGroupsOfSameOrder(moments):
+ """Returns a dictionary mapping the order (int) to a list of moments with that order."""
+ result = defaultdict(list)
+ for i, moment in enumerate(moments):
+ order = getOrder(moment)
+ result[order].append(moment)
+ return result
+
# -------------------- Common Function working with exponent tuples and polynomial moments -----------------------------
diff --git a/serialscenario.py b/serialscenario.py
index 73d342e29d799c8dd4d944727593e44e40c5ce13..c30d443629820e0099a0613217138bb69cd6e6e4 100644
--- a/serialscenario.py
+++ b/serialscenario.py
@@ -9,7 +9,7 @@ from lbmpy.stencils import getStencil
def createScenario(domainSize, boundarySetupFunction, methodParameters, optimizationParams, lbmKernel=None,
- initialVelocity=None, preUpdateFunctions=[]):
+ initialVelocity=None, preUpdateFunctions=[], kernelParams={}):
if 'target' not in optimizationParams:
optimizationParams['target'] = 'cpu'
@@ -65,7 +65,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
f(pdfArrays[0])
if boundaryHandling is not None:
boundaryHandling(pdfs=pdfArrays[0])
- lbmKernel(src=pdfArrays[0], dst=pdfArrays[1])
+ lbmKernel(src=pdfArrays[0], dst=pdfArrays[1], **kernelParams)
pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0]
getMacroscopic(pdfs=pdfArrays[0], density=densityArr[0], velocity=velocityArr[0])
@@ -81,7 +81,7 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
f(pdfGpuArrays[0])
if boundaryHandling is not None:
boundaryHandling(pdfs=pdfGpuArrays[0])
- lbmKernel(src=pdfGpuArrays[0], dst=pdfGpuArrays[1])
+ lbmKernel(src=pdfGpuArrays[0], dst=pdfGpuArrays[1], **kernelParams)
pdfGpuArrays[0], pdfGpuArrays[1] = pdfGpuArrays[1], pdfGpuArrays[0]
@@ -98,7 +98,8 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza
return gpuTimeLoop if optimizationParams['target'] == 'gpu' else cpuTimeLoop
-def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None, **kwargs):
+def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, lbmKernel=None,
+ kernelParams={}, **kwargs):
def boundarySetupFunction(boundaryHandling, method):
myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:method.dim])
myUbb.name = 'ubb'
@@ -106,11 +107,12 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={},
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, optimizationParams, lbmKernel=lbmKernel)
+ return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel,
+ kernelParams=kernelParams)
def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None,
- lbmKernel=None, optimizationParams={}, **kwargs):
+ lbmKernel=None, optimizationParams={}, kernelParams={}, **kwargs):
assert dim in (2, 3)
if radius is not None:
@@ -149,14 +151,13 @@ def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None
boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim))
assert domainSize is not None
- if 'forceModel' not in kwargs:
- kwargs['forceModel'] = 'guo'
-
- return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel)
+ return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel,
+ kernelParams=kernelParams)
def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=None, lbmKernel=None,
- optimizationParams={}, initialVelocity=None, boundarySetupFunctions=[], **kwargs):
+ optimizationParams={}, initialVelocity=None, boundarySetupFunctions=[],
+ kernelParams={}, **kwargs):
assert dim in (2, 3)
kwargs['force'] = tuple([force, 0, 0][:dim])
@@ -198,5 +199,6 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No
kwargs['forceModel'] = 'guo'
return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel,
- initialVelocity=initialVelocity, preUpdateFunctions=[periodicity])
+ initialVelocity=initialVelocity, preUpdateFunctions=[periodicity],
+ kernelParams=kernelParams)