From 3efcd2835d4caf853f00c1fc114251902ccdbd95 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Thu, 2 Mar 2017 13:05:17 +0100
Subject: [PATCH] lbmpy: bugfix in pressure/density boundary for incompressible
methods
---
boundaries/boundaryconditions.py | 12 +++++++++---
1 file changed, 9 insertions(+), 3 deletions(-)
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 26785bb8..a5ead12f 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -56,7 +56,8 @@ def fixedDensity(pdfField, direction, lbMethod, density):
neighbor = offsetFromDir(direction, lbMethod.dim)
inverseDir = invDir(direction)
- velocity = lbMethod.conservedQuantityComputation.definedSymbols()['velocity']
+ cqc = lbMethod.conservedQuantityComputation
+ velocity = cqc.definedSymbols()['velocity']
symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
@@ -64,12 +65,17 @@ def fixedDensity(pdfField, direction, lbMethod, density):
simplification = createSimplificationStrategy(lbMethod)
symmetricEq = simplification(symmetricEq)
- densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
+ densitySymbol = cqc.definedSymbols()['density']
+
+ densityEq = cqc.equilibriumInputEquationsFromInitValues(density=density).insertSubexpressions().mainEquations[0]
+ assert densityEq.lhs == densitySymbol
+ transformedDensity = densityEq.rhs
conditions = [(eq_i.rhs, sp.Equality(direction, i))
for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
eq_component = sp.Piecewise(*conditions)
- subExprs = [sp.Eq(eq.lhs, density if eq.lhs == densitySymbol else eq.rhs) for eq in symmetricEq.subexpressions]
+ subExprs = [sp.Eq(eq.lhs, transformedDensity if eq.lhs == densitySymbol else eq.rhs)
+ for eq in symmetricEq.subexpressions]
return subExprs + [sp.Eq(pdfField[neighbor](inverseDir), 2 * eq_component - pdfField(direction))]
--
GitLab