From 6f0c299eedaa455ec9ee8aef9050fa1dbc5c755b Mon Sep 17 00:00:00 2001
From: "AMATERASU\\geier" <geier@irmb.tu-bs.de>
Date: Wed, 11 Oct 2023 13:28:15 +0200
Subject: [PATCH] Bug fix regarding the interface condition at a solid wall. Is
 here implemented only for BB without delay

---
 .../MultiphaseScaleDistributionLBMKernel.cpp  | 24 +++++++++++++++----
 1 file changed, 19 insertions(+), 5 deletions(-)

diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
index e791f0e7e..5d7500c69 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
@@ -627,7 +627,9 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 							real sumRho = 0;
 							real sumWeight = 1.e-100;
 							for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
-                                if ((phi[fdir] <= phiLim)) {
+                                if ((phi[fdir] <= phiLim) &&
+                                    !bcArray->isSolid(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+                                                      x3 + D3Q27System::DX3[fdir])) {
 									sumRho += WEIGTH[fdir] * (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]);
 									sumWeight += WEIGTH[fdir];
 								}
@@ -1191,7 +1193,9 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 								// real sumVy = 0;
 								// real sumVz = 0;
 								for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
-                                    if ((phi[fdir] > phiLim)) {
+                                    if ((phi[fdir] > phiLim) &&
+                                        !bcArray->isSolid(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+                                                          x3 + D3Q27System::DX3[fdir])) {
 
 										sumRho += WEIGTH[fdir] * (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]);// * tempRho;
 										// sumVx += WEIGTH[fdir] * (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]);
@@ -1301,9 +1305,19 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
                                         real fG = distribution->getDistributionInvForDirection(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
 										real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
 										real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoL, (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-                                        real fBC = feqNew + (fG - feqOLD) * (c1o1 / collFactorL - c1o1) / (c1o1 / collFactorG - c1o1);
-										//not adding laplace pressure as it should be included in rhoL;
-                                        ff[D3Q27System::INVDIR[fdir]] = fBC;
+                                       //11.10.2023 correction for BB BC assuming no delay BB
+                                        if (bcArray->isSolid(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+                                                             x3 + D3Q27System::DX3[fdir])){
+										feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 , x2 , x3 ),vx,vy,vz);
+										feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoL,vx,vy,vz);
+										
+										}//!11.10.2023
+
+
+                                        real fBC = feqNew +
+                                                  (fG - feqOLD) * (c1o1 / collFactorL - c1o1) / (c1o1 / collFactorG - c1o1);
+ 
+										ff[D3Q27System::INVDIR[fdir]] = fBC;
 										//15.5.23
 										//ff[D3Q27System::INVDIR[fdir]] = (fGEQ - WEIGTH[fdir] * dvDir * (c1o1 / collFactorL)) - c6o1 * WEIGTH[fdir] * (vDir);
 										//real fBC = (distribution->getDistributionInvForDirection(x1, x2, x3, D3Q27System::INVDIR[fdir]) - c6o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorL - c1o1)) ;
-- 
GitLab