lbmpy: bugfixes in macroscopic value computation

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)
 

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):

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)

macroscopicValueKernels.py

-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:

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':