From 75b1d916b7b970083a0f2a70edd58acb7ffefa68 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Wed, 1 Nov 2017 18:47:59 +0100
Subject: [PATCH] Simple functions to compute force on boundary
- added functions to compute force on boundary for generic boundary and
NoSlip
- simple calculations for serial scenarios
- to distinguish between different obstacles, NoSlip boundaries can now
be named
---
boundaries/boundaryconditions.py | 20 ++++++++++--
boundaries/forceevaluation.py | 52 ++++++++++++++++++++++++++++++++
boundaries/handlinginterface.py | 4 +++
scenarios.py | 17 +++++++++++
4 files changed, 90 insertions(+), 3 deletions(-)
create mode 100644 boundaries/forceevaluation.py
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 19d2d690..e7146e3f 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -51,6 +51,11 @@ class Boundary(object):
class NoSlip(Boundary):
+
+ def __init__(self, name=None):
+ """Set an optional name here, to mark boundaries, for example for force evaluations"""
+ self._name = name
+
"""No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
@@ -58,11 +63,20 @@ class NoSlip(Boundary):
return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(directionSymbol))]
def __hash__(self):
- # All boundaries of these class behave equal -> should also be equal
- return hash("NoSlip")
+ # All boundaries of these class behave equal -> should also be equal (as long as name is equal)
+ return hash(self.name)
def __eq__(self, other):
- return type(other) == NoSlip
+ if not isinstance(other, NoSlip):
+ return False
+ return self.name == other.name
+
+ @property
+ def name(self):
+ if self._name:
+ return self._name
+ else:
+ return type(self).__name__
class NoSlipFullWay(Boundary):
diff --git a/boundaries/forceevaluation.py b/boundaries/forceevaluation.py
new file mode 100644
index 00000000..3e741d73
--- /dev/null
+++ b/boundaries/forceevaluation.py
@@ -0,0 +1,52 @@
+import numpy as np
+from lbmpy.stencils import inverseDirection
+
+
+def calculateForceOnNoSlipBoundary(boundaryObject, boundaryHandling, pdfArray):
+ indArr = boundaryHandling.getBoundaryIndexArray(boundaryObject)
+ method = boundaryHandling.lbMethod
+
+ stencil = np.array(method.stencil)
+
+ if method.dim == 2:
+ x, y = indArr['x'], indArr['y']
+ values = 2 * pdfArray[x, y, indArr['dir']]
+ else:
+ assert method.dim == 3
+ x, y, z = indArr['x'], indArr['y'], indArr['z']
+ values = 2 * pdfArray[x, y, z, indArr['dir']]
+
+ forces = stencil[indArr['dir']] * values[:, np.newaxis]
+ return forces.sum(axis=0)
+
+
+def calculateForceOnBoundary(boundaryObject, boundaryHandling, pdfArray):
+ indArr = boundaryHandling.getBoundaryIndexArray(boundaryObject)
+ method = boundaryHandling.lbMethod
+
+ stencil = np.array(method.stencil)
+ invDirection = np.array([method.stencil.index(inverseDirection(d))
+ for d in method.stencil])
+
+ if method.dim == 2:
+ x, y = indArr['x'], indArr['y']
+ offsets = stencil[indArr['dir']]
+
+ fluidValues = pdfArray[x, y, indArr['dir']]
+ boundaryValues = pdfArray[x + offsets[:, 0],
+ y + offsets[:, 1],
+ invDirection[indArr['dir']]]
+ else:
+ assert method.dim == 3
+ x, y, z = indArr['x'], indArr['y'], indArr['z']
+ offsets = stencil[indArr['dir']]
+
+ fluidValues = pdfArray[x, y, z, indArr['dir']]
+ boundaryValues = pdfArray[x + offsets[:, 0],
+ y + offsets[:, 1],
+ z + offsets[:, 2],
+ invDirection[indArr['dir']]]
+
+ values = fluidValues + boundaryValues
+ forces = stencil[indArr['dir']] * values[:, np.newaxis]
+ return forces.sum(axis=0)
diff --git a/boundaries/handlinginterface.py b/boundaries/handlinginterface.py
index de3f28fb..ba768e4c 100644
--- a/boundaries/handlinginterface.py
+++ b/boundaries/handlinginterface.py
@@ -113,6 +113,10 @@ class GenericBoundaryHandling(object):
if target not in ('cpu', 'gpu'):
raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
+ @property
+ def lbMethod(self):
+ return self._lbMethod
+
def getBoundaryNameToFlagDict(self):
result = {self._flagFieldInterface.getName(o): self._flagFieldInterface.getFlag(o) for o in self._boundaryInfos}
result['fluid'] = self._flagFieldInterface.getFlag('fluid')
diff --git a/scenarios.py b/scenarios.py
index d1ab830f..0687bb27 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -26,6 +26,7 @@ that defines the viscosity of the fluid (valid values being between 0 and 2).
import numpy as np
import sympy as sp
+from lbmpy.boundaries.forceevaluation import calculateForceOnBoundary, calculateForceOnNoSlipBoundary
from lbmpy.geometry import setupChannelWalls, addParabolicVelocityInflow
from pystencils.field import getLayoutOfArray, createNumpyArrayWithLayout
from pystencils.slicing import sliceFromDirection, addGhostLayers, removeGhostLayers, normalizeSlice, makeSlice
@@ -346,6 +347,14 @@ class Scenario(object):
"""Equation collection defining the LBM update rule (already in simplified form)"""
return self._lbmKernel.updateRule
+ def calculateForceOnBoundary(self, boundaryObject):
+ """Computes force on boundary using simple momentum exchange method"""
+ if isinstance(boundaryObject, NoSlip):
+ return calculateForceOnNoSlipBoundary(boundaryObject, self.boundaryHandling, self._pdfArrays[0])
+ else:
+ self.runBoundaryHandlingOnly()
+ return calculateForceOnBoundary(boundaryObject, self.boundaryHandling, self._pdfArrays[0])
+
def animateVelocity(self, steps=10, **kwargs):
import lbmpy.plot2d as plt
@@ -415,3 +424,11 @@ class Scenario(object):
self._getMacroscopic(pdfs=self._pdfArrays[0], density=self._density, velocity=self._velocity,
**self.kernelParams)
+ def runBoundaryHandlingOnly(self):
+ isGpuSimulation = len(self._pdfGpuArrays) > 0
+ if isGpuSimulation:
+ self._pdfGpuArrays[0].set(self._pdfArrays[0])
+ self._boundaryHandling(pdfs=self._pdfGpuArrays[0], **self.kernelParams)
+ self._pdfGpuArrays[0].get(self._pdfArrays[0])
+ else:
+ self._boundaryHandling(pdfs=self._pdfArrays[0], **self.kernelParams)
--
GitLab