From a52302dbe0d9d1b175c73ac0d55825b0b30531cf Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Mon, 13 Feb 2017 15:52:49 +0100
Subject: [PATCH] lbmpy: bugfixes in macroscopic value computation
---
boundaries/boundaryconditions.py | 8 ++++++--
boundaries/boundaryhandling.py | 1 +
forcemodels.py | 24 ++++++++++++++----------
macroscopicValueKernels.py | 17 +++++++++++++----
methods/conservedquantitycomputation.py | 2 +-
5 files changed, 35 insertions(+), 17 deletions(-)
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 459c8799..26785bb8 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -12,7 +12,7 @@ def noSlip(pdfField, direction, lbMethod):
return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(direction))]
-def ubb(pdfField, direction, lbMethod, velocity):
+def ubb(pdfField, direction, lbMethod, velocity, adaptVelocityToForce=False):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
assert len(velocity) == lbMethod.dim, \
@@ -20,7 +20,11 @@ def ubb(pdfField, direction, lbMethod, velocity):
neighbor = offsetFromDir(direction, lbMethod.dim)
inverseDir = invDir(direction)
- # TODO adapt velocity to force
+ if adaptVelocityToForce:
+ cqc = lbMethod.conservedQuantityComputation
+ shiftedVelEqs = cqc.equilibriumInputEquationsFromInitValues(velocity=velocity)
+ velocity = [eq.rhs for eq in shiftedVelEqs.extract(cqc.firstOrderMomentSymbols).mainEquations]
+
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)
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 14e40146..213fc1f6 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -140,6 +140,7 @@ def generateBoundaryHandling(pdfField, indexArr, lbMethod, boundaryFunctor):
cellLoop = LoopOverCoordinate(cellLoopBody, coordinateToLoopOver=0, start=0, stop=indexArr.shape[0])
indexField = Field.createFromNumpyArray("indexField", indexArr, indexDimensions=1)
+ indexField.isIndexField = True
coordinateSymbols = [TypedSymbol(name, "int") for name in ['x', 'y', 'z']]
for d in range(dim):
diff --git a/forcemodels.py b/forcemodels.py
index 5b7fe308..ca0c4203 100644
--- a/forcemodels.py
+++ b/forcemodels.py
@@ -16,10 +16,14 @@ class Simple:
def __init__(self, force):
self._force = force
- def __call__(self, lbmMethod, **kwargs):
- assert len(self._force) == lbmMethod.dim
- return [3 * w_i * sum([d_i * f_i for d_i, f_i in zip(direction, self._force)])
- for direction, w_i in zip(lbmMethod.stencil, lbmMethod.weights)]
+ def __call__(self, lbMethod, **kwargs):
+ assert len(self._force) == lbMethod.dim
+
+ def scalarProduct(a, b):
+ return sum(a_i * b_i for a_i, b_i in zip(a, b))
+
+ return [3 * w_i * scalarProduct(self._force, direction)
+ for direction, w_i in zip(lbMethod.stencil, lbMethod.weights)]
class Luo:
@@ -35,12 +39,12 @@ class Luo:
def __init__(self, force):
self._force = force
- def __call__(self, lbmMethod):
- u = sp.Matrix(lbmMethod.firstOrderEquilibriumMomentSymbols)
+ def __call__(self, lbMethod):
+ u = sp.Matrix(lbMethod.firstOrderEquilibriumMomentSymbols)
force = sp.Matrix(self._force)
result = []
- for direction, w_i in zip(lbmMethod.stencil, lbmMethod.weights):
+ for direction, w_i in zip(lbMethod.stencil, lbMethod.weights):
direction = sp.Matrix(direction)
result.append(3 * w_i * force.dot(direction - u + 3 * direction * direction.dot(u)))
return result
@@ -62,11 +66,11 @@ class Guo:
def __init__(self, force):
self._force = force
- def __call__(self, lbmMethod):
+ def __call__(self, lbMethod):
luo = Luo(self._force)
- shearRelaxationRate = lbmMethod.getShearRelaxationRate()
+ shearRelaxationRate = lbMethod.getShearRelaxationRate()
correctionFactor = (1 - sp.Rational(1, 2) * shearRelaxationRate)
- return [correctionFactor * t for t in luo(lbmMethod)]
+ return [correctionFactor * t for t in luo(lbMethod)]
def macroscopicVelocityShift(self, density):
return defaultVelocityShift(density, self._force)
diff --git a/macroscopicValueKernels.py b/macroscopicValueKernels.py
index 084a6f6a..4061adfe 100644
--- a/macroscopicValueKernels.py
+++ b/macroscopicValueKernels.py
@@ -1,4 +1,3 @@
-import sympy as sp
from copy import deepcopy
from pystencils.field import Field, getLayoutFromNumpyArray
from lbmpy.simplificationfactory import createSimplificationStrategy
@@ -30,14 +29,24 @@ def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fiel
pdfField = Field.createGeneric('pdfs', lbMethod.dim, indexDimensions=1, layout=fieldLayout)
else:
pdfField = Field.createFromNumpyArray('pdfs', pdfArr, indexDimensions=1)
- fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape)-1])
outputMapping = {}
for outputQuantity in outputQuantities:
numberOfElements = cqc.conservedQuantities[outputQuantity]
assert numberOfElements >= 1
- outputField = Field.createGeneric(outputQuantity, lbMethod.dim, layout=fieldLayout,
- indexDimensions=0 if numberOfElements <= 1 else 1)
+
+ indDims = 0 if numberOfElements <= 1 else 1
+ if pdfArr is None:
+ fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape) - 1])
+ outputField = Field.createGeneric(outputQuantity, lbMethod.dim, layout=fieldLayout, indexDimensions=indDims)
+ else:
+ outputFieldShape = pdfArr.shape[:-1]
+ if indDims > 0:
+ outputFieldShape += (numberOfElements,)
+ fieldLayout = getLayoutFromNumpyArray(pdfArr)
+ else:
+ fieldLayout = getLayoutFromNumpyArray(pdfArr, indexDimensionIds=[len(pdfField.shape) - 1])
+ outputField = Field.createFixedSize(outputQuantity, outputFieldShape, indDims, pdfArr.dtype, fieldLayout)
outputMapping[outputQuantity] = [outputField(i) for i in range(numberOfElements)]
if len(outputMapping[outputQuantity]) == 1:
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 34060528..41f31460 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -90,7 +90,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
@property
def conservedQuantities(self):
return {'density': 1,
- 'velocity': 3}
+ 'velocity': len(self._stencil[0])}
def definedSymbols(self, order='all'):
if order == 'all':
--
GitLab