From 403c658e154711b5c39d294fa7e2509313a38eff Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Wed, 9 Dec 2020 16:18:55 +0100
Subject: [PATCH] fix shear rate again

---
 src/cpu/VirtualFluidsCore/LBM/D3Q27System.h   |  3 +-
 .../LBM/ThixotropyModelLBMKernel.cpp          |  2 +-
 .../LBM/ThixotropyModelLBMKernel2.cpp         | 32 ++++++++++++-------
 3 files changed, 22 insertions(+), 15 deletions(-)

diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index 1477a226e..a15a7a82b 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -31,7 +31,6 @@
 //! \author Konstantin Kutscher, Sebastian Geller, Soeren Freudiger
 //=======================================================================================
 
-
 #ifndef D3Q27SYSTEM_H
 #define D3Q27SYSTEM_H
 
@@ -1174,7 +1173,7 @@ usage: ...
 
          
          //TODO: may be factor 2
-         return sqrt(c2 * dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+         return sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
       }
    }
 #endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
index df0d9c457..d0d0542b6 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
@@ -480,7 +480,7 @@ void ThixotropyModelLBMKernel::calculate(int step)
 						LBMReal Dyz = -three * collFactorF * mfabb;
 						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 						//non Newtonian fluid collision factor
-						LBMReal shearRate = sqrt(c2 * dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+						LBMReal shearRate = sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
 						collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
 						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
index cafded555..da74d7a67 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
@@ -480,22 +480,30 @@ void ThixotropyModelLBMKernel2::calculate(int step)
 						LBMReal Dyz = -three * collFactorF * mfabb;
 						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 						//non Newtonian fluid collision factor
-						LBMReal shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
-						//collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
+						LBMReal shearRate = sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+						collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
 						////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 						//relax
+						//mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+						//
+						//LBMReal collFactorFyy = getThyxotropyCollFactor(collFactorF, std::sqrt(dxux*dxux + dyuy*dyuy) / (rho + one), rho);
+						//mxxMyy += collFactorFyy * (-mxxMyy) - 3. * (1. - c1o2 * collFactorFyy) * (vx2 * dxux - vy2 * dyuy);
+						//
+						//LBMReal collFactorFzz = getThyxotropyCollFactor(collFactorF, std::sqrt(dxux*dxux + dzuz*dzuz) / (rho + one), rho);
+						//mxxMzz += collFactorFzz * (-mxxMzz) - 3. * (1. - c1o2 * collFactorFzz) * (vx2 * dxux - vz2 * dzuz);
+
+						//mfabb += getThyxotropyCollFactor(collFactorF, std::abs(Dyz) / (rho + one), rho) * (-mfabb);
+						//mfbab += getThyxotropyCollFactor(collFactorF, std::abs(Dxz) / (rho + one), rho) * (-mfbab);
+						//mfbba += getThyxotropyCollFactor(collFactorF, std::abs(Dxy) / (rho + one), rho) * (-mfbba);
+
 						mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
-						
-						LBMReal collFactorFyy = getThyxotropyCollFactor(collFactorF, std::sqrt(dxux*dxux + dyuy*dyuy) / (rho + one), rho);
-						mxxMyy += collFactorFyy * (-mxxMyy) - 3. * (1. - c1o2 * collFactorFyy) * (vx2 * dxux - vy2 * dyuy);
-						
-						LBMReal collFactorFzz = getThyxotropyCollFactor(collFactorF, std::sqrt(dxux*dxux + dzuz*dzuz) / (rho + one), rho);
-						mxxMzz += collFactorFzz * (-mxxMzz) - 3. * (1. - c1o2 * collFactorFzz) * (vx2 * dxux - vz2 * dzuz);
-
-						mfabb += getThyxotropyCollFactor(collFactorF, std::abs(Dyz) / (rho + one), rho) * (-mfabb);
-						mfbab += getThyxotropyCollFactor(collFactorF, std::abs(Dxz) / (rho + one), rho) * (-mfbab);
-						mfbba += getThyxotropyCollFactor(collFactorF, std::abs(Dxy) / (rho + one), rho) * (-mfbba);
+						mxxMyy += collFactorF * (-mxxMyy) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vy2 * dyuy);
+						mxxMzz += collFactorF * (-mxxMzz) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vz2 * dzuz);
+
+						mfabb += collFactorF * (-mfabb);
+						mfbab += collFactorF * (-mfbab);
+						mfbba += collFactorF * (-mfbba);
 
 						// linear combinations back
 						mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
-- 
GitLab