From bb81981d9af7a81270e06285130cf9c9de3d5f63 Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Tue, 26 Oct 2021 15:13:38 +0200
Subject: [PATCH] fix viscosity correction and optimize muParser variable
definition in MultiphaseTwoPhaseFieldsPressureFilterLBMKernel
---
.../cpu/MultiphaseDropletTest/DropletTest.cfg | 4 +-
apps/cpu/MultiphaseDropletTest/droplet.cpp | 41 +++----
...eTwoPhaseFieldsPressureFilterLBMKernel.cpp | 107 ++++++++++++------
...aseTwoPhaseFieldsPressureFilterLBMKernel.h | 17 ++-
4 files changed, 109 insertions(+), 60 deletions(-)
diff --git a/apps/cpu/MultiphaseDropletTest/DropletTest.cfg b/apps/cpu/MultiphaseDropletTest/DropletTest.cfg
index 174066574..a93c46436 100644
--- a/apps/cpu/MultiphaseDropletTest/DropletTest.cfg
+++ b/apps/cpu/MultiphaseDropletTest/DropletTest.cfg
@@ -1,5 +1,5 @@
#pathname = d:/temp/MultiphaseDropletTest
-pathname = E:/Multiphase/DropletTestImpVel-Hesam2
+pathname = E:/Multiphase/DropletTest_new_corr5
numOfThreads = 4
availMem = 10e9
@@ -13,7 +13,7 @@ dx = 1
refineLevel = 0
#Simulation
-uLB = 0.01#0.005#0.005
+uLB = 0 #0.001#0.005#0.005
Re = 10
nuL = 1e-2 #1e-5# 1.0e-5 #!1e-2
nuG = 0.015811388300841892 #5e-2 #1e-4 # 1e-8 # 1.16e-4 #!1e-2
diff --git a/apps/cpu/MultiphaseDropletTest/droplet.cpp b/apps/cpu/MultiphaseDropletTest/droplet.cpp
index 269544600..451133022 100644
--- a/apps/cpu/MultiphaseDropletTest/droplet.cpp
+++ b/apps/cpu/MultiphaseDropletTest/droplet.cpp
@@ -109,11 +109,11 @@ void run(string configname)
LBMReal mu_h = rho_h * nu_h;
//gravity
- LBMReal g_y = Re*Re*mu_h*mu_h/(rho_h*(rho_h-rho_l)*D*D*D);
+ LBMReal g_y = Re* Re* mu_h* mu_h / (rho_h * (rho_h - rho_l) * D * D * D);
//Eotvos number
LBMReal Eo = 100;
//surface tension
- sigma = rho_h*g_y*D*D/Eo;
+ sigma = rho_h* g_y* D* D / Eo;
//g_y = 0;
@@ -125,7 +125,7 @@ void run(string configname)
//UBLOG(logINFO, "rho = " << rhoLB);
UBLOG(logINFO, "D = " << D);
UBLOG(logINFO, "nuL = " << nuL);
- UBLOG(logINFO, "nuL = " << nuG);
+ UBLOG(logINFO, "nuG = " << nuG);
UBLOG(logINFO, "Re = " << Re);
UBLOG(logINFO, "Eo = " << Eo);
UBLOG(logINFO, "g_y = " << g_y);
@@ -141,18 +141,18 @@ void run(string configname)
SPtr<LBMKernel> kernel;
//kernel = SPtr<LBMKernel>(new MultiphaseScratchCumulantLBMKernel());
- kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
- //kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
+ kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
- //mu::Parser fgr;
- //fgr.SetExpr("-(rho-rho_l)*g_y");
- //fgr.DefineConst("rho_l", rho_l);
- //fgr.DefineConst("g_y", g_y);
+ mu::Parser fgr;
+ fgr.SetExpr("-(rho-rho_l)*g_y");
+ fgr.DefineConst("rho_l", rho_l);
+ fgr.DefineConst("g_y", g_y);
- //kernel->setWithForcing(true);
- //kernel->setForcingX1(0.0);
- //kernel->setForcingX2(fgr);
- //kernel->setForcingX3(0.0);
+ kernel->setWithForcing(true);
+ kernel->setForcingX1(0.0);
+ kernel->setForcingX2(fgr);
+ kernel->setForcingX3(0.0);
kernel->setPhiL(phiL);
kernel->setPhiH(phiH);
@@ -184,7 +184,7 @@ void run(string configname)
grid->setPeriodicX1(true);
grid->setPeriodicX2(true);
grid->setPeriodicX3(true);
- grid->setGhostLayerWidth(1);
+ grid->setGhostLayerWidth(2);
SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::RECURSIVE));
@@ -298,6 +298,7 @@ void run(string configname)
mu::Parser fct2;
fct2.SetExpr("0.5*uLB-uLB*0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
+ //fct2.SetExpr("uLB");
fct2.DefineConst("uLB", uLB);
fct2.DefineConst("x1c", x1c);
fct2.DefineConst("x2c", x2c);
@@ -305,8 +306,8 @@ void run(string configname)
fct2.DefineConst("radius", radius);
fct2.DefineConst("interfaceThickness", interfaceThickness);
- MultiphaseInitDistributionsBlockVisitor initVisitor(densityRatio);
- //MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
+ //MultiphaseInitDistributionsBlockVisitor initVisitor(densityRatio);
+ MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
initVisitor.setPhi(fct1);
initVisitor.setVx1(fct2);
grid->accept(initVisitor);
@@ -344,12 +345,12 @@ void run(string configname)
grid->accept(bcVisitor);
- TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
- grid->accept(setConnsVisitor);
-
- //ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ //TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
//grid->accept(setConnsVisitor);
+ ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ grid->accept(setConnsVisitor);
+
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
SPtr<WriteMultiphaseQuantitiesCoProcessor> pp(new WriteMultiphaseQuantitiesCoProcessor(
grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
index 846dd8272..e72f69d04 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
@@ -50,14 +50,21 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::initDataSet()
SPtr<DistributionArray3D> f(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9));
SPtr<DistributionArray3D> h(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9)); // For phase-field
SPtr<DistributionArray3D> h2(new D3Q27EsoTwist3DSplittedVector(nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9)); // For phase-field
- SPtr<PhaseFieldArray3D> divU(new PhaseFieldArray3D( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
+ //SPtr<PhaseFieldArray3D> divU(new PhaseFieldArray3D( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
+ SPtr<PhaseFieldArray3D> divU1(new PhaseFieldArray3D( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure(new CbArray3D<LBMReal, IndexerX3X2X1>( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
pressureOld = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
dataSet->setFdistributions(f);
dataSet->setHdistributions(h); // For phase-field
dataSet->setH2distributions(h2); // For phase-field
- dataSet->setPhaseField(divU);
+ //dataSet->setPhaseField(divU);
+ dataSet->setPhaseField(divU1);
dataSet->setPressureField(pressure);
+
+ phaseField = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.0));
+ phaseField2 = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.0));
+ divU = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
+
}
//////////////////////////////////////////////////////////////////////////
SPtr<LBMKernel> MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::clone()
@@ -83,6 +90,7 @@ SPtr<LBMKernel> MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::clone()
kernel->setIndex(ix1, ix2, ix3);
kernel->setDeltaT(deltaT);
kernel->setGhostLayerWidth(2);
+ dynamicPointerCast<MultiphaseTwoPhaseFieldsPressureFilterLBMKernel>(kernel)->initForcing();
return kernel;
}
@@ -170,14 +178,16 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
+
+
//TODO
//very expensive !!!!!
- CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField(
- new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, -999.0));
- CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField2(
- new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, -999.0));
- CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr divU(
- new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
+ //CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField(
+ // new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, -999.0));
+ // CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField2(
+ // new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, -999.0));
+ // CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr divU(
+ // new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
//#pragma omp parallel for
for (int x3 = minX3-ghostLayerWidth; x3 < maxX3+ghostLayerWidth; x3++) {
@@ -713,7 +723,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
//+WEIGTH[NE] * (((phi2[TE] - phi2[BW]) - (phi2[BE] - phi2[TW])) + ((phi2[TS] - phi2[BN]) + (phi2[TN] - phi2[BS])))) +
//+WEIGTH[N] * (phi2[T] - phi2[B]));
- if (withForcing) {
+ //if (withForcing) {
// muX1 = static_cast<double>(x1-1+ix1*maxX1);
// muX2 = static_cast<double>(x2-1+ix2*maxX2);
// muX3 = static_cast<double>(x3-1+ix3*maxX3);
@@ -722,11 +732,13 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
// forcingX2 = muForcingX2.Eval() + c1o3*drho*dX2_phi * rhoToPhi / rho;//-gradPy/rho;
//forcingX3 = muForcingX3.Eval() + c1o3*drho*dX3_phi * rhoToPhi / rho;//-gradPz/rho;
- muForcingX1.DefineVar("rho",&muRho);
- muForcingX2.DefineVar("rho",&muRho);
- muForcingX3.DefineVar("rho",&muRho);
muRho = rho;
+
+ //muForcingX1.DefineConst("rho", rho);
+ //muForcingX2.DefineConst("rho", rho);
+ //muForcingX3.DefineConst("rho", rho);
+
forcingX1 = muForcingX1.Eval()/rho - gradPx/rho;
forcingX2 = muForcingX2.Eval()/rho - gradPy/rho;
@@ -740,7 +752,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
vvy += forcingX2 * deltaT * 0.5; // Y
vvz += forcingX3 * deltaT * 0.5; // Z
- }
+ //}
///surface tension force
@@ -958,20 +970,20 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
//forcing
///////////////////////////////////////////////////////////////////////////////////////////
- if (withForcing)
- {
- muX1 = static_cast<double>(x1 - 1 + ix1 * maxX1);
- muX2 = static_cast<double>(x2 - 1 + ix2 * maxX2);
- muX3 = static_cast<double>(x3 - 1 + ix3 * maxX3);
-
- //forcingX1 = muForcingX1.Eval();
- //forcingX2 = muForcingX2.Eval();
- //forcingX3 = muForcingX3.Eval();
-
- //vvx += forcingX1 * deltaT * 0.5; // X
- //vvy += forcingX2 * deltaT * 0.5; // Y
- //vvz += forcingX3 * deltaT * 0.5; // Z
- }
+ //if (withForcing)
+ //{
+ // muX1 = static_cast<double>(x1 - 1 + ix1 * maxX1);
+ // muX2 = static_cast<double>(x2 - 1 + ix2 * maxX2);
+ // muX3 = static_cast<double>(x3 - 1 + ix3 * maxX3);
+
+ // //forcingX1 = muForcingX1.Eval();
+ // //forcingX2 = muForcingX2.Eval();
+ // //forcingX3 = muForcingX3.Eval();
+
+ // //vvx += forcingX1 * deltaT * 0.5; // X
+ // //vvy += forcingX2 * deltaT * 0.5; // Y
+ // //vvz += forcingX3 * deltaT * 0.5; // Z
+ //}
LBMReal vx2;
LBMReal vy2;
@@ -1296,9 +1308,9 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
+ (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3 * oMdrho) + c1o27 * oMdrho;
///storing pre collision second moments
- LBMReal mbxx = mfcaa;
- LBMReal mbyy = mfaca;
- LBMReal mbzz = mfaac;
+ LBMReal mbxx = mfcaa - c1o3 * mfaaa;
+ LBMReal mbyy = mfaca - c1o3 * mfaaa;
+ LBMReal mbzz = mfaac - c1o3 * mfaaa;
LBMReal mbxy = mfbba;
LBMReal mbxz = mfbab;
LBMReal mbyz = mfabb;
@@ -1468,9 +1480,12 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
mfaba = -mfaba;
mfaab = -mfaab;
//////////////////////////////////////////////////////////////////////////////////////
- mfbaa += -rho * rhoToPhi * c1o2 * ((mbxx + mfcaa) * dX1_phi + (mbxy + mfbba) * dX2_phi + (mbxz + mfbab) * dX3_phi);
- mfaba += -rho * rhoToPhi * c1o2 * ((mbxy + mfbba) * dX1_phi + (mbyy + mfaca) * dX2_phi + (mbyz + mfabb) * dX3_phi);
- mfaab += -rho * rhoToPhi * c1o2 * ((mbxz + mfbab) * dX1_phi + (mbyz + mfabb) * dX2_phi + (mbzz + mfaac) * dX3_phi);
+ //mfbaa += -rho * rhoToPhi * c1o2 * ((mbxx + mfcaa) * dX1_phi + (mbxy + mfbba) * dX2_phi + (mbxz + mfbab) * dX3_phi);
+ //mfaba += -rho * rhoToPhi * c1o2 * ((mbxy + mfbba) * dX1_phi + (mbyy + mfaca) * dX2_phi + (mbyz + mfabb) * dX3_phi);
+ //mfaab += -rho * rhoToPhi * c1o2 * ((mbxz + mfbab) * dX1_phi + (mbyz + mfabb) * dX2_phi + (mbzz + mfaac) * dX3_phi);
+ mfbaa += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (2 * dxux * dX1_phi + Dxy * dX2_phi + Dxz * dX3_phi) / (rho);
+ mfaba += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (Dxy * dX1_phi + 2 * dyuy * dX2_phi + Dyz * dX3_phi) / (rho);
+ mfaab += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (Dxz * dX1_phi + Dyz * dX2_phi + 2 * dyuy * dX3_phi) / (rho);
////////////////////////////////////////////////////////////////////////////////////
//back
////////////////////////////////////////////////////////////////////////////////////
@@ -3538,4 +3553,28 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::swapDistributions()
LBMKernel::swapDistributions();
dataSet->getHdistributions()->swap();
dataSet->getH2distributions()->swap();
-}
\ No newline at end of file
+}
+
+void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::initForcing()
+{
+ muForcingX1.DefineVar("x1", &muX1); muForcingX1.DefineVar("x2", &muX2); muForcingX1.DefineVar("x3", &muX3);
+ muForcingX2.DefineVar("x1", &muX1); muForcingX2.DefineVar("x2", &muX2); muForcingX2.DefineVar("x3", &muX3);
+ muForcingX3.DefineVar("x1", &muX1); muForcingX3.DefineVar("x2", &muX2); muForcingX3.DefineVar("x3", &muX3);
+
+ muDeltaT = deltaT;
+
+ muForcingX1.DefineVar("dt", &muDeltaT);
+ muForcingX2.DefineVar("dt", &muDeltaT);
+ muForcingX3.DefineVar("dt", &muDeltaT);
+
+ muNu = (1.0 / 3.0) * (1.0 / collFactor - 1.0 / 2.0);
+
+ muForcingX1.DefineVar("nu", &muNu);
+ muForcingX2.DefineVar("nu", &muNu);
+ muForcingX3.DefineVar("nu", &muNu);
+
+ muForcingX1.DefineVar("rho",&muRho);
+ muForcingX2.DefineVar("rho",&muRho);
+ muForcingX3.DefineVar("rho",&muRho);
+
+}
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.h
index 302aaf508..7d20f8210 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.h
@@ -51,10 +51,7 @@ public:
virtual ~MultiphaseTwoPhaseFieldsPressureFilterLBMKernel(void) = default;
void calculate(int step) override;
SPtr<LBMKernel> clone() override;
- void forwardInverseChimeraWithKincompressible(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K, LBMReal oneMinusRho);
- void backwardInverseChimeraWithKincompressible(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K, LBMReal oneMinusRho);
- void forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
- void backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
+
///refactor
//CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure;
@@ -64,6 +61,14 @@ public:
protected:
virtual void initDataSet();
void swapDistributions() override;
+
+ void initForcing();
+
+ void forwardInverseChimeraWithKincompressible(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K, LBMReal oneMinusRho);
+ void backwardInverseChimeraWithKincompressible(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K, LBMReal oneMinusRho);
+ void forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
+ void backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
+
LBMReal f1[D3Q27System::ENDF+1];
CbArray4D<LBMReal,IndexerX4X3X2X1>::CbArray4DPtr localDistributionsF;
@@ -82,6 +87,10 @@ protected:
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressureOld;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField2;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr divU;
+
LBMReal h [D3Q27System::ENDF+1];
LBMReal h2[D3Q27System::ENDF + 1];
LBMReal g [D3Q27System::ENDF+1];
--
GitLab