From 71219dbb204e7cfc3a04c09090c95f3248281ce0 Mon Sep 17 00:00:00 2001
From: Soeren Peters <peters.soeren@gmx.net>
Date: Thu, 11 Feb 2021 10:53:45 +0100
Subject: [PATCH] Bugfix: fReturn should be calculated with collFactorF.

---
 .../BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp  | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
index 9f6ca522a..4ea802eb4 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
@@ -26,7 +26,10 @@ void RheologicalVelocityBCAlgorithm::applyBC()
    calcMacrosFct(f, drho, vx1, vx2, vx3);
    calcFeqFct(feq, drho, vx1, vx2, vx3);
 
-   rho = 1.0+drho*compressibleFactor;
+    LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
+    LBMReal collFactorF = getThyxotropyCollFactor(collFactor, shearRate, drho);
+
+    rho = 1.0+drho*compressibleFactor;
 
    for (int fdir = D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++)
    {
@@ -35,7 +38,7 @@ void RheologicalVelocityBCAlgorithm::applyBC()
          const int invDir = D3Q27System::INVDIR[fdir];
          LBMReal q = bcPtr->getQ(invDir);// m+m q=0 stabiler
          LBMReal velocity = bcPtr->getBoundaryVelocity(invDir);
-         LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity*rho)/(1.0+q));
+         LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactorF)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity*rho)/(1.0+q));
          distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
       }
    }
-- 
GitLab